From 7594dc59c079904022161d5030a3c2985263e4b6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20Tschumperl=C3=A9?= Date: Tue, 28 Jul 2020 08:50:31 +0200 Subject: [PATCH] CImg::deriche() : do calculus in 'double' precision by default. --- CImg.h | 40 ++++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/CImg.h b/CImg.h index a589218f..0ef265a1 100644 --- a/CImg.h +++ b/CImg.h @@ -38392,52 +38392,52 @@ namespace cimg_library_suffixed { CImg& deriche(const float sigma, const unsigned int order=0, const char axis='x', const bool boundary_conditions=true) { #define _cimg_deriche_apply \ - CImg Y(N); \ - Tfloat *ptrY = Y._data, yb = 0, yp = 0; \ + CImg 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=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;