// Subdivider.cpp: Subdivider クラスのインプリメンテーション
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "MeshEditor.h"
#include "Subdivider.h"
#include <math.h>

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

//////////////////////////////////////////////////////////////////////
// 構築/消滅
//////////////////////////////////////////////////////////////////////

Subdivider::Subdivider()
{
	vertex_index = NULL;
}

Subdivider::~Subdivider()
{
	if(vertex_index != NULL)
		delete[] vertex_index;
}

void Subdivider::setMeshData(MeshData *mesh)
{
	this->mesh = mesh;
}

void Subdivider::linear()
{
	int (*face)[3] = mesh->face;
	float (*vertex)[3] = mesh->vertex;
	int face_N = mesh->face_N;
	int vertex_N = mesh->vertex_N;
	int (*link)[3] = mesh->face_link_E;

	BOOL* flag = new BOOL[face_N];
	for(int i=0; i<face_N; i++){
		//if(vertex[face[i][0]][2] > 0)
			flag[i] = true;
		//else
			//flag[i] = false;
	}

	this->subdivideTopology(flag);

	for(i=0; i<vertex_N; i++){
		new_vertex[i][0] = vertex[i][0];
		new_vertex[i][1] = vertex[i][1];
		new_vertex[i][2] = vertex[i][2];
	}

	int counter= 0;
	for(i=0; i<face_N; i++){
		if(!flag[i])
			continue;
		for(int j=0; j<3; j++){
			int pair = link[i][j];
			if(pair < i || !flag[pair]){
				int index = vertex_index[i][j];
				int i1 = face[i][j];
				int i2 = face[i][(j+1)%3];
				new_vertex[index][0] = 0.5f*(vertex[i1][0] + vertex[i2][0]);
				new_vertex[index][1] = 0.5f*(vertex[i1][1] + vertex[i2][1]);
				new_vertex[index][2] = 0.5f*(vertex[i1][2] + vertex[i2][2]);		
			}
		}
	}

	delete[] flag;
	delete[] vertex_index;
	vertex_index = NULL;

	this->updateMesh();
}

void Subdivider::LOOP()
{
	int (*face)[3] = mesh->face;
	float (*vertex)[3] = mesh->vertex;
	int face_N = mesh->face_N;
	int vertex_N = mesh->vertex_N;
	int (*link_f)[3] = mesh->face_link_E;
	BOOL* isBound = mesh->isBound;

	BOOL* flag = new BOOL[face_N];
	for(int i=0; i<face_N; i++)
		flag[i] = true;

	this->subdivideTopology(flag);

	for(i=0; i<vertex_N; i++){
		if(isBound[i] == INNER){
			int *link, link_N;
			mesh->getConsistent1Ring(i, link, link_N);
			if(link_N == 0){
				new_vertex[i][0] = vertex[i][0];
				new_vertex[i][1] = vertex[i][1];
				new_vertex[i][2] = vertex[i][2];
				continue;
			}
			double b = (0.625 - (0.375 + 0.25*cos(2.0*PI/(double)link_N))
				*(0.375 + 0.25*cos(2.0*PI/(double)link_N)))/(double)link_N;
			double sum[3];
			sum[0] = sum[1] = sum[2] = 0;
			for(int j=0; j<link_N; j++){
				sum[0] += vertex[link[j]][0];
				sum[1] += vertex[link[j]][1];
				sum[2] += vertex[link[j]][2];
			}
			new_vertex[i][0] = (float)((1.0-(double)link_N*b)*vertex[i][0] + b*sum[0]);
			new_vertex[i][1] = (float)((1.0-(double)link_N*b)*vertex[i][1] + b*sum[1]);
			new_vertex[i][2] = (float)((1.0-(double)link_N*b)*vertex[i][2] + b*sum[2]);
		}
		else if(isBound[i] == BOUNDARY){
			int *link, link_N;
			mesh->getConsistent1Ring(i, link, link_N);
			int i1 = link[0];
			int i2 = link[link_N-1];
			new_vertex[i][0] = (float)(0.125*(vertex[i1][0] + vertex[i2][0] + 6.0*vertex[i][0]));
			new_vertex[i][1] = (float)(0.125*(vertex[i1][1] + vertex[i2][1] + 6.0*vertex[i][1]));
			new_vertex[i][2] = (float)(0.125*(vertex[i1][2] + vertex[i2][2] + 6.0*vertex[i][2]));
		}
		else{
			new_vertex[i][0] = vertex[i][0];
			new_vertex[i][1] = vertex[i][1];
			new_vertex[i][2] = vertex[i][2];
		}
	}

	for(i=0; i<face_N; i++){
		if(!flag[i])
			continue;
		for(int j=0; j<3; j++){
			int pair = link_f[i][j];
			if(pair < i || !flag[pair]){
				int index = vertex_index[i][j];
				int i1 = face[i][j];
				int i2 = face[i][(j+1)%3];

				if(pair < 0){
					new_vertex[index][0] = 0.5f*(vertex[i1][0] + vertex[i2][0]);
					new_vertex[index][1] = 0.5f*(vertex[i1][1] + vertex[i2][1]);
					new_vertex[index][2] = 0.5f*(vertex[i1][2] + vertex[i2][2]);
					continue;
				}
			
				int i3 = face[i][(j+2)%3];
				int i4;
				if(link_f[pair][0] == i)
					i4 = face[pair][2];
				else if(link_f[pair][1] == i)
					i4 = face[pair][0];
				else
					i4 = face[pair][1];
				if(isBound[i1] == INNER || isBound[i2] == INNER ||
					(i3 >= 0 && i4 >= 0)){
					new_vertex[index][0] = (float)(0.125*(3.0*(vertex[i1][0] + vertex[i2][0]) 
														+ vertex[i3][0] + vertex[i4][0]));
					new_vertex[index][1] = (float)(0.125*(3.0*(vertex[i1][1] + vertex[i2][1]) 
														+ vertex[i3][1] + vertex[i4][1]));
					new_vertex[index][2] = (float)(0.125*(3.0*(vertex[i1][2] + vertex[i2][2]) 
														+ vertex[i3][2] + vertex[i4][2]));
				}
				else{
					new_vertex[index][0] = 0.5f*(vertex[i1][0] + vertex[i2][0]);
					new_vertex[index][1] = 0.5f*(vertex[i1][1] + vertex[i2][1]);
					new_vertex[index][2] = 0.5f*(vertex[i1][2] + vertex[i2][2]);
				}
			}
		}
	}

	delete[] flag;
	delete[] vertex_index;
	vertex_index = NULL;
	
	this->updateMesh();

	
}

