热传播测地线距离的计算

这里要跟大家分享的是2013年Siggraph上面的一篇paper,名为《Geodesics in Heat:A New Approach to Computing Distance Based on Heat Flow》,这篇paper没有提供源代码,但是因为算法的思想相当新颖,如果你之前有研究过其它的测地三角网格曲面上的测地距离算法,那么看到这篇paper后,你会非常的激动,觉得这个算法相当神奇,网格曲面上测地距离的计算方法又有了新的突破。因为看到这篇paper非常激动,以至于我兴奋得马上去把代码写了一遍,虽然测地距离的计算方法,对于我来说没什么用,但它的想法,它的思想值得我们学习。

《Geodesics in Heat:A New Approach to Computing Distance Based on Heat Flow》算法主要是通过向量场的方法,通过求解热流,求解泊松方程。这篇paper我觉得应该把它归类为向量场类型的文献,曲面上向量场的应用非常广泛,可以用于网格优化重建、参数化纹理映射、网格变形……等,如今这篇paper又让我明白向量场也可以求解测地距离。通过热量传播的方法,去求解测地距离,真是大牛啊


一、理论知识

在很早之前对于测地距离,大牛们就推导出了测地距离的求解归结为求解eikonal equation(程函方程):


且满足边界约束条件为:

上面符号Φ就是测地距离。


也就是不管是欧式空间还是曲面空间,其测地距离的梯度模长恒等于1,然后边界条件:源点的测地距离值为0。

因此很多大牛,都是针对上面的程函方程的求解进行研究,然而所提出的算法都是非线性的方法,因为上面的方程就是非线性方程,在网格曲面上计算量非常大。而这篇paper作者的真正贡献在于把非线性问题转换为线性问题,到最后只需要求解一个泊松方程就可以了。开始这个算法之前,建议先看一下《Mesh Editing with Poisson-Based Gradient Field Manipulation》这篇利用泊松方程进行三角网格模型编辑,其实这篇paper的思想应该是借用了《Poisson Image Editing》,如果你已经很熟悉《Mesh Editing with Poisson-Based Gradient Field Manipulation》那么看这篇paper将事半功倍。

1、时间离散化,热传播方程时间离散化为:


2、空间离散化,在网格曲面上,离散的拉普拉斯坐标定义为:



Ai表示三角面片的面积的三份之一,j表示顶点i的邻接顶点。对于有V个顶点的网格模型,我们可以列出上面的V个方程,写出矩阵形式为:

其中,Lc为拉普拉矩阵,A为包含每个顶点面积的对角矩阵。

代入公式(3)可得:


在网格曲面上,对于一个给定的三角面片平面上,标量场的梯度计算公式如下:



Af表示三角形的面积,N表示三角面片的法矢,ui就是标量场,我们可以把网格曲面每个顶点的测地距离看成是一个标量场。

顶点i的散度离散形式为:



懒得多说废话了,总之到最后归结为求解如下方程:


其中b我们需要先通过求解标量场u,然后根据散度的计算公式,计算出顶点测地距离的散度。

二、算法实现

算法总流程:


示意图:


其中,Φ就是我们要求的测地距离,t是一个无穷小的数,趋近于0

算法具体实现,先说明一下,我下面自己写的代码很乱,懒得整理,因为我只是为了学习:

1、求解u

根据公式:


求解u,因此我们需要先算出A,t,Lc,还有δ。接着我将逐渐讲解着四个参数的求解:

a、时间t计算:见文献3.2.4,时间理论上来说是一个非常小的数,文中作者给出了其合适的值为:网格平均边长的平方


代码实现如下:

  1. //计算时间步长,文献中时间步长为:网格模型的边长平均,然后平方  
  2. void CHeatGeodesics::Set_Time()  
  3. {  
  4.     m_BaseMesh->need_edge();  
  5.     int en=m_BaseMesh->m_edges.size();  
  6.     float sumlength=0;  
  7.     for (int i=0;i<en;i++)  
  8.     {  
  9.         sumlength+=m_BaseMesh->m_edges[i].length();  
  10.     }  
  11.     sumlength=sumlength/en;//边长的平均长度  
  12.     sumlength=sumlength*sumlength;  
  13.     m_Time_Step=sumlength;  
  14. }  
