diff --git a/CMakeLists.txt b/CMakeLists.txt
index 8cd43729f..7152ca841 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -117,12 +117,12 @@ if (WITH_RAWZOR)
     endif (WIN32)
 endif (WITH_RAWZOR)
 
-if (OPTION_OMP)
-	find_package(OpenMP)
-	if (OPENMP_FOUND)
-		set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
-	endif (OPENMP_FOUND)
-endif (OPTION_OMP)
+#if (OPTION_OMP)
+#	find_package(OpenMP)
+#	if (OPENMP_FOUND)
+#		set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+#	endif (OPENMP_FOUND)
+#endif (OPTION_OMP)
 
 add_subdirectory (rtexif)
 add_subdirectory (rtengine)
diff --git a/rtdata/languages/English (UK) b/rtdata/languages/English (UK)
index d992de66b..37b52708d 100644
--- a/rtdata/languages/English (UK)	
+++ b/rtdata/languages/English (UK)	
@@ -363,6 +363,8 @@ PREFERENCES_DATEFORMAT;Date Format
 PREFERENCES_DATEFORMATHINT;<i>You can use the following formatting strings:</i>\n<b>%y</b><i> : year</i>\n<b>%m</b><i> : month</i>\n<b>%d</b><i> : day</i>\n<i>\nFor example, the hungarian date format is:</i>\n<b>%y/%m/%d</b>
 PREFERENCES_DCBENHANCE;Apply DCB enhancment step
 PREFERENCES_DCBITERATIONS;Number of DCB iterations
+#Emil's CA autocorrection
+PREFERENCES_CACORRECTION;Apply CA auto correction
 PREFERENCES_DEFAULTLANG;Default language
 PREFERENCES_DEFAULTTHEME;Default theme
 PREFERENCES_DEMOSAICINGALGO;Demosaicing Algorithm
diff --git a/rtdata/languages/English (US) b/rtdata/languages/English (US)
index b87063e3e..492a1b35c 100644
--- a/rtdata/languages/English (US)	
+++ b/rtdata/languages/English (US)	
@@ -365,6 +365,8 @@ PREFERENCES_DATEFORMAT;Date Format
 PREFERENCES_DATEFORMATHINT;<i>You can use the following formatting strings:</i>\n<b>%y</b><i> : year</i>\n<b>%m</b><i> : month</i>\n<b>%d</b><i> : day</i>\n<i>\nFor example, the hungarian date format is:</i>\n<b>%y/%m/%d</b>
 PREFERENCES_DCBENHANCE;Apply DCB enhancement step
 PREFERENCES_DCBITERATIONS;Number of DCB iterations
+#Emil's CA autocorrection
+PREFERENCES_CACORRECTION;Apply CA auto correction
 PREFERENCES_DEFAULTLANG;Default language
 PREFERENCES_DEFAULTTHEME;Default theme
 PREFERENCES_DEMOSAICINGALGO;Demosaicing Algorithm