void Subdivider::subdivideTopology(BOOL *flag)
{
	int face_N = mesh->face_N;
	int vertex_N = mesh->vertex_N;
	int (*link)[3] = mesh->face_link_E;
	float (*vertex)[3] = mesh->vertex;
	int (*face)[3] = mesh->face;
	BOOL grow = true;
	while(grow){
		grow = false;
		for(int i=0; i<face_N; i++){
			if(flag[i])
				continue;
			int count = 0;
			for(int j=0; j<3; j++){
				if(flag[link[i][j]])
					count++;
			}
			if(count > 1){
				flag[i] = true;
				grow = true;
			}
		}
	}

	if(vertex_index != NULL)
		delete[] vertex_index;
	vertex_index = new int[face_N][3];
	for(int i=0; i<face_N; i++)
		vertex_index[i][0] = vertex_index[i][1] = vertex_index[i][2] = -1;

	new_vertex_N = vertex_N;
	for(i=0; i<face_N; i++){
		if(!flag[i])
			continue;
		for(int j=0; j<3; j++){
			int pair = link[i][j];
			if(pair < i || !flag[pair]){
				vertex_index[i][j] = new_vertex_N;
				if(pair >= 0){
					if(link[pair][0] == i)
						vertex_index[pair][0] = new_vertex_N;
					else if(link[pair][1] == i)
						vertex_index[pair][1] = new_vertex_N;
					else
						vertex_index[pair][2] = new_vertex_N;
				}
				new_vertex_N++;
			}
		}
	}

	new_vertex = new float[new_vertex_N][3];

	new_face_N = face_N;
	for(i=0; i<face_N; i++){
		if(flag[i])
			new_face_N += 3;
		else{
			int count = 0;
			for(int j=0; j<3; j++){
				if(flag[link[i][j]])
					count++;
			}
			if(count != 0)
				new_face_N++;
		}
	}
	new_face = new int[new_face_N][3];
	int face_id = 0;
	for(i=0; i<face_N; i++){
		if(flag[i]){
			new_face[face_id][0] = face[i][0];
			new_face[face_id][1] = vertex_index[i][0];
			new_face[face_id][2] = vertex_index[i][2];
			face_id++;

			new_face[face_id][0] = face[i][1];
			new_face[face_id][1] = vertex_index[i][1];
			new_face[face_id][2] = vertex_index[i][0];
			face_id++;

			new_face[face_id][0] = face[i][2];
			new_face[face_id][1] = vertex_index[i][2];
			new_face[face_id][2] = vertex_index[i][1];
			face_id++;

			new_face[face_id][0] = vertex_index[i][0];
			new_face[face_id][1] = vertex_index[i][1];
			new_face[face_id][2] = vertex_index[i][2];
			face_id++;
		}
		else{
			int id = -1;
			for(int j=0; j<3; j++){
				if(flag[link[i][j]])
					id = j;
			}
			switch(id){
			case -1:
				new_face[face_id][0] = face[i][0];
				new_face[face_id][1] = face[i][1];
				new_face[face_id][2] = face[i][2];
				face_id++;
				break;

			case 0:
				new_face[face_id][0] = face[i][0];
				new_face[face_id][1] = vertex_index[i][0];
				new_face[face_id][2] = face[i][2];
				face_id++;

				new_face[face_id][0] = face[i][1];
				new_face[face_id][1] = face[i][2];
				new_face[face_id][2] = vertex_index[i][0];
				face_id++;
				break;

			case 1:
				new_face[face_id][0] = face[i][0];
				new_face[face_id][1] = face[i][1];
				new_face[face_id][2] = vertex_index[i][1];
				face_id++;

				new_face[face_id][0] = face[i][2];
				new_face[face_id][1] = face[i][0];
				new_face[face_id][2] = vertex_index[i][1];
				face_id++;
				break;

			case 2:
				new_face[face_id][0] = face[i][0];
				new_face[face_id][1] = face[i][1];
				new_face[face_id][2] = vertex_index[i][2];
				face_id++;

				new_face[face_id][0] = face[i][1];
				new_face[face_id][1] = face[i][2];
				new_face[face_id][2] = vertex_index[i][2];
				face_id++;
				break;
			}	
		}
	}
}

void Subdivider::updateMesh()
{
	mesh->deleteLinkData();
	mesh->vertex_N = new_vertex_N;
	mesh->vertex = new_vertex;
	mesh->face_N = new_face_N;
	mesh->face = new_face;
	
	mesh->generateVertexLink();
	mesh->generateFaceLink(); 
	mesh->computeFaceNormal();
	mesh->computeNormal();
}

