图形处理（十一）Stroke Parameterization

Stroke Parameterization

原文地址：/article/7649387.html

作者：hjimce

转眼已经过去了好几年，最近开始写技术博客，是为了回顾。《Stroke Parameterization》这篇paper是我人生写的第一篇作者没有提供代码的文章，也是初次学会阅读外文文献的开始。三年前菜鸟一只，连如何通过paper写代码都还不懂，然而没想到这篇paper花了两周的时间，竟然被我搞定了，有了信心，从此菜鸟开始学习起飞……

这篇paper的作者也是一大牛，发了好多篇Siggraph的文章，所以自然这篇paper的质量还是挺不错的。参数化算法的好坏一般是通过纹理贴图的方法，进行验证的。

一、相关理论





Stroke Parameterization顾名思义就是沿着曲线进行参数化的意思，在我的另外一篇博文中《离散指数映射Decal》是以一个点为源点，进行参数化，参数化结果为一圆形参数域。然而在网格曲面上，可能有的时候我们并不紧紧是想要圆形Decal，而是希望沿着曲线进行参数化，比如上面文字贴图中，我们的图片是一张长方形图片，这个时候如果用固定边界的参数化方法，或者用离散指数Decal,它们的参数域一般都类似于圆形，用于上面的贴图肯定不行。



或者又如上图，我们给定的一张鱼的图片是矩形的，把鱼图片贴到那个鱼缸上，这个时候，我们就要用到沿着曲线进行参数化的方法了。

算法原理：



在网格曲面上，参数化无非就是要求解网格顶点的(u,v)坐标，如上图所示，已知曲线C(tx)，我们的目标便是要求出曲线附近区域的每个顶点的(u,v)坐标，也就是我们要求出tx，dx，然后就可以得到二维的参数坐标：





dx表示网格顶点x到曲线的最短测地距离。



算法以Dijkstra算法为遍历依据，根据加权平均的方法，通过已Frozen的邻接顶点更新计算未知的X点的相关信息。在参数化的过程中通过遍历的方法，逐个计算顶点的局部坐标系，参数化坐标。

1、顶点x局部标价的更新



顶点x的坐标基底e1更新公式：



其中Nu(x)表示已经被Dijkstra算法遍历，且标记为Fronzen的顶点(进入队列，并且又从队列中删除的点，也就是已经确定最短距离的顶点)，且其为X点的一邻接顶点。然后x点局部坐标系的n轴为顶点的法矢，这个直接通过邻接三角面片的加权平均就可以计算了。

然后，在已知n,e1轴后，我们可以直接用右手法则确定e2轴，也就是直接通过叉乘的方法确定e2:



2、参数化坐标(u,v)更新

顶点x的参数化坐标更新：





其中权重w(qi,x)的计算公式为：





ε是一个非常小的数，以防分母为0，说的简单一点就是以邻接边长的倒数作为权重。

二、算法实现

后面我将结合我写的代码，进行算法实现讲解，因为这个算法是我还是菜鸟的时候写的代码，然后后面也没有经过整理，只是把效果显示出来，得出结果，所以代码很粗糙，将就一下。

Alogrithm：

1、初始化部分：

初始曲线s={pi}上的点pi相关参数初始化：

a、建立pi的局部标价e1为曲线pi点处的切矢，n为顶底法矢，以此根据右手法则计算出e2

vector<CVector3D>m_e1(m_SeedID.size());
point pt;
for(int i=0;i<m_SeedID.size();i++)
{
//计算种子点数组的e1基底：曲线的切向量作为e1
if (i==0)
{
pt=Tmesh->vertices[m_SeedID[i+1]]-Tmesh->vertices[m_SeedID[i]];
m_nodes[m_SeedID[i]].m_e1=CVector3D(pt[0],pt[1],pt[2]);
m_nodes[m_SeedID[i]].m_e1.Normalize();
}
else if (i<m_SeedID.size()-1)
{
pt=Tmesh->vertices[m_SeedID[i+1]]-Tmesh->vertices[m_SeedID[i-1]];
m_nodes[m_SeedID[i]].m_e1=CVector3D(pt[0],pt[1],pt[2]);
m_nodes[m_SeedID[i]].m_e1.Normalize();
}
else
{
pt=Tmesh->vertices[m_SeedID[i]]-Tmesh->vertices[m_SeedID[i-1]];
m_nodes[m_SeedID[i]].m_e1=CVector3D(pt[0],pt[1],pt[2]);
m_nodes[m_SeedID[i]].m_e1.Normalize();
}
}

