我已经实现了测地线的加热方法：

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

当我运行它时，我得到的解决方案在视觉上似乎是正确的：

在这张图片中，箭头是负渐变，颜色是距离，这意味着当一个点的颜色趋向白色时，该点离源顶点更远。

我看到的图像在视觉上是我所期望的，但是我注意到负值。如果我打印每个点的距离的绝对值，而不是距离，我得到：

本质上，靠近源顶点的点正在获得负值。

我已经检查了多次计算这个的代码，据我所知，我做的事情是正确的。这篇论文需要求解 2 个线性系统，这两个系统都可以在论文第 5 页的第一列中看到。

我不完全确定为什么我的一些距离是负的。

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

为清楚起见，我试图解决的计算在链接论文的第 5 页，即这些：

特别是，求解方程组：

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

$L_C\phi=b$