Improvements for highlight reconstruction modes 'Blend' and 'Color propagation'.

rtengine/hilite_recon.cc

@@ -41,6 +41,7 @@
 #define MAX(a,b) ((a)<(b)?(b):(a))
 #define MIN(a,b) ((a)>(b)?(b):(a))
 
+#define FOREACHCOLOR for (int c=0; c < ColorCount; c++)				
 
 //#include "RGBdefringe.cc"
 
@@ -198,11 +199,21 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b
 	
 	static const int numdirs = 4;
 
-	static const float threshpct = 0.25;
+	//%%%%%%%%%%%%%%%%%%%%
-	static const float fixthreshpct = 0.7;
+	//for blend algorithm:
-	static const float maxpct = 0.95;
+	static const float clipthresh = 0.95, blendthresh=1.0;
+	const int ColorCount=3;
+	// Transform matrixes rgb>lab and back
+	static const float trans[2][ColorCount][ColorCount] =
+	{ { { 1,1,1 }, { 1.7320508,-1.7320508,0 }, { -1,-1,2 } },
+		{ { 1,1,1 }, { 1,-1,1 }, { 1,1,-1 } } };
+	static const float itrans[2][ColorCount][ColorCount] =
+	{ { { 1,0.8660254,-0.5 }, { 1,-0.8660254,-0.5 }, { 1,0,1 } },
+		{ { 1,1,1 }, { 1,-1,1 }, { 1,1,-1 } } };
+	//%%%%%%%%%%%%%%%%%%%%
 
-	float max[3], thresh[3], fixthresh[3], norm[3];
+	
+	float max[3], thresh[3], fixthresh[3];
 		
 	//float red1, green1, blue1;//diagnostic
 	float chmaxalt[4]={0,0,0,0};//diagnostic
@@ -229,31 +240,39 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b
 	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 	
+	static const float threshpct = 0.25;//percentage of channel max above which we have highlights
+	static const float fixthreshpct = 0.7;//percentage of channel max above which we estimate reconstructed highlights
+	static const float maxpct = 0.95;//percentage of channel max above which we consider data to be in clipping range
+	static const float satthreshpct = 0.5;//threshold pct. of clipping above which we desaturate the reconstructed estimate
+	
 	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;
+		max[c] = chmax[c]*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 medpt   = max[0]+max[1]+max[2]-whitept-clippt;
+	float maxave=(max[0]+max[1]+max[2])/3;
 	
+	float clip[3];
+	FOREACHCOLOR clip[c]=MIN(maxave,chmax[c]);
 		
-	float sat = ri->get_white();
+	// Determine the maximum level (clip) of all channels
-	float black = ri->get_black();
+	const float fixpt = fixthreshpct*clippt;
-	sat -= black;
+	const float desatpt = satthreshpct*maxave+(1-satthreshpct)*clippt;
-	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,3> channelblur(width,height,ARRAY2D_CLEAR_DATA);
 	
 	multi_array2D<float,5> hilite_full(width,height,ARRAY2D_CLEAR_DATA);
 		
 	// blur RGB channels
-	boxblur2(red  ,channelblur[0],height,width,4);
+	boxblur2(red  ,channelblur[0],height,width,3);
-	boxblur2(green,channelblur[1],height,width,4);
+	boxblur2(green,channelblur[1],height,width,3);
-	boxblur2(blue ,channelblur[2],height,width,4);
+	boxblur2(blue ,channelblur[2],height,width,3);
 	
 	float hipass_sum=0, hipass_norm=0.01;
 	
@@ -469,6 +488,7 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b
 	
 	const float Yclip = (0.299*max[0] + 0.587*max[1] + 0.114*max[2]);
 	//const float Yclip = 0.3333*(max[0] + max[1] + max[2]);
+	//float sumwt=0, counts=0;
 	
 #ifdef _OPENMP
 #pragma omp parallel for
@@ -482,19 +502,95 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b
 			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
 			