void Subdivider::butterfly()
{
	int (*face)[3] = mesh->face;
	float (*vertex)[3] = mesh->vertex;
	int face_N = mesh->face_N;
	int vertex_N = mesh->vertex_N;
	int (*link_f)[3] = mesh->face_link_E;
	BOOL* isBound = mesh->isBound;
	int* degree = mesh->degree_v;

	BOOL* flag = new BOOL[face_N];
	for(int i=0; i<face_N; i++){
		//if(vertex[face[i][0]][2] > 0)
			flag[i] = true;
		//else
			//flag[i] = false;
	}

	this->subdivideTopology(flag);

	for(i=0; i<vertex_N; i++){
		new_vertex[i][0] = vertex[i][0];
		new_vertex[i][1] = vertex[i][1];
		new_vertex[i][2] = vertex[i][2];
	}

	for(i=0; i<face_N; i++){
		if(!flag[i])
			continue;
		for(int j=0; j<3; j++){
			int pair = link_f[i][j];
			if(pair < i || !flag[pair]){
				int index = vertex_index[i][j];
				int i1 = face[i][j];
				int i2 = face[i][(j+1)%3];

				int *link1, *link2, link1_N, link2_N;
				mesh->getConsistent1Ring(i1, link1, link1_N);
				mesh->getConsistent1Ring(i2, link2, link2_N);

				if(pair < 0 && isBound[i1] == BOUNDARY && isBound[i2] == BOUNDARY){
					int i3 = link1[0];
					if(i3 == i2)
						i3 = link1[link1_N-1];
					int i4 = link2[0];
					if(i4 == i1)
						i4 = link2[link2_N-1];
					new_vertex[index][0] = (float)(0.0625*(9.0*(vertex[i1][0] + vertex[i2][0]) 
						- (vertex[i3][0] + vertex[i4][0])));
					new_vertex[index][1] = (float)(0.0625*(9.0*(vertex[i1][1] + vertex[i2][1]) 
						- (vertex[i3][1] + vertex[i4][1])));
					new_vertex[index][2] = (float)(0.0625*(9.0*(vertex[i1][2] + vertex[i2][2]) 
						- (vertex[i3][2] + vertex[i4][2])));
					continue;
				}

				int i3, i4;
				i3 = face[i][(j+2)%3];
				if(link_f[pair][0] == i)
					i4 = face[pair][2];
				else if(link_f[pair][1] == i)
					i4 = face[pair][0];
				else
					i4 = face[pair][1];

				if(isBound[i1] == INNER && isBound[i2] == INNER){
					if(degree[i1] == 6 && degree[i2] == 6){
						int index1, index2;
						for(int k=0; k<link1_N; k++)
							if(link1[k] == i2)
								index1 = k;
						for(k=0; k<link2_N; k++)
							if(link2[k] == i1)
								index2 = k;
						int i5, i6, i7, i8;
						i5 = link2[(index2+link1_N-2)%link2_N];
						i6 = link1[(index1+2)%link1_N];
						i7 = link1[(index1+link1_N-2)%link1_N];
						i8 = link2[(index2+2)%link2_N];

						new_vertex[index][0] = (float)(0.0625*(8.0*(vertex[i1][0] + vertex[i2][0]) 
							+ 2.0*(vertex[i3][0] + vertex[i4][0])
							- (vertex[i5][0] + vertex[i6][0] + vertex[i7][0] + vertex[i8][0])));
						new_vertex[index][1] = (float)(0.0625*(8.0*(vertex[i1][1] + vertex[i2][1]) 
							+ 2.0*(vertex[i3][1] + vertex[i4][1])
							- (vertex[i5][1] + vertex[i6][1] + vertex[i7][1] + vertex[i8][1])));
						new_vertex[index][2] = (float)(0.0625*(8.0*(vertex[i1][2] + vertex[i2][2]) 
							+ 2.0*(vertex[i3][2] + vertex[i4][2])
							- (vertex[i5][2] + vertex[i6][2] + vertex[i7][2] + vertex[i8][2])));
					}
					else{
						if(degree[i1] == 6 && degree[i2] != 6){
							int tmp = i1;
							i1 = i2;
							i2 = tmp;
							link1 = link2;
							link1_N = link2_N;
						}
						if(degree[i1] != 6 && degree[i2] == 6){
							int index1;
							for(int k=0; k<link1_N; k++){
								if(link1[k] == i2)
									index1 = k;
							}
							double sum[3];
							
							if(degree[i1] == 3){
								sum[0] = 0.25*(5.0*vertex[i2][0] 
									- (vertex[link1[(index1+1)%3]][0] + vertex[link1[(index1+2)%3]][0]));
								sum[1] = 0.25*(5.0*vertex[i2][1] 
									- (vertex[link1[(index1+1)%3]][1] + vertex[link1[(index1+2)%3]][1]));
								sum[2] = 0.25*(5.0*vertex[i2][2] 
									- (vertex[link1[(index1+1)%3]][2] + vertex[link1[(index1+2)%3]][2]));
							}
							else if(degree[i1] == 4){
								sum[0] = 0.5*(3.0*vertex[i2][0] - vertex[link1[(index1+2)%4]][0]);
								sum[1] = 0.5*(3.0*vertex[i2][1] - vertex[link1[(index1+2)%4]][1]);
								sum[2] = 0.5*(3.0*vertex[i2][2] - vertex[link1[(index1+2)%4]][2]);
							}
							else{
								sum[0] = sum[1] = sum[2] = 0;
								for(int l=0; l<link1_N; l++){
									int m = link1[(index1+l)%link1_N];
									double s = (0.25 + cos(2.0*PI*(double)l/(double)link1_N) 
										+ 0.5*cos(4.0*PI*(double)l/(double)link1_N));
									sum[0] += s*vertex[m][0];
									sum[1] += s*vertex[m][1];
									sum[2] += s*vertex[m][2];
								}
							}
							new_vertex[index][0] = (float)(0.75*vertex[i1][0] + sum[0]/(double)link1_N);
							new_vertex[index][1] = (float)(0.75*vertex[i1][1] + sum[1]/(double)link1_N);
							new_vertex[index][2] = (float)(0.75*vertex[i1][2] + sum[2]/(double)link1_N);
						}
						else{
							int index1;
							for(int k=0; k<link1_N; k++){
								if(link1[k] == i2)
									index1 = k;
							}
							double sum[3];
							
							if(degree[i1] == 3){
								sum[0] = 0.25*(5.0*vertex[i2][0] 
									- (vertex[link1[(index1+1)%3]][0] + vertex[link1[(index1+2)%3]][0]));
								sum[1] = 0.25*(5.0*vertex[i2][1] 
									- (vertex[link1[(index1+1)%3]][1] + vertex[link1[(index1+2)%3]][1]));
								sum[2] = 0.25*(5.0*vertex[i2][2] 
									- (vertex[link1[(index1+1)%3]][2] + vertex[link1[(index1+2)%3]][2]));
							}
							else if(degree[i1] == 4){
								sum[0] = 0.5*(3.0*vertex[i2][0] - vertex[link1[(index1+2)%4]][0]);
								sum[1] = 0.5*(3.0*vertex[i2][1] - vertex[link1[(index1+2)%4]][1]);
								sum[2] = 0.5*(3.0*vertex[i2][2] - vertex[link1[(index1+2)%4]][2]);
							}
							else{
								sum[0] = sum[1] = sum[2] = 0;
								for(int l=0; l<link1_N; l++){
									int m = link1[(index1+l)%link1_N];
									double s = (0.25 + cos(2.0*PI*(double)l/(double)link1_N) 
										+ 0.5*cos(4.0*PI*(double)l/(double)link1_N));
									sum[0] += s*vertex[m][0];
									sum[1] += s*vertex[m][1];
									sum[2] += s*vertex[m][2];
								}
							}
							double tmp1[3];
							tmp1[0] = 0.75*vertex[i1][0] + sum[0]/(double)link1_N;
							tmp1[1] = 0.75*vertex[i1][1] + sum[1]/(double)link1_N;
							tmp1[2] = 0.75*vertex[i1][2] + sum[2]/(double)link1_N;

							int tmp = i1;
							i1 = i2;
							i2 = tmp;
							link1 = link2;
							link1_N = link2_N;

							for(k=0; k<link1_N; k++){
								if(link1[k] == i2)
									index1 = k;
							}
							
							if(degree[i1] == 3){
								sum[0] = 0.25*(5.0*vertex[i2][0] 
									- (vertex[link1[(index1+1)%3]][0] + vertex[link1[(index1+2)%3]][0]));
								sum[1] = 0.25*(5.0*vertex[i2][1] 
									- (vertex[link1[(index1+1)%3]][1] + vertex[link1[(index1+2)%3]][1]));
								sum[2] = 0.25*(5.0*vertex[i2][2] 
									- (vertex[link1[(index1+1)%3]][2] + vertex[link1[(index1+2)%3]][2]));
							}
							else if(degree[i1] == 4){
								sum[0] = 0.5*(3.0*vertex[i2][0] - vertex[link1[(index1+2)%4]][0]);
								sum[1] = 0.5*(3.0*vertex[i2][1] - vertex[link1[(index1+2)%4]][1]);
								sum[2] = 0.5*(3.0*vertex[i2][2] - vertex[link1[(index1+2)%4]][2]);
							}
							else{
								sum[0] = sum[1] = sum[2] = 0;
								for(int l=0; l<link1_N; l++){
									int m = link1[(index1+l)%link1_N];
									double s = (0.25 + cos(2.0*PI*(double)l/(double)link1_N) 
										+ 0.5*cos(4.0*PI*(double)l/(double)link1_N));
									sum[0] += s*vertex[m][0];
									sum[1] += s*vertex[m][1];
									sum[2] += s*vertex[m][2];
								}
							}
							double tmp2[3];
							tmp2[0] = 0.75*vertex[i1][0] + sum[0]/(double)link1_N;
							tmp2[1] = 0.75*vertex[i1][1] + sum[1]/(double)link1_N;
							tmp2[2] = 0.75*vertex[i1][2] + sum[2]/(double)link1_N;

							new_vertex[index][0] = (float)(0.5*(tmp1[0] + tmp2[0]));
							new_vertex[index][1] = (float)(0.5*(tmp1[1] + tmp2[1]));
							new_vertex[index][2] = (float)(0.5*(tmp1[2] + tmp2[2]));
						}
					}
				}
				else{
					new_vertex[index][0] = 0.5f*(vertex[i1][0] + vertex[i2][0]);
					new_vertex[index][1] = 0.5f*(vertex[i1][1] + vertex[i2][1]);
					new_vertex[index][2] = 0.5f*(vertex[i1][2] + vertex[i2][2]);
				}
			}
		}
	}

	delete[] flag;
	delete[] vertex_index;
	vertex_index = NULL;
	
	this->updateMesh();

	
}

