图形处理（十一）Stroke Parameterization

转眼已经过去了好几年，最近开始写技术博客，是为了回顾。《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

[cpp] view plain copy

1. vector<CVector3D>m_e1(m_SeedID.size());  
2.    point pt;  
3. for(int i=0;i<m_SeedID.size();i++)  
4. {  
5. //计算种子点数组的e1基底：曲线的切向量作为e1  
6.     if (i==0)  
7.     {  
8.         pt=Tmesh->vertices[m_SeedID[i+1]]-Tmesh->vertices[m_SeedID[i]];  
9.         m_nodes[m_SeedID[i]].m_e1=CVector3D(pt[0],pt[1],pt[2]);  
10.         m_nodes[m_SeedID[i]].m_e1.Normalize();  
11.     }  
12.     else if (i<m_SeedID.size()-1)  
13.     {  
14.         pt=Tmesh->vertices[m_SeedID[i+1]]-Tmesh->vertices[m_SeedID[i-1]];  
15.         m_nodes[m_SeedID[i]].m_e1=CVector3D(pt[0],pt[1],pt[2]);  
16.         m_nodes[m_SeedID[i]].m_e1.Normalize();  
17.     }  
18.     else   
19.     {  
20.         pt=Tmesh->vertices[m_SeedID[i]]-Tmesh->vertices[m_SeedID[i-1]];  
21.         m_nodes[m_SeedID[i]].m_e1=CVector3D(pt[0],pt[1],pt[2]);  
22.         m_nodes[m_SeedID[i]].m_e1.Normalize();  
23.     }  
24. }  

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

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

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

[cpp] view plain copy

1. vec ab;  
2. vector<double>arccoordate(m_SeedID.size(),0.0);  
3.   
4. for(int i=0;i<m_SeedID.size()-1;i++)//弧长参数化  
5. {  
6.     ab=Tmesh->vertices[m_SeedID[i+1]]-Tmesh->vertices[m_SeedID[i]];   
7.     sumlength+=len(ab);  
8.     arccoordate[i+1]=sumlength;  
9. }  

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

[cpp] view plain copy

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

2、Dijkstra算法更新邻域点

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

[cpp] view plain copy

