Solving geodesics on triangular meshes gives negative distances

Question

I have implemented the heat method for geodesics:

https://www.cs.cmu.edu/~kmcrane/Projects/HeatMethod/paperCACM.pdf

When I run it I am getting a solution that, visually, seems correct:

In this image, the arrows are the negated gradient, and the color is the distance, meaning that as a point's color goes towards white that point is farther away from the source vertex.

The image I am seeing is, visually, what I would expect, however I am noticing negative values. If I print the absolute value of the distance at each point, rather than the distance, I get:

Essentially, points close to the source vertex are getting negative values.

I have checked the code that is calculating this multiple times and as far as I can tell I am doing things correctly. The paper requires solving 2 linear systems, both visible in the first column of page 5 in the paper.

I am not entirely sure why some of my distances are negative.

vector<Eigen::Vector3f> gradient;
Eigen::VectorXf HeatGeodesics(
    const HMesh<VertexData>& mesh, const Eigen::VectorXf& sources)
{
    Eigen::DiagonalMatrix<float, Eigen::Dynamic> areas(sources.size());
    Eigen::SparseMatrix<float> laplace_operator(sources.size(), sources.size());

    for(uint i = 0; i < sources.size(); i++)
        laplace_operator.insert(i, i) = 0;

    areas.setZero();

    for(uint e = 0; e < mesh.Edges().size(); e += 2)
    {
        const auto& edge = mesh.Edges()[e];
        const uint i = edge.Vert().ID();
        const uint j = edge.Pair().Vert().ID();

        const Eigen::Vector3f e1 = -edge.Next().Dir().normalized();
        const Eigen::Vector3f e2 = edge.Prev().Dir().normalized();

        const float alpha = acos(e1.dot(e2));

        const Eigen::Vector3f p1 = -edge.Pair().Next().Dir().normalized();
        const Eigen::Vector3f p2 = edge.Pair().Prev().Dir().normalized();

        const float beta = acos(p1.dot(p2));

        const float laplace_coeff = -0.5f * ((1.f / tan(alpha)) + (1.f / tan(beta)));

        laplace_operator.insert(i, j) = laplace_coeff;
        laplace_operator.insert(j, i) = laplace_coeff;

        laplace_operator.coeffRef(i, i) -= laplace_coeff;
        laplace_operator.coeffRef(j, j) -= laplace_coeff;
    }

    uint count = 0;
    for(auto& v: mesh.Verts())
    {
        std::vector<HMesh<VertexData>::MFace*> faces = v.ContainingFaces();
        for(auto f: faces)
        {
            areas.diagonal()[count] += f->Area();
        }

        areas.diagonal()[count] /= 3.d;

        count++;
    }

    SparseLU<SparseMatrix<float>, COLAMDOrdering<int> > solver;

    const Eigen::SparseMatrix<float> a =
        Eigen::SparseMatrix<float>(areas) + 0.1 * laplace_operator;
    solver.analyzePattern(a);
    solver.factorize(a);

    Eigen::VectorXf heat = solver.solve(sources);
    gradient = CalculateSimplifiedVertexGradient(mesh, heat);

    Eigen::VectorXf integrated_gradient(gradient.size());
    integrated_gradient.setZero();

    for(const auto face: mesh.Faces())
    {
        using Edge = HMesh<VertexData>::MEdge;
        const Edge& e1 = face.Edge();
        const Edge& e2 = e1.Next();
        const Edge& e3 = e2.Next();

        const Eigen::Vector3f d1 = e1.Dir();
        const Eigen::Vector3f d2 = e2.Dir();
        const Eigen::Vector3f d3 = e3.Dir();

        const Eigen::Vector3f nd1 = d1.normalized();
        const Eigen::Vector3f nd2 = d2.normalized();
        const Eigen::Vector3f nd3 = d3.normalized();

        const float cot_theta1 = 1.f / tan(acos(nd1.dot(-nd3)));
        const float cot_theta2 = 1.f / tan(acos(nd2.dot(-nd1)));
        const float cot_theta3 = 1.f / tan(acos(nd3.dot(-nd2)));

        const uint v1 = e1.Vert().ID();
        const uint v2 = e2.Vert().ID();
        const uint v3 = e3.Vert().ID();

        const auto g1 = -gradient[v1].normalized();
        const auto g2 = -gradient[v2].normalized();
        const auto g3 = -gradient[v3].normalized();

        integrated_gradient[v1] += cot_theta3 * d1.dot(g1) + cot_theta2 * (-d3.dot(g1));
        integrated_gradient[v2] += cot_theta1 * d2.dot(g2) + cot_theta3 * (-d1.dot(g2));
        integrated_gradient[v3] += cot_theta2 * d3.dot(g3) + cot_theta1 * (-d2.dot(g3));
    }

    solver.analyzePattern(-laplace_operator);
    solver.factorize(-laplace_operator);

    return solver.solve(integrated_gradient);
}

For clarity, the calculations I am trying to solve are on the 5th page of the linked paper, i.e these ones:

In particular, solving the systems of equations:

$(M-tL_c)u=\delta_\gamma$

$L_C\phi=b$

Answer 1

The solution to the final Poisson equation is defined only up to an additive constant. So you just need to shift the solution vector so the smallest value is zero.