设为首页收藏本站

中国膜结构网

 找回密码
 立即注册

QQ登录

只需一步,快速开始

膜结构车棚
膜结构车棚膜结构资质国产膜材 膜结构网中国膜结构协会
查看: 67|回复: 7

热传播测地线距离的计算

[复制链接]
  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

    发表于 2021-5-22 16:08 | 显示全部楼层 |阅读模式
    这里要跟大家分享的是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又让我明白向量场也可以求解测地距离。通过热量传播的方法,去求解测地距离,真是大牛啊
    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. }
    复制代码
    回复


    http://www.mjgw.org/ 专业从事膜结构设计、制作加工、施工安装的膜结构工程服务,能够为客户提供专业的膜结构整体解决方案。做中国最好的膜结构综合服务平台。欢迎大家联系电话:198-7840-1958,QQ:463017170.
    相关关键词:膜结构车棚,膜结构车棚覆盖,膜结构车棚公司,膜结构车棚多少钱,膜结构车棚厂家,膜结构车棚价格,社区膜结构车棚,膜结构车棚膜布厂家 ,膜结构车棚哪家好,膜结构车棚多少钱一米,膜结构车棚报价,膜结构车棚哪里有,膜结构车棚定制,膜结构车棚安装,膜结构车棚设计,膜结构车棚电话,膜结构车棚加工,膜结构车棚膜布价格,膜结构车棚批发,膜结构车棚制造商,膜结构车棚生产厂家,膜结构车棚设计,膜结构车棚施工,膜结构车棚多少钱一平米,膜结构车棚订制,张拉膜车棚,张拉膜车棚覆盖,张拉膜车棚公司,张拉膜车棚多少钱,张拉膜车棚厂家,张拉膜车棚价格,社区张拉膜车棚,张拉膜车棚膜布厂家 ,张拉膜车棚哪家好,张拉膜车棚多少钱一米,张拉膜车棚报价,张拉膜车棚哪里有,张拉膜车棚定制,张拉膜车棚安装,张拉膜车棚设计,张拉膜车棚电话,张拉膜车棚加工,张拉膜车棚膜布价格,张拉膜车棚批发,张拉膜车棚制造商,张拉膜车棚生产厂家,张拉膜车棚设计,张拉膜车棚施工,张拉膜车棚多少钱一平米,张拉膜车棚订制,常用膜材品牌:德国杜肯、法国法拉利、德国海德斯、德国米乐、日本平岗、韩国秀博、比利时希运、美国赫虏伯、中国科宝、上海慧遥。

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2021-5-22 16:08 | 显示全部楼层
    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. }  
    复制代码
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2021-5-22 16:09 | 显示全部楼层
    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. }  
    复制代码
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2021-5-22 16:10 | 显示全部楼层
    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. }  
    复制代码
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2021-5-22 16:10 | 显示全部楼层
    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. }  
    复制代码
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2021-5-22 16:11 | 显示全部楼层
    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. }  
    复制代码
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2021-5-22 16:11 | 显示全部楼层
    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.     }
    复制代码
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    开心
    2021-6-20 09:04
  • 签到天数: 1540 天

    [LV.Master]伴坛终老

     楼主| 发表于 2021-5-22 16:11 | 显示全部楼层
    1. https://blog.csdn.net/yrred_rock/article/details/77622969?utm_medium=distribute.pc_relevant.none-task-blog-baidujs_title-1&spm=1001.2101.3001.4242
    复制代码
    回复 支持 反对

    使用道具 举报

    您需要登录后才可以回帖 登录 | 立即注册

    本版积分规则

    关闭

    推荐膜材品牌上一条 /6 下一条

    进口膜材 国产膜材 pvdf膜材ptfe膜材ETFE膜材
    最好的膜结构公司 一级膜结构资质 膜结构一级资质
    膜结构设计-膜结构十大品牌-etfe设计-充气膜结构
    诺科膜结构
    遨都膜结构设计
    中国膜结构网
    中国空间膜结构

    QQ|申请友链|手机版|中国膜结构论坛