mesh.cpp

此代码主要用于汽车车身的优化

/*		
		method using : Beta/k, Gama/k 
		if( k > 4 && k != 0 )
		{
			for( int j = 0; j < k; j ++ )
				q += ((*m_pSubdivVertexes)[ getFace( (*m_pVertexes)[i].m_adjFaceIdxSet[j]).m_newFacePointIdx ]).m_v \
						* getFace( (*m_pVertexes)[i].m_adjFaceIdxSet[j] ).m_vertexIdxSet.size();
			q = q / k / k / 4.0;		// else q just = 0,0,0
			// R:  midpoints
			//   what if edge number != adj face num?? use edge num as k??
			
			//k = (*m_pVertexes)[i].m_adjEdgePosSet.size();
			for( j = 0; j < k; j ++ )
				r += ((*m_pVertexes)[i].m_adjEdgePosSet[j])->m_midPoint;
			r = r / k / k;
			ASSERT(		((*m_pVertexes)[i].m_adjFaceIdxSet.size() == 0)  \
					==	((*m_pVertexes)[i].m_adjEdgePosSet.size() == 0) );
			// S:  original vertex point : (*m_pVertexes)[i].m_v

			int n = m_pVertexes->size();
			// Q + 2R + S / 4
			Vertex vx( (q + r + (*m_pVertexes)[i].m_v * (1-3./k) ) );	// 9 = Beta + Gama + Gama
			m_pSubdivVertexes->push_back( vx );
			(*m_pVertexes)[i].m_newVertexPointIdx = m_pSubdivVertexes->size() - 1;
		}
		else	*/
/*
		if ( k != 0 )		// ordinary point
		{
			// Q:  average of the new face points of all faces adjacent to the original vertex point
			for( int j = 0; j < (*m_pVertexes)[i].m_adjFaceIdxSet.size(); j ++ )
				q += ((*m_pSubdivVertexes)[ getFace( (*m_pVertexes)[i].m_adjFaceIdxSet[j] ).m_newFacePointIdx ]).m_v;
			q /= (*m_pVertexes)[i].m_adjFaceIdxSet.size();		// else q just = 0,0,0
			
			// R:  average of the midpoints of all original edges incident on the original vertex point
			int edgen = (*m_pVertexes)[i].m_adjEdgePosSet.size();
			for( j = 0; j < edgen; j ++ )
				r += ((*m_pVertexes)[i].m_adjEdgePosSet[j])->m_midPoint;
			if ( edgen != 0 )
				r /= edgen;
			else
			{
				cout << "No Use Vertex Exist:";
				(*m_pVertexes)[i].dump();
				continue;
			}
			
			int n = m_pVertexes->size();
			// Q + 2R + S / 4
			Vector3d v = (q + r * 2 + (*m_pVertexes)[i].m_v * (n-3) / n )  / 4;
			(*m_pVertexes)[i].m_newVertexPointIdx = addSubdivVertex(v);

		}
		else
		{	// 可能存在面内，棱上，或者独立的点。即没有被用作顶点的点
			cout << "No Use Vertex Exist:";
			(*m_pVertexes)[i].dump();
			continue;
		}
		
	}
}

*/