-			float pixref[3]={MIN(Yclip,hfsize[0][i1][j1]),MIN(Yclip,hfsize[1][i1][j1]),MIN(Yclip,hfsize[2][i1][j1])};
+			//%%%%%%%%%%%%%%%%%%%%%%%
+			//estimate recovered values using modified HLRecovery_blend algorithm
+			
+			float rgb[ColorCount], rgb_blend[ColorCount], cam[2][ColorCount], lab[2][ColorCount], sum[2], chratio;
+			float L,C,H,Lfrac;
+			
+			// Copy input pixel to rgb so it's easier to access in loops
+			rgb[0] = pixel[0]; rgb[1] = pixel[1]; rgb[2] = pixel[2];
+			
+			// Initialize cam with raw input [0] and potentially clipped input [1]
+			FOREACHCOLOR {
+				cam[0][c] = rgb[c];
+				cam[1][c] = MIN(cam[0][c],clippt);
+			}
+			
+			// Calculate the lightness correction ratio (chratio)
+			for (int i=0; i<2; i++) {
+				FOREACHCOLOR {
+					lab[i][c]=0;
+					for (int j=0; j < ColorCount; j++)
+						lab[i][c] += trans[ColorCount-3][c][j] * cam[i][j];
+				}
+				
+				sum[i]=0;
+				for (int c=1; c < ColorCount; c++)
+					sum[i] += SQR(lab[i][c]);
+			}
+			chratio = sqrt(sum[1]/sum[0]);
+			
+			
+			// Apply ratio to lightness in lab space
+			for (int c=1; c < ColorCount; c++) 
+				lab[0][c] *= chratio;
+			
+			// Transform back from lab to RGB
+			FOREACHCOLOR {
+				cam[0][c]=0;
+				for (int j=0; j < ColorCount; j++) {
+					cam[0][c] += itrans[ColorCount-3][c][j] * lab[0][j];
+				}
+			}
+			FOREACHCOLOR rgb[c] = cam[0][c] / ColorCount;
+			
+			// Copy converted pixel back			
+			if (pixel[0] > fixpt) {
+				float rfrac = SQR((MIN(clip[0],pixel[0])-fixpt)/(clip[0]-fixpt));
+				pixel[0]= MIN(maxave,rfrac*rgb[0]+(1-rfrac)*pixel[0]);
+			}
+			if (pixel[1] > fixpt) {
+				float gfrac = SQR((MIN(clip[1],pixel[1])-fixpt)/(clip[1]-fixpt));
+				pixel[1]= MIN(maxave,gfrac*rgb[1]+(1-gfrac)*pixel[1]); 
+			}
+			if (pixel[2] > fixpt) {
+				float bfrac = SQR((MIN(clip[2],pixel[2])-fixpt)/(clip[2]-fixpt));
+				pixel[2]= MIN(maxave,bfrac*rgb[2]+(1-bfrac)*pixel[2]);
+			}
+			
+			if ((L=(pixel[0]+pixel[1]+pixel[2])/3) > desatpt) {
+				Lfrac = MAX(0,(maxave-L)/(maxave-desatpt));
+				C = Lfrac * 1.732050808 * (pixel[0] - pixel[1]);
+				H = Lfrac * (2 * pixel[2] - pixel[0] - pixel[1]);
+				pixel[0] = L - H / 6.0 + C / 3.464101615;
+				pixel[1] = L - H / 6.0 - C / 3.464101615;
+				pixel[2] = L + H / 3.0;
+			}
+			
+			// Copy converted pixel back
+			/*float rfrac = SQR((pixel[0]-clippt)/(chmax[0]-clippt));
+			float gfrac = SQR((pixel[1]-clippt)/(chmax[1]-clippt));
+			float bfrac = SQR((pixel[2]-clippt)/(chmax[2]-clippt));
+			if (pixel[0] > blendthresh*clippt) rgb_blend[0]= rfrac*rgb[0]+(1-rfrac)*pixel[0]; 
+			if (pixel[1] > blendthresh*clippt) rgb_blend[1]= gfrac*rgb[1]+(1-gfrac)*pixel[1]; 
+			if (pixel[2] > blendthresh*clippt) rgb_blend[2]= bfrac*rgb[2]+(1-bfrac)*pixel[2];*/
+			
+			//end of HLRecovery_blend estimation
+			//%%%%%%%%%%%%%%%%%%%%%%%
 						
 			//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 dirwt, factor, Y, I, Q;
