/************************************************************************
Yutaka Ohtake 

Simplifier.cpp
Simplification methods

Copyright (c) 1999-2001 
The University of Aizu. All Rights Reserved.
************************************************************************/

#include "stdafx.h"
#include "MeshEditor.h"
#include "Simplifier.h"

/*
#include "vtkDelaunay2D.h"
#include "vtkPolygon.h"
#include "vtkPoints.h"
#include "vtkCellType.h"
*/

#define PI 3.14159265

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// �\�z/����
//////////////////////////////////////////////////////////////////////

Simplifier::Simplifier()
{
	mesh = NULL;
	tag = NULL;
	//ridge_edge = NULL;
	//ravine_edge = NULL;
	both_edge = NULL;
}

Simplifier::~Simplifier()
{
	if(tag != NULL)
		delete[] tag;
	if(both_edge != NULL){
		for(int i=0; i<mesh->vertex_N; i++)
			delete both_edge[i];
		delete[] both_edge;
	}
}

BOOL Simplifier::vertexRemove(int index)
{
	/*
	Node** v_link = mesh->vertex_link;
	int (*f_link)[3] = mesh->face_link_E;
	float (*vertex)[3] = mesh->vertex;
	int (*face)[3] = mesh->face;

	int degree;
	int *star;
	int *deleted_face;

	mesh->getConsistent1Ring(index, star, deleted_face, degree);

	for(int i=0; i<degree; i++){
		mesh->deleteFace(deleted_face[i]);
	}

	for(i=0; i<degree; i++){
		Node *tmp_adjacent = new Node;
		tmp_adjacent->next = NULL;

		for(Node* current = v_link[star[i]]; current->next != NULL; current = current->next->next){
			if(current->v != index && current->next->v != index){
				tmp_adjacent->append(current->v, current->f);
				tmp_adjacent->append(current->next->v, current->next->f);
			}
		}
		delete v_link[star[i]];
		v_link[star[i]] = tmp_adjacent;
	}

	float total_angle = 0;
	if(mesh->isBound[index])
		for(i=0; i<degree-1; i++){
			float v1[3], v2[3];
			MeshData::VEC(v1, vertex[index], vertex[star[i]]);
			MeshData::VEC(v2, vertex[index], vertex[star[i+1]]);
			float l1 = MeshData::LENGTH(v1);
			float l2 = MeshData::LENGTH(v2);
			if(l1 != 0 && l2 != 0)
				total_angle += (float)acos(MeshData::DOT(v1, v2)/(l1*l2));
		}
	else
		for(i=0; i<degree; i++){
			float v1[3], v2[3];
			MeshData::VEC(v1, vertex[index], vertex[star[i]]);
			MeshData::VEC(v2, vertex[index], vertex[star[(i+1)%degree]]);
			float l1 = MeshData::LENGTH(v1);
			float l2 = MeshData::LENGTH(v2);
			if(l1 != 0 && l2 != 0)
				total_angle += (float)acos(MeshData::DOT(v1, v2)/(l1*l2));
		}

	float (*z)[2] = new float[degree][2];
	float S = 0;
	float a;
	if(!mesh->isBound[index])
		a = 2.0*PI/(double)total_angle;
	else{
		if(total_angle < 2.0*PI)
			a = 1.0;
		else
			a = 2.0*PI/total_angle;
	}
	float v1[3], v2[3];
	float r1, r2, max_r = -1;
	MeshData::VEC(v1, vertex[index], vertex[star[0]]);
	r1 = MeshData::LENGTH(v1);
	for(i=0; i<degree-1; i++){
		if(max_r < r1)
			max_r = r1;

		z[i][0] = (float)(r1*cos(S*a));
		z[i][1] = (float)(r1*sin(S*a));

		MeshData::VEC(v2, vertex[index], vertex[star[i+1]]);
		r2 = MeshData::LENGTH(v2);
		
		if(r1 != 0 && r2 != 0)
			S += (float)acos((v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2])/(r1*r2));

		v1[0] = v2[0];
		v1[1] = v2[1];
		v1[2] = v2[2];
		r1 = r2;
	}
	if(max_r < r1)
			max_r = r1;
	z[i][0] = (float)(r1*cos(S*a));
	z[i][1] = (float)(r1*sin(S*a));

	for(i=0; i<degree; i++){
		z[i][0] = z[i][0]*500.0f/max_r;
		z[i][1] = z[i][1]*500.0f/max_r;
	}
	
	vtkDelaunay2D* d2d = vtkDelaunay2D::New();
	vtkPoints* ps = vtkPoints::New();
	vtkPolyData* poly_data = vtkPolyData::New();

	ps->Initialize();
	int *id = new int[degree];
	for(i=0; i<degree; i++){
		id[i] = i;
		ps->InsertPoint(i, z[i][0], z[i][1], 0.0f);
	}
	poly_data->Initialize();
	poly_data->SetPoints(ps);
	poly_data->Allocate(1,1);
	poly_data->InsertNextCell(VTK_POLYGON, degree, id); 
	d2d->SetInput(poly_data);
	d2d->SetSource(poly_data);
	d2d->SetTolerance(0.00001);
	d2d->SetOffset(500.0);
	d2d->SetAlpha(0);
	d2d->BoundingTriangulationOff();
	d2d->Update();
	vtkCellArray* tri_data = d2d->GetOutput()->GetPolys();

	int tri_N = tri_data->GetNumberOfCells();
	int (*tri)[3] = new int[tri_N][3];
	for(i = 0; i < tri_N; i++){
		int dummy;
		int *dummy_tri;
		tri_data->GetNextCell(dummy, dummy_tri);
		tri[i][0] = dummy_tri[0];
		tri[i][1] = dummy_tri[1];
		tri[i][2] = dummy_tri[2];
	}

	d2d->Delete();
	poly_data->Delete();
	ps->Delete();

	int m = 0;
	for(i=0; i<tri_N; i++){
		face[deleted_face[m]][0] = star[tri[i][0]];
		face[deleted_face[m]][1] = star[tri[i][1]];
		face[deleted_face[m]][2] = star[tri[i][2]];

		v_link[star[tri[i][0]]]->append(star[tri[i][1]], deleted_face[m]);
		v_link[star[tri[i][0]]]->append(star[tri[i][2]], deleted_face[m]);

		v_link[star[tri[i][1]]]->append(star[tri[i][2]], deleted_face[m]);
		v_link[star[tri[i][1]]]->append(star[tri[i][0]], deleted_face[m]);

		v_link[star[tri[i][2]]]->append(star[tri[i][0]], deleted_face[m]);
		v_link[star[tri[i][2]]]->append(star[tri[i][1]], deleted_face[m]);

		m++;
	}
	delete[] deleted_face;
	delete[] id;
	//for(i=0; i<tri_N; i++)
		//delete tri[i];
	delete[] tri;
	delete[] z;

	for(i=0; i<degree; i++){
		mesh->isBound[star[i]] = 0;

		for(Node* current = v_link[star[i]]; current->next != NULL; current = current->next->next){
			int pairFace = -1;
			for(Node* search = v_link[star[i]]->next; search != NULL; search = search->next->next){
				if(current->v == search->v){
					pairFace = search->f;
					break;
				}
			}
			if(pairFace < 0){
				mesh->isBound[star[i]] = current->f;
			}
			else{
				int *face1 = face[current->f];
				if(star[i] == face1[0])
					f_link[current->f][0] = pairFace;
				else if(star[i] == face1[1])
					f_link[current->f][1] = pairFace;
				else 
					f_link[current->f][2] = pairFace;
			}
		}
	}
	delete[] star;
	delete v_link[index];
	v_link[index] = new Node;
	*/
	return true;
}

BOOL Simplifier::edgeCollapse(int index1, int index2)
{
	return true;
}

BOOL Simplifier::DK(int iter)
{
	int vertex_N = mesh->vertex_N;
	float (*vertex)[3] = mesh->vertex;
	float (*normal)[3] = mesh->normal;

	BOOL *mark = new BOOL[vertex_N];
	BOOL *isValid = new BOOL[vertex_N];
	float *cur = new float[vertex_N];
	float max_cur = -1.0f;
	float *area = new float[vertex_N];
	float max_area = -1.0f;
	
	for(int i=0; i<vertex_N; i++){
		area[i] = 0;
		cur[i] = 0;
		mark[i] = false;
		isValid[i] = true;

		int nei_N;
		int *nei;

		mesh->getConsistent1Ring(i, nei, nei_N);
		if(nei_N == 0){
			isValid[i] = false;
			mark[i] = true;	
			continue;
		}
		
		/*for(int j=0; j<nei_N; j++){
			float v1[3];
			MeshData::VEC(v1, vertex[i], vertex[nei[j]]);
			float l1 = MeshData::LENGTH(v1);
			if(l1 != 0)
				cur[i] += (float)fabs((v1[0]*normal[i][0] + v1[1]*normal[i][1] + v1[2]*normal[i][2])/l1);
		}*/
		if(mesh->isBound[i]){
			for(int j=0; j<nei_N-1; j++){
				int s = nei[j];
				int t = nei[j+1];
				area[i] += MeshData::AREA(vertex[i], vertex[s], vertex[t]);
			}
		}
		else{
			for(int j=0; j<nei_N; j++){
				int s = nei[j];
				int t = nei[(j+1)%nei_N];
				area[i] += MeshData::AREA(vertex[i], vertex[s], vertex[t]);
			}
		}
		if(area[i] > max_area)
			max_area = area[i];
		//if(cur[i] > max_cur)
			//max_cur = cur[i];
	}

	int* index = new int[vertex_N];
	float* w = new float[vertex_N];
	int valid_N = 0;
	for(i=0; i<vertex_N; i++){
		if(isValid[i]){
			index[valid_N] = i;
			w[valid_N] = area[i]; //cur[i]/max_cur;
			valid_N ++;
		}
	}
	delete[] cur;
	delete[] area;

	quickSort(index, w, 0, valid_N-1);
	delete[] w;

	BOOL update = false;
	for(i=0; i<valid_N; i++){
		if(!mark[index[i]]){
			int nei_N;
			int *nei;
			mesh->getConsistent1Ring(index[i], nei, nei_N);

			mark[index[i]] = true;
			
			BOOL flag = this->vertexRemoveCarefully(index[i]);
			update |= flag; 

			if(flag){
				mark[index[i]] = true;
				for(int j=0; j<nei_N; j++)
					mark[nei[j]] = true;
			}
		}
	}
	delete[] index;
	delete[] mark;

	return update;
}

void Simplifier::setMeshData(MeshData *mesh)
{
	this->mesh = mesh;
	int vertex_N = mesh->vertex_N;
	tag = new BOOL[vertex_N];
	for(int i=0; i<vertex_N; i++)
		tag[i] = false;
}

void Simplifier::quickSort(int *index, float *w, int start, int end)
{
	if(start < end){
		float v = w[end];
		int i = start-1;
		int j = end;
		while(j > i){
			for(i = i+1; w[i] < v; i++);
			for(j = j-1; w[j] > v; j--);
			float t = w[i]; 
			w[i] = w[j]; 
			w[j] = t;

			int tmp = index[i];
			index[i] = index[j];
			index[j] = tmp;
		}
		float t = w[j];
		w[j] = w[i];
		w[i] = w[end];
		w[end] = t;

		int tmp = index[j];
		index[j] = index[i];
		index[i] = index[end];
		index[end] = tmp;

		quickSort(index, w, start, i-1);
		quickSort(index, w, i+1, end);
	}
	else
		return;
}

void Simplifier::setRidgeAsTag()
{
	int vertex_N = mesh->vertex_N;
	BOOL* isRidge = mesh->isRidge;

	if(tag != NULL)
		delete[] tag;
	tag = new BOOL[vertex_N];
	for(int i=0; i<vertex_N; i++)
		if(isRidge[i])
			tag[i] = true;
		else
			tag[i] = false;
}

