图形处理（二）固定边界参数化

固定边界参数化方法是参数化方法中的一类经典算法，至今还有很广泛的用途。这类算法可以说是我读研阶段写的第一个算法。当年年少无知，连外文文献怎么阅读都不懂，导师发给了我好几篇paper，没有一篇看得懂，就连三角网格模型的拓扑邻接关系都不懂。参数化国内相关的硕士、博士论文非常多，所以我就从国内文献开始看起，看了20篇paper才完全知道代码要怎么写，要怎么把文献转换成代码，看了将近一个月的文献，又花了一周的时间把代码实现出来。

现在回想起来，其实这个算法只要一天的时间就可以搞定，因为算法就只只需要求解一个方程组AX=b，因此我们只需要把A、b的计算方法，写出来，然后求解方程组，就可以求出每个顶点的参数化坐标。

1、边界点的参数化。边界的参数化形状常用的有两种，当然其实任意凸多边形都是可以的。

(1)参数化到正方形上：对边界点按逆时针方向弦长参数化到正方形上，参数化后的坐标作为求解内点坐标的已知条件。

(2)参数化到圆周上：对边界点按同样的步骤，参数化到圆周上。



固定边界的参数化方法，根据权值的不同又有好个不同的方法，但是万变不离其中，只是更改一下权值而已。

2、内点的参数化（非边界点）。

内点坐标的计算是根据如下计算公式：



(1)均匀参数化所采用的权值：



均匀参数化是最早的固定边界参数化算法了。均匀参数化的基本思想是把网格的边界点映射到一个平面上的凸多边形,而内点坐标为其一环邻近点的凸组合，其权值取以它为中心点的一环的平均值。



(2)保形参数化。保形参数化是在Tutte的基础上根据凸组合的思想提出的一种方法。此方法保持了每个径向弧长不变，要求所有的三角片都是非退化的。在满足以下两个式子的前提下，其各邻接顶点的权值计算公式如下：









从上面三个公式，其实说白了就是保证边长长度不变，同时角度比例不变的性质。

(3)中值坐标参数化。其权值λ的计算公式为：



(4)调和映射。其权值计算公式如下：



2、算法代码实现

这个代码是我写的第一个网格曲面的代码，已经是3年前的事，所以代码非常粗糙，一大堆bug，所以将就着看看实现：

#include "Parameterization.h"
#include <algorithm>
#include <assert.h>
#include <math.h>

CParameterization::CParameterization(int m_SelectType,bool ParameterRegionType,TriMesh *tmesh)
{
m_ParameterType=m_SelectType;//参数化类型
m_mesh=tmesh;
seedInnerv=-1;
m_SumBorderLength=0.0;
m_ParameterRegionType=ParameterRegionType;
m_ParameterResult.resize(tmesh->vertices.size());
}