-			float clipfix[3]={0,0,0};
+			float totwt=0;//0.0003;
+			float clipfix[3]={0,0,0};//={totwt*rgb_blend[0],totwt*rgb_blend[1],totwt*rgb_blend[2]};
 			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+(pixref[0]+pixref[1]+pixref[2]);
+				float Y   = 0.001+(rgb_blend[0]+rgb_blend[1]+rgb_blend[2]);
 				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(pixref[0]/Y-hilite_dir[dir*4+0][i1][j1]/Yhi) + \
+					dirwt = invfn[65535*(SQR(rgb_blend[0]/Y-hilite_dir[dir*4+0][i1][j1]/Yhi) + \
-										 SQR(pixref[1]/Y-hilite_dir[dir*4+1][i1][j1]/Yhi) + \
+										 SQR(rgb_blend[1]/Y-hilite_dir[dir*4+1][i1][j1]/Yhi) + \
-										 SQR(pixref[2]/Y-hilite_dir[dir*4+2][i1][j1]/Yhi))];
+										 SQR(rgb_blend[2]/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);
@@ -504,6 +600,8 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b
 			clipfix[0] /= totwt;
 			clipfix[1] /= totwt;
 			clipfix[2] /= totwt;
+			//sumwt += totwt;
+			//counts ++;
 					
 			//now correct clipped channels
 			if (pixel[0]>max[0] && pixel[1]>max[1] && pixel[2]>max[2]) {
@@ -561,6 +659,7 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b
 	}
 	
 	if(plistener) plistener->setProgress(1.00);
+	//printf("ave wt=%f\n",sumwt/counts);
 
 	
 	// diagnostic output

rtengine/rawimagesource.cc

