www.pudn.com > 3dSimplifier.zip > Mesh.cpp, change:2000-01-12,size:33363b


// Mesh.cpp: implementation of the CMesh class. 
// 
////////////////////////////////////////////////////////////////////// 
 
#include "stdafx.h" 
#include <gl/gl.h> 
#include "Mesh.h" 
#include "Heap.h" 
 
#ifdef _DEBUG 
#undef THIS_FILE 
static char THIS_FILE[]=__FILE__; 
#define new DEBUG_NEW 
#endif 
 
////////////////////////////////////////////////////////////////////// 
// Construction/Destruction 
////////////////////////////////////////////////////////////////////// 
 
CMesh::CMesh() 
{ 
	m_Triangle = NULL; 
	m_sTriangle = NULL; 
	m_Vertex = NULL; 
	m_sVertex = NULL; 
	m_nTriangle = 0 ; 
	m_nVertex = 0; 
	m_VNormal = NULL; 
	m_TNormal = NULL; 
	m_sVNormal = NULL; 
	m_sTNormal = NULL; 
 
	m_pProgress = NULL; 
} 
 
CMesh::~CMesh() 
{ 
	if (m_Vertex) 
		delete []m_Vertex; 
	if (m_sVertex) 
		delete []m_sVertex; 
	if (m_Triangle) 
		delete []m_Triangle; 
	if (m_sTriangle) 
		delete []m_sTriangle; 
	if (m_VNormal) 
		delete []m_VNormal; 
	if (m_TNormal) 
		delete []m_TNormal; 
	if (m_sVNormal) 
		delete []m_sVNormal; 
	if (m_sTNormal) 
		delete []m_sTNormal; 
} 
 
int CMesh::ReadData(FILE *fp) 
{ 
	int i,j; 
	float box[2][3]; 
	float center[3]; 
 
	fscanf(fp,"%d%d",&m_nVertex,&m_nTriangle); 
 
	if (m_pProgress) 
		m_pProgress->SetRange(0, m_nVertex+m_nTriangle); 
	m_Vertex = new Coord3[m_nVertex]; 
	m_Triangle = new Triangle[m_nTriangle]; 
 
	for (i=0;i<m_nVertex;i++) 
	{ 
		fscanf(fp,"%f%f%f", &(m_Vertex[i][0]), &(m_Vertex[i][1]), &(m_Vertex[i][2])); 
		m_pProgress->SetPos(i); 
	} 
	for (i=0;i<m_nTriangle;i++) 
	{ 
		fscanf(fp,"%d%d%d", &(m_Triangle[i][0]), &(m_Triangle[i][1]), &(m_Triangle[i][2])); 
		m_pProgress->SetPos(i+m_nVertex); 
	} 
 
	box[0][0] = box[0][1] = box[0][2] = MAXFLOAT; 
	box[1][0] = box[1][1] = box[1][2] = -MAXFLOAT; 
	for (i=0;i<m_nVertex;i++) 
		for (j=0;j<3;j++) 
		{ 
			if (box[0][j]>m_Vertex[i][j]) 
				box[0][j] = m_Vertex[i][j]; 
			if (box[1][j]<m_Vertex[i][j]) 
				box[1][j] = m_Vertex[i][j]; 
		} 
	center[0] = (box[0][0]+box[1][0])/2.0; 
	center[1] = (box[0][1]+box[1][1])/2.0; 
	center[2] = (box[0][2]+box[1][2])/2.0; 
	float scale; 
	scale = FMAX(box[1][0]-box[0][0],FMAX(box[1][1]-box[0][1],box[1][2]-box[0][2])); 
	scale /=2.0; 
	for (i=0;i<m_nVertex;i++) 
	{ 
		VEC3_VOPV_OP_S(m_Vertex[i],m_Vertex[i],-,center,/,scale); 
	} 
	VEC3_VOPV_OP_S(box[0],box[0],-,center,/,scale); 
	VEC3_VOPV_OP_S(box[1],box[1],-,center,/,scale); 
	VEC3_ZERO(center); 
 
	m_VNormal = new Coord3[m_nVertex]; 
	m_TNormal = new Coord3[m_nTriangle]; 
	m_nsTriangle = m_nTriangle; 
	m_nsVertex = m_nVertex; 
	m_sTriangle = new Triangle[m_nTriangle]; 
	m_sVertex = new Coord3[m_nVertex]; 
 
	for (i=0;i<m_nVertex;i++) 
	{ 
		VEC3_ASN_OP(m_sVertex[i],=,m_Vertex[i]); 
	} 
	for (i=0;i<m_nTriangle;i++) 
	{ 
		VEC3_ASN_OP(m_sTriangle[i],=,m_Triangle[i]); 
	} 
 
	for (i=0;i<m_nVertex;i++) 
	{ 
		VEC3_ZERO(m_VNormal[i]); 
	} 
	for (i=0;i<m_nTriangle;i++) 
	{ 
		Coord3 vect[3]; 
		VEC3_V_OP_V(vect[0],m_Vertex[m_Triangle[i][1]],-,m_Vertex[m_Triangle[i][0]]); 
		VEC3_V_OP_V(vect[1],m_Vertex[m_Triangle[i][2]],-,m_Vertex[m_Triangle[i][0]]); 
		CROSSPROD3(vect[2],vect[0],vect[1]); 
		V3Standarize(vect[2]); 
		VEC3_ASN_OP(m_TNormal[i],=,vect[2]); 
		for (j=0;j<3;j++) 
		{ 
			VEC3_V_OP_V(m_VNormal[m_Triangle[i][j]],m_VNormal[m_Triangle[i][j]],+,vect[2]); 
		} 
	} 
	for (i=0;i<m_nVertex;i++) 
		V3Standarize(m_VNormal[i]); 
 
	m_sVNormal = new Coord3[m_nVertex]; 
	m_sTNormal = new Coord3[m_nTriangle]; 
	for (i=0;i<m_nVertex;i++) 
	{ 
		VEC3_ASN_OP(m_sVNormal[i],=,m_VNormal[i]); 
	} 
 
	for (i=0;i<m_nTriangle;i++) 
	{ 
		VEC3_ASN_OP(m_sTNormal[i],=,m_TNormal[i]); 
	} 
 
	return m_nTriangle; 
} 
 
void CMesh::SaveData(FILE *fp) 
{ 
	int i; 
 
	fprintf(fp,"%d %d\n",m_nsVertex,m_nsTriangle); 
 
	if (m_pProgress) 
		m_pProgress->SetRange(0, m_nsVertex+m_nsTriangle); 
 
	for (i=0;i<m_nsVertex;i++) 
	{ 
		fprintf(fp,"%f %f %f\n", m_sVertex[i][0], m_sVertex[i][1], m_sVertex[i][2]); 
		m_pProgress->SetPos(i); 
	} 
	for (i=0;i<m_nsTriangle;i++) 
	{ 
		fprintf(fp,"%d %d %d\n", m_sTriangle[i][0], m_sTriangle[i][1], m_sTriangle[i][2]); 
		m_pProgress->SetPos(i+m_nsVertex); 
	} 
} 
 