void Mesh::prepareCatmullSubdivData()//catmull clark细分预处理
{
	if(m_pSubdivVertexes)
	{
		m_pSubdivVertexes->clear();
		delete m_pSubdivVertexes;
	}
	m_pSubdivVertexes = new VERTEXSET;

	// face central point:
	FOR_EACH_FACE(fitr)//新面点
	{
		Vector3d v;
		for( int j = 0; j < fitr->m_vertexIdxSet.size(); j ++ )
			v += getVertex(fitr->m_vertexIdxSet[j]).m_v;

		fitr->m_newFacePointIdx = addSubdivVertex(v / fitr->m_vertexIdxSet.size());
	}

	// new edgepoint: two adj face point + 2 * midpoint / 4
//	Vector3d v2;
	FOR_EACH_EDGE(vp)//新边点
	{
		ASSERT(vp->m_adjFaceIdxSet.size() <= 2);		// deals with manifold only
		ASSERT(vp->m_adjFaceIdxSet.size() > 0);			// all edges are built from faces
		
		// view otherside as a very slim face, that is, a line covering the edge
		Vector3d v;
		if( vp->m_adjFaceIdxSet.size() == 2 && ( getVertex(vp->m_bgVertexIdx).isOnBoundary() == getVertex(vp->m_edVertexIdx).isOnBoundary()) )	
		{
			// when quadrangle: 1/16 1/16 3/8 3/8 <==> 
			//                  == sum of edge points / 2 + sum of face point / 4
			//                     this is applicable to other polygons too
			v =	(getVertex(vp->m_bgVertexIdx).m_v + getVertex(vp->m_edVertexIdx).m_v) / 4 + \
				 ( ( (*m_pSubdivVertexes)[ getFace( vp->m_adjFaceIdxSet[1] ).m_newFacePointIdx] ).m_v +
				   ( (*m_pSubdivVertexes)[ getFace( vp->m_adjFaceIdxSet[0] ).m_newFacePointIdx] ).m_v )  / 4;
			vp->m_newEdgePointIdx = addSubdivVertex(v);
		}
		else if ( vp->m_adjFaceIdxSet.size() == 2 && getVertex(vp->m_bgVertexIdx).isOnBoundary() )
		{	// 1/16 1/16 5/8 1/8
			v =	getVertex(vp->m_bgVertexIdx).m_v / 2 + \
				 ( ( (*m_pSubdivVertexes)[ getFace( vp->m_adjFaceIdxSet[1] ).m_newFacePointIdx] ).m_v +
				   ( (*m_pSubdivVertexes)[ getFace( vp->m_adjFaceIdxSet[0] ).m_newFacePointIdx] ).m_v )  / 4;
			vp->m_newEdgePointIdx = addSubdivVertex(v);
			
		}
		else if ( vp->m_adjFaceIdxSet.size() == 2 && getVertex(vp->m_edVertexIdx).isOnBoundary() )
		{
			v =	getVertex(vp->m_edVertexIdx).m_v / 2 + \
				 ( ( (*m_pSubdivVertexes)[ getFace( vp->m_adjFaceIdxSet[1] ).m_newFacePointIdx] ).m_v +
				   ( (*m_pSubdivVertexes)[ getFace( vp->m_adjFaceIdxSet[0] ).m_newFacePointIdx] ).m_v )  / 4;
			vp->m_newEdgePointIdx = addSubdivVertex(v);

		}
		else	// vp->m_adjFaceIdxSet.size() == 1
		{
			vp->m_newEdgePointIdx = addSubdivVertex((getVertex(vp->m_bgVertexIdx).m_v + getVertex(vp->m_edVertexIdx).m_v) / 2);
		}
	}
	// new vertex point:  Beta/k * all edge + Gama/k * all anti vertex + (1-Beta-Gama)old vertex 
	//计算新点的位置
	FOR_EACH_VERTEX(vitr)//新顶点
	{
		Vector3d v;							// initially zero
		int k = vitr->m_adjFaceIdxSet.size();
		if ( ! vitr->isOnBoundary() )		// ordinary point: k must equal to number of adj edges
		{
			float beta = 3. / 2. / k;
			float gama = 1. / 4. / k;
			Vector3d v;
			
			// add all point on other side of the edge ( * beta / k )
			for( int i = 0; i < k; i ++ )
			{
				v += (getVertex( getEdge(vitr->m_adjEdgePosSet[i]).getCounterVertex(vitr->m_index) 
								).m_v * beta / k);
			}
			// add all face point not directly connected to this point but on the same face with this point ( * game / k )
			for( i = 0; i < k; i ++ )
			{
				Face &f = getFace(vitr->m_adjFaceIdxSet[i]);
				// if quadrangle: += the 4th point on that face * gama / k
				// else           += all other points(not directly connected to vitr) on that face / (pointNum-3) * gama / k
				//   what if triangles?????
				int fk = f.m_vertexIdxSet.size();
				for( int vi = 0; vi < fk;vi ++ )
					if( f.m_vertexIdxSet[vi] == vitr->m_index )
					{
						vi += 2;	// skip to the first other points' index ( whatever the dir is, clockwise or reversely)
						break;
					}
				for( int j = 0; j < fk - 3; j ++ )
				{
					v += ( getVertex(f.m_vertexIdxSet[(vi+j) % fk]).m_v * gama / k / (fk - 3) );
				}
			}
			v += ( vitr->m_v * (1. - beta - gama) );
			vitr->m_newVertexPointIdx = addSubdivVertex(v);

		}
		else if ( k > 0 && vitr->m_adjEdgePosSet.size() >= 3 )//( k == 1 && (*m_pVertexes)[i].m_adjEdgePosSet.size() == 2 )
		{
			for (int i = 0; i < vitr->m_adjEdgePosSet.size(); i ++ )
			{
				if( getEdge(vitr->m_adjEdgePosSet[i]).isOnBoundary() )
					v += getVertex(vitr->m_adjEdgePosSet[i]->getCounterVertex(vitr->m_index)).m_v;
			}
			v = (vitr->m_v * 6 + v) / 8;
			vitr->m_newVertexPointIdx = addSubdivVertex(v);
		}
		else if (vitr->m_adjEdgePosSet.size() == 2 && vitr->m_adjFaceIdxSet.size())	// keep the vertex on corner unchanged
		{
			vitr->m_newVertexPointIdx = addSubdivVertex(vitr->m_v);
		}
		else
		{	// 可能存在面内，棱上，或者独立的点。即没有被用作顶点的点
			cout << "No Use Vertex Exist: (in prepareCatmullSubdivData )";
			vitr->dump();
			continue;
		}
		
	}
}
















