gistfile1.cpp

    uint w = m_image.width();
    uint h = m_image.height();
    uint m = m_image.size();
    const Scalar de = 1.0/gamma;
    const Scalar dc = -4.0/gamma + 1.0;
    printf("w h %i %i\n", w, h);
    using Matrix = Eigen::SparseMatrix<Scalar>;
    Matrix A(m, m);

    Eigen::VectorXd b(m);

    using T = Eigen::Triplet<double>;
    std::vector<T> coefficients;
    for (uint i = 0; i < h; i++) {
      for (uint j = 0; j < w; j++) {
        uint index = i*w + j;

        b(index) = m_image.data()[index];

        // boundary condition on land.
        //if (b(index) > 240.0) { b(index) -= b(index); }

        if (i > 0) coefficients.push_back(T(index-w, index, de));
        if (i < h-1) coefficients.push_back(T(index+w, index, de));
        if (j > 0) coefficients.push_back(T(index-1, index, de));
        if (j < w-1) coefficients.push_back(T(index+1, index, de));
        coefficients.push_back(T(index, index, dc));
      }
    }

    A.setFromTriplets(coefficients.begin(), coefficients.end());

    Eigen::SparseQR<Eigen::SparseMatrix<Scalar>, Eigen::COLAMDOrdering<int>> luA(A);
    int rank = luA.rank();
    std::cout << "Rank is " << rank << std::endl;

    Eigen::ConjugateGradient<Eigen::SparseMatrix<Scalar>> solver(A);
    Eigen::VectorXd x = solver.solveWithGuess(b, b);
    std::cout << "Solved." << std::endl;
    std::cout << "#iterations:     " << solver.iterations() << std::endl;
    std::cout << "estimated error: " << solver.error()      << std::endl;

    //x = A * x;

    Scalar min=std::numeric_limits<float>::max();
    Scalar max=std::numeric_limits<float>::lowest();
    for (uint i = 0; i < m; i++) {
      min = glm::min(min, x(i));
      max = glm::max(min, x(i));
    }
    for (uint i = 0; i < m; i++) {
      x(i) = (x(i) - min)/(max - min);
    }

    unsigned char* data = new unsigned char[3 * m];
    for (auto i = 0; i < w; i++) {
      for (auto j = 0; j < h; j++) {
        auto o = ((w-i-1)*w + j)*3;
        auto oo = (i*w + j);
        printf("X is %f\n", x(oo));
        data[o+0] = x(oo)*255.0;
        data[o+1] = x(oo)*255.0;
        data[o+2] = x(oo)*255.0;
      }
    }

    saveTGA("/tmp/processed", data, w, h);
    delete []data;