CImg<Tfloat> get_structure_tensors(const bool is_fwbw_scheme=false) const {
      if (is_empty()) return *this;
      CImg<Tfloat> res;
      if (_depth>1) { // 3D
        res.assign(_width,_height,_depth,6,0);
        if (!is_fwbw_scheme) { // Classical central finite differences
          cimg_pragma_openmp(parallel for cimg_openmp_if(_width*_height*_depth>=(cimg_openmp_sizefactor)*1048576 &&
                                                         _spectrum>=2))
          cimg_forC(*this,c) {
            Tfloat
              *ptrd0 = res.data(0,0,0,0), *ptrd1 = res.data(0,0,0,1), *ptrd2 = res.data(0,0,0,2),
              *ptrd3 = res.data(0,0,0,3), *ptrd4 = res.data(0,0,0,4), *ptrd5 = res.data(0,0,0,5);
            CImg_3x3x3(I,Tfloat);
            cimg_for3x3x3(*this,x,y,z,c,I,Tfloat) {
              const Tfloat
                ix = (Incc - Ipcc)/2,
                iy = (Icnc - Icpc)/2,
                iz = (Iccn - Iccp)/2;
              cimg_pragma_openmp(atomic) *(ptrd0++)+=ix*ix;
              cimg_pragma_openmp(atomic) *(ptrd1++)+=ix*iy;
              cimg_pragma_openmp(atomic) *(ptrd2++)+=ix*iz;
              cimg_pragma_openmp(atomic) *(ptrd3++)+=iy*iy;
              cimg_pragma_openmp(atomic) *(ptrd4++)+=iy*iz;
              cimg_pragma_openmp(atomic) *(ptrd5++)+=iz*iz;
            }
          }
        } else { // Forward/backward finite differences
          cimg_pragma_openmp(parallel for cimg_openmp_if(_width*_height*_depth>=(cimg_openmp_sizefactor)*1048576 &&
                                                         _spectrum>=2))
          cimg_forC(*this,c) {
            Tfloat
              *ptrd0 = res.data(0,0,0,0), *ptrd1 = res.data(0,0,0,1), *ptrd2 = res.data(0,0,0,2),
              *ptrd3 = res.data(0,0,0,3), *ptrd4 = res.data(0,0,0,4), *ptrd5 = res.data(0,0,0,5);
            CImg_3x3x3(I,Tfloat);
            cimg_for3x3x3(*this,x,y,z,c,I,Tfloat) {
              const Tfloat
                ixf = Incc - Iccc, ixb = Iccc - Ipcc, ixc = (Incc - Ipcc)/2,
                iyf = Icnc - Iccc, iyb = Iccc - Icpc, iyc = (Icnc - Icpc)/2,
                izf = Iccn - Iccc, izb = Iccc - Iccp, izc = (Iccn - Iccp)/2;
              cimg_pragma_openmp(atomic) *(ptrd0++)+=(ixf*ixf + ixb*ixb)/2;
              cimg_pragma_openmp(atomic) *(ptrd1++)+=ixc*iyc;
              cimg_pragma_openmp(atomic) *(ptrd2++)+=ixc*izc;
              cimg_pragma_openmp(atomic) *(ptrd3++)+=(iyf*iyf + iyb*iyb)/2;
              cimg_pragma_openmp(atomic) *(ptrd4++)+=iyc*izc;
              cimg_pragma_openmp(atomic) *(ptrd5++)+=(izf*izf + izb*izb)/2;
            }
          }
        }
      } else { // 2D
        res.assign(_width,_height,_depth,3,0);
        if (!is_fwbw_scheme) { // Classical central finite differences
          cimg_pragma_openmp(parallel for cimg_openmp_if(_width*_height>=(cimg_openmp_sizefactor)*1048576 &&
                                                         _depth*_spectrum>=2))
          cimg_forC(*this,c) {
            Tfloat *ptrd0 = res.data(0,0,0,0), *ptrd1 = res.data(0,0,0,1), *ptrd2 = res.data(0,0,0,2);
            CImg_3x3(I,Tfloat);
            cimg_for3x3(*this,x,y,0,c,I,Tfloat) {
              const Tfloat
                ix = (Inc - Ipc)/2,
                iy = (Icn - Icp)/2;
              cimg_pragma_openmp(atomic) *(ptrd0++)+=ix*ix;
              cimg_pragma_openmp(atomic) *(ptrd1++)+=ix*iy;
              cimg_pragma_openmp(atomic) *(ptrd2++)+=iy*iy;
            }
          }
        } else { // Forward/backward finite differences (version 2)
          cimg_pragma_openmp(parallel for cimg_openmp_if(_width*_height>=(cimg_openmp_sizefactor)*1048576 &&
                                                         _depth*_spectrum>=2))
          cimg_forC(*this,c) {
            Tfloat *ptrd0 = res.data(0,0,0,0), *ptrd1 = res.data(0,0,0,1), *ptrd2 = res.data(0,0,0,2);
            CImg_3x3(I,Tfloat);
            cimg_for3x3(*this,x,y,0,c,I,Tfloat) {
              const Tfloat
                ixf = Inc - Icc, ixb = Icc - Ipc, ixc = (Inc - Ipc)/2,
                iyf = Icn - Icc, iyb = Icc - Icp, iyc = (Icn - Icp)/2;
              cimg_pragma_openmp(atomic) *(ptrd0++)+=(ixf*ixf + ixb*ixb)/2;
              cimg_pragma_openmp(atomic) *(ptrd1++)+=ixc*iyc;
              cimg_pragma_openmp(atomic) *(ptrd2++)+=(iyf*iyf + iyb*iyb)/2;
            }
          }
        }
      }
      return res;
    }