图形处理（十一）Stroke Parameterization

Stroke Parameterization

原文地址：http://blog.csdn.net/hjimce/article/details/46489913

作者：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