void Subdivider::normalAverage()
{
	
	float (*vertex)[3] = mesh->vertex;
	int face_N = mesh->face_N;
	int (*face)[3] = new int[face_N][3];
	int vertex_N = mesh->vertex_N;
	int (*link)[3] = new int[face_N][3];
	double (*L)[3] = new double[vertex_N][3];

	for(int i=0; i<vertex_N; i++)
		mesh->laplacian(i, L[i]);

	for(i=0; i<face_N; i++){
		link[i][0] = mesh->face_link_E[i][0];
		link[i][1] = mesh->face_link_E[i][1];
		link[i][2] = mesh->face_link_E[i][2];

		face[i][0] = mesh->face[i][0];
		face[i][1] = mesh->face[i][1];
		face[i][2] = mesh->face[i][2];
	}

	BOOL* flag = new BOOL[face_N];
	for(i=0; i<face_N; i++){
		//if(vertex[face[i][0]][2] > 0)
			flag[i] = true;
		//else
			//flag[i] = false;
	}

	this->subdivideTopology(flag);

	for(i=0; i<vertex_N; i++){
		new_vertex[i][0] = vertex[i][0];
		new_vertex[i][1] = vertex[i][1];
		new_vertex[i][2] = vertex[i][2];
	}

	int counter= 0;
	for(i=0; i<face_N; i++){
		if(!flag[i])
			continue;
		for(int j=0; j<3; j++){
			int pair = link[i][j];
			if(pair < i || !flag[pair]){
				int index = vertex_index[i][j];
				int i1 = face[i][j];
				int i2 = face[i][(j+1)%3];
				int i3 = face[i][(j+2)%3];
				
				double a = MeshData::DIST(vertex[i1], vertex[i2]);
				double b = MeshData::DIST(vertex[i3], vertex[i1]);
				double c = MeshData::DIST(vertex[i2], vertex[i3]);
				double s1 = (a+b-c)/(2.0*a);

				if(face[pair][0] == i1)
					i3 = face[pair][1];
				else if(face[pair][1] == i1)
					i3 = face[pair][2];
				else
					i3 = face[pair][0];

				b = MeshData::DIST(vertex[i3], vertex[i1]);
				c = MeshData::DIST(vertex[i2], vertex[i3]);
				double s2 = (a+b-c)/(2.0*a);

				s1 = s2 = 0.5;

				new_vertex[index][0] = (float)(0.5*((s1+s2)*vertex[i2][0] + (2.0-s1-s2)*vertex[i1][0]));
				new_vertex[index][1] = (float)(0.5*((s1+s2)*vertex[i2][1] + (2.0-s1-s2)*vertex[i1][1]));
				new_vertex[index][2] = (float)(0.5*((s1+s2)*vertex[i2][2] + (2.0-s1-s2)*vertex[i1][2]));
			}
		}
	}
	
	this->updateMesh();
	
	mesh->computeFaceNormal();
	float (*normal)[3] = mesh->normal_f;
	double (*new_normal)[3] = new double[new_face_N][3]; 
		
	//face = mesh->face;
	vertex = mesh->vertex;
	int **link_f = mesh->vertex_link_f;
	int **link_v = mesh->vertex_link_v;
	int *degree_f = mesh->degree_f;
	int *degree_v = mesh->degree_v;

	for(i=0; i<vertex_N; i++){

		int deg = degree_f[i];
		int* l = link_f[i];
		int* l_v = link_v[i];

		double *length = new double[deg];
		double total_length = 0;
		for(int j=0; j<deg; j++){
			int f1 = l[j];
			int f2 = l[(j+1)%deg];
			int k = l_v[j];
			double c1[3], c2[3], m[3];
			mesh->faceCenter(c1, f1);
			mesh->faceCenter(c2, f2);
			/*
			m[0] = 0.5*(vertex[i][0] + vertex[k][0]);
			m[1] = 0.5*(vertex[i][1] + vertex[k][1]);
			m[2] = 0.5*(vertex[i][2] + vertex[k][2]);
			length[j] = MeshData::DIST(c1,m) + MeshData::DIST(m,c2);
			*/
			
			double v1[3], v2[3];
			MeshData::VEC(v1, vertex[i], c1);
			MeshData::VEC(v2, vertex[i], c2);
			double dot = MeshData::DOT(v1, v2)/(MeshData::LENGTH(v1)*MeshData::LENGTH(v2));
			length[j] = acos(dot);
			//length[j] = mesh->faceArea(f1) + mesh->faceArea(f2);
			total_length += length[j];
		}
		for(j=0; j<deg; j++)
			length[j] /= total_length;

		for(j=0; j<deg; j++){
			double nj[3];
			nj[0] = normal[l[j]][0];
			nj[1] = normal[l[j]][1];
			nj[2] = normal[l[j]][2];

			double total_w = (deg+5.0)/(double)(4*deg);
			total_length = 0;

			for(int m=j+1; m<j+deg; m++){
				int k = m%deg;
				int d = abs(j-k);
				if(d == 0)
					continue;
				double w = (3.0 + 2.0*cos(2*PI*d/deg))/(double)(4*deg);
				//total_length += length[k];
				//double w = (3.0 + 2.0*cos(2*PI*total_length))/(double)(4*deg);

				double ro[3];
				MeshData::CROSS(ro, nj, normal[l[k]]);
				double len = MeshData::LENGTH(ro);
				if(len != 0){
					ro[0] /= len;
					ro[1] /= len;
					ro[2] /= len;

					double dot = MeshData::DOT(nj, normal[l[k]]);
					if(dot > 1)
						dot = 1;
					else if(dot < -1)
						dot = -1;

					double angle = w*acos(dot)/(total_w+w);

					double R[3][3];
					GENERATE_MAT(R, angle, ro);

					double tmp[3];
					tmp[0] = nj[0];
					tmp[1] = nj[1];
					tmp[2] = nj[2];
					MAT_VEC(nj, R, tmp);
				}
				total_w += w;
			}
			new_normal[l[j]][0] = nj[0];
			new_normal[l[j]][1] = nj[1];
			new_normal[l[j]][2] = nj[2];

			/*
			double U[3][3];
			double nj[3];
			nj[0] = normal[l[j]][0];
			nj[1] = normal[l[j]][1];
			nj[2] = normal[l[j]][2];
			double c[3];
			mesh->faceCenter(c, l[j]);
			double m[3];
			MeshData::VEC(m, vertex[i], c);
			ROTATION_MATRIX(U, nj, m);
			
			double test[3];
			MAT_VEC(test, U, m);

			double theta = 0;
			double phi = 0;	

			for(int k=0; k<deg; k++){
				int d = abs(j-k);
				if(d == 0)
					continue;
				double w = (3.0 + 2.0*cos(2*PI*d/deg))/(double)(4*deg);

				double nk[3];
				nk[0] = normal[l[k]][0];
				nk[1] = normal[l[k]][1];
				nk[2] = normal[l[k]][2];
				double n_tmp[3];
				MAT_VEC(n_tmp, U, nk);
				double theta_k, phi_k;
				ANGLE(theta_k, phi_k, n_tmp);

				theta += w*theta_k;
				phi += w*phi_k;
			}
			theta = 0;
			
			nj[0] = sin(theta)*sin(phi);
			nj[1] = cos(theta)*sin(phi);
			nj[2] = cos(phi);
			
			INVERSE_U(U);

			MAT_VEC(new_normal[l[j]], U, nj);
			*/
			/*
			double tmp_n[3] = {0,0,0};
			for(int k=0; k<deg; k++){
				int d = abs(j-k);
				double w;
				
				//DOO-SABIN
				if(d != 0){
					w = (3.0 + 2.0*cos(2*PI*d/deg))/(double)(4*deg);
				}
				else
					w = (deg+5.0)/(double)(4*deg);

				//CATMULL and CLARK
				
				//if(d == 0)
					//w = 0.5 + 1.0/(4.0*deg);
				//else if(d == 1 || d == deg-1)
					//w = 0.125 + 1.0/(4.0*deg);
				//else
					//w = 1.0/(4.0*deg);

				tmp_n[0] += w*normal[l[k]][0];
				tmp_n[1] += w*normal[l[k]][1];
				tmp_n[2] += w*normal[l[k]][2];
			}
			double len = MeshData::LENGTH(tmp_n);
			if(len != 0){
				new_normal[l[j]][0] = tmp_n[0]/len;
				new_normal[l[j]][1] = tmp_n[1]/len;
				new_normal[l[j]][2] = tmp_n[2]/len;
			}*/
		}
		delete[] length;
	}
	
	for(i=0; i<face_N; i++){
		if(!flag[i])
			continue;
		for(int j=0; j<3; j++){
			int pair = link[i][j];
			if(pair < i || !flag[pair]){
				int index = vertex_index[i][j];
				int i1 = face[i][j];
				int i2 = face[i][(j+1)%3];

				//double pro[3] = {0,0,0};

				int *l_v = link_v[i1];
				int *l_f = link_f[i1];
				int deg = degree_v[i1];
				int f[4];
				for(int k=0; k<deg; k++){
					if(index == l_v[k]){
						f[0] = l_f[k];
						f[1] = l_f[(k+deg-1)%deg];
						break;
					}
				}
				l_v = link_v[i2];
				l_f = link_f[i2];
				deg = degree_v[i2];

				for(k=0; k<deg; k++){
					if(index == l_v[k]){
						f[2] = l_f[k];
						f[3] = l_f[(k+deg-1)%deg];
						break;
					}
				}
	
				double A[6] = {0,0,0,0,0,0};
				double d[3] = {0,0,0};
				double n[3] = {0,0,0};
				double p[3] = {0,0,0};
				for(k=0; k<4; k++){
					int fi = f[k];
					double area = 1; //mesh->faceArea(fi);
					n[0] += area*normal[fi][0];
					n[1] += area*normal[fi][1];
					n[2] += area*normal[fi][2];

					A[0] += area*new_normal[fi][0] * new_normal[fi][0];
					A[1] += area*new_normal[fi][0] * new_normal[fi][1];
					A[2] += area*new_normal[fi][0] * new_normal[fi][2];

					A[3] += area*new_normal[fi][1] * new_normal[fi][1];
					A[4] += area*new_normal[fi][1] * new_normal[fi][2];

					A[5] += area*new_normal[fi][2] * new_normal[fi][2];
					
					double dot;
					if(k < 2)
						dot = MeshData::DOT(new_normal[fi], vertex[i1]);
					else
						dot = MeshData::DOT(new_normal[fi], vertex[i2]);

					d[0] += area*dot*new_normal[fi][0];
					d[1] += area*dot*new_normal[fi][1];
					d[2] += area*dot*new_normal[fi][2];

					/*
					float c[3];
					mesh->faceCenter(c, fi);
					p[0] += 0.25*c[0];
					p[1] += 0.25*c[1];
					p[2] += 0.25*c[2];*/
				}

				//n[0] = L[i1][0] + L[i2][0];
				//n[1] = L[i1][1] + L[i2][1];
				//n[2] = L[i1][2] + L[i2][2];

				//t = ((d-Ap).n)/((An).n)
				
				p[0] = vertex[index][0];
				p[1] = vertex[index][1];
				p[2] = vertex[index][2];
				/*
				for(int m=0; m<degree_v[index]; m++){
					p[0] += vertex[link_v[index][m]][0];
					p[1] += vertex[link_v[index][m]][1];
					p[2] += vertex[link_v[index][m]][2];
				}
				p[0] = 0.5*(vertex[index][0] + p[0]/m);
				p[1] = 0.5*(vertex[index][1] + p[1]/m);
				p[2] = 0.5*(vertex[index][2] + p[2]/m);
				*/
				double over = (d[0]-A[0]*p[0] - A[1]*p[1] - A[2]*p[2])*n[0] +
							  (d[1]-A[1]*p[0] - A[3]*p[1] - A[4]*p[2])*n[1] +
							  (d[2]-A[2]*p[0] - A[4]*p[1] - A[5]*p[2])*n[2] ;
				double under = (A[0]*n[0] + A[1]*n[1] + A[2]*n[2])*n[0] + 
							   (A[1]*n[0] + A[3]*n[1] + A[4]*n[2])*n[1] + 
							   (A[2]*n[0] + A[4]*n[1] + A[5]*n[2])*n[2] ;
				if(under != 0){
					double t = over/under;
					vertex[index][0] = (float)(p[0] + t*n[0]);
					vertex[index][1] = (float)(p[1] + t*n[1]);
					vertex[index][2] = (float)(p[2] + t*n[2]);
				}
				/*
				double pro[3] = {0,0,0};

				int *l_v = link_v[i1];
				int *l_f = link_f[i1];
				int deg = degree_v[i1];
				int f1, f2;
				for(int k=0; k<deg; k++){
					if(index == l_v[k]){
						f1 = l_f[k];
						f2 = l_f[(k+deg-1)%deg];
						break;
					}
				}

				float v[3];
				MeshData::VEC(v, vertex[index], vertex[i1]);

				double dot = MeshData::DOT(new_normal[f1], v);
				pro[0] += dot*new_normal[f1][0];
				pro[1] += dot*new_normal[f1][1];
				pro[2] += dot*new_normal[f1][2];

				dot = MeshData::DOT(new_normal[f2], v);
				pro[0] += dot*new_normal[f2][0];
				pro[1] += dot*new_normal[f2][1];
				pro[2] += dot*new_normal[f2][2];

				l_v = link_v[i2];
				l_f = link_f[i2];
				deg = degree_v[i2];

				for(k=0; k<deg; k++){
					if(index == l_v[k]){
						f1 = l_f[k];
						f2 = l_f[(k+deg-1)%deg];
						break;
					}
				}

				MeshData::VEC(v, vertex[index], vertex[i2]);

				dot = MeshData::DOT(new_normal[f1], v);
				pro[0] += dot*new_normal[f1][0];
				pro[1] += dot*new_normal[f1][1];
				pro[2] += dot*new_normal[f1][2];

				dot = MeshData::DOT(new_normal[f2], v);
				pro[0] += dot*new_normal[f2][0];
				pro[1] += dot*new_normal[f2][1];
				pro[2] += dot*new_normal[f2][2];

				vertex[index][0] += (float)(0.25*pro[0]);
				vertex[index][1] += (float)(0.25*pro[1]);
				vertex[index][2] += (float)(0.25*pro[2]);
				*/
			}
		}
	}

	delete[] face;
	delete[] link;
	delete[] new_normal;
	delete[] flag;
	delete[] vertex_index;
	vertex_index = NULL;
}