// Each new face point is connected to the new edge points of the old face's edges
// Each new vertex point is connected to the new edge points of all original adj-edges.

//  we can find that each face is divided into n faces, where n is the number of its edges or vertex
Mesh* Mesh::calcCatmallSubdivMesh()//catmull clark细分开始
{
	prepareCatmullSubdivData();

	bool blFound;
	Mesh *pSubMesh = new Mesh();
	// link(copy) vertexes data
	pSubMesh->m_pVertexes = this->m_pSubdivVertexes;

	for( int i = 0; i < m_faces.size(); i ++)			// foreach face: generate vertex_num new faces
	{
		Vector3d v;
		for( int j = 0; j < m_faces[i].m_vertexIdxSet.size(); j ++ )// foreach vertex: add 1 face
		{
			vector<int> newVertexIdxSet;
			int vi = ( getVertex(m_faces[i].m_vertexIdxSet[j]) ).m_newVertexPointIdx;	// new vi
			newVertexIdxSet.push_back(vi);				// first vertex;

			vi = m_faces[i].m_vertexIdxSet[j];			// old vertex index
			ASSERT( (m_faces[i].m_edgePosSet[j])->m_bgVertexIdx == vi || (m_faces[i].m_edgePosSet[j])->m_edVertexIdx == vi );
			
			int k = j;
			// face => edges => edge whose bgVertex = (*m_pSubdivVertexes)[j]

			newVertexIdxSet.push_back( (m_faces[i].m_edgePosSet[k])->m_newEdgePointIdx);// next edge point in clock wise direction
			newVertexIdxSet.push_back( m_faces[i].m_newFacePointIdx );	// next point: face point

			// get the anti-clock-wise edge.
			k = (k + m_faces[i].m_edgePosSet.size() - 1) % m_faces[i].m_edgePosSet.size();
			ASSERT(m_faces[i].m_edgePosSet[k]->m_edVertexIdx == vi || m_faces[i].m_edgePosSet[k]->m_bgVertexIdx == vi);

			newVertexIdxSet.push_back( (m_faces[i].m_edgePosSet[k])->m_newEdgePointIdx);// next point on edge in clock wise direction

			pSubMesh->addFace(newVertexIdxSet);
		}
	}
	// break old link
	this->m_pSubdivVertexes = NULL;
	pSubMesh->printMeshInfo();
	return pSubMesh;
}