1. void CStrokeParameterization::DecalGeodesicVectors(TriMesh *G1Mesh,double r)     
2. {  
3.   
4.     G1Mesh->need_normals();  
5.     G1Mesh->need_neighbors();  
6.     VertexNormals.clear();  
7.     m_VertexNormals.clear();  
8.     for (int i=0;i<G1Mesh->vertices.size();i++)  
9.     {  
10.         vec nor;  
11.         nor=G1Mesh->normals[i];  
12.         VertexNormals.push_back(nor);  
13.         nor=normalize(nor);  
14.         m_VertexNormals.push_back(CVector3D(nor[0],nor[1],nor[2]));  
15.     }  
16.   
17.     //数据结构转化  
18.     unsigned vn,fn;  
19.     vn=G1Mesh->vertices.size();   
20.     fn=G1Mesh->faces.size();  
21.     double *pts;  
22.     pts = new double[vn*3];  
23.     unsigned *fs;  
24.     fs = new unsigned[fn*3];  
25.   
26.     for (int i=0;i<vn;i++)  
27.     {  
28.         point p=G1Mesh->vertices[i];  
29.         int shift=i*3;  
30.         pts[shift]=p[0];  
31.         pts[shift+1]=p[1];  
32.         pts[shift+2]=p[2];  
33.     }  
34.     for (int i=0;i<fn;i++)  
35.     {  
36.         int  shift=i*3;  
37.         fs[shift]=G1Mesh->faces[i][0];  
38.         fs[shift+1]=G1Mesh->faces[i][1];  
39.         fs[shift+2]=G1Mesh->faces[i][2];  
40.     }  
41.     Gmesh.from_TriMeshData(vn,pts,fn,fs);  
42.     delete []pts;  
43.     delete []fs;  
44.     numOfVertices=G1Mesh->vertices.size();  
45.     numOfFaces=G1Mesh->faces.size();  
46.   
47.   
48.   
49.     m_nodes.resize(numOfVertices);  
50.     for(unsigned i=0; i<m_nodes.size(); ++i)  
51.     {  
52.         m_nodes[i].vertex() = &Gmesh.vertices()[i];  
53.         m_nodes[i].clear();  
54.         m_nodes[i].m_VertexID=i;  
55.           
56.     }  
57.     std::set<DEMNode_pointer, CDEMNode>  Queue0;   
58.     Queue0.clear();  
59.     InitializationSeed();  
60.     for(int i=0;i<m_CurveNeighbor.size();i++)  
61.     {  
62.       Queue0.insert(&m_nodes[m_SeedCurve[i]]);  
63.     }  
64.     std::vector<double> distances_between_nodes;  
65.     std::vector<DEMNode_pointer> neighbor_nodes;  
66.     while(!Queue0.empty())  
67.     {  
68.         DEMNode_pointer min_node = *Queue0.begin();  
69.         Queue0.erase(Queue0.begin());  
70.         assert(min_node->distance_from_source() < GEODESIC_INF);  
71.         min_node->m_VisitFlag=CDEMNode::Frozen;  
72.         vector<int>::iterator iter = find(m_CurveNeighbor.begin(),m_CurveNeighbor.end(),min_node->m_VertexID);  
73.         if (iter==m_CurveNeighbor.end())  
74.         {  
75.             min_node->m_e1=Average_e1(min_node->m_VertexID);  
76.             min_node->m_Normal=m_VertexNormals[min_node->m_VertexID];  
77.             min_node->m_e2=min_node->m_e1 * min_node->m_Normal;  
78.             min_node->m_UV=Average_GVector(min_node->m_VertexID);  
79.             min_node->m_GeodesicDistance=min_node->m_UV.GetLength();  
80.         }  
81.       
82.         neighbor_nodes.clear();  
83.         distances_between_nodes.clear();  
84.         list_neighbor_from_node(min_node, neighbor_nodes, distances_between_nodes);  
85.         for(unsigned i=0; i<neighbor_nodes.size(); ++i)        
86.         {  
87.             DEMNode_pointer next_node = neighbor_nodes[i];  
88.             if(next_node->distance_from_source() > min_node->distance_from_source() +  
89.                 distances_between_nodes[i])  
90.             {  
91.                 next_node->distance_from_source() = min_node->distance_from_source() +  
92.                     distances_between_nodes[i];  
93.                 next_node->previous() = min_node;  
94.             }  
95.         }  
96.       
97.         int neighbor_size_of_u=neighbor_size(min_node->vertex());  
98.         if((min_node->m_UV.dx<=(sumlength+r))&&(-r<(min_node->m_UV.dx))&&(abs(min_node->m_UV.dy)<=r))  
99.         {  
100.             vertex_pointer vertex_of_u=min_node->vertex();  
101.             face_pointer adjface_of_u;  
102.             for (int i=0;i<vertex_of_u->adjacent_faces().size();i++)  
103.             {  
104.                 adjface_of_u=vertex_of_u->adjacent_faces()[i];  
105.                 int face_id=adjface_of_u->id();  
106.                 G1Mesh->faces[face_id].beSelect=true;  
107.             }  
108.   
109.   
110.             for (int j=0;j<neighbor_size_of_u;j++)  
111.             {  
112.                 int neighbors_of_u;  
113.                 neighbors_of_u=neighbor_i(min_node->m_VertexID,j);  
114.                 if (m_nodes[neighbors_of_u].m_VisitFlag==CDEMNode::Inactive)  
115.                 {  
116.                     Queue0.insert(&m_nodes[neighbors_of_u]);  
117.                     m_nodes[neighbors_of_u].m_VisitFlag=CDEMNode::Active;             
118.                 }  
119.             }  
120.         }  
121.       
122.     }//while  
123.    DecalNormalization(r);  
124. }  
125.   
126. int CStrokeParameterization::neighbor_size(geodesic::vertex_pointer u)  
127. {  
128.     int number;  
129.     number=u->adjacent_edges().size();  
130.     return number;  
131. }  
132. int CStrokeParameterization::neighbor_i(geodesic::vertex_pointer u,int i)  
133. {  
134.     int neighbor_id;  
135.     edge_pointer adjacent_e_u=u->adjacent_edges()[i];  
136.     neighbor_id=adjacent_e_u->opposite_vertex(u)->id();  
137.     return neighbor_id;  
138. }  
139. int CStrokeParameterization::neighbor_i(int u,int i)  
140. {  
141.     int neighbor_id;  
142.     vertex_pointer vertex_u=&Gmesh.vertices()[u];  
143.     edge_pointer adjacent_e_u=vertex_u->adjacent_edges()[i];  
144.     vertex_pointer adjacent_v_u=adjacent_e_u->opposite_vertex(vertex_u);  
145.     neighbor_id=adjacent_v_u->id();  
146.     return neighbor_id;  
147. }  
148. vector<int> CStrokeParameterization::Co_neighbor(int u_id,int v_id)  
149. {  
150.     vector<int> co_neighbor;  
151.     co_neighbor.clear();  
152.     vertex_pointer u,v;  
153.     u=&Gmesh.vertices()[u_id];  
154.     v=&Gmesh.vertices()[v_id];    
155.     for (int i=0;i<neighbor_size(u);i++)  
156.     {  
157.         int u_nei,v_nei;  
158.         u_nei=neighbor_i(u_id,i);  
159.         for(int j=0;j<neighbor_size(v);j++)  
160.         {  
161.             v_nei=neighbor_i(v_id,j);  
162.             if (u_nei==v_nei)  
163.             {  
164.                 co_neighbor.push_back(v_nei);  
165.             }  
166.         }  
167.   
168.     }  
169.   
170.     return co_neighbor;  
171. }  
172. //归一化函数  
173. void CStrokeParameterization::DecalNormalization(double radius)  
174. {  
175.     double dScale= 1.0/ (3.0*radius);                   
176.     for (int i=0;i<m_nodes.size();i++)  
177.     {  
178.         m_nodes[i].m_UV.dx=dScale*m_nodes[i].m_UV.dx;                             
179.         m_nodes[i].m_UV.dy=dScale*m_nodes[i].m_UV.dy;   
180.     }  
181. }  
182.   
183.   
184.   
185. void CStrokeParameterization::list_neighbor_from_node(DEMNode_pointer node,   
186.                                                                     std::vector<DEMNode_pointer>& storage,  
187.                                                                     std::vector<double>& distances)  
188. {  
189.     vertex_pointer v = node->vertex();  
190.     assert(storage.size() == distances.size());  
191.   
192.     for(unsigned i=0; i<v->adjacent_edges().size(); ++i)  
193.     {  
194.         edge_pointer e = v->adjacent_edges()[i];  
195.         vertex_pointer new_v = e->opposite_vertex(v);  
196.         DEMNode_pointer DEMNew_node=&m_nodes[node_index(new_v)];  
197.         double l=e->length();  
198.        if(DEMNew_node->m_VisitFlag!=CDEMNode::Frozen)  
199.        {  
200.            storage.push_back(DEMNew_node);  
201.            distances.push_back(e->length());  
202.        }          
203.     }  
204. }  
205.   
206.   
207. unsigned CStrokeParameterization:: node_index(vertex_pointer v)       
208. {  
209.     return v->id();  
210. };  