void Subdivider::subdivideTopologySqrt3()
{
	decideCenterForSqrt3();

	int face_N = mesh->face_N;
	int vertex_N = mesh->vertex_N;
	int **link = mesh->vertex_link_f;
	float (*vertex)[3] = mesh->vertex;
	float (*center)[3] = mesh->center;
	int (*face)[3] = mesh->face;
	int* degree = mesh->degree_f;

	new_vertex_N = vertex_N + face_N;
	new_vertex = new float[new_vertex_N][3];

	new_face_N = face_N;
	for(int i=0; i<vertex_N; i++)
		new_face_N += degree[i];
	new_face = new int[new_face_N][3];

	int face_id = 0;
	for(i=0; i<vertex_N; i++){
		int deg = degree[i];
		int* l = link[i];
		for(int j=0; j<deg; j++){
			new_face[face_id][0] = i;
			new_face[face_id][1] = vertex_N + l[j];
			new_face[face_id][2] = vertex_N + l[(j+1)%deg];
			face_id++;
		}
	}
	
	for(i=0; i<vertex_N; i++){
		new_vertex[i][0] = vertex[i][0];
		new_vertex[i][1] = vertex[i][1];
		new_vertex[i][2] = vertex[i][2];
	}

	for(i=0; i<face_N; i++){
		int j = vertex_N + i;
		new_vertex[j][0] = center[i][0];
		new_vertex[j][1] = center[i][1];
		new_vertex[j][2] = center[i][2];
	}

	this->updateMesh();
}