diff --git a/rtengine/CA_correct_RT.cc b/rtengine/CA_correct_RT.cc
new file mode 100644
index 000000000..820edab40
--- /dev/null
+++ b/rtengine/CA_correct_RT.cc
@@ -0,0 +1,678 @@
+////////////////////////////////////////////////////////////////
+//
+//		Chromatic Aberration Auto-correction
+//
+//		copyright (c) 2008-2010  Emil Martinec <ejmartin@uchicago.edu>
+//
+//
+// code dated: June 2, 2010
+//
+//	CA_correct_RT.cc is free software: you can redistribute it and/or modify
+//	it under the terms of the GNU General Public License as published by
+//	the Free Software Foundation, either version 3 of the License, or
+//	(at your option) any later version.
+//
+//	This program is distributed in the hope that it will be useful,
+//	but WITHOUT ANY WARRANTY; without even the implied warranty of
+//	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//	GNU General Public License for more details.
+//
+//	You should have received a copy of the GNU General Public License
+//	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+////////////////////////////////////////////////////////////////
+
+
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+int RawImageSource::LinEqSolve(int c, int dir, int nDim, float* pfMatr, float* pfVect, float* pfSolution) 
+{
+//==============================================================================
+// return 1 if system not solving, 0 if system solved
+// nDim - system dimension
+// pfMatr - matrix with coefficients
+// pfVect - vector with free members
+// pfSolution - vector with system solution
+// pfMatr becames trianglular after function call
+// pfVect changes after function call
+//
+// Developer: Henry Guennadi Levkin
+//
+//==============================================================================
+
+	float fMaxElem;
+	float fAcc;
+	
+	int i, j, k, m;
+	
+	for(k=0; k<(nDim-1); k++) {// base row of matrix
+		// search of line with max element
+		fMaxElem = fabs( pfMatr[k*nDim + k] );
+		m = k;
+		for (i=k+1; i<nDim; i++) {
+			if(fMaxElem < fabs(pfMatr[i*nDim + k]) ) {
+				fMaxElem = pfMatr[i*nDim + k];
+				m = i;
+			}
+		}
+		
+		// permutation of base line (index k) and max element line(index m)
+		if(m != k) {
+			for(i=k; i<nDim; i++) {
+				fAcc               = pfMatr[k*nDim + i];
+				pfMatr[k*nDim + i] = pfMatr[m*nDim + i];
+				pfMatr[m*nDim + i] = fAcc;
+			}
+			fAcc = pfVect[k];
+			pfVect[k] = pfVect[m];
+			pfVect[m] = fAcc;
+		}
+		
+		if( pfMatr[k*nDim + k] == 0.) {
+			//linear system has no solution
+			return 1; // needs improvement !!!
+		}
+		
+		// triangulation of matrix with coefficients
+		for(j=(k+1); j<nDim; j++) {// current row of matrix
+			fAcc = - pfMatr[j*nDim + k] / pfMatr[k*nDim + k];
+			for(i=k; i<nDim; i++) {
+				pfMatr[j*nDim + i] = pfMatr[j*nDim + i] + fAcc*pfMatr[k*nDim + i];
+			}
+			pfVect[j] = pfVect[j] + fAcc*pfVect[k]; // free member recalculation
+		}
+	}
+	
+	for(k=(nDim-1); k>=0; k--) {
+		pfSolution[k] = pfVect[k];
+		for(i=(k+1); i<nDim; i++) {
+			pfSolution[k] -= (pfMatr[k*nDim + i]*pfSolution[i]);
+		}
+		pfSolution[k] = pfSolution[k] / pfMatr[k*nDim + k];
+	}
+	
+	return 0;
+}
+//end of linear equation solver
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+void RawImageSource::CA_correct_RT() { 
+
+#define TS 256		// Tile size
+//#define polyord	4	// max order of fit monomial
+//#define numpar 16	// number of fit parameters = SQR(polyord)
+#define border 8
+#define border2	16
+	
+	#define PIX_SORT(a,b) { if ((a)>(b)) {temp=(a);(a)=(b);(b)=temp;} }
+	#define SQR(x) ((x)*(x))
+	
+	// local variables
+	int width=W, height=H;
+	float (*Gtmp);
+	Gtmp = (float (*)) calloc ((height)*(width), sizeof *Gtmp);
+	
+	int polyord=4, numpar=16, numblox[3]={0,0,0};
+
+	//static const int border=8;
+	int rrmin, rrmax, ccmin, ccmax;
+	int top, bottom, left, right, row, col;
+	int rr, cc, rr1, cc1, c, indx, indx1, i, j, k, m, n, dir;
+	int GRBdir[2][3], offset[2][3];
+	int	shifthfloor[3], shiftvfloor[3], shifthceil[3], shiftvceil[3];
+	int vblsz, hblsz, vblock, hblock;
+	//int verbose=1;
+	int res;
+	int v1=TS, v2=2*TS, v3=3*TS, v4=4*TS;//, p1=-TS+1, p2=-2*TS+2, p3=-3*TS+3, m1=TS+1, m2=2*TS+2, m3=3*TS+3;
+	
+	float eps=1e-10;			//tolerance to avoid dividing by zero
+	
+	float	wtu, wtd, wtl, wtr;
+	float	coeff[2][3][3], CAshift[2][3];
+	float	polymat[3][2][1296], shiftmat[3][2][36], fitparams[3][2][36];
+	float	shifthfrac[3], shiftvfrac[3], temp, p[9];
+	float	gdiff, deltgrb, denom, Ginthfloor, Ginthceil, Gint, gradwt;
+	float	grbdiffinthfloor, grbdiffinthceil, grbdiffint;
+	float	blockgave, blockgsqave;
+	
+	static const float gaussg[5] = {0.171582, 0.15839, 0.124594, 0.083518, 0.0477063};
+	static const float gaussrb[3] = {0.332406, 0.241376, 0.0924212};
+	
+	//char		*buffer1;				// vblsz*hblsz*3*2
+	//float		(*blockshifts)[3][2];	// vblsz*hblsz*3*2
+	float blockshifts[10000][3][2]; //fixed memory allocation
+	float blockwt[10000]; //fixed memory allocation
+	
+	
+	char		*buffer;			// TS*TS*16
+	float         (*rgb)[3];		// TS*TS*12
+	float         (*grbdiff);		// TS*TS*4
+	float         (*gshift);		// TS*TS*4
+	float         (*ghpfh);		// TS*TS*4
+	float         (*ghpfv);		// TS*TS*4
+	float         (*rbhpfh);		// TS*TS*4
+	float         (*rbhpfv);		// TS*TS*4
+
+	
+	/* assign working space; this would not be necessary
+	 if the algorithm is part of the larger pre-interpolation processing */
+	buffer = (char *) malloc(36*TS*TS);
+	//merror(buffer,"CA_correct()");
+	memset(buffer,0,36*TS*TS);
+	
+	// rgb array
+	rgb         = (float (*)[3])		buffer;
+	grbdiff		= (float (*))			(buffer +	12*TS*TS);
+	gshift		= (float (*))			(buffer +	16*TS*TS);
+	ghpfh		= (float (*))			(buffer +	20*TS*TS);
+	ghpfv		= (float (*))			(buffer +	24*TS*TS);
+	rbhpfh		= (float (*))			(buffer +	28*TS*TS);
+	rbhpfv		= (float (*))			(buffer +	32*TS*TS);
+	
+	
+	vblsz=ceil((float)(height+border2)/(TS-border2))+2;
+	hblsz=ceil((float)(width+border2)/(TS-border2))+2;
+	
+	/*buffer1 = (char *) malloc(4*vblsz*hblsz*3*2);
+	 merror(buffer1,"CA_correct()");
+	 memset(buffer1,0,4*vblsz*hblsz*3*2);
+	 // block CA shifts
+	 blockshifts	= (float (*)[3][2])		buffer1;*/
+	
+	
+	
+	
+	// Main algorithm: Tile loop
+	//#pragma omp parallel for shared(image,height,width) private(top,left,indx,indx1) schedule(dynamic)
+	for (top=-border, vblock=1; top < height; top += TS-border2, vblock++)
+		for (left=-border, hblock=1; left < width; left += TS-border2, hblock++) {
+			bottom = MIN( top+TS,height+border);
+			right  = MIN(left+TS, width+border);
+			rr1 = bottom - top;
+			cc1 = right - left;
+			//t1_init = clock();
+			// rgb from input CFA data
+			// rgb values should be floating point number between 0 and 1 
+			// after white balance multipliers are applied 
+			if (top<0) {rrmin=border;} else {rrmin=0;}
+			if (left<0) {ccmin=border;} else {ccmin=0;}
+			if (bottom>height) {rrmax=height-top;} else {rrmax=rr1;}
+			if (right>width) {ccmax=width-left;} else {ccmax=cc1;}
+			
+			blockgave=blockgsqave=denom=0;
+			for (rr=rrmin; rr < rrmax; rr++)
+				for (row=rr+top, cc=ccmin; cc < ccmax; cc++) {
+					col = cc+left;
+					c = FC(rr,cc);
+					indx=row*width+col;
+					indx1=rr*TS+cc;
+					rgb[indx1][c] = (ri->data[row][col])/65535.0f;
+				}
+			blockwt[vblock*hblsz+hblock]=0;
+			//blockwt[vblock*hblsz+hblock]=1;
+			
+			// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+			//fill borders
+			if (rrmin>0) {
+				for (rr=0; rr<border; rr++) 
+					for (cc=ccmin; cc<ccmax; cc++) {
+						c = FC(rr,cc);
+						rgb[rr*TS+cc][c] = rgb[(border2-rr)*TS+cc][c];
+					}
+			}
+			if (rrmax<rr1) {
+				for (rr=0; rr<border; rr++) 
+					for (cc=ccmin; cc<ccmax; cc++) {
+						c=FC(rr,cc);
+						rgb[(rrmax+rr)*TS+cc][c] = (ri->data[(height-rr-2)][left+cc])/65535.0f;
+					}
+			}
+			if (ccmin>0) {
+				for (rr=rrmin; rr<rrmax; rr++) 
+					for (cc=0; cc<border; cc++) {
+						c=FC(rr,cc);
+						rgb[rr*TS+cc][c] = rgb[rr*TS+border2-cc][c];
+					}
+			}
+			if (ccmax<cc1) {
+				for (rr=rrmin; rr<rrmax; rr++)
+					for (cc=0; cc<border; cc++) {
+						c=FC(rr,cc);
+						rgb[rr*TS+ccmax+cc][c] = (ri->data[(top+rr)][(width-cc-2)])/65535.0f;
+					}
+			}
+			
+			//also, fill the image corners
+			if (rrmin>0 && ccmin>0) {
+				for (rr=0; rr<border; rr++) 
+					for (cc=0; cc<border; cc++) {
+						c=FC(rr,cc);
+						rgb[(rr)*TS+cc][c] = (ri->data[border2-rr][border2-cc])/65535.0f;
+					}
+			}
+			if (rrmax<rr1 && ccmax<cc1) {
+				for (rr=0; rr<border; rr++) 
+					for (cc=0; cc<border; cc++) {
+						c=FC(rr,cc);
+						rgb[(rrmax+rr)*TS+ccmax+cc][c] = (ri->data[(height-rr-2)][(width-cc-2)])/65535.0f;
+					}
+			}
+			if (rrmin>0 && ccmax<cc1) {
+				for (rr=0; rr<border; rr++) 
+					for (cc=0; cc<border; cc++) {
+						c=FC(rr,cc);
+						rgb[(rr)*TS+ccmax+cc][c] = (ri->data[(border2-rr)][(width-cc-2)])/65535.0f;
+					}
+			}
+			if (rrmax<rr1 && ccmin>0) {
+				for (rr=0; rr<border; rr++) 
+					for (cc=0; cc<border; cc++) {
+						c=FC(rr,cc);
+						rgb[(rrmax+rr)*TS+cc][c] = (ri->data[(height-rr-2)][(border2-cc)])/65535.0f;
+					}
+			}
+			
+			//end of border fill
+			// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+			
+				for (j=0; j<2; j++)
+					for (k=0; k<3; k++)
+						for (c=0; c<3; c+=2) coeff[j][k][c]=0;
+			
+			//end of initialization
+			
+						
+			for (rr=3; rr < rr1-3; rr++)
+				for (row=rr+top, cc=3, indx=rr*TS+cc; cc < cc1-3; cc++, indx++) {
+					col = cc+left;
+					c = FC(rr,cc);
+					
+					if (c!=1) {
+						//compute directional weights using image gradients
+						wtu=1/SQR(eps+fabs(rgb[(rr+1)*TS+cc][1]-rgb[(rr-1)*TS+cc][1])+fabs(rgb[(rr)*TS+cc][c]-rgb[(rr-2)*TS+cc][c])+fabs(rgb[(rr-1)*TS+cc][1]-rgb[(rr-3)*TS+cc][1]));
+						wtd=1/SQR(eps+fabs(rgb[(rr-1)*TS+cc][1]-rgb[(rr+1)*TS+cc][1])+fabs(rgb[(rr)*TS+cc][c]-rgb[(rr+2)*TS+cc][c])+fabs(rgb[(rr+1)*TS+cc][1]-rgb[(rr+3)*TS+cc][1]));
+						wtl=1/SQR(eps+fabs(rgb[(rr)*TS+cc+1][1]-rgb[(rr)*TS+cc-1][1])+fabs(rgb[(rr)*TS+cc][c]-rgb[(rr)*TS+cc-2][c])+fabs(rgb[(rr)*TS+cc-1][1]-rgb[(rr)*TS+cc-3][1]));
+						wtr=1/SQR(eps+fabs(rgb[(rr)*TS+cc-1][1]-rgb[(rr)*TS+cc+1][1])+fabs(rgb[(rr)*TS+cc][c]-rgb[(rr)*TS+cc+2][c])+fabs(rgb[(rr)*TS+cc+1][1]-rgb[(rr)*TS+cc+3][1]));
+
+						
+						//store in rgb array the interpolated G value at R/B grid points using directional weighted average
+						rgb[indx][1]=(wtu*rgb[indx-v1][1]+wtd*rgb[indx+v1][1]+wtl*rgb[indx-1][1]+wtr*rgb[indx+1][1])/(wtu+wtd+wtl+wtr);
+					}
+					if (row>-1 && row<height && col>-1 && col<width)
+						Gtmp[row*width + col] = rgb[indx][1];
+				}
+			
+			//color difference curvatures
+			for (rr=4; rr < rr1-4; rr++)
+				for (cc=4+(FC(rr,2)&1), indx=rr*TS+cc, c = FC(rr,cc); cc < cc1-4; cc+=2, indx+=2) {
+					
+					ghpfv[indx] = 2*rgb[indx][1]-rgb[indx+v2][1]-rgb[indx-v2][1];
+					ghpfh[indx] = 2*rgb[indx][1]-rgb[indx+2][1]-rgb[indx-2][1];
+					rbhpfv[indx] = ghpfv[indx]-(2*rgb[indx][c]-rgb[indx+v2][c]-rgb[indx-v2][c]);
+					rbhpfh[indx] = ghpfh[indx]-(2*rgb[indx][c]-rgb[indx+2][c]-rgb[indx-2][c]);
+					
+				}
+			
+			// along line segments, find the point along each segment that minimizes the color variance
+			// averaged over the tile; evaluate for up/down and left/right away from R/B grid point 
+			for (rr=4; rr < rr1-4; rr++)
+				for (cc=4+(FC(rr,2)&1), indx=rr*TS+cc, c = FC(rr,cc); cc < cc1-4; cc+=2, indx+=2) {
+					
+					//in linear interpolation, color differences are a quadratic function of interpolation position;
+					//solve for the interpolation position that minimizes color difference variance over the tile
+
+					//vertical
+					gdiff=0.625*(rgb[indx+TS][1]-rgb[indx-TS][1])+0.1875*(rgb[indx+TS+1][1]-rgb[indx-TS+1][1]+rgb[indx+TS-1][1]-rgb[indx-TS-1][1]);
+					gradwt=fabs(0.25*rbhpfv[indx]+0.125*(rbhpfv[indx+2]+rbhpfv[indx-2]) );
+					
+					deltgrb=(rgb[indx][c]-rgb[indx][1]);
+					
+					coeff[0][0][c] += gradwt*deltgrb*deltgrb;
+					coeff[0][1][c] += gradwt*gdiff*deltgrb;
+					coeff[0][2][c] += gradwt*gdiff*gdiff;
+
+					//horizontal
+					gdiff=0.625*(rgb[indx+1][1]-rgb[indx-1][1])+0.1875*(rgb[indx+1+TS][1]-rgb[indx-1+TS][1]+rgb[indx+1-TS][1]-rgb[indx-1-TS][1]);
+					gradwt=fabs(0.25*rbhpfh[indx]+0.125*(rbhpfh[indx+v2]+rbhpfh[indx-v2]) );
+					
+					deltgrb=(rgb[indx][c]-rgb[indx][1]);
+					
+					coeff[1][0][c] += gradwt*deltgrb*deltgrb;
+					coeff[1][1][c] += gradwt*gdiff*deltgrb;
+					coeff[1][2][c] += gradwt*gdiff*gdiff;
+
+					
+					//	In Mathematica,
+					//  f[x_]=Expand[Total[Flatten[
+					//  ((1-x) RotateLeft[Gint,shift1]+x RotateLeft[Gint,shift2]-cfapad)^2[[dv;;-1;;2,dh;;-1;;2]]]]];
+					//  extremum = -.5Coefficient[f[x],x]/Coefficient[f[x],x^2]	
+				}
+			for (c=0; c<3; c+=2){
+				for (j=0; j<2; j++) {// vert/hor
+					
+					CAshift[j][c]=coeff[j][1][c]/coeff[j][2][c];
+					blockwt[vblock*hblsz+hblock]= (rr1-8)*(cc1-8)/4;
+
+					//data structure = CAshift[vert/hor][color]
+					//j=0=vert, 1=hor
+					
+					if ((CAshift[j][c])<0) {
+						GRBdir[j][c]=-1;
+					} else {
+						GRBdir[j][c]=1;
+					}	
+					offset[j][c]=floor(CAshift[j][c]);
+					//offset gives NW corner of square containing the min; j=0=vert, 1=hor
+					
+				}//vert/hor
+			}//color
+			
+			
+						
+			/* CAshift[j][c] are the locations  
+			 that minimize color difference variances; 
+			 This is the approximate _optical_ location of the R/B pixels */
+			
+			
+			for (c=0; c<3; c+=2) {
+				//evaluate the shifts to the location that minimizes CA within the tile
+				blockshifts[(vblock)*hblsz+hblock][c][0]=(CAshift[0][c]); //vert CA shift for R/B
+				blockshifts[(vblock)*hblsz+hblock][c][1]=(CAshift[1][c]); //hor CA shift for R/B
+				//data structure: blockshifts[blocknum][R/B][v/h]
+			}
+		}
+	//end of diagnostic pass
+	
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+	//now prepare for CA correction pass
+	//first, fill border blocks of blockshift array
+	for (vblock=1; vblock<vblsz-1; vblock++) {//left and right sides
+		for (c=0; c<3; c+=2) {
+			for (i=0; i<2; i++) {
+				blockshifts[vblock*hblsz][c][i]=blockshifts[(vblock)*hblsz+2][c][i];
+				blockshifts[vblock*hblsz+hblsz-1][c][i]=blockshifts[(vblock)*hblsz+hblsz-3][c][i];
+			}
+		}
+	}
+	for (hblock=0; hblock<hblsz; hblock++) {//top and bottom sides
+		for (c=0; c<3; c+=2) {
+			for (i=0; i<2; i++) {
+				blockshifts[hblock][c][i]=blockshifts[2*hblsz+hblock][c][i];
+				blockshifts[(vblsz-1)*hblsz+hblock][c][i]=blockshifts[(vblsz-3)*hblsz+hblock][c][i];
+			}
+		}
+	}
+	//end of filling border pixels of blockshift array
+	
+	//initialize fit arrays
+	for (i=0; i<1296; i++) {polymat[0][0][i] = polymat[0][1][i] = polymat[2][0][i] = polymat[2][1][i] = 0;}
+	for (i=0; i<36; i++) {shiftmat[0][0][i] = shiftmat[0][1][i] = shiftmat[2][0][i] = shiftmat[2][1][i] = 0;}
+	
+	for (vblock=1; vblock<vblsz-1; vblock++)
+		for (hblock=1; hblock<hblsz-1; hblock++) {
+			// block 3x3 median of blockshifts for robustness
+			for (c=0; c<3; c+=2) {
+				for (dir=0; dir<2; dir++) {
+					p[0] = blockshifts[(vblock-1)*hblsz+hblock-1][c][dir];
+					p[1] = blockshifts[(vblock-1)*hblsz+hblock][c][dir];
+					p[2] = blockshifts[(vblock-1)*hblsz+hblock+1][c][dir];
+					p[3] = blockshifts[(vblock)*hblsz+hblock-1][c][dir];
+					p[4] = blockshifts[(vblock)*hblsz+hblock][c][dir];
+					p[5] = blockshifts[(vblock)*hblsz+hblock+1][c][dir];
+					p[6] = blockshifts[(vblock+1)*hblsz+hblock-1][c][dir];
+					p[7] = blockshifts[(vblock+1)*hblsz+hblock][c][dir];
+					p[8] = blockshifts[(vblock+1)*hblsz+hblock+1][c][dir];
+					PIX_SORT(p[1],p[2]); PIX_SORT(p[4],p[5]); PIX_SORT(p[7],p[8]);
+					PIX_SORT(p[0],p[1]); PIX_SORT(p[3],p[4]); PIX_SORT(p[6],p[7]);
+					PIX_SORT(p[1],p[2]); PIX_SORT(p[4],p[5]); PIX_SORT(p[7],p[8]);
+					PIX_SORT(p[0],p[3]); PIX_SORT(p[5],p[8]); PIX_SORT(p[4],p[7]);
+					PIX_SORT(p[3],p[6]); PIX_SORT(p[1],p[4]); PIX_SORT(p[2],p[5]);
+					PIX_SORT(p[4],p[7]); PIX_SORT(p[4],p[2]); PIX_SORT(p[6],p[4]);
+					PIX_SORT(p[4],p[2]);
+					blockshifts[(vblock)*hblsz+hblock][c][dir] = p[4];
+					if (p[4]<0) {GRBdir[dir][c]=-1;} else {GRBdir[dir][c]=1;}
+				}
+				
+				//now prepare coefficient matrix
+				if (fabs(blockshifts[(vblock)*hblsz+hblock][c][0])>1.0 || fabs(blockshifts[(vblock)*hblsz+hblock][c][1])>1.0) continue;
+				numblox[c] += 1;
+				for (dir=0; dir<2; dir++) {
+					for (i=0; i<polyord; i++)
+						for (j=0; j<polyord; j++) {
+							for (m=0; m<polyord; m++)
+								for (n=0; n<polyord; n++) {
+									polymat[c][dir][numpar*(polyord*i+j)+(polyord*m+n)] += (float)pow(vblock,i+m)*pow(hblock,j+n)*blockwt[vblock*hblsz+hblock];
+								}
+							shiftmat[c][dir][(polyord*i+j)] += (float)pow(vblock,i)*pow(hblock,j)*blockshifts[(vblock)*hblsz+hblock][c][dir]*blockwt[vblock*hblsz+hblock];
+						}//monomials	
+				}//dir 
+				
+			}//c
+		}//blocks
+	
+	numblox[1]=MIN(numblox[0],numblox[2]);
+	//if (numblox[1]<72) {
+	//	polyord=4; numpar=16;
+		if (numblox[1]<32) {
+			polyord=2; numpar=4;
+			if (numblox[1]< 10) return;
+		}
+	//}	
+	
+	//fit parameters to blockshifts
+	for (c=0; c<3; c+=2)
+		for (dir=0; dir<2; dir++) {
+			res = LinEqSolve(c, dir, numpar, polymat[c][dir], shiftmat[c][dir], fitparams[c][dir]);
+			if (res) {
+				//if (verbose) fprintf(stderr,_("problem fitting CA parameters for (c = %d dir = %d)\n"),c,dir);
+				for (i=0; i<numpar; i++) fitparams[c][dir][i]=0;
+			}
+		}
+	//fitparams[5*i+j] gives the coefficients of (vblock^i hblock^j) in a polynomial fit for i,j<=4
+	
+	//end of initialization for CA correction pass
+	
+
+	
+	// Main algorithm: Tile loop
+	//#pragma omp parallel for shared(image,height,width) private(top,left,indx,indx1) schedule(dynamic)
+	for (top=-border, vblock=1; top < height; top += TS-border2, vblock++)
+		for (left=-border, hblock=1; left < width; left += TS-border2, hblock++) {
+			bottom = MIN( top+TS,height+border);
+			right  = MIN(left+TS, width+border);
+			rr1 = bottom - top;
+			cc1 = right - left;
+			//t1_init = clock();
+			// rgb from input CFA data
+			// rgb values should be floating point number between 0 and 1 
+			// after white balance multipliers are applied 
+			if (top<0) {rrmin=border;} else {rrmin=0;}
+			if (left<0) {ccmin=border;} else {ccmin=0;}
+			if (bottom>height) {rrmax=height-top;} else {rrmax=rr1;}
+			if (right>width) {ccmax=width-left;} else {ccmax=cc1;}
+			
+			for (rr=rrmin; rr < rrmax; rr++)
+				for (row=rr+top, cc=ccmin; cc < ccmax; cc++) {
+					col = cc+left;
+					c = FC(rr,cc);
+					indx=row*width+col;
+					indx1=rr*TS+cc;	
+					//rgb[indx1][c] = image[indx][c]/65535.0f;
+					rgb[indx1][c] = (ri->data[row][col])/65535.0f;
+					if ((c&1)==0) rgb[indx1][1] = Gtmp[indx];
+				}
+			// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+			//fill borders
+			if (rrmin>0) {
+				for (rr=0; rr<border; rr++) 
+					for (cc=ccmin; cc<ccmax; cc++) {
+						c = FC(rr,cc);
+						rgb[rr*TS+cc][c] = rgb[(border2-rr)*TS+cc][c];
+						rgb[rr*TS+cc][1] = rgb[(border2-rr)*TS+cc][1];
+					}
+			}
+			if (rrmax<rr1) {
+				for (rr=0; rr<border; rr++) 
+					for (cc=ccmin; cc<ccmax; cc++) {
+						c=FC(rr,cc);
+						rgb[(rrmax+rr)*TS+cc][c] = (ri->data[(height-rr-2)][left+cc])/65535.0f;
+						rgb[(rrmax+rr)*TS+cc][1] = Gtmp[(height-rr-2)*width+left+cc];
+					}
+			}
+			if (ccmin>0) {
+				for (rr=rrmin; rr<rrmax; rr++) 
+					for (cc=0; cc<border; cc++) {
+						c=FC(rr,cc);
+						rgb[rr*TS+cc][c] = rgb[rr*TS+border2-cc][c];
+						rgb[rr*TS+cc][1] = rgb[rr*TS+border2-cc][1];
+					}
+			}
+			if (ccmax<cc1) {
+				for (rr=rrmin; rr<rrmax; rr++)
+					for (cc=0; cc<border; cc++) {
+						c=FC(rr,cc);
+						rgb[rr*TS+ccmax+cc][c] = (ri->data[(top+rr)][(width-cc-2)])/65535.0f;
+						rgb[rr*TS+ccmax+cc][1] = Gtmp[(top+rr)*width+(width-cc-2)];
+					}
+			}
+			
+			//also, fill the image corners
+			if (rrmin>0 && ccmin>0) {
+				for (rr=0; rr<border; rr++) 
+					for (cc=0; cc<border; cc++) {
+						c=FC(rr,cc);
+						rgb[(rr)*TS+cc][c] = (ri->data[border2-rr][border2-cc])/65535.0f;
+						rgb[(rr)*TS+cc][1] = Gtmp[(border2-rr)*width+border2-cc];
+					}
+			}
+			if (rrmax<rr1 && ccmax<cc1) {
+				for (rr=0; rr<border; rr++) 
+					for (cc=0; cc<border; cc++) {
+						c=FC(rr,cc);
+						rgb[(rrmax+rr)*TS+ccmax+cc][c] = (ri->data[(height-rr-2)][(width-cc-2)])/65535.0f;
+						rgb[(rrmax+rr)*TS+ccmax+cc][1] = Gtmp[(height-rr-2)*width+(width-cc-2)];
+					}
+			}
+			if (rrmin>0 && ccmax<cc1) {
+				for (rr=0; rr<border; rr++) 
+					for (cc=0; cc<border; cc++) {
+						c=FC(rr,cc);
+						rgb[(rr)*TS+ccmax+cc][c] = (ri->data[(border2-rr)][(width-cc-2)])/65535.0f;
+						rgb[(rr)*TS+ccmax+cc][1] = Gtmp[(border2-rr)*width+(width-cc-2)];
+					}
+			}
+			if (rrmax<rr1 && ccmin>0) {
+				for (rr=0; rr<border; rr++) 
+					for (cc=0; cc<border; cc++) {
+						c=FC(rr,cc);
+						rgb[(rrmax+rr)*TS+cc][c] = (ri->data[(height-rr-2)][(border2-cc)])/65535.0f;
+						rgb[(rrmax+rr)*TS+cc][1] = Gtmp[(height-rr-2)*width+(border2-cc)];
+					}
+			}
+			
+			//end of border fill
+			// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%			
+			
+			
+			blockshifts[(vblock)*hblsz+hblock][0][0] = blockshifts[(vblock)*hblsz+hblock][0][1] = 0;
+			blockshifts[(vblock)*hblsz+hblock][2][0] = blockshifts[(vblock)*hblsz+hblock][2][1] = 0;
+			for (i=0; i<polyord; i++)
+				for (j=0; j<polyord; j++) {
+					blockshifts[(vblock)*hblsz+hblock][0][0] += (float)pow(vblock,i)*pow(hblock,j)*fitparams[0][0][polyord*i+j];
+					blockshifts[(vblock)*hblsz+hblock][0][1] += (float)pow(vblock,i)*pow(hblock,j)*fitparams[0][1][polyord*i+j];
+					blockshifts[(vblock)*hblsz+hblock][2][0] += (float)pow(vblock,i)*pow(hblock,j)*fitparams[2][0][polyord*i+j];
+					blockshifts[(vblock)*hblsz+hblock][2][1] += (float)pow(vblock,i)*pow(hblock,j)*fitparams[2][1][polyord*i+j];
+				}
+			
+			for (c=0; c<3; c+=2) {
+				
+				//some parameters for the bilinear interpolation
+				shiftvfloor[c]=floor((float)blockshifts[(vblock)*hblsz+hblock][c][0]);
+				shiftvceil[c]=ceil((float)blockshifts[(vblock)*hblsz+hblock][c][0]);
+				shiftvfrac[c]=blockshifts[(vblock)*hblsz+hblock][c][0]-shiftvfloor[c];
+				
+				shifthfloor[c]=floor((float)blockshifts[(vblock)*hblsz+hblock][c][1]);
+				shifthceil[c]=ceil((float)blockshifts[(vblock)*hblsz+hblock][c][1]);
+				shifthfrac[c]=blockshifts[(vblock)*hblsz+hblock][c][1]-shifthfloor[c];
+			}
+			
+			
+			for (rr=4; rr < rr1-4; rr++)
+				for (cc=4+(FC(rr,2)&1), c = FC(rr,cc); cc < cc1-4; cc+=2) {
+					//perform CA correction using color ratios or color differences
+					
+					Ginthfloor=(1-shifthfrac[c])*rgb[(rr+shiftvfloor[c])*TS+cc+shifthfloor[c]][1]+(shifthfrac[c])*rgb[(rr+shiftvfloor[c])*TS+cc+shifthceil[c]][1];
+					Ginthceil=(1-shifthfrac[c])*rgb[(rr+shiftvceil[c])*TS+cc+shifthfloor[c]][1]+(shifthfrac[c])*rgb[(rr+shiftvceil[c])*TS+cc+shifthceil[c]][1];
+					//Gint is blinear interpolation of G at CA shift point
+					Gint=(1-shiftvfrac[c])*Ginthfloor+(shiftvfrac[c])*Ginthceil;
+					
+					//determine R/B at grid points using color differences at shift point plus interpolated G value at grid point
+					//but first we need to interpolate G-R/G-B to grid points...
+					grbdiff[(rr)*TS+cc]=Gint-rgb[(rr)*TS+cc][c];
+					
+				}
+			
+			for (rr=8; rr < rr1-8; rr++)
+				for (cc=8+(FC(rr,2)&1), c = FC(rr,cc); cc < cc1-8; cc+=2) {
+					//interpolate color difference from optical R/B locations to grid locations
+					
+					grbdiffinthfloor=(1-shifthfrac[c]/2)*grbdiff[(rr)*TS+cc]+(shifthfrac[c]/2)*grbdiff[(rr)*TS+cc-2*GRBdir[1][c]];
+					grbdiffinthceil=(1-shifthfrac[c]/2)*grbdiff[(rr-2*GRBdir[0][c])*TS+cc]+(shifthfrac[c]/2)*grbdiff[(rr-2*GRBdir[0][c])*TS+cc-2*GRBdir[1][c]];
+					//grbdiffint is bilinear interpolation of G-R/G-B at grid point
+					grbdiffint=(1-shiftvfrac[c]/2)*grbdiffinthfloor+(shiftvfrac[c]/2)*grbdiffinthceil;
+					
+					//now determine R/B at grid points using interpolated color differences and interpolated G value at grid point
+					Gint=rgb[(rr)*TS+cc][1]-grbdiffint;
+					
+					if (fabs(Gint-rgb[(rr)*TS+cc][c])<0.25*(Gint+rgb[(rr)*TS+cc][c])) {
+						rgb[(rr)*TS+cc][c]=Gint;
+					} else {
+						
+						//gradient weights using difference from G at CA shift points and G at grid points
+						p[0]=1/(eps+fabs(rgb[(rr)*TS+cc][1]-gshift[(rr)*TS+cc]));
+						p[1]=1/(eps+fabs(rgb[(rr)*TS+cc][1]-gshift[(rr)*TS+cc-2*GRBdir[1][c]]));
+						p[2]=1/(eps+fabs(rgb[(rr)*TS+cc][1]-gshift[(rr-2*GRBdir[0][c])*TS+cc]));
+						p[3]=1/(eps+fabs(rgb[(rr)*TS+cc][1]-gshift[(rr-2*GRBdir[0][c])*TS+cc-2*GRBdir[1][c]]));
+						
+						grbdiffint = (p[0]*grbdiff[(rr)*TS+cc]+p[1]*grbdiff[(rr)*TS+cc-2*GRBdir[1][c]]+ \
+									  p[2]*grbdiff[(rr-2*GRBdir[0][c])*TS+cc]+p[3]*grbdiff[(rr-2*GRBdir[0][c])*TS+cc-2*GRBdir[1][c]])/(p[0]+p[1]+p[2]+p[3]);
+						
+						//there is an assumption that the grid point is no more than 2 pixels from the optical point; perhaps should correct for this???
+						//now determine R/B at grid points using interpolated color differences and interpolated G value at grid point
+						rgb[(rr)*TS+cc][c]=rgb[(rr)*TS+cc][1]-grbdiffint;
+						
+					}
+				}
+			
+			// copy CA corrected results back to image matrix
+			for (rr=border; rr < rr1-border; rr++)
+				for (row=rr+top, cc=border+(FC(rr,2)&1); cc < cc1-border; cc+=2) {
+					col = cc + left;
+					indx = row*width + col;
+					c = FC(row,col);
+					 
+					ri->data[row][col] = CLIP((int)(65535.0f*rgb[(rr)*TS+cc][c] + 0.5f));
+
+				} 
+		}
+	
+	// clean up
+	free(buffer);
+	free(Gtmp);
+	//free(buffer1);
+	
+
+	
+#undef TS
+#undef polyord
+#undef numpar
+#undef border
+#undef PIX_SORT
+#undef SQR
+
+}
+
+
+
diff --git a/rtengine/rawimagesource.cc b/rtengine/rawimagesource.cc
index 4f6c23d36..8ef9b8a70 100644
--- a/rtengine/rawimagesource.cc
+++ b/rtengine/rawimagesource.cc
@@ -772,6 +772,20 @@ int RawImageSource::load (Glib::ustring fname) {
         delete corr_data;
 */                        
     }