m_SeedID为源曲线上点按顺序存储的索引。

b、计算pi点的参数坐标为：



其中α(pi)为沿着曲线s，pi点的累积弧长，就是相当于累积弧长参数化。

vec ab;
vector<double>arccoordate(m_SeedID.size(),0.0);

for(int i=0;i<m_SeedID.size()-1;i++)//弧长参数化
{
ab=Tmesh->vertices[m_SeedID[i+1]]-Tmesh->vertices[m_SeedID[i]];
sumlength+=len(ab);
arccoordate[i+1]=sumlength;
}

c、pi点的测地距离设置为0(Dijkstra算法源点集设置)。

for(int i=0;i<m_SeedID.size();i++)
{
m_nodes[m_SeedID[i]].m_UV.dx=arccoordate[i];//（u，v）坐标设置
m_nodes[m_SeedID[i]].m_UV.dy=0.0;
m_nodes[m_SeedID[i]].m_VisitFlag=CDEMNode::Active;//标记为活动，已加入队列，但还未从队列中删除
m_nodes[m_SeedID[i]].distance_from_source()=0.0;//距离设置为0
m_nodes[m_SeedID[i]].m_GeodesicDistance=m_nodes[m_SeedID[i]].m_UV.GetLength();
m_nodes[m_SeedID[i]].m_Normal=m_VertexNormals[m_SeedID[i]];
m_nodes[m_SeedID[i]].m_e2=m_nodes[m_SeedID[i]].m_e1 * m_nodes[m_SeedID[i]].m_Normal;//局部标价初始化
}

2、Dijkstra算法更新邻域点

根据前面所说的计算方法进行更新参数化坐标，及每个顶点的局部标价。这一步主要就是用到公式2和公式3，然后在结合Dijkstra算法就OK了

