    static void _cimg_recursive_apply(T *data, const double filter[], const int N, const ulongT off,
                                      const unsigned int order, const bool boundary_conditions) {
      double val[4] = { 0 };  // res[n,n - 1,n - 2,n - 3,..] or res[n,n + 1,n + 2,n + 3,..]
      const double
        sumsq = filter[0], sum = sumsq * sumsq,
        a1 = filter[1], a2 = filter[2], a3 = filter[3],
        scaleM = 1. / ( (1. + a1 - a2 + a3) * (1. - a1 - a2 - a3) * (1. + a2 + (a1 - a3) * a3) );
      double M[9]; // Triggs matrix
      M[0] = scaleM * (-a3 * a1 + 1. - a3 * a3 - a2);
      M[1] = scaleM * (a3 + a1) * (a2 + a3 * a1);
      M[2] = scaleM * a3 * (a1 + a3 * a2);
      M[3] = scaleM * (a1 + a3 * a2);
      M[4] = -scaleM * (a2 - 1.) * (a2 + a3 * a1);
      M[5] = -scaleM * a3 * (a3 * a1 + a3 * a3 + a2 - 1.);
      M[6] = scaleM * (a3 * a1 + a2 + a1 * a1 - a2 * a2);
      M[7] = scaleM * (a1 * a2 + a3 * a2 * a2 - a1 * a3 * a3 - a3 * a3 * a3 - a3 * a2 + a3);
      M[8] = scaleM * a3 * (a1 + a3 * a2);
      switch (order) {
      case 0 : {
        const double iplus = (boundary_conditions?data[(N - 1)*off]:(T)0);
        for (int pass = 0; pass<2; ++pass) {
          if (!pass) {
            for (int k = 1; k<4; ++k) val[k] = (boundary_conditions?*data/sumsq:0);
          } else {
            // 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;
            val[0] = (M[0] * unp + M[1] * unp1 + M[2] * unp2 + vplus) * sum;
            val[1] = (M[3] * unp + M[4] * unp1 + M[5] * unp2 + vplus) * sum;
            val[2] = (M[6] * unp + M[7] * unp1 + M[8] * unp2 + vplus) * sum;
            *data = (T)val[0];
            data -= off;
            for (int k = 3; k>0; --k) val[k] = val[k - 1];
          }
          for (int n = pass; n<N; ++n) {
            val[0] = (*data);
            if (pass) val[0] *= sum;
            for (int k = 1; k<4; ++k) val[0] += val[k] * filter[k];
            *data = (T)val[0];
            if (!pass) data += off; else data -= off;
            for (int k = 3; k>0; --k) val[k] = val[k - 1];
          }
          if (!pass) data -= off;
        }
      } break;
      case 1 : {
        double x[3]; // [front,center,back]
        for (int pass = 0; pass<2; ++pass) {
          if (!pass) {
            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
            const double
              unp  = val[1], unp1 = val[2], unp2 = val[3];
            val[0] = (M[0] * unp + M[1] * unp1 + M[2] * unp2) * sum;
            val[1] = (M[3] * unp + M[4] * unp1 + M[5] * unp2) * sum;
            val[2] = (M[6] * unp + M[7] * unp1 + M[8] * unp2) * sum;
            *data = (T)val[0];
            data -= off;
            for (int k = 3; k>0; --k) val[k] = val[k - 1];
          }
          for (int n = pass; n<N - 1; ++n) {
            if (!pass) {
              x[0] = *(data + off);
              val[0] = 0.5f * (x[0] - x[2]);
            } else val[0] = (*data) * sum;
            for (int k = 1; k<4; ++k) val[0] += val[k] * filter[k];
            *data = (T)val[0];
            if (!pass) {
              data += off;
              for (int k = 2; k>0; --k) x[k] = x[k - 1];
            } else { data-=off;}
            for (int k = 3; k>0; --k) val[k] = val[k - 1];
          }
          *data = (T)0;
        }
      } break;
      case 2: {
        double x[3]; // [front,center,back]
        for (int pass = 0; pass<2; ++pass) {
          if (!pass) {
            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
            const double
              unp  = val[1], unp1 = val[2], unp2 = val[3];
            val[0] = (M[0] * unp + M[1] * unp1 + M[2] * unp2) * sum;
            val[1] = (M[3] * unp + M[4] * unp1 + M[5] * unp2) * sum;
            val[2] = (M[6] * unp + M[7] * unp1 + M[8] * unp2) * sum;
            *data = (T)val[0];
            data -= off;
            for (int k = 3; k>0; --k) val[k] = val[k - 1];
          }
          for (int n = pass; n<N - 1; ++n) {
            if (!pass) { x[0] = *(data + off); val[0] = (x[1] - x[2]); }
            else { x[0] = *(data - off); val[0] = (x[2] - x[1]) * sum; }
            for (int k = 1; k<4; ++k) val[0] += val[k]*filter[k];
            *data = (T)val[0];
            if (!pass) data += off; else data -= off;
            for (int k = 2; k>0; --k) x[k] = x[k - 1];
            for (int k = 3; k>0; --k) val[k] = val[k - 1];
          }
          *data = (T)0;
        }
      } break;
      case 3: {
        double x[3]; // [front,center,back]
        for (int pass = 0; pass<2; ++pass) {
          if (!pass) {
            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
            const double
              unp = val[1], unp1 = val[2], unp2 = val[3];
            val[0] = (M[0] * unp + M[1] * unp1 + M[2] * unp2) * sum;
            val[1] = (M[3] * unp + M[4] * unp1 + M[5] * unp2) * sum;
            val[2] = (M[6] * unp + M[7] * unp1 + M[8] * unp2) * sum;
            *data = (T)val[0];
            data -= off;
            for (int k = 3; k>0; --k) val[k] = val[k - 1];
          }
          for (int n = pass; n<N - 1; ++n) {
            if (!pass) { x[0] = *(data + off); val[0] = (x[0] - 2*x[1] + x[2]); }
            else { x[0] = *(data - off); val[0] = 0.5f * (x[2] - x[0]) * sum; }
            for (int k = 1; k<4; ++k) val[0] += val[k] * filter[k];
            *data = (T)val[0];
            if (!pass) data += off; else data -= off;
            for (int k = 2; k>0; --k) x[k] = x[k - 1];
            for (int k = 3; k>0; --k) val[k] = val[k - 1];
          }
          *data = (T)0;
        }
      } break;
      }
    }
    static void _cimg_recursive_apply(T *data, const double filter[], const int N, const ulongT off,
                                      const unsigned int order, const bool boundary_conditions) {
      double val[4] = { 0 };   
      const double
        sumsq = filter[0], sum = sumsq * sumsq,
        a1 = filter[1], a2 = filter[2], a3 = filter[3],
        scaleM = 1. / ( (1. + a1 - a2 + a3) * (1. - a1 - a2 - a3) * (1. + a2 + (a1 - a3) * a3) );
      double M[9];  
      M[0] = scaleM * (-a3 * a1 + 1. - a3 * a3 - a2);
      M[1] = scaleM * (a3 + a1) * (a2 + a3 * a1);
      M[2] = scaleM * a3 * (a1 + a3 * a2);
      M[3] = scaleM * (a1 + a3 * a2);
      M[4] = -scaleM * (a2 - 1.) * (a2 + a3 * a1);
      M[5] = -scaleM * a3 * (a3 * a1 + a3 * a3 + a2 - 1.);
      M[6] = scaleM * (a3 * a1 + a2 + a1 * a1 - a2 * a2);
      M[7] = scaleM * (a1 * a2 + a3 * a2 * a2 - a1 * a3 * a3 - a3 * a3 * a3 - a3 * a2 + a3);
      M[8] = scaleM * a3 * (a1 + a3 * a2);
      switch (order) {
      case 0 : {
        const double iplus = (boundary_conditions?data[(N - 1)*off]:(T)0);
        for (int pass = 0; pass<2; ++pass) {
          if (!pass) {
            for (int k = 1; k<4; ++k) val[k] = (boundary_conditions?*data/sumsq:0);
          } else {
             
            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;
            val[0] = (M[0] * unp + M[1] * unp1 + M[2] * unp2 + vplus) * sum;
            val[1] = (M[3] * unp + M[4] * unp1 + M[5] * unp2 + vplus) * sum;
            val[2] = (M[6] * unp + M[7] * unp1 + M[8] * unp2 + vplus) * sum;
            *data = (T)val[0];
            data -= off;
            for (int k = 3; k>0; --k) val[k] = val[k - 1];
          }
          for (int n = pass; n<N; ++n) {
            val[0] = (*data);
            if (pass) val[0] *= sum;
            for (int k = 1; k<4; ++k) val[0] += val[k] * filter[k];
            *data = (T)val[0];
            if (!pass) data += off; else data -= off;
            for (int k = 3; k>0; --k) val[k] = val[k - 1];
          }
          if (!pass) data -= off;
        }
      } break;
      case 1 : {
        double x[3];  
        for (int pass = 0; pass<2; ++pass) {
          if (!pass) {
            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 {
             
            const double
              unp  = val[1], unp1 = val[2], unp2 = val[3];
            val[0] = (M[0] * unp + M[1] * unp1 + M[2] * unp2) * sum;
            val[1] = (M[3] * unp + M[4] * unp1 + M[5] * unp2) * sum;
            val[2] = (M[6] * unp + M[7] * unp1 + M[8] * unp2) * sum;
            *data = (T)val[0];
            data -= off;
            for (int k = 3; k>0; --k) val[k] = val[k - 1];
          }
          for (int n = pass; n<N - 1; ++n) {
            if (!pass) {
              x[0] = *(data + off);
              val[0] = 0.5f * (x[0] - x[2]);
            } else val[0] = (*data) * sum;
            for (int k = 1; k<4; ++k) val[0] += val[k] * filter[k];
            *data = (T)val[0];
            if (!pass) {
              data += off;
              for (int k = 2; k>0; --k) x[k] = x[k - 1];
            } else { data-=off;}
            for (int k = 3; k>0; --k) val[k] = val[k - 1];
          }
          *data = (T)0;
        }
      } break;
      case 2: {
        double x[3];  
        for (int pass = 0; pass<2; ++pass) {
          if (!pass) {
            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 {
             
            const double
              unp  = val[1], unp1 = val[2], unp2 = val[3];
            val[0] = (M[0] * unp + M[1] * unp1 + M[2] * unp2) * sum;
            val[1] = (M[3] * unp + M[4] * unp1 + M[5] * unp2) * sum;
            val[2] = (M[6] * unp + M[7] * unp1 + M[8] * unp2) * sum;
            *data = (T)val[0];
            data -= off;
            for (int k = 3; k>0; --k) val[k] = val[k - 1];
          }
          for (int n = pass; n<N - 1; ++n) {
            if (!pass) { x[0] = *(data + off); val[0] = (x[1] - x[2]); }
            else { x[0] = *(data - off); val[0] = (x[2] - x[1]) * sum; }
            for (int k = 1; k<4; ++k) val[0] += val[k]*filter[k];
            *data = (T)val[0];
            if (!pass) data += off; else data -= off;
            for (int k = 2; k>0; --k) x[k] = x[k - 1];
            for (int k = 3; k>0; --k) val[k] = val[k - 1];
          }
          *data = (T)0;
        }
      } break;
      case 3: {
        double x[3];  
        for (int pass = 0; pass<2; ++pass) {
          if (!pass) {
            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 {
             
            const double
              unp = val[1], unp1 = val[2], unp2 = val[3];
            val[0] = (M[0] * unp + M[1] * unp1 + M[2] * unp2) * sum;
            val[1] = (M[3] * unp + M[4] * unp1 + M[5] * unp2) * sum;
            val[2] = (M[6] * unp + M[7] * unp1 + M[8] * unp2) * sum;
            *data = (T)val[0];
            data -= off;
            for (int k = 3; k>0; --k) val[k] = val[k - 1];
          }
          for (int n = pass; n<N - 1; ++n) {
            if (!pass) { x[0] = *(data + off); val[0] = (x[0] - 2*x[1] + x[2]); }
            else { x[0] = *(data - off); val[0] = 0.5f * (x[2] - x[0]) * sum; }
            for (int k = 1; k<4; ++k) val[0] += val[k] * filter[k];
            *data = (T)val[0];
            if (!pass) data += off; else data -= off;
            for (int k = 2; k>0; --k) x[k] = x[k - 1];
            for (int k = 3; k>0; --k) val[k] = val[k - 1];
          }
          *data = (T)0;
        }
      } break;
      }
    }