CImg<T>::deriche() : do calculus in 'double' precision by default.

CImg.h

@@ -38392,52 +38392,52 @@ namespace cimg_library_suffixed {
     CImg<T>& deriche(const float sigma, const unsigned int order=0, const char axis='x',
                      const bool boundary_conditions=true) {
 #define _cimg_deriche_apply \
-  CImg<Tfloat> Y(N); \
-  Tfloat *ptrY = Y._data, yb = 0, yp = 0; \
+  CImg<doubleT> Y(N); \
+  double *ptrY = Y._data, yb = 0, yp = 0; \
   T xp = (T)0; \
-  if (boundary_conditions) { xp = *ptrX; yb = yp = (Tfloat)(coefp*xp); } \
+  if (boundary_conditions) { xp = *ptrX; yb = yp = (double)(coefp*xp); } \
   for (int m = 0; m<N; ++m) { \
     const T xc = *ptrX; ptrX+=off; \
-    const Tfloat yc = *(ptrY++) = (Tfloat)(a0*xc + a1*xp - b1*yp - b2*yb); \
+    const double yc = *(ptrY++) = (double)(a0*xc + a1*xp - b1*yp - b2*yb); \
     xp = xc; yb = yp; yp = yc; \
   } \
   T xn = (T)0, xa = (T)0; \
-  Tfloat yn = 0, ya = 0; \
-  if (boundary_conditions) { xn = xa = *(ptrX-off); yn = ya = (Tfloat)coefn*xn; } \
+  double yn = 0, ya = 0; \
+  if (boundary_conditions) { xn = xa = *(ptrX - off); yn = ya = (double)coefn*xn; } \
   for (int n = N - 1; n>=0; --n) { \
     const T xc = *(ptrX-=off); \
-    const Tfloat yc = (Tfloat)(a2*xn + a3*xa - b1*yn - b2*ya); \
+    const double yc = (double)(a2*xn + a3*xa - b1*yn - b2*ya); \
     xa = xn; xn = xc; ya = yn; yn = yc; \
     *ptrX = (T)(*(--ptrY)+yc); \
   }
       const char naxis = cimg::lowercase(axis);
-      const float nsigma = sigma>=0?sigma:-sigma*(naxis=='x'?_width:naxis=='y'?_height:naxis=='z'?_depth:_spectrum)/100;
+      const double nsigma = sigma>=0?sigma:-sigma*(naxis=='x'?_width:naxis=='y'?_height:naxis=='z'?_depth:_spectrum)/100;
       if (is_empty() || (nsigma<0.1f && !order)) return *this;
-      const float
+      const double
         nnsigma = nsigma<0.1f?0.1f:nsigma,
         alpha = 1.695f/nnsigma,
-        ema = (float)std::exp(-alpha),
-        ema2 = (float)std::exp(-2*alpha),
+        ema = std::exp(-alpha),
+        ema2 = std::exp(-2*alpha),
         b1 = -2*ema,
         b2 = ema2;
-      float a0 = 0, a1 = 0, a2 = 0, a3 = 0, coefp = 0, coefn = 0;
+      double a0 = 0, a1 = 0, a2 = 0, a3 = 0, coefp = 0, coefn = 0;
       switch (order) {
       case 0 : {
-        const float k = (1-ema)*(1-ema)/(1 + 2*alpha*ema-ema2);
+        const double k = (1-ema)*(1-ema)/(1 + 2*alpha*ema-ema2);
         a0 = k;
         a1 = k*(alpha - 1)*ema;
         a2 = k*(alpha + 1)*ema;
         a3 = -k*ema2;
       } break;
       case 1 : {
-        const float k = -(1-ema)*(1-ema)*(1-ema)/(2*(ema + 1)*ema);
+        const double k = -(1-ema)*(1-ema)*(1-ema)/(2*(ema + 1)*ema);
 	a0 = a3 = 0;
 	a1 = k*ema;
         a2 = -a1;
       } break;
       case 2 : {
-        const float
-          ea = (float)std::exp(-alpha),
+        const double
+          ea = std::exp(-alpha),
           k = -(ema2 - 1)/(2*alpha*ema),
           kn = (-2*(-1 + 3*ea - 3*ea*ea + ea*ea*ea)/(3*ea + 1 + 3*ea*ea + ea*ea*ea));
         a0 = kn;
@@ -38527,7 +38527,7 @@ namespace cimg_library_suffixed {
 	  if (!pass) {
 	    for (int k = 1; k<4; ++k) val[k] = (boundary_conditions?*data/sumsq:0);
 	  } else {
-	    // apply Triggs boundary conditions
+	    // Apply Triggs boundary conditions
 	    const double
 	      uplus = iplus/(1. - a1 - a2 - a3), vplus = uplus/(1. - a1 - a2 - a3),
 	      unp  = val[1] - uplus, unp1 = val[2] - uplus, unp2 = val[3] - uplus;
@@ -38556,7 +38556,7 @@ namespace cimg_library_suffixed {
 	    for (int k = 0; k<3; ++k) x[k] = (boundary_conditions?*data:(T)0);
 	    for (int k = 0; k<4; ++k) val[k] = 0;
 	  } else {
-	    // apply Triggs boundary conditions
+	    // Apply Triggs boundary conditions
 	    const double
 	      unp  = val[1], unp1 = val[2], unp2 = val[3];
 	    val[0] = (M[0] * unp + M[1] * unp1 + M[2] * unp2) * sum;
@@ -38589,7 +38589,7 @@ namespace cimg_library_suffixed {
 	    for (int k = 0; k<3; ++k) x[k] = (boundary_conditions?*data:(T)0);
 	    for (int k = 0; k<4; ++k) val[k] = 0;
 	  } else {
-	    // apply Triggs boundary conditions
+	    // Apply Triggs boundary conditions
 	    const double
 	      unp  = val[1], unp1 = val[2], unp2 = val[3];
 	    val[0] = (M[0] * unp + M[1] * unp1 + M[2] * unp2) * sum;
@@ -38618,7 +38618,7 @@ namespace cimg_library_suffixed {
 	    for (int k = 0; k<3; ++k) x[k] = (boundary_conditions?*data:(T)0);
 	    for (int k = 0; k<4; ++k) val[k] = 0;
 	  } else {
-	    // apply Triggs boundary conditions
+	    // Apply Triggs boundary conditions
 	    const double
 	      unp = val[1], unp1 = val[2], unp2 = val[3];
 	    val[0] = (M[0] * unp + M[1] * unp1 + M[2] * unp2) * sum;