# Solving geodesics on triangular meshes gives negative distances

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

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

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

}


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$$

• It is a bit difficult to help you since you are showing us a code instead of the calculations that you are doing. Apr 19, 2021 at 18:51
• Thank you for the feedback, I made an edit to try to make it easier to understand. Apr 19, 2021 at 19:33
• 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? Apr 19, 2021 at 20:09
• No, triangles are guaranteed to not be degenerate so alpha and beta are also never zero. Apr 19, 2021 at 20:11