//计算时间步长,文献中时间步长为:网格模型的边长平均,然后平方
void CHeatGeodesics::Set_Time()
{
    m_BaseMesh->need_edge();
    int en=m_BaseMesh->m_edges.size();
    float sumlength=0;
    for (int i=0;i<en;i++)
    {
        sumlength+=m_BaseMesh->m_edges[i].length();
    }
    sumlength=sumlength/en;//边长的平均长度
    sumlength=sumlength*sumlength;
    m_Time_Step=sumlength;
}
b、面积对角矩阵A的计算:

  1. //计算Geodesics in Heat 文献中的包含顶点Vi的面积矩阵  
  2. void CHeatGeodesics::Get_Matrix_A_areas()  
  3. {  
  4.     m_BaseMesh->need_adjacentfaces();  
  5.     gsi::SparseMatrix &A=m_A_areas;  
  6.     A.Resize(m_NumberV,m_NumberV);  
  7.     //计算每个三角面片的面积  
  8.     int fn=m_BaseMesh->faces.size();  
  9.     vector<float>&face_area=m_Faces_Area;  
  10.     face_area.resize(fn);  
  11.     for (int i=0;i<fn;i++)  
  12.     {  
  13.         TriMesh::Face &f=m_BaseMesh->faces[i];  
  14.         vec vij=m_BaseMesh->vertices[f[1]]-m_BaseMesh->vertices[f[0]];  
  15.         vec vik=m_BaseMesh->vertices[f[2]]-m_BaseMesh->vertices[f[0]];  
  16.         float areaf=0.5*len(vij CROSS vik);  
  17.         face_area[i]=areaf;  
  18.     }  
  19.     //计算拉普拉斯算子中每个顶点所占据的面积,即邻接三角面片面积总和的三分之一  
  20.     for (int i=0;i<m_NumberV;i++)  
  21.     {  
  22.         vector<int>&af=m_BaseMesh->adjacentfaces[i];  
  23.         int n_af=af.size();  
  24.         float sumarea=0.0;  
  25.         for (int j=0;j<n_af;j++)  
  26.         {  
  27.             sumarea+=face_area[af[j]];  
  28.         }  
  29.         //包含顶点的面积为:邻接三角面片面积总和的三分之一  
  30.         sumarea=sumarea/3.0;  
  31.         A(i,i) =sumarea;  
  32.     }  
  33. }  
//计算Geodesics in Heat 文献中的包含顶点Vi的面积矩阵
void CHeatGeodesics::Get_Matrix_A_areas()
{
    m_BaseMesh->need_adjacentfaces();
    gsi::SparseMatrix &A=m_A_areas;
    A.Resize(m_NumberV,m_NumberV);
    //计算每个三角面片的面积
    int fn=m_BaseMesh->faces.size();
    vector<float>&face_area=m_Faces_Area;
    face_area.resize(fn);
    for (int i=0;i<fn;i++)
    {
        TriMesh::Face &f=m_BaseMesh->faces[i];
        vec vij=m_BaseMesh->vertices[f[1]]-m_BaseMesh->vertices[f[0]];
        vec vik=m_BaseMesh->vertices[f[2]]-m_BaseMesh->vertices[f[0]];
        float areaf=0.5*len(vij CROSS vik);
        face_area[i]=areaf;
    }
    //计算拉普拉斯算子中每个顶点所占据的面积,即邻接三角面片面积总和的三分之一
    for (int i=0;i<m_NumberV;i++)
    {
        vector<int>&af=m_BaseMesh->adjacentfaces[i];
        int n_af=af.size();
        float sumarea=0.0;
        for (int j=0;j<n_af;j++)
        {
            sumarea+=face_area[af[j]];
        }
        //包含顶点的面积为:邻接三角面片面积总和的三分之一
        sumarea=sumarea/3.0;
        A(i,i) =sumarea;
    }
}
C、拉普拉斯矩阵Lc计算:

  1. //计算Geodesics in Heat 文献中的Lc矩阵,即拉普拉斯算子  
  2. void CHeatGeodesics::Get_Matrix_Lc()  
  3. {  
  4.     Ls.resize(m_NumberV,m_NumberV);  
  5.     int m_nEdges=10000 ;  
  6.     Ls.reserve(m_nEdges+m_NumberV);  
  7.     for (int i = 0;i<m_NumberV;++i)   
  8.     {  
  9.         VProperty & vi = m_vertices[i];  
  10.         int nNbrs = vi.VNeighbors.size();  
  11.         for (int k = 0;k<nNbrs;++k)   
  12.         {  
  13.             Ls.insert(i, vi.VNeighbors[k]) = vi.VNeiWeight[k];  
  14.         }  
  15.         Ls.insert(i, i) = -vi.VSumWeight;  
  16.     }  
  17.     Ls.finalize();  
  18.     gsi::SparseMatrix &A=m_Lc;  
  19.     A.Resize(m_NumberV,m_NumberV);  
  20.     for(int k=0;k<m_NumberV;++k)   
  21.     {  
  22.         for ( SparseMatrixType::InnerIterator it(Ls,k); it; ++it)  
  23.         {  
  24.             A.Set( it.row(), it.col(), it.value() );  
  25.         }     
  26.     }  
  27.   
  28. }  