int CMesh::DrawModel(int RenderMode, int RenderModel) 
{ 
	int i; 
 
	if (RenderMode == 0) 
	{ 
		glDisable(GL_LIGHTING); 
		if (RenderModel == 0) 
		{ 
			for (i=0; i<m_nTriangle; i++) 
			{ 
				glBegin(GL_LINE_LOOP); 
				glVertex3fv(m_Vertex[m_Triangle[i][0]]); 
				glVertex3fv(m_Vertex[m_Triangle[i][1]]); 
				glVertex3fv(m_Vertex[m_Triangle[i][2]]); 
				glEnd(); 
			} 
		} 
		else if (RenderModel == 1) 
		{ 
			for (i=0; i<m_nsTriangle; i++) 
			{ 
				glBegin(GL_LINE_LOOP); 
				glVertex3fv(m_sVertex[m_sTriangle[i][0]]); 
				glVertex3fv(m_sVertex[m_sTriangle[i][1]]); 
				glVertex3fv(m_sVertex[m_sTriangle[i][2]]); 
				glEnd(); 
			} 
		} 
		glEnable(GL_LIGHTING); 
	} 
	else if (RenderMode == 1) 
	{ 
		glBegin(GL_TRIANGLES); 
		if (RenderModel == 0) 
		{ 
			for (i=0; i<m_nTriangle; i++) 
			{ 
				glNormal3fv(m_TNormal[i]); 
				glVertex3fv(m_Vertex[m_Triangle[i][0]]); 
				glVertex3fv(m_Vertex[m_Triangle[i][1]]); 
				glVertex3fv(m_Vertex[m_Triangle[i][2]]); 
			} 
		} 
		else if (RenderModel == 1) 
		{ 
			for (i=0; i<m_nsTriangle; i++) 
			{ 
				glNormal3fv(m_sTNormal[i]); 
				glVertex3fv(m_sVertex[m_sTriangle[i][0]]); 
				glVertex3fv(m_sVertex[m_sTriangle[i][1]]); 
				glVertex3fv(m_sVertex[m_sTriangle[i][2]]); 
			} 
		} 
		glEnd(); 
	} 
 
	glEnable(GL_LIGHTING); 
	if (RenderModel == 0) 
		return m_nTriangle; 
	else if (RenderModel == 1) 
		return m_nsTriangle; 
	else 
		return 0; 
} 
 
