"Fossies" - the Fresh Open Source Software Archive

Member "rawtherapee-5.7/rtengine/diagonalcurves.cc" (10 Sep 2019, 16870 Bytes) of package /linux/misc/rawtherapee-5.7.tar.xz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "diagonalcurves.cc" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 5.6_vs_5.7.

    1 /*
    2  *  This file is part of RawTherapee.
    3  *
    4  *  Copyright (c) 2004-2010 Gabor Horvath <hgabor@rawtherapee.com>
    5  *
    6  *  RawTherapee is free software: you can redistribute it and/or modify
    7  *  it under the terms of the GNU General Public License as published by
    8  *  the Free Software Foundation, either version 3 of the License, or
    9  *  (at your option) any later version.
   10  *
   11  *  RawTherapee is distributed in the hope that it will be useful,
   12  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
   13  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   14  *  GNU General Public License for more details.
   15  *
   16  *  You should have received a copy of the GNU General Public License
   17  *  along with RawTherapee.  If not, see <https://www.gnu.org/licenses/>.
   18  */
   19 #include <glib.h>
   20 #include <glib/gstdio.h>
   21 #include "curves.h"
   22 #include <cmath>
   23 #include <vector>
   24 #include "mytime.h"
   25 #include <cstring>
   26 
   27 #define CLIPD(a) ((a)>0.0?((a)<1.0?(a):1.0):0.0)
   28 
   29 namespace rtengine
   30 {
   31 
   32 DiagonalCurve::DiagonalCurve (const std::vector<double>& p, int poly_pn)
   33 {
   34 
   35     ppn = poly_pn > 65500 ? 65500 : poly_pn;
   36 
   37     if (ppn < 500) {
   38         hashSize = 100;    // Arbitrary cut-off value, but multiple of 10
   39     }
   40 
   41     if (ppn < 50) {
   42         hashSize = 10;    // Arbitrary cut-off value, but multiple of 10
   43     }
   44 
   45     if (p.size() < 3) {
   46         kind = DCT_Empty;
   47     } else {
   48         bool identity = true;
   49         kind = (DiagonalCurveType)p[0];
   50 
   51         if (kind == DCT_Linear || kind == DCT_Spline || kind == DCT_NURBS || kind == DCT_CatumullRom) {
   52             N = (p.size() - 1) / 2;
   53             x = new double[N];
   54             y = new double[N];
   55             int ix = 1;
   56 
   57             for (int i = 0; i < N; i++) {
   58                 x[i] = p[ix++];
   59                 y[i] = p[ix++];
   60 
   61                 if (std::fabs(x[i] - y[i]) >= 0.000009) {
   62                     // the smallest possible difference between x and y curve point values is ~ 0.00001
   63                     // checking against >= 0.000009 is a bit saver than checking against >= 0.00001
   64                     identity = false;
   65                 }
   66             }
   67 
   68             if (x[0] != 0.0f || x[N - 1] != 1.0f)
   69                 // Special (and very rare) case where all points are on the identity line but
   70                 // not reaching the limits
   71             {
   72                 identity = false;
   73             }
   74 
   75             if(x[0] == 0.f && x[1] == 0.f)
   76                 // Avoid crash when first two points are at x = 0 (git Issue 2888)
   77             {
   78                 x[1] = 0.01f;
   79             }
   80 
   81             if(x[0] == 1.f && x[1] == 1.f)
   82                 // Avoid crash when first two points are at x = 1 (100 in gui) (git Issue 2923)
   83             {
   84                 x[0] = 0.99f;
   85             }
   86 
   87             if (!identity) {
   88                 if (kind == DCT_Spline && N > 2) {
   89                     spline_cubic_set ();
   90                 } else if (kind == DCT_NURBS && N > 2) {
   91                     NURBS_set ();
   92                     fillHash();
   93                 } else if (kind == DCT_CatumullRom && N > 2) {
   94                     catmull_rom_set();
   95                 } else {
   96                     kind = DCT_Linear;
   97                 }
   98             }
   99         } else if (kind == DCT_Parametric) {
  100             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)) {
  101                 identity = false;
  102 
  103                 x = new double[9];
  104                 x[0] = p[0];
  105 
  106                 for (int i = 1; i < 4; i++) {
  107                     x[i] = min(max(p[i], 0.001), 0.99);
  108                 }
  109 
  110                 for (int i = 4; i < 8; i++) {
  111                     x[i] = (p[i] + 100.0) / 200.0;
  112                 }
  113 
  114                 if (p.size() < 9) {
  115                     x[8] = 1.0;
  116                 } else {
  117                     x[8] = p[8] / 100.0;
  118                 }
  119 
  120                 mc = -xlog(2.0) / xlog(x[2]);
  121                 double mbase = pfull (0.5, x[8], x[6], x[5]);
  122                 mfc = mbase <= 1e-14 ? 0.0 : xexp(xlog(mbase) / mc);        // value of the curve at the center point
  123                 msc = -xlog(2.0) / xlog(x[1] / x[2]);
  124                 mhc = -xlog(2.0) / xlog((x[3] - x[2]) / (1 - x[2]));
  125             }
  126         }
  127 
  128         if (identity) {
  129             kind = DCT_Empty;
  130         }
  131     }
  132 }
  133 
  134 DiagonalCurve::~DiagonalCurve ()
  135 {
  136 
  137     delete [] x;
  138     delete [] y;
  139     delete [] ypp;
  140     poly_x.clear();
  141     poly_y.clear();
  142 }
  143 
  144 void DiagonalCurve::spline_cubic_set ()
  145 {
  146 
  147     double* u = new double[N - 1];
  148     delete [] ypp;      // TODO: why do we delete ypp here since it should not be allocated yet?
  149     ypp = new double [N];
  150 
  151     ypp[0] = u[0] = 0.0;    /* set lower boundary condition to "natural" */
  152 
  153     for (int i = 1; i < N - 1; ++i) {
  154         double sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
  155         double p = sig * ypp[i - 1] + 2.0;
  156         ypp[i] = (sig - 1.0) / p;
  157         u[i] = ((y[i + 1] - y[i])
  158                 / (x[i + 1] - x[i]) - (y[i] - y[i - 1]) / (x[i] - x[i - 1]));
  159         u[i] = (6.0 * u[i] / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p;
  160     }
  161 
  162     ypp[N - 1] = 0.0;
  163 
  164     for (int k = N - 2; k >= 0; --k) {
  165         ypp[k] = ypp[k] * ypp[k + 1] + u[k];
  166     }
  167 
  168     delete [] u;
  169 }
  170 
  171 void DiagonalCurve::NURBS_set ()
  172 {
  173 
  174     int nbSubCurvesPoints = N + (N - 3) * 2;
  175 
  176     std::vector<double> sc_x(nbSubCurvesPoints);  // X sub-curve points (  XP0,XP1,XP2,  XP2,XP3,XP4,  ...)
  177     std::vector<double> sc_y(nbSubCurvesPoints);  // Y sub-curve points (  YP0,YP1,YP2,  YP2,YP3,YP4,  ...)
  178     std::vector<double> sc_length(N + 2);         // Length of the subcurves
  179     double total_length = 0.;
  180 
  181     // Create the list of Bezier sub-curves
  182     // NURBS_set is called if N > 2 and non identity only
  183 
  184     int j = 0;
  185     int k = 0;
  186 
  187     for (int i = 0; i < N - 1;) {
  188         double length;
  189         double dx;
  190         double dy;
  191 
  192         // first point (on the curve)
  193         if (!i) {
  194             sc_x[j] = x[i];
  195             sc_y[j++] = y[i++];
  196         } else {
  197             sc_x[j] = (x[i - 1] + x[i]) / 2.;
  198             sc_y[j++] = (y[i - 1] + y[i]) / 2.;
  199         }
  200 
  201         // second point (control point)
  202         sc_x[j] = x[i];
  203         sc_y[j] = y[i++];
  204 
  205         dx = sc_x[j] - sc_x[j - 1];
  206         dy = sc_y[j] - sc_y[j - 1];
  207         length = sqrt(dx * dx + dy * dy);
  208         j++;
  209 
  210         // third point (on the curve)
  211         if (i == N - 1) {
  212             sc_x[j] = x[i];
  213             sc_y[j] = y[i];
  214         } else {
  215             sc_x[j] =  (x[i - 1] + x[i]) / 2.;
  216             sc_y[j] =  (y[i - 1] + y[i]) / 2.;
  217         }
  218 
  219         dx = sc_x[j] - sc_x[j - 1];
  220         dy = sc_y[j] - sc_y[j - 1];
  221         length += sqrt(dx * dx + dy * dy);
  222         j++;
  223 
  224         // Storing the length of all sub-curves and the total length (to have a better distribution
  225         // of the points along the curve)
  226         sc_length[k++] = length;
  227         total_length += length;
  228     }
  229 
  230     poly_x.clear();
  231     poly_y.clear();
  232     unsigned int sc_xsize = j - 1;
  233 
  234     // adding the initial horizontal segment, if any
  235     if (x[0] > 0.) {
  236         poly_x.push_back(0.);
  237         poly_y.push_back(y[0]);
  238     }
  239 
  240     // adding the initial horizontal segment, if any
  241     // create the polyline with the number of points adapted to the X range of the sub-curve
  242     for (unsigned int i = 0; i < sc_xsize /*sc_x.size()*/; i += 3) {
  243         // TODO: Speeding-up the interface by caching the polyline, instead of rebuilding it at each action on sliders !!!
  244         nbr_points = (int)(((double)(ppn + N - 2) * sc_length[i / 3] ) / total_length);
  245 
  246         if (nbr_points < 0) {
  247             for(unsigned int it = 0; it < sc_x.size(); it += 3) { // used unsigned int instead of size_t to avoid %zu in printf
  248                 printf("sc_length[%u/3]=%f \n", it, sc_length[it / 3]);
  249             }
  250 
  251             printf("NURBS diagonal curve: error detected!\n i=%u 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);
  252             exit(0);
  253         }
  254 
  255         // increment along the curve, not along the X axis
  256         increment = 1.0 / (double)(nbr_points - 1);
  257         x1 = sc_x[i];
  258         y1 = sc_y[i];
  259         x2 = sc_x[i + 1];
  260         y2 = sc_y[i + 1];
  261         x3 = sc_x[i + 2];
  262         y3 = sc_y[i + 2];
  263         firstPointIncluded = !i;
  264         AddPolygons ();
  265     }
  266 
  267     // adding the final horizontal segment, always (see under)
  268     poly_x.push_back(3.0);      // 3.0 is a hack for optimization purpose of the getVal method (the last value has to be beyond the normal range)
  269     poly_y.push_back(y[N - 1]);
  270 
  271     fillDyByDx();
  272 }
  273 
  274 
  275 /*****************************************************************************
  276  * Catmull Rom Spline
  277  * (https://en.wikipedia.org/wiki/Centripetal_Catmull%E2%80%93Rom_spline)
  278  *****************************************************************************/
  279 
  280 namespace {
  281 
  282 inline double pow2(double x)
  283 {
  284     return x*x;
  285 }
  286 
  287 
  288 inline double catmull_rom_tj(double ti,
  289                              double xi, double yi,
  290                              double xj, double yj)
  291 {
  292     // see https://github.com/Beep6581/RawTherapee/pull/4701#issuecomment-414054187
  293     static constexpr double alpha = 0.375;
  294     return pow(sqrt(pow2(xj-xi) + pow2(yj-yi)), alpha) + ti;
  295 }
  296 
  297 
  298 inline void catmull_rom_spline(int n_points,
  299                                double p0_x, double p0_y,
  300                                double p1_x, double p1_y,
  301                                double p2_x, double p2_y,
  302                                double p3_x, double p3_y,
  303                                std::vector<double> &res_x,
  304                                std::vector<double> &res_y)
  305 {
  306     res_x.reserve(n_points);
  307     res_y.reserve(n_points);
  308 
  309     double t0 = 0;
  310     double t1 = catmull_rom_tj(t0, p0_x, p0_y, p1_x, p1_y);
  311     double t2 = catmull_rom_tj(t1, p1_x, p1_y, p2_x, p2_y);
  312     double t3 = catmull_rom_tj(t2, p2_x, p2_y, p3_x, p3_y);
  313 
  314     double space = (t2-t1) / n_points;
  315 
  316     res_x.push_back(p1_x);
  317     res_y.push_back(p1_y);
  318 
  319     // special case, a segment at 0 or 1 is computed exactly
  320     if (p1_y == p2_y && (p1_y == 0 || p1_y == 1)) {
  321         for (int i = 1; i < n_points-1; ++i) {
  322             double t = p1_x + space * i;
  323             if (t >= p2_x) {
  324                 break;
  325             }
  326             res_x.push_back(t);
  327             res_y.push_back(p1_y);
  328         }
  329     } else {
  330         for (int i = 1; i < n_points-1; ++i) {
  331             double t = t1 + space * i;
  332         
  333             double c = (t1 - t)/(t1 - t0);
  334             double d = (t - t0)/(t1 - t0);
  335             double A1_x = c * p0_x + d * p1_x;
  336             double A1_y = c * p0_y + d * p1_y;
  337 
  338             c = (t2 - t)/(t2 - t1);
  339             d = (t - t1)/(t2 - t1);
  340             double A2_x = c * p1_x + d * p2_x;
  341             double A2_y = c * p1_y + d * p2_y;
  342 
  343             c = (t3 - t)/(t3 - t2);
  344             d = (t - t2)/(t3 - t2);
  345             double A3_x = c * p2_x + d * p3_x;
  346             double A3_y = c * p2_y + d * p3_y;
  347 
  348             c = (t2 - t)/(t2 - t0);
  349             d = (t - t0)/(t2 - t0);
  350             double B1_x = c * A1_x + d * A2_x;
  351             double B1_y = c * A1_y + d * A2_y;
  352 
  353             c = (t3 - t)/(t3 - t1);
  354             d = (t - t1)/(t3 - t1);
  355             double B2_x = c * A2_x + d * A3_x;
  356             double B2_y = c * A2_y + d * A3_y;        
  357 
  358             c = (t2 - t)/(t2 - t1);
  359             d = (t - t1)/(t2 - t1);
  360             double C_x = c * B1_x + d * B2_x;
  361             double C_y = c * B1_y + d * B2_y;
  362 
  363             res_x.push_back(C_x);
  364             res_y.push_back(C_y);
  365         }
  366     }
  367 
  368     res_x.push_back(p2_x);
  369     res_y.push_back(p2_y);
  370 }
  371 
  372 
  373 inline void catmull_rom_reflect(double px, double py, double cx, double cy,
  374                                 double &rx, double &ry)
  375 {
  376 #if 0
  377     double dx = px - cx;
  378     double dy = py - cy;
  379     rx = cx - dx;
  380     ry = cy - dy;
  381 #else
  382     // see https://github.com/Beep6581/RawTherapee/pull/4701#issuecomment-414054187
  383     static constexpr double epsilon = 1e-5;
  384     double dx = px - cx;
  385     double dy = py - cy;
  386     rx = cx - dx * 0.01;
  387     ry = dx > epsilon ? (dy / dx) * (rx - cx) + cy : cy;
  388 #endif    
  389 }
  390 
  391 
  392 void catmull_rom_chain(int n_points, int n_cp, double *x, double *y,
  393                        std::vector<double> &res_x, std::vector<double> &res_y)
  394 {
  395     double x_first, y_first;
  396     double x_last, y_last;
  397     catmull_rom_reflect(x[1], y[1], x[0], y[0], x_first, y_first);
  398     catmull_rom_reflect(x[n_cp-2], y[n_cp-2], x[n_cp-1], y[n_cp-1], x_last, y_last);
  399 
  400     int segments = n_cp - 1;
  401 
  402     res_x.reserve(n_points);
  403     res_y.reserve(n_points);
  404 
  405     for (int i = 0; i < segments; ++i) {
  406         int n = max(int(n_points * (x[i+1] - x[i]) + 0.5), 2);
  407         catmull_rom_spline(
  408             n, i == 0 ? x_first : x[i-1], i == 0 ? y_first : y[i-1],
  409             x[i], y[i], x[i+1], y[i+1],
  410             i == segments-1 ? x_last : x[i+2],
  411             i == segments-1 ? y_last : y[i+2],
  412             res_x, res_y);
  413     }
  414 }
  415 
  416 } // namespace
  417 
  418 
  419 void DiagonalCurve::catmull_rom_set()
  420 {
  421     int n_points = max(ppn * 65, 65000);
  422     poly_x.clear();
  423     poly_y.clear();
  424     catmull_rom_chain(n_points, N, x, y, poly_x, poly_y);
  425 }
  426 
  427 /*****************************************************************************/
  428 
  429 
  430 double DiagonalCurve::getVal (double t) const
  431 {
  432 
  433     switch (kind) {
  434 
  435     case DCT_Parametric : {
  436         if (t <= 1e-14) {
  437             return 0.0;
  438         }
  439 
  440         double tv = xexp(mc * xlog(t));
  441         double base = pfull (tv, x[8], x[6], x[5]);
  442         double stretched = base <= 1e-14 ? 0.0 : xexp(xlog(base) / mc);
  443 
  444         if (t < x[2]) {
  445             // add shadows effect:
  446             double stv = xexp(msc * xlog(stretched / mfc));
  447             double sbase = pfull (stv, x[8], x[7], 0.5);
  448             return mfc * (sbase <= 1e-14 ? 0.0 : xexp(xlog(sbase) / msc));
  449         } else {
  450             // add highlights effect:
  451             double htv = xexp(mhc * xlog((stretched - mfc) / (1 - mfc)));
  452             double hbase = pfull (htv, x[8], 0.5, x[4]);
  453             return mfc + (1 - mfc) * (hbase <= 1e-14 ? 0.0 : xexp(xlog(hbase) / mhc));
  454         }
  455 
  456         break;
  457     }
  458 
  459     case DCT_Linear :
  460     case DCT_Spline :
  461     {
  462         // values under and over the first and last point
  463         if (t > x[N - 1]) {
  464             return y[N - 1];
  465         } else if (t < x[0]) {
  466             return y[0];
  467         }
  468 
  469         // do a binary search for the right interval:
  470         unsigned int k_lo = 0, k_hi = N - 1;
  471 
  472         while (k_hi > 1 + k_lo) {
  473             unsigned int k = (k_hi + k_lo) / 2;
  474 
  475             if (x[k] > t) {
  476                 k_hi = k;
  477             } else {
  478                 k_lo = k;
  479             }
  480         }
  481 
  482         double h = x[k_hi] - x[k_lo];
  483 
  484         // linear
  485         if (kind == DCT_Linear) {
  486             return y[k_lo] + (t - x[k_lo]) * ( y[k_hi] - y[k_lo] ) / h;
  487         }
  488         // spline curve
  489         else { // if (kind==Spline) {
  490             double a = (x[k_hi] - t) / h;
  491             double b = (t - x[k_lo]) / h;
  492             double r = a * y[k_lo] + b * y[k_hi] + ((a * a * a - a) * ypp[k_lo] + (b * b * b - b) * ypp[k_hi]) * (h * h) * 0.1666666666666666666666666666666;
  493             return CLIPD(r);
  494         }
  495 
  496         break;
  497     }
  498 
  499     case DCT_CatumullRom: {
  500         auto it = std::lower_bound(poly_x.begin(), poly_x.end(), t);
  501         if (it == poly_x.end()) {
  502             return poly_y.back();
  503         }
  504         auto d = it - poly_x.begin();
  505         if (it+1 < poly_x.end() && t - *it > *(it+1) - t) {
  506             ++d;
  507         }
  508         return LIM01(*(poly_y.begin() + d));
  509     }
  510 
  511     case DCT_NURBS : {
  512         // get the hash table entry by rounding the value (previously multiplied by "hashSize")
  513         unsigned short int i = (unsigned short int)(t * hashSize);
  514 
  515         if (UNLIKELY(i > (hashSize + 1))) {
  516             //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.at(i), poly_x.size(), poly_x[hash.at(i)]);
  517             printf("\nOVERFLOW: hash #%d is used while seeking for value %.8f\n\n", i, t);
  518             return t;
  519         }
  520 
  521         unsigned int k_lo;
  522         unsigned int k_hi;
  523 
  524         k_lo = hash.at(i).smallerValue;
  525         k_hi = hash.at(i).higherValue;
  526 
  527         // do a binary search for the right interval :
  528         while (k_hi > 1 + k_lo) {
  529             unsigned int k = (k_hi + k_lo) / 2;
  530 
  531             if (poly_x[k] > t) {
  532                 k_hi = k;
  533             } else {
  534                 k_lo = k;
  535             }
  536         }
  537 
  538         return poly_y[k_lo] + (t - poly_x[k_lo]) * dyByDx[k_lo];
  539     }
  540 
  541     case DCT_Empty :
  542     default:
  543         // all other (unknown) kind
  544         return t;
  545     }
  546 }
  547 
  548 void DiagonalCurve::getVal (const std::vector<double>& t, std::vector<double>& res) const
  549 {
  550 
  551     res.resize (t.size());
  552 
  553     for (unsigned int i = 0; i < t.size(); i++) {
  554         res[i] = getVal(t[i]);
  555     }
  556 }
  557 
  558 }