[cpp] view plain copy

1. //求侧地矢量的公式  
2. CVector2D CStrokeParameterization::Average_GVector(int q)  
3. {     
4.   
5.     CVector2D AverageGvector;  
6.     vector<int> neighbor=Tmesh->neighbors[q];  
7.     CVector3D sume1;  
8.     double weight;  
9.     double  sumweight=0.0;  
10.     for(int i=0;i<neighbor.size();i++)  
11.     {     
12.         if (m_nodes[neighbor[i]].m_VisitFlag==CDEMNode::Frozen)  
13.        {   
14.         vec pq=Tmesh->vertices[q]-Tmesh->vertices[neighbor[i]];  
15.         weight=1.0/len(pq);  
16.         CVector3D pq0(pq[0],pq[1],pq[2]);  
17.         pq0=pq0*AlignNormal(q,neighbor[i]);  
18.         CVector2D uvofq(pq0|m_nodes[neighbor[i]].m_e1,pq0|m_nodes[neighbor[i]].m_e2);  
19.         uvofq=uvofq+m_nodes[neighbor[i]].m_UV;  
20.         uvofq=uvofq*weight;  
21.         AverageGvector=AverageGvector+uvofq;  
22.         sumweight+=weight;  
23.        }  
24.     }  
25.     AverageGvector=AverageGvector/sumweight;  
26.    
27.   
28.   return AverageGvector;  
29. }  

[cpp] view plain copy

1. //求p的法向量往q法向量的旋转矩阵  
2. CMatrix3D CStrokeParameterization::AlignNormal(int p,int q)  
3. {  
4.     double angle=0.0;  
5.     vec p_nor,q_nor;  
6.     CVector3D p_Nor,q_Nor,Cross_qN_pN;  
7.     CMatrix3D rotMtrx1;  
8.     p_nor=VertexNormals[p];  
9.     q_nor=VertexNormals[q];  
10.     p_Nor=CVector3D(p_nor[0],p_nor[1],p_nor[2]);  
11.     q_Nor=CVector3D(q_nor[0],q_nor[1],q_nor[2]);  
12.     p_Nor.Normalize();  
13.     q_Nor.Normalize();  
14.     angle=p_Nor|q_Nor;  
15.     angle=acos(angle);                                      
16.     Cross_qN_pN=p_Nor*q_Nor;  
17.     rotMtrx1=rotMtrx1.CreateRotateMatrix(angle,Cross_qN_pN);  
18.   
19. return rotMtrx1;  
20. }  
21. //求p的法向量往q法向量的旋转矩阵,参数为两个点的法向量  
22. CMatrix3D CStrokeParameterization::AlignNormaln(CVector3D np,CVector3D nq)  
23. {  
24.     double angle=0.0;  
25.     CVector3D p_Nor,q_Nor,Cross_qN_pN;  
26.     CMatrix3D rotMtrx1;  
27.     p_Nor=np;  
28.     q_Nor=nq;  
29.     p_Nor.Normalize();  
30.     q_Nor.Normalize();  
31.     angle=p_Nor|q_Nor;  
32.     angle=acos(angle);                                      
33.     Cross_qN_pN=p_Nor*q_Nor;  
34.     rotMtrx1=rotMtrx1.CreateRotateMatrix(angle,Cross_qN_pN);  
35.     return rotMtrx1;  
36. }  

Frame Field的更新显示结果：

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

贴图

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

本文地址：http://blog.csdn.net/hjimce/article/details/46489913     作者：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