void CStrokeParameterization::DecalGeodesicVectors(TriMesh *G1Mesh,double r)
{

G1Mesh->need_normals();
G1Mesh->need_neighbors();
VertexNormals.clear();
m_VertexNormals.clear();
for (int i=0;i<G1Mesh->vertices.size();i++)
{
vec nor;
nor=G1Mesh->normals[i];
VertexNormals.push_back(nor);
nor=normalize(nor);
m_VertexNormals.push_back(CVector3D(nor[0],nor[1],nor[2]));
}

//数据结构转化
unsigned vn,fn;
vn=G1Mesh->vertices.size();
fn=G1Mesh->faces.size();
double *pts;
pts = new double[vn*3];
unsigned *fs;
fs = new unsigned[fn*3];

for (int i=0;i<vn;i++)
{
point p=G1Mesh->vertices[i];
int shift=i*3;
pts[shift]=p[0];
pts[shift+1]=p[1];
pts[shift+2]=p[2];
}
for (int i=0;i<fn;i++)
{
int  shift=i*3;
fs[shift]=G1Mesh->faces[i][0];
fs[shift+1]=G1Mesh->faces[i][1];
fs[shift+2]=G1Mesh->faces[i][2];
}
Gmesh.from_TriMeshData(vn,pts,fn,fs);
delete []pts;
delete []fs;
numOfVertices=G1Mesh->vertices.size();
numOfFaces=G1Mesh->faces.size();

m_nodes.resize(numOfVertices);
for(unsigned i=0; i<m_nodes.size(); ++i)
{
m_nodes[i].vertex() = &Gmesh.vertices()[i];
m_nodes[i].clear();
m_nodes[i].m_VertexID=i;

}
std::set<DEMNode_pointer, CDEMNode>  Queue0;
Queue0.clear();
InitializationSeed();
for(int i=0;i<m_CurveNeighbor.size();i++)
{
Queue0.insert(&m_nodes[m_SeedCurve[i]]);
}
std::vector<double> distances_between_nodes;
std::vector<DEMNode_pointer> neighbor_nodes;
while(!Queue0.empty())
{
DEMNode_pointer min_node = *Queue0.begin();
Queue0.erase(Queue0.begin());
assert(min_node->distance_from_source() < GEODESIC_INF);
min_node->m_VisitFlag=CDEMNode::Frozen;
vector<int>::iterator iter = find(m_CurveNeighbor.begin(),m_CurveNeighbor.end(),min_node->m_VertexID);
if (iter==m_CurveNeighbor.end())
{
min_node->m_e1=Average_e1(min_node->m_VertexID);
min_node->m_Normal=m_VertexNormals[min_node->m_VertexID];
min_node->m_e2=min_node->m_e1 * min_node->m_Normal;
min_node->m_UV=Average_GVector(min_node->m_VertexID);
min_node->m_GeodesicDistance=min_node->m_UV.GetLength();
}

neighbor_nodes.clear();
distances_between_nodes.clear();
list_neighbor_from_node(min_node, neighbor_nodes, distances_between_nodes);
for(unsigned i=0; i<neighbor_nodes.size(); ++i)
{
DEMNode_pointer next_node = neighbor_nodes[i];
if(next_node->distance_from_source() > min_node->distance_from_source() +
distances_between_nodes[i])
{
next_node->distance_from_source() = min_node->distance_from_source() +
distances_between_nodes[i];
next_node->previous() = min_node;
}
}

int neighbor_size_of_u=neighbor_size(min_node->vertex());
if((min_node->m_UV.dx<=(sumlength+r))&&(-r<(min_node->m_UV.dx))&&(abs(min_node->m_UV.dy)<=r))
{
vertex_pointer vertex_of_u=min_node->vertex();
face_pointer adjface_of_u;
for (int i=0;i<vertex_of_u->adjacent_faces().size();i++)
{
adjface_of_u=vertex_of_u->adjacent_faces()[i];
int face_id=adjface_of_u->id();
G1Mesh->faces[face_id].beSelect=true;
}

for (int j=0;j<neighbor_size_of_u;j++)
{
int neighbors_of_u;
neighbors_of_u=neighbor_i(min_node->m_VertexID,j);
if (m_nodes[neighbors_of_u].m_VisitFlag==CDEMNode::Inactive)
{
Queue0.insert(&m_nodes[neighbors_of_u]);
m_nodes[neighbors_of_u].m_VisitFlag=CDEMNode::Active;
}
}
}

}//while
DecalNormalization(r);
}

int CStrokeParameterization::neighbor_size(geodesic::vertex_pointer u)
{
int number;
number=u->adjacent_edges().size();
return number;
}
int CStrokeParameterization::neighbor_i(geodesic::vertex_pointer u,int i)
{
int neighbor_id;
edge_pointer adjacent_e_u=u->adjacent_edges()[i];
neighbor_id=adjacent_e_u->opposite_vertex(u)->id();
return neighbor_id;
}
int CStrokeParameterization::neighbor_i(int u,int i)
{
int neighbor_id;
vertex_pointer vertex_u=&Gmesh.vertices()[u];
edge_pointer adjacent_e_u=vertex_u->adjacent_edges()[i];
vertex_pointer adjacent_v_u=adjacent_e_u->opposite_vertex(vertex_u);
neighbor_id=adjacent_v_u->id();
return neighbor_id;
}
vector<int> CStrokeParameterization::Co_neighbor(int u_id,int v_id)
{
vector<int> co_neighbor;
co_neighbor.clear();
vertex_pointer u,v;
u=&Gmesh.vertices()[u_id];
v=&Gmesh.vertices()[v_id];
for (int i=0;i<neighbor_size(u);i++)
{
int u_nei,v_nei;
u_nei=neighbor_i(u_id,i);
for(int j=0;j<neighbor_size(v);j++)
{
v_nei=neighbor_i(v_id,j);
if (u_nei==v_nei)
{
co_neighbor.push_back(v_nei);
}
}

}

return co_neighbor;
}
//归一化函数
void CStrokeParameterization::DecalNormalization(double radius)
{
double dScale= 1.0/ (3.0*radius);
for (int i=0;i<m_nodes.size();i++)
{
m_nodes[i].m_UV.dx=dScale*m_nodes[i].m_UV.dx;
m_nodes[i].m_UV.dy=dScale*m_nodes[i].m_UV.dy;
}
}