BOOL Simplifier::vertexRemoveCarefully(int index)
{
	/*
	Node** v_link = mesh->vertex_link;
	int (*f_link)[3] = mesh->face_link_E;
	float (*vertex)[3] = mesh->vertex;
	int (*face)[3] = mesh->face;
	BOOL *isBound = mesh->isBound;

	int degree;
	int *star;
	int *deleted_face;

	mesh->getConsistent1Ring(index, star, deleted_face, degree);

	//check degree
	if(degree < 3){
		delete[] star;
		delete[] deleted_face;
		return false;
	}

	//detecte tag edge 
	int tag_count = 0;
	int tag1 = -1;
	int tag2 = -1;
	for(int i=0; i<degree; i++)
		if(tag[star[i]]){
			tag_count++;
			if(tag1 < 0)
				tag1 = i;
			else
				tag2 = i;
		}
	if(tag_count == 1 && tag[index]){ 
		delete[] star;
		delete[] deleted_face;
		return false;
	}

	if(tag_count > 2 && tag2 - tag1 != tag_count - 1){
		delete[] star;
		delete[] deleted_face;
		return false;
	}

	BOOL isTagAdjacent = tag1+1 == tag2 || (tag1 == 0 && tag2 == degree-1);
	if(tag_count == 2 && !isTagAdjacent && !tag[index]){
		if(MeshData::DIST(vertex[index], vertex[star[tag1]]) 
			+ MeshData::DIST(vertex[index], vertex[star[tag2]]) > 1.0f ||
			fabs(MeshData::DOT(mesh->t_max[star[tag1]], mesh->t_max[star[tag2]])) < 0.5f ){
			delete[] star;
			delete[] deleted_face;
			return false;
		}
	}
	BOOL tag_edge = (tag_count == 2 && !isTagAdjacent);

	//check topology preservation
	for(i=0; i<degree; i++){
		int *nei, nei_N;
		mesh->getConsistent1Ring(star[i], nei, nei_N);
		int c = 0;
		for(int j=0; j<degree; j++)
			for(int k=0; k<nei_N; k++)
				if(star[j] == nei[k])
					c++;
		BOOL flag;
		if(isBound[index] != 0 && isBound[star[i]] != 0)
			flag = (c == 1);
		else
			flag = (c == 2);
		if(!flag){
			delete[] star;
			delete[] deleted_face;
			return false;
		}
	}
	
	//trianglization of the hole
	int tri_N;
	int (*tri)[3] = new int[degree-2][3];

	if(tag_edge){
		int id1_N = tag2-tag1+1;
		int* id1 = new int[id1_N];
		for(int i=0; i<id1_N; i++)
			id1[i] = tag1+i;

		int id2_N = degree-(tag2-tag1)+1;
		int *id2 = new int[id2_N]; 
		for(i=0; i<id2_N; i++)
			id2[i] = (tag2+i)%degree;

		float total_angle1 = 0;
		for(i=0; i<id1_N-1; i++){
			float v1[3], v2[3];
			MeshData::VEC(v1, vertex[index], vertex[star[id1[i]]]);
			MeshData::VEC(v2, vertex[index], vertex[star[id1[i+1]]]);
			float l1 = MeshData::LENGTH(v1);
			float l2 = MeshData::LENGTH(v2);
			if(l1 != 0 && l2 != 0){
				float arg = (float)acos(MeshData::DOT(v1, v2)/(l1*l2));
				if(arg != 0)
					total_angle1 += arg;
				else
					total_angle1 += 0.1f;
			}
			else
				total_angle1 += 0.1f;
		}

		float total_angle2 = 0;
		for(i=0; i<id2_N-1; i++){
			float v1[3], v2[3];
			MeshData::VEC(v1, vertex[index], vertex[star[id2[i]]]);
			MeshData::VEC(v2, vertex[index], vertex[star[id2[i+1]]]);
			float l1 = MeshData::LENGTH(v1);
			float l2 = MeshData::LENGTH(v2);
			if(l1 != 0 && l2 != 0){
				float arg = (float)acos(MeshData::DOT(v1, v2)/(l1*l2));
				if(arg != 0)
					total_angle2 += arg;
				else
					total_angle2 += 0.1f;
			}
			else
				total_angle2 += 0.1f;
		}

		float (*z)[2] = new float[degree][2];

		float S = 0;
		float a = PI/total_angle1;
		float v1[3], v2[3];
		float r1, r2, max_r = -1;

		MeshData::VEC(v1, vertex[index], vertex[star[id1[0]]]);
		r1 = MeshData::LENGTH(v1);
		for(i=0; i<id1_N-1; i++){
			if(max_r < r1)
				max_r = r1;

			if(r1 != 0){
				z[id1[i]][0] = (float)(r1*cos(S*a));
				z[id1[i]][1] = (float)(r1*sin(S*a));
			}
			else{
				z[id1[i]][0] = (float)(cos(S*a));
				z[id1[i]][1] = (float)(sin(S*a));
			}

			MeshData::VEC(v2, vertex[index], vertex[star[id1[i+1]]]);
			r2 = MeshData::LENGTH(v2);
			
			if(r1 != 0 && r2 != 0){
				float arg = (float)acos((v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2])/(r1*r2));
				if(arg != 0)
					S += arg;
				else
					S += 0.1f;
			}
			else
				S += 0.1f;

			v1[0] = v2[0];
			v1[1] = v2[1];
			v1[2] = v2[2];
			r1 = r2;
		}
		
		S = 0;
		a = PI/total_angle2;
		MeshData::VEC(v1, vertex[index], vertex[star[id2[0]]]);
		r1 = MeshData::LENGTH(v1);
		for(i=0; i<id2_N-1; i++){
			if(max_r < r1)
				max_r = r1;

			if(r1 != 0){
				z[id2[i]][0] = (float)(r1*cos(S*a + PI));
				z[id2[i]][1] = (float)(r1*sin(S*a + PI));
			}
			else{
				z[id2[i]][0] = (float)(cos(S*a + PI));
				z[id2[i]][1] = (float)(sin(S*a + PI));
			}

			MeshData::VEC(v2, vertex[index], vertex[star[id2[i+1]]]);
			r2 = MeshData::LENGTH(v2);
			
			if(r1 != 0 && r2 != 0){
				float arg = (float)acos((v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2])/(r1*r2));
				if(arg != 0)
					S += arg;
				else
					S += 0.1f;
			}
			else
				S += 0.1f;

			v1[0] = v2[0];
			v1[1] = v2[1];
			v1[2] = v2[2];
			r1 = r2;
		}
		if(max_r < r1)
				max_r = r1;
		
		vtkDelaunay2D* d2d = vtkDelaunay2D::New();
		d2d->SetTolerance(0.0000);
		d2d->SetOffset(500.0);
		d2d->SetAlpha(0);
		d2d->BoundingTriangulationOff();

		vtkPoints* ps = vtkPoints::New();
		vtkPolyData* poly_data = vtkPolyData::New();


		ps->Initialize();
		for(i=0; i<degree; i++)
			ps->InsertPoint(i, z[i][0], z[i][1], 0.0f);
		
		poly_data->Initialize();
		poly_data->SetPoints(ps);
		d2d->SetInput(poly_data);
		d2d->SetSource(poly_data);

		poly_data->Allocate(1);
		poly_data->InsertNextCell(VTK_POLYGON, id1_N, id1);
		d2d->Update();
		vtkCellArray* tri_data = d2d->GetOutput()->GetPolys();
		tri_N = tri_data->GetNumberOfCells();

		if(tri_N != id1_N-2){
			delete[] star;
			delete[] deleted_face;
			delete[] id1;
			delete[] id2;
			delete[] z;
			delete[] tri;
			d2d->Delete();
			poly_data->Delete();
			ps->Delete();
			return false;
		}

		for(i = 0; i < tri_N; i++){
			int dummy;
			int *dummy_tri;
			tri_data->GetNextCell(dummy, dummy_tri);
			tri[i][0] = dummy_tri[0];
			tri[i][1] = dummy_tri[1];
			tri[i][2] = dummy_tri[2];
		}

		poly_data->Initialize();
		poly_data->SetPoints(ps);
		poly_data->Allocate(1);
		poly_data->InsertNextCell(VTK_POLYGON, id2_N, id2);
		d2d->Update();
		tri_data = d2d->GetOutput()->GetPolys();
		tri_data->InitTraversal();
		tri_N += tri_data->GetNumberOfCells();
		
		if(tri_N != degree-2){
			delete[] star;
			delete[] deleted_face;
			delete[] id1;
			delete[] id2;
			delete[] z;
			delete[] tri;
			d2d->Delete();
			poly_data->Delete();
			ps->Delete();
			return false;
		}

		for(; i < tri_N; i++){
			int dummy;
			int *dummy_tri;
			tri_data->GetNextCell(dummy, dummy_tri);
			tri[i][0] = dummy_tri[0];
			tri[i][1] = dummy_tri[1];
			tri[i][2] = dummy_tri[2];
		}

		d2d->Delete();
		poly_data->Delete();
		ps->Delete();
		delete[] z;
		delete[] id1;
		delete[] id2;
	}
	//normal case
	else{
		float total_angle = 0;
		if(mesh->isBound[index])
			for(i=0; i<degree-1; i++){
				float v1[3], v2[3];
				MeshData::VEC(v1, vertex[index], vertex[star[i]]);
				MeshData::VEC(v2, vertex[index], vertex[star[i+1]]);
				float l1 = MeshData::LENGTH(v1);
				float l2 = MeshData::LENGTH(v2);
				if(l1 != 0 && l2 != 0)
					total_angle += (float)acos(MeshData::DOT(v1, v2)/(l1*l2));
		}
		else
			for(i=0; i<degree; i++){
				float v1[3], v2[3];
				MeshData::VEC(v1, vertex[index], vertex[star[i]]);
				MeshData::VEC(v2, vertex[index], vertex[star[(i+1)%degree]]);
				float l1 = MeshData::LENGTH(v1);
				float l2 = MeshData::LENGTH(v2);
				if(l1 != 0 && l2 != 0){
					float arg = (float)acos(MeshData::DOT(v1, v2)/(l1*l2));
					if(arg != 0)
						total_angle += arg;
					else
						total_angle += 0.1f;
				}
				else
					total_angle += 0.1f;
			}

		float (*z)[2] = new float[degree][2];
		float S = 0;
		float a;
		if(!mesh->isBound[index])
			a = 2.0*PI/total_angle;
		else{
			a = PI/total_angle;
		}
		float v1[3], v2[3];
		float r1, r2, max_r = -1;
		MeshData::VEC(v1, vertex[index], vertex[star[0]]);
		r1 = MeshData::LENGTH(v1);
		for(i=0; i<degree-1; i++){
			if(max_r < r1)
				max_r = r1;

			if(r1 != 0){
				z[i][0] = (float)(r1*cos(S*a));
				z[i][1] = (float)(r1*sin(S*a));
			}
			else{
				z[i][0] = (float)(cos(S*a));
				z[i][1] = (float)(sin(S*a));
			}

			MeshData::VEC(v2, vertex[index], vertex[star[i+1]]);
			r2 = MeshData::LENGTH(v2);
			
			if(r1 != 0 && r2 != 0){
				float arg = (float)acos((v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2])/(r1*r2));
				if(arg != 0)
					S += arg;
				else
					S += 0.1f;
			}
			else
				S += 0.1f;

			v1[0] = v2[0];
			v1[1] = v2[1];
			v1[2] = v2[2];
			r1 = r2;
		}
		if(max_r < r1)
				max_r = r1;

		if(r1 != 0){
			z[i][0] = (float)(r1*cos(S*a));
			z[i][1] = (float)(r1*sin(S*a));
		}
		else{
			z[i][0] = (float)(cos(S*a));
			z[i][1] = (float)(sin(S*a));
		}

		for(i=0; i<degree; i++){
			z[i][0] = z[i][0]/max_r;
			z[i][1] = z[i][1]/max_r;
		}
		
		vtkDelaunay2D* d2d = vtkDelaunay2D::New();
		d2d->SetTolerance(0.0000);
		d2d->SetOffset(500.0);
		d2d->SetAlpha(0);
		d2d->BoundingTriangulationOff();

		vtkPoints* ps = vtkPoints::New();
		vtkPolyData* poly_data = vtkPolyData::New();


		ps->Initialize();
		int *id = new int[degree];
		for(i=0; i<degree; i++){
			id[i] = i;
			ps->InsertPoint(i, z[i][0], z[i][1], 0.0f);
		}
		poly_data->Initialize();
		poly_data->SetPoints(ps);
		d2d->SetInput(poly_data);
		d2d->SetSource(poly_data);

		poly_data->Allocate(1);
		poly_data->InsertNextCell(VTK_POLYGON, degree, id); 

		d2d->Update();
		vtkCellArray* tri_data = d2d->GetOutput()->GetPolys();

		tri_N = tri_data->GetNumberOfCells();

		if(tri_N != degree-2){
			delete[] star;
			delete[] deleted_face;
			delete[] id;
			delete[] z;
			delete[] tri;
			d2d->Delete();
			poly_data->Delete();
			ps->Delete();
			return false;
		}

		for(i = 0; i < tri_N; i++){
			int dummy;
			int *dummy_tri;
			tri_data->GetNextCell(dummy, dummy_tri);
			tri[i][0] = dummy_tri[0];
			tri[i][1] = dummy_tri[1];
			tri[i][2] = dummy_tri[2];
		}

		d2d->Delete();
		poly_data->Delete();
		ps->Delete();
		delete[] id;
		delete[] z;
	}

	for(i=0; i<tri_N; i++){
		if(tri[i][0] > tri[i][1]){
			int tmp = tri[i][0];
			tri[i][0] = tri[i][1];
			tri[i][1] = tmp;
		}
		if(tri[i][1] > tri[i][2]){
			int tmp = tri[i][1];
			tri[i][1] = tri[i][2];
			tri[i][2] = tmp;
		}
		if(tri[i][0] > tri[i][1]){
			int tmp = tri[i][0];
			tri[i][0] = tri[i][1];
			tri[i][1] = tmp;
		}
	}

	if(isBound[index])
		for(i=0; i<degree-1; i++){
			mesh->deleteFace(deleted_face[i]);
		}
	else
		for(i=0; i<degree; i++){
			mesh->deleteFace(deleted_face[i]);
		}

	for(i=0; i<degree; i++){
		Node *tmp_adjacent = new Node;
		tmp_adjacent->next = NULL;

		for(Node* current = v_link[star[i]]; current->next != NULL; current = current->next->next){
			if(current->v != index && current->next->v != index){
				tmp_adjacent->append(current->v, current->f);
				tmp_adjacent->append(current->next->v, current->next->f);
			}
		}
		delete v_link[star[i]];
		v_link[star[i]] = tmp_adjacent;
	}

	int m = 0;
	for(i=0; i<tri_N; i++){
		face[deleted_face[m]][0] = star[tri[i][0]];
		face[deleted_face[m]][1] = star[tri[i][1]];
		face[deleted_face[m]][2] = star[tri[i][2]];

		v_link[star[tri[i][0]]]->append(star[tri[i][1]], deleted_face[m]);
		v_link[star[tri[i][0]]]->append(star[tri[i][2]], deleted_face[m]);

		v_link[star[tri[i][1]]]->append(star[tri[i][2]], deleted_face[m]);
		v_link[star[tri[i][1]]]->append(star[tri[i][0]], deleted_face[m]);

		v_link[star[tri[i][2]]]->append(star[tri[i][0]], deleted_face[m]);
		v_link[star[tri[i][2]]]->append(star[tri[i][1]], deleted_face[m]);

		m++;
	}
	delete[] deleted_face;
	delete[] tri;

	for(i=0; i<degree; i++){
		mesh->isBound[star[i]] = 0;

		for(Node* current = v_link[star[i]]; current->next != NULL; current = current->next->next){
			int pairFace = -1;
			for(Node* search = v_link[star[i]]->next; search != NULL; search = search->next->next){
				if(current->v == search->v){
					pairFace = search->f;
					break;
				}
			}
			if(pairFace < 0){
				mesh->isBound[star[i]] = current->f;
			}
			else{
				int *face1 = face[current->f];
				if(star[i] == face1[0])
					f_link[current->f][0] = pairFace;
				else if(star[i] == face1[1])
					f_link[current->f][1] = pairFace;
				else 
					f_link[current->f][2] = pairFace;
			}
		}
	}
	delete[] star;
	delete v_link[index];
	v_link[index] = new Node;
	*/
	return true;
}

void Simplifier::setRavineAsTag()
{
	int vertex_N = mesh->vertex_N;
	BOOL* isRavine = mesh->isRavine;
	if(tag != NULL)
		delete[] tag;
	tag = new BOOL[vertex_N];
	for(int i=0; i<vertex_N; i++)
		if(isRavine[i])
			tag[i] = true;
		else
			tag[i] = false;
}

void Simplifier::setTagEdges()
{
	int vertex_N = mesh->vertex_N;
	BOOL *isRidge = mesh->isRidge;
	BOOL *isRavine = mesh->isRavine;

	Node** ridge_edge = mesh->ridge_edge;
	Node** ravine_edge = mesh->ravine_edge;

	if(ridge_edge != NULL){
		for(int i=0; i<vertex_N; i++)
			delete ridge_edge[i];
		delete[] ridge_edge;
	}
	
	if(ravine_edge != NULL){
		for(int i=0; i<vertex_N; i++)
			delete ravine_edge[i];
		delete[] ravine_edge;
	}

	if(both_edge != NULL){
		for(int i=0; i<vertex_N; i++)
			delete both_edge[i];
		delete[] both_edge;
	}

	ridge_edge = mesh->ridge_edge = new Node*[vertex_N];
	ravine_edge = mesh->ravine_edge = new Node*[vertex_N];
	both_edge = new Node*[vertex_N];
	for(int i=0; i<vertex_N; i++){
		ridge_edge[i] = new Node;
		ravine_edge[i] = new Node;
		both_edge[i] = new Node;
	}

	for(i=0; i<vertex_N; i++){
		if(!isRidge[i] && !isRavine[i])
			continue;
		int *nei, nei_N;
		mesh->getConsistent1Ring(i, nei, nei_N);
		
		if(isRidge[i]){
			for(int j=0; j<nei_N; j++)
				if(isRidge[nei[j]]){
					ridge_edge[i]->append(nei[j], -1);
					both_edge[i]->append(nei[j], -1);
				}
		}

		if(isRavine[i]){
			for(int j=0; j<nei_N; j++)
				if(isRavine[nei[j]]){
					ravine_edge[i]->append(nei[j], -1);
					if(!isRidge[i] || !isRidge[nei[j]])
						both_edge[i]->append(nei[j], -1);
				}
		}
	}
}

