CImgList<Tfloat> get_hessian(const char *const axes=0) const {
      CImgList<Tfloat> res;
      char __axes[12] = { 0 };
      const char *_axes = axes?axes:__axes;
      if (!axes) {
        unsigned int k = 0;
        if (_width>1) { __axes[k++] = 'x'; __axes[k++] = 'x'; }
        if (_width>1 && _height>1) { __axes[k++] = 'x'; __axes[k++] = 'y'; }
        if (_width>1 && _depth>1) { __axes[k++] = 'x'; __axes[k++] = 'z'; }
        if (_height>1) { __axes[k++] = 'y'; __axes[k++] = 'y'; }
        if (_height>1 && _depth>1) { __axes[k++] = 'y'; __axes[k++] = 'z'; }
        if (_depth>1) { __axes[k++] = 'z'; __axes[k++] = 'z'; }
      }
      const unsigned int len = (unsigned int)std::strlen(_axes);
      if (len%2)
        throw CImgArgumentException(_cimg_instance
                                    "get_hessian(): Invalid specified axes '%s'.",
                                    cimg_instance,
                                    axes);
      CImg<Tfloat> hess;
      for (unsigned int k = 0; k<len; k+=2) {
        const char
          _axis1 = cimg::lowercase(_axes[k]),
          _axis2 = cimg::lowercase(_axes[k + 1]),
          axis1 = std::min(_axis1,_axis2),
          axis2 = std::max(_axis2,_axis2);
        if (axis1!='x' && axis1!='y' && axis1!='z' &&
            axis2!='x' && axis2!='y' && axis2!='z')
          throw CImgArgumentException(_cimg_instance
                                      "get_hessian(): Invalid specified axes '%s'.",
                                      cimg_instance,
                                      axes);
        const longT off = axis1=='x'?1:axis1=='y'?_width:_width*_height;

        hess.assign(_width,_height,_depth,_spectrum);

        if (((axis1=='x' || axis2=='x') && _width==1) ||
            ((axis1=='y' || axis2=='y') && _height==1) ||
            ((axis1=='z' || axis2=='z') && _depth==1)) {
          hess.fill(0).move_to(res);
          continue;
        }

        if (axis1==axis2) // Ixx, Iyy, Izz
          cimg_pragma_openmp(parallel for cimg_openmp_collapse(3) cimg_openmp_if_size(size(),16384))
          cimg_forXYZC(*this,x,y,z,c) {
            const ulongT pos = offset(x,y,z,c);
            if ((axis1=='x' && !x) || (axis1=='y' && !y) || (axis1=='z' && !z))
              hess[pos] = (Tfloat)_data[pos + off] - _data[pos];
            else if ((axis1=='x' && x==width() - 1) ||
                     (axis1=='y' && y==height() - 1) ||
                     (axis1=='z' && z==depth() - 1))
              hess[pos] = (Tfloat)_data[pos - off] - _data[pos];
            else
              hess[pos] = (Tfloat)_data[pos + off] + _data[pos - off] - 2*_data[pos];
          }
        else if (axis1=='x' && axis2=='y') // Ixy
          cimg_pragma_openmp(parallel for cimg_openmp_collapse(2)
                             cimg_openmp_if(_width*_height>=(cimg_openmp_sizefactor)*16384 &&
                                            _depth*_spectrum>=2))
          cimg_forZC(*this,z,c) {
            CImg_3x3(I,Tfloat);
            cimg_for3x3(*this,x,y,z,c,I,Tfloat) hess(x,y,z,c) = (Inn + Ipp - Inp - Ipn)/4;
          }
        else if (axis1=='x' && axis2=='z') // Ixz
          cimg_pragma_openmp(parallel for cimg_openmp_if(_width*_height*_depth>=(cimg_openmp_sizefactor)*16384 &&
                                                         _spectrum>=2))
          cimg_forC(*this,c) {
            CImg_3x3x3(I,Tfloat);
            cimg_for3x3x3(*this,x,y,z,c,I,Tfloat) hess(x,y,z,c) = (Incn + Ipcp - Incp - Ipcn)/4;
          }
        else // Iyz
          cimg_pragma_openmp(parallel for cimg_openmp_if(_width*_height*_depth>=(cimg_openmp_sizefactor)*16384 &&
                                                         _spectrum>=2))
          cimg_forC(*this,c) {
            CImg_3x3x3(I,Tfloat);
            cimg_for3x3x3(*this,x,y,z,c,I,Tfloat) hess(x,y,z,c) = (Icnn + Icpp - Icnp - Icpn)/4;
          }
        hess.move_to(res);
      }
      return res;
    }