From 1d4faecdfdbaaff2d64d2d19bdaee942f512ca29 Mon Sep 17 00:00:00 2001 From: Emil Martinec Date: Sun, 8 May 2011 08:47:39 -0500 Subject: [PATCH] For some reason improcfun.cc didn't get merged. --- rtengine/improcfun.cc | 916 ++++++++++++++++++++++++++---------------- 1 file changed, 559 insertions(+), 357 deletions(-) diff --git a/rtengine/improcfun.cc b/rtengine/improcfun.cc index 8da8614a1..30c476e14 100644 --- a/rtengine/improcfun.cc +++ b/rtengine/improcfun.cc @@ -32,15 +32,16 @@ #include #include #include +#include #ifdef _OPENMP #include #endif namespace rtengine { - -using namespace procparams; - + + using namespace procparams; + #undef MAXVAL #undef CMAXVAL #undef MAXL @@ -51,7 +52,7 @@ using namespace procparams; #undef CLIPS #undef CLIPC #undef CLIPTO - + #define MAXVAL 0xffff #define CMAXVAL 0xffff #define MAXL 0xffff @@ -62,70 +63,46 @@ using namespace procparams; #define CLIPS(a) ((a)>-32768?((a)<32767?(a):32767):-32768) #define CLIPC(a) ((a)>-32000?((a)<32000?(a):32000):-32000) #define CLIPTO(a,b,c) ((a)>(b)?((a)<(c)?(a):(c)):(b)) - +#define CLIP2(a) ((a)0.0?((a)<65535.5?(a):65535.5):0.0) + +#define D50x 0.96422 +#define D50z 0.82521 +#define u0 4.0*D50x/(D50x+15+3*D50z) +#define v0 9.0/(D50x+15+3*D50z) + +#define eps_max 580.40756 //(MAXVAL* 216.0f/24389.0); +#define kappa 903.29630 //24389.0/27.0; + + extern const Settings* settings; -int* ImProcFunctions::cacheL = 0; -int* ImProcFunctions::cachea = 0; -int* ImProcFunctions::cacheb = 0; -int* ImProcFunctions::xcache = 0; -int* ImProcFunctions::ycache = 0; -int* ImProcFunctions::zcache = 0; -unsigned short* ImProcFunctions::gamma2curve = 0; +LUTf ImProcFunctions::cachef ; +LUTf ImProcFunctions::gamma2curve = 0; void ImProcFunctions::initCache () { - const int maxindex = 2*65536; - cacheL = new int[maxindex]; - cachea = new int[maxindex]; - cacheb = new int[maxindex]; - gamma2curve = new unsigned short[65536]; + int maxindex = 65536; + cachef(maxindex,0/*LUT_CLIP_BELOW*/); - int threshold = (int)(0.008856*CMAXVAL); - for (int i=0; ithreshold) { - cacheL[i] = (int)round(655.35 * (116.0 * exp(1.0/3.0 * log((double)i / CMAXVAL)) - 16.0)); - cachea[i] = (int)round(32768.0 * 500.0 * exp(1.0/3.0 * log((double)i / CMAXVAL))); - cacheb[i] = (int)round(32768.0 * 200.0 * exp(1.0/3.0 * log((double)i / CMAXVAL))); + gamma2curve(maxindex); + + for (int i=0; ieps_max) { + cachef[i] = 327.68*( exp(1.0/3.0 * log((double)i / MAXVAL) )); } else { - cacheL[i] = (int)round(9033.0 * (double)i / 1000.0); // assuming CMAXVAL = 65535 - cachea[i] = (int)round(32768.0 * 500.0 * (7.787*i/CMAXVAL+16.0/116.0)); - cacheb[i] = (int)round(32768.0 * 200.0 * (7.787*i/CMAXVAL+16.0/116.0)); + cachef[i] = 327.68*((kappa*i/MAXVAL+16.0)/116.0); } + } - double fY; - ycache = new int[0x10000]; - for (int i=0; i<0x10000; i++) - ycache[i] = (int)round(65536.0 * ((fY=((double)i/655.35+16)/116) > 2.0689655172413793e-1 ? fY*fY*fY : 1.107056459879453852e-3*(double)i/655.35)); - for (int i=0; i<0x10000; i++) - ycache[i] = CLIP(ycache[i]); - xcache = new int[369621]; - for (int i=-141556; i<228064; i++) - xcache[i+141556] = (int)round(65536.0 * (i > 15728 ? ((double)i/76021)*((double)i/76021)*((double)i/76021)*0.96422 : (1.2841854934601665e-1*(double)i/76021-1.7712903358071262e-2)*0.96422)); - for (int i=0; i<369620; i++) - xcache[i] = CLIP(xcache[i]); - zcache = new int[825747]; - for (int i=-369619; i<456128; i++) - zcache[i+369619] = (int)round(65536.0 * (i > 15728 ? ((double)i/76021)*((double)i/76021)*((double)i/76021)*0.82521 : (1.2841854934601665e-1*(double)i/76021-1.7712903358071262e-2)*0.82521)); - for (int i=0; i<825747; i++) - zcache[i] = CLIP(zcache[i]); - - for (int i=0; i<65536; i++) { - int g = (int)(CurveFactory::gamma2(i/65535.0) * 65535.0); - gamma2curve[i] = CLIP(g); + for (int i=0; iworkingSpaceMatrix (wprofile); - lumimul[0] = wprof[0][1]; + lumimul[0] = wprof[1][0]; lumimul[1] = wprof[1][1]; - lumimul[2] = wprof[2][1]; + lumimul[2] = wprof[1][2]; int W = original->width; - int cradius = 1; for (int i=row_from; ig[i][j]; int b = original->b[i][j]; - int x = (toxyz[0][0] * r + toxyz[1][0] * g + toxyz[2][0] * b) >> 15; - int y = (toxyz[0][1] * r + toxyz[1][1] * g + toxyz[2][1] * b) >> 15; - int z = (toxyz[0][2] * r + toxyz[1][2] * g + toxyz[2][2] * b) >> 15; + int y = CLIP((int)(lumimul[0] * r + lumimul[1] * g + lumimul[2] * b)) ; - x = CLIPTO(x,0,2*65536-1); - y = CLIPTO(y,0,2*65536-1); - z = CLIPTO(z,0,2*65536-1); - - int oa = cachea[x] - cachea[y]; - int ob = cacheb[y] - cacheb[z]; - - if (oa<0) oa = -oa; - if (ob<0) ob = -ob; - - if (oa > cradius) - cradius = oa; - if (ob > cradius) - cradius = ob; - - if (histogram) { - int hval = CLIP(y); //(306 * original->r[i][j] + 601 * original->g[i][j] + 117 * original->b[i][j]) >> 10; - histogram[hval]++; + if (histogram) { + histogram[y]++; } } } - *chroma_radius = cradius; } -void ImProcFunctions::firstAnalysis (Image16* original, const ProcParams* params, unsigned int* histogram, double gamma) { +void ImProcFunctions::firstAnalysis (Imagefloat* original, const ProcParams* params, LUTu & histogram, double gamma) { // set up monitor transform - TMatrix wprof = iccStore->workingSpaceMatrix (params->icm.working); + Glib::ustring wprofile = params->icm.working; if (monitorTransform) cmsDeleteTransform (monitorTransform); monitorTransform = NULL; - Glib::ustring monitorProfile=settings->monitorProfile; - if (settings->autoMonitorProfile) monitorProfile=iccStore->defaultMonitorProfile; - //if (settings->verbose) printf("Using monitor profile: %s\n", monitorProfile.c_str()); + Glib::ustring monitorProfile=settings->monitorProfile; + if (settings->autoMonitorProfile) monitorProfile=iccStore->defaultMonitorProfile; cmsHPROFILE monitor = iccStore->getProfile ("file:"+monitorProfile); if (monitor) { cmsHPROFILE iprof = iccStore->getXYZProfile (); - cmsHPROFILE oprof = iccStore->getProfile (params->icm.output); - if (!oprof) - oprof = iccStore->getsRGBProfile (); lcmsMutex->lock (); - monitorTransform = cmsCreateTransform (iprof, TYPE_RGB_16, monitor, TYPE_RGB_8, settings->colorimetricIntent, 0); + monitorTransform = cmsCreateTransform (iprof, TYPE_RGB_FLT, monitor, TYPE_RGB_8, settings->colorimetricIntent, + settings->LCMSSafeMode ? cmsFLAGS_NOOPTIMIZE : cmsFLAGS_NOOPTIMIZE | cmsFLAGS_NOCACHE ); // NOCACHE is important for thread safety lcmsMutex->unlock (); } + + //chroma_scale = 1; - // calculate chroma radius and histogram of the y channel needed for exposure curve calculation + // calculate histogram of the y channel needed for contrast curve calculation in exposure adjustments #ifdef _OPENMP int T = omp_get_max_threads(); @@ -223,10 +171,8 @@ void ImProcFunctions::firstAnalysis (Image16* original, const ProcParams* params int T = 1; #endif - int* cr = new int [T]; unsigned int** hist = new unsigned int* [T]; for (int i=0; iheight); + firstAnalysisThread (original, wprofile, hist[0], 0, original->height); #endif - chroma_radius = cr[0]; - for (int i=0; ichroma_radius) - chroma_radius = cr[i]; - memset (histogram, 0, 65536*sizeof(int)); + histogram.clear(); for (int i=0; i<65536; i++) for (int j=0; jsh.localcontrast / 200.0; TMatrix wprof = iccStore->workingSpaceMatrix (params->icm.working); - int toxyz[3][3] = { + + float toxyz[3][3] = { { - floor(32768.0 * wprof[0][0] / 0.96422), - floor(32768.0 * wprof[0][1]), - floor(32768.0 * wprof[0][2] / 0.82521) + ( wprof[0][0] / D50x), + ( wprof[0][1] / D50x), + ( wprof[0][2] / D50x) },{ - floor(32768.0 * wprof[1][0] / 0.96422), - floor(32768.0 * wprof[1][1]), - floor(32768.0 * wprof[1][2] / 0.82521) + ( wprof[1][0] ), + ( wprof[1][1] ), + ( wprof[1][2] ) },{ - floor(32768.0 * wprof[2][0] / 0.96422), - floor(32768.0 * wprof[2][1]), - floor(32768.0 * wprof[2][2] / 0.82521) + ( wprof[2][0] / D50z), + ( wprof[2][1] / D50z), + ( wprof[2][2] / D50z) } }; - bool mixchannels = params->chmixer.red[0]!=100 || params->chmixer.red[1]!=0 || params->chmixer.red[2]!=0 || params->chmixer.green[0]!=0 || params->chmixer.green[1]!=100 || params->chmixer.green[2]!=0 || params->chmixer.blue[0]!=0 || params->chmixer.blue[1]!=0 || params->chmixer.blue[2]!=100; - int mapval; - double factor; + bool mixchannels = (params->chmixer.red[0]!=100 || params->chmixer.red[1]!=0 || params->chmixer.red[2]!=0 || \ + params->chmixer.green[0]!=0 || params->chmixer.green[1]!=100 || params->chmixer.green[2]!=0 || \ + params->chmixer.blue[0]!=0 || params->chmixer.blue[1]!=0 || params->chmixer.blue[2]!=100); + int tW = working->width; int tH = working->height; - int r, g, b; - float h, s, v; double pi = M_PI; FlatCurve* hCurve; FlatCurve* sCurve; FlatCurve* vCurve; - float* cossq = new float [8093]; - for (int i=0; i<8093; i++) - cossq[i] = SQR(cos(pi*(float)i/16384)); + + float* cossq = new float [8192]; + for (int i=0; i<8192; i++) + cossq[i] = SQR(cos(pi*(float)i/16384.0)); FlatCurveType hCurveType = (FlatCurveType)params->hsvequalizer.hcurve.at(0); FlatCurveType sCurveType = (FlatCurveType)params->hsvequalizer.scurve.at(0); @@ -323,28 +264,34 @@ void ImProcFunctions::rgbProc (Image16* working, LabImage* lab, float* hltonecur if (hCurveEnabled) hCurve = new FlatCurve(params->hsvequalizer.hcurve); if (sCurveEnabled) sCurve = new FlatCurve(params->hsvequalizer.scurve); if (vCurveEnabled) vCurve = new FlatCurve(params->hsvequalizer.vcurve); - -#pragma omp parallel for private(r, g, b,factor,mapval,h,s,v) if (multiThread) + + 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 hlrange = 65536.0-shoulder; + + +#pragma omp parallel for if (multiThread) for (int i=0; ir[i][j]; - g = working->g[i][j]; - b = working->b[i][j]; + float r = working->r[i][j]; + float g = working->g[i][j]; + float b = working->b[i][j]; + + //if (i==100 & j==100) printf("rgbProc input R= %f G= %f B= %f \n",r,g,b); if (mixchannels) { - int newr = (r*params->chmixer.red[0] + g*params->chmixer.red[1] + b*params->chmixer.red[2]) / 100; - int newg = (r*params->chmixer.green[0] + g*params->chmixer.green[1] + b*params->chmixer.green[2]) / 100; - int newb = (r*params->chmixer.blue[0] + g*params->chmixer.blue[1] + b*params->chmixer.blue[2]) / 100; - r = CLIP(newr); - g = CLIP(newg); - b = CLIP(newb); + r = (r*params->chmixer.red[0] + g*params->chmixer.red[1] + b*params->chmixer.red[2]) / 100; + g = (r*params->chmixer.green[0] + g*params->chmixer.green[1] + b*params->chmixer.green[2]) / 100; + b = (r*params->chmixer.blue[0] + g*params->chmixer.blue[1] + b*params->chmixer.blue[2]) / 100; + } if (processSH || processLCE) { - mapval = shmap->map[i][j]; - factor = 1.0; + double mapval = 1.0 + shmap->map[i][j]; + double factor = 1.0; if (processSH) { if (mapval > h_th) @@ -354,37 +301,41 @@ void ImProcFunctions::rgbProc (Image16* working, LabImage* lab, float* hltonecur } if (processLCE) { double sub = lceamount*(mapval-factor*(r*lumimul[0] + g*lumimul[1] + b*lumimul[2])); - r = CLIP((int)(factor*r-sub)); - g = CLIP((int)(factor*g-sub)); - b = CLIP((int)(factor*b-sub)); + r = CLIP(factor*r-sub); + g = CLIP(factor*g-sub); + b = CLIP(factor*b-sub); } else { - r = CLIP((int)(factor*r)); - g = CLIP((int)(factor*g)); - b = CLIP((int)(factor*b)); + r = CLIP(factor*r); + g = CLIP(factor*g); + b = CLIP(factor*b); } } - - float tonefactor=(hltonecurve[r]+hltonecurve[g]+hltonecurve[b])/3; - + //TODO: proper treatment of out-of-gamut colors + //float tonefactor = hltonecurve[(0.299f*r+0.587f*g+0.114f*b)]; + float tonefactor=((r0 ? (float)shtonecurve[Y]/Y : 1); + float Y = (0.299*r + 0.587*g + 0.114*b); + tonefactor = shtonecurve[Y]; r *= tonefactor; g *= tonefactor; b *= tonefactor; //brightness/contrast and user tone curve - r = tonecurve[CLIP(r)]; - g = tonecurve[CLIP(g)]; - b = tonecurve[CLIP(b)]; + r = tonecurve[r]; + g = tonecurve[g]; + b = tonecurve[b]; if (abs(sat)>0.5 || hCurveEnabled || sCurveEnabled || vCurveEnabled) { + float h,s,v; rgb2hsv(r,g,b,h,s,v); if (sat > 0.5) { s = (1-(float)sat/100)*s+(float)sat/100*(1-SQR(SQR(1-s))); @@ -409,12 +360,12 @@ void ImProcFunctions::rgbProc (Image16* working, LabImage* lab, float* hltonecur if (satparam < -0.00001) s *= 1+satparam; } - + /*s = sCurve->getVal((double)s); - if (s > 1.0) - s -= 1.0; - else if (s < 0.0) - s += 1.0;*/ + if (s > 1.0) + s -= 1.0; + else if (s < 0.0) + s += 1.0;*/ } if (vCurveEnabled) { //shift value @@ -426,30 +377,55 @@ void ImProcFunctions::rgbProc (Image16* working, LabImage* lab, float* hltonecur if (valparam < -0.00001) v *= (1+valparam); } - + /*v = vCurve->getVal((double)v); - if (v > 1.0) - v -= 1.0; - else if (v < 0.0) - v += 1.0;*/ + if (v > 1.0) + v -= 1.0; + else if (v < 0.0) + v += 1.0;*/ } hsv2rgb(h,s,v,r,g,b); } - //hsv2rgb(h,s,v,r,g,b); + + //r=FCLIP(r); + //g=FCLIP(g); + //b=FCLIP(b); + + float x = (toxyz[0][0] * r + toxyz[0][1] * g + toxyz[0][2] * b) ; + float y = (toxyz[1][0] * r + toxyz[1][1] * g + toxyz[1][2] * b) ; + float z = (toxyz[2][0] * r + toxyz[2][1] * g + toxyz[2][2] * b) ; + + /*lab->L[i][j] = CLIP2(116.0 * (CurveFactory::flinterp(cachef,y)) - 5242.88); //5242.88=16.0*327.68; + lab->a[i][j] = CLIPS(500.0 * (((CurveFactory::flinterp(cachef,x) - CurveFactory::flinterp(cachef,y)) ) )); + lab->b[i][j] = CLIPS(200.0 * (((CurveFactory::flinterp(cachef,y) - CurveFactory::flinterp(cachef,z)) ) ));*/ + + x = (x<65535.0 ? cachef[x] : (327.68*exp(log(x/MAXVAL)/3.0 ))); + y = (y<65535.0 ? cachef[y] : (327.68*exp(log(y/MAXVAL)/3.0 ))); + z = (z<65535.0 ? cachef[z] : (327.68*exp(log(z/MAXVAL)/3.0 ))); + lab->L[i][j] = (116.0 * y - 5242.88); //5242.88=16.0*327.68; + lab->a[i][j] = (500.0 * (x - y) ); + lab->b[i][j] = (200.0 * (y - z) ); - int x = (toxyz[0][0] * r + toxyz[1][0] * g + toxyz[2][0] * b) >> 15; - int y = (toxyz[0][1] * r + toxyz[1][1] * g + toxyz[2][1] * b) >> 15; - int z = (toxyz[0][2] * r + toxyz[1][2] * g + toxyz[2][2] * b) >> 15; + //float L1 = lab->L[i][j];//for testing + //float a1 = lab->a[i][j]; + //float b1 = lab->b[i][j]; + //float xxx=1; + + //test for color accuracy + /*float fy = (0.00862069 * lab->L[i][j])/327.68 + 0.137932; // (L+16)/116 + float fx = (0.002 * lab->a[i][j])/327.68 + fy; + float fz = fy - (0.005 * lab->b[i][j])/327.68; + + float x_ = 65535*Lab2xyz(fx)*D50x; + float y_ = 65535*Lab2xyz(fy); + float z_ = 65535*Lab2xyz(fz)*D50z; + + int R,G,B; + xyz2srgb(x_,y_,z_,R,G,B); + r=(float)R; g=(float)G; b=(float)B; + float xxx=1;*/ - x = CLIPTO(x,0,2*65536-1); - y = CLIPTO(y,0,2*65536-1); - z = CLIPTO(z,0,2*65536-1); - - int L = cacheL[y]; - lab->L[i][j] = L; - lab->a[i][j] = CLIPC(((cachea[x] - cachea[y]) * chroma_scale) >> 15); - lab->b[i][j] = CLIPC(((cacheb[y] - cacheb[z]) * chroma_scale) >> 15); } } @@ -457,71 +433,128 @@ void ImProcFunctions::rgbProc (Image16* working, LabImage* lab, float* hltonecur if (sCurveEnabled) delete sCurve; if (vCurveEnabled) delete vCurve; delete [] cossq; - //delete [] my_tonecurve; } -void ImProcFunctions::luminanceCurve (LabImage* lold, LabImage* lnew, int* curve) { +void ImProcFunctions::luminanceCurve (LabImage* lold, LabImage* lnew, LUTf & curve) { int W = lold->W; int H = lold->H; for (int i=0; iL[i][j] = curve[lold->L[i][j]]; + for (int j=0; jL[i][j]; + //if (Lin>0 && Lin<65535) + lnew->L[i][j] = curve[Lin]; + } } -void ImProcFunctions::chrominanceCurve (LabImage* lold, LabImage* lnew, float* acurve, float* bcurve) { +void ImProcFunctions::chrominanceCurve (LabImage* lold, LabImage* lnew, LUTf & acurve, LUTf & bcurve, LUTf & satcurve/*, double sat*/) { int W = lold->W; int H = lold->H; + /*for (int i=0; ia[i][j] = acurve[lold->a[i][j]+32768]-32768; - lnew->b[i][j] = bcurve[lold->b[i][j]+32768]-32768; - }*/ + float ain=lold->a[i][j]; + float bin=lold->b[i][j]; + //if (fabs(ain)<32767 && fabs(bin)<32767) + lnew->a[i][j] = acurve[ain+32768.0f]-32768.0; + lnew->b[i][j] = bcurve[bin+32768.0f]-32768.0; + }*/ + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + TMatrix wprof = iccStore->workingSpaceMatrix (params->icm.working); + + double wp[3][3] = { + {wprof[0][0],wprof[0][1],wprof[0][2]}, + {wprof[1][0],wprof[1][1],wprof[1][2]}, + {wprof[2][0],wprof[2][1],wprof[2][2]}}; + #pragma omp parallel for if (multiThread) for (int i=0; ia[i][j]; - int ob = lold->b[i][j]; + float atmp = acurve[lold->a[i][j]+32768.0f]-32768.0f; + float btmp = bcurve[lold->b[i][j]+32768.0f]-32768.0f; - float atmp = acurve[oa+32768]-32768; - float btmp = bcurve[ob+32768]-32768; - - double real_c = 1.0; + if (params->labCurve.saturation) { + float chroma = sqrt(SQR(atmp)+SQR(btmp)+0.001); + float satfactor = (satcurve[chroma+32768.0f]-32768.0f)/chroma; + atmp *= satfactor; + btmp *= satfactor; + } + + //double real_c = 1.0; if (params->labCurve.avoidclip) { - double Lclip = MIN(lnew->L[i][j]/655.35,100.0); - double cr = tightestroot (Lclip, (double)atmp/chroma_scale, (double)btmp/chroma_scale, 3.079935, -1.5371515, -0.54278342); - double cg = tightestroot (Lclip, (double)atmp/chroma_scale, (double)btmp/chroma_scale, -0.92123418, 1.87599, 0.04524418); - double cb = tightestroot (Lclip, (double)atmp/chroma_scale, (double)btmp/chroma_scale, 0.052889682, -0.20404134, 1.15115166); + //Luv limiter + float Y,u,v; + Lab2Yuv(lnew->L[i][j],atmp,btmp,Y,u,v); + //Yuv2Lab includes gamut restriction map + Yuv2Lab(Y,u,v,lnew->L[i][j],lnew->a[i][j],lnew->b[i][j], wp); + + //Gabor's limiter -- needs adapting to float branch + /*double Lclip = MIN(lnew->L[i][j]/655.35,100.0); + double cr = tightestroot (Lclip, (double)atmp, (double)btmp, 3.079935, -1.5371515, -0.54278342); + double cg = tightestroot (Lclip, (double)atmp, (double)btmp, -0.92123418, 1.87599, 0.04524418); + double cb = tightestroot (Lclip, (double)atmp, (double)btmp, 0.052889682, -0.20404134, 1.15115166); if (cr>0 && cr0 && cg0 && cb0 && cba[i][j] = atmp; + lnew->b[i][j] = btmp; + } + /* + //Gabor again int nna = (int)((atmp) * real_c ); int nnb = (int)((btmp) * real_c ); lnew->a[i][j] = CLIPTO(nna,-32000,32000); - lnew->b[i][j] = CLIPTO(nnb,-32000,32000); + lnew->b[i][j] = CLIPTO(nnb,-32000,32000);*/ } //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + /* + float slope = (sat>0) ? (125.0+sat)/(125.0-sat) : (1+sat/100.0); + + TMatrix wprof = iccStore->workingSpaceMatrix (params->icm.working); + + double wp[3][3] = { + {wprof[0][0],wprof[0][1],wprof[0][2]}, + {wprof[1][0],wprof[1][1],wprof[1][2]}, + {wprof[2][0],wprof[2][1],wprof[2][2]}}; + + int W = lold->W; + int H = lold->H; + for (int i=0; ia[i][j]+32768.0f]-32768.0; + float bin = bcurve[lold->b[i][j]+32768.0f]-32768.0; + Lab2Yuv(lold->L[i][j],ain,bin,Y,u,v); + u *= slope; + v *= slope; + Yuv2Lab(Y,u,v,lnew->L[i][j],lnew->a[i][j],lnew->b[i][j], wp); + //float ain=lold->a[i][j]; + //float bin=lold->b[i][j]; + //if (fabs(ain)<32767 && fabs(bin)<32767) + //lnew->a[i][j] = acurve[ain+32768.0f]-32768.0; + //lnew->b[i][j] = bcurve[bin+32768.0f]-32768.0; + }*/ + } -} #include "cubic.cc" void ImProcFunctions::colorCurve (LabImage* lold, LabImage* lnew) { - /*double* cmultiplier = new double [181021]; +/* LUT cmultiplier(181021); - double boost_a = (params->colorBoost.amount + 100.0) / 100.0; - double boost_b = (params->colorBoost.amount + 100.0) / 100.0; + double boost_a = ((float)params->colorBoost.amount + 100.0) / 100.0; + double boost_b = ((float)params->colorBoost.amount + 100.0) / 100.0; double c, amul = 1.0, bmul = 1.0; if (boost_a > boost_b) { @@ -535,14 +568,14 @@ void ImProcFunctions::colorCurve (LabImage* lold, LabImage* lnew) { amul = boost_a / boost_b; } - if (params->colorBoost.enable_saturationlimiter && c>1) { + if (params->colorBoost.enable_saturationlimiter && c>1.0) { // re-generate color multiplier lookup table - double d = params->colorBoost.saturationlimit * chroma_scale / 3.0; + double d = params->colorBoost.saturationlimit / 3.0; double alpha = 0.5; double threshold1 = alpha * d; double threshold2 = c*d*(alpha+1.0) - d; for (int i=0; i<=181020; i++) { // lookup table stores multipliers with a 0.25 chrominance resolution - double chrominance = (double)i/4; + double chrominance = (double)i/4.0; if (chrominance < threshold1) cmultiplier[i] = c; else if (chrominance < d) @@ -555,10 +588,10 @@ void ImProcFunctions::colorCurve (LabImage* lold, LabImage* lnew) { } float eps = 0.001; - double shift_a = params->colorShift.a * chroma_scale + eps, shift_b = params->colorShift.b * chroma_scale + eps; + double shift_a = params->colorShift.a + eps, shift_b = params->colorShift.b + eps; - short** oa = lold->a; - short** ob = lold->b; + float** oa = lold->a; + float** ob = lold->b; #pragma omp parallel for if (multiThread) for (int i=0; iH; i++) @@ -566,40 +599,40 @@ void ImProcFunctions::colorCurve (LabImage* lold, LabImage* lnew) { double wanted_c = c; if (params->colorBoost.enable_saturationlimiter && c>1) { - int chroma = (int)(4.0 * sqrt((oa[i][j]+shift_a)*(oa[i][j]+shift_a) + (ob[i][j]+shift_b)*(ob[i][j]+shift_b))); - wanted_c = cmultiplier [MIN(chroma,181020)]; + float chroma = (float)(4.0 * sqrt((oa[i][j]+shift_a)*(oa[i][j]+shift_a) + (ob[i][j]+shift_b)*(ob[i][j]+shift_b))); + wanted_c = cmultiplier [chroma]; } double real_c = wanted_c; if (wanted_c >= 1.0 && params->colorBoost.avoidclip) { - double cclip = 100000; - double cr = tightestroot ((double)lnew->L[i][j]/655.35, (double)(oa[i][j]+shift_a)/chroma_scale*amul, (double)(ob[i][j]+shift_b)/chroma_scale*bmul, 3.079935, -1.5371515, -0.54278342); - double cg = tightestroot ((double)lnew->L[i][j]/655.35, (double)(oa[i][j]+shift_a)/chroma_scale*amul, (double)(ob[i][j]+shift_b)/chroma_scale*bmul, -0.92123418, 1.87599, 0.04524418); - double cb = tightestroot ((double)lnew->L[i][j]/655.35, (double)(oa[i][j]+shift_a)/chroma_scale*amul, (double)(ob[i][j]+shift_b)/chroma_scale*bmul, 0.052889682, -0.20404134, 1.15115166); + double cclip = 100000.0; + double cr = tightestroot ((double)lnew->L[i][j]/655.35, (double)(oa[i][j]+shift_a)*amul, (double)(ob[i][j]+shift_b)*bmul, 3.079935, -1.5371515, -0.54278342); + double cg = tightestroot ((double)lnew->L[i][j]/655.35, (double)(oa[i][j]+shift_a)*amul, (double)(ob[i][j]+shift_b)*bmul, -0.92123418, 1.87599, 0.04524418); + double cb = tightestroot ((double)lnew->L[i][j]/655.35, (double)(oa[i][j]+shift_a)*amul, (double)(ob[i][j]+shift_b)*bmul, 0.052889682, -0.20404134, 1.15115166); if (cr>1.0 && cr1.0 && cg1.0 && cba[i][j] = CLIPTO(nna,-32000,32000); - lnew->b[i][j] = CLIPTO(nnb,-32000,32000); + float nna = ((oa[i][j]+shift_a) * real_c * amul); + float nnb = ((ob[i][j]+shift_b) * real_c * bmul); + lnew->a[i][j] = CLIPTO(nna,-32000.0f,32000.0f); + lnew->b[i][j] = CLIPTO(nnb,-32000.0f,32000.0f); } - - delete [] cmultiplier;*/ +*/ + //delete [] cmultiplier; } void ImProcFunctions::impulsedenoise (LabImage* lab) { if (params->impulseDenoise.enabled && lab->W>=8 && lab->H>=8) - impulse_nr (lab->L, lab->L, lab->W, lab->H, (float)params->impulseDenoise.thresh/20.0 ); + impulse_nr (lab, (float)params->impulseDenoise.thresh/20.0 ); } void ImProcFunctions::defringe (LabImage* lab) { @@ -613,7 +646,7 @@ void ImProcFunctions::colorCurve (LabImage* lold, LabImage* lnew) { if (params->dirpyrDenoise.enabled && lab->W>=8 && lab->H>=8) - dirpyrLab_denoise(lab, lab, params->dirpyrDenoise.luma, params->dirpyrDenoise.chroma, params->dirpyrDenoise.gamma/3.0 ); + dirpyrLab_denoise(lab, lab, params->dirpyrDenoise ); } void ImProcFunctions::dirpyrequalizer (LabImage* lab) { @@ -622,145 +655,314 @@ void ImProcFunctions::colorCurve (LabImage* lold, LabImage* lnew) { //dirpyrLab_equalizer(lab, lab, params->dirpyrequalizer.mult); dirpyr_equalizer(lab->L, lab->L, lab->W, lab->H, params->dirpyrequalizer.mult); - + } } - -void ImProcFunctions::lumadenoise (LabImage* lab, int** b2) { - - /*if (params->lumaDenoise.enabled && lab->W>=8 && lab->H>=8) + + void ImProcFunctions::lumadenoise (LabImage* lab, int** b2) { + + if (params->lumaDenoise.enabled && lab->W>=8 && lab->H>=8) #ifdef _OPENMP #pragma omp parallel #endif - bilateral (lab->L, lab->L, (unsigned short**)b2, lab->W, lab->H, \ - params->lumaDenoise.radius / scale, params->lumaDenoise.edgetolerance, multiThread); - */ -} - -void ImProcFunctions::colordenoise (LabImage* lab, int** b2) { - /* - if (params->colorDenoise.enabled && lab->W>=8 && lab->H>=8) { + bilateral (lab->L, lab->L, (float**)b2, lab->W, lab->H, params->lumaDenoise.radius / scale, params->lumaDenoise.edgetolerance, multiThread); + } + + void ImProcFunctions::colordenoise (LabImage* lab, int** b2) { + + if (params->colorDenoise.enabled && lab->W>=8 && lab->H>=8) { #ifdef _OPENMP #pragma omp parallel #endif - { - AlignedBuffer* buffer = new AlignedBuffer (MAX(lab->W,lab->H)); - gaussHorizontal (lab->a, lab->a, buffer, lab->W, lab->H, params->colorDenoise.amount / 10.0 / scale, multiThread); - gaussHorizontal (lab->b, lab->b, buffer, lab->W, lab->H, params->colorDenoise.amount / 10.0 / scale, multiThread); - gaussVertical (lab->a, lab->a, buffer, lab->W, lab->H, params->colorDenoise.amount / 10.0 / scale, multiThread); - gaussVertical (lab->b, lab->b, buffer, lab->W, lab->H, params->colorDenoise.amount / 10.0 / scale, multiThread); - - delete buffer; - } - } - */ -} - -void ImProcFunctions::getAutoExp (unsigned int* histogram, int histcompr, double expcomp, double clip, double& br, int& bl) { - - double sum = 0; - for (int i=0; i<65536>>histcompr; i++) - sum += histogram[i]; - - // compute clipping points based on the original histograms (linear, without exp comp.) - int clippable = (int)(sum * clip); - int clipped = 0; - int aw = (65536>>histcompr) - 1; - while (aw>1 && histogram[aw]+clipped <= clippable) { - clipped += histogram[aw]; - aw--; - } - - clipped = 0; - int shc = 0; - while (shc>histcompr; i++) - gavg += histogram[i] * CurveFactory::gamma2((int)(corr*(i<* buffer = new AlignedBuffer (MAX(lab->W,lab->H)); + gaussHorizontal (lab->a, lab->a, buffer, lab->W, lab->H, params->colorDenoise.amount / 10.0 / scale, multiThread); + gaussHorizontal (lab->b, lab->b, buffer, lab->W, lab->H, params->colorDenoise.amount / 10.0 / scale, multiThread); + gaussVertical (lab->a, lab->a, buffer, lab->W, lab->H, params->colorDenoise.amount / 10.0 / scale, multiThread); + gaussVertical (lab->b, lab->b, buffer, lab->W, lab->H, params->colorDenoise.amount / 10.0 / scale, multiThread); + + delete buffer; + } + } + } - awg = CurveFactory::igamma2 ((float)(awg/65535.0)) * 65535.0; //need to inverse gamma transform to get correct exposure compensation parameter - - bl = (int)((65535*bl)/awg); - br = log(65535.0*corr / (awg)) / log(2.0); - if (br<0) br = 0; - if (br>10) br=10; -} + void ImProcFunctions::getAutoExp (LUTu & histogram, int histcompr, double expcomp, double clip, double& br, int& bl) { + + double sum = 0; + for (int i=0; i<65536>>histcompr; i++) + sum += histogram[i]; + + // compute clipping points based on the original histograms (linear, without exp comp.) + int clippable = (int)(sum * clip); + int clipped = 0; + int aw = (65536>>histcompr) - 1; + while (aw>1 && histogram[aw]+clipped <= clippable) { + clipped += histogram[aw]; + aw--; + } + + clipped = 0; + int shc = 0; + while (shc>histcompr; i++) + gavg += histogram[i] * CurveFactory::gamma2((int)(corr*(i<10.0) br=10.0; + } -#include "calc_distort.h" - -double ImProcFunctions::getAutoDistor (const Glib::ustring &fname, int thumb_size) { - if (fname != "") { - rtengine::RawMetaDataLocation ri; - int w_raw=-1, h_raw=thumb_size; - int w_thumb=-1, h_thumb=thumb_size; - - Thumbnail* thumb = rtengine::Thumbnail::loadQuickFromRaw (fname, ri, w_thumb, h_thumb, 1, FALSE); - if (thumb == NULL) - return 0.0; - - Thumbnail* raw = rtengine::Thumbnail::loadFromRaw (fname, ri, w_raw, h_raw, 1, FALSE); - if (raw == NULL) { - delete thumb; - return 0.0; - } - if (h_thumb != h_raw) { - delete thumb; - delete raw; - return 0.0; - } - - int width; - - if (w_thumb > w_raw) - width = w_raw; - else - width = w_thumb; - - unsigned char* thumbGray; - unsigned char* rawGray; - thumbGray = thumb->getGrayscaleHistEQ (width); - rawGray = raw->getGrayscaleHistEQ (width); - - if (thumbGray == NULL || rawGray == NULL) { - if (thumbGray) delete thumbGray; - if (rawGray) delete rawGray; - delete thumb; - delete raw; - return 0.0; - } - - double dist_amount = calcDistortion (thumbGray, rawGray, width, h_thumb); - delete thumbGray; - delete rawGray; - delete thumb; - delete raw; - return dist_amount; - } - else - return 0.0; -} + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + #include "calc_distort.h" + + double ImProcFunctions::getAutoDistor (const Glib::ustring &fname, int thumb_size) { + if (fname != "") { + rtengine::RawMetaDataLocation ri; + int w_raw=-1, h_raw=thumb_size; + int w_thumb=-1, h_thumb=thumb_size; + + Thumbnail* thumb = rtengine::Thumbnail::loadQuickFromRaw (fname, ri, w_thumb, h_thumb, 1, FALSE); + if (thumb == NULL) + return 0.0; + + Thumbnail* raw = rtengine::Thumbnail::loadFromRaw (fname, ri, w_raw, h_raw, 1, FALSE); + if (raw == NULL) { + delete thumb; + return 0.0; + } + + if (h_thumb != h_raw) { + delete thumb; + delete raw; + return 0.0; + } + + int width; + + if (w_thumb > w_raw) + width = w_raw; + else + width = w_thumb; + + unsigned char* thumbGray; + unsigned char* rawGray; + thumbGray = thumb->getGrayscaleHistEQ (width); + rawGray = raw->getGrayscaleHistEQ (width); + + if (thumbGray == NULL || rawGray == NULL) { + if (thumbGray) delete thumbGray; + if (rawGray) delete rawGray; + delete thumb; + delete raw; + return 0.0; + } + + double dist_amount = calcDistortion (thumbGray, rawGray, width, h_thumb); + delete thumbGray; + delete rawGray; + delete thumb; + delete raw; + return dist_amount; + } + else + return 0.0; + } + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + void ImProcFunctions::rgb2hsv (float r, float g, float b, float &h, float &s, float &v) { + + double var_R = r / 65535.0; + double var_G = g / 65535.0; + double var_B = b / 65535.0; + + double var_Min = MIN(MIN(var_R,var_G),var_B); + double var_Max = MAX(MAX(var_R,var_G),var_B); + double del_Max = var_Max - var_Min; + v = var_Max; + if (fabs(del_Max)<0.00001) { + h = 0; + s = 0; + } + else { + s = del_Max/var_Max; + + if ( var_R == var_Max ) h = (var_G - var_B)/del_Max; + else if ( var_G == var_Max ) h = 2.0 + (var_B - var_R)/del_Max; + else if ( var_B == var_Max ) h = 4.0 + (var_R - var_G)/del_Max; + h /= 6.0; + + if ( h < 0 ) h += 1; + if ( h > 1 ) h -= 1; + } + } + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + void ImProcFunctions::hsv2rgb (float h, float s, float v, float &r, float &g, float &b) { + + float h1 = h*6; // sector 0 to 5 + int i = floor( h1 ); + float f = h1 - i; // fractional part of h + + float p = v * ( 1 - s ); + float q = v * ( 1 - s * f ); + float t = v * ( 1 - s * ( 1 - f ) ); + + float r1,g1,b1; + + if (i==0) {r1 = v; g1 = t; b1 = p;} + if (i==1) {r1 = q; g1 = v; b1 = p;} + if (i==2) {r1 = p; g1 = v; b1 = t;} + if (i==3) {r1 = p; g1 = q; b1 = v;} + if (i==4) {r1 = t; g1 = p; b1 = v;} + if (i==5) {r1 = v; g1 = p; b1 = q;} + + r = ((r1)*65535.0); + g = ((g1)*65535.0); + b = ((b1)*65535.0); + } + + void ImProcFunctions::xyz2srgb (float x, float y, float z, float &r, float &g, float &b) { + + //Transform to output color. Standard sRGB is D65, but internal representation is D50 + //Note that it is only at this point that we should have need of clipping color data + + /*float x65 = d65_d50[0][0]*x + d65_d50[0][1]*y + d65_d50[0][2]*z ; + float y65 = d65_d50[1][0]*x + d65_d50[1][1]*y + d65_d50[1][2]*z ; + float z65 = d65_d50[2][0]*x + d65_d50[2][1]*y + d65_d50[2][2]*z ; + + r = sRGB_xyz[0][0]*x65 + sRGB_xyz[0][1]*y65 + sRGB_xyz[0][2]*z65; + g = sRGB_xyz[1][0]*x65 + sRGB_xyz[1][1]*y65 + sRGB_xyz[1][2]*z65; + b = sRGB_xyz[2][0]*x65 + sRGB_xyz[2][1]*y65 + sRGB_xyz[2][2]*z65;*/ + + /*r = sRGBd65_xyz[0][0]*x + sRGBd65_xyz[0][1]*y + sRGBd65_xyz[0][2]*z ; + g = sRGBd65_xyz[1][0]*x + sRGBd65_xyz[1][1]*y + sRGBd65_xyz[1][2]*z ; + b = sRGBd65_xyz[2][0]*x + sRGBd65_xyz[2][1]*y + sRGBd65_xyz[2][2]*z ;*/ + + r = ((sRGB_xyz[0][0]*x + sRGB_xyz[0][1]*y + sRGB_xyz[0][2]*z)) ; + g = ((sRGB_xyz[1][0]*x + sRGB_xyz[1][1]*y + sRGB_xyz[1][2]*z)) ; + b = ((sRGB_xyz[2][0]*x + sRGB_xyz[2][1]*y + sRGB_xyz[2][2]*z)) ; + + } + + + void ImProcFunctions::xyz2rgb (float x, float y, float z, float &r, float &g, float &b, float rgb_xyz[3][3]) { + + //Transform to output color. Standard sRGB is D65, but internal representation is D50 + //Note that it is only at this point that we should have need of clipping color data + + /*float x65 = d65_d50[0][0]*x + d65_d50[0][1]*y + d65_d50[0][2]*z ; + float y65 = d65_d50[1][0]*x + d65_d50[1][1]*y + d65_d50[1][2]*z ; + float z65 = d65_d50[2][0]*x + d65_d50[2][1]*y + d65_d50[2][2]*z ; + + r = sRGB_xyz[0][0]*x65 + sRGB_xyz[0][1]*y65 + sRGB_xyz[0][2]*z65; + g = sRGB_xyz[1][0]*x65 + sRGB_xyz[1][1]*y65 + sRGB_xyz[1][2]*z65; + b = sRGB_xyz[2][0]*x65 + sRGB_xyz[2][1]*y65 + sRGB_xyz[2][2]*z65;*/ + + /*r = sRGBd65_xyz[0][0]*x + sRGBd65_xyz[0][1]*y + sRGBd65_xyz[0][2]*z ; + g = sRGBd65_xyz[1][0]*x + sRGBd65_xyz[1][1]*y + sRGBd65_xyz[1][2]*z ; + b = sRGBd65_xyz[2][0]*x + sRGBd65_xyz[2][1]*y + sRGBd65_xyz[2][2]*z ;*/ + + + + r = ((rgb_xyz[0][0]*x + rgb_xyz[0][1]*y + rgb_xyz[0][2]*z)) ; + g = ((rgb_xyz[1][0]*x + rgb_xyz[1][1]*y + rgb_xyz[1][2]*z)) ; + b = ((rgb_xyz[2][0]*x + rgb_xyz[2][1]*y + rgb_xyz[2][2]*z)) ; + + } + + void ImProcFunctions::Lab2XYZ(float L, float a, float b, float &x, float &y, float &z) { + float fy = (0.00862069 * L) + 0.137932; // (L+16)/116 + float fx = (0.002 * a) + fy; + float fz = fy - (0.005 * b); + + x = 65535*f2xyz(fx)*D50x; + y = 65535*f2xyz(fy); + z = 65535*f2xyz(fz)*D50z; + } + + void ImProcFunctions::XYZ2Lab(float X, float Y, float Z, float &L, float &a, float &b) { + + float fx = (X<65535.0 ? cachef[X] : (327.68*exp(log(X/MAXVAL)/3.0 ))); + float fy = (Y<65535.0 ? cachef[Y] : (327.68*exp(log(Y/MAXVAL)/3.0 ))); + float fz = (Z<65535.0 ? cachef[Z] : (327.68*exp(log(Z/MAXVAL)/3.0 ))); + + L = (116.0 * fy - 5242.88); //5242.88=16.0*327.68; + a = (500.0 * (fx - fy) ); + b = (200.0 * (fy - fz) ); + + } + + void ImProcFunctions::Lab2Yuv(float L, float a, float b, float &Y, float &u, float &v) { + float fy = (0.00862069 * L/327.68) + 0.137932; // (L+16)/116 + float fx = (0.002 * a/327.68) + fy; + float fz = fy - (0.005 * b/327.68); + + float X = 65535.0*f2xyz(fx)*D50x; + Y = 65535.0*f2xyz(fy); + float Z = 65535.0*f2xyz(fz)*D50z; + + u = 4.0*X/(X+15*Y+3*Z)-u0; + v = 9.0*Y/(X+15*Y+3*Z)-v0; + + /*float u0 = 4*D50x/(D50x+15+3*D50z); + float v0 = 9/(D50x+15+3*D50z); + u -= u0; + v -= v0;*/ + + } + + void ImProcFunctions::Yuv2Lab(float Yin, float u, float v, float &L, float &a, float &b, double wp[3][3]) { + + float u1 = u + u0; + float v1 = v + v0; + + float Y = Yin; + float X = (9*u1*Y)/(4*v1*D50x); + float Z = (12 - 3*u1 - 20*v1)*Y/(4*v1*D50z); + + gamutmap(X,Y,Z,wp); + + float fx = (X<65535.0 ? cachef[X] : (327.68*exp(log(X/MAXVAL)/3.0 ))); + float fy = (Y<65535.0 ? cachef[Y] : (327.68*exp(log(Y/MAXVAL)/3.0 ))); + float fz = (Z<65535.0 ? cachef[Z] : (327.68*exp(log(Z/MAXVAL)/3.0 ))); + + L = (116.0 * fy - 5242.88); //5242.88=16.0*327.68; + a = (500.0 * (fx - fy) ); + b = (200.0 * (fy - fz) ); + + } + +#include "gamutbdy.cc" } +#undef eps_max +#undef kappa