diff --git a/rtengine/curves.cc b/rtengine/curves.cc index ef313ffbc..6dcc221d0 100644 --- a/rtengine/curves.cc +++ b/rtengine/curves.cc @@ -38,8 +38,15 @@ namespace rtengine { x = 0; y = 0; ypp = 0; + hash = NULL; + hashSize = 1000; // has to be initiallised to the maximum value } + Curve::~Curve () { + if (hash) + delete [] hash; + } + void Curve::AddPolygons () { if (firstPointIncluded) { @@ -62,12 +69,71 @@ namespace rtengine { poly_y.push_back(y3); } + void Curve::fillHash() { + hash = new unsigned short int[hashSize+2]; + + unsigned int polyIter = 0; + double const increment = 1./hashSize; + double milestone = 0.; + + for (unsigned int i=0; i<(hashSize+1); i++) { + while(poly_x[polyIter] <= milestone) ++polyIter; + hash[i] = polyIter-1; + milestone += increment; + } + hash[hashSize+1] = poly_x.size()-1; + + /* + // Debug output to file + FILE* f = fopen ("hash.txt", "wt"); + for (int i=0; i<(hashSize+2); i++) + fprintf (f, "%d: %d > %.6f, %.6f\n", i, hash[i], poly_x[hash[i]], poly_y[hash[i]]); + fprintf (f, "\nppn: %d\npoly_x: %d\n", ppn, poly_x.size()); + fclose (f); + */ + } + // Wikipedia sRGB: Unlike most other RGB color spaces, the sRGB gamma cannot be expressed as a single numerical value. // The overall gamma is approximately 2.2, consisting of a linear (gamma 1.0) section near black, and a non-linear section elsewhere involving a 2.4 exponent // and a gamma (slope of log output versus log input) changing from 1.0 through about 2.3. const double CurveFactory::sRGBGamma = 2.2; const double CurveFactory::sRGBGammaCurve = 2.4; + void fillCurveArray(DiagonalCurve* diagCurve, LUTf &outCurve, int skip, bool needed) { + if (needed) { + LUTf lutCurve (65536); + + for (int i=0; i<=0xffff; i+= i<0xffff-skip ? skip : 1 ) { + // change to [0,1] range + double val = (double)i / 65535.0; + // apply custom/parametric/NURBS curve, if any + val = diagCurve->getVal (val); + // store result in a temporary array + lutCurve[i] = (val); + } + + // if skip>1, let apply linear interpolation in the skipped points of the curve + if (skip > 1) { + int prev = 0; + for (int i=1; i<=0xffff-skip; i++) { + if (i%skip==0) { + prev+=skip; + continue; + } + lutCurve[i] = ( lutCurve[prev] * (skip - i%skip) + lutCurve[prev+skip] * (i%skip) ) / skip; + } + } + + for (int i=0; i<=0xffff; i++) { + outCurve[i] = (65535.0 * lutCurve[i]); + } + } + else { + for (int i=0; i<=0xffff; i++) { + outCurve[i] = (float)i; + } + } + } //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% void CurveFactory::complexsgnCurve (double saturation, bool satlimit, double satlimthresh, \ @@ -76,118 +142,99 @@ namespace rtengine { //colormult = chroma_scale for Lab manipulations + //----------------------------------------------------- + + bool needed; + DiagonalCurve* dCurve = NULL; + // check if contrast curve is needed - bool needsaturation = (saturation<-0.0001 || saturation>0.0001); - - // curve without contrast - LUTf dacurve (65536); - LUTf dbcurve (65536); + needed = (saturation<-0.0001 || saturation>0.0001); - LUTf dscurve (65536); - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - std::vector satcurvePoints; - satcurvePoints.push_back((double)DCT_NURBS); - if (saturation>0) { - double satslope = (0.5+2*saturation/500.0)/(0.5-2*saturation/500.0); - double scale = (satlimthresh/100.1); - if (!satlimit) scale=100/100.1; - - satcurvePoints.push_back(0); //black point. Value in [0 ; 1] range - satcurvePoints.push_back(0); //black point. Value in [0 ; 1] range - - //if (satlimit) { - satcurvePoints.push_back(0.5-0.5*scale); //toe point - satcurvePoints.push_back(0.5-0.5*scale); //value at toe point - - satcurvePoints.push_back(0.5-(0.5/satslope)*scale); //toe point - satcurvePoints.push_back(0.5-0.5*scale); //value at toe point - - satcurvePoints.push_back(0.5+(0.5/satslope)*scale); //shoulder point - satcurvePoints.push_back(0.5+0.5*scale); //value at shoulder point - - satcurvePoints.push_back(0.5+0.5*scale); //shoulder point - satcurvePoints.push_back(0.5+0.5*scale); //value at shoulder point - /*} else { - satcurvePoints.push_back(0.25+saturation/500.0); //toe point - satcurvePoints.push_back(0.25-saturation/500.0); //value at toe point - - satcurvePoints.push_back(0.75-saturation/500.0); //shoulder point - satcurvePoints.push_back(0.75+saturation/500.0); //value at shoulder point - }*/ - - satcurvePoints.push_back(1); // white point - satcurvePoints.push_back(1); // value at white point - } else { - satcurvePoints.push_back(0); - satcurvePoints.push_back(-(saturation/200.0)); - - satcurvePoints.push_back(1); - satcurvePoints.push_back(1+saturation/200.0); + // Filling the curve if needed + if (needed) { + + //%%%%%%%%%%%%%%%%% Saturation curve's control points %%%%%%%%%%%%%%%%% + std::vector satcurvePoints; + satcurvePoints.push_back((double)DCT_NURBS); + if (saturation>0) { + double satslope = (0.5+2*saturation/500.0)/(0.5-2*saturation/500.0); + double scale = (satlimthresh/100.1); + if (!satlimit) scale=100/100.1; + + satcurvePoints.push_back(0); //black point. Value in [0 ; 1] range + satcurvePoints.push_back(0); //black point. Value in [0 ; 1] range + + //if (satlimit) { + satcurvePoints.push_back(0.5-0.5*scale); //toe point + satcurvePoints.push_back(0.5-0.5*scale); //value at toe point + + satcurvePoints.push_back(0.5-(0.5/satslope)*scale); //toe point + satcurvePoints.push_back(0.5-0.5*scale); //value at toe point + + satcurvePoints.push_back(0.5+(0.5/satslope)*scale); //shoulder point + satcurvePoints.push_back(0.5+0.5*scale); //value at shoulder point + + satcurvePoints.push_back(0.5+0.5*scale); //shoulder point + satcurvePoints.push_back(0.5+0.5*scale); //value at shoulder point + /*} else { + satcurvePoints.push_back(0.25+saturation/500.0); //toe point + satcurvePoints.push_back(0.25-saturation/500.0); //value at toe point + + satcurvePoints.push_back(0.75-saturation/500.0); //shoulder point + satcurvePoints.push_back(0.75+saturation/500.0); //value at shoulder point + }*/ + + satcurvePoints.push_back(1); // white point + satcurvePoints.push_back(1); // value at white point + } else { + satcurvePoints.push_back(0); + satcurvePoints.push_back(-(saturation/200.0)); + + satcurvePoints.push_back(1); + satcurvePoints.push_back(1+saturation/200.0); + } + dCurve = new DiagonalCurve (satcurvePoints, CURVES_MIN_POLY_POINTS/skip); + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + fillCurveArray(dCurve, satCurve, skip, needed); + + delete dCurve; + dCurve = NULL; } - DiagonalCurve* satcurve = new DiagonalCurve (satcurvePoints, CURVES_MIN_POLY_POINTS/skip); // Actually, CURVES_MIN_POLY_POINTS = 1000, - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - + else { + fillCurveArray(NULL, satCurve, skip, needed); + } + + //----------------------------------------------------- + + needed = false; // create a curve if needed - DiagonalCurve* tacurve = NULL; - if (acurvePoints.size()>0 && acurvePoints[0]!=0) - tacurve = new DiagonalCurve (acurvePoints, CURVES_MIN_POLY_POINTS/skip); - DiagonalCurve* tbcurve = NULL; - if (bcurvePoints.size()>0 && bcurvePoints[0]!=0) - tbcurve = new DiagonalCurve (bcurvePoints, CURVES_MIN_POLY_POINTS/skip); - - for (int i=0; i<=0xffff; i+= i<0xffff-skip ? skip : 1 ) { - - // change to [0,1] range - double aval = (double)i / 65535.0; - double bval = (double)i / 65535.0; - double sval = (double)i / 65535.0; - - - // apply saturation curve - if (needsaturation) - sval = satcurve->getVal (sval); - - // apply custom/parametric/NURBS curve, if any - if (tacurve) { - aval = tacurve->getVal (aval); - } - // apply custom/parametric/NURBS curve, if any - if (tbcurve) { - bval = tbcurve->getVal (bval); - } - - // store result in a temporary array - dacurve[i] = (aval); - dbcurve[i] = (bval); - dscurve[i] = (sval); + if (acurvePoints.size()>0 && acurvePoints[0]!=0) { + dCurve = new DiagonalCurve (acurvePoints, CURVES_MIN_POLY_POINTS/skip); + if (dCurve && !dCurve->isIdentity()) + needed = true; } - - delete tacurve; - delete tbcurve; - - // if skip>1, let apply linear interpolation in the skipped points of the curve - int prev = 0; - for (int i=1; i<=0xffff-skip; i++) { - if (i%skip==0) { - prev+=skip; - continue; - } - dacurve[i] = ( dacurve[prev] * (skip - i%skip) + dacurve[prev+skip] * (i%skip) ) / skip; - dbcurve[i] = ( dbcurve[prev] * (skip - i%skip) + dbcurve[prev+skip] * (i%skip) ) / skip; - dscurve[i] = ( dscurve[prev] * (skip - i%skip) + dscurve[prev+skip] * (i%skip) ) / skip; - } - - for (int i=0; i<=0xffff; i++) { - aoutCurve[i] = (65535.0 * dacurve[i]); - boutCurve[i] = (65535.0 * dbcurve[i]); - satCurve[i] = (65535.0 * dscurve[i]); + fillCurveArray(dCurve, aoutCurve, skip, needed); + if (dCurve) { + delete dCurve; + dCurve = NULL; } - delete satcurve; + //----------------------------------------------------- + + needed = false; + if (bcurvePoints.size()>0 && bcurvePoints[0]!=0) { + dCurve = new DiagonalCurve (bcurvePoints, CURVES_MIN_POLY_POINTS/skip); + if (dCurve && !dCurve->isIdentity()) + needed = true; + } + fillCurveArray(dCurve, boutCurve, skip, needed); + if (dCurve) { + delete dCurve; + dCurve = NULL; + } } - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -223,17 +270,9 @@ namespace rtengine { // curve without contrast LUTf dcurve(0x10000); - // check if contrast curve is needed - bool needcontrast = contr>0.00001 || contr<-0.00001; - // check if inverse gamma is needed at the end bool needigamma = igamma_ && gamma_>0; - // create a curve if needed - DiagonalCurve* tcurve = NULL; - if (curvePoints.size()>0 && curvePoints[0]!=0) - tcurve = new DiagonalCurve (curvePoints, CURVES_MIN_POLY_POINTS/skip); - // clear array that stores histogram valid before applying the custom curve outBeforeCCurveHistogram.clear(); @@ -241,31 +280,37 @@ namespace rtengine { // tone curve base. a: slope (from exp.comp.), b: black, def_mul: max. x value (can be>1), hr,sr: highlight,shadow recovery //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - std::vector brightcurvePoints; - brightcurvePoints.push_back((double)DCT_NURBS); - - brightcurvePoints.push_back(0); //black point. Value in [0 ; 1] range - brightcurvePoints.push_back(0); //black point. Value in [0 ; 1] range - - if(br>0) { - brightcurvePoints.push_back(0.1); //toe point - brightcurvePoints.push_back(0.1+br/150.0); //value at toe point + DiagonalCurve* brightcurve = NULL; + + // check if brightness curve is needed + if (br>0.00001 || br<-0.00001) { + + std::vector brightcurvePoints; + brightcurvePoints.push_back((double)DCT_NURBS); + + brightcurvePoints.push_back(0.); //black point. Value in [0 ; 1] range + brightcurvePoints.push_back(0.); //black point. Value in [0 ; 1] range - brightcurvePoints.push_back(0.7); //shoulder point - brightcurvePoints.push_back(MIN(1.0,0.7+br/300.0)); //value at shoulder point - } else { - brightcurvePoints.push_back(0.1-br/150.0); //toe point - brightcurvePoints.push_back(0.1); //value at toe point + if(br>0) { + brightcurvePoints.push_back(0.1); //toe point + brightcurvePoints.push_back(0.1+br/150.0); //value at toe point + + brightcurvePoints.push_back(0.7); //shoulder point + brightcurvePoints.push_back(MIN(1.0,0.7+br/300.0)); //value at shoulder point + } else { + brightcurvePoints.push_back(0.1-br/150.0); //toe point + brightcurvePoints.push_back(0.1); //value at toe point + + brightcurvePoints.push_back(MIN(1.0,0.7-br/300.0)); //shoulder point + brightcurvePoints.push_back(0.7); //value at shoulder point + } + brightcurvePoints.push_back(1.); // white point + brightcurvePoints.push_back(1.); // value at white point - brightcurvePoints.push_back(MIN(1.0,0.7-br/300.0)); //shoulder point - brightcurvePoints.push_back(0.7); //value at shoulder point + brightcurve = new DiagonalCurve (brightcurvePoints, CURVES_MIN_POLY_POINTS/skip); } - brightcurvePoints.push_back(1); // white point - brightcurvePoints.push_back(1); // value at white point - - DiagonalCurve* brightcurve = new DiagonalCurve (brightcurvePoints, CURVES_MIN_POLY_POINTS/skip); // Actually, CURVES_MIN_POLY_POINTS = 1000, //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - + float exp_scale = a; float scale = 65536.0; float comp = (MAX(0,ecomp) + 1.0)*hlcompr/100.0; @@ -307,21 +352,27 @@ namespace rtengine { val = (double)i / 65535.0; // gamma correction - if (gamma_>1) + if (gamma_>1) val = gamma (val, gamma_, start, slope, mul, add); // apply brightness curve - val = brightcurve->getVal (val); - + if (brightcurve) + val = brightcurve->getVal (val); // TODO: getVal(double) is very slow! Optimize with a LUTf + // store result in a temporary array dcurve[i] = CLIPD(val); } - + + if (brightcurve) + delete brightcurve; + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - if (needcontrast) { + + // check if contrast curve is needed + if (contr>0.00001 || contr<-0.00001) { + // compute mean luminance of the image with the curve applied int sum = 0; float avg = 0; @@ -352,49 +403,64 @@ namespace rtengine { contrastcurvePoints.push_back(1); // white point contrastcurvePoints.push_back(1); // value at white point - DiagonalCurve* contrastcurve = new DiagonalCurve (contrastcurvePoints, CURVES_MIN_POLY_POINTS/skip); // Actually, CURVES_MIN_POLY_POINTS = 1000, + DiagonalCurve* contrastcurve = new DiagonalCurve (contrastcurvePoints, CURVES_MIN_POLY_POINTS/skip); //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% // apply contrast enhancement for (int i=0; i<=0xffff; i++) { dcurve[i] = contrastcurve->getVal (dcurve[i]); } + delete contrastcurve; } + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + // create a curve if needed + bool histNeeded = false; + DiagonalCurve* tcurve = NULL; + if (curvePoints.size()>0 && curvePoints[0]!=0) { + tcurve = new DiagonalCurve (curvePoints, CURVES_MIN_POLY_POINTS/skip); + if (outBeforeCCurveHistogram && histogramCropped) + histNeeded = true; + } + if (tcurve && tcurve->isIdentity()) { + delete tcurve; + tcurve = NULL; + } + for (int i=0; i<=0xffff; i++) { float val; - + + if (histNeeded) { + float fi=i; + float hval = dcurve[shCurve[hlCurve[i]*fi]*fi]; + //if (needigamma) + // hval = igamma2 (hval); + int hi = (int)(255.0*(hval)); + outBeforeCCurveHistogram[hi] += histogramCropped[i] ; + } + // apply custom/parametric/NURBS curve, if any if (tcurve) { - if (outBeforeCCurveHistogram && histogramCropped) { - float fi=i; - float hval = dcurve[shCurve[hlCurve[i]*fi]*fi]; - //if (needigamma) - // hval = igamma2 (hval); - int hi = (int)(255.0*(hval)); - outBeforeCCurveHistogram[hi] += histogramCropped[i] ; - } - val = tcurve->getVal (dcurve[i]); + val = tcurve->getVal (dcurve[i]); // TODO: getVal(double) is very slow! Optimize with a LUTf } else { val = (dcurve[i]); } - + // if inverse gamma is needed, do it (standard sRGB inverse gamma is applied) if (needigamma) val = igamma (val, gamma_, start, slope, mul, add); - + outCurve[i] = (65535.0 * val); } - - - delete tcurve; - delete brightcurve; + + if (tcurve) + delete tcurve; + /*if (outBeforeCCurveHistogram) { for (int i=0; i<256; i++) printf("i= %d bchist= %d \n",i,outBeforeCCurveHistogram[i]); }*/ - } @@ -409,15 +475,6 @@ namespace rtengine { // curve without contrast LUTf dcurve(65536,0); - // check if contrast curve is needed - bool needcontrast = contr>0.00001 || contr<-0.00001; - - // create a curve if needed - DiagonalCurve* tcurve = NULL; - if (curvePoints.size()>0 && curvePoints[0]!=0) - tcurve = new DiagonalCurve (curvePoints, CURVES_MIN_POLY_POINTS/skip); - - // clear array that stores histogram valid before applying the custom curve if (outBeforeCCurveHistogram) outBeforeCCurveHistogram.clear(); @@ -426,48 +483,62 @@ namespace rtengine { // tone curve base. a: slope (from exp.comp.), b: black, def_mul: max. x value (can be>1), hr,sr: highlight,shadow recovery //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - std::vector brightcurvePoints; - brightcurvePoints.push_back((double)((CurveType)DCT_NURBS)); - - brightcurvePoints.push_back(0); // black point. Value in [0 ; 1] range - brightcurvePoints.push_back(0); // black point. Value in [0 ; 1] range - - if (br>0) { - brightcurvePoints.push_back(0.1); // toe point - brightcurvePoints.push_back(0.1+br/150.0); //value at toe point + // check if brightness curve is needed + if (br>0.00001 || br<-0.00001) { + + std::vector brightcurvePoints; + brightcurvePoints.push_back((double)((CurveType)DCT_NURBS)); + + brightcurvePoints.push_back(0.); // black point. Value in [0 ; 1] range + brightcurvePoints.push_back(0.); // black point. Value in [0 ; 1] range - brightcurvePoints.push_back(0.7); // shoulder point - brightcurvePoints.push_back(MIN(1.0,0.7+br/300.0)); //value at shoulder point - } else { - brightcurvePoints.push_back(0.1-br/150.0); // toe point - brightcurvePoints.push_back(0.1); // value at toe point + if (br>0) { + brightcurvePoints.push_back(0.1); // toe point + brightcurvePoints.push_back(0.1+br/150.0); //value at toe point + + brightcurvePoints.push_back(0.7); // shoulder point + brightcurvePoints.push_back(MIN(1.0,0.7+br/300.0)); //value at shoulder point + } else { + brightcurvePoints.push_back(0.1-br/150.0); // toe point + brightcurvePoints.push_back(0.1); // value at toe point + + brightcurvePoints.push_back(MIN(1.0,0.7-br/300.0)); // shoulder point + brightcurvePoints.push_back(0.7); // value at shoulder point + } + brightcurvePoints.push_back(1.); // white point + brightcurvePoints.push_back(1.); // value at white point - brightcurvePoints.push_back(MIN(1.0,0.7-br/300.0)); // shoulder point - brightcurvePoints.push_back(0.7); // value at shoulder point + DiagonalCurve* brightcurve = new DiagonalCurve (brightcurvePoints, CURVES_MIN_POLY_POINTS/skip); + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + // Applying brightness curve + for (int i=0; i<32768; i++) { // L values range up to 32767, higher values are for highlight overflow + + // change to [0,1] range + float val = (float)i / 32767.0; + + // apply brightness curve + val = brightcurve->getVal (val); + + // store result in a temporary array + dcurve[i] = CLIPD(val); + } + delete brightcurve; } - brightcurvePoints.push_back(1); // white point - brightcurvePoints.push_back(1); // value at white point - - DiagonalCurve* brightcurve = new DiagonalCurve (brightcurvePoints, CURVES_MIN_POLY_POINTS/skip); // Actually, CURVES_MIN_POLY_POINTS = 1000, - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - for (int i=0; i<32768; i++) { // L values range up to 32767, higher values are for highlight overflow - - // change to [0,1] range - float val = (float)i / 32767.0; - - // apply brightness curve - val = brightcurve->getVal (val); - - // store result in a temporary array - dcurve[i] = CLIPD(val); + else { + for (int i=0; i<32768; i++) { // L values range up to 32767, higher values are for highlight overflow + // set the identity curve in the temporary array + dcurve[i] = (float)i / 32767.0; + } } //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - if (needcontrast) { + // check if contrast curve is needed + if (contr>0.00001 || contr<-0.00001) { + // compute mean luminance of the image with the curve applied int sum = 0; float avg = 0; @@ -485,8 +556,8 @@ namespace rtengine { std::vector contrastcurvePoints; contrastcurvePoints.push_back((double)((CurveType)DCT_NURBS)); - contrastcurvePoints.push_back(0); // black point. Value in [0 ; 1] range - contrastcurvePoints.push_back(0); // black point. Value in [0 ; 1] range + contrastcurvePoints.push_back(0.); // black point. Value in [0 ; 1] range + contrastcurvePoints.push_back(0.); // black point. Value in [0 ; 1] range contrastcurvePoints.push_back(avg-avg*(0.6-contr/250.0)); // toe point contrastcurvePoints.push_back(avg-avg*(0.6+contr/250.0)); // value at toe point @@ -494,10 +565,10 @@ namespace rtengine { contrastcurvePoints.push_back(avg+(1-avg)*(0.6-contr/250.0)); // shoulder point contrastcurvePoints.push_back(avg+(1-avg)*(0.6+contr/250.0)); // value at shoulder point - contrastcurvePoints.push_back(1); // white point - contrastcurvePoints.push_back(1); // value at white point + contrastcurvePoints.push_back(1.); // white point + contrastcurvePoints.push_back(1.); // value at white point - DiagonalCurve* contrastcurve = new DiagonalCurve (contrastcurvePoints, CURVES_MIN_POLY_POINTS/skip); // Actually, CURVES_MIN_POLY_POINTS = 1000, + DiagonalCurve* contrastcurve = new DiagonalCurve (contrastcurvePoints, CURVES_MIN_POLY_POINTS/skip); //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% // apply contrast enhancement @@ -508,39 +579,64 @@ namespace rtengine { delete contrastcurve; } //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - for (int i=0; i<32768; i++) { // L values go up to 32767, last stop is for highlight overflow - float val; - - // apply custom/parametric/NURBS curve, if any - if (tcurve) { - if (outBeforeCCurveHistogram) { + + // create a curve if needed + DiagonalCurve* tcurve = NULL; + bool histNeeded = false; + if (curvePoints.size()>0 && curvePoints[0]!=0) { + tcurve = new DiagonalCurve (curvePoints, CURVES_MIN_POLY_POINTS/skip); + if (outBeforeCCurveHistogram && histogramCropped) + histNeeded = true; + } + if (tcurve && tcurve->isIdentity()) { + delete tcurve; + tcurve = NULL; + } + + if (tcurve) { + // L values go up to 32767, last stop is for highlight overflow + for (int i=0; i<32768; i++) { + float val; + + if (histNeeded) { float hval = dcurve[i]; int hi = (int)(255.0*CLIPD(hval)); outBeforeCCurveHistogram[hi]+=histogramCropped[i] ; } + + // apply custom/parametric/NURBS curve, if any val = tcurve->getVal (dcurve[i]); - } else { - val = (dcurve[i]); + + outCurve[i] = (32767.0 * val); } - - outCurve[i] = (32767.0 * val); } - for (int i=32768; i<65536; i++) outCurve[i]=i; - - - delete tcurve; - delete brightcurve; + else { + // Skip the slow getval method if no curve is used (or an identity curve) + // L values go up to 32767, last stop is for highlight overflow + for (int i=0; i<32768; i++) { + if (histNeeded) { + float hval = dcurve[i]; + int hi = (int)(255.0*CLIPD(hval)); + outBeforeCCurveHistogram[hi]+=histogramCropped[i] ; + } + + outCurve[i] = 32767.0*dcurve[i]; + } + } + for (int i=32768; i<65536; i++) outCurve[i]=(float)i; + + if (tcurve) + delete tcurve; + /*if (outBeforeCCurveHistogram) { for (int i=0; i<256; i++) printf("i= %d bchist= %d \n",i,outBeforeCCurveHistogram[i]); }*/ - } - - + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - + + LUTf CurveFactory::gammatab; LUTf CurveFactory::igammatab_srgb; diff --git a/rtengine/curves.h b/rtengine/curves.h index 017569bd8..5ebcf4af6 100644 --- a/rtengine/curves.h +++ b/rtengine/curves.h @@ -35,7 +35,6 @@ #define CLIPI(a) ((a)>0?((a)<65534?(a):65534):0) - namespace rtengine { class CurveFactory { @@ -211,7 +210,10 @@ class Curve { double* y; std::vector poly_x; // X points of the faceted curve std::vector poly_y; // Y points of the faceted curve - double* ypp; + unsigned short int* hash; + unsigned int hashSize; // hash table's size, between [10, 100, 1000] + + double* ypp; // Fields for the elementary curve polygonisation double x1, y1, x2, y2, x3, y3; @@ -225,12 +227,16 @@ class Curve { static inline double p10 (double x, double prot) { return x<=0.5 ? CurveFactory::cupper (x*2, 2.0, prot)/2.0 : 0.5 + CurveFactory::clower ((x-0.5)*2, 2.0, prot)/2.0; } static inline double pfull (double x, double prot, double sh, double hl) { return (1-sh)*(1-hl)*p00(x,prot) + sh*hl*p11(x,prot) + (1-sh)*hl*p01(x,prot) + sh*(1-hl)*p10(x,prot); } + void fillHash(); + public: Curve (); + ~Curve (); void AddPolygons (); virtual double getVal (double t) = 0; virtual void getVal (const std::vector& t, std::vector& res) = 0; + virtual bool isIdentity () = 0; }; class DiagonalCurve : public Curve { @@ -238,6 +244,10 @@ class DiagonalCurve : public Curve { protected: DiagonalCurveType kind; + unsigned int minSearch; // a effacer!!! + unsigned int maxSearch; // a effacer!!! + unsigned int searchArray[21]; // a effacer!!! + void spline_cubic_set (); void NURBS_set (); @@ -245,8 +255,9 @@ class DiagonalCurve : public Curve { DiagonalCurve (const std::vector& points, int ppn=CURVES_MIN_POLY_POINTS); ~DiagonalCurve (); - double getVal (double t); - void getVal (const std::vector& t, std::vector& res); + double getVal (double t); + void getVal (const std::vector& t, std::vector& res); + bool isIdentity () { return kind==DCT_Empty; }; }; class FlatCurve : public Curve { @@ -264,8 +275,9 @@ class FlatCurve : public Curve { FlatCurve (const std::vector& points, bool isPeriodic = true, int ppn=CURVES_MIN_POLY_POINTS); ~FlatCurve (); - double getVal (double t); - void getVal (const std::vector& t, std::vector& res); + double getVal (double t); + void getVal (const std::vector& t, std::vector& res); + bool isIdentity () { return kind==FCT_Empty; }; }; } diff --git a/rtengine/diagonalcurves.cc b/rtengine/diagonalcurves.cc index f07b8f0a7..5f1632c84 100644 --- a/rtengine/diagonalcurves.cc +++ b/rtengine/diagonalcurves.cc @@ -30,7 +30,11 @@ namespace rtengine { DiagonalCurve::DiagonalCurve (const std::vector& p, int poly_pn) { - ppn = poly_pn; + ppn = poly_pn; + bool identity = true; + + if (ppn < 500) hashSize = 100; // Arbitrary cut-off value + if (ppn < 50) hashSize = 10; // Arbitrary cut-off value if (p.size()<3) { kind = DCT_Empty; @@ -45,17 +49,28 @@ DiagonalCurve::DiagonalCurve (const std::vector& p, int poly_pn) { for (int i=0; i 2) + spline_cubic_set (); + else if (kind==DCT_NURBS && N > 2) { + NURBS_set (); + fillHash(); + } + else kind=DCT_Linear; } - if (kind==DCT_Spline) - spline_cubic_set (); - else if (kind==DCT_NURBS && N > 2) - NURBS_set (); - else kind=DCT_Linear; } else if (kind==DCT_Parametric) { - if (p.size()!=8 && p.size()!=9) - kind = DCT_Empty; - else { + if ((p.size()==8 || p.size()==9) && (p.at(4)!=0.0f || p.at(5)!=0.0f || p.at(6)!=0.0f || p.at(7)!=0.0f)) { + identity = false; + x = new double[9]; for (int i=0; i<4; i++) x[i] = p[i]; @@ -67,6 +82,9 @@ DiagonalCurve::DiagonalCurve (const std::vector& p, int poly_pn) { x[8] = p[8]/100.0; } } + if (identity) { + kind = DCT_Empty; + } } } @@ -82,7 +100,7 @@ DiagonalCurve::~DiagonalCurve () { void DiagonalCurve::spline_cubic_set () { double* u = new double[N-1]; - delete [] ypp; + delete [] ypp; // TODO: why do we delete ypp here since it should not be allocated yet? ypp = new double [N]; ypp[0] = u[0] = 0.0; /* set lower boundary condition to "natural" */ @@ -113,7 +131,7 @@ void DiagonalCurve::NURBS_set () { double total_length=0.; // Create the list of Bezier sub-curves - // NURBS_set is called if N > 2 only + // NURBS_set is called if N > 2 and non identity only int j = 0; int k = 0; @@ -122,76 +140,84 @@ void DiagonalCurve::NURBS_set () { double dx; double dy; - // first point (on the curve) - if (!i) { - sc_x[j] = x[i]; - sc_y[j++] = y[i++]; - } - else { - sc_x[j] = (x[i-1] + x[i]) / 2.; - sc_y[j++] = (y[i-1] + y[i]) / 2.; - } + // first point (on the curve) + if (!i) { + sc_x[j] = x[i]; + sc_y[j++] = y[i++]; + } + else { + sc_x[j] = (x[i-1] + x[i]) / 2.; + sc_y[j++] = (y[i-1] + y[i]) / 2.; + } - // second point (control point) - sc_x[j] = x[i]; - sc_y[j] = y[i++]; + // second point (control point) + sc_x[j] = x[i]; + sc_y[j] = y[i++]; - dx = sc_x[j] - sc_x[j-1]; - dy = sc_y[j] - sc_y[j-1]; - length = sqrt(dx*dx + dy*dy); - j++; + dx = sc_x[j] - sc_x[j-1]; + dy = sc_y[j] - sc_y[j-1]; + length = sqrt(dx*dx + dy*dy); + j++; - // third point (on the curve) - if (i==N-1) { - sc_x[j] = x[i]; - sc_y[j] = y[i]; - } - else { - sc_x[j] = (x[i-1] + x[i]) / 2.; - sc_y[j] = (y[i-1] + y[i]) / 2.; - } - dx = sc_x[j] - sc_x[j-1]; - dy = sc_y[j] - sc_y[j-1]; - length += sqrt(dx*dx + dy*dy); - j++; + // third point (on the curve) + if (i==N-1) { + sc_x[j] = x[i]; + sc_y[j] = y[i]; + } + else { + sc_x[j] = (x[i-1] + x[i]) / 2.; + sc_y[j] = (y[i-1] + y[i]) / 2.; + } + dx = sc_x[j] - sc_x[j-1]; + dy = sc_y[j] - sc_y[j-1]; + length += sqrt(dx*dx + dy*dy); + j++; - // Storing the length of all sub-curves and the total length (to have a better distribution - // of the points along the curve) - sc_length[k++] = length; - total_length += length; + // Storing the length of all sub-curves and the total length (to have a better distribution + // of the points along the curve) + sc_length[k++] = length; + total_length += length; } poly_x.clear(); - poly_y.clear(); - unsigned int sc_xsize=j-1; + poly_y.clear(); + unsigned int sc_xsize=j-1; j = 0; + + // adding the initial horizontal segment, if any + if (x[0] > 0.) { + poly_x.push_back(0.); + poly_y.push_back(y[0]); + } + + // adding the initial horizontal segment, if any // create the polyline with the number of points adapted to the X range of the sub-curve for (unsigned int i=0; i < sc_xsize /*sc_x.size()*/; i+=3) { - // TODO: Speeding-up the interface by caching the polyline, instead of rebuilding it at each action on sliders !!! - nbr_points = (int)(((double)(ppn+N-2) * sc_length[i/3] )/ total_length); - if (nbr_points<0){ - for(int it=0;it < sc_x.size(); it+=3) printf("sc_length[%d/3]=%f \n",it,sc_length[it/3]); - printf("NURBS: error detected!\n i=%d nbr_points=%d ppn=%d N=%d sc_length[i/3]=%f total_length=%f",i,nbr_points,ppn,N,sc_length[i/3],total_length); - exit(0); - } - // increment along the curve, not along the X axis - increment = 1.0 / (double)(nbr_points-1); - x1 = sc_x[i]; y1 = sc_y[i]; - x2 = sc_x[i+1]; y2 = sc_y[i+1]; - x3 = sc_x[i+2]; y3 = sc_y[i+2]; - firstPointIncluded = !i; - AddPolygons (); + // TODO: Speeding-up the interface by caching the polyline, instead of rebuilding it at each action on sliders !!! + nbr_points = (int)(((double)(ppn+N-2) * sc_length[i/3] )/ total_length); + if (nbr_points<0){ + for(int it=0;it < sc_x.size(); it+=3) printf("sc_length[%d/3]=%f \n",it,sc_length[it/3]); + printf("NURBS diagonal curve: error detected!\n i=%d nbr_points=%d ppn=%d N=%d sc_length[i/3]=%f total_length=%f",i,nbr_points,ppn,N,sc_length[i/3],total_length); + exit(0); + } + // increment along the curve, not along the X axis + increment = 1.0 / (double)(nbr_points-1); + x1 = sc_x[i]; y1 = sc_y[i]; + x2 = sc_x[i+1]; y2 = sc_y[i+1]; + x3 = sc_x[i+2]; y3 = sc_y[i+2]; + firstPointIncluded = !i; + AddPolygons (); } + + // adding the final horizontal segment, always (see under) + poly_x.push_back(1.1); // 1.1 is a hack for optimization purpose of the getVal method (the last value has to be beyond the normal range) + poly_y.push_back(y[N-1]); } double DiagonalCurve::getVal (double t) { switch (kind) { - case DCT_Empty : - return t; - break; - case DCT_Parametric : { if (t<=1e-14) return 0.0; @@ -222,7 +248,7 @@ double DiagonalCurve::getVal (double t) { } case DCT_Linear : case DCT_Spline : { - // values under and over the first and last point + // values under and over the first and last point if (t>x[N-1]) return y[N-1]; else if (tx[N-1]) - return y[N-1]; - else if (t 1){ - int k = (k_hi + k_lo) / 2; - if (poly_x[k] > t) - k_hi = k; - else - k_lo = k; - } + if (i > (hashSize+1)) { + //printf("\nOVERFLOW: hash #%d is used while seeking for value %.8f, corresponding polygon's point #%d (out of %d point) x value: %.8f\n\n", i, t, hash[i], poly_x.size(), poly_x[hash[i]]); + printf("\nOVERFLOW: hash #%d is used while seeking for value %.8f\n\n", i, t); + return t; + } - double h = poly_x[k_hi] - poly_x[k_lo]; - return poly_y[k_lo] + (t - poly_x[k_lo]) * ( poly_y[k_hi] - poly_y[k_lo] ) / h; - break; + unsigned int k_lo = 0; + unsigned int k_hi = 0; + + k_lo = hash[i]; + k_hi = hash[i+1]; + + // do a binary search for the right interval : + while (k_hi - k_lo > 1){ + unsigned int k = (k_hi + k_lo) / 2; + if (poly_x[k] > t) + k_hi = k; + else + k_lo = k; + } + if (k_lo == k_hi) + k_hi = k_lo+1; + + double dx = poly_x[k_hi] - poly_x[k_lo]; + double dy = poly_y[k_hi] - poly_y[k_lo]; + return poly_y[k_lo] + (t - poly_x[k_lo]) * ( dy ) / dx; + break; } + case DCT_Empty : default: - // all other (unknown) kind - return t; + // all other (unknown) kind + return t; } } void DiagonalCurve::getVal (const std::vector& t, std::vector& res) { -// TODO!!!! can be made much faster!!! Binary search of getVal(double) at each point can be avoided - res.resize (t.size()); for (unsigned int i=0; i& p, bool isPeriodic, int poly_pn kind = FCT_Empty; periodic = isPeriodic; + bool identity = true; + if (p.size()>4) { kind = (FlatCurveType)p[0]; if (kind==FCT_MinMaxCPoints) { - int oneMorePoint = periodic ? 1:0; + int oneMorePoint = periodic ? 1:0; N = (p.size()-1)/4; x = new double[N+oneMorePoint]; y = new double[N+oneMorePoint]; @@ -53,22 +55,28 @@ FlatCurve::FlatCurve (const std::vector& p, bool isPeriodic, int poly_pn y[i] = p[ix++]; leftTangent[i] = p[ix++]; rightTangent[i] = p[ix++]; + if (y[i] != 0.5f) + identity = false; } // The first point is copied to the end of the point list, to handle the curve periodicity if (periodic) { - x[N] = p[1]+1.0; - y[N] = p[2]; - leftTangent[N] = p[3]; - rightTangent[N] = p[4]; + x[N] = p[1]+1.0; + y[N] = p[2]; + leftTangent[N] = p[3]; + rightTangent[N] = p[4]; } else { - N--; + N--; } - if (N > 0+(periodic?1:0) ) + if (!identity && N > 0+(periodic?1:0) ) { CtrlPoints_set (); + fillHash(); + } } /*else if (kind==FCT_Parametric) { }*/ + if (identity) + kind = FCT_Empty; } } @@ -275,11 +283,11 @@ void FlatCurve::CtrlPoints_set () { } // increment along the curve, not along the X axis increment = 1.0 / (double)(nbr_points-1); - x1 = sc_x[j]; y1 = sc_y[j++]; + x1 = sc_x[j]; y1 = sc_y[j++]; x2 = sc_x[j]; y2 = sc_y[j++]; x3 = sc_x[j]; y3 = sc_y[j++]; AddPolygons (); - } + } } // adding an final horizontal line if necessary @@ -306,18 +314,16 @@ double FlatCurve::getVal (double t) { switch (kind) { case FCT_MinMaxCPoints : { - /* To be updated if we have to handle non periodic flat curves - // values under and over the first and last point + /* To be updated if we have to handle non periodic flat curves + // values under and over the first and last point if (t>x[N-1]) return y[N-1]; else if (t& t, std::vector& res) { -// TODO!!!! can be made much faster!!! Binary search of getVal(double) at each point can be avoided - res.resize (t.size()); for (unsigned int i=0; i