BOOL Simplifier::vertexRemoveTagEdge(int index, float T, float A)
{
	/*
	Node** v_link = mesh->vertex_link;
	int (*f_link)[3] = mesh->face_link_E;
	float (*vertex)[3] = mesh->vertex;
	int (*face)[3] = mesh->face;
	BOOL *isBound = mesh->isBound;

	Node** ridge_edge = mesh->ridge_edge;
	Node** ravine_edge = mesh->ravine_edge;

	int degree;
	int *star;
	int *deleted_face;

	mesh->getConsistent1Ring(index, star, deleted_face, degree);

	//check degree
	if(degree < 3){
		delete[] star;
		delete[] deleted_face;
		return false;
	}

	BOOL tag_edge = false;
	int tag1, tag2;
	BOOL is_ridge_tag = false;
	if(ridge_edge[index]->next == NULL && ravine_edge[index]->next == NULL){
		////////////////////////comment start
		int ridge_count = 0;
		int ridge_tag1 = -1;
		int ridge_tag2 = -1;
		for(int i=0; i<degree; i++){
			if(ridge_edge[star[i]]->next != NULL){
				ridge_count++;
				if(ridge_tag1 < 0)
					ridge_tag1 = i;
				else
					ridge_tag2 = i;
			}
		}

		int ravine_count = 0;
		int ravine_tag1 = -1;
		int ravine_tag2 = -1;
		for(i=0; i<degree; i++){
			if(ravine_edge[star[i]]->next != NULL){
				ravine_count++;
				if(ravine_tag1 < 0)
					ravine_tag1 = i;
				else
					ravine_tag2 = i;
			}
		}

		if(ravine_count == 0 && ridge_count == 2){
			tag_edge = true;
			is_ridge_tag = true;
			tag1 = ridge_tag1;
			tag2 = ridge_tag2;
		}
		else if(ridge_count == 0 && ravine_count == 2){
			tag_edge = true;
			tag1 = ravine_tag1;
			tag2 = ravine_tag2;
		}
		
		if(tag_edge){
			if(MeshData::DIST(vertex[index], vertex[star[tag1]]) 
				+ MeshData::DIST(vertex[index], vertex[star[tag2]]) > T){
				tag_edge = false;
			}
			else{
				float v1[3], v2[3];
				MeshData::VEC(v1, vertex[index], vertex[star[tag1]]);
				MeshData::VEC(v2, vertex[index], vertex[star[tag2]]);
				if(MeshData::DOT(v1, v2) > A*MeshData::LENGTH(v1)*MeshData::LENGTH(v2))
					tag_edge = false;
			}
		}/////////////////////////////comment end
	}
	else{
		int ridge_count = 0;
		int ridge_tag1 = -1;
		int ridge_tag2 = -1;
		for(Node* current = ridge_edge[index]; current->next!=NULL; current=current->next){
			ridge_count++;
			if(ridge_tag1 < 0)
				ridge_tag1 = current->v;
			else
				ridge_tag2 = current->v;
		}

		int ravine_count = 0;
		int ravine_tag1 = -1;
		int ravine_tag2 = -1;
		for(current = ravine_edge[index]; current->next!=NULL; current=current->next){
			ravine_count++;
			if(ravine_tag1 < 0)
				ravine_tag1 = current->v;
			else
				ravine_tag2 = current->v;
		}

		if(ravine_count == 0){
			if(ridge_count == 2){
				tag_edge = true;
				is_ridge_tag = true;
				tag1 = ridge_tag1;
				tag2 = ridge_tag2;
			}
			if(ridge_count == 1){
				int c = 0;
				int i1;
				for(int i=0; i<degree; i++){
					if(ridge_edge[star[i]]->next != NULL){
						if(star[i] != ridge_tag1){
							c++;
							i1 = star[i];
						}
					}
				}
				if(c == 1){
					if(MeshData::DIST(vertex[index], vertex[i1]) < 0.5*T){
						tag_edge = true;
						is_ridge_tag = true;
						tag1 = ridge_tag1;
						tag2 = i1;
					}
				}
			}
		}
		else if(ridge_count == 0){
			if(ravine_count == 2){
				tag_edge = true;
				tag1 = ravine_tag1;
				tag2 = ravine_tag2;
			}
			if(ravine_count == 1){
				int c = 0;
				int i1;
				for(int i=0; i<degree; i++){
					if(ravine_edge[star[i]]->next != NULL){
						if(star[i] != ravine_tag1){
							c++;
							i1 = star[i];
						}
					}
				}
				if(c == 1){
					if(MeshData::DIST(vertex[index], vertex[i1]) < 0.5*T){
						tag_edge = true;
						tag1 = ravine_tag1;
						tag2 = i1;
					}
				}
			}
		}
		/////////////////comment start
		if(tag_edge){
			float v1[3], v2[3];
			MeshData::VEC(v1, vertex[index], vertex[tag1]);
			MeshData::VEC(v2, vertex[index], vertex[tag2]);
			if(MeshData::DOT(v1, v2) > A*MeshData::LENGTH(v1)*MeshData::LENGTH(v2))
				tag_edge = false;
		}//////////////////comment end
		
		if(!tag_edge){
			delete[] star;
			delete[] deleted_face;
			return false;
		}

		for(int i=0; i<degree; i++){
			if(star[i] == tag1)
				tag1 = i;
			else if(star[i] == tag2)
				tag2 = i;
		}
		if(tag2 < tag1){
			int tmp = tag1;
			tag1 = tag2;
			tag2 = tmp;
		}
	}

	//check topology preservation
	for(int i=0; i<degree; i++){
		int *nei, nei_N;
		mesh->getConsistent1Ring(star[i], nei, nei_N);
		int c = 0;
		for(int j=0; j<degree; j++)
			for(int k=0; k<nei_N; k++)
				if(star[j] == nei[k])
					c++;
		BOOL flag;
		if(isBound[index] != 0 && isBound[star[i]] != 0)
			flag = (c == 1);
		else
			flag = (c == 2);
		if(!flag){
			delete[] star;
			delete[] deleted_face;
			return false;
		}
	}
	
	//trianglization of the hole
	int tri_N;
	int (*tri)[3] = new int[degree-2][3];

	if(tag_edge && !(tag2 - tag1 == 1 || (tag2+1)%degree == tag1)){
		int id1_N = tag2-tag1+1;
		int* id1 = new int[id1_N];
		for(int i=0; i<id1_N; i++)
			id1[i] = tag1+i;

		int id2_N = degree-(tag2-tag1)+1;
		int *id2 = new int[id2_N]; 
		for(i=0; i<id2_N; i++)
			id2[i] = (tag2+i)%degree;

		float total_angle1 = 0;
		for(i=0; i<id1_N-1; i++){
			float v1[3], v2[3];
			MeshData::VEC(v1, vertex[index], vertex[star[id1[i]]]);
			MeshData::VEC(v2, vertex[index], vertex[star[id1[i+1]]]);
			float l1 = MeshData::LENGTH(v1);
			float l2 = MeshData::LENGTH(v2);
			if(l1 != 0 && l2 != 0){
				float arg = (float)acos(MeshData::DOT(v1, v2)/(l1*l2));
				if(arg != 0)
					total_angle1 += arg;
				else
					total_angle1 += 0.1f;
			}
			else
				total_angle1 += 0.1f;
		}

		float total_angle2 = 0;
		for(i=0; i<id2_N-1; i++){
			float v1[3], v2[3];
			MeshData::VEC(v1, vertex[index], vertex[star[id2[i]]]);
			MeshData::VEC(v2, vertex[index], vertex[star[id2[i+1]]]);
			float l1 = MeshData::LENGTH(v1);
			float l2 = MeshData::LENGTH(v2);
			if(l1 != 0 && l2 != 0){
				float arg = (float)acos(MeshData::DOT(v1, v2)/(l1*l2));
				if(arg != 0)
					total_angle2 += arg;
				else
					total_angle2 += 0.1f;
			}
			else
				total_angle2 += 0.1f;
		}

		float (*z)[2] = new float[degree][2];

		float S = 0;
		float a = PI/total_angle1;
		float v1[3], v2[3];
		float r1, r2, max_r = -1;

		MeshData::VEC(v1, vertex[index], vertex[star[id1[0]]]);
		r1 = MeshData::LENGTH(v1);
		for(i=0; i<id1_N-1; i++){
			if(max_r < r1)
				max_r = r1;

			if(r1 != 0){
				z[id1[i]][0] = (float)(r1*cos(S*a));
				z[id1[i]][1] = (float)(r1*sin(S*a));
			}
			else{
				z[id1[i]][0] = (float)(cos(S*a));
				z[id1[i]][1] = (float)(sin(S*a));
			}

			MeshData::VEC(v2, vertex[index], vertex[star[id1[i+1]]]);
			r2 = MeshData::LENGTH(v2);
			
			if(r1 != 0 && r2 != 0){
				float arg = (float)acos((v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2])/(r1*r2));
				if(arg != 0)
					S += arg;
				else
					S += 0.1f;
			}
			else
				S += 0.1f;

			v1[0] = v2[0];
			v1[1] = v2[1];
			v1[2] = v2[2];
			r1 = r2;
		}
		
		S = 0;
		a = PI/total_angle2;
		MeshData::VEC(v1, vertex[index], vertex[star[id2[0]]]);
		r1 = MeshData::LENGTH(v1);
		for(i=0; i<id2_N-1; i++){
			if(max_r < r1)
				max_r = r1;

			if(r1 != 0){
				z[id2[i]][0] = (float)(r1*cos(S*a + PI));
				z[id2[i]][1] = (float)(r1*sin(S*a + PI));
			}
			else{
				z[id2[i]][0] = (float)(cos(S*a + PI));
				z[id2[i]][1] = (float)(sin(S*a + PI));
			}

			MeshData::VEC(v2, vertex[index], vertex[star[id2[i+1]]]);
			r2 = MeshData::LENGTH(v2);
			
			if(r1 != 0 && r2 != 0){
				float arg = (float)acos((v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2])/(r1*r2));
				if(arg != 0)
					S += arg;
				else
					S += 0.1f;
			}
			else
				S += 0.1f;

			v1[0] = v2[0];
			v1[1] = v2[1];
			v1[2] = v2[2];
			r1 = r2;
		}
		if(max_r < r1)
				max_r = r1;
		
		vtkDelaunay2D* d2d = vtkDelaunay2D::New();
		d2d->SetTolerance(0.0000);
		d2d->SetOffset(500.0);
		d2d->SetAlpha(0);
		d2d->BoundingTriangulationOff();

		vtkPoints* ps = vtkPoints::New();
		vtkPolyData* poly_data = vtkPolyData::New();


		ps->Initialize();
		for(i=0; i<degree; i++)
			ps->InsertPoint(i, z[i][0], z[i][1], 0.0f);
		
		poly_data->Initialize();
		poly_data->SetPoints(ps);
		d2d->SetInput(poly_data);
		d2d->SetSource(poly_data);

		poly_data->Allocate(1);
		poly_data->InsertNextCell(VTK_POLYGON, id1_N, id1);
		d2d->Update();
		vtkCellArray* tri_data = d2d->GetOutput()->GetPolys();
		tri_N = tri_data->GetNumberOfCells();

		if(tri_N != id1_N-2){
			delete[] star;
			delete[] deleted_face;
			delete[] id1;
			delete[] id2;
			delete[] z;
			delete[] tri;
			d2d->Delete();
			poly_data->Delete();
			ps->Delete();
			return false;
		}

		for(i = 0; i < tri_N; i++){
			int dummy;
			int *dummy_tri;
			tri_data->GetNextCell(dummy, dummy_tri);
			tri[i][0] = dummy_tri[0];
			tri[i][1] = dummy_tri[1];
			tri[i][2] = dummy_tri[2];
		}

		poly_data->Initialize();
		poly_data->SetPoints(ps);
		poly_data->Allocate(1);
		poly_data->InsertNextCell(VTK_POLYGON, id2_N, id2);
		d2d->Update();
		tri_data = d2d->GetOutput()->GetPolys();
		tri_data->InitTraversal();
		tri_N += tri_data->GetNumberOfCells();
		
		if(tri_N != degree-2){
			delete[] star;
			delete[] deleted_face;
			delete[] id1;
			delete[] id2;
			delete[] z;
			delete[] tri;
			d2d->Delete();
			poly_data->Delete();
			ps->Delete();
			return false;
		}

		for(; i < tri_N; i++){
			int dummy;
			int *dummy_tri;
			tri_data->GetNextCell(dummy, dummy_tri);
			tri[i][0] = dummy_tri[0];
			tri[i][1] = dummy_tri[1];
			tri[i][2] = dummy_tri[2];
		}

		d2d->Delete();
		poly_data->Delete();
		ps->Delete();
		delete[] z;
		delete[] id1;
		delete[] id2;
	}
	//normal case
	else{
		float total_angle = 0;
		if(mesh->isBound[index])
			for(i=0; i<degree-1; i++){
				float v1[3], v2[3];
				MeshData::VEC(v1, vertex[index], vertex[star[i]]);
				MeshData::VEC(v2, vertex[index], vertex[star[i+1]]);
				float l1 = MeshData::LENGTH(v1);
				float l2 = MeshData::LENGTH(v2);
				if(l1 != 0 && l2 != 0)
					total_angle += (float)acos(MeshData::DOT(v1, v2)/(l1*l2));
		}
		else
			for(i=0; i<degree; i++){
				float v1[3], v2[3];
				MeshData::VEC(v1, vertex[index], vertex[star[i]]);
				MeshData::VEC(v2, vertex[index], vertex[star[(i+1)%degree]]);
				float l1 = MeshData::LENGTH(v1);
				float l2 = MeshData::LENGTH(v2);
				if(l1 != 0 && l2 != 0){
					float arg = (float)acos(MeshData::DOT(v1, v2)/(l1*l2));
					if(arg != 0)
						total_angle += arg;
					else
						total_angle += 0.1f;
				}
				else
					total_angle += 0.1f;
			}

		float (*z)[2] = new float[degree][2];
		float S = 0;
		float a;
		if(!mesh->isBound[index])
			a = 2.0*PI/total_angle;
		else{
			a = PI/total_angle;
		}
		float v1[3], v2[3];
		float r1, r2, max_r = -1;
		MeshData::VEC(v1, vertex[index], vertex[star[0]]);
		r1 = MeshData::LENGTH(v1);
		for(i=0; i<degree-1; i++){
			if(max_r < r1)
				max_r = r1;

			if(r1 != 0){
				z[i][0] = (float)(r1*cos(S*a));
				z[i][1] = (float)(r1*sin(S*a));
			}
			else{
				z[i][0] = (float)(cos(S*a));
				z[i][1] = (float)(sin(S*a));
			}

			MeshData::VEC(v2, vertex[index], vertex[star[i+1]]);
			r2 = MeshData::LENGTH(v2);
			
			if(r1 != 0 && r2 != 0){
				float arg = (float)acos((v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2])/(r1*r2));
				if(arg != 0)
					S += arg;
				else
					S += 0.1f;
			}
			else
				S += 0.1f;

			v1[0] = v2[0];
			v1[1] = v2[1];
			v1[2] = v2[2];
			r1 = r2;
		}
		if(max_r < r1)
				max_r = r1;

		if(r1 != 0){
			z[i][0] = (float)(r1*cos(S*a));
			z[i][1] = (float)(r1*sin(S*a));
		}
		else{
			z[i][0] = (float)(cos(S*a));
			z[i][1] = (float)(sin(S*a));
		}

		for(i=0; i<degree; i++){
			z[i][0] = z[i][0]/max_r;
			z[i][1] = z[i][1]/max_r;
		}
		
		vtkDelaunay2D* d2d = vtkDelaunay2D::New();
		d2d->SetTolerance(0.0000);
		d2d->SetOffset(500.0);
		d2d->SetAlpha(0);
		d2d->BoundingTriangulationOff();

		vtkPoints* ps = vtkPoints::New();
		vtkPolyData* poly_data = vtkPolyData::New();


		ps->Initialize();
		int *id = new int[degree];
		for(i=0; i<degree; i++){
			id[i] = i;
			ps->InsertPoint(i, z[i][0], z[i][1], 0.0f);
		}
		poly_data->Initialize();
		poly_data->SetPoints(ps);
		d2d->SetInput(poly_data);
		d2d->SetSource(poly_data);

		poly_data->Allocate(1);
		poly_data->InsertNextCell(VTK_POLYGON, degree, id); 

		d2d->Update();
		vtkCellArray* tri_data = d2d->GetOutput()->GetPolys();

		tri_N = tri_data->GetNumberOfCells();

		if(tri_N != degree-2){
			delete[] star;
			delete[] deleted_face;
			delete[] id;
			delete[] z;
			delete[] tri;
			d2d->Delete();
			poly_data->Delete();
			ps->Delete();
			return false;
		}

		for(i = 0; i < tri_N; i++){
			int dummy;
			int *dummy_tri;
			tri_data->GetNextCell(dummy, dummy_tri);
			tri[i][0] = dummy_tri[0];
			tri[i][1] = dummy_tri[1];
			tri[i][2] = dummy_tri[2];
		}

		d2d->Delete();
		poly_data->Delete();
		ps->Delete();
		delete[] id;
		delete[] z;
	}

	for(i=0; i<tri_N; i++){
		if(tri[i][0] > tri[i][1]){
			int tmp = tri[i][0];
			tri[i][0] = tri[i][1];
			tri[i][1] = tmp;
		}
		if(tri[i][1] > tri[i][2]){
			int tmp = tri[i][1];
			tri[i][1] = tri[i][2];
			tri[i][2] = tmp;
		}
		if(tri[i][0] > tri[i][1]){
			int tmp = tri[i][0];
			tri[i][0] = tri[i][1];
			tri[i][1] = tmp;
		}
	}

	if(isBound[index])
		for(i=0; i<degree-1; i++){
			mesh->deleteFace(deleted_face[i]);
		}
	else
		for(i=0; i<degree; i++){
			mesh->deleteFace(deleted_face[i]);
		}

	for(i=0; i<degree; i++){
		Node *tmp_adjacent = new Node;
		tmp_adjacent->next = NULL;

		for(Node* current = v_link[star[i]]; current->next != NULL; current = current->next->next){
			if(current->v != index && current->next->v != index){
				tmp_adjacent->append(current->v, current->f);
				tmp_adjacent->append(current->next->v, current->next->f);
			}
		}
		delete v_link[star[i]];
		v_link[star[i]] = tmp_adjacent;
	}

	int m = 0;
	for(i=0; i<tri_N; i++){
		face[deleted_face[m]][0] = star[tri[i][0]];
		face[deleted_face[m]][1] = star[tri[i][1]];
		face[deleted_face[m]][2] = star[tri[i][2]];

		v_link[star[tri[i][0]]]->append(star[tri[i][1]], deleted_face[m]);
		v_link[star[tri[i][0]]]->append(star[tri[i][2]], deleted_face[m]);

		v_link[star[tri[i][1]]]->append(star[tri[i][2]], deleted_face[m]);
		v_link[star[tri[i][1]]]->append(star[tri[i][0]], deleted_face[m]);

		v_link[star[tri[i][2]]]->append(star[tri[i][0]], deleted_face[m]);
		v_link[star[tri[i][2]]]->append(star[tri[i][1]], deleted_face[m]);

		m++;
	}
	delete[] deleted_face;
	delete[] tri;

	for(i=0; i<degree; i++){
		mesh->isBound[star[i]] = 0;

		for(Node* current = v_link[star[i]]; current->next != NULL; current = current->next->next){
			int pairFace = -1;
			for(Node* search = v_link[star[i]]->next; search != NULL; search = search->next->next){
				if(current->v == search->v){
					pairFace = search->f;
					break;
				}
			}
			if(pairFace < 0){
				mesh->isBound[star[i]] = current->f;
			}
			else{
				int *face1 = face[current->f];
				if(star[i] == face1[0])
					f_link[current->f][0] = pairFace;
				else if(star[i] == face1[1])
					f_link[current->f][1] = pairFace;
				else 
					f_link[current->f][2] = pairFace;
			}
		}
	}

	//update tag edge information
	if(tag_edge){
		if(is_ridge_tag){
			Node *tmp = new Node;
			BOOL detect_pair = false;
			for(Node* current = ridge_edge[star[tag1]]; current->next != NULL; current=current->next){
				if(current->v != index)
					tmp->append(current->v, -1);
				if(current->v == star[tag2])
					detect_pair = true;
			}
			if(!detect_pair)
				tmp->append(star[tag2], -1);
			delete ridge_edge[star[tag1]];
			ridge_edge[star[tag1]] = tmp;

			tmp = new Node;
			detect_pair = false;
			for(current = ridge_edge[star[tag2]]; current->next != NULL; current=current->next){
				if(current->v != index)
					tmp->append(current->v, -1);
				if(current->v == star[tag1])
					detect_pair = true;
			}
			if(!detect_pair)
				tmp->append(star[tag1], -1);
			delete ridge_edge[star[tag2]];
			ridge_edge[star[tag2]] = tmp;

			delete ridge_edge[index];
			ridge_edge[index] = new Node;
		}
		else{
			Node *tmp = new Node;
			BOOL detect_pair = false;
			for(Node* current = ravine_edge[star[tag1]]; current->next != NULL; current=current->next){
				if(current->v != index)
					tmp->append(current->v, -1);
				if(current->v == star[tag2])
					detect_pair = true;
			}
			if(!detect_pair)
				tmp->append(star[tag2], -1);
			delete ravine_edge[star[tag1]];
			ravine_edge[star[tag1]] = tmp;

			tmp = new Node;
			detect_pair = false;
			for(current = ravine_edge[star[tag2]]; current->next != NULL; current=current->next){
				if(current->v != index)
					tmp->append(current->v, -1);
				if(current->v == star[tag1])
					detect_pair = true;
			}
			if(!detect_pair)
				tmp->append(star[tag1], -1);
			delete ravine_edge[star[tag2]];
			ravine_edge[star[tag2]] = tmp;

			delete ravine_edge[index];
			ravine_edge[index] = new Node;
		}
		//if(tag2 - tag1 == 1 || (tag2+1)%degree == tag1)
			//tag_edge = false;	
	}

	delete[] star;
	delete v_link[index];
	v_link[index] = new Node;
	*/
	return true;
}

