From 3248e6b0035fd6fff88d8f228990ba51ca16bbc6 Mon Sep 17 00:00:00 2001
From: Jonathan Bieler <jonathan.bieler@alumni.epfl.ch>
Date: Sat, 12 Jun 2021 21:07:44 +0200
Subject: [PATCH] Reduce artifacts in parametric curve calculation (#6219)

---
 rtengine/curves.h          | 15 +++++++++++++++
 rtengine/diagonalcurves.cc | 28 ++++++++++++++++------------
 2 files changed, 31 insertions(+), 12 deletions(-)

diff --git a/rtengine/curves.h b/rtengine/curves.h
index 8800e54db..02503aff2 100644
--- a/rtengine/curves.h
+++ b/rtengine/curves.h
@@ -142,11 +142,19 @@ protected:
         double lx = xlog(x);
         return m2 * x + (1.0 - m2) * (2.0 - xexp(k * lx)) * xexp(l * lx);
     }
+    static inline double basel_alt(double x)
+    {
+        return (2.0 - x) * x * x * x;
+    }
     // basic concave function between (0,0) and (1,1). m1 and m2 controls the slope at the start and end point
     static inline double baseu(double x, double m1, double m2)
     {
         return 1.0 - basel(1.0 - x, m1, m2);
     }
+    static inline double baseu_alt(double x)
+    {
+        return x * (2.0 + (x - 2.0) * x * x);
+    }
     // convex curve between (0,0) and (1,1) with slope m at (0,0). hr controls the highlight recovery
     static inline double cupper(double x, double m, double hr)
     {
@@ -448,6 +456,13 @@ protected:
     {
         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);
     }
+    static inline double pfull_alt(double x, double sh, double hl)
+    {
+        double t = (1.0 - sh) * (1.0 - hl) * CurveFactory::basel_alt(x) + sh * hl * CurveFactory::baseu_alt(x);
+        return x <= 0.5
+            ? t + (1.0 - sh) * hl * CurveFactory::basel_alt(2.0 * x) * 0.5 + sh * (1.0 - hl) * CurveFactory::baseu_alt(2.0 * x) * 0.5
+            : t + (1.0 - sh) * hl * (0.5 + CurveFactory::baseu_alt(2.0 * x - 1.0) * 0.5) + sh * (1.0 - hl) * (0.5 + CurveFactory::basel_alt(2.0 * x - 1.0) * 0.5);
+    }
 
     void fillHash();
     void fillDyByDx();
diff --git a/rtengine/diagonalcurves.cc b/rtengine/diagonalcurves.cc
index 622683ea2..8a1c6ccbb 100644
--- a/rtengine/diagonalcurves.cc
+++ b/rtengine/diagonalcurves.cc
@@ -116,8 +116,8 @@ DiagonalCurve::DiagonalCurve (const std::vector<double>& p, int poly_pn)
                 }
 
                 mc = -xlog(2.0) / xlog(x[2]);
-                double mbase = pfull (0.5, x[8], x[6], x[5]);
-                mfc = mbase <= 1e-14 ? 0.0 : xexp(xlog(mbase) / mc);        // value of the curve at the center point
+                double mbase = pfull_alt (0.5, x[6], x[5]);
+                mfc = xexp(xlog(mbase) / mc);        // value of the curve at the center point
                 msc = -xlog(2.0) / xlog(x[1] / x[2]);
                 mhc = -xlog(2.0) / xlog((x[3] - x[2]) / (1 - x[2]));
             }
@@ -424,7 +424,6 @@ void DiagonalCurve::catmull_rom_set()
 
 /*****************************************************************************/
 
-
 double DiagonalCurve::getVal (double t) const
 {
 
@@ -435,20 +434,25 @@ double DiagonalCurve::getVal (double t) const
             return 0.0;
         }
 
-        double tv = xexp(mc * xlog(t));
-        double base = pfull (tv, x[8], x[6], x[5]);
-        double stretched = base <= 1e-14 ? 0.0 : xexp(xlog(base) / mc);
+        double tv = xexp(max(mc * xlog(t),-236.0)); // prevents numerical issues when calling pfull, at the cost of minor artifacts
+        double base = pfull_alt (tv, x[6], x[5]);
+        double stretched = xexp(xlog(base) / mc);
 
         if (t < x[2]) {
             // add shadows effect:
-            double stv = xexp(msc * xlog(stretched / mfc));
-            double sbase = pfull (stv, x[8], x[7], 0.5);
-            return mfc * (sbase <= 1e-14 ? 0.0 : xexp(xlog(sbase) / msc));
+            double stv = xexp(max(msc * xlog(stretched / mfc),-236.0));
+            double sbase = pfull_alt (stv, x[7], 0.5);
+            return mfc * xexp(xlog(sbase) / msc);
         } else {
             // add highlights effect:
-            double htv = xexp(mhc * xlog((stretched - mfc) / (1 - mfc)));
-            double hbase = pfull (htv, x[8], 0.5, x[4]);
-            return mfc + (1 - mfc) * (hbase <= 1e-14 ? 0.0 : xexp(xlog(hbase) / mhc));
+            double htv = xexp(max(mhc * xlog((stretched - mfc) / (1.0 - mfc)),-236.0));
+            double hbase = pfull_alt (htv, 0.5, x[4]);
+            //this part of the curve isn't affected by highlight, return the base curve
+            if (hbase < 1e-6 ){ 
+                return stretched;
+            } else {
+                return mfc + (1.0 - mfc) * xexp(xlog(hbase) / mhc);
+            }
         }
 
         break;