int CMesh::Simplify(float ratio) 
{ 
  Matrix4x4 *tri_Q; 
  Matrix4x4 *vex_Q; 
  Coord3 *Vertex; 
  Triangle *Mesh; 
   
  int **vex_tri; 
  int **vRing; 
  int *Vex_map; 
  int i,j,k; 
  int tmp,tmp1,tmp2; 
 
  int **cad_pair; 
  float **pair_cost; 
  Coord3 **pair_pos; 
 
  BOOL *bMark; 
  int  *nMap; 
  int  *mMap; 
 
  if (m_nTriangle == 0) 
	  return 0; 
 
  delete []m_sVertex; 
  delete []m_sTriangle; 
  delete []m_sVNormal; 
  delete []m_sTNormal; 
 
  Vertex= new Coord3[m_nVertex]; 
  for (i=0;i<m_nVertex;i++) 
  { 
	  VEC3_ASN_OP(Vertex[i],=,m_Vertex[i]); 
  } 
  Mesh = new Triangle[m_nTriangle]; 
  for (i=0; i<m_nTriangle; i++) 
  { 
	  VEC3_ASN_OP(Mesh[i],=,m_Triangle[i]); 
  } 
   
  vex_tri=(int **)malloc(m_nVertex*sizeof(int *)); 
  vRing=(int **)malloc(m_nVertex*sizeof(int *)); 
  for (i=0;i<m_nVertex;i++) { 
	vex_tri[i]=(int *)malloc(6*Size_int); 
	vex_tri[i][0]=0; 
	vRing[i]=(int *)malloc(6*Size_int); 
	vRing[i][0]=0; 
  } 
  cad_pair=(int **)malloc(m_nVertex*sizeof(int *)); 
  for (i=0;i<m_nVertex;i++) { 
	cad_pair[i]=(int *)malloc(6*Size_int); 
	cad_pair[i][0]=0; 
  } 
  pair_cost=(float **)malloc(m_nVertex*sizeof(float *)); 
  pair_pos=(Coord3 **)malloc(m_nVertex*sizeof(Coord3 *)); 
  Vex_map= new int[m_nVertex]; 
  for (i=0;i<m_nVertex;i++) 
	Vex_map[i]=i; 
 
  for (k=0;k<m_nTriangle;k++) 
  { 
	for (i=0;i<3;i++) { 
		tmp=tmp1=Mesh[k][i]; 
		tmp2=Mesh[k][(i+2)%3]; 
		for (j=1;j<vRing[tmp][0]+1;j++) 
			if (vRing[tmp][j] == tmp2) 
				break; 
		if (j==vRing[tmp][0]+1) 
		{ 
			vRing[tmp][j]=tmp2; 
			vRing[tmp][0] += 1; 
			if (!(vRing[tmp][0]%5)) 
				vRing[tmp] = (int *)realloc(vRing[tmp],(vRing[tmp][0]+6)*sizeof(int)); 
		} 
		tmp2=Mesh[k][(i+1)%3]; 
		for (j=1;j<vRing[tmp][0]+1;j++) 
			if (vRing[tmp][j] == tmp2) 
				break; 
		if (j==vRing[tmp][0]+1) 
		{ 
			vRing[tmp][j]=tmp2; 
			vRing[tmp][0] += 1; 
			if (!(vRing[tmp][0]%5)) 
				vRing[tmp] = (int *)realloc(vRing[tmp],(vRing[tmp][0]+6)*sizeof(int)); 
		} 
		if (vex_tri[tmp]==NULL) { 
			vex_tri[tmp]=(int *)malloc(6*sizeof(int)); 
			vex_tri[tmp][0]=1; 
			vex_tri[tmp][1]=k; 
		} 
		else { 
			if (!(vex_tri[tmp][0]%5)) 
				vex_tri[tmp]=(int *)realloc(vex_tri[tmp],(vex_tri[tmp][0]+6)*sizeof(int)); 
			vex_tri[tmp][0]=vex_tri[tmp][0]+1; 
			vex_tri[tmp][vex_tri[tmp][0]]=k; 
		} 
		if (tmp1>tmp2) { 
			tmp=tmp1; 
			tmp1=tmp2; 
			tmp2=tmp; 
		} 
		for (j=1;j=cad_pair[tmp1][0];j++) 
			if (cad_pair[tmp1][j]==tmp2) 
				break; 
		if (j>cad_pair[tmp1][0]) { 
			if (!(cad_pair[tmp1][0]%5)) 
				cad_pair[tmp1]=(int *)realloc(cad_pair[tmp1],(cad_pair[tmp1][0]+6)*sizeof(int)); 
			cad_pair[tmp1][0]=j; 
			cad_pair[tmp1][j]=tmp2; 
		} 
	} 
 
  } 
 
  BOOL *border; 
  border = new BOOL[m_nVertex]; 
  for (i=0;i<m_nVertex;i++) 
  { 
	  if (vex_tri[i][0]==vRing[i][0]) 
		  border[i]=FALSE; 
	  else 
		  border[i]=TRUE; 
	  free(vRing[i]); 
  } 
  free(vRing); 
 
  tri_Q=(Matrix4x4 *)malloc(m_nTriangle*Size_Matrix4x4); 
  cal_tri_Q(m_nTriangle,Mesh,Vertex,m_TNormal,tri_Q); 
  vex_Q=(Matrix4x4 *)malloc(m_nVertex*Size_Matrix4x4); 
  cal_vex_Q(m_nVertex,border,vex_tri,Mesh,Vertex,m_TNormal,tri_Q,vex_Q); 
  free(tri_Q); 
 
 
  bMark = new BOOL[m_nVertex]; 
  for (i=0;i<m_nVertex;i++) 
	  bMark[i]=FALSE; 
  int nvex = 0; 
  for (i=0;i<m_nTriangle;i++) 
	  for (j=0;j<3;j++) 
		  if (!bMark[Mesh[i][j]]) 
		  { 
			  nvex++; 
			  bMark[Mesh[i][j]] = TRUE; 
		  } 
 
  for (i=0;i<m_nVertex;i++) 
	  if (vex_tri[i]) 
		free(vex_tri[i]); 
  free(vex_tri); 
 
  tmp=(int)(m_nTriangle*ratio); 
 
  if (m_pProgress) 
	  m_pProgress->SetRange(0,tmp); 
   
  pair_contract(m_nVertex,border,cad_pair,pair_cost,pair_pos, 
	        vex_Q,Vertex,Vex_map,tmp); 
 
  free(vex_Q); 
  delete []border; 
 
  nMap = new int[m_nVertex]; 
  mMap = new int[m_nVertex]; 
  for (i=0;i<m_nVertex;i++) 
	  bMark[i]=FALSE; 
  for (i=0;i<m_nVertex;i++) 
  { 
	tmp=i; 
	while (Vex_map[tmp]!=tmp) 
		tmp=Vex_map[tmp]; 
	Vex_map[i]=tmp; 
  } 
  m_nsTriangle=0; 
  for (i=0;i<m_nTriangle;i++) 
	  if ((Vex_map[Mesh[i][0]]!=Vex_map[Mesh[i][1]])&& 
		  (Vex_map[Mesh[i][0]]!=Vex_map[Mesh[i][2]])&& 
		  (Vex_map[Mesh[i][1]]!=Vex_map[Mesh[i][2]])) 
	  { 
		  bMark[Vex_map[Mesh[i][0]]] = TRUE; 
		  bMark[Vex_map[Mesh[i][1]]] = TRUE; 
		  bMark[Vex_map[Mesh[i][2]]] = TRUE; 
		  m_nsTriangle++; 
	  } 
  m_nsVertex=0; 
  for (i=0;i<m_nVertex;i++) 
	  if (bMark[i]) 
	  { 
		  nMap[i] = m_nsVertex; 
		  mMap[m_nsVertex] = i; 
		  m_nsVertex++; 
	  } 
  m_sVertex = new Coord3[m_nsVertex]; 
  for (i=0;i<m_nsVertex;i++) 
  { 
	  VEC3_ASN_OP(m_sVertex[i],=,Vertex[mMap[i]]); 
  } 
  m_sTriangle = new Triangle[m_nsTriangle]; 
  m_nsTriangle = 0; 
  for (i=0;i<m_nTriangle;i++) 
	  if ((Vex_map[Mesh[i][0]]!=Vex_map[Mesh[i][1]])&& 
		  (Vex_map[Mesh[i][0]]!=Vex_map[Mesh[i][2]])&& 
		  (Vex_map[Mesh[i][1]]!=Vex_map[Mesh[i][2]])) 
	  { 
		  m_sTriangle[m_nsTriangle][0] = nMap[Vex_map[Mesh[i][0]]]; 
		  m_sTriangle[m_nsTriangle][1] = nMap[Vex_map[Mesh[i][1]]]; 
		  m_sTriangle[m_nsTriangle][2] = nMap[Vex_map[Mesh[i][2]]]; 
		  m_nsTriangle++; 
	  } 
 
	m_sVNormal = new Coord3[m_nsVertex]; 
	m_sTNormal = new Coord3[m_nsTriangle]; 
	for (i=0;i<m_nsVertex;i++) 
	{ 
		VEC3_ZERO(m_sVNormal[i]); 
	} 
	for (i=0;i<m_nsTriangle;i++) 
	{ 
		Coord3 vect[3]; 
		VEC3_V_OP_V(vect[0],m_sVertex[m_sTriangle[i][1]],-,m_sVertex[m_sTriangle[i][0]]); 
		VEC3_V_OP_V(vect[1],m_sVertex[m_sTriangle[i][2]],-,m_sVertex[m_sTriangle[i][0]]); 
		CROSSPROD3(vect[2],vect[0],vect[1]); 
		V3Standarize(vect[2]); 
		VEC3_ASN_OP(m_sTNormal[i],=,vect[2]); 
		for (j=0;j<3;j++) 
		{ 
			VEC3_V_OP_V(m_sVNormal[m_sTriangle[i][j]],m_sVNormal[m_sTriangle[i][j]],+,vect[2]); 
		} 
	} 
	for (i=0;i<m_nsVertex;i++) 
		V3Standarize(m_sVNormal[i]); 
 
  delete []Vex_map; 
  delete []Mesh; 
  delete []Vertex; 
  delete []mMap; 
  delete []nMap; 
  delete []bMark; 
 
  return m_nsTriangle; 
} 
 