void Mesh::prepareDooSabinSubdivMesh()//DooSabin细分前处理
{
	// 对每个面的每个顶点生成一个相应的新顶点，加入pNewVertexs, 保存index

	// For each vertex Pi of each face of the object, generate a new point P0i as average of the vertex, the
	// two edge points and the face point of the face.
	if(m_pSubdivVertexes)
	{
		m_pSubdivVertexes->clear();
		delete m_pSubdivVertexes;
	}
	m_pSubdivVertexes = new VERTEXSET;

/*	FOR_EACH_EDGE(vp)
	{
		vp->setMidPoint( (getVertex(vp->m_bgVertexIdx).getVertex3d() + 
							getVertex(vp->m_edVertexIdx).getVertex3d()) / 2 );
	}
*/
	Mesh::EDGEPOS outEdgePos, inEdgePos;//定义内部边和外部边
	FOR_EACH_FACE(fitr)
	{
		Vector3d v, fPoint;			// initialize to be zero
		// face central point:
		/*
		for( int j = 0; j < fitr->m_vertexIdxSet.size(); j ++ )
			fPoint += getVertex(fitr->m_vertexIdxSet[j]).getVertex3d();
		fPoint /= fitr->m_vertexIdxSet.size();
		*/
		double mask[10];
		mask[0] = 0.25 + 1.25/fitr->m_vertexIdxSet.size();//(n+5)/4/n,n是面的边数
		for( int j = 1; j < fitr->m_vertexIdxSet.size(); j++ )//求K值
			mask[j] = (3+2*cos(2*j*3.1415926/fitr->m_vertexIdxSet.size()))/4/fitr->m_vertexIdxSet.size();

		int k = fitr->m_vertexIdxSet.size();
		// for each vertex, generate a new VertexPoint in that face
		for( j = 0; j < k; j ++ )
		{
			Vector3d newv;
			
			if( ! getVertex(fitr->m_vertexIdxSet[j]).isOnBoundary() )
			{
				// calc new vertex, push into vertex vector and save the dooSabin vertex index vector
				for( int k = 0; k < fitr->m_vertexIdxSet.size(); k ++ )
					newv += getVertex(fitr->m_vertexIdxSet[(k+j) % fitr->m_vertexIdxSet.size()]).m_v * mask[k];
			}
			else if( getVertex(fitr->m_vertexIdxSet[j]).m_adjEdgePosSet.size() > 2 ) // has two adj face
			{
				// find the boundary edge on face:
				if( getEdge(fitr->m_edgePosSet[j]).isOnBoundary() )
					newv = getVertex(fitr->m_edgePosSet[j]->getCounterVertex(fitr->m_vertexIdxSet[j])).m_v;
				else 
					newv = getVertex(fitr->m_edgePosSet[(j+k-1) % k]->getCounterVertex(fitr->m_vertexIdxSet[j])).m_v;

				newv = (getVertex(fitr->m_vertexIdxSet[j]).m_v * 3 + newv) / 4;
			}
			else	// has one adj face & two adjEdge只有两条边
			{
				ASSERT(    getVertex(fitr->m_vertexIdxSet[j]).m_adjEdgePosSet.size() == 2 \
						&& getVertex(fitr->m_vertexIdxSet[j]).m_adjFaceIdxSet.size() == 1 );
				newv = getVertex(fitr->m_vertexIdxSet[j]).m_v;
				
			}
/*			// get two adj edge idx
			outEdgePos = fitr->m_edgePosSet[j];
			inEdgePos = fitr->m_edgePosSet[ (j+fitr->m_vertexIdxSet.size()-1) % fitr->m_vertexIdxSet.size()];
			Vertex v( (	fPoint * (fitr->m_vertexIdxSet.size() + 3) / fitr->m_vertexIdxSet.size() + getVertex(fitr->m_vertexIdxSet[j]).m_v + \
						outEdgePos->m_midPoint + inEdgePos->m_midPoint)	 / 4, m_pSubdivVertexes->size()); */
			fitr->m_dooSabinVertexIdxSet.push_back( addSubdivVertex(newv) );	
		}
	}
	FOR_EACH_VERTEX(vitr)
	{
		if(vitr->isOnBoundary())
			continue;
		vitr->dump();
		sortVertexAdjFace(vitr->m_index);
	}
}

