diff --git a/rtengine/PF_correct_RT.cc b/rtengine/PF_correct_RT.cc
index c9f28d6c5..a084a9578 100644
--- a/rtengine/PF_correct_RT.cc
+++ b/rtengine/PF_correct_RT.cc
@@ -39,7 +39,7 @@
 
 namespace rtengine {
 
-void ImProcFunctions::PF_correct_RT(LabImage * src, LabImage * dst, double radius, int thresh, bool edges) { 
+void ImProcFunctions::PF_correct_RT(LabImage * src, LabImage * dst, double radius, int thresh) { 
 	
 	float threshsqr = SQR(thresh);
 	int halfwin = ceil(2*radius)+1;
diff --git a/rtengine/curves.cc b/rtengine/curves.cc
index 301d0664f..08113c625 100644
--- a/rtengine/curves.cc
+++ b/rtengine/curves.cc
@@ -262,8 +262,8 @@ namespace rtengine {
 		
 		float exp_scale = a;
 		float scale = 65536.0;
-		float comp = (ecomp+1.0)*hlcompr/100.0;
-		float shoulder = ((scale/exp_scale)*(hlcomprthresh/200.0))+0.1;
+		float comp = (MAX(0,ecomp) + 1.0)*hlcompr/100.0;
+		float shoulder = ((scale/MAX(1,exp_scale))*(hlcomprthresh/200.0))+0.1;
 		//printf("shoulder = %e\n",shoulder);
 		//printf ("exp_scale= %f comp= %f def_mul=%f a= %f \n",exp_scale,comp,def_mul,a);
 		
@@ -272,14 +272,6 @@ namespace rtengine {
 			// change to [0,1] range
 			float val = (float)i-shoulder;
 			
-			// apply default multiplier (that is >1 if highlight recovery is on)
-			// val *= def_mul;
-			
-			// apply base curve, thus, exposure compensation and black point with shadow and highlight protection
-			//val = basecurve (val*def_mul, a, 0, def_mul, hlcompr/100.0, 0);
-			
-			//hlCurve[i] = (65535.0 * CLIPD(val));
-			
 			if (comp>0.0)
 			{
 				if (val>0.0) {
@@ -313,7 +305,6 @@ namespace rtengine {
 				val = gamma (val, gamma_, start, slope, mul, add);
 			
 			// apply brightness curve
-			//val = brightness (val, br/100.0);
 			val = brightcurve->getVal (val);
 			
 			// store result in a temporary array
diff --git a/rtengine/curves.h b/rtengine/curves.h
index 75597e14a..eb3a5fd51 100644
--- a/rtengine/curves.h
+++ b/rtengine/curves.h
@@ -92,6 +92,7 @@ class CurveFactory {
         return 1.0 - hr + hr*baseu((x-x1)/hr, m, 0.3*hr);
     }
 	static inline double clower2 (double x, double m, double sr) {
+		//curve for b<0; starts with positive slope and then rolls over toward straight line to x=y=1
 		float x1 = sr/1.5 + 0.00001;
 		float y1 = 1-(1-x1)*m;
 		if (x>x1 || sr<0.001) 
@@ -99,14 +100,15 @@ class CurveFactory {
 		else
 			return y1+m*(x-x1)-(1-m)*SQR(SQR(1-x/x1));
 	}
-    // tone curve base. a: slope (from exp.comp.), b: black, D: max. x value (can be>1), hr,sr: highlight,shadow recovery
+    // tone curve base. a: slope (from exp.comp.), b: black point normalized by 65535, 
+	// D: max. x value (can be>1), hr,sr: highlight,shadow recovery
     static inline double basecurve (double x, double a, double b, double D, double hr, double sr) { 
         if (b<0) {
-			double m = 0.5;
-			double slope = 1.0+b;
-			double y = -b+m*slope;
+			double m = 0.5;//midpoint
+			double slope = 1.0+b;//slope of straight line between (0,-b) and (1,1)
+			double y = -b+m*slope;//value at midpoint
 			if (x>m) 
-				return y + (x - m)*slope;
+				return y + (x - m)*slope;//value on straight line between (m,y) and (1,1)
 			else 
 				return y*clower2(x/m, slope*m/y, 2.0-sr);
 		} else {
@@ -121,22 +123,6 @@ class CurveFactory {
 				return y+(x-m)*slope;
 		}
     }
-    // brightness curve at point x, only positive amount it supported
-    static inline double brightnessbase (double x, double amount) {
-        if (x<0.5)
-            return x + amount*cupper(2.0*x, 4.5, 0.0)/3.0;
-        else
-            return x + amount*cupper(2.0-2.0*x, 1.5, 0.0)/3.0;
-    }
-    // brightness curve at point x, positive negative and zero amount are supported
-    static inline double brightness (double x, double amount) {
-        if (amount==0.0)
-            return x;
-        else if (amount>0.0)
-            return brightnessbase (x, amount);
-        else 
-            return 1.0 - brightnessbase (1.0-x, -amount);
-    }
 	
 
   public:
diff --git a/rtengine/expo_before_b.cc b/rtengine/expo_before_b.cc
index cab685114..371a76ac2 100644
--- a/rtengine/expo_before_b.cc
+++ b/rtengine/expo_before_b.cc
@@ -26,96 +26,107 @@
 //   	     use fast-demo(provisional) from Emil Martinec
 //			 inspired from work Guillermo Luijk and Manuel LLorens(Perfectraw)
 //
-  // This function uses parameters:
+// This function uses parameters:
 //      exposure (linear): 2^(-8..0..8): currently 0.5 +3
-  //      preserve (log)   : 0..8 : currently 0.1 1
+//      preserve (log)   : 0..8 : currently 0.1 1
 
 void RawImageSource::processRawWhitepoint(float expos, float preser) {
-      	MyTime t1e,t2e;
-    	t1e.set();
-
+	MyTime t1e,t2e;
+	t1e.set();
+	
 	int width=W, height=H;
-
+	
     // exposure correction inspired from G.Luijk
     if (fabs(preser)<0.001) {	
         // No highlight protection - simple mutiplication
+		for (int c=0; c<4; c++) {
+			chmax[c] *= expos;
+		}
 #pragma omp parallel for shared(expos)
         for (int row=0;row<height;row++)
-            for (int col=0;col<width;col++)
+            for (int col=0;col<width;col++) {
                 rawData[row][col] *= expos;
+			}
     } else {
         // save old image as it's overwritten by demosaic
         float** imgd = allocArray< float >(W,H); 
-	 
+		
         // with memcpy it's faster than for (...)
         for (int i=0; i<H; i++) memcpy (imgd[i], rawData[i], W*sizeof(**imgd));
-
+		
         // Demosaic to calc luminosity
         fast_demosaic (0,0,W,H);
-
-// calculate CIE luminosity
+		
+		// calculate CIE luminosity
         float* luminosity = (float *) new float[width*height];
-
+		
 #pragma omp parallel default(shared)  
-{
+		{
             // CIE luminosity
 #pragma omp for  
-	 for(int row=0;row<height;row++)
-		for(int col=0;col<width;col++)
+			for(int row=0;row<height;row++)
+				for(int col=0;col<width;col++)
                     luminosity[row*width+col] = 
                     0.299*(float)red[row][col] + 0.587*(float)green[row][col] + 0.114*(float)blue[row][col]; 
-}	
-
+		}	
+		
         // restore image destroyed by demosaic
         for (int i=0; i<H; i++) memcpy (rawData[i], imgd[i], W*sizeof(**imgd));
         freeArray<float>(imgd, H);
-			
+		
         // Find maximum to adjust LUTs. New float engines clips only at the very end
         int maxVal=0;
-	 for(int row=0;row<height;row++)
-		for(int col=0;col<width;col++)
+		for(int row=0;row<height;row++)
+			for(int col=0;col<width;col++)
                 if (rawData[row][col]>maxVal) maxVal = rawData[row][col];
-
-    // Exposure correction with highlight preservation
+		
+		// Exposure correction with highlight preservation
         LUTf lut(maxVal+1);
-    if(expos>1){
+		if(expos>1){
             float K = (float) maxVal / expos*exp(-preser*log(2.0));
             for (int j=0;j<=maxVal;j++) 
                 lut[(int)j]=(((float)maxVal-K*expos)/((float)maxVal-K)*(j-maxVal)+(float) maxVal) / j;
-
+			
+			for (int c=0; c<4; c++) {
+				chmax[c] *= expos;
+			}
+			
 #pragma omp parallel for  shared(expos)
-	 for(int row=0;row<height;row++)
-		for(int col=0;col<width;col++){
+			for(int row=0;row<height;row++)
+				for(int col=0;col<width;col++){
                     if (luminosity[row*width + col] < K) {
                         rawData[row][col] *= expos;
                     } else {
                         rawData[row][col] *= lut[luminosity[row*width+col]];
-			}
+					}
                 }
         } else {
             // Negative exposure
-      float EV=log(expos)/log(2.0);                              // Convert exp. linear to EV
+			float EV=log(expos)/log(2.0);                              // Convert exp. linear to EV
             float K = (float)maxVal * exp(-preser * log(2.0));
-
+			
             for (int j=0;j<=maxVal;j++) 
                 lut[(int)j] = exp(EV*((float)maxVal-j) / ((float)maxVal-K) * log(2.0));
-
+			
 #pragma omp parallel for  shared(expos)	  
-	 for(int row=0;row<height;row++)
-		for(int col=0;col<width;col++){
+			for(int row=0;row<height;row++)
+				for(int col=0;col<width;col++){
                     if (luminosity[row*width+col]<K) {
                         rawData[row][col] *= expos;
                     } else {
                         rawData[row][col] *= lut[luminosity[row*width+col]];
-    }	
-}
-
+					}	
+				}
+			for (int c=0; c<4; c++) {
+				chmax[c] *= expos;
+			}
         }	
-
+		
+		
         delete[] luminosity;
     }
-		        t2e.set();
-        if( settings->verbose )
-           printf("Exposure before  %d usec\n", t2e.etime(t1e));
-
+	t2e.set();
+	if( settings->verbose )
+		printf("Exposure before  %d usec\n", t2e.etime(t1e));
+	
 }
\ No newline at end of file
diff --git a/rtengine/hilite_recon.cc b/rtengine/hilite_recon.cc
new file mode 100644
index 000000000..0df494ea1
--- /dev/null
+++ b/rtengine/hilite_recon.cc
@@ -0,0 +1,602 @@
+////////////////////////////////////////////////////////////////
+//
+//		Highlight reconstruction
+//
+//		copyright (c) 2008-2011  Emil Martinec <ejmartin@uchicago.edu>
+//
+//
+// code dated: June 16, 2011
+//
+//	hilite_recon.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/>.
+//
+////////////////////////////////////////////////////////////////
+
+//#include <rtengine.h>
+#include <cstddef>
+#include <math.h>
+#include <curves.h>
+#include <array2D.h>
+#include <improcfun.h>
+#include <rawimagesource.h>
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+#define SQR(x) ((x)*(x))
+
+#define MAX(a,b) ((a)<(b)?(b):(a))
+#define MIN(a,b) ((a)>(b)?(b):(a))
+
+
+//#include "RGBdefringe.cc"
+
+//namespace rtengine {
+
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+void RawImageSource::boxblur2(float** src, float** dst, int H, int W, int box ) 
+{
+	
+	array2D<float> temp(W,H,ARRAY2D_CLEAR_DATA);
+	
+	//box blur image channel; box size = 2*box+1
+	//horizontal blur
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+	for (int row = 0; row < H; row++) {
+		int len = box + 1;
+		temp[row][0] = src[row][0]/len;
+		for (int j=1; j<=box; j++) {
+			temp[row][0] += src[row][j]/len;
+		}
+		for (int col=1; col<=box; col++) {
+			temp[row][col] = (temp[row][col-1]*len + src[row][col+box])/(len+1);
+			len ++;
+		}
+		for (int col = box+1; col < W-box; col++) {
+			temp[row][col] = temp[row][col-1] + (src[row][col+box] - src[row][col-box-1])/len;
+		}
+		for (int col=W-box; col<W; col++) {
+			temp[row][col] = (temp[row][col-1]*len - src[row][col-box-1])/(len-1);
+			len --;
+		}
+	}
+	
+	//vertical blur
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+	for (int col = 0; col < W; col++) {
+		int len = box + 1;
+		dst[0][col] = temp[0][col]/len;
+		for (int i=1; i<=box; i++) {
+			dst[0][col] += temp[i][col]/len;
+		}
+		for (int row=1; row<=box; row++) {
+			dst[row][col] = (dst[(row-1)][col]*len + temp[(row+box)][col])/(len+1);
+			len ++;
+		}
+		for (int row = box+1; row < H-box; row++) {
+			dst[row][col] = dst[(row-1)][col] + (temp[(row+box)][col] - temp[(row-box-1)][col])/len;
+		}
+		for (int row=H-box; row<H; row++) {
+			dst[row][col] = (dst[(row-1)][col]*len - temp[(row-box-1)][col])/(len-1);
+			len --;
+		}
+	}
+	
+}
+
+void RawImageSource::boxblur_resamp(float **src, float **dst, float & max, int H, int W, int box, int samp ) 
+{
+	
+	array2D<float> temp(W,H,ARRAY2D_CLEAR_DATA);
+	array2D<float> temp1(W,H,ARRAY2D_CLEAR_DATA);
+	
+	float maxtmp=0;
+
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+	//box blur image channel; box size = 2*box+1
+	//horizontal blur
+	for (int row = 0; row < H; row++) {
+		int len = box + 1;
+		temp[row][0] = src[row][0]/len;
+		maxtmp = MAX(maxtmp,src[row][0]);
+		for (int j=1; j<=box; j++) {
+			temp[row][0] += src[row][j]/len;
+			maxtmp = MAX(maxtmp,src[row][j]);
+		}
+		for (int col=1; col<=box; col++) {
+			temp[row][col] = (temp[row][col-1]*len + src[row][col+box])/(len+1);
+			maxtmp = MAX(maxtmp,src[row][col]);
+			len ++;
+		}
+		for (int col = box+1; col < W-box; col++) {
+			temp[row][col] = temp[row][col-1] + (src[row][col+box] - src[row][col-box-1])/len;
+			maxtmp = MAX(maxtmp,src[row][col]);
+		}
+		for (int col=W-box; col<W; col++) {
+			temp[row][col] = (temp[row][col-1]*len - src[row][col-box-1])/(len-1);
+			maxtmp = MAX(maxtmp,src[row][col]);
+			len --;
+		}
+	}
+	
+	max = maxtmp;
+
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+	//vertical blur
+	for (int col = 0; col < W; col+=samp) {
+		int len = box + 1;
+		temp1[0][col] = temp[0][col]/len;
+		for (int i=1; i<=box; i++) {
+			temp1[0][col] += temp[i][col]/len;
+		}
+		for (int row=1; row<=box; row++) {
+			temp1[row][col] = (temp1[(row-1)][col]*len + temp[(row+box)][col])/(len+1);
+			len ++;
+		}
+		for (int row = box+1; row < H-box; row++) {
+			temp1[row][col] = temp1[(row-1)][col] + (temp[(row+box)][col] - temp[(row-box-1)][col])/len;
+		}
+		for (int row=H-box; row<H; row++) {
+			temp1[row][col] = (temp1[(row-1)][col]*len - temp[(row-box-1)][col])/(len-1);
+			len --;
+		}
+	}
+	
+	for (int row = 0; row < (H-(H%samp))/samp; row++)
+		for (int col = 0; col < (W-(W%samp))/samp; col++)
+			 dst[row][col] = temp1[samp*row][samp*col];
+	
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+	
+void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** blue)
+{
+	
+	if (plistener) {
+		plistener->setProgressStr ("HL reconstruction...");
+		plistener->setProgress (0.0);
+	}
+	float progress = 0.0;
+	
+	
+	int height = H;
+	int width = W;
+	
+	const int range = 2;
+	const int pitch = 4;
+
+	
+	int hfh = (height-(height%pitch))/pitch;
+	int hfw = (width-(width%pitch))/pitch;
+	
+	static const int numdirs = 4;				
+	
+	static const float threshpct = 0.25;
+	static const float fixthreshpct = 0.7;
+	static const float maxpct = 0.95;
+	
+	float max[3], thresh[3], fixthresh[3], norm[3];
+		
+	//float red1, green1, blue1;//diagnostic
+	float chmaxalt[4]={0,0,0,0};//diagnostic
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	//halfsize demosaic
+	
+	
+	multi_array2D<float,3> hfsize (hfw,hfh,ARRAY2D_CLEAR_DATA);
+	
+	boxblur_resamp(red,hfsize[0],chmaxalt[0],height,width,range,pitch);
+	boxblur_resamp(green,hfsize[1],chmaxalt[1],height,width,range,pitch);
+	boxblur_resamp(blue,hfsize[2],chmaxalt[2],height,width,range,pitch);
+	
+	if(plistener) plistener->setProgress(0.10);
+
+
+	
+	//blur image
+	//for (int m=0; m<3; m++)
+	//	boxblur2(hfsize[m],hfsizeblur[m],hfh,hfw,3);
+	
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+	for (int c=0; c<3; c++) {
+		thresh[c] = chmax[c]*threshpct;
+		fixthresh[c] = chmax[c]*fixthreshpct;
+		max[c] = chmax[c]*maxpct;//MIN(MIN(chmax[0],chmax[1]),chmax[2])*maxpct;
+		norm[c] = 1.0/(max[c]-fixthresh[c]);
+	}
+	float whitept = MAX(MAX(max[0],max[1]),max[2]);
+	float clippt  = MIN(MIN(max[0],max[1]),max[2]);
+	
+	
+	float sat = ri->get_white();
+	float black = ri->get_black();
+	sat -= black;
+	float camwb[4];
+	for (int c=0; c<4; c++) camwb[c]=ri->get_cam_mul(c);
+	float min = MIN(MIN(camwb[0],camwb[1]),camwb[2]);
+
+	
+	multi_array2D<float,5> hilite_full(width,height,ARRAY2D_CLEAR_DATA);
+
+	multi_array2D<float,4> hilite(hfw,hfh,ARRAY2D_CLEAR_DATA);
+	
+	multi_array2D<float,4*numdirs> hilite_dir(hfw,hfh,ARRAY2D_CLEAR_DATA);
+		
+	
+	// set up which pixels are clipped or near clipping
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+	for (int i=0; i<height; i++) {
+		for (int j=0; j<width; j++) {
+			
+			//if one or more channels is highlight but none are blown, add to highlight accumulator
+			
+			if ((red[i][j]>thresh[0] || green[i][j]>thresh[1] || blue[i][j]>thresh[2]) && \
+				(red[i][j]<max[0] && green[i][j]<max[1] && blue[i][j]<max[2])) {
+
+				hilite_full[0][i][j] = red[i][j];
+				hilite_full[1][i][j] = green[i][j];
+				hilite_full[2][i][j] = blue[i][j];
+				hilite_full[3][i][j] = 1;
+				hilite_full[4][i][j] = 1;
+			}
+			//if (i%100==0 && j%100==0)
+			//	printf("row=%d  col=%d  r=%f  g=%f  b=%f hilite=%f  \n",i,j,hilite_full[0][i][j],hilite_full[1][i][j],hilite_full[2][i][j],hilite_full[3][i][j]);
+		}
+	}//end of filling highlight array
+	
+	if(plistener) plistener->setProgress(0.25);
+
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	//blur highlight data
+	boxblur2(hilite_full[4],hilite_full[4],height,width,1);
+	
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+	for (int i=0; i<height; i++) {
+		for (int j=0; j<width; j++) {
+						
+			if (hilite_full[4][i][j]>0.00001 && hilite_full[4][i][j]<0.95) {
+				//too near an edge, could risk using CA affected pixels, therefore omit
+				hilite_full[0][i][j] = hilite_full[1][i][j] = hilite_full[2][i][j] = hilite_full[3][i][j] = 0;
+			}
+		}
+	}
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	// blur and resample highlight data; range=size of blur, pitch=sample spacing
+	for (int m=0; m<4; m++)
+		boxblur_resamp(hilite_full[m],hilite[m],chmaxalt[m],height,width,range,pitch);
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+	//blur highlights
+	//for (int m=0; m<4; m++)
+	//	boxblur2(hilite[m],hilite[m],hfh,hfw,4);
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+	if(plistener) plistener->setProgress(0.50);
+	
+	LUTf invfn(0x10000);
+	
+	for (int i=0; i<0x10000; i++)
+		invfn[i] = 1.0/(1.0+i);
+	
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+	
+	//fill gaps in highlight map by directional extension
+	//raster scan from four corners
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+	for (int j=1; j<hfw-1; j++) 
+		for (int i=2; i<hfh-2; i++) {
+			//from left 
+			if (hilite[3][i][j]>0.01) {
+				for (int c=0; c<4; c++) {
+					hilite_dir[c][i][j] = hilite[c][i][j]/hilite[3][i][j];
+				}
+			} else {
+				for (int c=0; c<4; c++) {
+					hilite_dir[c][i][j] = 0.1*((hilite_dir[0+c][i-2][j-1]+hilite_dir[0+c][i-1][j-1]+hilite_dir[0+c][i][j-1]+hilite_dir[0+c][i+1][j-1]+hilite_dir[0+c][i+2][j-1])/ \
+											  (hilite_dir[0+3][i-2][j-1]+hilite_dir[0+3][i-1][j-1]+hilite_dir[0+3][i][j-1]+hilite_dir[0+3][i+1][j-1]+hilite_dir[0+3][i+2][j-1]+0.00001));
+					hilite_dir[4+c][i][j+1]  += hilite_dir[c][i][j];
+					hilite_dir[8+c][i-2][j]  += hilite_dir[c][i][j];
+					hilite_dir[12+c][i+2][j] += hilite_dir[c][i][j];
+				}
+			}
+		}
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+	for (int j=hfw-2; j>0; j--)  
+		for (int i=2; i<hfh-2; i++) {
+			//from right 
+			if (hilite[3][i][j]>0.01) {
+				for (int c=0; c<4; c++) {
+					hilite_dir[4+c][i][j] = hilite[c][i][j]/hilite[3][i][j];
+				}
+			} else {
+				for (int c=0; c<4; c++) {
+					hilite_dir[4+c][i][j] = 0.1*((hilite_dir[4+c][(i-2)][(j+1)]+hilite_dir[4+c][(i-1)][(j+1)]+hilite_dir[4+c][(i)][(j+1)]+hilite_dir[4+c][(i+1)][(j+1)]+hilite_dir[4+c][(i+2)][(j+1)])/ \
+														  (hilite_dir[4+3][(i-2)][(j+1)]+hilite_dir[4+3][(i-1)][(j+1)]+hilite_dir[4+3][(i)][(j+1)]+hilite_dir[4+3][(i+1)][(j+1)]+hilite_dir[4+3][(i+2)][(j+1)]+0.00001));
+					hilite_dir[8+c][i-2][j] += hilite_dir[4+c][i][j];
+					hilite_dir[12+c][i+2][j] += hilite_dir[4+c][i][j];
+				}
+			}
+		}
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+	for (int i=1; i<hfh-1; i++) 
+		for (int j=2; j<hfw-2; j++) {
+			//if (i%100==0 && j%100==0)
+			//	printf("row=%d  col=%d  r=%f  g=%f  b=%f hilite=%f  \n",i,j,hilite[0][i][j],hilite[1][i][j],hilite[2][i][j],hilite[3][i][j]);
+			
+			//from top 
+			if (hilite[3][i][j]>0.01) {
+				for (int c=0; c<4; c++) {
+					hilite_dir[8+c][i][j] = hilite[c][i][j]/hilite[3][i][j];
+				}
+			} else {
+				for (int c=0; c<4; c++) {
+					hilite_dir[8+c][i][j] = 0.1*((hilite_dir[8+c][i-1][j-2]+hilite_dir[8+c][i-1][j-1]+hilite_dir[8+c][i-1][j]+hilite_dir[8+c][i-1][j+1]+hilite_dir[8+c][i-1][j+2])/ \
+											  (hilite_dir[8+3][i-1][j-2]+hilite_dir[8+3][i-1][j-1]+hilite_dir[8+3][i-1][j]+hilite_dir[8+3][i-1][j+1]+hilite_dir[8+3][i-1][j+2]+0.00001));
+					hilite_dir[12+c][i+1][j] += hilite_dir[8+c][i][j];
+				}
+			}
+		}
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+	for (int i=hfh-2; i>0; i--) 
+		for (int j=2; j<hfw-2; j++) {
+			//from bottom 
+			if (hilite[3][i][j]>0.01) {
+				for (int c=0; c<4; c++) {
+					hilite_dir[12+c][i][j] = hilite[c][i][j]/hilite[3][i][j];
+				}
+			} else {
+				for (int c=0; c<4; c++) {
+					hilite_dir[12+c][i][j] = 0.1*((hilite_dir[12+c][(i+1)][(j-2)]+hilite_dir[12+c][(i+1)][(j-1)]+hilite_dir[12+c][(i+1)][(j)]+hilite_dir[12+c][(i+1)][(j+1)]+hilite_dir[12+c][(i+1)][(j+2)])/ \
+														  (hilite_dir[12+3][(i+1)][(j-2)]+hilite_dir[12+3][(i+1)][(j-1)]+hilite_dir[12+3][(i+1)][(j)]+hilite_dir[12+3][(i+1)][(j+1)]+hilite_dir[12+3][(i+1)][(j+2)]+0.00001));
+				}
+			}
+	
+		}
+	
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+	//fill in edges
+	for (int dir=0; dir<numdirs; dir++) {
+		for (int i=1; i<hfh-1; i++) 
+			for (int c=0; c<4; c++) {
+				hilite_dir[dir*4+c][i][0] = hilite_dir[dir*4+c][i][1];
+				hilite_dir[dir*4+c][i][hfw-1] = hilite_dir[dir*4+c][i][hfw-2];
+			}
+		
+		for (int j=1; j<hfw-1; j++) 
+			for (int c=0; c<4; c++) {
+				hilite_dir[dir*4+c][0][j] = hilite_dir[dir*4+c][1][j];
+				hilite_dir[dir*4+c][hfh-1][j] = hilite_dir[dir*4+c][hfh-2][j];
+			}
+		
+		for (int c=0; c<4; c++) {
+			hilite_dir[dir*4+c][0][0] = hilite_dir[dir*4+c][1][0] = hilite_dir[dir*4+c][0][1] = hilite_dir[dir*4+c][1][1] = hilite_dir[dir*4+c][2][2];
+			hilite_dir[dir*4+c][0][hfw-1] = hilite_dir[dir*4+c][1][hfw-1] = hilite_dir[dir*4+c][0][hfw-2] = hilite_dir[dir*4+c][1][hfw-2] = hilite_dir[dir*4+c][2][hfw-3];
+			hilite_dir[dir*4+c][hfh-1][0] = hilite_dir[dir*4+c][hfh-2][0] = hilite_dir[dir*4+c][hfh-1][1] = hilite_dir[dir*4+c][hfh-2][1] = hilite_dir[dir*4+c][hfh-3][2];
+			hilite_dir[dir*4+c][hfh-1][hfw-1] = hilite_dir[dir*4+c][hfh-2][hfw-1] = hilite_dir[dir*4+c][hfh-1][hfw-2] = hilite_dir[dir*4+c][hfh-2][hfw-2] = hilite_dir[dir*4+c][hfh-3][hfw-3];
+		}
+	}
+	
+	if(plistener) plistener->setProgress(0.75);
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+	/*for (int m=0; m<4*numdirs; m++) {
+	 boxblur2(hilite_dir[m],hilite_dir[m],hfh,hfw,4);
+	 }*/
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	//now we have highlight data extended in various directions
+	//next step is to build back local data by averaging
+	
+
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	// now reconstruct clipped channels using color ratios
+	
+	//const float Yclip = (0.299*max[0] + 0.587*max[0] + 0.114*max[0]);
+	const float Yclip = 0.3333*(max[0] + max[1] + max[2]);
+	
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+	for (int i=0; i<height; i++) {
+		int i1 = MIN((i-(i%pitch))/pitch,hfh-1);
+		for (int j=0; j<width; j++) {
+			int j1 = MIN((j-(j%pitch))/pitch,hfw-1);
+						
+			float pixel[3]={red[i][j],green[i][j],blue[i][j]};
+			if (pixel[0]<max[0] && pixel[1]<max[1] && pixel[2]<max[2]) continue;//pixel not clipped
+			//if (pixel[0]<fixthresh[0] && pixel[1]<fixthresh[1] && pixel[2]<fixthresh[2]) continue;//pixel not clipped
+			
+			//there are clipped highlights
+			//first, determine weighted average of unclipped extensions (weighting is by 'hue' proximity)
+			float dirwt, totwt=0, factor, Y, I, Q;
+			float clipfix[3]={0,0,0};
+			for (int dir=0; dir<numdirs; dir++) {
+				float Yhi = 0.001+(hilite_dir[dir*4+0][i1][j1] + hilite_dir[dir*4+1][i1][j1] + hilite_dir[dir*4+2][i1][j1]);
+				float Y   = 0.001+(hfsize[0][i1][j1] + hfsize[1][i1][j1] + hfsize[2][i1][j1]);
+				if (hilite_dir[dir*4+0][i1][j1]+hilite_dir[dir*4+1][i1][j1]+hilite_dir[dir*4+2][i1][j1]>0.5) {
+					dirwt = invfn[65535*(SQR(hfsize[0][i1][j1]/Y-hilite_dir[dir*4+0][i1][j1]/Yhi) + \
+										 SQR(hfsize[1][i1][j1]/Y-hilite_dir[dir*4+1][i1][j1]/Yhi) + \
+										 SQR(hfsize[2][i1][j1]/Y-hilite_dir[dir*4+2][i1][j1]/Yhi))];
+					totwt += dirwt;
+					clipfix[0] += dirwt*hilite_dir[dir*4+0][i1][j1]/(hilite_dir[dir*4+3][i1][j1]+0.00001);
+					clipfix[1] += dirwt*hilite_dir[dir*4+1][i1][j1]/(hilite_dir[dir*4+3][i1][j1]+0.00001);
+					clipfix[2] += dirwt*hilite_dir[dir*4+2][i1][j1]/(hilite_dir[dir*4+3][i1][j1]+0.00001);
+				}
+			}
+			clipfix[0] /= totwt;
+			clipfix[1] /= totwt;
+			clipfix[2] /= totwt;
+					
+			//now correct clipped channels
+			if (pixel[0]>max[0] && pixel[1]>max[1] && pixel[2]>max[2]) {
+				//all channels clipped
+				float Y = (0.299*clipfix[0] + 0.587*clipfix[1] + 0.114*clipfix[2]);
+				//float Y = (clipfix[0] + clipfix[1] + clipfix[2]);
+				factor = whitept/Y;
+				red[i][j]   = clipfix[0]*factor;
+				green[i][j] = clipfix[1]*factor;
+				blue[i][j]  = clipfix[2]*factor;
+			} else {//some channels clipped
+				int notclipped[3] = {pixel[0]<max[0] ? 1 : 0, pixel[1]<max[1] ? 1 : 0, pixel[2]<max[2] ? 1 : 0};
+				
+				if (notclipped[0]==0) {//red clipped
+					red[i][j]  = (clipfix[0]*MAX(1,((notclipped[1]*pixel[1] + notclipped[2]*pixel[2])/ \
+																(notclipped[1]*clipfix[1] + notclipped[2]*clipfix[2]))));
+				}
+				if (notclipped[1]==0) {//green clipped
+					green[i][j] = (clipfix[1]*MAX(1,((notclipped[2]*pixel[2] + notclipped[0]*pixel[0])/ \
+																   (notclipped[2]*clipfix[2] + notclipped[0]*clipfix[0]))));
+				}
+				if (notclipped[2]==0) {//blue clipped
+					blue[i][j]  = (clipfix[2]*MAX(1,((notclipped[0]*pixel[0] + notclipped[1]*pixel[1])/ \
+																 (notclipped[0]*clipfix[0] + notclipped[1]*clipfix[1]))));
+				}
+			}
+			
+			/*if (hilite[3][i1][j1]>0.01) {
+				red[i][j]	= (red[i][j]  + hilite[0][i1][j1])/(1+hilite[3][i1][j1]);
+				green[i][j] = (green[i][j]+ hilite[1][i1][j1])/(1+hilite[3][i1][j1]);
+				blue[i][j]	= (blue[i][j] + hilite[2][i1][j1])/(1+hilite[3][i1][j1]);
+			}*/
+			
+			Y = (0.299 * red[i][j] + 0.587 * green[i][j] + 0.114 * blue[i][j]);
+			
+			if (Y>whitept) {
+				factor = whitept/Y;
+				
+				/*I = (0.596 * red[i][j] - 0.275 * green[i][j] - 0.321 * blue[i][j]);
+				Q = (0.212 * red[i][j] - 0.523 * green[i][j] + 0.311 * blue[i][j]);
+				
+				Y *= factor;
+				I *= factor;//MAX(0,MIN(1,(whitept-Y)/(whitept-clippt)));
+				Q *= factor;//MAX(0,MIN(1,(whitept-Y)/(whitept-clippt)));
+				
+				red[i][j]   = Y + 0.956*I + 0.621*Q;
+				green[i][j] = Y - 0.272*I - 0.647*Q;
+				blue[i][j]  = Y - 1.105*I + 1.702*Q;*/
+				
+				red[i][j]   *= factor;
+				green[i][j] *= factor;
+				blue[i][j]  *= factor;
+			}
+		}
+	}
+	
+	if(plistener) plistener->setProgress(1.00);
+
+	
+	// diagnostic output
+	/*for (int i=0; i<height; i++) {
+		int i1 = (i-(i%pitch))/pitch;
+		for (int j=0; j<width; j++) {
+			int j1 = (j-(j%pitch))/pitch;
+			
+			//red[i][j]  =hfsize[0][i1][j1];
+			//green[i][j]=hfsize[1][i1][j1];
+			//blue[i][j] =hfsize[2][i1][j1];
+			
+			//red[i][j]=  hilite[0][i1][j1]/(hilite[3][i1][j1]+0.001);
+			//green[i][j]=hilite[1][i1][j1]/(hilite[3][i1][j1]+0.001);
+			//blue[i][j]= hilite[2][i1][j1]/(hilite[3][i1][j1]+0.001);
+			
+			red[i][j]=  hilite_dir[8+0][i1][j1]/hilite_dir[8+3][i1][j1];
+			green[i][j]=hilite_dir[8+1][i1][j1]/hilite_dir[8+3][i1][j1];
+			blue[i][j]= hilite_dir[8+2][i1][j1]/hilite_dir[8+3][i1][j1];
+		
+			//red[i][j]=  clipfix[0][i1][j1];
+			//green[i][j]=clipfix[1][i1][j1];
+			//blue[i][j]= clipfix[2][i1][j1];
+			
+		}
+	}*/
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+}// end of HLReconstruction
+	
+/*
+	void RawImageSource::halfsize()
+	{
+		int ex,ey;
+		//determine GRBG coset; (ey,ex) is the offset of the R subarray
+		if (FC(0,0)==1) {//first pixel is G
+			if (FC(0,1)==0) {ey=0; ex=1;} else {ey=1; ex=0;}
+		} else {//first pixel is R or B
+			if (FC(0,0)==0) {ey=0; ex=0;} else {ey=1; ex=1;}
+		}
+		
+		for (int i=0; i<(H-(H&1)); i++) {
+			for (int j=0; j<(W-(W&1)); j++){
+				red[i][j] = rawData[i+ey-(i&1)][j+ex-(j&1)];
+				green[i][j] = (rawData[i+(1-ey)-(i&1)][j+ex-(j&1)]+rawData[i+ey-(i&1)][j+(1-ex)-(j&1)])/2;
+				blue[i][j] = rawData[i+(1-ey)-(i&1)][j+(1-ex)-(j&1)];
+			}
+			if (W&1) {
+				red[i][W-1]=red[i][W-2];
+				green[i][W-1]=green[i][W-2];
+				blue[i][W-1]=blue[i][W-2];
+			}
+		}
+		if (H&1) {
+			for (int j=0; j<W; j++){
+				red[H-1][j] = red[H-2][j];
+				green[H-1][j] = green[H-2][j];
+				blue[H-1][j] = blue[H-2][j];
+			}
+		}
+	}
+*/
+	
+//}
+
diff --git a/rtengine/imagesource.h b/rtengine/imagesource.h
index 05a65bdcb..17d47b225 100644
--- a/rtengine/imagesource.h
+++ b/rtengine/imagesource.h
@@ -68,8 +68,10 @@ class ImageSource : public InitialImage {
 
         virtual ~ImageSource            () {}
         virtual int         load        (Glib::ustring fname, bool batch = false) =0;
-        virtual void        preprocess  (const RAWParams &raw){};
-        virtual void        demosaic    (const RAWParams &raw){};
+        virtual void        preprocess  (const RAWParams &raw, HRecParams hrp){};
+        virtual void        demosaic    (const RAWParams &raw, HRecParams hrp){};
+	virtual void		HLRecovery_inpaint (float** red, float** green, float** blue){};
+
         virtual void        getImage    (ColorTemp ctemp, int tran, Imagefloat* image, PreviewProps pp, HRecParams hlp, ColorManagementParams cmp, RAWParams raw) {}
         virtual ColorTemp   getWB       () =0;
         virtual ColorTemp   getAutoWB   () =0;
diff --git a/rtengine/improccoordinator.cc b/rtengine/improccoordinator.cc
index 3de8a472e..bb65e0a71 100644
--- a/rtengine/improccoordinator.cc
+++ b/rtengine/improccoordinator.cc
@@ -112,13 +112,13 @@ void ImProcCoordinator::updatePreviewImage (int todo, Crop* cropCall) {
 
     progress ("Applying white balance, color correction & sRBG conversion...",100*readyphase/numofphases);
     if ( todo & M_PREPROC) {
-    	imgsrc->preprocess( rp );
+    	imgsrc->preprocess( rp, params.hlrecovery );
         imgsrc->getRAWHistogram( histRedRaw, histGreenRaw, histBlueRaw );
     }
 
     if( todo & M_RAW){
     	fineDetailsProcessed = highDetailNeeded;
-    	imgsrc->demosaic( rp );
+    	imgsrc->demosaic( rp, params.hlrecovery );
     }
 
     if (todo & M_INIT) {
@@ -268,7 +268,7 @@ void ImProcCoordinator::updatePreviewImage (int todo, Crop* cropCall) {
             progress ("Pyramid equalizer...",100*readyphase/numofphases);
             ipf.dirpyrequalizer (nprevl);
 
-            // Superset by pyreq
+            // Superseded by dirpyreq
             //progress ("Wavelet...",100*readyphase/numofphases);
             //ipf.waveletEqualizer (nprevl, true, true);
         }
@@ -523,8 +523,8 @@ void ImProcCoordinator::saveInputICCReference (const Glib::ustring& fname) {
 	ppar.icm.input = "(none)";
 	Imagefloat* im = new Imagefloat (fW, fH);
 	Image16* im16 = new Image16 (fW, fH);
-	imgsrc->preprocess( ppar.raw );
-	imgsrc->demosaic(ppar.raw );
+	imgsrc->preprocess( ppar.raw, ppar.hlrecovery );
+	imgsrc->demosaic(ppar.raw, ppar.hlrecovery );
 	//imgsrc->getImage (imgsrc->getWB(), 0, im, pp, ppar.hlrecovery, ppar.icm, ppar.raw);
 	ColorTemp currWB = ColorTemp (params.wb.temperature, params.wb.green);
 	if (params.wb.method=="Camera")
diff --git a/rtengine/improcfun.cc b/rtengine/improcfun.cc
index 0baa62c56..26810d5f8 100644
--- a/rtengine/improcfun.cc
+++ b/rtengine/improcfun.cc
@@ -266,8 +266,8 @@ void ImProcFunctions::rgbProc (Imagefloat* working, LabImage* lab, LUTf & hltone
 	if (vCurveEnabled) vCurve = new FlatCurve(params->hsvequalizer.vcurve);
 	
 	const float exp_scale = pow (2.0, params->toneCurve.expcomp);
-	const float comp = (params->toneCurve.expcomp + 1.0)*params->toneCurve.hlcompr/100.0;
-	const float shoulder = ((65536.0/exp_scale)*(params->toneCurve.hlcomprthresh/200.0))+0.1;
+	const float comp = (MAX(0,params->toneCurve.expcomp) + 1.0)*params->toneCurve.hlcompr/100.0;
+	const float shoulder = ((65536.0/MAX(1,exp_scale))*(params->toneCurve.hlcomprthresh/200.0))+0.1;
 	const float hlrange = 65536.0-shoulder;
 	
 	
@@ -591,7 +591,7 @@ void ImProcFunctions::colorCurve (LabImage* lold, LabImage* lnew) {
 		
 		if (params->defringe.enabled && lab->W>=8 && lab->H>=8)
 			
-			PF_correct_RT(lab, lab, params->defringe.radius, params->defringe.threshold, false /*edges only*/ );
+			PF_correct_RT(lab, lab, params->defringe.radius, params->defringe.threshold);
 	}
 	
 	void ImProcFunctions::dirpyrdenoise (LabImage* lab) {
diff --git a/rtengine/improcfun.h b/rtengine/improcfun.h
index a5c5969d7..44b724247 100644
--- a/rtengine/improcfun.h
+++ b/rtengine/improcfun.h
@@ -110,7 +110,7 @@ class ImProcFunctions {
 		void idirpyr_eq_channel  (float ** data_coarse, float ** data_fine, float ** buffer, int width, int height, int level, const double * mult );
 
 		void defringe         (LabImage* lab);
-		void PF_correct_RT    (LabImage * src, LabImage * dst, double radius, int thresh, bool edges);
+		void PF_correct_RT    (LabImage * src, LabImage * dst, double radius, int thresh);
 
 		Image8*     lab2rgb   (LabImage* lab, int cx, int cy, int cw, int ch, Glib::ustring profile);
 		Image16*    lab2rgb16b (LabImage* lab, int cx, int cy, int cw, int ch, Glib::ustring profile, Glib::ustring profi, Glib::ustring gam, bool freegamma, double gampos, double slpos);// for gamma output		
diff --git a/rtengine/rawimagesource.cc b/rtengine/rawimagesource.cc
index 4a350d8db..da0e9fea0 100644
--- a/rtengine/rawimagesource.cc
+++ b/rtengine/rawimagesource.cc
@@ -82,7 +82,7 @@ RawImageSource::RawImageSource ()
     hrmap[0] = NULL;
     hrmap[1] = NULL;
     hrmap[2] = NULL;
-	needhr = NULL;
+	//needhr = NULL;
     //hpmap = NULL;
 	camProfile = NULL;
 	embProfile = NULL;
@@ -113,8 +113,8 @@ RawImageSource::~RawImageSource () {
         freeArray<float>(hrmap[1], dh);
         freeArray<float>(hrmap[2], dh);
     }
-    if (needhr)
-        freeArray<char>(needhr, H);
+    //if (needhr)
+    //    freeArray<char>(needhr, H);
     //if (hpmap)
     //    freeArray<char>(hpmap, H);
     if (camProfile)
@@ -253,8 +253,8 @@ void RawImageSource::getImage (ColorTemp ctemp, int tran, Imagefloat* image, Pre
 
 
 
-	if (hrp.enabled==true && hrp.method=="Color" && hrmap[0]==NULL) 
-        updateHLRecoveryMap_ColorPropagation ();
+	//if (hrp.enabled==true && hrp.method=="Color" && hrmap[0]==NULL) 
+    //    updateHLRecoveryMap_ColorPropagation ();
 
     // compute image area to render in order to provide the requested part of the image
     int sx1, sy1, imwidth, imheight, fw;
@@ -278,6 +278,7 @@ void RawImageSource::getImage (ColorTemp ctemp, int tran, Imagefloat* image, Pre
         imwidth = maximwidth;
     if (!fuji && imheight>maximheight)
         imheight = maximheight;
+	
     int maxx=this->W,maxy=this->H,skip=pp.skip;
 
     //if (sx1+skip*imwidth>maxx) imwidth --; // very hard to fix this situation without an 'if' in the loop.
@@ -351,8 +352,9 @@ void RawImageSource::getImage (ColorTemp ctemp, int tran, Imagefloat* image, Pre
 				
             }
         }
-                
-        if (hrp.enabled)
+		
+		//process all highlight recovery other than "Color"
+        if (hrp.enabled && hrp.method!="Color")
 			hlRecovery (hrp.method, line_red, line_grn, line_blue, i, sx1, imwidth, skip, raw);
 
         transLine (line_red, line_grn, line_blue, ix, image, tran, imwidth, imheight, fw);
@@ -407,6 +409,7 @@ void RawImageSource::getImage (ColorTemp ctemp, int tran, Imagefloat* image, Pre
         }
     }
 
+
 	
     // Flip if needed
     if (tran & TR_HFLIP)
@@ -945,7 +948,7 @@ int RawImageSource::load (Glib::ustring fname, bool batch) {
 	
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-void RawImageSource::preprocess  (const RAWParams &raw)
+void RawImageSource::preprocess  (const RAWParams &raw, HRecParams hrp)
 {
 	MyTime t1,t2;
 	t1.set();
@@ -1092,7 +1095,7 @@ void RawImageSource::preprocess  (const RAWParams &raw)
 
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 	
-void RawImageSource::demosaic(const RAWParams &raw)
+void RawImageSource::demosaic(const RAWParams &raw, HRecParams hrp )
 {
     if (ri->isBayer()) {
     	MyTime t1,t2;
@@ -1122,6 +1125,10 @@ void RawImageSource::demosaic(const RAWParams &raw)
         plistener->setProgressStr ("Ready.");
         plistener->setProgress (1.0);
     }
+	
+	//color propagation highlight recovery 
+	if (hrp.enabled && hrp.method=="Color")
+		HLRecovery_inpaint (red,green,blue);
 
 }
 	
@@ -1280,10 +1287,12 @@ void RawImageSource::copyOriginalPixels(const RAWParams &raw, RawImage *src, Raw
 	}
 	
 
-/* Scale original pixels into the range 0 65535 using black offsets and multipliers */
+// Scale original pixels into the range 0 65535 using black offsets and multipliers 
 void RawImageSource::scaleColors(int winx,int winy,int winw,int winh, const RAWParams &raw)
 {
 float black_lev[4];//black level
+	
+	chmax[0]=chmax[1]=chmax[2]=chmax[3]=0;//channel maxima
 
 //adjust black level  (eg Canon)
 // cblack Bayer
@@ -1341,6 +1350,7 @@ black_lev[3]= raw.blackzero;
 				val*=scale_mul[2];}	
 				
 				rawData[row][col] = (val);
+				chmax[c] = MAX(chmax[c],val);
 			}
 		}
 	}else{
@@ -1352,21 +1362,26 @@ black_lev[3]= raw.blackzero;
 					val -= cblack[0];
 					val *= scale_mul[0];
 					rawData[row][3*col+0] = (val);
+					chmax[0] = MAX(chmax[0],val);
 				}
 				val = rawData[row][3*col+1];
 				if (val){
 					val -= cblack[1];
 					val *= scale_mul[1];
 					rawData[row][3*col+1] = (val);
+					chmax[1] = MAX(chmax[1],val);
 				}
 				val = rawData[row][3*col+2];
 				if (val){
 					val -= cblack[2];
 					val *= scale_mul[2];
 					rawData[row][3*col+2] = (val);
+					chmax[2] = MAX(chmax[2],val);
 				}
 			}
 		}
+		
+		chmax[3]=chmax[1];
 	}
 
 }
@@ -1405,11 +1420,11 @@ void RawImageSource::processFalseColorCorrectionThread  (Imagefloat* im, int row
  
   int W = im->width;
 
-  float** rbconv_Y = allocArray<float>(W,3);
-  float** rbconv_I = allocArray<float>(W,3);
-  float** rbconv_Q = allocArray<float>(W,3);
-  float** rbout_I = allocArray<float>(W,3);
-  float** rbout_Q = allocArray<float>(W,3);
+  array2D<float> rbconv_Y (W,3);
+  array2D<float> rbconv_I (W,3);
+  array2D<float> rbconv_Q (W,3);
+  array2D<float> rbout_I (W,3);
+  array2D<float> rbout_Q (W,3);
 
   float* row_I = new float[W];
   float* row_Q = new float[W];
@@ -1511,11 +1526,6 @@ void RawImageSource::processFalseColorCorrectionThread  (Imagefloat* im, int row
     row_Q[W-1] = rbout_Q[cx][W-1];
     convert_row_to_RGB (im->r[row_to-1], im->g[row_to-1], im->b[row_to-1], rbconv_Y[cx], row_I, row_Q, W);
 
-  freeArray<float>(rbconv_Y, 3);
-  freeArray<float>(rbconv_I, 3);
-  freeArray<float>(rbconv_Q, 3);
-  freeArray<float>(rbout_I, 3);
-  freeArray<float>(rbout_Q, 3);
   delete [] row_I;
   delete [] row_Q;
   delete [] pre1_I;
@@ -1951,8 +1961,8 @@ void RawImageSource::hlRecovery (std::string method, float* red, float* green, f
         HLRecovery_Luminance (red, green, blue, red, green, blue, width, 65535.0);
     else if (method=="CIELab blending")
         HLRecovery_CIELab (red, green, blue, red, green, blue, width, 65535.0, xyz_cam, cam_xyz);
-    else if (method=="Color")
-        HLRecovery_ColorPropagation (red, green, blue, i, sx1, width, skip);
+    /*else if (method=="Color")
+        HLRecovery_ColorPropagation (red, green, blue, i, sx1, width, skip);*/
 	else if (method=="Blend")	// derived from Dcraw
 			{	float pre_mul[4];
 				for(int c=0;c<4;c++) pre_mul[c]=ri->get_pre_mul(c);
@@ -2324,12 +2334,14 @@ void RawImageSource::inverse33 (double (*rgb_cam)[3], double (*cam_rgb)[3]) {
 //#include "demosaic_algos.cc"
 	
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-//Emil's code for AMaZE
+//Emil's code 
 #include "fast_demo.cc"//fast demosaic	
 #include "amaze_demosaic_RT.cc"//AMaZE demosaic	
 #include "CA_correct_RT.cc"//Emil's CA auto correction
 #include "cfa_linedn_RT.cc"//Emil's line denoise
 #include "green_equil_RT.cc"//Emil's green channel equilibration
+#include "hilite_recon.cc"//Emil's highlight reconstruction
+
 #include "expo_before_b.cc"//Jacques's exposure before interpolation
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
diff --git a/rtengine/rawimagesource.h b/rtengine/rawimagesource.h
index 52a04a71d..ad24ba0de 100644
--- a/rtengine/rawimagesource.h
+++ b/rtengine/rawimagesource.h
@@ -21,6 +21,7 @@
 
 #include <imagesource.h>
 #include <lcms2.h>
+#include <array2D.h>
 
 #define HR_SCALE 2
 
@@ -82,10 +83,10 @@ class RawImageSource : public ImageSource {
         //char** hpmap;
         float** hrmap[3];   // for color propagation
         char** needhr;      // for color propagation
-        int max[3], g;
+		int max[3];
+		float chmax[4];
         double initialGain; // initial gain calculated after scale_colors
         double defGain;
-        //int blcode[16][16][32];  // Looks like it's an unused variable...
         bool full;
         cmsHPROFILE camProfile;
         cmsHPROFILE embProfile;
@@ -128,8 +129,8 @@ class RawImageSource : public ImageSource {
         ~RawImageSource ();
 
         int         load        (Glib::ustring fname, bool batch = false);
-        void        preprocess  (const RAWParams &raw);
-        void        demosaic    (const RAWParams &raw);
+        void        preprocess  (const RAWParams &raw, HRecParams hrp);
+        void        demosaic    (const RAWParams &raw, HRecParams hrp);
         void        copyOriginalPixels(const RAWParams &raw, RawImage *ri, RawImage *riDark, RawImage *riFlatFile  );
         void        cfaboxblur  (RawImage *riFlatFile, float* cfablur, int boxH, int boxW );
         void        scaleColors (int winx,int winy,int winw,int winh, const RAWParams &raw);// raw for cblack
@@ -155,6 +156,14 @@ class RawImageSource : public ImageSource {
         static void colorSpaceConversion (Imagefloat* im, ColorManagementParams cmp, cmsHPROFILE embedded, cmsHPROFILE camprofile, double cam[3][3], double& defgain);
         static void inverse33 (double (*coeff)[3], double (*icoeff)[3]);
 
+	void boxblur2(float** src, float** dst, int H, int W, int box );
+	void boxblur_resamp(float **src, float **dst, float & max, int H, int W, int box, int samp ); 
+
+	//void boxblur_resamp(float **red, float **green, float **blue, int H, int W, float thresh[3], float max[3], \
+						multi_array2D<float,3> & hfsize, multi_array2D<float,3> & hilite, int box );
+	void HLRecovery_inpaint (float** red, float** green, float** blue);
+	//void HLRecovery_inpaint ();
+	
         static void HLRecovery_Luminance (float* rin, float* gin, float* bin, float* rout, float* gout, float* bout, int width, float maxval);
         static void HLRecovery_CIELab (float* rin, float* gin, float* bin, float* rout, float* gout, float* bout, int width, float maxval, double cam[3][3], double icam[3][3]);
 		static void HLRecovery_blend (float* rin, float* gin, float* bin, int width, float maxval, float* pre_mul, const RAWParams &raw);
diff --git a/rtengine/simpleprocess.cc b/rtengine/simpleprocess.cc
index 5d881e84d..df320e616 100644
--- a/rtengine/simpleprocess.cc
+++ b/rtengine/simpleprocess.cc
@@ -98,9 +98,9 @@ IImage16* processImage (ProcessingJob* pjob, int& errorCode, ProgressListener* p
         currWB = imgsrc->getAutoWB ();
 
     PreviewProps pp (0, 0, fw, fh, 1);
-    imgsrc->preprocess( params.raw );
+    imgsrc->preprocess( params.raw, params.hlrecovery );
 	if (pl) pl->setProgress (0.20);
-    imgsrc->demosaic( params.raw );
+    imgsrc->demosaic( params.raw, params.hlrecovery );
 	if (pl) pl->setProgress (0.40);
     Imagefloat* baseImg = new Imagefloat (fw, fh);
     imgsrc->getImage (currWB, tr, baseImg, pp, params.hlrecovery, params.icm, params.raw);
diff --git a/rtgui/tonecurve.cc b/rtgui/tonecurve.cc
index 0ac88d03e..b262bdd61 100644
--- a/rtgui/tonecurve.cc
+++ b/rtgui/tonecurve.cc
@@ -51,7 +51,7 @@ ToneCurve::ToneCurve () : Gtk::VBox(), FoldableToolPanel(this) {
 //----------- Exposure Compensation ------------------------
   expcomp   = Gtk::manage (new Adjuster (M("TP_EXPOSURE_EXPCOMP"), -5, 10, 0.05, 0));
   pack_start (*expcomp);
-  hlcompr = Gtk::manage (new Adjuster (M("TP_EXPOSURE_COMPRHIGHLIGHTS"), 0, 100, 1, 70));
+  hlcompr = Gtk::manage (new Adjuster (M("TP_EXPOSURE_COMPRHIGHLIGHTS"), 0, 500, 1, 70));
   pack_start (*hlcompr);
   hlcomprthresh = Gtk::manage (new Adjuster (M("TP_EXPOSURE_COMPRHIGHLIGHTSTHRESHOLD"), 0, 100, 1, 0));
   pack_start (*hlcomprthresh);