@@ -287,6 +287,11 @@ void RawImageSource::getImage (ColorTemp ctemp, int tran, Imagefloat* image, Pre
     gm/=area;
     bm/=area;
 	
+	hlmax[0]=chmax[0]*rm*area;
+	hlmax[1]=chmax[1]*gm*area;
+	hlmax[2]=chmax[2]*bm*area;
+
+
 #ifdef _OPENMP
 #pragma omp parallel
     {
@@ -295,7 +300,7 @@ void RawImageSource::getImage (ColorTemp ctemp, int tran, Imagefloat* image, Pre
     float* line_red  = new float[imwidth];
     float* line_grn  = new float[imwidth];
     float* line_blue = new float[imwidth];
-
+	//printf("clip[0]=%f  clip[1]=%f  clip[2]=%f\n",hlmax[0],hlmax[1],hlmax[2]);
 
 
 #ifdef _OPENMP
@@ -355,7 +360,7 @@ void RawImageSource::getImage (ColorTemp ctemp, int tran, Imagefloat* image, Pre
 
 		//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);
+			hlRecovery (hrp.method, line_red, line_grn, line_blue, i, sx1, imwidth, skip, raw, hlmax);
 
         transLine (line_red, line_grn, line_blue, ix, image, tran, imwidth, imheight, fw);
 		
@@ -1771,10 +1776,9 @@ bool RawImageSource::findInputProfile(Glib::ustring inProfile, cmsHPROFILE embed
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 // derived from Dcraw "blend_highlights()"
 //  very effective to reduce (or remove) the magenta, but with levels of grey !
-void RawImageSource::HLRecovery_blend(float* rin, float* gin, float* bin, int width, float maxval, float* pre_mul, const RAWParams &raw)
+void RawImageSource::HLRecovery_blend(float* rin, float* gin, float* bin, int width, float maxval, float* pre_mul, const RAWParams &raw, float* hlmax)
 	{
 		const int ColorCount=3;
-    int clip=INT_MAX;
 
 		// Transform matrixes rgb>lab and back
 		static const float trans[2][ColorCount][ColorCount] =
@@ -1787,29 +1791,47 @@ void RawImageSource::HLRecovery_blend(float* rin, float* gin, float* bin, int wi
 #define FOREACHCOLOR for (int c=0; c < ColorCount; c++)
 #define SQR(x) ((x)*(x))
 		
+		
+		float minpt=MIN(MIN(hlmax[0],hlmax[1]),hlmax[2]);
+		//float maxpt=MAX(MAX(hlmax[0],hlmax[1]),hlmax[2]);
+		//float medpt=hlmax[0]+hlmax[1]+hlmax[2]-minpt-maxpt;
+		float maxave=(hlmax[0]+hlmax[1]+hlmax[2])/3;
+		const float clipthresh = 0.95;
+		const float fixthresh = 0.5;
+		const float satthresh = 0.5;
+
+		float clip[3];
+		FOREACHCOLOR clip[c]=MIN(maxave,hlmax[c]);
+		
 		// Determine the maximum level (clip) of all channels
-    int i;
+		const float clippt = clipthresh*maxval;
-    FOREACHCOLOR if (clip > (i = (int) maxval * pre_mul[c] * raw.expos)) clip = i;	
+		const float fixpt = fixthresh*minpt;
+		const float desatpt = satthresh*maxave+(1-satthresh)*maxval;
+
 
 #pragma omp parallel for
 		for (int col=0; col<width; col++) {
 			float rgb[ColorCount], cam[2][ColorCount], lab[2][ColorCount], sum[2], chratio;
+			float L,C,H,Lfrac;
 			
 			// Copy input pixel to rgb so it's easier to access in loops
 			rgb[0] = rin[col]; rgb[1] = gin[col]; rgb[2] = bin[col];
 			
 			// If no channel is clipped, do nothing on pixel
 			int c;
-        for (c=0; c<ColorCount; c++) { if (rgb[c] > clip) break; }
+			for (c=0; c<ColorCount; c++) { if (rgb[c] > clippt) break; }
 			if (c == ColorCount) continue;
 			
+			float clipfrac[3];
+			FOREACHCOLOR clipfrac[c] = MIN(1,rgb[c]/maxave);
+			
 			// Initialize cam with raw input [0] and potentially clipped input [1]
 			FOREACHCOLOR {
 				cam[0][c] = rgb[c];
-            cam[1][c] = MIN(cam[0][c],clip);
+				cam[1][c] = MIN(cam[0][c],maxval);
 			}
 			
-        // Calculate the lightness correction ration (chratio)
+			// Calculate the lightness correction ratio (chratio)
 			for (int i=0; i<2; i++) {
 				FOREACHCOLOR {
 					lab[i][c]=0;
@@ -1821,10 +1843,11 @@ void RawImageSource::HLRecovery_blend(float* rin, float* gin, float* bin, int wi
 				for (int c=1; c < ColorCount; c++)
 					sum[i] += SQR(lab[i][c]);
 			}
-        chratio = sqrt(sum[1]/sum[0]);
+			chratio = sqrt(sqrt(sum[1]/sum[0]));
 			
 			// Apply ratio to lightness in lab space
-        for (int c=1; c < ColorCount; c++) lab[0][c] *= chratio;
+			for (int c=1; c < ColorCount; c++) 
+				lab[0][c] *= chratio;
 			
 			// Transform back from lab to RGB
 			FOREACHCOLOR {
@@ -1836,7 +1859,27 @@ void RawImageSource::HLRecovery_blend(float* rin, float* gin, float* bin, int wi
 			FOREACHCOLOR rgb[c] = cam[0][c] / ColorCount;
 
 			// Copy converted pixel back			
-        rin[col]=rgb[0]; gin[col]=rgb[1]; bin[col]=rgb[2];
+			if (rin[col] > fixpt) {
+				float rfrac = SQR((MIN(clip[0],rin[col])-fixpt)/(clip[0]-fixpt));
+				rin[col]= MIN(maxave,rfrac*rgb[0]+(1-rfrac)*rin[col]);
+			}
+			if (gin[col] > fixpt) {
+				float gfrac = SQR((MIN(clip[1],gin[col])-fixpt)/(clip[1]-fixpt));
+				gin[col]= MIN(maxave,gfrac*rgb[1]+(1-gfrac)*gin[col]); 
+			}
+			if (bin[col] > fixpt) {
+				float bfrac = SQR((MIN(clip[2],bin[col])-fixpt)/(clip[2]-fixpt));
+				bin[col]= MIN(maxave,bfrac*rgb[2]+(1-bfrac)*bin[col]);
+			}
+			
+			if ((L=(rin[col]+gin[col]+bin[col])/3) > desatpt) {
+				Lfrac = MAX(0,(maxave-L)/(maxave-desatpt));
+				C = Lfrac * 1.732050808 * (rin[col] - gin[col]);
+				H = Lfrac * (2 * bin[col] - rin[col] - gin[col]);
+				rin[col] = L - H / 6.0 + C / 3.464101615;
+				gin[col] = L - H / 6.0 - C / 3.464101615;
+				bin[col] = L + H / 3.0;
+			}
 		}
 	}
 
@@ -1932,7 +1975,7 @@ void RawImageSource::HLRecovery_CIELab (float* rin, float* gin, float* bin, floa
 	
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-void RawImageSource::hlRecovery (std::string method, float* red, float* green, float* blue, int i, int sx1, int width, int skip,const RAWParams &raw ) {
+void RawImageSource::hlRecovery (std::string method, float* red, float* green, float* blue, int i, int sx1, int width, int skip,const RAWParams &raw, float* hlmax ) {
 	
     if (method=="Luminance")
         HLRecovery_Luminance (red, green, blue, red, green, blue, width, 65535.0);
@@ -1943,7 +1986,7 @@ void RawImageSource::hlRecovery (std::string method, float* red, float* green, f
 	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);
-				HLRecovery_blend(red, green, blue, width, 65535.0, pre_mul, raw );}
+				HLRecovery_blend(red, green, blue, width, 65535.0, pre_mul, raw, hlmax );}
 
 }
 	

rtengine/rawimagesource.h

@@ -87,7 +87,7 @@ class RawImageSource : public ImageSource {
         float** hrmap[3];   // for color propagation
         char** needhr;      // for color propagation
 		int max[3];
-		float chmax[4];
+		float chmax[4],hlmax[4];
         double initialGain; // initial gain calculated after scale_colors
         double defGain;
         bool full;
@@ -115,7 +115,7 @@ class RawImageSource : public ImageSource {
         void hphd_horizontal     (float** hpmap, int row_from, int row_to);
         void hphd_green          (float** hpmap);
         void processFalseColorCorrectionThread (Imagefloat* im, int row_from, int row_to);
-        void hlRecovery          (std::string method, float* red, float* green, float* blue, int i, int sx1, int width, int skip, const RAWParams &raw);
+        void hlRecovery          (std::string method, float* red, float* green, float* blue, int i, int sx1, int width, int skip, const RAWParams &raw, float* hlmax);
         int  defTransform        (int tran);
         void rotateLine          (float* line, float** channel, int tran, int i, int w, int h);
         void transformRect       (PreviewProps pp, int tran, int &sx1, int &sy1, int &width, int &height, int &fw);
@@ -170,7 +170,7 @@ class RawImageSource : public ImageSource {
 	
         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);
+		static void HLRecovery_blend (float* rin, float* gin, float* bin, int width, float maxval, float* pre_mul, const RAWParams &raw, float* hlmax);
 
     protected:
         typedef unsigned short ushort;