//计算Geodesics in Heat 文献中的Lc矩阵,即拉普拉斯算子
void CHeatGeodesics::Get_Matrix_Lc()
{
    Ls.resize(m_NumberV,m_NumberV);
    int m_nEdges=10000 ;
    Ls.reserve(m_nEdges+m_NumberV);
    for (int i = 0;i<m_NumberV;++i) 
    {
        VProperty & vi = m_vertices[i];
        int nNbrs = vi.VNeighbors.size();
        for (int k = 0;k<nNbrs;++k) 
        {
            Ls.insert(i, vi.VNeighbors[k]) = vi.VNeiWeight[k];
        }
        Ls.insert(i, i) = -vi.VSumWeight;
    }
    Ls.finalize();
    gsi::SparseMatrix &A=m_Lc;
    A.Resize(m_NumberV,m_NumberV);
    for(int k=0;k<m_NumberV;++k) 
    {
        for ( SparseMatrixType::InnerIterator it(Ls,k); it; ++it)
        {
            A.Set( it.row(), it.col(), it.value() );
        }   
    }

}

邻接顶点的余切cot权重计算:

  1. //邻接顶点的余切权重计算  
  2. void CHeatGeodesics::CotangentWeights(TriMesh*TMesh,int vIndex,vector<float>&vweight,float &WeightSum,bool bNormalize)//计算一阶邻近点的各自cottan权重  
  3. {     
  4.     int NeighborNumber=TMesh->neighbors[vIndex].size();  
  5.     vweight.resize(NeighborNumber);  
  6.     WeightSum=0;  
  7.     vector<int>&NeiV=TMesh->neighbors[vIndex];  
  8.     for (int i=0;i<NeighborNumber;i++)  
  9.     {  
  10.         int j_nei=NeiV[i];  
  11.         vector<int>tempnei;  
  12.         Co_neighbor(TMesh,vIndex,j_nei,tempnei);  
  13.         float cotsum=0.0;  
  14.         for (int j=0;j<tempnei.size();j++)  
  15.         {  
  16.             vec vivo=TMesh->vertices[vIndex]-TMesh->vertices[tempnei[j]];  
  17.             vec vjvo=TMesh->vertices[j_nei]-TMesh->vertices[tempnei[j]];  
  18.             float dotvector=vivo DOT vjvo;  
  19.             dotvector=dotvector/sqrt(len2(vivo)*len2(vjvo)-dotvector*dotvector);  
  20.             cotsum+=dotvector;  
  21.         }  
  22.         vweight[i]=cotsum/2.0;  
  23.         WeightSum+=vweight[i];  
  24.     }  
  25.   
  26.     if ( bNormalize )   
  27.     {  
  28.         for (int k=0;k<NeighborNumber;++k)  
  29.         {  
  30.             vweight[k]/=WeightSum;  
  31.         }  
  32.         WeightSum=1.0;  
  33.     }  
  34. }  
  35.   
  36. void CHeatGeodesics:: UniformWeights(TriMesh*TMesh,int vIndex,vector<float>&vweight,float &WeightSum,bool bNormalize)  
  37. {  
  38.     int NeighborNumber=TMesh->neighbors[vIndex].size();  
  39.     vweight.resize(NeighborNumber);  
  40.     WeightSum = 0;  
  41.     for (int j = 0; j <NeighborNumber; ++j )  
  42.     {  
  43.         vweight[j] = 1;  
  44.         WeightSum += vweight[j];  
  45.     }  
  46.   
  47.     if ( bNormalize )  
  48.     {  
  49.         for ( int k = 0; k < NeighborNumber; ++k )  
  50.             vweight[k] /= WeightSum;  
  51.         WeightSum=1.0;  
  52.     }  
  53. }  
  54. //获取两顶点的共同邻接顶点  
  55. void CHeatGeodesics::Co_neighbor(TriMesh *Tmesh,int u_id,int v_id,vector<int>&co_neiv)  
  56. {  
  57.     Tmesh->need_adjacentedges();  
  58.     vector<int>&u_id_ae=Tmesh->adjancetedge[u_id];   
  59.     int en=u_id_ae.size();  
  60.     Tedge Co_Edge;  
  61.     for (int i=0;i<en;i++)  
  62.     {  
  63.         Tedge &ae=Tmesh->m_edges[u_id_ae[i]];  
  64.         int opsi=ae.opposite_vertex(u_id);  
  65.         if (opsi==v_id)  
  66.         {  
  67.             Co_Edge=ae;  
  68.             break;  
  69.         }  
  70.     }  
  71.     for (int i=0;i<Co_Edge.m_adjacent_faces.size();i++)  
  72.     {  
  73.         TriMesh::Face af=Co_Edge.m_adjacent_faces[i];  
  74.         for (int j=0;j<3;j++)  
  75.         {  
  76.             if((af[j]!=u_id)&&(af[j]!=v_id))  
  77.             {  
  78.                 co_neiv.push_back(af[j]);  
  79.             }  
  80.         }  
  81.     }  
  82. }  