int CMesh::Simplify(int nFaceNo) 
{ 
  Matrix4x4 *tri_Q; 
  Matrix4x4 *vex_Q; 
  Coord3 *Vertex; 
  Triangle *Mesh; 
   
  int **vex_tri; 
  int **vRing; 
  int *Vex_map; 
  int i,j,k; 
  int tmp,tmp1,tmp2; 
 
  int **cad_pair; 
  float **pair_cost; 
  Coord3 **pair_pos; 
 
  BOOL *bMark; 
  int  *nMap; 
  int  *mMap; 
 
  if (m_nTriangle == 0) 
	  return 0; 
 
  delete []m_sVertex; 
  delete []m_sTriangle; 
  delete []m_sVNormal; 
  delete []m_sTNormal; 
 
  Vertex= new Coord3[m_nVertex]; 
  for (i=0;i<m_nVertex;i++) 
  { 
	  VEC3_ASN_OP(Vertex[i],=,m_Vertex[i]); 
  } 
  Mesh = new Triangle[m_nTriangle]; 
  for (i=0; i<m_nTriangle; i++) 
  { 
	  VEC3_ASN_OP(Mesh[i],=,m_Triangle[i]); 
  } 
   
  vex_tri=(int **)malloc(m_nVertex*sizeof(int *)); 
  vRing=(int **)malloc(m_nVertex*sizeof(int *)); 
  for (i=0;i<m_nVertex;i++) { 
	vex_tri[i]=(int *)malloc(6*Size_int); 
	vex_tri[i][0]=0; 
	vRing[i]=(int *)malloc(6*Size_int); 
	vRing[i][0]=0; 
  } 
  cad_pair=(int **)malloc(m_nVertex*sizeof(int *)); 
  for (i=0;i<m_nVertex;i++) { 
	cad_pair[i]=(int *)malloc(6*Size_int); 
	cad_pair[i][0]=0; 
  } 
  pair_cost=(float **)malloc(m_nVertex*sizeof(float *)); 
  pair_pos=(Coord3 **)malloc(m_nVertex*sizeof(Coord3 *)); 
  Vex_map= new int[m_nVertex]; 
  for (i=0;i<m_nVertex;i++) 
	Vex_map[i]=i; 
 
  for (k=0;k<m_nTriangle;k++) 
  { 
	for (i=0;i<3;i++) { 
		tmp=tmp1=Mesh[k][i]; 
		tmp2=Mesh[k][(i+2)%3]; 
		for (j=1;j<vRing[tmp][0]+1;j++) 
			if (vRing[tmp][j] == tmp2) 
				break; 
		if (j==vRing[tmp][0]+1) 
		{ 
			vRing[tmp][j]=tmp2; 
			vRing[tmp][0] += 1; 
			if (!(vRing[tmp][0]%5)) 
				vRing[tmp] = (int *)realloc(vRing[tmp],(vRing[tmp][0]+6)*sizeof(int)); 
		} 
		tmp2=Mesh[k][(i+1)%3]; 
		for (j=1;j<vRing[tmp][0]+1;j++) 
			if (vRing[tmp][j] == tmp2) 
				break; 
		if (j==vRing[tmp][0]+1) 
		{ 
			vRing[tmp][j]=tmp2; 
			vRing[tmp][0] += 1; 
			if (!(vRing[tmp][0]%5)) 
				vRing[tmp] = (int *)realloc(vRing[tmp],(vRing[tmp][0]+6)*sizeof(int)); 
		} 
		if (vex_tri[tmp]==NULL) { 
			vex_tri[tmp]=(int *)malloc(6*sizeof(int)); 
			vex_tri[tmp][0]=1; 
			vex_tri[tmp][1]=k; 
		} 
		else { 
			if (!(vex_tri[tmp][0]%5)) 
				vex_tri[tmp]=(int *)realloc(vex_tri[tmp],(vex_tri[tmp][0]+6)*sizeof(int)); 
			vex_tri[tmp][0]=vex_tri[tmp][0]+1; 
			vex_tri[tmp][vex_tri[tmp][0]]=k; 
		} 
		if (tmp1>tmp2) { 
			tmp=tmp1; 
			tmp1=tmp2; 
			tmp2=tmp; 
		} 
		for (j=1;j=cad_pair[tmp1][0];j++) 
			if (cad_pair[tmp1][j]==tmp2) 
				break; 
		if (j>cad_pair[tmp1][0]) { 
			if (!(cad_pair[tmp1][0]%5)) 
				cad_pair[tmp1]=(int *)realloc(cad_pair[tmp1],(cad_pair[tmp1][0]+6)*sizeof(int)); 
			cad_pair[tmp1][0]=j; 
			cad_pair[tmp1][j]=tmp2; 
		} 
	} 
 
  } 
 
  BOOL *border; 
  border = new BOOL[m_nVertex]; 
  for (i=0;i<m_nVertex;i++) 
  { 
	  if (vex_tri[i][0]==vRing[i][0]) 
		  border[i]=FALSE; 
	  else 
		  border[i]=TRUE; 
	  free(vRing[i]); 
  } 
  free(vRing); 
 
  tri_Q=(Matrix4x4 *)malloc(m_nTriangle*Size_Matrix4x4); 
  cal_tri_Q(m_nTriangle,Mesh,Vertex,m_TNormal,tri_Q); 
  vex_Q=(Matrix4x4 *)malloc(m_nVertex*Size_Matrix4x4); 
  cal_vex_Q(m_nVertex,border,vex_tri,Mesh,Vertex,m_TNormal,tri_Q,vex_Q); 
  free(tri_Q); 
 
 
  bMark = new BOOL[m_nVertex]; 
  for (i=0;i<m_nVertex;i++) 
	  bMark[i]=FALSE; 
  int nvex = 0; 
  for (i=0;i<m_nTriangle;i++) 
	  for (j=0;j<3;j++) 
		  if (!bMark[Mesh[i][j]]) 
		  { 
			  nvex++; 
			  bMark[Mesh[i][j]] = TRUE; 
		  } 
 
  for (i=0;i<m_nVertex;i++) 
	  if (vex_tri[i]) 
		free(vex_tri[i]); 
  free(vex_tri); 
 
  tmp = m_nTriangle - nFaceNo; 
 
  if (m_pProgress) 
	  m_pProgress->SetRange(0,tmp); 
 
  pair_contract(m_nVertex,border,cad_pair,pair_cost,pair_pos, 
	        vex_Q,Vertex,Vex_map,tmp); 
 
  free(vex_Q); 
  delete []border; 
 
  nMap = new int[m_nVertex]; 
  mMap = new int[m_nVertex]; 
  for (i=0;i<m_nVertex;i++) 
	  bMark[i]=FALSE; 
  for (i=0;i<m_nVertex;i++) 
  { 
	tmp=i; 
	while (Vex_map[tmp]!=tmp) 
		tmp=Vex_map[tmp]; 
	Vex_map[i]=tmp; 
  } 
  m_nsTriangle=0; 
  for (i=0;i<m_nTriangle;i++) 
	  if ((Vex_map[Mesh[i][0]]!=Vex_map[Mesh[i][1]])&& 
		  (Vex_map[Mesh[i][0]]!=Vex_map[Mesh[i][2]])&& 
		  (Vex_map[Mesh[i][1]]!=Vex_map[Mesh[i][2]])) 
	  { 
		  bMark[Vex_map[Mesh[i][0]]] = TRUE; 
		  bMark[Vex_map[Mesh[i][1]]] = TRUE; 
		  bMark[Vex_map[Mesh[i][2]]] = TRUE; 
		  m_nsTriangle++; 
	  } 
  m_nsVertex=0; 
  for (i=0;i<m_nVertex;i++) 
	  if (bMark[i]) 
	  { 
		  nMap[i] = m_nsVertex; 
		  mMap[m_nsVertex] = i; 
		  m_nsVertex++; 
	  } 
  m_sVertex = new Coord3[m_nsVertex]; 
  for (i=0;i<m_nsVertex;i++) 
  { 
	  VEC3_ASN_OP(m_sVertex[i],=,Vertex[mMap[i]]); 
  } 
  m_sTriangle = new Triangle[m_nsTriangle]; 
  m_nsTriangle = 0; 
  for (i=0;i<m_nTriangle;i++) 
	  if ((Vex_map[Mesh[i][0]]!=Vex_map[Mesh[i][1]])&& 
		  (Vex_map[Mesh[i][0]]!=Vex_map[Mesh[i][2]])&& 
		  (Vex_map[Mesh[i][1]]!=Vex_map[Mesh[i][2]])) 
	  { 
		  m_sTriangle[m_nsTriangle][0] = nMap[Vex_map[Mesh[i][0]]]; 
		  m_sTriangle[m_nsTriangle][1] = nMap[Vex_map[Mesh[i][1]]]; 
		  m_sTriangle[m_nsTriangle][2] = nMap[Vex_map[Mesh[i][2]]]; 
		  m_nsTriangle++; 
	  } 
 
	m_sVNormal = new Coord3[m_nsVertex]; 
	m_sTNormal = new Coord3[m_nsTriangle]; 
	for (i=0;i<m_nsVertex;i++) 
	{ 
		VEC3_ZERO(m_sVNormal[i]); 
	} 
	for (i=0;i<m_nsTriangle;i++) 
	{ 
		Coord3 vect[3]; 
		VEC3_V_OP_V(vect[0],m_sVertex[m_sTriangle[i][1]],-,m_sVertex[m_sTriangle[i][0]]); 
		VEC3_V_OP_V(vect[1],m_sVertex[m_sTriangle[i][2]],-,m_sVertex[m_sTriangle[i][0]]); 
		CROSSPROD3(vect[2],vect[0],vect[1]); 
		V3Standarize(vect[2]); 
		VEC3_ASN_OP(m_sTNormal[i],=,vect[2]); 
		for (j=0;j<3;j++) 
		{ 
			VEC3_V_OP_V(m_sVNormal[m_sTriangle[i][j]],m_sVNormal[m_sTriangle[i][j]],+,vect[2]); 
		} 
	} 
	for (i=0;i<m_nsVertex;i++) 
		V3Standarize(m_sVNormal[i]); 
 
  delete []Vex_map; 
  delete []Mesh; 
  delete []Vertex; 
  delete []mMap; 
  delete []nMap; 
  delete []bMark; 
 
   return m_nsTriangle; 
} 
 