CParameterization::~CParameterization(void)
{

}
//提取框选中三角面片的边界顶点、内点索引，并存于m_BorderVertex、m_InnerVertex中
void CParameterization::FindBorderInnerVertex()
{
int i,j,k;
m_mesh->need_adjacentfaces();
m_mesh->need_neighbors();
int n=m_mesh->faces.size();
//调试专用
/*m_mesh->faces[229].beSelect=true;
m_mesh->faces[248].beSelect=true;
m_mesh->faces[249].beSelect=true;
m_mesh->faces[262].beSelect=true;
m_mesh->faces[263].beSelect=true;
m_mesh->faces[283].beSelect=true;*/

/*m_mesh->faces[229].beSelect=true;
m_mesh->faces[248].beSelect=true;
m_mesh->faces[249].beSelect=true;
m_mesh->faces[262].beSelect=true;
m_mesh->faces[263].beSelect=true;
m_mesh->faces[268].beSelect=true;
m_mesh->faces[269].beSelect=true;
m_mesh->faces[282].beSelect=true;
m_mesh->faces[283].beSelect=true;
m_mesh->faces[303].beSelect=true;*/

//寻找到一个内点，作为种子点
for (i=0;i<n;i++)
if (m_mesh->faces[i].beSelect)
{
for (j=0;j<3;j++)
{
const vector<int>&a=m_mesh->adjacentfaces[m_mesh->faces[i][j]];
for(k=0;k<a.size();k++) if (m_mesh->faces[a[k]].beSelect==false)  break;

if(k==a.size()) {seedInnerv=m_mesh->faces[i][j];break;}

}
if (seedInnerv!=-1)break;//内点种子点找到
}

if(i==n)  {AfxMessageBox("未选中任何区域");return;}

m_InteriorVertex.push_back(seedInnerv);
vector<int> nearby;
vector<bool> beSearch(m_mesh->vertices.size(),false);
beSearch[seedInnerv]=TRUE;
// add all the neighbors of the starting vertex into nearby
for(j=0;j< m_mesh->neighbors[seedInnerv].size();j++)
{
int nid=m_mesh->neighbors[seedInnerv][j];
nearby.push_back(nid);
}
//广度优先遍历
while(nearby.size()>0)
{

int iNearby,iNeighbor;
for(i=0; i<nearby.size(); i++)
{
iNearby = nearby[i];
if(beSearch[iNearby])
{
vector<int>::iterator iter;
iter=find(nearby.begin(),nearby.end(),iNearby);
nearby.erase(iter);
continue;
}
vector<int> CountAdjacentSelFaces;
for(j=0;j<m_mesh->adjacentfaces[iNearby].size();j++)
{
if(m_mesh->faces[m_mesh->adjacentfaces[iNearby][j]].beSelect==TRUE)
CountAdjacentSelFaces.push_back(m_mesh->adjacentfaces[iNearby][j]);
/*  if(m_mesh->faces[m_mesh->adjacentfaces[iNearby][j]].beSelect==false)
{
m_BorderVertex.push_back(iNearby);
break;
}//边界点
else if (m_mesh->is_bdy(iNearby))
{
m_BorderVertex.push_back(iNearby);
break;
}*/
}
if(CountAdjacentSelFaces.size()==1)
{
m_mesh->faces[CountAdjacentSelFaces[0]].beSelect=false;
}
else if(CountAdjacentSelFaces.size()<m_mesh->adjacentfaces[iNearby].size())
{
m_BorderVertex.push_back(iNearby);
}
else if(CountAdjacentSelFaces.size()==m_mesh->adjacentfaces[iNearby].size())
{
if(m_mesh->is_bdy(iNearby))
{
m_BorderVertex.push_back(iNearby);
}
else m_InteriorVertex.push_back(iNearby);
}//内点
beSearch[iNearby]=TRUE;
vector<int>::iterator iter;
iter=find(nearby.begin(),nearby.end(),iNearby);
nearby.erase(iter);
for(k=0;k<m_mesh->neighbors[iNearby].size();k++)
{
int h1,h2;
h2=0;
iNeighbor = m_mesh->neighbors[iNearby][k];
if(beSearch[iNeighbor]) continue;
for (h1=0;h1<m_mesh->adjacentfaces[iNeighbor].size();h1++)
{
if(m_mesh->faces[m_mesh->adjacentfaces[iNeighbor][h1]].beSelect==false)
h2++;
}
if (h2<m_mesh->adjacentfaces[iNeighbor].size())nearby.push_back(iNeighbor);
}
}//for
}//while
//边界点排序

vector<int>BorderVertex;
vector<int>::iterator iter;
vector<int>::iterator iter0;
int fid;
if (m_BorderVertex.size()<3) AfxMessageBox("边界点出错");
BorderVertex.push_back(m_BorderVertex[0]);
iter=find(m_BorderVertex.begin(),m_BorderVertex.end(),m_BorderVertex[0]);
m_BorderVertex.erase(iter);
while(m_BorderVertex.size()>0)
{

for (j=0;j<m_mesh->neighbors[BorderVertex[BorderVertex.size()-1]].size();j++)
{
fid=m_mesh->neighbors[BorderVertex[BorderVertex.size()-1]][j];
iter=find(m_BorderVertex.begin(),m_BorderVertex.end(),fid);
if (iter!=m_BorderVertex.end()) break;

}
if (j==m_mesh->neighbors[BorderVertex[BorderVertex.size()-1]].size())
{
int SpecialVertex=BorderVertex[BorderVertex.size()-1];
for(int k=0;k<m_mesh->adjacentfaces[SpecialVertex].size();k++)
{
m_mesh->faces[m_mesh->adjacentfaces[SpecialVertex][k]].beSelect=false;
}
BorderVertex.pop_back();
//AfxMessageBox("参数化失败");
continue;
}
BorderVertex.push_back(m_mesh->neighbors[BorderVertex[BorderVertex.size()-1]][j]);
m_BorderVertex.erase(iter);
}

m_BorderVertex=BorderVertex;

}
//边界顶点参数化
void CParameterization::ParameterBorderVertex()
{

//计算边界总长度
if (m_BorderVertex.size()<3)
{
return;
}
vec v1v2;
vector<float> VertexU(m_BorderVertex.size(),0.0);//边界顶点的弧长参数化坐标
VertexU[0]=0.0;
for(int i=1;i<m_BorderVertex.size();i++)
{
v1v2=m_mesh->vertices[m_BorderVertex[i]]-m_mesh->vertices[m_BorderVertex[i-1]];
m_SumBorderLength+=len(v1v2);
VertexU[i]=m_SumBorderLength;
}
v1v2=m_mesh->vertices[m_BorderVertex[m_BorderVertex.size()-1]]-m_mesh->vertices[m_BorderVertex[0]];
m_SumBorderLength+=len(v1v2);

for (int i=0;i<VertexU.size();i++)
{
VertexU[i]=VertexU[i]/m_SumBorderLength;
}

//分为正方形和圆形参数域

//暂时未写圆形参数化

//正方形参数化
float unionscal=0.25;
for (int i=0;i<VertexU.size();i++)
{
if(unionscal>VertexU[i])m_ParameterResult[m_BorderVertex[i]]=vec2(VertexU[i]*4.0,0.0);
else if(unionscal*2.0>VertexU[i])m_ParameterResult[m_BorderVertex[i]]=vec2(1.0,(VertexU[i]-unionscal)*4.0);
else if(unionscal*3.0>VertexU[i])m_ParameterResult[m_BorderVertex[i]]=vec2((unionscal*3.0-VertexU[i])*4.0,1.0);
else if(unionscal*4.0>VertexU[i])m_ParameterResult[m_BorderVertex[i]]=vec2(0.0,(1.0-VertexU[i])*4.0);
}
}
void CParameterization::ParameterInnerVertex()
{
if (m_BorderVertex.size()<3)
{
return;
}
int n=m_InteriorVertex.size()+m_BorderVertex.size();
FyMatrix CoefficientMatrix(n,n),b(n,2),pt(n,2);//系数矩阵
CoefficientMatrix.EqualEye();//方程组系数矩阵单位化
b.AllZero();//系数矩阵
for (int i=m_InteriorVertex.size();i<n;i++)
for (int j=0;j<2;j++)
b(i,j)=m_ParameterResult[m_BorderVertex[i-m_InteriorVertex.size()]][j];

float Coefficient;//参数化类型的系数 均匀参数化为度的倒数
int ineighbors;
for (int i=0;i<m_InteriorVertex.size();i++)
{

int  NeighborNumber=m_mesh->neighbors[m_InteriorVertex[i]].size();
vector<float>NeighborCoefficient;
if(m_ParameterType==0) NeighborCoefficient=ShapePreserveParameterization(m_InteriorVertex[i]);
else if((m_ParameterType==2)||(m_ParameterType==3)) NeighborCoefficient=MeanValueHarmonyParameterization(m_InteriorVertex[i]);

for (int j=0;j<NeighborNumber;j++)
{

if (m_ParameterType==1)
Coefficient=1.0/NeighborNumber;//均匀参数
else if (m_ParameterType==0)//
{
Coefficient=NeighborCoefficient[j];
}//保形 else if

else if ((m_ParameterType==2)||(m_ParameterType==3))
{
Coefficient=NeighborCoefficient[j];
}//中值 else if

ineighbors=m_mesh->neighbors[m_InteriorVertex[i]][j];
int Index=0;
vector<int>::iterator iter;
iter=find(m_InteriorVertex.begin(),m_InteriorVertex.end(),ineighbors);
if(iter!=m_InteriorVertex.end())  Index=distance(m_InteriorVertex.begin(),iter);
else
{
iter=find(m_BorderVertex.begin(),m_BorderVertex.end(),ineighbors);
if (iter!=m_BorderVertex.end()) Index=m_InteriorVertex.size()+distance(m_BorderVertex.begin(),iter);
else {AfxMessageBox("参数化失败1");/*assert(0)*/;/* return;*/}

}
CoefficientMatrix(i,Index)=-Coefficient;

}//for

}//for

CoefficientMatrix=CoefficientMatrix.Inverse();
pt=CoefficientMatrix*b;
for (int i=0;i<m_InteriorVertex.size();i++)
for (int j=0;j<2;j++)
{
m_ParameterResult[m_InteriorVertex[i]][j]=pt(i,j);
}

for (int i=0;i<m_BorderVertex.size();i++)
if( (m_ParameterResult[m_BorderVertex[i]][0]!=pt(i+m_InteriorVertex.size(),0))||(m_ParameterResult[m_BorderVertex[i]][1]!=pt(i+m_InteriorVertex.size(),1))  )
{
AfxMessageBox("参数化失败2");
}

}
//均匀参数化
float CParameterization::UniformParameterization()
{

return 0.3;
}
//保形参数化  参数p是三维模型的顶点索引
vector<float> CParameterization::ShapePreserveParameterization(int p)
{
int PNeighborNumber=m_mesh->neighbors[p].size();
//P的邻接顶点排序
vector<int> PNeighbor=m_mesh->neighbors[p];
vector<int> SortNeighbor;
SortNeighbor.push_back(PNeighbor[0]);
vector<int>::iterator iter;
iter=find(PNeighbor.begin(),PNeighbor.end(),PNeighbor[0]);
PNeighbor.erase(iter);
int fid;
while(PNeighbor.size()>0)
{

int i;
for (i=0;i<m_mesh->neighbors[SortNeighbor[SortNeighbor.size()-1]].size();i++)
{
fid=m_mesh->neighbors[SortNeighbor[SortNeighbor.size()-1]][i];
iter=find(PNeighbor.begin(),PNeighbor.end(),fid);
if (iter!=PNeighbor.end()) break;
}
if (i==m_mesh->neighbors[SortNeighbor[SortNeighbor.size()-1]].size())  AfxMessageBox("参数化失败3");
SortNeighbor.push_back(fid);
PNeighbor.erase(iter);
}

/*for (int k=0;k<m_mesh->adjacentfaces[p].size();k++)
{
int q=m_mesh->neighbors[p][k];
TriMesh::Face AdjacentFaces;
AdjacentFaces=m_mesh->faces[m_mesh->adjacentfaces[p][k]];

vector<int>v1v2v3;
for (int a=0;a<3;a++)v1v2v3.push_back(AdjacentFaces[a]);

vector<int>::iterator v;
v=find(v1v2v3.begin(),v1v2v3.end(),p);
v1v2v3.erase(v);//剩余的两个顶点组成一条边

for (int a=0;a<NeighborNumber;a++)
{
if ((v1v2v3[0]==m_mesh->neighbors[p][a])||(v1v2v3[1]==m_mesh->neighbors[p][a])) continue;

}
}//for*/
//把p点的邻接顶点展平到二维平面,P点为坐标原点
//计算每个点的角度
vector<float> Angle2D(PNeighborNumber,0.0);//极角
Angle2D[0]=0.0;
float SumAngle=0.0;

for (int i=1;i<SortNeighbor.size();i++)
{
point pa=m_mesh->vertices[SortNeighbor[i-1]]-m_mesh->vertices[p];
point pb=m_mesh->vertices[SortNeighbor[i]]-m_mesh->vertices[p];
float angle=(pa DOT pb)/(len(pa)*len(pb));
SumAngle+=acos(angle);
Angle2D[i]=SumAngle;
}
point pa=m_mesh->vertices[SortNeighbor[PNeighborNumber-1]]-m_mesh->vertices[p];
point pb=m_mesh->vertices[SortNeighbor[0]]-m_mesh->vertices[p];
float angle=(pa DOT pb)/(len(pa)*len(pb));
SumAngle+=acos(angle);
for (int i=0;i<PNeighborNumber;i++)
{
Angle2D[i]=Angle2D[i]*2.0*3.1415926/SumAngle;
}
//计算极坐标的径向长度
vector<float> Radius(PNeighborNumber,0.0);
for(int i=0;i<PNeighborNumber;i++)
{
point pq=m_mesh->vertices[SortNeighbor[i]]-m_mesh->vertices[p];
Radius[i]=len(pq);
}

//计算二维坐标
vector<vec3> P2Dcoordinate(PNeighborNumber);
vec3  pSeed(0.0,0.0,0.0);
for(int i=0;i<PNeighborNumber;i++)
{
P2Dcoordinate[i][0]=10*Radius[i]*cos(Angle2D[i]);
P2Dcoordinate[i][1]=10*Radius[i]*sin(Angle2D[i]);
P2Dcoordinate[i][2]=0.0;
}
//计算个点邻接顶点的权值
int p1,p2,p3;
vec3 pp1,pp2,pp3;

vector<vector<float>>Weight;
Weight.resize(PNeighborNumber);
for (int i=0;i<PNeighborNumber;i++)
{
Weight[i].reserve(PNeighborNumber+2);
}
for (int i=0;i<PNeighborNumber;i++)
{
p3=i;
pp3=P2Dcoordinate[i]-pSeed;
float Inaccuracy=10000.0;
int accuratep1,accuratep2;
float u=0.0;
float v=0.0;
float w=0.0;
for (int j=0;j<PNeighborNumber;j++)
{

if (j<PNeighborNumber-1)
{
p1=j ;
pp1=P2Dcoordinate[j]-pSeed;
p2=j+1;
pp2=P2Dcoordinate[j+1]-pSeed;
}
else
{
p1=j;
pp1=P2Dcoordinate[j]-pSeed;
p2=0;
pp2=P2Dcoordinate[0]-pSeed;
}
if ((p1!=p3)&&(p2!=p3))
{

//判断P点是否在三角形p1p2p3内部
vec3  p1p2=P2Dcoordinate[p2]-P2Dcoordinate[p1];
vec3  p1p3=P2Dcoordinate[p3]-P2Dcoordinate[p1];
float AreaP1PP2=0.5*len(pp1 CROSS pp2);
float AreaP2PP3=0.5*len(pp2 CROSS pp3);
float AreaP3PP1=0.5*len(pp3 CROSS pp1);
float AreaP1P2P3=0.5*len(p1p2 CROSS p1p3);
float AreaSum=AreaP1PP2+AreaP2PP3+AreaP3PP1;

if(Inaccuracy>(AreaSum-AreaP1P2P3))
{
Inaccuracy=AreaSum-AreaP1P2P3;
accuratep1=p1;
accuratep2=p2;
u=AreaP2PP3/AreaSum;
v=AreaP3PP1/AreaSum;
w=AreaP1PP2/AreaSum;
}
}//if
}//for
Weight[accuratep1].push_back(u);
Weight[accuratep2].push_back(v);
Weight[p3].push_back(w);
}//for
//计算各顶点的平均权值,但此权值是根据排序后顶点的顺序排序的
vector<float> AverageWeight0(PNeighborNumber,0.0);
for (int i=0;i<PNeighborNumber;i++)
{
for (int j=0;j<Weight[i].size();j++)
{
AverageWeight0[i]+=Weight[i][j];
}
AverageWeight0[i]=AverageWeight0[i]/PNeighborNumber;
}
//回归未排序的权值  即最终的权值
vector<float> AverageWeight(PNeighborNumber,0.0);
PNeighbor=m_mesh->neighbors[p];
for (int i=0;i<PNeighborNumber;i++)
{
iter=find(PNeighbor.begin(),PNeighbor.end(),SortNeighbor[i]);
if (iter==PNeighbor.end()) AfxMessageBox("参数化失败4");//Index=distance(m_InteriorVertex.begin(),iter);
int dis=distance(PNeighbor.begin(),iter);
AverageWeight[dis]=AverageWeight0[i];
}

return AverageWeight;
}

