diff --git a/rtengine/hilite_recon.cc b/rtengine/hilite_recon.cc index d5937cfdf..722469f9e 100644 --- a/rtengine/hilite_recon.cc +++ b/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" @@ -196,13 +197,23 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b int hfh = (height-(height%pitch))/pitch; int hfw = (width-(width%pitch))/pitch; - static const int numdirs = 4; + static const int numdirs = 4; + + //%%%%%%%%%%%%%%%%%%%% + //for blend algorithm: + 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 } } }; + //%%%%%%%%%%%%%%%%%%%% + - 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 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; - norm[c] = 1.0/(max[c]-fixthresh[c]); + max[c] = chmax[c]*maxpct; } 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]); + + // Determine the maximum level (clip) of all channels + const float fixpt = fixthreshpct*clippt; + const float desatpt = satthreshpct*maxave+(1-satthreshpct)*clippt; - 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 channelblur(width,height,ARRAY2D_CLEAR_DATA); multi_array2D hilite_full(width,height,ARRAY2D_CLEAR_DATA); // blur RGB channels - boxblur2(red ,channelblur[0],height,width,4); - boxblur2(green,channelblur[1],height,width,4); - boxblur2(blue ,channelblur[2],height,width,4); + boxblur2(red ,channelblur[0],height,width,3); + boxblur2(green,channelblur[1],height,width,3); + 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 @@ -477,24 +497,100 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b int i1 = MIN((i-(i%pitch))/pitch,hfh-1); for (int j=0; j 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 clipfix[3]={0,0,0}; + float dirwt, factor, Y, I, Q; + 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; dir0.5) { - dirwt = invfn[65535*(SQR(pixref[0]/Y-hilite_dir[dir*4+0][i1][j1]/Yhi) + \ - SQR(pixref[1]/Y-hilite_dir[dir*4+1][i1][j1]/Yhi) + \ - SQR(pixref[2]/Y-hilite_dir[dir*4+2][i1][j1]/Yhi))]; + dirwt = invfn[65535*(SQR(rgb_blend[0]/Y-hilite_dir[dir*4+0][i1][j1]/Yhi) + \ + SQR(rgb_blend[1]/Y-hilite_dir[dir*4+1][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 diff --git a/rtengine/rawimagesource.cc b/rtengine/rawimagesource.cc index 75f6adb12..8f45796bb 100644 --- a/rtengine/rawimagesource.cc +++ b/rtengine/rawimagesource.cc @@ -286,6 +286,11 @@ void RawImageSource::getImage (ColorTemp ctemp, int tran, Imagefloat* image, Pre rm/=area; 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 @@ -352,10 +357,10 @@ 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); @@ -1220,27 +1225,27 @@ void RawImageSource::copyOriginalPixels(const RAWParams &raw, RawImage *src, Raw } // flatfield - } else { - // No bayer pattern - // TODO: Is there a flat field correction possible? - if (!rawData) rawData = allocArray(3*W,H); - - if (riDark && W == riDark->get_width() && H == riDark->get_height()) { - for (int row = 0; row < H; row++) { - for (int col = 0; col < W; col++) { - rawData[row][3*col+0] = MAX (src->data[row][3*col+0]+ri->get_black() - riDark->data[row][3*col+0], 0); - rawData[row][3*col+1] = MAX (src->data[row][3*col+1]+ri->get_black() - riDark->data[row][3*col+1], 0); - rawData[row][3*col+2] = MAX (src->data[row][3*col+2]+ri->get_black() - riDark->data[row][3*col+2], 0); - } - } - } else { - for (int row = 0; row < H; row++) { - for (int col = 0; col < W; col++) { - rawData[row][3*col+0] = src->data[row][3*col+0]; - rawData[row][3*col+1] = src->data[row][3*col+1]; - rawData[row][3*col+2] = src->data[row][3*col+2]; - } - } + } else { + // No bayer pattern + // TODO: Is there a flat field correction possible? + if (!rawData) rawData = allocArray(3*W,H); + + if (riDark && W == riDark->get_width() && H == riDark->get_height()) { + for (int row = 0; row < H; row++) { + for (int col = 0; col < W; col++) { + rawData[row][3*col+0] = MAX (src->data[row][3*col+0]+ri->get_black() - riDark->data[row][3*col+0], 0); + rawData[row][3*col+1] = MAX (src->data[row][3*col+1]+ri->get_black() - riDark->data[row][3*col+1], 0); + rawData[row][3*col+2] = MAX (src->data[row][3*col+2]+ri->get_black() - riDark->data[row][3*col+2], 0); + } + } + } else { + for (int row = 0; row < H; row++) { + for (int col = 0; col < W; col++) { + rawData[row][3*col+0] = src->data[row][3*col+0]; + rawData[row][3*col+1] = src->data[row][3*col+1]; + rawData[row][3*col+2] = src->data[row][3*col+2]; + } + } } } } @@ -1771,74 +1776,112 @@ 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) -{ - const int ColorCount=3; - int clip=INT_MAX; - - // 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 } } }; +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; + // 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 } } }; + #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 + const float clippt = clipthresh*maxval; + const float fixpt = fixthresh*minpt; + const float desatpt = satthresh*maxave+(1-satthresh)*maxval; - // Determine the maximum level (clip) of all channels - int i; - FOREACHCOLOR if (clip > (i = (int) maxval * pre_mul[c] * raw.expos)) clip = i; #pragma omp parallel for - for (int col=0; col 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],maxval); + } + + // 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(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 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 clip) break; } - if (c == ColorCount) continue; - - // 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); - } - - // Calculate the lightness correction ration (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 - rin[col]=rgb[0]; gin[col]=rgb[1]; bin[col]=rgb[2]; - } -} + // Copy converted pixel back + 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; + } + } + } void RawImageSource::HLRecovery_Luminance (float* rin, float* gin, float* bin, float* rout, float* gout, float* bout, int width, float maxval) { @@ -1932,8 +1975,8 @@ 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); else if (method=="CIELab blending") @@ -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 );} } diff --git a/rtengine/rawimagesource.h b/rtengine/rawimagesource.h index 94ebd994b..855b963ff 100644 --- a/rtengine/rawimagesource.h +++ b/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;