BOOL Simplifier::DK2(float T, float A)
{
	int vertex_N = mesh->vertex_N;
	float (*vertex)[3] = mesh->vertex;
	float (*normal)[3] = mesh->normal;

	BOOL *mark = new BOOL[vertex_N];
	BOOL *isValid = new BOOL[vertex_N];
	float *cur = new float[vertex_N];
	float max_cur = -1.0f;
	float *area = new float[vertex_N];
	float max_area = -1.0f;
	
	for(int i=0; i<vertex_N; i++){
		area[i] = 0;
		cur[i] = 0;
		mark[i] = false;
		isValid[i] = true;

		int nei_N;
		int *nei;

		mesh->getConsistent1Ring(i, nei, nei_N);
		if(nei_N == 0){
			isValid[i] = false;
			mark[i] = true;	
			continue;
		}
		
		/*for(int j=0; j<nei_N; j++){
			float v1[3];
			MeshData::VEC(v1, vertex[i], vertex[nei[j]]);
			float l1 = MeshData::LENGTH(v1);
			if(l1 != 0)
				cur[i] += (float)fabs((v1[0]*normal[i][0] + v1[1]*normal[i][1] + v1[2]*normal[i][2])/l1);
		}*/
		if(mesh->isBound[i]){
			for(int j=0; j<nei_N-1; j++){
				int s = nei[j];
				int t = nei[j+1];
				area[i] += MeshData::AREA(vertex[i], vertex[s], vertex[t]);
			}
		}
		else{
			for(int j=0; j<nei_N; j++){
				int s = nei[j];
				int t = nei[(j+1)%nei_N];
				area[i] += MeshData::AREA(vertex[i], vertex[s], vertex[t]);
			}
		}

		if(area[i] > max_area)
			max_area = area[i];
		//if(cur[i] > max_cur)
			//max_cur = cur[i];
	}

	int* index = new int[vertex_N];
	float* w = new float[vertex_N];
	int valid_N = 0;
	for(i=0; i<vertex_N; i++){
		if(isValid[i]){
			index[valid_N] = i;
			w[valid_N] = area[i]; //cur[i]/max_cur;
			valid_N ++;
		}
	}
	delete[] cur;
	delete[] area;

	quickSort(index, w, 0, valid_N-1);
	delete[] w;

	BOOL update = false;
	for(i=0; i<valid_N; i++){
		if(!mark[index[i]]){
			int nei_N;
			int *nei;
			mesh->getConsistent1Ring(index[i], nei, nei_N);

			mark[index[i]] = true;
			
			BOOL flag = this->vertexRemoveTagEdge(index[i], T, A);
			update |= flag; 

			if(flag){
				mark[index[i]] = true;
				for(int j=0; j<nei_N; j++)
					mark[nei[j]] = true;
			}
		}
	}
	delete[] index;
	delete[] mark;

	return update;
}

void Simplifier::setTagEdgesAngle(float T)
{
	int vertex_N = mesh->vertex_N;
	int face_N = mesh->face_N;
	int (*face)[3] = mesh->face;
	float (*vertex)[3] = mesh->vertex;
	float (*normal)[3] = mesh->normal_f;
	int (*adjacent)[3] = mesh->face_link_E;

	Node** ridge_edge = mesh->ridge_edge;
	Node** ravine_edge = mesh->ravine_edge;

	if(ridge_edge != NULL){
		for(int i=0; i<vertex_N; i++)
			delete ridge_edge[i];
	}
	if(ravine_edge != NULL){
		for(int i=0; i<vertex_N; i++)
			delete ravine_edge[i];
	}
	if(both_edge != NULL){
		for(int i=0; i<vertex_N; i++)
			delete both_edge[i];
	}

	ridge_edge = mesh->ridge_edge = new Node*[vertex_N];
	ravine_edge = mesh->ravine_edge = new Node*[vertex_N];
	both_edge = new Node*[vertex_N];
	for(int i=0; i<vertex_N; i++){
		ridge_edge[i] = new Node;
		ravine_edge[i] = new Node;
		both_edge[i] = new Node;
	}

	for(i=0; i<face_N; i++){
		if(face[i][0] < 0)
			continue;

		for(int j=0; j<3; j++){
			int pair = adjacent[i][j];
			BOOL flag = false;
			BOOL ridge = true;
			if(pair < 0)
				flag = true;
			else if(i > pair && MeshData::DOT(normal[i], normal[pair]) < T){
				flag = true;

				float c1[3];
				c1[0] = (vertex[face[i][0]][0] + vertex[face[i][1]][0] + vertex[face[i][2]][0])/3;
				c1[1] = (vertex[face[i][0]][1] + vertex[face[i][1]][1] + vertex[face[i][2]][1])/3;
				c1[2] = (vertex[face[i][0]][2] + vertex[face[i][1]][2] + vertex[face[i][2]][2])/3;

				float c2[3];
				c2[0] = (vertex[face[pair][0]][0] + vertex[face[pair][1]][0] + vertex[face[pair][2]][0])/3;
				c2[1] = (vertex[face[pair][0]][1] + vertex[face[pair][1]][1] + vertex[face[pair][2]][1])/3;
				c2[2] = (vertex[face[pair][0]][2] + vertex[face[pair][1]][2] + vertex[face[pair][2]][2])/3;

				float v[3];
				MeshData::VEC(v, c1, c2);
				float len = MeshData::LENGTH(v);
				float dot1 = MeshData::DOT(v, normal[i])/len;
				float dot2 = -MeshData::DOT(v, normal[pair])/len;
				if(acos(dot1) + acos(dot2) < PI)
					ridge = false;
			}
			if(flag){
				if(ridge){
					ridge_edge[face[i][j]]->append(face[i][(j+1)%3], -1);
					ridge_edge[face[i][(j+1)%3]]->append(face[i][j], -1);
				}
				else{
					ravine_edge[face[i][j]]->append(face[i][(j+1)%3], -1);
					ravine_edge[face[i][(j+1)%3]]->append(face[i][j], -1);
				}
				both_edge[face[i][j]]->append(face[i][(j+1)%3], -1);
				both_edge[face[i][(j+1)%3]]->append(face[i][j], -1);
			}
		}
	}
}

void Simplifier::deleteShortFeature(float T)
{
	int vertex_N = mesh->vertex_N;
	float (*vertex)[3] = mesh->vertex;
	int *degree = new int[vertex_N];

	Node** ridge_edge = mesh->ridge_edge;
	Node** ravine_edge = mesh->ravine_edge;

	for(int i=0; i<vertex_N; i++){
		degree[i] = 0;
		for(Node* current = ridge_edge[i]; current->next != NULL; current = current->next)
			degree[i]++;
	}
	BOOL is_changed = true;
	while(is_changed){
		is_changed = false;
		for(i=0; i<vertex_N; i++){
			if(degree[i] != 1)
				continue;

			int i1 = i;
			int i2 = ridge_edge[i]->v;
			float len = 0;
			while(true){
				len += MeshData::DIST(vertex[i1], vertex[i2]);

				if(degree[i2] != 2)
					break;

				int next;
				if(ridge_edge[i2]->v != i1)
					next = ridge_edge[i2]->v;
				else
					next = ridge_edge[i2]->next->v;
				i1 = i2;
				i2 = next;
			}
			if(len < T){
				is_changed = true;
				i1 = i;
				i2 = ridge_edge[i]->v;
				while(true){
					degree[i1] = 0;
					delete ridge_edge[i1];
					ridge_edge[i1] = new Node;
					if(degree[i2] != 2)
						break;
					int next;
					if(ridge_edge[i2]->v != i1)
						next = ridge_edge[i2]->v;
					else
						next = ridge_edge[i2]->next->v;
					
					i1 = i2;
					i2 = next;
				}
				degree[i2]--;
				Node* tmp = new Node;
				for(Node* c = ridge_edge[i2]; c->next != NULL; c = c->next){
					if(c->v != i1)
						tmp->append(c->v, -1);
				}
				delete ridge_edge[i2];
				ridge_edge[i2] = tmp;
			}
		}
	}

	for(i=0; i<vertex_N; i++){
		degree[i] = 0;
		for(Node* current = ravine_edge[i]; current->next != NULL; current = current->next)
			degree[i]++;
	}
	is_changed = true;
	while(is_changed){
		is_changed = false;
		for(i=0; i<vertex_N; i++){
			if(degree[i] != 1)
				continue;

			int i1 = i;
			int i2 = ravine_edge[i]->v;
			float len = 0;
			while(true){
				len += MeshData::DIST(vertex[i1], vertex[i2]);

				if(degree[i2] != 2)
					break;

				int next;
				if(ravine_edge[i2]->v != i1)
					next = ravine_edge[i2]->v;
				else
					next = ravine_edge[i2]->next->v;
				i1 = i2;
				i2 = next;
			}
			if(len < T){
				is_changed = true;
				i1 = i;
				i2 = ravine_edge[i]->v;
				while(true){
					degree[i1] = 0;
					delete ravine_edge[i1];
					ravine_edge[i1] = new Node;
					if(degree[i2] != 2)
						break;
					int next;
					if(ravine_edge[i2]->v != i1)
						next = ravine_edge[i2]->v;
					else
						next = ravine_edge[i2]->next->v;
					
					i1 = i2;
					i2 = next;
				}
				degree[i2]--;
				Node* tmp = new Node;
				for(Node* c = ravine_edge[i2]; c->next != NULL; c = c->next){
					if(c->v != i1)
						tmp->append(c->v, -1);
				}
				delete ravine_edge[i2];
				ravine_edge[i2] = tmp;
			}
		}
	}

	delete[] degree;
}

void Simplifier::deleteSmallCycle(float T)
{
	int vertex_N = mesh->vertex_N;
	float (*vertex)[3] = mesh->vertex;
	int *degree = new int[vertex_N];

	Node** ridge_edge = mesh->ridge_edge;
	Node** ravine_edge = mesh->ravine_edge;

	for(int i=0; i<vertex_N; i++){
		degree[i] = 0;
		for(Node* current = ridge_edge[i]; current->next != NULL; current = current->next)
			degree[i]++;
	}
	for(i=0; i<vertex_N; i++){
		if(degree[i] == 1)
			continue;

		for(Node* current = ridge_edge[i]; current->next != NULL; current = current->next){
			int index = checkSmallRidgeCycle(i, i, current->v, MeshData::DIST(vertex[i], vertex[current->v]), T);
			if(index >= 0 && index != current->v)
				break;
		}
		if(current->next != NULL){
			int i1 = i;
			int i2 = current->v;
			degree[i1]--;
			degree[i2]--;

			Node* tmp = new Node;
			for(Node* c = ridge_edge[i1]; c->next != NULL; c = c->next){
				if(c->v != i2)
					tmp->append(c->v, -1);
			}
			delete ridge_edge[i1];
			ridge_edge[i1] = tmp;

			tmp = new Node;
			for(c = ridge_edge[i2]; c->next != NULL; c = c->next){
				if(c->v != i1)
					tmp->append(c->v, -1);
			}
			delete ridge_edge[i2];
			ridge_edge[i2] = tmp;
		}
	}

	for(i=0; i<vertex_N; i++){
		degree[i] = 0;
		for(Node* current = ravine_edge[i]; current->next != NULL; current = current->next)
			degree[i]++;
	}
	for(i=0; i<vertex_N; i++){
		if(degree[i] == 1)
			continue;

		for(Node* current = ravine_edge[i]; current->next != NULL; current = current->next){
			int index = checkSmallRavineCycle(i, i, current->v, MeshData::DIST(vertex[i], vertex[current->v]), T);
			if(index >= 0 && index != current->v)
				break;
		}
		if(current->next != NULL){
			int i1 = i;
			int i2 = current->v;
			degree[i1]--;
			degree[i2]--;

			Node* tmp = new Node;
			for(Node* c = ravine_edge[i1]; c->next != NULL; c = c->next){
				if(c->v != i2)
					tmp->append(c->v, -1);
			}
			delete ravine_edge[i1];
			ravine_edge[i1] = tmp;

			tmp = new Node;
			for(c = ravine_edge[i2]; c->next != NULL; c = c->next){
				if(c->v != i1)
					tmp->append(c->v, -1);
			}
			delete ravine_edge[i2];
			ravine_edge[i2] = tmp;
		}
	}

	delete[] degree;
}

int Simplifier::checkSmallRidgeCycle(int start, int previous, int next, float len, float T)
{
	Node** ridge_edge = mesh->ridge_edge;

	if(len > T)
		return -1;
	if(next == start)
		return previous;
	for(Node* current = ridge_edge[next]; current->next!=NULL; current=current->next){
		int index;
		if(current->v != previous && 
			(index = checkSmallRidgeCycle(start, next, current->v, len+MeshData::DIST(mesh->vertex[next], mesh->vertex[current->v]),T)) >= 0)
			return index;
	}
	return -1;
}

int Simplifier::checkSmallRavineCycle(int start, int previous, int next, float len, float T)
{
	Node** ravine_edge = mesh->ravine_edge;

	if(len > T)
		return -1;
	if(next == start)
		return previous;
	for(Node* current = ravine_edge[next]; current->next!=NULL; current=current->next){
		int index;
		if(current->v != previous && 
			(index = checkSmallRavineCycle(start, next, current->v, len+MeshData::DIST(mesh->vertex[next], mesh->vertex[current->v]),T)) >= 0)
			return index;
	}
	return -1;
}

void Simplifier::connectNearEdge(float T)
{
	/*
	int vertex_N = mesh->vertex_N;
	float (*vertex)[3] = mesh->vertex;
	int *degree = new int[vertex_N];

	for(int i=0; i<vertex_N; i++){
		degree[i] = 0;
		for(Node* current = ridge_edge[i]; current->next != NULL; current = current->next)
			degree[i]++;
	}
	for(i=0; i<vertex_N; i++){
		if(degree[i] != 1)
			continue;

		Node* que = new Node;
		que->append(i, -1);

		Node* visit = new Node;
		visit->append(i, 0);

		while(que->next!=NULL){
			int now = que->v;
			Node* next = que->next;
			next->tail = que->tail;
			que->next = NULL;
			delete que;
			que = next;

			int *nei, nei_N;
			mesh->getConsistent1Ring(now, nei, nei_N);

			for(int j=0; j<nei_n; j++){
				if(ridge_edge[now]->doesInclude(nei[j]) == NULL)
					continue;

			}
		}

		for(Node* current = ridge_edge[i]; current->next != NULL; current = current->next){
			if(checkSmallRidgeCycle(i, i, current->v, MeshData::DIST(vertex[i], vertex[current->v]), T))
				break;
		}
		if(current->next != NULL){
			int i1 = i;
			int i2 = current->v;
			degree[i1]--;
			degree[i2]--;

			Node* tmp = new Node;
			for(Node* c = ridge_edge[i1]; c->next != NULL; c = c->next){
				if(c->v != i2)
					tmp->append(c->v, -1);
			}
			delete ridge_edge[i1];
			ridge_edge[i1] = tmp;

			tmp = new Node;
			for(c = ridge_edge[i2]; c->next != NULL; c = c->next){
				if(c->v != i1)
					tmp->append(c->v, -1);
			}
			delete ridge_edge[i2];
			ridge_edge[i2] = tmp;
		}
	}

	for(i=0; i<vertex_N; i++){
		degree[i] = 0;
		for(Node* current = ravine_edge[i]; current->next != NULL; current = current->next)
			degree[i]++;
	}
	for(i=0; i<vertex_N; i++){
		if(degree[i] == 1)
			continue;

		for(Node* current = ravine_edge[i]; current->next != NULL; current = current->next){
			if(checkSmallRavineCycle(i, i, current->v, MeshData::DIST(vertex[i], vertex[current->v]), T))
				break;
		}
		if(current->next != NULL){
			int i1 = i;
			int i2 = current->v;
			degree[i1]--;
			degree[i2]--;

			Node* tmp = new Node;
			for(Node* c = ravine_edge[i1]; c->next != NULL; c = c->next){
				if(c->v != i2)
					tmp->append(c->v, -1);
			}
			delete ravine_edge[i1];
			ravine_edge[i1] = tmp;

			tmp = new Node;
			for(c = ravine_edge[i2]; c->next != NULL; c = c->next){
				if(c->v != i1)
					tmp->append(c->v, -1);
			}
			delete ravine_edge[i2];
			ravine_edge[i2] = tmp;
		}
	}

	delete[] degree;
	*/
}

