    CImg<Tfloat> get_pseudoinvert(const bool use_LU=false) const {

      // LU-based method.
      if (use_LU) {
        CImg<Tfloat> AtA(width(),width());
        cimg_pragma_openmp(parallel for cimg_openmp_if_size(_width*_height,128*128))
        cimg_forY(AtA,i)
          for (int j = 0; j<=i; ++j) {
            double res = 0;
            cimg_forY(*this,k) res+=(*this)(i,k)*(*this)(j,k);
            AtA(j,i) = AtA(i,j) = (Tfloat)res;
          }
        AtA.invert(true);
        return AtA*get_transpose();
      }

      // SVD-based method.
      CImg<Tfloat> U, S, V;
      SVD(U,S,V,false);
      const Tfloat epsilon = (sizeof(Tfloat)<=4?5.96e-8f:1.11e-16f)*std::max(_width,_height)*S.max();
      cimg_forX(V,x) {
        const Tfloat s = S(x), invs = s>epsilon?1/s:0;
        cimg_forY(V,y) V(x,y)*=invs;
      }
      return V*U.transpose();
    }
    CImg<Tfloat> get_pseudoinvert(const bool use_LU=false) const {

       
      if (use_LU) {
        CImg<Tfloat> AtA(width(),width());
        cimg_pragma_openmp(parallel for cimg_openmp_if_size(_width*_height,128*128))
        cimg_forY(AtA,i)
          for (int j = 0; j<=i; ++j) {
            double res = 0;
            cimg_forY(*this,k) res+=(*this)(i,k)*(*this)(j,k);
            AtA(j,i) = AtA(i,j) = (Tfloat)res;
          }
        AtA.invert(true);
        return AtA*get_transpose();
      }

       
      CImg<Tfloat> U, S, V;
      SVD(U,S,V,false);
      const Tfloat epsilon = (sizeof(Tfloat)<=4?5.96e-8f:1.11e-16f)*std::max(_width,_height)*S.max();
      cimg_forX(V,x) {
        const Tfloat s = S(x), invs = s>epsilon?1/s:0;
        cimg_forY(V,y) V(x,y)*=invs;
      }
      return V*U.transpose();
    }