void CMesh::cal_tri_Q(int Num_triangle,Triangle *mesh,Coord3 *vertex,Vector3 *tri_normal,Matrix4x4 *tri_Q) 
{ 
  int i; 
  float tmp; 
 
  for (i=0;i<Num_triangle;i++) { 
	tmp=-DOTPROD3(tri_normal[i],vertex[mesh[i][0]]); 
	tri_Q[i][0][0]=tri_normal[i][0]*tri_normal[i][0]; 
	tri_Q[i][1][1]=tri_normal[i][1]*tri_normal[i][1]; 
	tri_Q[i][2][2]=tri_normal[i][2]*tri_normal[i][2]; 
	tri_Q[i][3][3]=tmp*tmp; 
	tri_Q[i][0][1]=tri_Q[i][1][0]=tri_normal[i][0]*tri_normal[i][1]; 
	tri_Q[i][0][2]=tri_Q[i][2][0]=tri_normal[i][0]*tri_normal[i][2]; 
	tri_Q[i][0][3]=tri_Q[i][3][0]=tri_normal[i][0]*tmp; 
	tri_Q[i][1][2]=tri_Q[i][2][1]=tri_normal[i][1]*tri_normal[i][2]; 
	tri_Q[i][1][3]=tri_Q[i][3][1]=tri_normal[i][1]*tmp; 
	tri_Q[i][2][3]=tri_Q[i][3][2]=tri_normal[i][2]*tmp; 
  } 
} 
 
void CMesh::cal_vex_Q(int Num_vertex,BOOL *border,int **vex_tri,Triangle *mesh,Coord3 *vertex,Vector3 *tri_normal,Matrix4x4 *tri_Q,Matrix4x4 *vex_Q) 
{ 
  int i,j,k,m,p; 
  Matrix4x4 tmp_Q; 
  float tmp; 
  Vector3 vect[2]; 
 
  for (i=0;i<Num_vertex;i++) { 
	for (j=0;j<4;j++) 
		VEC4_ZERO(vex_Q[i][j]); 
	for (j=1;j=vex_tri[i][0];j++) { 
		p=vex_tri[i][j]; 
		for (k=0;k<4;k++) 
			VEC4_V_OP_V(vex_Q[i][k],vex_Q[i][k],+,tri_Q[p][k]); 
		for (k=0;k<3;k++) 
			if (mesh[p][k]==i) 
				break; 
		int a,b; 
		if (mesh[p][k]>mesh[p][(k+2)%3]) 
		{ 
			a=mesh[p][k]; 
			b=mesh[p][(k+2)%3]; 
		} 
		else 
		{ 
			b=mesh[p][k]; 
			a=mesh[p][(k+2)%3]; 
		} 
		if (border[a]&&border[b]) { 
			VEC3_V_OP_V(vect[0],vertex[mesh[p][(k+2)%3]],-,vertex[i]); 
			V3Standarize(vect[0]); 
			CROSSPROD3(vect[1],vect[0],tri_normal[p]); 
			V3Standarize(vect[1]); 
			tmp=-DOTPROD3(vect[1],vertex[i]); 
			tmp_Q[0][0]=vect[1][0]*vect[1][0]; 
			tmp_Q[1][1]=vect[1][1]*vect[1][1]; 
			tmp_Q[2][2]=vect[1][2]*vect[1][2]; 
			tmp_Q[3][3]=tmp*tmp; 
			tmp_Q[0][1]=tmp_Q[1][0]=vect[1][0]*vect[1][1]; 
			tmp_Q[0][2]=tmp_Q[2][0]=vect[1][0]*vect[1][2]; 
			tmp_Q[0][3]=tmp_Q[3][0]=vect[1][0]*tmp; 
			tmp_Q[1][2]=tmp_Q[2][1]=vect[1][1]*vect[1][2]; 
			tmp_Q[1][3]=tmp_Q[3][1]=vect[1][1]*tmp; 
			tmp_Q[2][3]=tmp_Q[3][2]=vect[1][2]*tmp; 
			for (m=0;m<4;m++) 
				VEC4_V_OP_V(vex_Q[i][m],vex_Q[i][m],+,tmp_Q[m]); 
		} 
		if (mesh[p][k]>mesh[p][(k+1)%3]) 
		{ 
			a=mesh[p][k]; 
			b=mesh[p][(k+1)%3]; 
		} 
		else 
		{ 
			b=mesh[p][k]; 
			a=mesh[p][(k+1)%3]; 
		} 
		if (border[a]&&border[b]) { 
			VEC3_V_OP_V(vect[0],vertex[mesh[p][(k+1)%3]],-,vertex[i]); 
			V3Standarize(vect[0]); 
			CROSSPROD3(vect[1],vect[0],tri_normal[p]); 
			V3Standarize(vect[1]); 
			tmp=-DOTPROD3(vect[1],vertex[i]); 
			tmp_Q[0][0]=vect[1][0]*vect[1][0]; 
			tmp_Q[1][1]=vect[1][1]*vect[1][1]; 
			tmp_Q[2][2]=vect[1][2]*vect[1][2]; 
			tmp_Q[3][3]=tmp*tmp; 
			tmp_Q[0][1]=tmp_Q[1][0]=vect[1][0]*vect[1][1]; 
			tmp_Q[0][2]=tmp_Q[2][0]=vect[1][0]*vect[1][2]; 
			tmp_Q[0][3]=tmp_Q[3][0]=vect[1][0]*tmp; 
			tmp_Q[1][2]=tmp_Q[2][1]=vect[1][1]*vect[1][2]; 
			tmp_Q[1][3]=tmp_Q[3][1]=vect[1][1]*tmp; 
			tmp_Q[2][3]=tmp_Q[3][2]=vect[1][2]*tmp; 
			for (m=0;m<4;m++) 
				VEC4_V_OP_V(vex_Q[i][m],vex_Q[i][m],+,tmp_Q[m]); 
		} 
	} 
  } 
} 
 
