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

#include "stdafx.h"
#include "MeshEditor.h"
#include "FutureDetector.h"

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

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

FutureDetector::FutureDetector()
{

}

FutureDetector::~FutureDetector()
{

}

void FutureDetector::ridge(MeshData *mesh, BOOL threshold_on)
{
	/*
	int i;
	int vertex_N = mesh->vertex_N;
	//allocate memory
	if(mesh->isRidge != NULL)
		delete[] mesh->isRidge;
	mesh->isRidge = new BOOL[vertex_N];
	BOOL *isRidge = mesh->isRidge;

	if(mesh->isRavine != NULL)
		delete[] mesh->isRavine;
	mesh->isRavine = new BOOL[vertex_N];
	BOOL *isRavine = mesh->isRavine;
	
	for(i=0; i<vertex_N; i++){
		if(mesh->isBound[i]){
			isRidge[i] = isRavine[i] = false;
			continue;
		}

		if(k_max[i] > 0)
			isRidge[i] = true;
		else
			isRidge[i] = false;

		
		if(k_min[i] < 0)
			isRavine[i] = true;
		else
			isRavine[i] = false;
	
		int degree, *nei;
		mesh->get1Ring(i, neu, degree);
	
		float k_max_max1, k_max_max2, k_max_min1, k_max_min2;
		float k_min_max1, k_min_max2, k_min_min1, k_min_min2;
		float l_max1, l_max2, l_min1, l_min2;

		float n1[3];
		CROSS(n1, mesh->normal[i], mesh->t_max[i]);
		float n2[3];
		CROSS(n2, mesh->normal[i], mesh->t_min[i]);

		for(int j=0; j<degree; j++){
		}

		for(int face_id = start_face[i];(face_id >= 0)&&(m < degree);face_id = nextFace(face_id, i)){
			int j = nextVertex(face_id, i);
			int k = nextVertex(face_id, j);

			float t1[3];
			VEC(t1, vertex[j], vertex[i]);
			float t2[3];
			VEC(t2, vertex[k], vertex[i]);

			
			float in1 = DOT(n1, t1);
			float in2 = DOT(n1, t2);
			if(in1 * in2 <= 0){
				float v[3];
				BOOL isOne = true;
				if(in1 < 0)
					in1 = -in1;
				if(in2 < 0){
					in2 = -in2;
					isOne = false;
				}
				v[0] = (in2*t1[0] + in1*t2[0])/(in1 + in2);
				v[1] = (in2*t1[1] + in1*t2[1])/(in1 + in2);
				v[2] = (in2*t1[2] + in1*t2[2])/(in1 + in2);

				float o[3] = {0,0,0};
				float area1 = AREA(o, t1, v);
				float area2 = AREA(o, t2, v);
				float k_max_in = (area2*k_max[j] + area1*k_max[k])/(area1+area2);
				float k_min_in = (area2*k_min[j] + area1*k_min[k])/(area1+area2);

				if(isOne){
					k_max_max1 = k_max_in;
					k_min_max1 = k_min_in;
					l_max1 = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
				}
				else{
					k_max_max2 = k_max_in;
					k_min_max2 = k_min_in;
					l_max2 = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
				}
			}

			in1 = DOT(n2, t1);
			in2 = DOT(n2, t2);
			if(in1 * in2 <= 0){
				float v[3];
				BOOL isOne = true;
				if(in1 < 0)
					in1 = -in1;
				if(in2 < 0){
					in2 = -in2;
					isOne = false;
				}
				v[0] = (in2*t1[0] + in1*t2[0])/(in1 + in2);
				v[1] = (in2*t1[1] + in1*t2[1])/(in1 + in2);
				v[2] = (in2*t1[2] + in1*t2[2])/(in1 + in2);

				float o[3] = {0,0,0};
				float area1 = AREA(o, t1, v);
				float area2 = AREA(o, t2, v);
				float k_max_in = (area2*k_max[j] + area1*k_max[k])/(area1+area2);
				float k_min_in = (area2*k_min[j] + area1*k_min[k])/(area1+area2);

				if(isOne){
					k_max_min1 = k_max_in;
					k_min_min1 = k_min_in;
					l_min1 = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
				}
				else{
					k_max_min2 = k_max_in;
					k_min_min2 = k_min_in;
					l_min2 = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
				}	
			}
			m++;
		}
		double d1 = sqrt(l_min1);
		double d2 = sqrt(l_min2);
		
		float b = (k_max_min2-k_max[i])/d2 + (k_max[i]-k_max_min1)/d1 
			- (k_max_min2-k_max_min1)/(d1+d2);

		float K1 = k_max_max1 + 0.5*l_max1*b*b/(k_max[i]-k_min[i]);
		float K2 = k_max_max2 + 0.5*l_max2*b*b/(k_max[i]-k_min[i]);

		isRidge[i] &= (k_max[i] > K1) && (k_max[i] > K2);
	
		d1 = sqrt(l_max1);
		d2 = sqrt(l_max2);
	
		b = (k_min_max2-k_min[i])/d2 + (k_min[i]-k_min_max1)/d1 
			- (k_min_max2-k_min_max1)/(d1+d2);

		K1 = k_min_min1 - 0.5*l_min1*b*b/(k_max[i]-k_min[i]);
		K2 = k_min_min2 - 0.5*l_min2*b*b/(k_max[i]-k_min[i]);

		isRavine[i] &= (k_min[i] < K1) && (k_min[i] < K2);
	}

	BOOL* weak_ridge = new BOOL[vertex_N];
	BOOL* weak_ravine = new BOOL[vertex_N];
	for(i=0; i<vertex_N; i++){
		weak_ridge[i] = false;
		weak_ravine[i] = false;
	}

	float k_max_ave = 0;
	float k_min_ave = 0;
	for(i=0; i<vertex_N; i++){
		if(adjacent[i]->next != NULL){
			k_max_ave += k_max[i];
			k_min_ave += k_min[i];
		}
	}
	k_max_ave /= (float)vertex_N;
	k_min_ave /= (float)vertex_N;

	float k_max_var = 0;
	float k_min_var = 0;
	for(i=0; i<vertex_N; i++){
		if(adjacent[i]->next != NULL){
			k_max_var += k_max[i]*k_max[i];
			k_min_var += k_min[i]*k_min[i];
		}
	}
	k_max_var = k_max_var/((float)vertex_N) - k_max_ave*k_max_ave;
	k_min_var = k_min_var/((float)vertex_N) - k_min_ave*k_min_ave;
	
	float ridge_hi = k_max_ave + 0.25*sqrt(k_max_var);
	float ridge_ro = k_max_ave - 0.53*sqrt(k_max_var);

	float ravine_hi = k_min_ave - 0.25*sqrt(k_min_var);
	float ravine_ro = k_min_ave + 0.53*sqrt(k_min_var);

	for(i=0; i<vertex_N; i++){
		if(isBound[i])
			continue;

		if(isRidge[i]){
			if(k_max[i] > ridge_hi)
				isRidge[i] = true;
			else
				isRidge[i] = false;
			if(k_max[i] > ridge_ro)
				weak_ridge[i] = true;
			else
				weak_ridge[i] = false;
		}
		else
			weak_ridge[i] = false;

		if(isRavine[i]){
			if(k_min[i] < ravine_hi)
				isRavine[i] = true;
			else
				isRavine[i] = false;
			if(k_min[i] < ravine_ro)
				weak_ravine[i] = true;
			else
				weak_ravine[i] = false;
		}
		else
			weak_ravine[i] = false;
	}

	for(i=0; i<vertex_N; i++){
		if(isBound[i])
			continue;

		if(isRidge[i]){
			BOOL flag_ridge = false;
			for(node* current=adjacent[i]; current->next!=NULL; current=current->next->next)
				flag_ridge |= isRidge[current->e];
			if(!flag_ridge)
				isRidge[i] = false;
		}
	
		if(isRavine[i]){
			BOOL flag_ravine = false;
			for(node* current=adjacent[i]; current->next!=NULL; current=current->next->next)
				flag_ravine |= isRavine[current->e];
			if(!flag_ravine)
				isRavine[i] = false;
		}
	}

	BOOL isGrow = true;
	while(isGrow){
		isGrow = false;
		for(i=0; i<vertex_N; i++){
			if(isBound[i])
				continue;
*/
			/*************************
			if(isRidge[i])
				for(node* current=adjacent[i]; current->next!=NULL; current=current->next->next)
					if(weak_ridge[current->e] && !isRidge[current->e]){
						isRidge[current->e] = true;
						isGrow = true;
					}

			if(isRavine[i])
				for(node* current=adjacent[i]; current->next!=NULL; current=current->next->next)
					if(weak_ravine[current->e] && !isRavine[current->e]){
						isRavine[current->e] = true;
						isGrow = true;
					}
			***********************/

			/*********************
			if(!isRidge[i] && !isRavine[i])
				continue;

			int degree = 0;
			for(node *current=adjacent[i]; current->next!=NULL; current=current->next->next)
				degree++;
		
			int m = 0;
			for(int face_id = start_face[i];(face_id >= 0)&&(m < degree);face_id = nextFace(face_id, i)){
				int j = nextVertex(face_id, i);
				int k = nextVertex(face_id, j);

				float v[3];
				VEC(v, vertex[j], vertex[i]);
				float t1[3];
				t1[0] = (1.0f-normal[i][0]*normal[i][0])*v[0] 
					- normal[i][0]*normal[i][1]*v[1] - normal[i][0]*normal[i][2]*v[2];
				t1[1] = -normal[i][0]*normal[i][1]*v[0] + (1.0f-normal[i][1]*normal[i][1])*v[1] 
					- normal[i][1]*normal[i][2]*v[2];
				t1[2] = -normal[i][0]*normal[i][2]*v[0] - normal[i][1]*normal[i][2]*v[1] 
					+ (1.0f-normal[i][2]*normal[i][2])*v[2];

				VEC(v, vertex[k], vertex[i]);
				float t2[3];
				t2[0] = (1.0f-normal[i][0]*normal[i][0])*v[0] 
					- normal[i][0]*normal[i][1]*v[1] - normal[i][0]*normal[i][2]*v[2];
				t2[1] = -normal[i][0]*normal[i][1]*v[0] + (1.0f-normal[i][1]*normal[i][1])*v[1] 
					- normal[i][1]*normal[i][2]*v[2];
				t2[2] = -normal[i][0]*normal[i][2]*v[0] - normal[i][1]*normal[i][2]*v[1] 
					+ (1.0f-normal[i][2]*normal[i][2])*v[2];

				float n1[3];
				CROSS(n1, normal[i], t1);
				float n2[3];
				CROSS(n2, t2, normal[i]);

				if(isRidge[i] && DOT(n1,t_min[i])*DOT(n2,t_min[i])>=0){
					float o[3] = {0,0,0};
					float area1 = AREA(o, t_min[i], t1);
					float area2 = AREA(o, t_min[i], t2);
					if(!isRidge[j] && weak_ridge[j] && area1 < area2){
						isGrow = true;
						isRidge[j] = true;
					}
					else if(!isRidge[k] && weak_ridge[k]){
						isGrow = true;
						isRidge[k] = true;
					}
				}

				if(isRavine[i] && DOT(n1,t_max[i])*DOT(n2,t_max[i])>=0){
					float o[3] = {0,0,0};
					float area1 = AREA(o, t_max[i], t1);
					float area2 = AREA(o, t_max[i], t2);
					if(!isRavine[j] && weak_ravine[j] && area1 < area2){
						isGrow = true;
						isRavine[j] = true;
					}
					else if(!isRavine[k] && weak_ravine[k]){
						isGrow = true;
						isRavine[k] = true;
					}
				}
				m++;
			}
			*********************/
		//}
	//}
	
	/******************************
	for(int N = 0; N<0; N++){
	for(i=0; i<vertex_N; i++){
		weak_ridge[i] = false;
		weak_ravine[i] = false;
	}

	for(i=0; i<vertex_N; i++){
		if(isBound[i]){
			continue;
		}

		if(!isRidge[i] && !isRavine[i])
			continue;

		int degree = 0;
		for(node *current=adjacent[i]; current->next!=NULL; current=current->next->next)
			degree++;

		int m = 0;
		for(int face_id = start_face[i];(face_id >= 0)&&(m < degree);face_id = nextFace(face_id, i)){
			int j = nextVertex(face_id, i);
			int k = nextVertex(face_id, j);

			float v[3];
			VEC(v, vertex[j], vertex[i]);
			float t1[3];
			t1[0] = (1.0f-normal[i][0]*normal[i][0])*v[0] 
				- normal[i][0]*normal[i][1]*v[1] - normal[i][0]*normal[i][2]*v[2];
			t1[1] = -normal[i][0]*normal[i][1]*v[0] + (1.0f-normal[i][1]*normal[i][1])*v[1] 
				- normal[i][1]*normal[i][2]*v[2];
			t1[2] = -normal[i][0]*normal[i][2]*v[0] - normal[i][1]*normal[i][2]*v[1] 
				+ (1.0f-normal[i][2]*normal[i][2])*v[2];

			VEC(v, vertex[k], vertex[i]);
			float t2[3];
			t2[0] = (1.0f-normal[i][0]*normal[i][0])*v[0] 
				- normal[i][0]*normal[i][1]*v[1] - normal[i][0]*normal[i][2]*v[2];
			t2[1] = -normal[i][0]*normal[i][1]*v[0] + (1.0f-normal[i][1]*normal[i][1])*v[1] 
				- normal[i][1]*normal[i][2]*v[2];
			t2[2] = -normal[i][0]*normal[i][2]*v[0] - normal[i][1]*normal[i][2]*v[1] 
				+ (1.0f-normal[i][2]*normal[i][2])*v[2];

			float n1[3];
			CROSS(n1, normal[i], t1);
			float n2[3];
			CROSS(n2, t2, normal[i]);

			if(isRidge[i] && DOT(n1,t_min[i])*DOT(n2,t_min[i])>=0){
				float o[3] = {0,0,0};
				float area1 = AREA(o, t_min[i], t1);
				float area2 = AREA(o, t_min[i], t2);
				if(area1 < area2)
					weak_ridge[j] = true;
				else
					weak_ridge[k] = true;
			}

			if(isRavine[i] && DOT(n1,t_max[i])*DOT(n2,t_max[i])>=0){
				float o[3] = {0,0,0};
				float area1 = AREA(o, t_max[i], t1);
				float area2 = AREA(o, t_max[i], t2);
				if(area1 < area2)
					weak_ravine[j] = true;
				else
					weak_ravine[k] = true;
				
			}
			m++;
		}
	}
	for(i=0; i<vertex_N; i++){
		isRidge[i] |= weak_ridge[i];
		isRavine[i] |= weak_ravine[i];
	}
	}
	***************************/
	
	//delete[] weak_ridge;
	//delete[] weak_ravine;
}