//邻接顶点的余切权重计算
void CHeatGeodesics::CotangentWeights(TriMesh*TMesh,int vIndex,vector<float>&vweight,float &WeightSum,bool bNormalize)//计算一阶邻近点的各自cottan权重
{   
    int NeighborNumber=TMesh->neighbors[vIndex].size();
    vweight.resize(NeighborNumber);
    WeightSum=0;
    vector<int>&NeiV=TMesh->neighbors[vIndex];
    for (int i=0;i<NeighborNumber;i++)
    {
        int j_nei=NeiV[i];
        vector<int>tempnei;
        Co_neighbor(TMesh,vIndex,j_nei,tempnei);
        float cotsum=0.0;
        for (int j=0;j<tempnei.size();j++)
        {
            vec vivo=TMesh->vertices[vIndex]-TMesh->vertices[tempnei[j]];
            vec vjvo=TMesh->vertices[j_nei]-TMesh->vertices[tempnei[j]];
            float dotvector=vivo DOT vjvo;
            dotvector=dotvector/sqrt(len2(vivo)*len2(vjvo)-dotvector*dotvector);
            cotsum+=dotvector;
        }
        vweight[i]=cotsum/2.0;
        WeightSum+=vweight[i];
    }

    if ( bNormalize ) 
    {
        for (int k=0;k<NeighborNumber;++k)
        {
            vweight[k]/=WeightSum;
        }
        WeightSum=1.0;
    }
}