void Subdivider::decideCenterForSqrt3()
{
	mesh->computeCenter();
	mesh->computeFaceNormal();

	int face_N = mesh->face_N;
	int (*link)[3] = mesh->face_link_E;
	float (*vertex)[3] = mesh->vertex;
	float (*center)[3] = mesh->center;
	int (*face)[3] = mesh->face;
	float (*normal)[3] = mesh->normal_f;

	double (*tmp_normal)[3][3] = new double[face_N][3][3];

	/*
	for(int i=0; i<face_N; i++){
		float *n = normal[i];
		double ni[3][3];
		for(int j=0; j<3; j++){
			int index = link[i][j];
			float *m = normal[index];

			double ro[3];
			MeshData::CROSS(ro, n, m);
			double len = MeshData::LENGTH(ro);
			if(len != 0){
				ro[0] /= len;
				ro[1] /= len;
				ro[2] /= len;

				double dot = MeshData::DOT(n, m);
				if(dot > 1)
					dot = 1;
				else if(dot < -1)
					dot = -1;

				double angle = 1.0*acos(dot)/2.0;

				double R[3][3];
				GENERATE_MAT(R, angle, ro);

				double tmp[3];
				tmp[0] = n[0];
				tmp[1] = n[1];
				tmp[2] = n[2];
				MAT_VEC(ni[j], R, tmp);
			}
			else{
				ni[j][0] = n[0];
				ni[j][1] = n[1];
				ni[j][2] = n[2];
			}
		}

		for(j=0; j<3; j++){
			double ro[3];
			double* n1 = ni[j];
			double* n2 = ni[(j+2)%3];
			MeshData::CROSS(ro, n1, n2);
			double len = MeshData::LENGTH(ro);
			if(len != 0){
				ro[0] /= len;
				ro[1] /= len;
				ro[2] /= len;

				double dot = MeshData::DOT(n1, n2);
				if(dot > 1)
					dot = 1;
				else if(dot < -1)
					dot = -1;

				double angle = 0.0*acos(dot);

				double R[3][3];
				GENERATE_MAT(R, angle, ro);

				MAT_VEC(tmp_normal[i][2*j], R, n1);
			
				angle = -angle;

				GENERATE_MAT(R, angle, ro);

				MAT_VEC(tmp_normal[i][2*j+1], R, n2);
			}
			else{
				tmp_normal[i][2*j][0] = n1[0];
				tmp_normal[i][2*j][1] = n1[1];
				tmp_normal[i][2*j][2] = n1[2];

				tmp_normal[i][2*j+1][0] = n2[0];
				tmp_normal[i][2*j+1][1] = n2[1];
				tmp_normal[i][2*j+1][2] = n2[2];
			}
		}
	}*/
	
	int **link_f = mesh->vertex_link_f;
	int **link_v = mesh->vertex_link_v;
	int vertex_N = mesh->vertex_N;
	int *degree = mesh->degree_f;
	for(int i=0; i<vertex_N; i++){
		int deg = degree[i];
		int* l = link_f[i];
		int* l_v = link_v[i];
	
		double *length = new double[deg];
		double total_length = 0;
		for(int j=0; j<deg; j++){
			int f1 = l[j];
			int f2 = l[(j+1)%deg];
			int k = l_v[j];
			double c1[3], c2[3], m[3];
			mesh->faceCenter(c1, f1);
			mesh->faceCenter(c2, f2);
			/*
			m[0] = 0.5*(vertex[i][0] + vertex[k][0]);
			m[1] = 0.5*(vertex[i][1] + vertex[k][1]);
			m[2] = 0.5*(vertex[i][2] + vertex[k][2]);
			length[j] = MeshData::DIST(c1,m) + MeshData::DIST(m,c2);
			*/
			
			double v1[3], v2[3];
			MeshData::VEC(v1, vertex[i], c1);
			MeshData::VEC(v2, vertex[i], c2);
			double dot = MeshData::DOT(v1, v2)/(MeshData::LENGTH(v1)*MeshData::LENGTH(v2));
			length[j] = acos(dot);
			//length[j] = mesh->faceArea(f1) + mesh->faceArea(f2);
			total_length += length[j];
		}
		for(j=0; j<deg; j++)
			length[j] /= total_length;

		for(j=0; j<deg; j++){
			double nj[3];
			nj[0] = normal[l[j]][0];
			nj[1] = normal[l[j]][1];
			nj[2] = normal[l[j]][2];

			double total_w = (deg+5.0)/(double)(4*deg);
			total_length = 0;

			for(int m=j+1; m<j+deg; m++){
				int k = m%deg;
				int d = abs(j-k);
				if(d == 0)
					continue;
				double w = (3.0 + 2.0*cos(2*PI*d/deg))/(double)(4*deg);
				//total_length += length[k];
				//double w = (3.0 + 2.0*cos(2*PI*total_length))/(double)(4*deg);

				double ro[3];
				MeshData::CROSS(ro, nj, normal[l[k]]);
				double len = MeshData::LENGTH(ro);
				if(len != 0){
					ro[0] /= len;
					ro[1] /= len;
					ro[2] /= len;

					double dot = MeshData::DOT(nj, normal[l[k]]);
					if(dot > 1)
						dot = 1;
					else if(dot < -1)
						dot = -1;

					double angle = w*acos(dot)/(total_w+w);

					double R[3][3];
					GENERATE_MAT(R, angle, ro);

					double tmp[3];
					tmp[0] = nj[0];
					tmp[1] = nj[1];
					tmp[2] = nj[2];
					MAT_VEC(nj, R, tmp);
				}
				total_w += w;
			}
			int f = l[j];
			int index;
			if(face[f][0] == i)
				index = 0;
			else if(face[f][1] == i)
				index = 1;
			else if(face[f][2] == i)
				index = 2;
			else
				int a = 0;

			tmp_normal[f][index][0] = nj[0];
			tmp_normal[f][index][1] = nj[1];
			tmp_normal[f][index][2] = nj[2];
		}
	}

	float (*tmp_center)[3] = new float[face_N][3];
	//for(int k=0; k<10; k++){
	for(i=0;i<face_N; i++){
		double A[6] = {0,0,0,0,0,0};
		double d[3] = {0,0,0};
		float *n = normal[i];
		//float n[] = {0,0,0};
		for(int j=0; j<3; j++){
			double *m = tmp_normal[i][j];

			//n[0] += m[0];
			//n[1] += m[1];
			//n[2] += m[2];

			A[0] += m[0] * m[0];
			A[1] += m[0] * m[1];
			A[2] += m[0] * m[2];

			A[3] += m[1] * m[1];
			A[4] += m[1] * m[2];

			A[5] += m[2] * m[2];
			int i1 = face[i][j];
			double dot = MeshData::DOT(m, vertex[i1]);
				
			d[0] += dot*m[0];
			d[1] += dot*m[1];
			d[2] += dot*m[2];
		}			
		float* p = center[i];
		/*
		float p[] = {0,0,0};
		for(j=0; j<3; j++){
			int index = link[i][j];
			p[0] += center[index][0];
			p[1] += center[index][1];
			p[2] += center[index][2];
		}
		p[0] /= 3;
		p[1] /= 3;
		p[2] /= 3;

		n[0] = center[i][0] - p[0];
		n[1] = center[i][1] - p[1];
		n[2] = center[i][2] - p[2];

		p[0] = center[i][0];
		p[1] = center[i][1];
		p[2] = center[i][2];*/

		double over = (d[0]-A[0]*p[0] - A[1]*p[1] - A[2]*p[2])*n[0] +
					  (d[1]-A[1]*p[0] - A[3]*p[1] - A[4]*p[2])*n[1] +
					  (d[2]-A[2]*p[0] - A[4]*p[1] - A[5]*p[2])*n[2] ;
		double under = (A[0]*n[0] + A[1]*n[1] + A[2]*n[2])*n[0] + 
					   (A[1]*n[0] + A[3]*n[1] + A[4]*n[2])*n[1] + 
					   (A[2]*n[0] + A[4]*n[1] + A[5]*n[2])*n[2] ;
		if(under != 0){
			double t = over/under;
			tmp_center[i][0] = (float)(p[0] + t*n[0]);
			tmp_center[i][1] = (float)(p[1] + t*n[1]); 
			tmp_center[i][2] = (float)(p[2] + t*n[2]);
		}
		else{
			tmp_center[i][0] = p[0];
			tmp_center[i][1] = p[1];
			tmp_center[i][2] = p[2];
		}
		/*
		double B[6];
		if(INVERSE(B, A)){
			double v[3];
			MAT_BY_VEC(v, B, d);
			tmp_center[i][0] = v[0];
			tmp_center[i][1] = v[1];
			tmp_center[i][2] = v[2];
		}
		else{
			tmp_center[i][0] = p[0];
			tmp_center[i][1] = p[1];
			tmp_center[i][2] = p[2];
		}*/
	}

	for(i=0;i<face_N; i++){
		center[i][0] = tmp_center[i][0];
		center[i][1] = tmp_center[i][1];
		center[i][2] = tmp_center[i][2];
	}
	//}

	/*
	mesh->computeNormal();
	for(i=0;i<vertex_N; i++){
		int deg = degree[i];
		int *l = link_f[i];

		double A[6] = {0,0,0,0,0,0};
		double d[3] = {0,0,0};
		
		//float *n = mesh->normal[i];
		for(int j=0; j<deg; j++){
			int f = l[j];
			int index;
			if(face[f][0] == i)
				index = 0;
			else if(face[f][1] == i)
				index = 1;
			else
				index = 2;

			double *m = tmp_normal[f][index];
			A[0] += m[0] * m[0];
			A[1] += m[0] * m[1];
			A[2] += m[0] * m[2];

			A[3] += m[1] * m[1];
			A[4] += m[1] * m[2];

			A[5] += m[2] * m[2];
			
			double dot = MeshData::DOT(m, center[f]);
				
			d[0] += dot*m[0];
			d[1] += dot*m[1];
			d[2] += dot*m[2];
		}
		double n[3];
		mesh->laplacian(i, n);
		float *p = vertex[i];
		double over = (d[0]-A[0]*p[0] - A[1]*p[1] - A[2]*p[2])*n[0] +
					  (d[1]-A[1]*p[0] - A[3]*p[1] - A[4]*p[2])*n[1] +
					  (d[2]-A[2]*p[0] - A[4]*p[1] - A[5]*p[2])*n[2] ;
		double under = (A[0]*n[0] + A[1]*n[1] + A[2]*n[2])*n[0] + 
					   (A[1]*n[0] + A[3]*n[1] + A[4]*n[2])*n[1] + 
					   (A[2]*n[0] + A[4]*n[1] + A[5]*n[2])*n[2] ;

		if(under == 0)
			continue;

		double t = over/under;
		
		vertex[i][0] += (float)(t*n[0]);
		vertex[i][1] += (float)(t*n[1]);
		vertex[i][2] += (float)(t*n[2]);
		
		//double B[6];
		//if(INVERSE(B, A)){
			//double v[3];
			//MAT_BY_VEC(v, B, d);
			//vertex[i][0] = v[0];
			//vertex[i][1] = v[1];
			//vertex[i][2] = v[2];
		//}
	}*/

	delete[] tmp_normal;
	delete[] tmp_center;
	
}