+	
+	
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+//Emil's CA auto correction
+	
+	if (settings->ca_autocorrect) {
+		if (plistener) {
+			plistener->setProgressStr ("CA Auto Correction...");
+			plistener->setProgress (0.0);
+		}
+		
+		CA_correct_RT();
+	}
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
     if (ri->filters) {
         // demosaic
@@ -3127,10 +3141,11 @@ void RawImageSource::dcb_demosaic(int iterations, int dcb_enhance)
         free(image);
 }
 	
-	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-	//Emil's code for AMaZE
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+//Emil's code for AMaZE
 #include "amaze_interpolate_RT.cc"//AMaZE demosaic	
-	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+#include "CA_correct_RT.cc"//Emil's CA auto correction
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
 } /* namespace */
diff --git a/rtengine/rawimagesource.h b/rtengine/rawimagesource.h
index 0a9c7d31e..bab0bafd4 100644
--- a/rtengine/rawimagesource.h
+++ b/rtengine/rawimagesource.h
@@ -135,6 +135,9 @@ class RawImageSource : public ImageSource {
         inline  void interpolate_row_rb     (unsigned short* ar, unsigned short* ab, unsigned short* pg, unsigned short* cg, unsigned short* ng, int i);
         inline  void interpolate_row_rb_mul_pp (unsigned short* ar, unsigned short* ab, unsigned short* pg, unsigned short* cg, unsigned short* ng, int i, double r_mul, double g_mul, double b_mul, int x1, int width, int skip);
 
+		int	LinEqSolve			(int c, int dir, int nDim, float* pfMatr, float* pfVect, float* pfSolution);//Emil's CA auto correction
+		void CA_correct_RT		();
+	
         void    eahd_demosaic   ();
         void    hphd_demosaic   ();
         void    vng4_demosaic   ();
diff --git a/rtengine/settings.h b/rtengine/settings.h
index bfe53173e..cd2dffa07 100644
--- a/rtengine/settings.h
+++ b/rtengine/settings.h
@@ -33,6 +33,7 @@ namespace rtengine {
             bool            verbose;
 	    int dcb_iterations;	// number of dcb iterations
 	    bool dcb_enhance;	// whether to do image refinment
+		bool ca_autocorrect;	// Emil's auto CA correction
 
         /** Creates a new instance of Settings.
           * @return a pointer to the new Settings instance. */
diff --git a/rtgui/options.cc b/rtgui/options.cc
index 01d7e36b8..c28542616 100644
--- a/rtgui/options.cc
+++ b/rtgui/options.cc
@@ -120,6 +120,7 @@ void Options::setDefaults () {
     
     rtSettings.dualThreadEnabled = true;
     rtSettings.demosaicMethod = "amaze";//Emil's code for AMaZE
+	rtSettings.ca_autocorrect = false;//Emil's CA correction
     rtSettings.colorCorrectionSteps = 0;
     rtSettings.dcb_iterations = 2;
     rtSettings.dcb_enhance = true;
@@ -384,6 +385,7 @@ int Options::saveToFile (Glib::ustring fname) {
     keyFile.set_integer ("Algorithms", "ColorCorrection", rtSettings.colorCorrectionSteps);
     keyFile.set_integer ("Algorithms", "DCBIterations", rtSettings.dcb_iterations);
     keyFile.set_boolean ("Algorithms", "DCBEnhance", rtSettings.dcb_enhance);
+	keyFile.set_boolean ("Algorithms", "CACorrect", rtSettings.ca_autocorrect);//Emil's CA correction
     
     keyFile.set_integer ("Crop Settings", "DPI", cropDPI);
 
diff --git a/rtgui/preferences.cc b/rtgui/preferences.cc
index acf7ffd3b..1a3ed53b0 100644
--- a/rtgui/preferences.cc
+++ b/rtgui/preferences.cc
@@ -286,11 +286,14 @@ Gtk::Widget* Preferences::getProcParamsPanel () {
     hb13->pack_start (*dcbIterations);
 
     dcbEnhance = Gtk::manage(new Gtk::CheckButton((M("PREFERENCES_DCBENHANCE"))));
+	
+	caAutoCorrect = Gtk::manage(new Gtk::CheckButton((M("PREFERENCES_CACORRECTION"))));//Emil's CA correction
 
     fdb->pack_start (*hb11, Gtk::PACK_SHRINK, 4);
     fdb->pack_start (*hb12, Gtk::PACK_SHRINK, 4);
     fdb->pack_start (*hb13, Gtk::PACK_SHRINK, 4);
     fdb->pack_start (*dcbEnhance, Gtk::PACK_SHRINK, 4);
+	fdb->pack_start (*caAutoCorrect, Gtk::PACK_SHRINK, 4);//Emil's CA correction
     mvbpp->pack_start (*fdem, Gtk::PACK_SHRINK, 4);
     mvbpp->set_border_width (4);
     //  drlab->set_size_request (drimg->get_width(), -1);
@@ -718,6 +721,7 @@ void Preferences::storePreferences () {
         moptions.rtSettings.demosaicMethod = "ahd";
     moptions.rtSettings.dcb_iterations=(int)dcbIterations->get_value();
     moptions.rtSettings.dcb_enhance=dcbEnhance->get_active();
+	moptions.rtSettings.ca_autocorrect=caAutoCorrect->get_active();//Emil's CA correction
 
     if (sdcurrent->get_active ()) 
         moptions.startupDir = STARTUPDIR_CURRENT;
@@ -816,6 +820,7 @@ void Preferences::fillPreferences () {
     dcbEnhance->set_sensitive(moptions.rtSettings.demosaicMethod=="dcb");
     dcbIterations->set_sensitive(moptions.rtSettings.demosaicMethod=="dcb");
     dcbIterationsLabel->set_sensitive(moptions.rtSettings.demosaicMethod=="dcb");
+	caAutoCorrect->set_active(moptions.rtSettings.ca_autocorrect);//Emil's CA Auto Correction
 
     if (moptions.startupDir==STARTUPDIR_CURRENT) 
         sdcurrent->set_active ();
diff --git a/rtgui/preferences.h b/rtgui/preferences.h
index ebf2d9956..b721cb01e 100644
--- a/rtgui/preferences.h
+++ b/rtgui/preferences.h
@@ -75,6 +75,7 @@ class Preferences : public Gtk::Dialog {
     Gtk::Label* dcbIterationsLabel;
     Gtk::SpinButton*  dcbIterations;
     Gtk::CheckButton* dcbEnhance;
+	Gtk::CheckButton* caAutoCorrect;//Emil's CA correction
     Gtk::FileChooserButton* iccDir;
     Gtk::FileChooserButton* monProfile;