CImg<Tfloat> get_blur_patch(const CImg<t>& guide,
                                const float sigma_s, const float sigma_r, const unsigned int patch_size=3,
                                const unsigned int lookup_size=4, const float smoothness=0,
                                const bool is_fast_approx=true) const {

#define _cimg_blur_patch3d_fast(N) { \
      cimg_pragma_openmp(parallel for cimg_openmp_collapse(2) \
                         cimg_openmp_if(res._width>=(cimg_openmp_sizefactor)*32 && res._height*res._depth>=4) \
                         firstprivate(P,Q)) \
      cimg_forXYZ(res,x,y,z) _cimg_abort_try_openmp2 { \
        cimg_abort_test2; \
        cimg_def##N##x##N##x##N(res,x,y,z); \
        tfloat *pP = P._data; cimg_forC(_guide,c) { cimg_get##N##x##N##x##N(_guide,x,y,z,c,pP,tfloat); pP+=N3; } \
        const int x0 = x - rsize1, y0 = y - rsize1, z0 = z - rsize1, \
                  x1 = x + rsize2, y1 = y + rsize2, z1 = z + rsize2; \
        tfloat sum_weights = 0; \
        cimg_for_in##N##XYZ(res,x0,y0,z0,x1,y1,z1,p,q,r) \
          if (cimg::abs(_guide(x,y,z,0) - _guide(p,q,r,0))<sigma_r3) { \
            tfloat *pQ = Q._data; cimg_forC(_guide,c) { cimg_get##N##x##N##x##N(_guide,p,q,r,c,pQ,tfloat); pQ+=N3; } \
            tfloat distance2 = 0; \
            pQ = Q._data; cimg_for(P,_pP,tfloat) { const tfloat dI = *_pP - *(pQ++); distance2+=dI*dI; } \
            distance2/=Pnorm; \
            const tfloat dx = (tfloat)p - x, dy = (tfloat)q - y, dz = (tfloat)r - z, \
              alldist = distance2 + (dx*dx + dy*dy + dz*dz)/sigma_s2, weight = alldist>3?0:1; \
            sum_weights+=weight; \
            cimg_forC(res,c) res(x,y,z,c)+=(Tfloat)weight*(*this)(p,q,r,c); \
          } \
        if (sum_weights>1e-10) cimg_forC(res,c) res(x,y,z,c)/=(Tfloat)sum_weights; \
        else cimg_forC(res,c) res(x,y,z,c) = (Tfloat)((*this)(x,y,z,c)); \
    } _cimg_abort_catch_openmp2 }

#define _cimg_blur_patch3d(N) { \
      cimg_pragma_openmp(parallel for cimg_openmp_collapse(2) \
                         cimg_openmp_if(res._width>=(cimg_openmp_sizefactor)*32 && res._height*res._depth>=4) \
                         firstprivate(P,Q)) \
      cimg_forXYZ(res,x,y,z) _cimg_abort_try_openmp2 { \
        cimg_abort_test2; \
        cimg_def##N##x##N##x##N(res,x,y,z); \
        tfloat *pP = P._data; cimg_forC(_guide,c) { cimg_get##N##x##N##x##N(_guide,x,y,z,c,pP,tfloat); pP+=N3; } \
        const int x0 = x - rsize1, y0 = y - rsize1, z0 = z - rsize1, \
                  x1 = x + rsize2, y1 = y + rsize2, z1 = z + rsize2; \
        tfloat sum_weights = 0, weight_max = 0; \
        cimg_for_in##N##XYZ(res,x0,y0,z0,x1,y1,z1,p,q,r) if (p!=x || q!=y || r!=z) { \
          tfloat *pQ = Q._data; cimg_forC(_guide,c) { cimg_get##N##x##N##x##N(_guide,p,q,r,c,pQ,tfloat); pQ+=N3; } \
          tfloat distance2 = 0; \
          pQ = Q._data; cimg_for(P,_pP,tfloat) { const tfloat dI = *_pP - *(pQ++); distance2+=dI*dI; } \
          distance2/=Pnorm; \
          const tfloat dx = (tfloat)p - x, dy = (tfloat)q - y, dz = (tfloat)r - z, \
            alldist = distance2 + (dx*dx + dy*dy + dz*dz)/sigma_s2, weight = std::exp(-alldist); \
          if (weight>weight_max) weight_max = weight; \
          sum_weights+=weight; \
          cimg_forC(res,c) res(x,y,z,c)+=(Tfloat)weight*(*this)(p,q,r,c); \
        } \
        sum_weights+=weight_max; cimg_forC(res,c) res(x,y,z,c)+=(Tfloat)weight_max*(*this)(x,y,z,c); \
        if (sum_weights>1e-10) cimg_forC(res,c) res(x,y,z,c)/=(Tfloat)sum_weights; \
        else cimg_forC(res,c) res(x,y,z,c) = (Tfloat)((*this)(x,y,z,c)); \
      } _cimg_abort_catch_openmp2 }

#define _cimg_blur_patch2d_fast(N) { \
        cimg_pragma_openmp(parallel for cimg_openmp_if(res._width>=(cimg_openmp_sizefactor)*32 && res._height>=4) \
                           firstprivate(P,Q)) \
        cimg_forXY(res,x,y) _cimg_abort_try_openmp2 { \
          cimg_abort_test2; \
          cimg_def##N##x##N(res,x,y); \
          tfloat *pP = P._data; cimg_forC(_guide,c) { cimg_get##N##x##N(_guide,x,y,0,c,pP,tfloat); pP+=N2; } \
          const int x0 = x - rsize1, y0 = y - rsize1, x1 = x + rsize2, y1 = y + rsize2; \
          tfloat sum_weights = 0; \
          cimg_for_in##N##XY(res,x0,y0,x1,y1,p,q) \
            if (cimg::abs(_guide(x,y,0,0) - _guide(p,q,0,0))<sigma_r3) { \
              tfloat *pQ = Q._data; cimg_forC(_guide,c) { cimg_get##N##x##N(_guide,p,q,0,c,pQ,tfloat); pQ+=N2; } \
              tfloat distance2 = 0; \
              pQ = Q._data; cimg_for(P,_pP,tfloat) { const tfloat dI = *_pP - *(pQ++); distance2+=dI*dI; } \
              distance2/=Pnorm; \
              const tfloat dx = (tfloat)p - x, dy = (tfloat)q - y, \
                alldist = distance2 + (dx*dx+dy*dy)/sigma_s2, weight = alldist>3?0:1; \
              sum_weights+=weight; \
              cimg_forC(res,c) res(x,y,c)+=(Tfloat)weight*(*this)(p,q,c); \
            } \
          if (sum_weights>1e-10) cimg_forC(res,c) res(x,y,c)/=(Tfloat)sum_weights; \
          else cimg_forC(res,c) res(x,y,c) = (Tfloat)((*this)(x,y,c)); \
        } _cimg_abort_catch_openmp2 }

#define _cimg_blur_patch2d(N) { \
        cimg_pragma_openmp(parallel for cimg_openmp_if(res._width>=(cimg_openmp_sizefactor)*32 && res._height>=4) \
                           firstprivate(P,Q)) \
        cimg_forXY(res,x,y) _cimg_abort_try_openmp2 { \
          cimg_abort_test2; \
          cimg_def##N##x##N(res,x,y); \
          tfloat *pP = P._data; cimg_forC(_guide,c) { cimg_get##N##x##N(_guide,x,y,0,c,pP,tfloat); pP+=N2; } \
          const int x0 = x - rsize1, y0 = y - rsize1, x1 = x + rsize2, y1 = y + rsize2; \
          tfloat sum_weights = 0, weight_max = 0; \
          cimg_for_in##N##XY(res,x0,y0,x1,y1,p,q) if (p!=x || q!=y) { \
            tfloat *pQ = Q._data; cimg_forC(_guide,c) { cimg_get##N##x##N(_guide,p,q,0,c,pQ,tfloat); pQ+=N2; } \
            tfloat distance2 = 0; \
            pQ = Q._data; cimg_for(P,_pP,tfloat) { const tfloat dI = *_pP - *(pQ++); distance2+=dI*dI; } \
            distance2/=Pnorm; \
            const tfloat dx = (tfloat)p - x, dy = (tfloat)q - y, \
              alldist = distance2 + (dx*dx+dy*dy)/sigma_s2, weight = std::exp(-alldist); \
            if (weight>weight_max) weight_max = weight; \
            sum_weights+=weight; \
            cimg_forC(res,c) res(x,y,c)+=(Tfloat)weight*(*this)(p,q,c); \
          } \
          sum_weights+=weight_max; cimg_forC(res,c) res(x,y,c)+=(Tfloat)weight_max*(*this)(x,y,c); \
          if (sum_weights>1e-10) cimg_forC(res,c) res(x,y,c)/=(Tfloat)sum_weights; \
          else cimg_forC(res,c) res(x,y,c) = (Tfloat)((*this)(x,y,c)); \
    } _cimg_abort_catch_openmp2 }

      typedef _cimg_tfloat tfloat;
      if (!is_sameXYZ(guide))
        throw CImgArgumentException(_cimg_instance
                                    "blur_patch(): Invalid size for specified guide image (%u,%u,%u,%u,%p).",
                                    cimg_instance,
                                    guide._width,guide._height,guide._depth,guide._spectrum,guide._data);
      if (is_empty() || !patch_size || !lookup_size) return +*this;
      Tfloat val_min, val_max = (Tfloat)max_min(val_min);
      _cimg_abort_init_openmp;
      cimg_abort_init;

      CImg<Tfloat> res(_width,_height,_depth,_spectrum,0);
      const CImg<tfloat>
        __guide = guide?CImg<tfloat>(guide,guide.pixel_type()==cimg::type<tfloat>::string()):
                        CImg<tfloat>(*this,pixel_type()==cimg::type<tfloat>::string()),
        _guide = smoothness>0?__guide.get_blur(smoothness):__guide.get_shared();
      CImg<tfloat> P(_guide._spectrum*patch_size*patch_size*(_depth>1?patch_size:1)), Q(P);

      t guide_min = (t)0, guide_max = (t)0;
      if (sigma_r<0) guide_max = guide.max_min(guide_min);
      const float
        guide_delta = (float)(guide_max - guide_min),
        _sigma_s = sigma_s>=0?sigma_s:-sigma_s*cimg::max(_width,_height,_depth)/100,
        _sigma_r = sigma_r>=0?sigma_r:-sigma_r*guide_delta/100,
        sigma_s2 = _sigma_s*_sigma_s,
        sigma_r2 = _sigma_r*_sigma_r,
        sigma_r3 = 3*_sigma_r,
        Pnorm = P.size()*sigma_r2;
      const int rsize2 = (int)lookup_size/2, rsize1 = (int)lookup_size - rsize2 - 1;
      const unsigned int N2 = patch_size*patch_size, N3 = N2*patch_size;
      cimg::unused(N2,N3);
      if (_depth>1) switch (patch_size) { // 3D
        case 2 : if (is_fast_approx) _cimg_blur_patch3d_fast(2) else _cimg_blur_patch3d(2) break;
        case 3 : if (is_fast_approx) _cimg_blur_patch3d_fast(3) else _cimg_blur_patch3d(3) break;
        default : {
          const int psize2 = (int)patch_size/2, psize1 = (int)patch_size - psize2 - 1;
          if (is_fast_approx) {
            cimg_pragma_openmp(parallel for cimg_openmp_collapse(2)
                               cimg_openmp_if(res._width>=(cimg_openmp_sizefactor)*32 && res._height*res._depth>=4)
                               firstprivate(P,Q))
            cimg_forXYZ(res,x,y,z) _cimg_abort_try_openmp2 { // Fast
              cimg_abort_test2;
              P = _guide.get_crop(x - psize1,y - psize1,z - psize1,x + psize2,y + psize2,z + psize2,true);
              const int x0 = x - rsize1, y0 = y - rsize1, z0 = z - rsize1,
                x1 = x + rsize2, y1 = y + rsize2, z1 = z + rsize2;
              tfloat sum_weights = 0;
              cimg_for_inXYZ(res,x0,y0,z0,x1,y1,z1,p,q,r)
                if (cimg::abs(_guide(x,y,z,0) - _guide(p,q,r,0))<sigma_r3) {
                  (Q = _guide.get_crop(p - psize1,q - psize1,r - psize1,p + psize2,q + psize2,r + psize2,true))-=P;
                  const tfloat
                    dx = (tfloat)x - p, dy = (tfloat)y - q, dz = (tfloat)z - r,
                    distance2 = (tfloat)(Q.pow(2).sum()/Pnorm + (dx*dx + dy*dy + dz*dz)/sigma_s2),
                    weight = distance2>3?0:1;
                  sum_weights+=weight;
                  cimg_forC(res,c) res(x,y,z,c)+=(Tfloat)weight*(*this)(p,q,r,c);
                }
              if (sum_weights>1e-10) cimg_forC(res,c) res(x,y,z,c)/=(Tfloat)sum_weights;
              else cimg_forC(res,c) res(x,y,z,c) = (Tfloat)((*this)(x,y,z,c));
            } _cimg_abort_catch_openmp2
          } else {
            cimg_pragma_openmp(parallel for cimg_openmp_collapse(2)
                               cimg_openmp_if(res._width>=(cimg_openmp_sizefactor)*32 && res._height*res._depth>=4)
                               firstprivate(P,Q))
            cimg_forXYZ(res,x,y,z) _cimg_abort_try_openmp2 { // Exact
              cimg_abort_test2;
              P = _guide.get_crop(x - psize1,y - psize1,z - psize1,x + psize2,y + psize2,z + psize2,true);
              const int x0 = x - rsize1, y0 = y - rsize1, z0 = z - rsize1,
                        x1 = x + rsize2, y1 = y + rsize2, z1 = z + rsize2;
              tfloat sum_weights = 0, weight_max = 0;
              cimg_for_inXYZ(res,x0,y0,z0,x1,y1,z1,p,q,r) if (p!=x || q!=y || r!=z) {
                (Q = _guide.get_crop(p - psize1,q - psize1,r - psize1,p + psize2,q + psize2,r + psize2,true))-=P;
                const tfloat
                  dx = (tfloat)x - p, dy = (tfloat)y - q, dz = (tfloat)z - r,
                  distance2 = (tfloat)(Q.pow(2).sum()/Pnorm + (dx*dx + dy*dy + dz*dz)/sigma_s2),
                  weight = std::exp(-distance2);
                if (weight>weight_max) weight_max = weight;
                sum_weights+=weight;
                cimg_forC(res,c) res(x,y,z,c)+=(Tfloat)weight*(*this)(p,q,r,c);
              }
              sum_weights+=weight_max; cimg_forC(res,c) res(x,y,z,c)+=(Tfloat)weight_max*(*this)(x,y,z,c);
              if (sum_weights>1e-10) cimg_forC(res,c) res(x,y,z,c)/=(Tfloat)sum_weights;
              else cimg_forC(res,c) res(x,y,z,c) = (Tfloat)((*this)(x,y,z,c));
            } _cimg_abort_catch_openmp2
          }
        }
        } else switch (patch_size) { // 2D
        case 2 : if (is_fast_approx) _cimg_blur_patch2d_fast(2) else _cimg_blur_patch2d(2) break;
        case 3 : if (is_fast_approx) _cimg_blur_patch2d_fast(3) else _cimg_blur_patch2d(3) break;
        case 4 : if (is_fast_approx) _cimg_blur_patch2d_fast(4) else _cimg_blur_patch2d(4) break;
        case 5 : if (is_fast_approx) _cimg_blur_patch2d_fast(5) else _cimg_blur_patch2d(5) break;
        case 6 : if (is_fast_approx) _cimg_blur_patch2d_fast(6) else _cimg_blur_patch2d(6) break;
        case 7 : if (is_fast_approx) _cimg_blur_patch2d_fast(7) else _cimg_blur_patch2d(7) break;
        case 8 : if (is_fast_approx) _cimg_blur_patch2d_fast(8) else _cimg_blur_patch2d(8) break;
        case 9 : if (is_fast_approx) _cimg_blur_patch2d_fast(9) else _cimg_blur_patch2d(9) break;
        default : { // Fast
          const int psize2 = (int)patch_size/2, psize1 = (int)patch_size - psize2 - 1;
          if (is_fast_approx) {
            cimg_pragma_openmp(parallel for cimg_openmp_if(res._width>=(cimg_openmp_sizefactor)*32 && res._height>=4)
                               firstprivate(P,Q))
            cimg_forXY(res,x,y) _cimg_abort_try_openmp2 { // Fast
              cimg_abort_test2;
              P = _guide.get_crop(x - psize1,y - psize1,x + psize2,y + psize2,true);
              const int x0 = x - rsize1, y0 = y - rsize1, x1 = x + rsize2, y1 = y + rsize2;
              tfloat sum_weights = 0;
              cimg_for_inXY(res,x0,y0,x1,y1,p,q)
                if (cimg::abs(_guide(x,y,0) - _guide(p,q,0))<sigma_r3) {
                  (Q = _guide.get_crop(p - psize1,q - psize1,p + psize2,q + psize2,true))-=P;
                  const tfloat
                    dx = (tfloat)x - p, dy = (tfloat)y - q,
                    distance2 = (tfloat)(Q.pow(2).sum()/Pnorm + (dx*dx + dy*dy)/sigma_s2),
                    weight = distance2>3?0:1;
                  sum_weights+=weight;
                  cimg_forC(res,c) res(x,y,c)+=(Tfloat)weight*(*this)(p,q,c);
                }
              if (sum_weights>1e-10) cimg_forC(res,c) res(x,y,c)/=(Tfloat)sum_weights;
              else cimg_forC(res,c) res(x,y,c) = (Tfloat)((*this)(x,y,c));
            } _cimg_abort_catch_openmp2
          } else {
            cimg_pragma_openmp(parallel for cimg_openmp_if(res._width>=(cimg_openmp_sizefactor)*32 && res._height>=4)
                               firstprivate(P,Q))
            cimg_forXY(res,x,y) _cimg_abort_try_openmp2 { // Exact
              cimg_abort_test2;
              P = _guide.get_crop(x - psize1,y - psize1,x + psize2,y + psize2,true);
              const int x0 = x - rsize1, y0 = y - rsize1, x1 = x + rsize2, y1 = y + rsize2;
              tfloat sum_weights = 0, weight_max = 0;
              cimg_for_inXY(res,x0,y0,x1,y1,p,q) if (p!=x || q!=y) {
                (Q = _guide.get_crop(p - psize1,q - psize1,p + psize2,q + psize2,true))-=P;
                const tfloat
                  dx = (tfloat)x - p, dy = (tfloat)y - q,
                  distance2 = (tfloat)(Q.pow(2).sum()/Pnorm + (dx*dx + dy*dy)/sigma_s2),
                  weight = std::exp(-distance2);
                if (weight>weight_max) weight_max = weight;
                sum_weights+=weight;
                cimg_forC(res,c) res(x,y,c)+=(Tfloat)weight*(*this)(p,q,c);
              }
              sum_weights+=weight_max; cimg_forC(res,c) res(x,y,c)+=(Tfloat)weight_max*(*this)(x,y,c);
              if (sum_weights>1e-10) cimg_forC(res,c) res(x,y,c)/=(Tfloat)sum_weights;
              else cimg_forC(res,c) res(x,y,c) = (Tfloat)((*this)(x,y,c));
            } _cimg_abort_catch_openmp2
          }
        }
        }
      cimg_abort_test;
      return res.cut(val_min,val_max);
    }