CImg<Tfloat> get_project_matrix(const CImg<t>& dictionary, const unsigned int method=0,
                                    const unsigned int max_iter=0, const double max_residual=1e-6) const {
      if (_depth!=1 || _spectrum!=1)
        throw CImgInstanceException(_cimg_instance
                                    "project_matrix(): Instance image is not a matrix.",
                                    cimg_instance);
      if (dictionary._height!=_height || dictionary._depth!=1 || dictionary._spectrum!=1)
        throw CImgArgumentException(_cimg_instance
                                    "project_matrix(): Specified dictionary (%u,%u,%u,%u) has an invalid size.",
                                    cimg_instance,
                                    dictionary._width,dictionary._height,dictionary._depth,dictionary._spectrum);

      if (!method) return get_solve(dictionary,true);
      CImg<Tfloat> W(_width,dictionary._width,1,1,0);

      // Compute dictionary norm and normalize it.
      CImg<Tfloat> D(dictionary,false), Dnorm(D._width);
      cimg_pragma_openmp(parallel for cimg_openmp_if(_width>=2 && _width*_height>=32))
      cimg_forX(Dnorm,d) {
        Tfloat norm = 0;
        cimg_forY(D,y) norm+=cimg::sqr(D(d,y));
        Dnorm[d] = std::max((Tfloat)1e-8,std::sqrt(norm));
      }
      cimg_forXY(D,d,y) D(d,y)/=Dnorm[d];

      // Matching pursuit.
      const unsigned int proj_step = method<3?1:method - 2;
      bool is_orthoproj = false;

      cimg_pragma_openmp(parallel for cimg_openmp_if(_width>=2 && _width*_height>=32))
        cimg_forX(*this,x) {
        CImg<Tfloat> S = get_column(x);
        const CImg<Tfloat> S0 = method<2?CImg<Tfloat>():S;
        Tfloat residual = S.magnitude()/S._height;
        const unsigned int nmax = max_iter?max_iter:D._width;

        for (unsigned int n = 0; n<nmax && residual>max_residual; ++n) {

          // Find best matching column in D.
          int dmax = 0;
          Tfloat absdotmax = 0, dotmax = 0;
          cimg_pragma_openmp(parallel for cimg_openmp_if(D._width>=2 && D._width*D._height>=32))
          cimg_forX(D,d) {
            Tfloat _dot = 0;
            cimg_forY(D,y) _dot+=S[y]*D(d,y);
            Tfloat absdot = cimg::abs(_dot);
            cimg_pragma_openmp(critical(get_project_matrix)) {
              if (absdot>absdotmax) {
                absdotmax = absdot;
                dotmax = _dot;
                dmax = d;
              }
            }
          }

          if (!n || method<3 || n%proj_step) {
            // Matching Pursuit: Subtract component to signal.
            W(x,dmax)+=dotmax;
            residual = 0;
            cimg_forY(S,y) {
              S[y]-=dotmax*D(dmax,y);
              residual+=cimg::sqr(S[y]);
            }
            residual = std::sqrt(residual)/S._height;
            is_orthoproj = false;

          } else {
            // Orthogonal Matching Pursuit: Orthogonal projection step.
            W(x,dmax) = 1; // Used as a marker only.
            unsigned int nbW = 0;
            cimg_forY(W,d) if (W(x,d)) ++nbW;
            CImg<Tfloat> sD(nbW,D._height);
            CImg<uintT> inds(nbW);
            int sd = 0;
            cimg_forY(W,d) if (W(x,d)) {
              cimg_forY(sD,y) sD(sd,y) = D(d,y);
              inds[sd++] = d;
            }
            S0.get_solve(sD,true).move_to(sD); // sD is now a one-column vector of weights

            // Recompute residual signal.
            S = S0;
            cimg_forY(sD,k) {
              const Tfloat weight = sD[k];
              const unsigned int ind = inds[k];
              W(x,ind) = weight;
              cimg_forY(S,y) S[y]-=weight*D(ind,y);
            }
            residual = S.magnitude()/S._height;
            is_orthoproj = true;
          }
        }

        // Perform last orthoprojection step if needed.
        if (method>=2 && !is_orthoproj) {
          unsigned int nbW = 0;
          cimg_forY(W,d) if (W(x,d)) ++nbW;
          if (nbW) { // Avoid degenerated case where 0 coefs are used
            CImg<Tfloat> sD(nbW,D._height);
            CImg<uintT> inds(nbW);
            int sd = 0;
            cimg_forY(W,d) if (W(x,d)) {
              cimg_forY(sD,y) sD(sd,y) = D(d,y);
              inds[sd++] = d;
            }
            S0.get_solve(sD,true).move_to(sD);
            cimg_forY(sD,k) W(x,inds[k]) = sD[k];
          }
        }
      }

      // Normalize resulting coefficients according to initial (non-normalized) dictionary.
      cimg_forXY(W,x,y) W(x,y)/=Dnorm[y];
      return W;
    }