void CHeatGeodesics:: UniformWeights(TriMesh*TMesh,int vIndex,vector<float>&vweight,float &WeightSum,bool bNormalize)
{
    int NeighborNumber=TMesh->neighbors[vIndex].size();
    vweight.resize(NeighborNumber);
    WeightSum = 0;
    for (int j = 0; j <NeighborNumber; ++j )
    {
        vweight[j] = 1;
        WeightSum += vweight[j];
    }

    if ( bNormalize )
    {
        for ( int k = 0; k < NeighborNumber; ++k )
            vweight[k] /= WeightSum;
        WeightSum=1.0;
    }
}
//获取两顶点的共同邻接顶点
void CHeatGeodesics::Co_neighbor(TriMesh *Tmesh,int u_id,int v_id,vector<int>&co_neiv)
{
    Tmesh->need_adjacentedges();
    vector<int>&u_id_ae=Tmesh->adjancetedge[u_id]; 
    int en=u_id_ae.size();
    Tedge Co_Edge;
    for (int i=0;i<en;i++)
    {
        Tedge &ae=Tmesh->m_edges[u_id_ae[i]];
        int opsi=ae.opposite_vertex(u_id);
        if (opsi==v_id)
        {
            Co_Edge=ae;
            break;
        }
    }
    for (int i=0;i<Co_Edge.m_adjacent_faces.size();i++)
    {
        TriMesh::Face af=Co_Edge.m_adjacent_faces[i];
        for (int j=0;j<3;j++)
        {
            if((af[j]!=u_id)&&(af[j]!=v_id))
            {
                co_neiv.push_back(af[j]);
            }
        }
    }
}

最后求解方程组,就可以把u求出来了。

2、求解三角网格曲面的热量场▽u (Heat flow):

这一步直接根据公式:


求解就可以了。然后对▽u进行归一化,并取热量场的反方向,即求测地距离的梯度场:

  1. for (int i=0;i<fn;i++)  
  2. {  
  3.     TriMesh::Face &f=m_BaseMesh->faces[i];  
  4.     for (int j=0;j<3;j++)  
  5.     {  
  6.         vec ei=m_BaseMesh->vertices[f[(j+2)%3]]-m_BaseMesh->vertices[f[(j+1)%3]];  
  7.         vec FNormal=normalize(m_BaseMesh->FaceNormal[i]);  
  8.         vec gradient=float(m_vertices[f[j]].VU*0.5/m_Faces_Area[i])*(FNormal CROSS ei);  
  9.         FGradient[i]=FGradient[i]+gradient;  
  10.     }  
  11.     normalize(FGradient[i]);//文献的重要两步 即梯度单位化和梯度反向  
  12.     FGradient[i]=float(-1.0)*FGradient[i];  
  13.     //length0[i]=len(FGradient[i]);  
  14. }  
    for (int i=0;i<fn;i++)
    {
        TriMesh::Face &f=m_BaseMesh->faces[i];
        for (int j=0;j<3;j++)
        {
            vec ei=m_BaseMesh->vertices[f[(j+2)%3]]-m_BaseMesh->vertices[f[(j+1)%3]];
            vec FNormal=normalize(m_BaseMesh->FaceNormal[i]);
            vec gradient=float(m_vertices[f[j]].VU*0.5/m_Faces_Area[i])*(FNormal CROSS ei);
            FGradient[i]=FGradient[i]+gradient;
        }
        normalize(FGradient[i]);//文献的重要两步 即梯度单位化和梯度反向
        FGradient[i]=float(-1.0)*FGradient[i];
        //length0[i]=len(FGradient[i]);
    }

3、对测地距离的梯度场X,求取散度。

根据公式:


