CImg<T>& deriche(const float sigma, const unsigned int order=0, const char axis='x',
                     const unsigned int boundary_conditions=1) {
#define _cimg_deriche_apply \
  CImg<doubleT> Y(N); \
  double *ptrY = Y._data, yb = 0, yp = 0; \
  T xp = (T)0; \
  if (boundary_conditions) { xp = *ptrX; yb = yp = (double)(coefp*xp); } \
  for (int m = 0; m<N; ++m) { \
    const T xc = *ptrX; ptrX+=off; \
    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; \
  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 double yc = (double)(a2*xn + a3*xa - b1*yn - b2*ya); \
    xa = xn; xn = xc; ya = yn; yn = yc; \
    *ptrX = (T)(*(--ptrY)+yc); \
  }

      if (order>2)
        throw CImgArgumentException(_cimg_instance
                                    "deriche(): Invalid specified order '%d' "
                                    "('order' can be { 0=smoothing | 1=1st-derivative | 2=2nd-derivative }).",
                                    cimg_instance,
                                    order);

      const char naxis = cimg::lowercase(axis);
      if (naxis!='x' && naxis!='y' && naxis!='z' && naxis!='c')
        throw CImgArgumentException(_cimg_instance
                                    "deriche(): Invalid specified axis '%c'.",
                                    cimg_instance,
                                    axis);
      const double
        nsigma = sigma>=0?sigma:-sigma*(naxis=='x'?_width:
                                        naxis=='y'?_height:
                                        naxis=='z'?_depth:_spectrum)/100,
        nnsigma = nsigma<0.1f?0.1f:nsigma;

      if (is_empty() || (nsigma<0.1f && !order)) return *this;
      if (boundary_conditions>1) {
        const int w = width(), h = height(), d = depth(), s = spectrum(), border = (int)cimg::round(1 + 3*nnsigma);
        switch (naxis) {
        case 'x' :
          return draw_image(get_resize(w + 2*border,h,d,s,0,boundary_conditions,0.5).
                            deriche(nnsigma,order,naxis,1).columns(border,w - 1 + border));
        case 'y' :
          return draw_image(get_resize(w,h + 2*border,d,s,0,boundary_conditions,0,0.5).
                            deriche(nnsigma,order,naxis,1).rows(border,h - 1 + border));
        case 'z' :
          return draw_image(get_resize(w,h,d + 2*border,s,0,boundary_conditions,0,0,0.5).
                            deriche(nnsigma,order,naxis,1).slices(border,d - 1 + border));
        default :
          return draw_image(get_resize(w,h,d,s + 2*border,0,boundary_conditions,0,0,0,0.5).
                            deriche(nnsigma,order,naxis,1).channels(border,d - 1 + border));
        }
      }

      const double
        alpha = 1.695f/nnsigma,
        ema = std::exp(-alpha),
        ema2 = std::exp(-2*alpha),
        b1 = -2*ema,
        b2 = ema2;
      double a0 = 0, a1 = 0, a2 = 0, a3 = 0, coefp = 0, coefn = 0;
      switch (order) {
      case 0 : {
        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 double k = -(1-ema)*(1-ema)*(1-ema)/(2*(ema + 1)*ema);
        a0 = a3 = 0;
        a1 = k*ema;
        a2 = -a1;
      } break;
      default : {
        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;
        a1 = -kn*(1 + k*alpha)*ema;
        a2 = kn*(1 - k*alpha)*ema;
        a3 = -kn*ema2;
      } break;
      }

      coefp = (a0 + a1)/(1 + b1 + b2);
      coefn = (a2 + a3)/(1 + b1 + b2);
      switch (naxis) {
      case 'x' : {
        const int N = width();
        const ulongT off = 1U;
        cimg_pragma_openmp(parallel for cimg_openmp_collapse(3) cimg_openmp_if(_width>=(cimg_openmp_sizefactor)*256 &&
                                                                   _height*_depth*_spectrum>=16))
        cimg_forYZC(*this,y,z,c) { T *ptrX = data(0,y,z,c); _cimg_deriche_apply; }
      } break;
      case 'y' : {
        const int N = height();
        const ulongT off = (ulongT)_width;
        cimg_pragma_openmp(parallel for cimg_openmp_collapse(3) cimg_openmp_if(_width>=(cimg_openmp_sizefactor)*256 &&
                                                                   _height*_depth*_spectrum>=16))
        cimg_forXZC(*this,x,z,c) { T *ptrX = data(x,0,z,c); _cimg_deriche_apply; }
      } break;
      case 'z' : {
        const int N = depth();
        const ulongT off = (ulongT)_width*_height;
        cimg_pragma_openmp(parallel for cimg_openmp_collapse(3) cimg_openmp_if(_width>=(cimg_openmp_sizefactor)*256 &&
                                                                   _height*_depth*_spectrum>=16))
        cimg_forXYC(*this,x,y,c) { T *ptrX = data(x,y,0,c); _cimg_deriche_apply; }
      } break;
      default : {
        const int N = spectrum();
        const ulongT off = (ulongT)_width*_height*_depth;
        cimg_pragma_openmp(parallel for cimg_openmp_collapse(3) cimg_openmp_if(_width>=(cimg_openmp_sizefactor)*256 &&
                                                                   _height*_depth*_spectrum>=16))
        cimg_forXYZ(*this,x,y,z) { T *ptrX = data(x,y,z,0); _cimg_deriche_apply; }
      }
      }
      return *this;
    }