void Simplifier::cutMeddleTagEdge()
{
	/*
	int vertex_N = mesh->vertex_N;
	float (*vertex)[3] = mesh->vertex;
	int *degree = new int[vertex_N];

	for(int i=0; i<vertex_N; i++){
		degree[i] = 0;
		for(Node* current = ridge_edge[i]; current->next != NULL; current = current->next)
			degree[i]++;
	}
	for(i=0; i<vertex_N; i++){
		if(degree[i] == 2)
			continue;
			
		for(Node* current=ridge_edge[i]; current->next!=NULL; current=current->next){
			int i1 = i;
			int i2 = current->v;
		
			while(true){
				if(degree[i2] != 2 || i2 == i)
					break;

				int next;
				if(ridge_edge[i2]->v != i1)
					next = ridge_edge[i2]->v;
				else
					next = ridge_edge[i2]->next->v;
				i1 = i2;
				i2 = next;
			}
			if(len < T){
				is_changed = true;
				i1 = i;
				i2 = ridge_edge[i]->v;
				while(true){
					degree[i1] = 0;
					delete ridge_edge[i1];
					ridge_edge[i1] = new Node;
					if(degree[i2] != 2)
						break;
					int next;
					if(ridge_edge[i2]->v != i1)
						next = ridge_edge[i2]->v;
					else
						next = ridge_edge[i2]->next->v;
					
					i1 = i2;
					i2 = next;
				}
				degree[i2]--;
				Node* tmp = new Node;
				for(Node* c = ridge_edge[i2]; c->next != NULL; c = c->next){
					if(c->v != i1)
						tmp->append(c->v, -1);
				}
				delete ridge_edge[i2];
				ridge_edge[i2] = tmp;
			}
		}



		for(Node* current = ridge_edge[i]; current->next != NULL; current = current->next){
			if(checkSmallRidgeCycle(i, i, current->v, MeshData::DIST(vertex[i], vertex[current->v]), T))
				break;
		}
		if(current->next != NULL){
			int i1 = i;
			int i2 = current->v;
			degree[i1]--;
			degree[i2]--;

			Node* tmp = new Node;
			for(Node* c = ridge_edge[i1]; c->next != NULL; c = c->next){
				if(c->v != i2)
					tmp->append(c->v, -1);
			}
			delete ridge_edge[i1];
			ridge_edge[i1] = tmp;

			tmp = new Node;
			for(c = ridge_edge[i2]; c->next != NULL; c = c->next){
				if(c->v != i1)
					tmp->append(c->v, -1);
			}
			delete ridge_edge[i2];
			ridge_edge[i2] = tmp;
		}
	}

	for(i=0; i<vertex_N; i++){
		degree[i] = 0;
		for(Node* current = ravine_edge[i]; current->next != NULL; current = current->next)
			degree[i]++;
	}
	for(i=0; i<vertex_N; i++){
		if(degree[i] == 1)
			continue;

		for(Node* current = ravine_edge[i]; current->next != NULL; current = current->next){
			if(checkSmallRavineCycle(i, i, current->v, MeshData::DIST(vertex[i], vertex[current->v]), T))
				break;
		}
		if(current->next != NULL){
			int i1 = i;
			int i2 = current->v;
			degree[i1]--;
			degree[i2]--;

			Node* tmp = new Node;
			for(Node* c = ravine_edge[i1]; c->next != NULL; c = c->next){
				if(c->v != i2)
					tmp->append(c->v, -1);
			}
			delete ravine_edge[i1];
			ravine_edge[i1] = tmp;

			tmp = new Node;
			for(c = ravine_edge[i2]; c->next != NULL; c = c->next){
				if(c->v != i1)
					tmp->append(c->v, -1);
			}
			delete ravine_edge[i2];
			ravine_edge[i2] = tmp;
		}
	}

	delete[] degree;*/
}
/*
BOOL Simplifier::contractEdge(int i1, int i2, float new_vertex[], int &i3, int &i4)
{
	float (*vertex)[3] = mesh->vertex;
	int **link_v = mesh->vertex_link_v;
	int **link_f = mesh->vertex_link_f;
	int *degree_v = mesh->degree_v;
	int *degree_f = mesh->degree_f;
	int (*face)[3] = mesh->face;
	float (*normal)[3] = mesh->normal_f;
	int *isBound = mesh->isBound;

	if(isBound[i1] == NON_MANIFOLD || isBound[i2] == NON_MANIFOLD)
		return false;

	//adjacent data for i1
	int *nei1, nei1_N, *face1, nei1f_N;
	mesh->getConsistent1Ring(i1, nei1, face1, nei1_N);
	nei1f_N = degree_f[i1];
	
	//adjacent data for i2
	int *nei2, nei2_N, *face2, nei2f_N;
	mesh->getConsistent1Ring(i2, nei2, face2, nei2_N);
	nei2f_N = degree_f[i2];

	//degree check
	if(isBound[i1] == BOUNDARY && isBound[i2] == BOUNDARY){
		if(nei1[0] == i2 || nei1[nei1_N-1] == i2 
			|| nei2[0] == i1 || nei2[nei2_N-1] == i1){
			if(nei1_N +  nei2_N < 6){
				return false;
			}
		}
		else if(nei1_N +  nei2_N < 7){
			return false;
		}
	}
	else if(nei1_N +  nei2_N < 7){
		return false;
	}

	//opposite vertcies to the contracted edge;
	i3 = -1;
	i4 = -1;
	//degenerated face index;
	int f1 = -1;
	int f2 = -1;

	for(int i=0; i<nei1f_N; i++){
		if(nei1[i] == i2){
			i3 = nei1[(i+1)%nei1_N];
			f1 = face1[i];
		}
	}
	if(i3 >= 0 && mesh->isBound[i3] == NON_MANIFOLD)
		return false;

	for(i=0; i<nei2f_N; i++){
		if(nei2[i] == i1){
			i4 = nei2[(i+1)%nei2_N];
			f2 = face2[i];
		}
	}
	if(i4 >= 0 && mesh->isBound[i4] == NON_MANIFOLD)
		return false;

	if(i3 >= 0 && i4 >=0 && i3 == i4)
		return false;

	for(i=0; i<nei1_N; i++){
		if(nei1[i] == i2 || nei1[i] == i3 || nei1[i] == i4)
			continue;
		for(int j=0; j<nei2_N; j++)
			if(nei2[j] != i1 && nei1[i] == nei2[j])
				return false;
	}

	int *nei3, nei3_N, *face3, nei3f_N;
	if(i3 >= 0){
		mesh->getConsistent1Ring(i3, nei3, face3, nei3_N);
		nei3f_N = degree_f[i3];
	}

	int *nei4, nei4_N, *face4, nei4f_N;
	if(i4 >= 0){
		mesh->getConsistent1Ring(i4, nei4, face4,  nei4_N);
		nei4f_N = degree_f[i4];
	}

	if((i3 >= 0 && nei3_N < 4) || (i4 >= 0 && nei4_N < 4))
		return false;

	//flip check by "zone" technique
	BOOL flip = false;
	for(i=0; i<nei1f_N; i++){
		if(face1[i] == f1 || face1[i] == f2)
			continue;
		float v[3];
		MeshData::VEC(v, vertex[nei1[i]], vertex[nei1[(i+1)%nei1_N]]);
		double n[3];
		MeshData::CROSS(n, v, normal[face1[i]]);
		double len = MeshData::LENGTH(n);
		if(len == 0)
			continue;
		n[0] /= len;
		n[1] /= len;
		n[2] /= len;
		MeshData::VEC(v, vertex[i1], vertex[nei1[i]]);
		double side1 = MeshData::DOT(n, v);
		MeshData::VEC(v, new_vertex, vertex[nei1[i]]);
		double side2 = MeshData::DOT(n, v);
		if(side1 * side2 <= 0){
			return false;
		}
	}
	for(i=0; i<nei2f_N; i++){
		if(face2[i] == f1 || face2[i] == f2)
			continue;
		float v[3];
		MeshData::VEC(v, vertex[nei2[i]], vertex[nei2[(i+1)%nei2_N]]);
		double n[3];
		MeshData::CROSS(n, v, normal[face2[i]]);
		double len = MeshData::LENGTH(n);
		if(len == 0)
			continue;
		n[0] /= len;
		n[1] /= len;
		n[2] /= len;
		MeshData::VEC(v, vertex[i2], vertex[nei2[i]]);
		double side1 = MeshData::DOT(n, v);
		MeshData::VEC(v, new_vertex, vertex[nei2[i]]);
		double side2 = MeshData::DOT(n, v);
		if(side1 * side2 <= 0){
			return false;
		}
	}

	//update vertex position;
	vertex[i1][0] = new_vertex[0];
	vertex[i1][1] = new_vertex[1];
	vertex[i1][2] = new_vertex[2];

	//Updated adjacent vertex data;
	//for vertex i1 and i2;
	int link_N = nei1f_N + nei2f_N - 4;
	if(f1 < 0)
		link_N += 2;
	else if(f2 < 0)
		link_N += 2;
	int* link = new int[3*link_N];
	link_N = 0;
	for(i=0; i<nei1f_N; i++){
		if(face1[i] != f1 && face1[i] != f2){
			link[3*link_N] = nei1[i];
			link[3*link_N+1] = nei1[(i+1)%nei1_N];
			link[3*link_N+2] = face1[i];
			link_N++;
		}
	}
	for(i=0; i<nei2f_N; i++){
		if(face2[i] != f1 && face2[i] != f2){
			link[3*link_N] = nei2[i];
			link[3*link_N+1] = nei2[(i+1)%nei2_N];
			link[3*link_N+2] = face2[i];

			if(face[face2[i]][0] == i2)
				face[face2[i]][0] = i1;
			else if(face[face2[i]][1] == i2)
				face[face2[i]][1] = i1;
			else 
				face[face2[i]][2] = i1;

			link_N++;
		}
		else
			face[face2[i]][0] = face[face2[i]][1] = face[face2[i]][2] = -1;
	}
	delete[] link_v[i1];
	delete[] link_f[i1];
	isBound[i1] = mesh->sortVertexLink(link, link_N, link_v[i1], link_f[i1], degree_v[i1], degree_f[i1]);
	delete[] link;

	//upate link infomation incident to i2
	for(i=0; i<nei2_N; i++){
		if(nei2[i] == i1)
			continue;
		int* l = link_v[nei2[i]];
		int N = degree_v[nei2[i]];
		for(int j=0; j<N; j++)
			if(l[j] == i2)
				l[j] = i1;
	}
	delete[] link_v[i2];
	link_v[i2] = NULL;
	delete[] link_f[i2];
	link_f[i2] = NULL;
	degree_v[i2] = 0;
	degree_f[i2] = 0;
	isBound[i2] = NON_MANIFOLD;

	//upadate face normal incident to i1
	int* new_nei = link_v[i1];
	int* new_face = link_f[i1];
	int new_nei_N = degree_v[i1];
	int new_neif_N = degree_f[i1];
	for(i=0; i<new_neif_N; i++){
		float v1[3], v2[3];
		MeshData::VEC(v1, vertex[i1], vertex[new_nei[i]]);
		MeshData::VEC(v2, vertex[i1], vertex[new_nei[(i+1)%new_nei_N]]);
		double n[3];
		MeshData::CROSS(n, v1, v2);
		double len = MeshData::LENGTH(n);

		int f = new_face[i];
		if((float)len != 0){
			normal[f][0] = (float)(n[0]/len);
			normal[f][1] = (float)(n[1]/len);
			normal[f][2] = (float)(n[2]/len);
		}
		else
			normal[f][0] = normal[f][1] = normal[f][2] = 0;
	}

	//for vertex i3;
	if(i3 >= 0){
		link_N = nei3f_N - 1;
		link = new int[3*link_N];
		link_N = 0;
		for(i=0; i<nei3f_N; i++){
			if(face3[i] != f1){
				link[3*link_N] = nei3[i];
				link[3*link_N+1] = nei3[(i+1)%nei3_N];
				link[3*link_N+2] = face3[i];
				link_N++;
			}
		}
		delete[] link_v[i3];
		delete[] link_f[i3];
		isBound[i3] = mesh->sortVertexLink(link, link_N, link_v[i3], link_f[i3], 
				degree_v[i3], degree_f[i3]);
		delete[] link;
	}

	//for vertex i4
	if(i4 >= 0){
		link_N = nei4f_N - 1;
		link = new int[3*link_N];
		link_N = 0;
		for(i=0; i<nei4f_N; i++){
			if(face4[i] != f2){
				link[3*link_N] = nei4[i];
				link[3*link_N+1] = nei4[(i+1)%nei4_N];
				link[3*link_N+2] = face4[i];
				link_N++;
			}
		}
		delete[] link_v[i4];
		delete[] link_f[i4];
		isBound[i4] = mesh->sortVertexLink(link, link_N, link_v[i4], link_f[i4], 
				degree_v[i4], degree_f[i4]);
		delete[] link;
	}

	return true;
}*/

void Simplifier::halfEdgeCollapse(float T, BOOL is_ridge)
{
	/*
	float (*vertex)[3] = mesh->vertex;
	Node** v_link = mesh->vertex_link;
	int vertex_N = mesh->vertex_N;
	Node** e_table = (Node**)new Node[vertex_N];

	Node **tag_edge;
	Node **cluster;
	if(is_ridge){
		tag_edge = mesh->ridge_edge;
		ridge_cluster = (Node**)new Node[vertex_N];
		for(int i=0; i<vertex_N; i++){
			ridge_cluster[i] = new Node;
			ridge_cluster[i]->append(i, -1);
		}
		cluster = ridge_cluster;
	}
	else{
		tag_edge = mesh->ravine_edge;
		ravine_cluster = (Node**)new Node[vertex_N];
		for(int i=0; i<vertex_N; i++){
			ravine_cluster[i] = new Node;
			ravine_cluster[i]->append(i, -1);
		}
		cluster = ravine_cluster;
	}

	int *degree = new int[vertex_N];
	for(int i=0; i<vertex_N; i++){
		degree[i] = 0;
		for(Node* current = tag_edge[i]; current->next != NULL; current = current->next)
			degree[i]++;
	}

	int edge_N = 0;
	for(i=0; i<vertex_N; i++){
		e_table[i] = new Node;
		for(Node* current = v_link[i]; current->next!=NULL; current = current->next->next){
			if(i > current->v){
				edge_N++;
			}
		}
	}

	int (*edge)[2] = new int[edge_N][2];
	float* dist = new float[edge_N];

	int counter = 0;
	for(i=0; i<vertex_N; i++){
		for(Node* current = v_link[i]; current->next!=NULL; current = current->next->next){
			if(i > current->v){
				edge[counter][0] = i;
				edge[counter][1] = current->v;
				dist[counter] = MeshData::DIST(vertex[i], vertex[current->v]);
				counter++;
			}
		}
	}

	int* heap = new int[edge_N+1];
	int* index = new int[edge_N];
	int last_heap = 0;

	for(i=0; i<edge_N; i++){
		e_table[edge[i][0]]->append(i, -1);
		e_table[edge[i][1]]->append(i, -1);

		//insert;
		heap[++last_heap] = i;
		index[i] = last_heap;
		upheap(dist, last_heap, last_heap, heap, index);
	}

	while(last_heap != 0){
		int min = heap[1];
		index[min] = -1;
		//remove
		heap[1] = heap[last_heap--];
		index[heap[1]] = 1;
		downheap(dist, last_heap, 1, heap, index);

		if(edge[min][0] < 0)
			continue;

		if(dist[min] > T)
			break;

		
		BOOL is_tag_edge = false;
		for(Node* current=tag_edge[edge[min][0]]; current->next!=NULL; current=current->next){
			if(current->v == edge[min][1]){
				is_tag_edge = true;
				break;
			}
		}

		int con_v;
		if(degree[edge[min][0]] == 0 && degree[edge[min][1]] == 0){
			float sum1 = 0;
			for(current=e_table[edge[min][0]]; current->next!=NULL; current=current->next){
				if(edge[current->v][0] < 0)
					continue;
				sum1 += MeshData::DIST(vertex[edge[current->v][0]], vertex[edge[current->v][1]]);
			}
			float sum2 = 0;
			for(current=e_table[edge[min][1]]; current->next!=NULL; current=current->next){
				if(edge[current->v][0] < 0)
					continue;
				sum2 += MeshData::DIST(vertex[edge[current->v][0]], vertex[edge[current->v][1]]);
			}
			if(sum1 > sum2)
				con_v = edge[min][0];
			else
				con_v = edge[min][1];
		}
		else if(degree[edge[min][0]] == 0){
			con_v = edge[min][1];
		}
		else if(degree[edge[min][1]] == 0){
			con_v = edge[min][0];
		}
		else if(degree[edge[min][0]] == 1 && degree[edge[min][1]] == 1){
			if(is_tag_edge)
				continue;
			else{
				float sum1 = 0;
				for(current=e_table[edge[min][0]]; current->next!=NULL; current=current->next){
					if(edge[current->v][0] < 0)
						continue;
					sum1 += MeshData::DIST(vertex[edge[current->v][0]], vertex[edge[current->v][1]]);
				}
				float sum2 = 0;
				for(current=e_table[edge[min][1]]; current->next!=NULL; current=current->next){
					if(edge[current->v][0] < 0)
						continue;
					sum2 += MeshData::DIST(vertex[edge[current->v][0]], vertex[edge[current->v][1]]);
				}
				if(sum1 > sum2)
					con_v = edge[min][0];
				else
					con_v = edge[min][1];
			}
		}
		else if(degree[edge[min][0]] == 1){
			if(is_tag_edge){
				if(degree[edge[min][1]] == 2)
					con_v = edge[min][0];
				else
					continue;
			}
			else 
				con_v = edge[min][1];
		}
		else if(degree[edge[min][1]] == 1){
			if(is_tag_edge){
				if(degree[edge[min][1]] == 2)
					con_v = edge[min][1];
				else
					continue;
			}
			else 
				con_v = edge[min][0];
		}
		else if(degree[edge[min][0]] == 2 && degree[edge[min][1]] == 2){
			if(is_tag_edge){
				float sum1 = 0;
				for(current=e_table[edge[min][0]]; current->next!=NULL; current=current->next){
					if(edge[current->v][0] < 0)
						continue;
					sum1 += MeshData::DIST(vertex[edge[current->v][0]], vertex[edge[current->v][1]]);
				}
				float sum2 = 0;
				for(current=e_table[edge[min][1]]; current->next!=NULL; current=current->next){
					if(edge[current->v][0] < 0)
						continue;
					sum2 += MeshData::DIST(vertex[edge[current->v][0]], vertex[edge[current->v][1]]);
				}
				if(sum1 > sum2)
					con_v = edge[min][0];
				else
					con_v = edge[min][1];
			}
			else 
				continue;
		}
		else if(degree[edge[min][0]] == 2){
			if(is_tag_edge){
				con_v = edge[min][1];
			}
			else 
				continue;
		}
		else if(degree[edge[min][1]] == 2){
			if(is_tag_edge){
				con_v = edge[min][0];
			}
			else 
				continue;
		}
		else{
			continue;
		}
		if(con_v == edge[min][1]){
			int t = edge[min][0];
			edge[min][0] = edge[min][1];
			edge[min][1] = t;
		}

		

		if(!contractEdge(edge[min][0], edge[min][1], vertex[con_v]))
			continue;

		//update edge infomation
		Node* tmp = new Node;
		int c = 0;
		for(current = e_table[edge[min][0]]; current->next!=NULL; current=current->next){
			if(edge[current->v][0] < 0)
				continue;
			if(current->v != min){
				c++;
				tmp->append(current->v, -1);
			}
		}
		for(current = e_table[edge[min][1]]; current->next!=NULL; current=current->next){
			if(current->v == min || edge[current->v][0] < 0)
				continue;
			int pair1;
			if(edge[current->v][1] == edge[min][1]){
				pair1 = edge[current->v][0];
				edge[current->v][1] = edge[min][0];
			}
			else{
				pair1 = edge[current->v][1];
				edge[current->v][0] = edge[min][0];
			}

			i = 0;
			BOOL flag = true;
			for(Node* search = tmp; i < c; search=search->next){
				int pair2;
				if(edge[search->v][0] != edge[min][0])
					pair2 = edge[search->v][0];
				else
					pair2 = edge[search->v][1];
				if(pair1 == pair2){
					flag = false;
					break;
				}
				i++;
			}
			if(flag)
				tmp->append(current->v, -1);
			else{
				//delete
				dist[current->v] = -1;
				if(index[current->v] > 0){
					upheap(dist, last_heap, index[current->v], heap, index);
					index[current->v] = -1;
					heap[1] = heap[last_heap--];
					index[heap[1]] = 1;
					downheap(dist, last_heap, 1, heap, index);
					edge[current->v][0] = edge[current->v][1] = -1;
				}
			}
		}
		delete e_table[edge[min][0]];
		e_table[edge[min][0]] = tmp;
		delete e_table[edge[min][1]];
		e_table[edge[min][1]] = new Node;

		for(current = e_table[edge[min][0]]; current->next!=NULL; current=current->next){
			int v = current->v;
			dist[v] = MeshData::DIST(vertex[edge[v][0]], vertex[edge[v][1]]);
			if(index[v] < 0){
				//insert
				heap[++last_heap] = v;
				index[v] = last_heap;
				upheap(dist, last_heap, last_heap, heap, index);
			}
			else{
				//change;
				if(index[v] != 1 && dist[v] < dist[heap[index[v]/2]])
					upheap(dist, last_heap, index[v], heap, index);
				else
					downheap(dist, last_heap, index[v], heap, index);
			}
		}

		//update vertex cluster
		tmp = new Node;
		for(current=cluster[edge[min][0]]; current->next!=NULL; current=current->next){
			tmp->append(current->v, -1);
		}
		for(current=cluster[edge[min][1]]; current->next!=NULL; current=current->next){
			tmp->append(current->v, -1);
		}
		delete cluster[edge[min][0]];
		delete cluster[edge[min][1]];
		cluster[edge[min][0]] = tmp;
		cluster[edge[min][1]] = new Node;

		//upate tag infomarion
		if(degree[edge[min][0]] == 0 && degree[edge[min][1]] == 0)
			continue;

		tmp = new Node;
		c = 0;
		degree[edge[min][0]] = 0;
		degree[edge[min][1]] = 0;
		for(current = tag_edge[edge[min][0]]; current->next!=NULL; current=current->next){
			if(current->v != edge[min][1]){
				tmp->append(current->v, -1);
				degree[edge[min][0]]++;
				c++;
			}
		}
		for(current = tag_edge[edge[min][1]]; current->next!=NULL; current=current->next){
			if(current->v == edge[min][0])
				continue;
			i = 0;
			BOOL flag = true;
			for(Node* search = tmp; i<c; search=search->next){
				if(current->v == search->v){
					flag = false;
					break;
				}
				i++;
			}
			if(flag){
				tmp->append(current->v, -1);
				degree[edge[min][0]]++;
				for(Node* search = tag_edge[current->v]; search->next!=NULL; search=search->next){
					if(search->v == edge[min][1]){
						search->v = edge[min][0];
						break;
					}
				}
			}
			else{
				flag = false;
				for(Node* search = tag_edge[current->v]; search->next!=NULL; search=search->next){
					if(search->v == edge[min][0]){
						flag = true;
						break;
					}
				}
				for(search = tag_edge[current->v]; search->next!=NULL; search=search->next){
					if(search->v == edge[min][1]){
						if(flag){
							Node *tmp2 = new Node;
							for(Node* current2 = tag_edge[current->v]; current2->next!=NULL; current2=current2->next){
								if(current2->v != edge[min][1])
									tmp2->append(current2->v, -1);
							}
							delete tag_edge[current->v];
							tag_edge[current->v] = tmp2;
							degree[current->v]--;
						}
						else
							search->v = edge[min][0];
						break;
					}
				}
			}
		}
		delete tag_edge[edge[min][0]];
		tag_edge[edge[min][0]] = tmp;
		delete tag_edge[edge[min][1]];
		tag_edge[edge[min][1]] = new Node;
	}

	delete[] edge;
	*/
}