void CMesh::pair_contract(int Num_vertex,BOOL *border,int **cad_pair,float **pair_cost,Coord3 **pair_pos,Matrix4x4 *vex_Q,Coord3 *vertex,int *vex_map,int require) 
{ 
  int i,j,k,m,n,min_i,min_j; 
  float **tmp_Q; 
  float M3x3[3][3]; 
  float vec_b[4]; 
  int remain,deletenum; 
  float min_cost; 
 
  tmp_Q=(float **)malloc(4*sizeof(float *)); 
  for (i=0;i<4;i++) 
	tmp_Q[i]=(float *)malloc(4*Size_float); 
 
  for (i=0;i<Num_vertex;i++) { 
	pair_cost[i]=(float *)malloc(cad_pair[i][0]*Size_float); 
	pair_pos[i]=(Coord3 *)malloc(cad_pair[i][0]*Size_Coord3); 
	for (k=1;k=cad_pair[i][0];k++) { 
		j=cad_pair[i][k]; 
		for (m=0;m<3;m++) 
			VEC4_V_OP_V(tmp_Q[m],vex_Q[i][m],+,vex_Q[j][m]); 
		tmp_Q[3][0]=tmp_Q[3][1]=tmp_Q[3][2]=0.0; 
		tmp_Q[3][3]=1.0; 
		vec_b[0]=vec_b[1]=vec_b[2]=0.0; 
		vec_b[3]=1.0; 
		for (m=0;m<3;m++) 
			VEC3_ASN_OP(M3x3[m],=,tmp_Q[m]); 
		if ((fabs(deta3(M3x3))>=0.01)&&(Gauss(4,tmp_Q,vec_b))) { 
			VEC3_ASN_OP(pair_pos[i][k-1],=,vec_b); 
			for (m=0;m<4;m++) 
				VEC4_V_OP_V(tmp_Q[m],vex_Q[i][m],+,vex_Q[j][m]); 
			pair_cost[i][k-1]=tmp_Q[0][0]*vec_b[0]*vec_b[0]+ 
			                  2*tmp_Q[0][1]*vec_b[0]*vec_b[1]+ 
			                  2*tmp_Q[0][2]*vec_b[0]*vec_b[2]+ 
			                  2*tmp_Q[0][3]*vec_b[0]+ 
			                  tmp_Q[1][1]*vec_b[1]*vec_b[1]+ 
			                  2*tmp_Q[1][2]*vec_b[1]*vec_b[2]+ 
			                  2*tmp_Q[1][3]*vec_b[1]+ 
			                  tmp_Q[2][2]*vec_b[2]*vec_b[2]+ 
			                  2*tmp_Q[2][3]*vec_b[2]+ 
					  tmp_Q[3][3]; 
		} 
		else { 
 
			for (m=0;m<4;m++) 
				VEC4_V_OP_V(tmp_Q[m],vex_Q[i][m],+,vex_Q[j][m]); 
			VEC3_ASN_OP(vec_b,=,vertex[i]); 
			VEC3_ASN_OP(pair_pos[i][k-1],=,vec_b); 
			vec_b[3]=tmp_Q[0][0]*vec_b[0]*vec_b[0]+ 
			                  2*tmp_Q[0][1]*vec_b[0]*vec_b[1]+ 
			                  2*tmp_Q[0][2]*vec_b[0]*vec_b[2]+ 
			                  2*tmp_Q[0][3]*vec_b[0]+ 
			                  tmp_Q[1][1]*vec_b[1]*vec_b[1]+ 
			                  2*tmp_Q[1][2]*vec_b[1]*vec_b[2]+ 
			                  2*tmp_Q[1][3]*vec_b[1]+ 
			                  tmp_Q[2][2]*vec_b[2]*vec_b[2]+ 
			                  2*tmp_Q[2][3]*vec_b[2]+ 
					  tmp_Q[3][3]; 
			min_cost=vec_b[3]; 
			VEC3_ASN_OP(vec_b,=,vertex[j]); 
			vec_b[3]=tmp_Q[0][0]*vec_b[0]*vec_b[0]+ 
			                  2*tmp_Q[0][1]*vec_b[0]*vec_b[1]+ 
			                  2*tmp_Q[0][2]*vec_b[0]*vec_b[2]+ 
			                  2*tmp_Q[0][3]*vec_b[0]+ 
			                  tmp_Q[1][1]*vec_b[1]*vec_b[1]+ 
			                  2*tmp_Q[1][2]*vec_b[1]*vec_b[2]+ 
			                  2*tmp_Q[1][3]*vec_b[1]+ 
			                  tmp_Q[2][2]*vec_b[2]*vec_b[2]+ 
			                  2*tmp_Q[2][3]*vec_b[2]+ 
					  tmp_Q[3][3]; 
			if (min_cost>vec_b[3]) { 
				min_cost=vec_b[3]; 
				VEC3_ASN_OP(pair_pos[i][k-1],=,vec_b); 
			} 
			VEC3_V_OP_V(vec_b,vertex[i],+,vertex[j]); 
			VEC3_V_OP_S(vec_b,vec_b,/,2.0f); 
			vec_b[3]=tmp_Q[0][0]*vec_b[0]*vec_b[0]+ 
			                  2*tmp_Q[0][1]*vec_b[0]*vec_b[1]+ 
			                  2*tmp_Q[0][2]*vec_b[0]*vec_b[2]+ 
			                  2*tmp_Q[0][3]*vec_b[0]+ 
			                  tmp_Q[1][1]*vec_b[1]*vec_b[1]+ 
			                  2*tmp_Q[1][2]*vec_b[1]*vec_b[2]+ 
			                  2*tmp_Q[1][3]*vec_b[1]+ 
			                  tmp_Q[2][2]*vec_b[2]*vec_b[2]+ 
			                  2*tmp_Q[2][3]*vec_b[2]+ 
					  tmp_Q[3][3]; 
			if (min_cost>vec_b[3]) { 
				min_cost=vec_b[3]; 
				VEC3_ASN_OP(pair_pos[i][k-1],=,vec_b); 
			} 
			pair_cost[i][k-1]=min_cost; 
		} 
	} 
  } 
 
 
  // 99,6,1 
  int numEdge = 0; 
  Heap *heap; 
  Heapable *pair; 
  int **VexEdge; 
  //VexEdge = new int51[Num_vertex]; 
  VexEdge = new Pint[Num_vertex]; 
  for (i=0;i<Num_vertex;i++) 
  { 
	  VexEdge[i] = (int *)malloc(51*sizeof(int)); 
	  VexEdge[i][0] = 0; 
	  numEdge += cad_pair[i][0]; 
  } 
  heap = new Heap(numEdge); 
  pair = new Heapable[numEdge]; 
 
  k = 0; 
  for (i=0;i<Num_vertex;i++) 
	for (j=1; j=cad_pair[i][0]; j++) 
	{ 
 
		pair[k].vert[0] = i; 
		m = pair[k].vert[1] = cad_pair[i][j]; 
		pair[k].cost = -pair_cost[i][j-1]; 
		VEC3_ASN_OP(pair[k].newposition,=,pair_pos[i][j-1]) 
		heap->insert(&(pair[k]), pair[k].cost); 
		if (!(pair[k].isInHeap())) 
			break; 
		n = VexEdge[i][0]+1; 
		if (n>50) 
		{ 
			VexEdge[i] = (int *)realloc(VexEdge[i],(n+50)*sizeof(int)); 
		} 
		VexEdge[i][n] = k; 
		VexEdge[i][0] = n; 
		n = VexEdge[m][0]+1; 
		if (n>50) 
		{ 
			VexEdge[m] = (int *)realloc(VexEdge[m],(n+50)*sizeof(int)); 
		} 
		VexEdge[m][n] = k; 
		VexEdge[m][0] = n; 
		k++; 
	} 
 
  for (i=0;i<Num_vertex;i++) 
  { 
	  if (pair_cost[i]) 
		  free(pair_cost[i]); 
	  if (pair_pos[i]) 
		  free(pair_pos[i]); 
	  if (cad_pair[i]) 
		  free(cad_pair[i]); 
  } 
  free(pair_cost); 
  free(pair_pos); 
  free(cad_pair); 
 
  // 99,6,1 
 
  deletenum=0; 
  remain=Num_vertex;   
  int delEdge = 0; 
   
  while (1) { 
 
	heap_node *top; 
	Heapable *dpair; 
 
	top = heap->extract(); 
	if (!top) 
		break; 
	dpair = top->obj; 
 
	min_i = dpair->vert[0]; 
	min_j = dpair->vert[1]; 
	if (border[min_i]&&border[min_j]) 
		deletenum ++; 
	else 
		deletenum +=2; 
	border[min_i] = border[min_i]||border[min_j]; 
 
	VEC3_ASN_OP(vertex[min_i],=,dpair->newposition); 
	VEC3_ASN_OP(vertex[min_j],=,dpair->newposition); 
	vex_map[min_j]=min_i; 
 
	delEdge++; 
 
	for (i=1;i=VexEdge[min_i][0];i++) 
	{ 
		m = VexEdge[min_i][i]; 
		if ((pair[m].vert[0] == min_j)||(pair[m].vert[1] == min_j)) 
		{ 
			for (j=i; j<VexEdge[min_i][0]; j++) 
				VexEdge[min_i][j] = VexEdge[min_i][j+1]; 
			VexEdge[min_i][0] = VexEdge[min_i][0]-1; 
			break; 
		} 
	} 
 
	for (i=1;i=VexEdge[min_j][0];i++) 
	{ 
		int p,q; 
 
		m = VexEdge[min_j][i]; 
		if ((pair[m].vert[0] == min_i)||(pair[m].vert[1] == min_i)) 
			continue; 
		if (pair[m].vert[0] == min_j) 
		{ 
			pair[m].vert[0] = min_i; 
			for (j=1;j=VexEdge[min_i][0]; j++) 
			{ 
				n = VexEdge[min_i][j]; 
				if (pair[n].vert[0] == pair[m].vert[1]) 
				{ 
					k = pair[n].vert[0]; 
					heap->kill(pair[m].getHeapPos()); 
					for (p=1; p=VexEdge[k][0]; p++) 
					{ 
						if (VexEdge[k][p] == m) 
						{ 
							for (q=p; q<VexEdge[k][0]; q++) 
								VexEdge[k][q]=VexEdge[k][q+1]; 
							VexEdge[k][0] = VexEdge[k][0]-1; 
							break; 
						} 
					} 
					break; 
				} 
				else if (pair[n].vert[1] == pair[m].vert[1]) 
				{ 
					k = pair[n].vert[1]; 
					heap->kill(pair[m].getHeapPos()); 
					for (p=1; p=VexEdge[k][0]; p++) 
					{ 
						if (VexEdge[k][p] == m) 
						{ 
							for (q=p; q<VexEdge[k][0]; q++) 
								VexEdge[k][q]=VexEdge[k][q+1]; 
							VexEdge[k][0] = VexEdge[k][0]-1; 
							break; 
						} 
					} 
					break; 
				} 
			} 
 
			if (j>VexEdge[min_i][0]) 
			{ 
				if (j>50) 
				{ 
					VexEdge[min_i] = (int *)realloc(VexEdge[min_i],(j+50)*sizeof(int)); 
				} 
				VexEdge[min_i][0] = j; 
				VexEdge[min_i][j] = m; 
			} 
		} 
		else if (pair[m].vert[1] == min_j) 
		{ 
			pair[m].vert[1] = min_i; 
			for (j=1;j=VexEdge[min_i][0]; j++) 
			{ 
				n = VexEdge[min_i][j]; 
				if (pair[n].vert[0] == pair[m].vert[0]) 
				{ 
					k = pair[n].vert[0]; 
					heap->kill(pair[m].getHeapPos()); 
					for (p=1; p=VexEdge[k][0]; p++) 
					{ 
						if (VexEdge[k][p] == m) 
						{ 
							for (q=p; q<VexEdge[k][0]; q++) 
								VexEdge[k][q]=VexEdge[k][q+1]; 
							VexEdge[k][0] = VexEdge[k][0]-1; 
							break; 
						} 
					} 
					break; 
				} 
				else if (pair[n].vert[1] == pair[m].vert[0]) 
				{ 
					k = pair[n].vert[1]; 
					heap->kill(pair[m].getHeapPos()); 
					for (p=1; p=VexEdge[k][0]; p++) 
					{ 
						if (VexEdge[k][p] == m) 
						{ 
							for (q=p; q<VexEdge[k][0]; q++) 
								VexEdge[k][q]=VexEdge[k][q+1]; 
							VexEdge[k][0] = VexEdge[k][0]-1; 
							break; 
						} 
					} 
					break; 
				} 
			} 
 
			if (j>VexEdge[min_i][0]) 
			{ 
				if (j>50) 
				{ 
					VexEdge[min_i] = (int *)realloc(VexEdge[min_i],(j+50)*sizeof(int)); 
				} 
				VexEdge[min_i][0] = j; 
				VexEdge[min_i][j] = m; 
			} 
		} 
	} 
 
	for (i=0;i<4;i++) 
		VEC4_V_OP_V(vex_Q[min_i][i],vex_Q[min_i][i],+,vex_Q[min_j][i]); 
 
	for (i=1; i=VexEdge[min_i][0];i++) 
	{ 
		n = VexEdge[min_i][i]; 
		j = pair[n].vert[0]; 
		if (j == min_i ) 
			j = pair[n].vert[1]; 
		for (m=0;m<3;m++) 
			VEC4_V_OP_V(tmp_Q[m],vex_Q[min_i][m],+,vex_Q[j][m]); 
		tmp_Q[3][0]=tmp_Q[3][1]=tmp_Q[3][2]=0.0; 
		tmp_Q[3][3]=1.0; 
		vec_b[0]=vec_b[1]=vec_b[2]=0.0; 
		vec_b[3]=1.0; 
		for (m=0;m<3;m++) 
			VEC3_ASN_OP(M3x3[m],=,tmp_Q[m]); 
		if ((fabs(deta3(M3x3))>=0.01)&&(Gauss(4,tmp_Q,vec_b))) { 
			VEC3_ASN_OP(pair[n].newposition,=,vec_b); 
			for (m=0;m<4;m++) 
				VEC4_V_OP_V(tmp_Q[m],vex_Q[min_i][m],+,vex_Q[j][m]); 
			pair[n].cost  =   tmp_Q[0][0]*vec_b[0]*vec_b[0]+ 
			                  2*tmp_Q[0][1]*vec_b[0]*vec_b[1]+ 
			                  2*tmp_Q[0][2]*vec_b[0]*vec_b[2]+ 
			                  2*tmp_Q[0][3]*vec_b[0]+ 
			                  tmp_Q[1][1]*vec_b[1]*vec_b[1]+ 
			                  2*tmp_Q[1][2]*vec_b[1]*vec_b[2]+ 
			                  2*tmp_Q[1][3]*vec_b[1]+ 
			                  tmp_Q[2][2]*vec_b[2]*vec_b[2]+ 
			                  2*tmp_Q[2][3]*vec_b[2]+ 
							  tmp_Q[3][3]; 
		} 
		else { 
 
			for (m=0;m<4;m++) 
				VEC4_V_OP_V(tmp_Q[m],vex_Q[min_i][m],+,vex_Q[j][m]); 
			VEC3_ASN_OP(vec_b,=,vertex[min_i]); 
			VEC3_ASN_OP(pair[n].newposition,=,vec_b); 
			vec_b[3]=tmp_Q[0][0]*vec_b[0]*vec_b[0]+ 
			                  2*tmp_Q[0][1]*vec_b[0]*vec_b[1]+ 
			                  2*tmp_Q[0][2]*vec_b[0]*vec_b[2]+ 
			                  2*tmp_Q[0][3]*vec_b[0]+ 
			                  tmp_Q[1][1]*vec_b[1]*vec_b[1]+ 
			                  2*tmp_Q[1][2]*vec_b[1]*vec_b[2]+ 
			                  2*tmp_Q[1][3]*vec_b[1]+ 
			                  tmp_Q[2][2]*vec_b[2]*vec_b[2]+ 
			                  2*tmp_Q[2][3]*vec_b[2]+ 
							  tmp_Q[3][3]; 
			min_cost=vec_b[3]; 
			VEC3_ASN_OP(vec_b,=,vertex[j]); 
			vec_b[3]=tmp_Q[0][0]*vec_b[0]*vec_b[0]+ 
			                  2*tmp_Q[0][1]*vec_b[0]*vec_b[1]+ 
			                  2*tmp_Q[0][2]*vec_b[0]*vec_b[2]+ 
			                  2*tmp_Q[0][3]*vec_b[0]+ 
			                  tmp_Q[1][1]*vec_b[1]*vec_b[1]+ 
			                  2*tmp_Q[1][2]*vec_b[1]*vec_b[2]+ 
			                  2*tmp_Q[1][3]*vec_b[1]+ 
			                  tmp_Q[2][2]*vec_b[2]*vec_b[2]+ 
			                  2*tmp_Q[2][3]*vec_b[2]+ 
							  tmp_Q[3][3]; 
			if (min_cost>vec_b[3]) { 
				min_cost=vec_b[3]; 
				VEC3_ASN_OP(pair[n].newposition,=,vec_b); 
			} 
			VEC3_V_OP_V(vec_b,vertex[min_i],+,vertex[j]); 
			VEC3_V_OP_S(vec_b,vec_b,/,2.0); 
			vec_b[3]=tmp_Q[0][0]*vec_b[0]*vec_b[0]+ 
			                  2*tmp_Q[0][1]*vec_b[0]*vec_b[1]+ 
			                  2*tmp_Q[0][2]*vec_b[0]*vec_b[2]+ 
			                  2*tmp_Q[0][3]*vec_b[0]+ 
			                  tmp_Q[1][1]*vec_b[1]*vec_b[1]+ 
			                  2*tmp_Q[1][2]*vec_b[1]*vec_b[2]+ 
			                  2*tmp_Q[1][3]*vec_b[1]+ 
			                  tmp_Q[2][2]*vec_b[2]*vec_b[2]+ 
			                  2*tmp_Q[2][3]*vec_b[2]+ 
							  tmp_Q[3][3]; 
			if (min_cost>vec_b[3]) { 
				min_cost=vec_b[3]; 
				VEC3_ASN_OP(pair[n].newposition,=,vec_b); 
			} 
			pair[n].cost=min_cost; 
		} 
		pair[n].cost = -pair[n].cost; 
		heap->update(&(pair[n]),pair[n].cost); 
	} 
 
	if (m_pProgress) 
		m_pProgress->SetPos(deletenum); 
 
	if (deletenum >= require) 
		break; 
 
  } /* while (1) */ 
 
  delete heap; 
  delete []pair; 
  for (i=0;i<Num_vertex;i++) 
  { 
	  free(VexEdge[i]); 
  } 
  delete []VexEdge; 
 
  for (i=0;i<4;i++) 
	free(tmp_Q[i]); 
  free(tmp_Q); 
} 
 
