计算机视觉
图像处理

图形处理(十一)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

  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点的累积弧长,就是相当于累积弧长参数化。

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

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

  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了

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

 

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

 

  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. return rotMtrx1;
  19. }
  20. //求p的法向量往q法向量的旋转矩阵,参数为两个点的法向量
  21. CMatrix3D CStrokeParameterization::AlignNormaln(CVector3D np,CVector3D nq)
  22. {
  23.     double angle=0.0;
  24.     CVector3D p_Nor,q_Nor,Cross_qN_pN;
  25.     CMatrix3D rotMtrx1;
  26.     p_Nor=np;
  27.     q_Nor=nq;
  28.     p_Nor.Normalize();
  29.     q_Nor.Normalize();
  30.     angle=p_Nor|q_Nor;
  31.     angle=acos(angle);
  32.     Cross_qN_pN=p_Nor*q_Nor;
  33.     rotMtrx1=rotMtrx1.CreateRotateMatrix(angle,Cross_qN_pN);
  34.     return rotMtrx1;
  35. }

Frame Field的更新显示结果:

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

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

参考文献:

1、Texture Brush: An Interactive Surface Texturing Interface

2、Stroke Parameterization

3、Interactive Decal Compositing with Discrete Exponential Maps

转载注明来源:CV视觉网 » 图形处理(十一)Stroke Parameterization

分享到:更多 ()
扫描二维码,给作者 打赏
pay_weixinpay_weixin

请选择你看完该文章的感受:

0不错 1超赞 0无聊 0扯淡 0不解 0路过

评论 3

评论前必须登录!