inline void Simplifier::upheap(float *a, int N, int k, int *p, int *q)
{
	int v;
	v = p[k];
	while(k > 1 && a[p[k/2]] >= a[v]){
		p[k] = p[k/2]; q[p[k/2]] = k; k = k/2;
	}
	p[k] = v; q[v] = k;
}

inline void Simplifier::downheap(float *a, int N, int k, int *p, int *q)
{
	int j, v;
	v = p[k];
	while(k <= N/2){
		j = k+k;
		if(j < N && a[p[j]] > a[p[j+1]]) j++;
		if(a[v] <= a[p[j]]) break;
		p[k] = p[j]; q[p[j]] = k; k = j;
	}
	p[k] = v; q[v] = k;
}

void Simplifier::halfEdgeCollapse(int until)
{
	float (*vertex)[3] = mesh->vertex;
	int** v_link = mesh->vertex_link_v;
	int vertex_N = mesh->vertex_N;
	int current_N = mesh->countValidVertex();
	//Node** e_table = (Node**)new Node[vertex_N];
	int** e_table = new int*[vertex_N];
	int* e_table_N = new int[vertex_N];
	int* degree_v = mesh->degree_v;
	float T = 0;
	int* isBound = mesh->isBound;
	//Node **tag_edge;

	/*
	if(both_edge == NULL){
		both_edge = new Node*[vertex_N];
		for(int i=0; i<vertex_N; i++)
			both_edge[i] = new Node;
	}
	tag_edge = both_edge;
	*/

	/*
	int *degree = new int[vertex_N];
	count tag degree
	*/

	int edge_N = 0;
	for(int i=0; i<vertex_N; i++){
		e_table[i] = new int[degree_v[i]];
		e_table_N[i] = 0;
		for(int j=0; j<degree_v[i]; j++)
			if(i > v_link[i][j] 
				&& (isBound[i] != NON_MANIFOLD && isBound[v_link[i][j]] != NON_MANIFOLD))
				edge_N++;
	}

	/*
	int edge_N = 0;
	for(i=0; i<vertex_N; i++){
		e_table[i] = new Node;
		for(Node* current = v_link[i]; current->next!=NULL; current = current->next->next){
			if(i > current->v){
				edge_N++;
			}
		}
	}
	*/

	int (*edge)[2] = new int[edge_N][2];
	float* dist = new float[edge_N];

	int counter = 0;
	for(i=0; i<vertex_N; i++){
		for(int j=0; j<degree_v[i]; j++)
			if(i > v_link[i][j]
				&& (isBound[i] != NON_MANIFOLD && isBound[v_link[i][j]] != NON_MANIFOLD)){
				edge[counter][0] = i;
				edge[counter][1] = v_link[i][j];
				dist[counter] = (float)MeshData::DIST(vertex[i], vertex[v_link[i][j]]);
				counter++;
			}
	}
/*
		for(Node* current = v_link[i]; current->next!=NULL; current = current->next->next){
			if(i > current->v){
				edge[counter][0] = i;
				edge[counter][1] = current->v;
				dist[counter] = MeshData::DIST(vertex[i], vertex[current->v]);
				counter++;
			}
		}
	}*/

	int* heap = new int[edge_N+1];
	int* index = new int[edge_N];
	int last_heap = 0;

	for(i=0; i<edge_N; i++){
		e_table[edge[i][0]][e_table_N[edge[i][0]]++] = i;
		e_table[edge[i][1]][e_table_N[edge[i][1]]++] = i;

		/*
		e_table[edge[i][0]]->append(i, -1);
		e_table[edge[i][1]]->append(i, -1);
		*/

		//insert;
		heap[++last_heap] = i;
		index[i] = last_heap;
		upheap(dist, last_heap, last_heap, heap, index);
	}

	while(last_heap != 0){
		int min = heap[1];
		index[min] = -1;
		//remove
		heap[1] = heap[last_heap--];
		index[heap[1]] = 1;
		downheap(dist, last_heap, 1, heap, index);

		int i1 = edge[min][0];
		int i2 = edge[min][1];

		if(i1 < 0 || i2 < 0)
			continue;
		/*
		BOOL is_tag_edge = false;
		for(Node* current=tag_edge[edge[min][0]]; current->next!=NULL; current=current->next){
			if(current->v == edge[min][1]){
				is_tag_edge = true;
				break;
			}
		}
		*/

		int con_v;
		//if(degree[edge[min][0]] == 0 && degree[edge[min][1]] == 0){
			float sum1 = 0;
			for(i=0; i<e_table_N[i1]; i++)
				sum1 += dist[e_table[i1][i]];
			/*
			for(current=e_table[edge[min][0]]; current->next!=NULL; current=current->next){
				if(edge[current->v][0] < 0)
					continue;
				sum1 += MeshData::DIST(vertex[edge[current->v][0]], vertex[edge[current->v][1]]);
			}
			*/
			float sum2 = 0;
			for(i=0; i<e_table_N[i2]; i++)
				sum1 += dist[e_table[i2][i]];
			/*
			for(current=e_table[edge[min][1]]; current->next!=NULL; current=current->next){
				if(edge[current->v][0] < 0)
					continue;
				sum2 += MeshData::DIST(vertex[edge[current->v][0]], vertex[edge[current->v][1]]);
			}
			*/
			if(sum1 > sum2)
				con_v = i1;
			else
				con_v = i2;
		//}
			/*
		else if(degree[edge[min][0]] == 0){
			con_v = edge[min][1];
		}
		else if(degree[edge[min][1]] == 0){
			con_v = edge[min][0];
		}
		else if(degree[edge[min][0]] == 1 && degree[edge[min][1]] == 1){
			if(is_tag_edge)
				continue;
			else{
				float sum1 = 0;
				for(current=e_table[edge[min][0]]; current->next!=NULL; current=current->next){
					if(edge[current->v][0] < 0)
						continue;
					sum1 += MeshData::DIST(vertex[edge[current->v][0]], vertex[edge[current->v][1]]);
				}
				float sum2 = 0;
				for(current=e_table[edge[min][1]]; current->next!=NULL; current=current->next){
					if(edge[current->v][0] < 0)
						continue;
					sum2 += MeshData::DIST(vertex[edge[current->v][0]], vertex[edge[current->v][1]]);
				}
				if(sum1 > sum2)
					con_v = edge[min][0];
				else
					con_v = edge[min][1];
			}
		}
		else if(degree[edge[min][0]] == 1){
			if(is_tag_edge){
				if(degree[edge[min][1]] == 2)
					con_v = edge[min][0];
				else
					continue;
			}
			else 
				con_v = edge[min][1];
		}
		else if(degree[edge[min][1]] == 1){
			if(is_tag_edge){
				if(degree[edge[min][1]] == 2)
					con_v = edge[min][1];
				else
					continue;
			}
			else 
				con_v = edge[min][0];
		}
		else if(degree[edge[min][0]] == 2 && degree[edge[min][1]] == 2){
			if(is_tag_edge){
				float sum1 = 0;
				for(current=e_table[edge[min][0]]; current->next!=NULL; current=current->next){
					if(edge[current->v][0] < 0)
						continue;
					sum1 += MeshData::DIST(vertex[edge[current->v][0]], vertex[edge[current->v][1]]);
				}
				float sum2 = 0;
				for(current=e_table[edge[min][1]]; current->next!=NULL; current=current->next){
					if(edge[current->v][0] < 0)
						continue;
					sum2 += MeshData::DIST(vertex[edge[current->v][0]], vertex[edge[current->v][1]]);
				}
				if(sum1 > sum2)
					con_v = edge[min][0];
				else
					con_v = edge[min][1];
			}
			else 
				continue;
		}
		else if(degree[edge[min][0]] == 2){
			if(is_tag_edge){
				con_v = edge[min][1];
			}
			else 
				continue;
		}
		else if(degree[edge[min][1]] == 2){
			if(is_tag_edge){
				con_v = edge[min][0];
			}
			else 
				continue;
		}
		else{
			continue;
		}*/

		if(con_v == i2){
			int t = edge[min][0];
			edge[min][0] = edge[min][1];
			edge[min][1] = t;

			i1 = edge[min][0];
			i2 = edge[min][1];
		}

		//half edge collapse
		int i3, i4;
		if(!contractEdge(i1, i2, vertex[i1], i3, i4))
			continue;

		edge[min][0] = edge[min][1] = -1;

		//update edge infomation
		int* new_table;
		int c;
		int e3 = -1;
		int e4 = -1;
		
		//for i3
		if(i3 >= 0){
			new_table = new int[degree_v[i3]];
			c = 0;
			for(i=0; i<e_table_N[i3] ;i++){
				int e = e_table[i3][i];
				if(edge[e][0] == i2 || edge[e][1] == i2){
					e3 = e;
					edge[e][0] = -1;
					edge[e][1] = -1;
					//delete from heap
					if(index[e] > 0){
						dist[e] = -1;
						upheap(dist, last_heap, index[e], heap, index);
						index[e] = -1;
						heap[1] = heap[last_heap--];
						index[heap[1]] = 1;
						downheap(dist, last_heap, 1, heap, index);
					}
				}
				else
					new_table[c++] = e;
			}
			e_table_N[i3] = c;
			delete[] e_table[i3];
			e_table[i3] = new_table;
		}
		//for i4
		if(i4 >= 0){
			new_table = new int[degree_v[i4]];
			c = 0;
			for(i=0; i<e_table_N[i4] ;i++){
				int e = e_table[i4][i];
				if(edge[e][0] == i2 || edge[e][1] == i2){
					e4 = e;
					edge[e][0] = -1;
					edge[e][1] = -1;
					//delete from heap
					if(index[e] > 0){
						dist[e] = -1;
						upheap(dist, last_heap, index[e], heap, index);
						index[e] = -1;
						heap[1] = heap[last_heap--];
						index[heap[1]] = 1;
						downheap(dist, last_heap, 1, heap, index);
					}
				}
				else
					new_table[c++] = e;
			}
			e_table_N[i4] = c;
			delete[] e_table[i4];
			e_table[i4] = new_table;
		}

		//for i1
		c = 0;
		new_table = new int[degree_v[i1]];
		for(i=0; i<e_table_N[i1]; i++){
			int e = e_table[i1][i];
			if(e != min)
				new_table[c++] = e;
		}
		for(i=0; i<e_table_N[i2]; i++){
			int e = e_table[i2][i];
			if(e != min && e != e3 && e != e4){
				new_table[c++] = e;
				if(edge[e][0] == i2)
					edge[e][0] = i1;
				else
					edge[e][1] = i1;
			}
		}
		e_table_N[i1] = c;
		delete[] e_table[i1];
		e_table[i1] = new_table;

		//for i2 (delete)
		e_table_N[i2] = 0;
		delete[] e_table[i2]; 
		e_table[i2] = NULL;

		//update heap 
		for(i=0; i<e_table_N[i1]; i++){
			int e = e_table[i1][i];
			dist[e] = (float)MeshData::DIST(vertex[edge[e][0]], vertex[edge[e][1]]);
			if(index[e] < 0){
				//insert
				heap[++last_heap] = e;
				index[e] = last_heap;
				upheap(dist, last_heap, last_heap, heap, index);
			}
			else{
				//change;
				if(index[e] != 1 && dist[e] < dist[heap[index[e]/2]])
					upheap(dist, last_heap, index[e], heap, index);
				else
					downheap(dist, last_heap, index[e], heap, index);
			}
		}

	/*
		Node* tmp = new Node;
		int c = 0;
		for(current = e_table[edge[min][0]]; current->next!=NULL; current=current->next){
			if(edge[current->v][0] < 0)
				continue;
			if(current->v != min){
				c++;
				tmp->append(current->v, -1);
			}
		}
		for(current = e_table[edge[min][1]]; current->next!=NULL; current=current->next){
			if(current->v == min || edge[current->v][0] < 0)
				continue;
			int pair1;
			if(edge[current->v][1] == edge[min][1]){
				pair1 = edge[current->v][0];
				edge[current->v][1] = edge[min][0];
			}
			else{
				pair1 = edge[current->v][1];
				edge[current->v][0] = edge[min][0];
			}

			i = 0;
			BOOL flag = true;
			for(Node* search = tmp; i < c; search=search->next){
				int pair2;
				if(edge[search->v][0] != edge[min][0])
					pair2 = edge[search->v][0];
				else
					pair2 = edge[search->v][1];
				if(pair1 == pair2){
					flag = false;
					break;
				}
				i++;
			}
			if(flag)
				tmp->append(current->v, -1);
			else{
				//delete
				dist[current->v] = -1;
				if(index[current->v] > 0){
					upheap(dist, last_heap, index[current->v], heap, index);
					index[current->v] = -1;
					heap[1] = heap[last_heap--];
					index[heap[1]] = 1;
					downheap(dist, last_heap, 1, heap, index);
					edge[current->v][0] = edge[current->v][1] = -1;
				}
			}
		}
		delete e_table[edge[min][0]];
		e_table[edge[min][0]] = tmp;
		delete e_table[edge[min][1]];
		e_table[edge[min][1]] = new Node;

		for(current = e_table[edge[min][0]]; current->next!=NULL; current=current->next){
			int v = current->v;
			dist[v] = MeshData::DIST(vertex[edge[v][0]], vertex[edge[v][1]]);
			if(index[v] < 0){
				//insert
				heap[++last_heap] = v;
				index[v] = last_heap;
				upheap(dist, last_heap, last_heap, heap, index);
			}
			else{
				//change;
				if(index[v] != 1 && dist[v] < dist[heap[index[v]/2]])
					upheap(dist, last_heap, index[v], heap, index);
				else
					downheap(dist, last_heap, index[v], heap, index);
			}
		}

		//update vertex cluster
		/////////////comment start
		/*
		tmp = new Node;
		for(current=cluster[edge[min][0]]; current->next!=NULL; current=current->next){
			tmp->append(current->v, -1);
		}
		for(current=cluster[edge[min][1]]; current->next!=NULL; current=current->next){
			tmp->append(current->v, -1);
		}
		delete cluster[edge[min][0]];
		delete cluster[edge[min][1]];
		cluster[edge[min][0]] = tmp;
		cluster[edge[min][1]] = new Node;
		//////////////////comment end
		*/

		current_N--;
		if(current_N <= until)
			break;

		//upate tag infomarion
		/*
		if(degree[edge[min][0]] == 0 && degree[edge[min][1]] == 0)
			continue;

		tmp = new Node;
		c = 0;
		degree[edge[min][0]] = 0;
		degree[edge[min][1]] = 0;
		for(current = tag_edge[edge[min][0]]; current->next!=NULL; current=current->next){
			if(current->v != edge[min][1]){
				tmp->append(current->v, -1);
				degree[edge[min][0]]++;
				c++;
			}
		}
		for(current = tag_edge[edge[min][1]]; current->next!=NULL; current=current->next){
			if(current->v == edge[min][0])
				continue;
			i = 0;
			BOOL flag = true;
			for(Node* search = tmp; i<c; search=search->next){
				if(current->v == search->v){
					flag = false;
					break;
				}
				i++;
			}
			if(flag){
				tmp->append(current->v, -1);
				degree[edge[min][0]]++;
				for(Node* search = tag_edge[current->v]; search->next!=NULL; search=search->next){
					if(search->v == edge[min][1]){
						search->v = edge[min][0];
						break;
					}
				}
			}
			else{
				flag = false;
				for(Node* search = tag_edge[current->v]; search->next!=NULL; search=search->next){
					if(search->v == edge[min][0]){
						flag = true;
						break;
					}
				}
				for(search = tag_edge[current->v]; search->next!=NULL; search=search->next){
					if(search->v == edge[min][1]){
						if(flag){
							Node *tmp2 = new Node;
							for(Node* current2 = tag_edge[current->v]; current2->next!=NULL; current2=current2->next){
								if(current2->v != edge[min][1])
									tmp2->append(current2->v, -1);
							}
							delete tag_edge[current->v];
							tag_edge[current->v] = tmp2;
							degree[current->v]--;
						}
						else
							search->v = edge[min][0];
						break;
					}
				}
			}
		}
		delete tag_edge[edge[min][0]];
		tag_edge[edge[min][0]] = tmp;
		delete tag_edge[edge[min][1]];
		tag_edge[edge[min][1]] = new Node;
		*/
	}

	for(i=0; i<vertex_N; i++)
		if(e_table[i] != NULL)
			delete[] e_table[i];
	delete[] e_table;
	delete[] e_table_N;
	delete[] dist;
	delete[] heap;
	delete[] index;
	delete[] edge;
}

