I have implemented the heat method for geodesics:


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

enter image description here

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:

enter image description here

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;


    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;


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

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

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

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

    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));


    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:

enter image description here

In particular, solving the systems of equations:



  • $\begingroup$ It is a bit difficult to help you since you are showing us a code instead of the calculations that you are doing. $\endgroup$
    – nicoguaro
    Apr 19, 2021 at 18:51
  • $\begingroup$ Thank you for the feedback, I made an edit to try to make it easier to understand. $\endgroup$
    – Makogan
    Apr 19, 2021 at 19:33
  • $\begingroup$ Just curious: are there situations where alpha or beta can be zero? If so, laplace_coeff could have a very funky value. Perhaps this occurs near the source vertex? $\endgroup$
    – Charlie S
    Apr 19, 2021 at 20:09
  • $\begingroup$ No, triangles are guaranteed to not be degenerate so alpha and beta are also never zero. $\endgroup$
    – Makogan
    Apr 19, 2021 at 20:11

1 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.