void Subdivider::subdivideTopologyBary()
{
	decideCenterForSqrt3();

	int face_N = mesh->face_N;
	int vertex_N = mesh->vertex_N;
	int **link = mesh->vertex_link_f;
	float (*vertex)[3] = mesh->vertex;
	float (*center)[3] = mesh->center;
	int (*face)[3] = mesh->face;
	int* degree = mesh->degree_f;

	new_vertex_N = vertex_N + face_N;
	new_vertex = new float[new_vertex_N][3];

	new_face_N = 3*face_N;
	new_face = new int[new_face_N][3];

	int face_id = 0;
	for(int i=0; i<face_N; i++){
		for(int j=0; j<3; j++){
			new_face[face_id][0] = face[i][j];
			new_face[face_id][1] = face[i][(j+1)%3];
			new_face[face_id][2] = vertex_N+i;
			face_id++;
		}
	}
	
	for(i=0; i<vertex_N; i++){
		new_vertex[i][0] = vertex[i][0];
		new_vertex[i][1] = vertex[i][1];
		new_vertex[i][2] = vertex[i][2];
	}

	for(i=0; i<face_N; i++){
		int j = vertex_N + i;
		new_vertex[j][0] = center[i][0];
		new_vertex[j][1] = center[i][1];
		new_vertex[j][2] = center[i][2];
	}

	this->updateMesh();
}