void Simplifier::Garland(float rate, BOOL mesure_err)
{
	float (*vertex)[3] = mesh->vertex;
	float (*normal)[3] = mesh->normal_f;
	int (*face)[3] = mesh->face;
	int** v_link = mesh->vertex_link_v;
	int** f_link = mesh->vertex_link_f;
	int vertex_N = mesh->vertex_N;
	int face_N = mesh->face_N;
	int** e_table = new int*[vertex_N];
	int* e_table_N = new int[vertex_N];
	int* degree_v = mesh->degree_v;
	int* degree_f = mesh->degree_f;
	float T = 0;
	int* isBound = mesh->isBound;
	int i;

	int v_N = mesh->countValidVertex();
	int until = v_N*(1.0f - rate);
	int current_N = v_N;

	//quadlic computation
	double (*Qf)[10] = new double[face_N][10];
	for(i=0; i<face_N; i++){
		double n[3];
		n[0] = normal[i][0];
		n[1] = normal[i][1];
		n[2] = normal[i][2];
		double area = mesh->faceArea(i);
		double d = -MeshData::DOT(n, vertex[face[i][0]]);
		MATRIX(Qf[i], n, d);
		MAT_TIMES(Qf[i], area);
	}

	double (*Q)[10] = new double[vertex_N][10];
	for(i=0; i<vertex_N; i++){
		int fn = degree_f[i];
		int* link = f_link[i];
		MAT_INIT(Q[i]);
		double total_area = 0;
		for(int j=0; j<fn; j++){
			MAT_SUM(Q[i], Qf[link[j]]);
			total_area += mesh->faceArea(link[j]);
		}
		if(total_area != 0)
			MAT_TIMES(Q[i], 1.0/total_area);
		if(isBound[i] == BOUNDARY){
			double dummy[10];
			float v1[3], v2[3];
            v1[0] = mesh->normal[i][0];
			v1[1] = mesh->normal[i][1];
			v1[2] = mesh->normal[i][2];
            MeshData::VEC(v2, vertex[v_link[i][0]], vertex[v_link[i][degree_v[i]-1]]);

			//MeshData::VEC(v1, vertex[i], vertex[v_link[i][0]]);
			//MeshData::VEC(v2, vertex[i], vertex[v_link[i][degree_v[i]-1]]);
			double n[3];
			MeshData::CROSS(n, v1, v2);
			double len = MeshData::LENGTH(n);
			if((float)len < 0.000001){
				n[0] = vertex[link[1]][0];
				n[1] = vertex[link[1]][1];
				n[2] = vertex[link[1]][2];
				len = MeshData::LENGTH(n);
			}
			n[0] /= len;
			n[1] /= len;
			n[2] /= len;

			double d = -MeshData::DOT(n, vertex[i]);
			MATRIX(dummy, n, d);
			MAT_TIMES(dummy, 10.0);
			MAT_SUM(Q[i], dummy);
		}
	}
	delete[] Qf;


	int edge_N = 0;
	for(i=0; i<vertex_N; i++){
		e_table[i] = new int[degree_v[i]];
		e_table_N[i] = 0;
		for(int j=0; j<degree_v[i]; j++)
			if(i > v_link[i][j] 
				&& (isBound[i] != NON_MANIFOLD && isBound[v_link[i][j]] != NON_MANIFOLD))
				edge_N++;
	}

	int (*edge)[2] = new int[edge_N][2];
	float* err = new float[edge_N];
	float (*opt_v)[3]  = new float[edge_N][3];

	int counter = 0;
	for(i=0; i<vertex_N; i++){
		for(int j=0; j<degree_v[i]; j++)
			if(i > v_link[i][j]
				&& (isBound[i] != NON_MANIFOLD && isBound[v_link[i][j]] != NON_MANIFOLD)){
				edge[counter][0] = i;
				edge[counter][1] = v_link[i][j];
				counter++;
			}
	}

	int* heap = new int[edge_N+1];
	int* index = new int[edge_N];
	int last_heap = 0;

	for(i=0; i<edge_N; i++){
		int i1 = edge[i][0];
		int i2 = edge[i][1];
		e_table[i1][e_table_N[i1]++] = i;
		e_table[i2][e_table_N[i2]++] = i;

		double v_min[3];
		double inv[6], Q12[10];
		MAT_PLUS(Q12, Q[i1], Q[i2]);

		if(INVERSE(inv, Q12)){
			double b[3];
			b[0] = -Q12[6];
			b[1] = -Q12[7];
			b[2] = -Q12[8];
			MAT_BY_VEC(v_min, inv, b);
		}
		else{
			v_min[0] = 0.5*(vertex[i1][0] + vertex[i2][0]);
			v_min[1] = 0.5*(vertex[i1][1] + vertex[i2][1]);
			v_min[2] = 0.5*(vertex[i1][2] + vertex[i2][2]);
			err[i] = Q_ERR(Q12, v_min);

			double v_tmp[3];
			v_tmp[0] = vertex[i1][0];
			v_tmp[1] = vertex[i1][1];
			v_tmp[2] = vertex[i1][2];
			double tmp_err = Q_ERR(Q12, v_tmp);
			if(err[i] > tmp_err){
				err[i] = tmp_err;
				v_min[0] = v_tmp[0];
				v_min[1] = v_tmp[1];
				v_min[2] = v_tmp[2];
			}

			v_tmp[0] = vertex[i2][0];
			v_tmp[1] = vertex[i2][1];
			v_tmp[2] = vertex[i2][2];
			tmp_err = Q_ERR(Q12, v_tmp);
			if(err[i] > tmp_err){
				err[i] = tmp_err;
				v_min[0] = v_tmp[0];
				v_min[1] = v_tmp[1];
				v_min[2] = v_tmp[2];
			}
		}

		opt_v[i][0] = v_min[0];
		opt_v[i][1] = v_min[1];
		opt_v[i][2] = v_min[2];

		err[i] = Q_ERR(Q12, v_min);

		//insert;
		heap[++last_heap] = i;
		index[i] = last_heap;
		upheap(err, last_heap, last_heap, heap, index);
	}

	while(last_heap != 0){
		int min = heap[1];
		index[min] = -1;
		//remove
		heap[1] = heap[last_heap--];
		index[heap[1]] = 1;
		downheap(err, last_heap, 1, heap, index);

		int i1 = edge[min][0];
		int i2 = edge[min][1];

		if(i1 < 0 || i2 < 0)
			continue;

		//if(opt_v[min][0] < 0.5)
			//continue;

		//half edge collapse
		int i3, i4;
		if(!contractEdge(i1, i2, opt_v[min], i3, i4))
			continue;

		err[min] = -1;

		MAT_SUM(Q[i1], Q[i2]);

		edge[min][0] = edge[min][1] = -1;

		//update edge infomation
		int* new_table;
		int c;
		int e3 = -1;
		int e4 = -1;
		
		//for i3
		if(i3 >= 0){
			new_table = new int[degree_v[i3]];
			c = 0;
			for(i=0; i<e_table_N[i3] ;i++){
				int e = e_table[i3][i];
				if(edge[e][0] == i2 || edge[e][1] == i2){
					e3 = e;
					edge[e][0] = -1;
					edge[e][1] = -1;
					//delete from heap
					if(index[e] > 0){
						err[e] = -1;
						upheap(err, last_heap, index[e], heap, index);
						index[e] = -1;
						heap[1] = heap[last_heap--];
						index[heap[1]] = 1;
						downheap(err, last_heap, 1, heap, index);
					}
				}
				else
					new_table[c++] = e;
			}
			e_table_N[i3] = c;
			delete[] e_table[i3];
			e_table[i3] = new_table;
		}
		//for i4
		if(i4 >= 0){
			new_table = new int[degree_v[i4]];
			c = 0;
			for(i=0; i<e_table_N[i4] ;i++){
				int e = e_table[i4][i];
				if(edge[e][0] == i2 || edge[e][1] == i2){
					e4 = e;
					edge[e][0] = -1;
					edge[e][1] = -1;
					//delete from heap
					if(index[e] > 0){
						err[e] = -1;
						upheap(err, last_heap, index[e], heap, index);
						index[e] = -1;
						heap[1] = heap[last_heap--];
						index[heap[1]] = 1;
						downheap(err, last_heap, 1, heap, index);
					}
				}
				else
					new_table[c++] = e;
			}
			e_table_N[i4] = c;
			delete[] e_table[i4];
			e_table[i4] = new_table;
		}

		//for i1
		c = 0;
		new_table = new int[degree_v[i1]];
		for(i=0; i<e_table_N[i1]; i++){
			int e = e_table[i1][i];
			if(e != min)
				new_table[c++] = e;
		}
		for(i=0; i<e_table_N[i2]; i++){
			int e = e_table[i2][i];
			if(e != min && e != e3 && e != e4){
				new_table[c++] = e;
				if(edge[e][0] == i2)
					edge[e][0] = i1;
				else
					edge[e][1] = i1;
			}
		}
		e_table_N[i1] = c;
		delete[] e_table[i1];
		e_table[i1] = new_table;

		//for i2 (delete)
		e_table_N[i2] = 0;
		delete[] e_table[i2]; 
		e_table[i2] = NULL;

		//update heap 
		for(i=0; i<e_table_N[i1]; i++){
			int e = e_table[i1][i];

			int ii1 = edge[e][0];
			int ii2 = edge[e][1];
			
			double v_min[3];
			double inv[6], Q12[10];
			MAT_PLUS(Q12, Q[ii1], Q[ii2]);

			if(INVERSE(inv, Q12)){
				double b[3];
				b[0] = -Q12[6];
				b[1] = -Q12[7];
				b[2] = -Q12[8];
				MAT_BY_VEC(v_min, inv, b);
			}
			else{
				v_min[0] = 0.5*(vertex[ii1][0] + vertex[ii2][0]);
				v_min[1] = 0.5*(vertex[ii1][1] + vertex[ii2][1]);
				v_min[2] = 0.5*(vertex[ii1][2] + vertex[ii2][2]);
				err[e] = Q_ERR(Q12, v_min);

				double v_tmp[3];
				v_tmp[0] = vertex[ii1][0];
				v_tmp[1] = vertex[ii1][1];
				v_tmp[2] = vertex[ii1][2];
				double tmp_err = Q_ERR(Q12, v_tmp);
				if(err[e] > tmp_err){
					err[e] = tmp_err;
					v_min[0] = v_tmp[0];
					v_min[1] = v_tmp[1];
					v_min[2] = v_tmp[2];
				}

				v_tmp[0] = vertex[ii2][0];
				v_tmp[1] = vertex[ii2][1];
				v_tmp[2] = vertex[ii2][2];
				tmp_err = Q_ERR(Q12, v_tmp);
				if(err[e] > tmp_err){
					err[e] = tmp_err;
					v_min[0] = v_tmp[0];
					v_min[1] = v_tmp[1];
					v_min[2] = v_tmp[2];
				}
			}
			opt_v[e][0] = v_min[0];  
			opt_v[e][1] = v_min[1]; 
			opt_v[e][2] = v_min[2];

			err[e] = Q_ERR(Q12, v_min);

			if(index[e] < 0){
				//insert
				heap[++last_heap] = e;
				index[e] = last_heap;
				upheap(err, last_heap, last_heap, heap, index);
			}
			else{
				//change;
				if(index[e] != 1 && err[e] < err[heap[index[e]/2]])
					upheap(err, last_heap, index[e], heap, index);
				else
					downheap(err, last_heap, index[e], heap, index);
			}
		}

		current_N--;
		if(current_N <= until)
			break;	
	}

	for(i=0; i<vertex_N; i++)
		if(e_table[i] != NULL)
			delete[] e_table[i];
	delete[] e_table;
	delete[] e_table_N;
	delete[] err;
	delete[] heap;
	delete[] index;
	delete[] edge;
	delete[] opt_v;
	delete[] Q;
}