求解测地距离标量场梯度场的散度。

  1. //计算顶点的散度  
  2. for (int i=0;i<vn;i++)  
  3. {     
  4.     vector<int>&adjacentface=m_BaseMesh->adjacentfaces[i];  
  5.     for (int j=0;j<adjacentface.size();j++)  
  6.     {  
  7.         TriMesh::Face &f=m_BaseMesh->faces[adjacentface[j]];  
  8.         for (int k=0;k<3;k++)  
  9.         {  
  10.             if (f[k]==i)  
  11.             {  
  12.                 vec ei=m_BaseMesh->vertices[f[(k+2)%3]]-m_BaseMesh->vertices[f[(k+1)%3]];  
  13.                 /*vec gradient=float(0.5/m_Faces_Area[adjacentface[j]])*(m_BaseMesh->FaceNormal[adjacentface[j]] CROSS ei); 
  14.                 //计算散度 
  15.                 m_vertices[i].VDivergence+=(FGradient[adjacentface[j]] DOT gradient)*m_Faces_Area[adjacentface[j]];*/  
  16.   
  17.                 vec e1=m_BaseMesh->vertices[f[(k+1)%3]]-m_BaseMesh->vertices[f[k]];  
  18.                 vec e2=m_BaseMesh->vertices[f[(k+2)%3]]-m_BaseMesh->vertices[f[k]];  
  19.                 float cot_angle1=Cot_angle(e2,ei);  
  20.                 float cot_angle2=Cot_angle(float(-1.0)*e1,ei);  
  21.             m_vertices[i].VDivergence+=0.5*(cot_angle1*(e1 DOT FGradient[adjacentface[j]])+cot_angle2*(e2 DOT FGradient[adjacentface[j]]));  
  22.                 break;  
  23.             }  
  24.         }  
  25.     }  
  26.   
  27. }  
   //计算顶点的散度
    for (int i=0;i<vn;i++)
    {   
        vector<int>&adjacentface=m_BaseMesh->adjacentfaces[i];
        for (int j=0;j<adjacentface.size();j++)
        {
            TriMesh::Face &f=m_BaseMesh->faces[adjacentface[j]];
            for (int k=0;k<3;k++)
            {
                if (f[k]==i)
                {
                    vec ei=m_BaseMesh->vertices[f[(k+2)%3]]-m_BaseMesh->vertices[f[(k+1)%3]];
                    /*vec gradient=float(0.5/m_Faces_Area[adjacentface[j]])*(m_BaseMesh->FaceNormal[adjacentface[j]] CROSS ei);
                    //计算散度
                    m_vertices[i].VDivergence+=(FGradient[adjacentface[j]] DOT gradient)*m_Faces_Area[adjacentface[j]];*/

                    vec e1=m_BaseMesh->vertices[f[(k+1)%3]]-m_BaseMesh->vertices[f[k]];
                    vec e2=m_BaseMesh->vertices[f[(k+2)%3]]-m_BaseMesh->vertices[f[k]];
                    float cot_angle1=Cot_angle(e2,ei);
                    float cot_angle2=Cot_angle(float(-1.0)*e1,ei);
                m_vertices[i].VDivergence+=0.5*(cot_angle1*(e1 DOT FGradient[adjacentface[j]])+cot_angle2*(e2 DOT FGradient[adjacentface[j]]));
                    break;
                }
            }
        }

    }
4、最后求解通过求解泊松方程,求取网格曲面上每个顶点的测地距离。

  1. gsi::SparseLinearSystem *pSystemPos2 = new gsi::SparseLinearSystem();  
  2.     gsi::Solver_TAUCS * pSolverPos2 = new gsi::Solver_TAUCS(pSystemPos2);  
  3.     pSystemPos2->Resize(m_NumberV, m_NumberV);  
  4.     pSystemPos2->ResizeRHS(1);  
  5.     //由于以下解方程组使用Solver_TAUCS::TAUCS_LLT ,即半正定矩阵模式  
  6.     m_Lc.Multiply(-1.0);  
  7.   
  8.     m_Lc(0,0) =m_Lc(0,0) +10;  
  9.   
  10.     pSystemPos2->SetMatrix(m_Lc);  
  11.     if ( ! pSystemPos2->Matrix().IsSymmetric() )assert(0);  
  12.     pSolverPos2->OnMatrixChanged();  
  13.     pSolverPos2->SetStoreFactorization(true);  
  14.     pSolverPos2->SetSolverMode( gsi::Solver_TAUCS::TAUCS_LLT );  
  15.     pSolverPos2->SetOrderingMode( gsi::Solver_TAUCS::TAUCS_METIS );  
  16.     gsi::Vector pRHSPos2 ;  
  17.     pRHSPos2.Resize(m_NumberV);  
  18.     // 右边项约束添加  
  19.     for (int i=0;i<m_NumberV;i++)  
  20.     {  
  21.         pRHSPos2[i]=float(-1.0)*m_vertices[i].VDivergence;    
  22.     }  
  23.   
  24.      //初始条件为方程热源的测地距离为0  
  25.     pRHSPos2[0]=pRHSPos2[0]+10;  
  26.       
  27.     pSystemPos2->SetRHS(0, pRHSPos2);  
  28.     bool boksoveLs2=pSolverPos2->Solve();  
  29.     if (!boksoveLs2)assert(0);  
  30.     m_GeodesicsDistance.resize(m_NumberV);  
  31.     for ( int i=0;i<m_NumberV;++i)   
  32.     {     
  33.        m_GeodesicsDistance[i]=(float)pSystemPos2->GetSolution(i,0);  
  34.     }  