//中值坐标参数化
vector<float> CParameterization::MeanValueHarmonyParameterization(int p)
{
int NeighborNumber=m_mesh->neighbors[p].size();
vector<float>w(NeighborNumber,0.0);//存储最后的系数
for (int k=0;k<NeighborNumber;k++)
{
int q=m_mesh->neighbors[p][k];
vector<TriMesh::Face> CommonFaces;
for(int a=0;a<m_mesh->adjacentfaces[p].size();a++)
{   //寻找边ij的两个三角面片

vector<int> &JAdjacentfaces=m_mesh->adjacentfaces[q];
vector<int>::iterator CommonFaceIter;
CommonFaceIter=find(JAdjacentfaces.begin(),JAdjacentfaces.end(),m_mesh->adjacentfaces[p][a]);
if (CommonFaceIter!=JAdjacentfaces.end()) CommonFaces.push_back(m_mesh->faces[m_mesh->adjacentfaces[p][a]]);
}
if(CommonFaces.size()!=2) AfxMessageBox("参数化失败5");//每条边必有两个三角面片
//寻找ij边的两个对顶点

for (int a=0;a<2;a++)
{
vector<int> v1v2v3;
for(int b=0;b<3;b++) v1v2v3.push_back(CommonFaces[a][b]);
vector<int>::iterator v;
v=find(v1v2v3.begin(),v1v2v3.end(),p);
v1v2v3.erase(v);
v=find(v1v2v3.begin(),v1v2v3.end(),q);
v1v2v3.erase(v);
if (v1v2v3.size()!=1)
{
AfxMessageBox("参数化失败6");
}
vec pa=m_mesh->vertices[p]-m_mesh->vertices[v1v2v3[0]];
vec pq=m_mesh->vertices[p]-m_mesh->vertices[q];
vec qa=m_mesh->vertices[p]-m_mesh->vertices[v1v2v3[0]];
//中值坐标参数化
if (m_ParameterType==3)
{
double angle=pa DOT pq;
angle=angle/(len(pa)*len(pq));
angle=tan(acos(angle)/2.0)/len(pq);
w[k]+=angle;
}
//调和映射
else if (m_ParameterType==2)
{
w[k]+=(len2(pa)+len2(qa)-len2(pq))*2.0/len(pa CROSS pq);
}

}
}//k for
float sumw=0.0;
for (int k=0;k<w.size();k++)sumw+=w[k];
for (int i=0;i<NeighborNumber;i++)
{
w[i]=w[i]/sumw;
}

return  w;
}

//判断参数化类型，并调用相应的参数化函数
vector<vec2> CParameterization::ComputeParameterResult()
{
FindBorderInnerVertex();
ParameterBorderVertex();
ParameterInnerVertex();
return m_ParameterResult;
}

往事不堪回首，当时菜鸟一枚，现在看到那样的代码，真的是很惭愧啊。

本文地址：/article/7649378.html
作者：hjimce 联系qq：1393852684 更多资源请关注我的博客：http://blog.csdn.net/hjimce
原创文章，版权所有，转载请保留本行信息。

参考文献：

1、Floater M S. Parameterization and smooth approximation of surface triang-ulations.Computer Aided Geometric Design, 1997,14(3): 231250.

2、 M.Floater and K.Hormann "Surface parameterization: A tutorial and surve- y," in Advances in Multiresolution for Georrzetric Modelling, pp. 157-186 2005.