BOOL Simplifier::contractEdge(int i1, int i2, float new_vertex[], int &i3, int &i4)
{
	float (*vertex)[3] = mesh->vertex;
	int **link_v = mesh->vertex_link_v;
	int **link_f = mesh->vertex_link_f;
	int *degree_v = mesh->degree_v;
	int *degree_f = mesh->degree_f;
	int (*face)[3] = mesh->face;
	float (*normal)[3] = mesh->normal_f;
	int *isBound = mesh->isBound;

	if(isBound[i1] == NON_MANIFOLD || isBound[i2] == NON_MANIFOLD)
		return false;

	//adjacent data for i1
	int *nei1, nei1_N, *face1, nei1f_N;
	mesh->getConsistent1Ring(i1, nei1, face1, nei1_N);
	nei1f_N = degree_f[i1];
	
	//adjacent data for i2
	int *nei2, nei2_N, *face2, nei2f_N;
	mesh->getConsistent1Ring(i2, nei2, face2, nei2_N);
	nei2f_N = degree_f[i2];

	//degree check
	if(isBound[i1] == BOUNDARY && isBound[i2] == BOUNDARY){
		if(nei1[0] == i2 || nei1[nei1_N-1] == i2 
			|| nei2[0] == i1 || nei2[nei2_N-1] == i1){
			if(nei1_N +  nei2_N < 6){
				return false;
			}
		}
		else if(nei1_N +  nei2_N < 7){
			return false;
		}
	}
	else if(nei1_N +  nei2_N < 7 || nei1_N +  nei2_N > 17){
		return false;
	}

	//opposite vertcies to the contracted edge;
	i3 = -1;
	i4 = -1;
	//degenerated face index;
	int f1 = -1;
	int f2 = -1;

	for(int i=0; i<nei1f_N; i++){
		if(nei1[i] == i2){
			i3 = nei1[(i+1)%nei1_N];
			f1 = face1[i];
		}
	}
	if(i3 >= 0 && mesh->isBound[i3] == NON_MANIFOLD)
		return false;

	for(i=0; i<nei2f_N; i++){
		if(nei2[i] == i1){
			i4 = nei2[(i+1)%nei2_N];
			f2 = face2[i];
		}
	}
	if(i4 >= 0 && mesh->isBound[i4] == NON_MANIFOLD)
		return false;

	if(i3 >= 0 && i4 >=0 && i3 == i4)
		return false;

	for(i=0; i<nei1_N; i++){
		if(nei1[i] == i2 || nei1[i] == i3 || nei1[i] == i4)
			continue;
		for(int j=0; j<nei2_N; j++)
			if(nei2[j] != i1 && nei1[i] == nei2[j])
				return false;
	}

	int *nei3, nei3_N, *face3, nei3f_N;
	if(i3 >= 0){
		mesh->getConsistent1Ring(i3, nei3, face3, nei3_N);
		nei3f_N = degree_f[i3];
	}

	int *nei4, nei4_N, *face4, nei4f_N;
	if(i4 >= 0){
		mesh->getConsistent1Ring(i4, nei4, face4,  nei4_N);
		nei4f_N = degree_f[i4];
	}

	if((i3 >= 0 && nei3_N < 4) || (i4 >= 0 && nei4_N < 4))
		return false;

	//flip check by "zone" technique
	BOOL flip = false;
	for(i=0; i<nei1f_N; i++){
		if(face1[i] == f1 || face1[i] == f2)
			continue;
		if(SHAPE(new_vertex, vertex[nei1[i]], vertex[nei1[(i+1)%nei1_N]]) < 0.000001){
			return false;
		}
		float v[3];
		MeshData::VEC(v, vertex[nei1[i]], vertex[nei1[(i+1)%nei1_N]]);
		double n[3];
		MeshData::CROSS(n, v, normal[face1[i]]);
		double len = MeshData::LENGTH(n);
		if(len == 0)
			continue;
		n[0] /= len;
		n[1] /= len;
		n[2] /= len;
		MeshData::VEC(v, vertex[i1], vertex[nei1[i]]);
		double side1 = MeshData::DOT(n, v);
		MeshData::VEC(v, new_vertex, vertex[nei1[i]]);
		double side2 = MeshData::DOT(n, v);
		if(side1 * side2 <= 0){
			return false;
		}
	}
	for(i=0; i<nei2f_N; i++){
		if(face2[i] == f1 || face2[i] == f2)
			continue;
		if(SHAPE(new_vertex, vertex[nei2[i]], vertex[nei2[(i+1)%nei2_N]]) < 0.000001){
			return false;
		}
		float v[3];
		MeshData::VEC(v, vertex[nei2[i]], vertex[nei2[(i+1)%nei2_N]]);
		double n[3];
		MeshData::CROSS(n, v, normal[face2[i]]);
		double len = MeshData::LENGTH(n);
		if(len == 0)
			continue;
		n[0] /= len;
		n[1] /= len;
		n[2] /= len;
		MeshData::VEC(v, vertex[i2], vertex[nei2[i]]);
		double side1 = MeshData::DOT(n, v);
		MeshData::VEC(v, new_vertex, vertex[nei2[i]]);
		double side2 = MeshData::DOT(n, v);
		if(side1 * side2 <= 0){
			return false;
		}
	}

	//update vertex position;
	vertex[i1][0] = new_vertex[0];
	vertex[i1][1] = new_vertex[1];
	vertex[i1][2] = new_vertex[2];

	//delete face
	if(f1 >= 0)
		face[f1][0] = face[f1][1] = face[f1][2] = -1;
	if(f2 >= 0)
		face[f2][0] = face[f2][1] = face[f2][2] = -1;

	//Updated adjacent vertex data;
	//for vertex i1 and i2;
	int link_N = nei1f_N + nei2f_N - 4;
	if(f1 < 0)
		link_N += 2;
	else if(f2 < 0)
		link_N += 2;
	int* link = new int[3*link_N];
	link_N = 0;
	for(i=0; i<nei1f_N; i++){
		if(face1[i] != f1 && face1[i] != f2){
			link[3*link_N] = nei1[i];
			link[3*link_N+1] = nei1[(i+1)%nei1_N];
			link[3*link_N+2] = face1[i];
			link_N++;
		}
	}
	for(i=0; i<nei2f_N; i++){
		if(face2[i] != f1 && face2[i] != f2){
			link[3*link_N] = nei2[i];
			link[3*link_N+1] = nei2[(i+1)%nei2_N];
			link[3*link_N+2] = face2[i];

			if(face[face2[i]][0] == i2)
				face[face2[i]][0] = i1;
			else if(face[face2[i]][1] == i2)
				face[face2[i]][1] = i1;
			else 
				face[face2[i]][2] = i1;

			link_N++;
		}
		else
			face[face2[i]][0] = face[face2[i]][1] = face[face2[i]][2] = -1;
	}
	delete[] link_v[i1];
	delete[] link_f[i1];
	isBound[i1] = mesh->sortVertexLink(link, link_N, link_v[i1], link_f[i1], degree_v[i1], degree_f[i1]);
	delete[] link;

	//upate link infomation incident to i2
	for(i=0; i<nei2_N; i++){
		if(nei2[i] == i1)
			continue;
		int* l = link_v[nei2[i]];
		int N = degree_v[nei2[i]];
		for(int j=0; j<N; j++)
			if(l[j] == i2)
				l[j] = i1;
	}
	delete[] link_v[i2];
	link_v[i2] = NULL;
	delete[] link_f[i2];
	link_f[i2] = NULL;
	degree_v[i2] = 0;
	degree_f[i2] = 0;
	isBound[i2] = NON_MANIFOLD;

	//upadate face normal incident to i1
	int* new_nei = link_v[i1];
	int* new_face = link_f[i1];
	int new_nei_N = degree_v[i1];
	int new_neif_N = degree_f[i1];
	for(i=0; i<new_neif_N; i++){
		float v1[3], v2[3];
		MeshData::VEC(v1, vertex[i1], vertex[new_nei[i]]);
		MeshData::VEC(v2, vertex[i1], vertex[new_nei[(i+1)%new_nei_N]]);
		double n[3];
		MeshData::CROSS(n, v1, v2);
		double len = MeshData::LENGTH(n);

		int f = new_face[i];
		if((float)len != 0){
			normal[f][0] = (float)(n[0]/len);
			normal[f][1] = (float)(n[1]/len);
			normal[f][2] = (float)(n[2]/len);
		}
		else
			normal[f][0] = normal[f][1] = normal[f][2] = 0;
	}

	//for vertex i3;
	if(i3 >= 0){
		link_N = nei3f_N - 1;
		link = new int[3*link_N];
		link_N = 0;
		for(i=0; i<nei3f_N; i++){
			if(face3[i] != f1){
				link[3*link_N] = nei3[i];
				link[3*link_N+1] = nei3[(i+1)%nei3_N];
				link[3*link_N+2] = face3[i];
				link_N++;
			}
		}
		delete[] link_v[i3];
		delete[] link_f[i3];
		isBound[i3] = mesh->sortVertexLink(link, link_N, link_v[i3], link_f[i3], 
				degree_v[i3], degree_f[i3]);
		delete[] link;
	}

	//for vertex i4
	if(i4 >= 0){
		link_N = nei4f_N - 1;
		link = new int[3*link_N];
		link_N = 0;
		for(i=0; i<nei4f_N; i++){
			if(face4[i] != f2){
				link[3*link_N] = nei4[i];
				link[3*link_N+1] = nei4[(i+1)%nei4_N];
				link[3*link_N+2] = face4[i];
				link_N++;
			}
		}
		delete[] link_v[i4];
		delete[] link_f[i4];
		isBound[i4] = mesh->sortVertexLink(link, link_N, link_v[i4], link_f[i4], 
				degree_v[i4], degree_f[i4]);
		delete[] link;
	}

	return true;
}

void Simplifier::GralandWithCluster(float rate)
{
	float (*vertex)[3] = mesh->vertex;
	float (*normal)[3] = mesh->normal_f;
	int (*face)[3] = mesh->face;
	int** v_link = mesh->vertex_link_v;
	int** f_link = mesh->vertex_link_f;
	int vertex_N = mesh->vertex_N;
	int face_N = mesh->face_N;
	int** e_table = new int*[vertex_N];
	int* e_table_N = new int[vertex_N];
	int* degree_v = mesh->degree_v;
	int* degree_f = mesh->degree_f;
	float T = 0;
	int* isBound = mesh->isBound;
	int i;

	int v_N = mesh->countValidVertex();
	int until = v_N*(1.0f - rate);
	int current_N = v_N;

	//quadlic computation
	double (*Qf)[10] = new double[face_N][10];
	for(i=0; i<face_N; i++){
		double n[3];
		n[0] = normal[i][0];
		n[1] = normal[i][1];
		n[2] = normal[i][2];
		double area = mesh->faceArea(i);
		double d = -MeshData::DOT(n, vertex[face[i][0]]);
		MATRIX(Qf[i], n, d);
		MAT_TIMES(Qf[i], area);
	}

	double (*Q)[10] = new double[vertex_N][10];
	for(i=0; i<vertex_N; i++){
		int fn = degree_f[i];
		int* link = f_link[i];
		MAT_INIT(Q[i]);
		double total_area = 0;
		for(int j=0; j<fn; j++){
			MAT_SUM(Q[i], Qf[link[j]]);
			total_area += mesh->faceArea(link[j]);
		}
		if(total_area != 0)
			MAT_TIMES(Q[i], 1.0/total_area);
		if(isBound[i] == BOUNDARY){
			double dummy[10];
			float v1[3], v2[3];
			MeshData::VEC(v1, vertex[i], vertex[v_link[i][0]]);
			MeshData::VEC(v2, vertex[i], vertex[v_link[i][degree_v[i]-1]]);
			double n[3];
			MeshData::CROSS(n, v1, v2);
			double len = MeshData::LENGTH(n);
			if((float)len < 0.000001){
				n[0] = vertex[link[1]][0];
				n[1] = vertex[link[1]][1];
				n[2] = vertex[link[1]][2];
				len = MeshData::LENGTH(n);
			}
			n[0] /= len;
			n[1] /= len;
			n[2] /= len;
			double d = -MeshData::DOT(n, vertex[i]);
			MATRIX(dummy, n, d);
			MAT_TIMES(dummy, 10.0);
			MAT_SUM(Q[i], dummy);
		}
	}
	delete[] Qf;


	int edge_N = 0;
	for(i=0; i<vertex_N; i++){
		e_table[i] = new int[degree_v[i]];
		e_table_N[i] = 0;
		for(int j=0; j<degree_v[i]; j++)
			if(i > v_link[i][j] 
				&& (isBound[i] != NON_MANIFOLD && isBound[v_link[i][j]] != NON_MANIFOLD))
				edge_N++;
	}

	int (*edge)[2] = new int[edge_N][2];
	float* err = new float[edge_N];
	float (*opt_v)[3]  = new float[edge_N][3];

	int counter = 0;
	for(i=0; i<vertex_N; i++){
		for(int j=0; j<degree_v[i]; j++)
			if(i > v_link[i][j]
				&& (isBound[i] != NON_MANIFOLD && isBound[v_link[i][j]] != NON_MANIFOLD)){
				edge[counter][0] = i;
				edge[counter][1] = v_link[i][j];
				counter++;
			}
	}

	int* heap = new int[edge_N+1];
	int* index = new int[edge_N];
	int last_heap = 0;

	for(i=0; i<edge_N; i++){
		int i1 = edge[i][0];
		int i2 = edge[i][1];
		e_table[i1][e_table_N[i1]++] = i;
		e_table[i2][e_table_N[i2]++] = i;

		double v_min[3];
		double inv[6], Q12[10];
		MAT_PLUS(Q12, Q[i1], Q[i2]);

		if(INVERSE(inv, Q12)){
			double b[3];
			b[0] = -Q12[6];
			b[1] = -Q12[7];
			b[2] = -Q12[8];
			MAT_BY_VEC(v_min, inv, b);
		}
		else{
			v_min[0] = 0.5*(vertex[i1][0] + vertex[i2][0]);
			v_min[1] = 0.5*(vertex[i1][1] + vertex[i2][1]);
			v_min[2] = 0.5*(vertex[i1][2] + vertex[i2][2]);
			err[i] = Q_ERR(Q12, v_min);

			double v_tmp[3];
			v_tmp[0] = vertex[i1][0];
			v_tmp[1] = vertex[i1][1];
			v_tmp[2] = vertex[i1][2];
			double tmp_err = Q_ERR(Q12, v_tmp);
			if(err[i] > tmp_err){
				err[i] = tmp_err;
				v_min[0] = v_tmp[0];
				v_min[1] = v_tmp[1];
				v_min[2] = v_tmp[2];
			}

			v_tmp[0] = vertex[i2][0];
			v_tmp[1] = vertex[i2][1];
			v_tmp[2] = vertex[i2][2];
			tmp_err = Q_ERR(Q12, v_tmp);
			if(err[i] > tmp_err){
				err[i] = tmp_err;
				v_min[0] = v_tmp[0];
				v_min[1] = v_tmp[1];
				v_min[2] = v_tmp[2];
			}
		}

		opt_v[i][0] = v_min[0];
		opt_v[i][1] = v_min[1];
		opt_v[i][2] = v_min[2];

		err[i] = Q_ERR(Q12, v_min);

		//insert;
		heap[++last_heap] = i;
		index[i] = last_heap;
		upheap(err, last_heap, last_heap, heap, index);
	}

	Node **cluster;
	ridge_cluster = (Node**)new Node[vertex_N];
	for(i=0; i<vertex_N; i++){
		ridge_cluster[i] = new Node;
		ridge_cluster[i]->append(i, -1);
	}
	cluster = ridge_cluster;

	while(last_heap != 0){
		int min = heap[1];
		index[min] = -1;
		//remove
		heap[1] = heap[last_heap--];
		index[heap[1]] = 1;
		downheap(err, last_heap, 1, heap, index);

		int i1 = edge[min][0];
		int i2 = edge[min][1];

		if(i1 < 0 || i2 < 0)
			continue;

		//half edge collapse
		int i3, i4;
		if(!contractEdge(i1, i2, opt_v[min], i3, i4))
			continue;

		err[min] = -1;

		MAT_SUM(Q[i1], Q[i2]);

		edge[min][0] = edge[min][1] = -1;

		//update edge infomation
		int* new_table;
		int c;
		int e3 = -1;
		int e4 = -1;
		
		//for i3
		if(i3 >= 0){
			new_table = new int[degree_v[i3]];
			c = 0;
			for(i=0; i<e_table_N[i3] ;i++){
				int e = e_table[i3][i];
				if(edge[e][0] == i2 || edge[e][1] == i2){
					e3 = e;
					edge[e][0] = -1;
					edge[e][1] = -1;
					//delete from heap
					if(index[e] > 0){
						err[e] = -1;
						upheap(err, last_heap, index[e], heap, index);
						index[e] = -1;
						heap[1] = heap[last_heap--];
						index[heap[1]] = 1;
						downheap(err, last_heap, 1, heap, index);
					}
				}
				else
					new_table[c++] = e;
			}
			e_table_N[i3] = c;
			delete[] e_table[i3];
			e_table[i3] = new_table;
		}
		//for i4
		if(i4 >= 0){
			new_table = new int[degree_v[i4]];
			c = 0;
			for(i=0; i<e_table_N[i4] ;i++){
				int e = e_table[i4][i];
				if(edge[e][0] == i2 || edge[e][1] == i2){
					e4 = e;
					edge[e][0] = -1;
					edge[e][1] = -1;
					//delete from heap
					if(index[e] > 0){
						err[e] = -1;
						upheap(err, last_heap, index[e], heap, index);
						index[e] = -1;
						heap[1] = heap[last_heap--];
						index[heap[1]] = 1;
						downheap(err, last_heap, 1, heap, index);
					}
				}
				else
					new_table[c++] = e;
			}
			e_table_N[i4] = c;
			delete[] e_table[i4];
			e_table[i4] = new_table;
		}

		//for i1
		c = 0;
		new_table = new int[degree_v[i1]];
		for(i=0; i<e_table_N[i1]; i++){
			int e = e_table[i1][i];
			if(e != min)
				new_table[c++] = e;
		}
		for(i=0; i<e_table_N[i2]; i++){
			int e = e_table[i2][i];
			if(e != min && e != e3 && e != e4){
				new_table[c++] = e;
				if(edge[e][0] == i2)
					edge[e][0] = i1;
				else
					edge[e][1] = i1;
			}
		}
		e_table_N[i1] = c;
		delete[] e_table[i1];
		e_table[i1] = new_table;

		//for i2 (delete)
		e_table_N[i2] = 0;
		delete[] e_table[i2]; 
		e_table[i2] = NULL;

		//update heap 
		for(i=0; i<e_table_N[i1]; i++){
			int e = e_table[i1][i];

			int ii1 = edge[e][0];
			int ii2 = edge[e][1];
			
			double v_min[3];
			double inv[6], Q12[10];
			MAT_PLUS(Q12, Q[ii1], Q[ii2]);

			if(INVERSE(inv, Q12)){
				double b[3];
				b[0] = -Q12[6];
				b[1] = -Q12[7];
				b[2] = -Q12[8];
				MAT_BY_VEC(v_min, inv, b);
			}
			else{
				v_min[0] = 0.5*(vertex[ii1][0] + vertex[ii2][0]);
				v_min[1] = 0.5*(vertex[ii1][1] + vertex[ii2][1]);
				v_min[2] = 0.5*(vertex[ii1][2] + vertex[ii2][2]);
				err[e] = Q_ERR(Q12, v_min);

				double v_tmp[3];
				v_tmp[0] = vertex[ii1][0];
				v_tmp[1] = vertex[ii1][1];
				v_tmp[2] = vertex[ii1][2];
				double tmp_err = Q_ERR(Q12, v_tmp);
				if(err[e] > tmp_err){
					err[e] = tmp_err;
					v_min[0] = v_tmp[0];
					v_min[1] = v_tmp[1];
					v_min[2] = v_tmp[2];
				}

				v_tmp[0] = vertex[ii2][0];
				v_tmp[1] = vertex[ii2][1];
				v_tmp[2] = vertex[ii2][2];
				tmp_err = Q_ERR(Q12, v_tmp);
				if(err[e] > tmp_err){
					err[e] = tmp_err;
					v_min[0] = v_tmp[0];
					v_min[1] = v_tmp[1];
					v_min[2] = v_tmp[2];
				}
			}
			opt_v[e][0] = v_min[0];  
			opt_v[e][1] = v_min[1]; 
			opt_v[e][2] = v_min[2];

			err[e] = Q_ERR(Q12, v_min);

			if(index[e] < 0){
				//insert
				heap[++last_heap] = e;
				index[e] = last_heap;
				upheap(err, last_heap, last_heap, heap, index);
			}
			else{
				//change;
				if(index[e] != 1 && err[e] < err[heap[index[e]/2]])
					upheap(err, last_heap, index[e], heap, index);
				else
					downheap(err, last_heap, index[e], heap, index);
			}
		}

		//update vertex cluster
		Node* tmp = new Node;
		for(Node* current=cluster[i1]; current->next!=NULL; current=current->next){
			tmp->append(current->v, -1);
		}
		for(current=cluster[i2]; current->next!=NULL; current=current->next){
			tmp->append(current->v, -1);
		}
		delete cluster[i1];
		delete cluster[i2];
		cluster[i1] = tmp;
		cluster[i2] = new Node;

		current_N--;
		if(current_N <= until)
			break;	
	}

	for(i=0; i<vertex_N; i++)
		if(e_table[i] != NULL)
			delete[] e_table[i];
	delete[] e_table;
	delete[] e_table_N;
	delete[] err;
	delete[] heap;
	delete[] index;
	delete[] edge;
	delete[] opt_v;
	delete[] Q;
}