gsi::SparseLinearSystem *pSystemPos2 = new gsi::SparseLinearSystem();
    gsi::Solver_TAUCS * pSolverPos2 = new gsi::Solver_TAUCS(pSystemPos2);
    pSystemPos2->Resize(m_NumberV, m_NumberV);
    pSystemPos2->ResizeRHS(1);
    //由于以下解方程组使用Solver_TAUCS::TAUCS_LLT ,即半正定矩阵模式
    m_Lc.Multiply(-1.0);

    m_Lc(0,0) =m_Lc(0,0) +10;

    pSystemPos2->SetMatrix(m_Lc);
    if ( ! pSystemPos2->Matrix().IsSymmetric() )assert(0);
    pSolverPos2->OnMatrixChanged();
    pSolverPos2->SetStoreFactorization(true);
    pSolverPos2->SetSolverMode( gsi::Solver_TAUCS::TAUCS_LLT );
    pSolverPos2->SetOrderingMode( gsi::Solver_TAUCS::TAUCS_METIS );
    gsi::Vector pRHSPos2 ;
    pRHSPos2.Resize(m_NumberV);
    // 右边项约束添加
    for (int i=0;i<m_NumberV;i++)
    {
        pRHSPos2[i]=float(-1.0)*m_vertices[i].VDivergence;  
    }

     //初始条件为方程热源的测地距离为0
    pRHSPos2[0]=pRHSPos2[0]+10;

    pSystemPos2->SetRHS(0, pRHSPos2);
    bool boksoveLs2=pSolverPos2->Solve();
    if (!boksoveLs2)assert(0);
    m_GeodesicsDistance.resize(m_NumberV);
    for ( int i=0;i<m_NumberV;++i) 
    {   
       m_GeodesicsDistance[i]=(float)pSystemPos2->GetSolution(i,0);
    }

可以说到了这里算法已经结束了,当然paper后面还有后续的处理,比如距离光顺,还有边界条件的等问题。但都属于后处理,你如果要实现跟paper一模一样的效果,还是得好好看一看后面的边界条件问题。


边界条件对于结果的影响还是蛮大的。

OK,做一下简单的总结:算法整个过程就说白了就是求解一个泊松方程,说的更白一点,就是求解一个大型的方程组AX=b,因此我们的目标就是先算出A、b,其中对于一个给定的网格曲面模型来说,A就是拉普拉斯矩阵,是固定的。b的求解就是整个算法成功的关键。本文地址:http://blog.csdn.net/hjimce/article/details/46415499     作者:hjimce     联系qq:1393852684   更多资源请关注我的博客:http://blog.csdn.net/hjimce                原创文章,转载请保留本行信息。

参考文献:

1、《Geodesics in Heat:A New Approach to Computing Distance Based on Heat Flow》

2、《Mesh Editing with Poisson-Based Gradient Field Manipulation》

3、《Poisson Image Editing》

  • 4
    点赞
  • 16
    收藏
    觉得还不错? 一键收藏
  • 1
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值