Graphics processing (vii) Calculation of geodesic distance based on thermal propagation-siggraph 2013_ graphic processing

Here to share with you is the 2013 siggraph above a paper, named "geodesics in heat:a New approach to Computing Distance Based" flow, this Heat does not Provide source code, but because the idea of the algorithm is quite novel, if you have studied other geodesic triangulation algorithms on other geodesic surface, then you will be very excited when you see this paper, and find this algorithm quite magical, and a new breakthrough has been made in the computational method of geodesic distance on the grid surface. Because see this paper is very excited, so I am excited to write the code once, although geodesic distance calculation method, for me, it is not useful, but its idea, its idea is worth us to learn.

The algorithm of geodesics in heat:a the New approach to Computing Distance Based is mainly through the method of vector field to solve the Poisson equation by solving the heat flux. This paper I think should classify it as the literature of vector field type, the application of vector field on surface is very extensive, it can be used for mesh optimization reconstruction, parametric texture mapping, mesh deformation ... And so on, now this article paper again let me know that the vector field can also solve the geodesic distance. Through the heat transmission method, to solve the geodesic distance, it's Daniel.

I. Theoretical knowledge

In the early time for geodesic distances, Daniel deduced the solution of geodesic distances to solve eikonal equation (process function equation):

And the boundary constraints are satisfied with the following conditions:

The symbol φ above is the geodesic distance.

That is, whether it is a European space or a curved surface space, its geodesic distance gradient modulus is equal to 1, and then boundary conditions: The source point of the geodesic distance value of 0.

So many Daniel are studying the solution of the process function equation above, however, the proposed algorithm is non-linear, because the above equation is nonlinear equation, the computational amount is very large on the mesh surface. The real contribution of this paper author is to transform the non-linear problem into a linear one, and finally we need to solve a Poisson equation. Before starting this algorithm, it is recommended to take a look at "mesh Editing with poisson-based gradient Field manipulation," the use of Poisson equation for triangular mesh model Editor, in fact, this paper thought should be borrowed from the " Poisson Image Editing, if you are already familiar with "Mesh Editing with poisson-based gradient Field manipulation" then look at this paper will be more effective.

1, time discretization, the thermal propagation equation time discretization:

2, Space discretization, on the grid surface, the discrete laplace coordinates are defined as:

AI represents one of three parts of the area of a triangular surface, and J represents the adjacency vertex of vertex i. For a mesh model with a V vertex, we can list the V equation above and write out the matrix form:

Where the LC is the Lapura matrix, and A is the diagonal matrix that contains each vertex area.

The surrogate formula (3) may be:

On the surface of a mesh, for a given triangular plane, the gradient formula for a scalar field is as follows:

AF represents the area of a triangle, n represents the normal vector of a triangular plane, and the UI is a scalar field, and we can consider the geodesic distance of each vertex of a mesh surface as a scalar field.

The divergence of vertex i is discrete in the form of:

Lazy to say a lot of nonsense, in short to finally boil down to solve the following equation:

b We need to solve the scalar field U First, then calculate the divergence of the vertex geodesic distance according to the calculation formula of divergence degree.

Second, the implementation of the algorithm

Algorithm Total Flow:

Diagram:

Where φ is the geodesic distance we require, T is an infinitesimal number, reaching 0

The implementation of the algorithm, first of all, I wrote the following code is very messy, lazy to tidy, because I just to learn:

1. Solving U

According to the formula:

To solve u, so we need to figure out A,T,LC, and δ. And then I'm going to explain the solution to four parameters:

A, time t calculation: see the literature 3.2.4, time is theoretically a very small number, the author gives the appropriate value is: the average length of the grid square

The code implementation is as follows:

[CPP] View plain copy//Time step, in the literature the time step is: The edge length average of the grid model, then the square void Cheatgeodesics::set_time () {m_basemesh->need_e       Dge ();       int en=m_basemesh->m_edges.size ();       float sumlength=0;       for (int i=0;i<en;i++) {sumlength+=m_basemesh->m_edges[i].length ();       The average length of the sumlength=sumlength/en;//side is sumlength=sumlength*sumlength;   M_time_step=sumlength; b, calculation of the area diagonal matrix A:

[cpp]  View Plain  copy///Compute the area matrix of Vertex VI in geodesics in heat  literature    void  Cheatgeodesics::get_matrix_a_areas ()    {       m_basemesh->need_ Adjacentfaces ();       gsi::SparseMatrix &A=m_A_areas;        a.resize (M_numberv,m_numberv);       //calculates the area of each triangular slice        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;        }       //calculates the area occupied by each vertex in the Laplace operator, That is, the total area of the adjacent triangular surface One-third        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++)            &nbsP {               sumarea+=face_area[af[j]] ;           }          The area of the   //contains vertices: One-third            sumarea= of the total area of the adjacent triangular patches sumarea/3.0;           a (i,i)  =sumarea;       }  }   C, Laplace matrix LC calculation:

[cpp]  View Plain  copy//The LC matrix in geodesics in heat  literature, i.e. Laplace operator    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 = &NBSP;0;K&LT;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 ()  );            }          }     }  

The cotangent cot weights of adjacent vertices are calculated:

[cpp]  View Plain  copy//adjacency vertex cotangent weight calculation    void cheatgeodesics::cotangentweights (trimesh* Tmesh,int vindex,vector<float>&vweight,float &weightsum,bool bnormalize)// Calculate the respective Cottan weights    {          int neighbornumber= of the first order neighboring points 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;       }  }  // Gets the common adjacency vertex    Void cheatgeodesics::co_neighbor of two vertices (trimesh *tmesh,int u_id,int v_id, VECTOR&LT;INT&GT;&AMP;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]);                }           }        }  }  

Finally solving the equations, you can find the U.

2. Solving the thermal field ▽u (Heat flow) of triangular mesh surfaces:

This step is directly based on the formula:

The solution is OK. Then the ▽u is normalized and the opposite direction of the heat field is obtained, that is, the gradient field of the geodesic distance is obtained:

[cpp]  View Plain  copy 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]);//The important two-step   of the literature is gradient unit and gradient reverse        fgradient[i]=float ( -1.0) *fgradient[i];       //length0[i]=len (Fgradient[i]);  }  
3. To the gradient field X of geodesic distance, the divergence degree is obtained.

According to the formula:

The divergence of the scalar field gradient field is solved by the geodesic distance.

[cpp]  View plain  copy//compute vertex divergence    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);                   //Calculation divergence                    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]];       &Nbsp;           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, the Final solution is to obtain the geodesic distance of each vertex on the mesh surface by solving the Poisson equation.

[cpp]  View Plain  copy Gsi::sparselinearsystem *psystempos2 = new gsi:: Sparselinearsystem ();       gsi::solver_taucs * psolverpos2 =  new gsi::solver_taucs (PSYSTEMPOS2);       psystempos2->resize (m_ Numberv, m_numberv);   &NBSP;&NBSP;&NBSP;&NBSP;PSYSTEMPOS2-&GT;RESIZERHS (1);        //because of the following solution equations the use of solver_taucs::taucs_llt&nbsp, that is, semidefinite matrix mode        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);   &NBsp;   psolverpos2->setsolvermode ( gsi::Solver_TAUCS::TAUCS_LLT );       psolverpos2->setorderingmode ( gsi::Solver_TAUCS::TAUCS_METIS );        gsi::Vector pRHSPos2 ;        Prhspos2.resize (m_numberv);       //  right constraint add         for  (int i=0;i<m_numberv;i++)        { &