Mesh* Mesh::calcDooSabinSubdivMesh()
{
	prepareDooSabinSubdivMesh();
	bool blFound;
	Mesh *pSubMesh = new Mesh();
	// link(copy) vertexes data连接顶点数据
	pSubMesh->m_pVertexes = this->m_pSubdivVertexes;

	// 对每个面：连接新顶点，生成一个新面


	// 对每个顶点，连接其在相邻面的各个新顶点，生成一个面
	FOR_EACH_VERTEX(vitr)
	{
		int adjFaceNum = vitr->m_adjFaceIdxSet.size();
		if( adjFaceNum > 0 )
		{
			if( vitr->isOnBoundary() )	// keep the vertex on corner uncut
				continue;
			
			// 逆时针顺序：
			vector<int> vFace;
			for(int i = adjFaceNum - 1; i >= 0; i --)
			{	// 获得邻面
				blFound = false;
				Face &f = getFace(vitr->m_adjFaceIdxSet[i]);
	/*
				cout << vitr->m_adjFaceIdxSet[i];
				f.dump();
				vitr->dump();
	*/		
				// 获得vitr顶点在邻面的上对应的新顶点
				for( int j = 0; j < f.m_vertexIdxSet.size(); j ++)
					if( f.m_vertexIdxSet[j] == vitr->m_index)		// j is the index of this vertex in face m_adjFaceIdxSet[i]
					{
						vFace.push_back(f.m_dooSabinVertexIdxSet[j]);// corresponding dooSabinVertex[j] is new vertex idx
						blFound = true;
						break;
					}
				ASSERT(blFound);
			}

			pSubMesh->addFace(vFace);	// 由其在相邻面的各个新定点生成一个面
			vFace.clear();
			//if (vitr->m_adjFaceIdxSet.size() == 1 && vitr->m_adjEdgeIdxSet.size() == 2)	// this vertex will be unchangely copied when preparing

		}
		else
		{
			cout << "Unused Vertex: ";
			vitr->dump();
		}
	}
	// 对每条棱，连接其定点在其相邻面的分别两个新顶点，共4点，生成一个面
	FOR_EACH_EDGE(eitr)
	{
		if (eitr->isOnBoundary())	// keep the edge on boundary uncut（未删减）
			continue;

		vector<int> eFace;
		// 找到棱在邻面的idx
		blFound = false;
		Face &f1 = getFace(eitr->m_adjFaceIdxSet[0]);
		for (int i = 0; i < f1.m_edgePosSet.size(); ++ i)
		{
			if(f1.m_edgePosSet[i] == eitr)	// from face's vertex i to i+1 is this edge
			{	//  reversely add face's vertex to new edge-face's vertex 
				eFace.push_back(f1.m_dooSabinVertexIdxSet[(i+1) % f1.m_edgePosSet.size()]);	// add i+1 vertex
				eFace.push_back(f1.m_dooSabinVertexIdxSet[i]);		// add i vertex
				blFound = true;
				break;
			}
		}
		ASSERT(blFound);
		blFound = false;
		Face &f2 = getFace(eitr->m_adjFaceIdxSet[1]);
		for ( i = 0; i < f2.m_edgePosSet.size(); ++ i)
		{
			if(f2.m_edgePosSet[i] == eitr)	// from face's vertex i to i+1 is this edge
			{	//  reversely add face's vertex to new edge-face's vertex 
				eFace.push_back(f2.m_dooSabinVertexIdxSet[(i+1) % f2.m_edgePosSet.size()]);	// add i+1 vertex
				eFace.push_back(f2.m_dooSabinVertexIdxSet[i]);		// add i vertex
				blFound = true;
				break;
			}
		}
		ASSERT(blFound);

		pSubMesh->addFace(eFace);	
		eFace.clear();
	}
	FOR_EACH_FACE(fitr)
	{
		pSubMesh->addFace(fitr->m_dooSabinVertexIdxSet);
	}
	this->m_pSubdivVertexes = NULL;
	pSubMesh->printMeshInfo();
	return pSubMesh;
}
//以下是根3细分函数
void Mesh::preparesqrt3SubdivMesh()//以下是对新顶点位置的计算
{//首先求aerfa的值，同loop细分求法相同
	if(m_pSubdivVertexes)
	{
		m_pSubdivVertexes->clear();
		delete m_pSubdivVertexes;
	}
	m_pSubdivVertexes = new VERTEXSET;
	
	// for each vertex, calc a new Vertex对于每一个顶点产生一个新顶点
	double aerfa[20];
	for( int i = 1; i < 20; i ++ )
	{
		aerfa[i] = (4.0-2*cos(2 * 3.1415926 / i))/9.0;
	//aerfa[i] =(3-2*cos(2*3.1415/i))/9/i;
	}
	//proposed by Warren
   //求新面点
	FOR_EACH_FACE(fitr)
	{   //Mesh::EDGEPOS vp;
		Vector3d v;
		for( int j = 0; j < fitr->m_vertexIdxSet.size(); j ++ )
			v += getVertex(fitr->m_vertexIdxSet[j]).m_v;

		fitr->m_newFacePointIdx = addSubdivVertex(v / fitr->m_vertexIdxSet.size());
	
	}
	//求新顶点
	FOR_EACH_VERTEX(vitr)