void CStrokeParameterization::list_neighbor_from_node(DEMNode_pointer node,
std::vector<DEMNode_pointer>& storage,
std::vector<double>& distances)
{
vertex_pointer v = node->vertex();
assert(storage.size() == distances.size());

for(unsigned i=0; i<v->adjacent_edges().size(); ++i)
{
edge_pointer e = v->adjacent_edges()[i];
vertex_pointer new_v = e->opposite_vertex(v);
DEMNode_pointer DEMNew_node=&m_nodes[node_index(new_v)];
double l=e->length();
if(DEMNew_node->m_VisitFlag!=CDEMNode::Frozen)
{
storage.push_back(DEMNew_node);
distances.push_back(e->length());
}
}
}

unsigned CStrokeParameterization:: node_index(vertex_pointer v)
{
return v->id();
};

//求侧地矢量的公式
CVector2D CStrokeParameterization::Average_GVector(int q)
{

CVector2D AverageGvector;
vector<int> neighbor=Tmesh->neighbors[q];
CVector3D sume1;
double weight;
double  sumweight=0.0;
for(int i=0;i<neighbor.size();i++)
{
if (m_nodes[neighbor[i]].m_VisitFlag==CDEMNode::Frozen)
{
vec pq=Tmesh->vertices[q]-Tmesh->vertices[neighbor[i]];
weight=1.0/len(pq);
CVector3D pq0(pq[0],pq[1],pq[2]);
pq0=pq0*AlignNormal(q,neighbor[i]);
CVector2D uvofq(pq0|m_nodes[neighbor[i]].m_e1,pq0|m_nodes[neighbor[i]].m_e2);
uvofq=uvofq+m_nodes[neighbor[i]].m_UV;
uvofq=uvofq*weight;
AverageGvector=AverageGvector+uvofq;
sumweight+=weight;
}
}
AverageGvector=AverageGvector/sumweight;

return AverageGvector;
}

//求p的法向量往q法向量的旋转矩阵
CMatrix3D CStrokeParameterization::AlignNormal(int p,int q)
{
double angle=0.0;
vec p_nor,q_nor;
CVector3D p_Nor,q_Nor,Cross_qN_pN;
CMatrix3D rotMtrx1;
p_nor=VertexNormals[p];
q_nor=VertexNormals[q];
p_Nor=CVector3D(p_nor[0],p_nor[1],p_nor[2]);
q_Nor=CVector3D(q_nor[0],q_nor[1],q_nor[2]);
p_Nor.Normalize();
q_Nor.Normalize();
angle=p_Nor|q_Nor;
angle=acos(angle);
Cross_qN_pN=p_Nor*q_Nor;
rotMtrx1=rotMtrx1.CreateRotateMatrix(angle,Cross_qN_pN);

return rotMtrx1;
}
//求p的法向量往q法向量的旋转矩阵,参数为两个点的法向量
CMatrix3D CStrokeParameterization::AlignNormaln(CVector3D np,CVector3D nq)
{
double angle=0.0;
CVector3D p_Nor,q_Nor,Cross_qN_pN;
CMatrix3D rotMtrx1;
p_Nor=np;
q_Nor=nq;
p_Nor.Normalize();
q_Nor.Normalize();
angle=p_Nor|q_Nor;
angle=acos(angle);
Cross_qN_pN=p_Nor*q_Nor;
rotMtrx1=rotMtrx1.CreateRotateMatrix(angle,Cross_qN_pN);
return rotMtrx1;
}

Frame Field的更新显示结果：



这种沿着曲线进行参数化的paper较少，还有另外一篇paper：《Texture Brush: An Interactive Surface Texturing Interface》也是沿着曲线参数化，不过效率速度都感觉没有这个爽。不过那篇paper的纹理贴图的效果看起来倒是挺漂亮的：



贴图
这篇paper就讲到这里吧。参数化在Siggraph上面的paper还是很多的，每一年都有Parameterization相关的模块，所以还有很多paper等着我们去学习

***************作者：hjimce 联系qq：1393852684 更多资源请关注我的博客：http://blog.csdn.net/hjimce
原创文章，转载请保留本行信息***************

参考文献：

1、Texture Brush: An Interactive Surface Texturing Interface

2、Stroke Parameterization

3、Interactive Decal Compositing with Discrete Exponential Maps