void CMesh::V3Standarize(Vector3 v) 
{ 
	float len; 
	len=sqrt(DOTPROD3(v,v)); 
	if (len!=0.0) 
	{ 
		v[X]=v[X]/len; 
		v[Y]=v[Y]/len; 
		v[Z]=v[Z]/len; 
	} 
} 
 
float CMesh::deta3(float matrix[3][3]) 
{ 
  float deta; 
  deta=matrix[0][0]*(matrix[1][1]*matrix[2][2]-matrix[1][2]*matrix[2][1])- 
       matrix[0][1]*(matrix[1][0]*matrix[2][2]-matrix[1][2]*matrix[2][0])+ 
       matrix[0][2]*(matrix[1][0]*matrix[2][1]-matrix[1][1]*matrix[2][0]); 
  return deta; 
} 
 
int CMesh::Gauss(int n,float **A,float *b) 
{ 
  int *js,i,j,k,is; 
  float d,t; 
 
  js=(int *)malloc(n*sizeof(int)); 
  for (k=0;k<n;k++) { 
	d=0.0; 
	for (i=k;i<n;i++) 
		for (j=k;j<n;j++) 
			if (A[i][j]!=0.0) { 
				t=fabs(A[i][j]); 
				if (t>d) { 
					d=t; 
					js[k]=j; 
					is=i; 
				} 
			} 
 
	if (d+1.0==1.0) { 
		free(js); 
		return 0; 
	} 
 
	if (is!=k) { 
		for (j=k;j<n;j++) { 
			t=A[k][j]; 
			A[k][j]=A[is][j]; 
			A[is][j]=t; 
						 
		} 
		t=b[k]; 
		b[k]=b[is]; 
		b[is]=t; 
	} 
 
	if (js[k]!=k) 
		for (i=0;i<n;i++) { 
			t=A[i][k]; 
			A[i][k]=A[i][js[k]]; 
			A[i][js[k]]=t; 
		} 
 
	t=A[k][k]; 
	for (j=k+1;j<n;j++) 
		if (A[k][j]!=0.0) 
			A[k][j]/=t; 
 
	b[k]=b[k]/t; 
	for (j=k+1;j<n;j++) 
		if (A[k][j]!=0.0) 
			for (i=0;i<n;i++) 
				if ((i!=k)&&(A[i][k]!=0.0)) 
					A[i][j]=A[i][j]-A[i][k]*A[k][j]; 
 
	for (i=0;i<n;i++) 
		if ((i!=k)&&(A[i][k]!=0.0)) 
			b[i]=b[i]-A[i][k]*b[k]; 
 
  } 
 
  for (k=n-1;k>=0;k--) 
	if (k!=js[k]) { 
		t=b[k]; 
		b[k]=b[js[k]]; 
		b[js[k]]=t; 
	} 
 
  free(js); 
 
  return 1; 
}