// Subdivider.h: Subdivider クラスのインターフェイス
//
//////////////////////////////////////////////////////////////////////

#if !defined(AFX_SUBDIVIDER_H__1168D90D_29F3_455D_8C32_4302007876E6__INCLUDED_)
#define AFX_SUBDIVIDER_H__1168D90D_29F3_455D_8C32_4302007876E6__INCLUDED_

#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000

#include "MeshData.h"

#define PI 3.14159265358979
#define EPSILON2 0.000000001
class Subdivider  
{
public:
	MeshData* mesh;
	int new_vertex_N;
	int new_face_N;
	int (*new_face)[3];
	float (*new_vertex)[3];
	int (*vertex_index)[3];

public:
	void subdivideTopologyBary();
	void decideCenterForSqrt3();
	void subdivideTopologySqrt3();
	void normalAverage();
	void butterfly();
	void LOOP();
	void linear();
	void setMeshData(MeshData* mesh);
	Subdivider();
	virtual ~Subdivider();

protected:
	void updateMesh();
	void subdivideTopology(BOOL* flag);

	//Vector v must be a unit vector, and v.x = 0
	static inline void ROTATION_MATRIX(double R[3][3], double v[3], double x[3]){
		double l = sqrt(v[0]*v[0] + v[1]*v[1]);
		double a;
		if(l == 0)
			a = 0;
		else
			a = acos(v[1]/l);
		if(v[0] < 0)
			a = -a;
		double b = acos(v[2]/MeshData::LENGTH(v));

		double Rz[3][3], Rx[3][3], tmp[3][3];
		Rz[0][0] = Rz[1][1] = cos(a);
		Rz[0][1] = -sin(a);
		Rz[1][0] = -Rz[0][1];
		Rz[2][0] = Rz[2][1] = Rz[0][2] = Rz[1][2] = 0;
		Rz[2][2] = 1;

		Rx[0][0] = 1;
		Rx[0][1] = Rx[0][2] = Rx[1][0] = Rx[2][0] = 0;
		Rx[1][1] = Rx[2][2] = cos(b);
		Rx[1][2] = -sin(b);
		Rx[2][1] = -Rx[1][2];

		MAT_TIMES(tmp, Rx, Rz);

		double x0[3];
		MAT_VEC(x0, tmp, x);
		double len = sqrt(x0[0]*x0[0] + x0[1]*x0[1]);
		if(len != 0){
			a = acos(x0[1]/len);
			if(x0[0] < 0)
				a = -a;
		}
		else
			a = 0;

		Rz[0][0] = Rz[1][1] = cos(a);
		Rz[0][1] = -sin(a);
		Rz[1][0] = -Rz[0][1];
		Rz[2][0] = Rz[2][1] = Rz[0][2] = Rz[1][2] = 0;
		Rz[2][2] = 1;

		MAT_TIMES(R, Rz, tmp);
	}

	//C = A.B
	static inline void MAT_TIMES(double C[3][3], double A[3][3], double B[3][3]){
		C[0][0] = A[0][0]*B[0][0] + A[0][1]*B[1][0] + A[0][2]*B[2][0];
		C[0][1] = A[0][0]*B[0][1] + A[0][1]*B[1][1] + A[0][2]*B[2][1];
		C[0][2] = A[0][0]*B[0][2] + A[0][1]*B[1][2] + A[0][2]*B[2][2];

		C[1][0] = A[1][0]*B[0][0] + A[1][1]*B[1][0] + A[1][2]*B[2][0];
		C[1][1] = A[1][0]*B[0][1] + A[1][1]*B[1][1] + A[1][2]*B[2][1];
		C[1][2] = A[1][0]*B[0][2] + A[1][1]*B[1][2] + A[1][2]*B[2][2];

		C[2][0] = A[2][0]*B[0][0] + A[2][1]*B[1][0] + A[2][2]*B[2][0];
		C[2][1] = A[2][0]*B[0][1] + A[2][1]*B[1][1] + A[2][2]*B[2][1];
		C[2][2] = A[2][0]*B[0][2] + A[2][1]*B[1][2] + A[2][2]*B[2][2];
	}

	//y = A.x
	static inline void MAT_VEC(double y[3], double A[3][3], double x[3]){
		y[0] = A[0][0]*x[0] + A[0][1]*x[1] + A[0][2]*x[2];
		y[1] = A[1][0]*x[0] + A[1][1]*x[1] + A[1][2]*x[2];
		y[2] = A[2][0]*x[0] + A[2][1]*x[1] + A[2][2]*x[2];
	}

	//Vector v must be a unit vector
	static inline void ANGLE(double &theta, double &phi, double v[3]){
		double l = sqrt(v[0]*v[0] + v[1]*v[1]);
		if(l == 0)
			theta = 0;
		else
			theta = acos(v[1]/l);
		if(v[0] < 0)
			theta = -theta;
		phi = acos(v[2]/MeshData::LENGTH(v));
	}

	static inline void INVERSE_U(double U[3][3]){
		double tmp;

		tmp = U[0][1]; U[0][1] = U[1][0]; U[1][0] = tmp;
		tmp = U[0][2]; U[0][2] = U[2][0]; U[2][0] = tmp;
		tmp = U[1][2]; U[1][2] = U[2][1]; U[2][1] = tmp;
	}

	static inline void GENERATE_MAT(double R[3][3], double angle, double v[3]){
		double l = sqrt(v[0]*v[0] + v[1]*v[1]);
		double a;
		if(l == 0)
			a = 0;
		else
			a = acos(v[1]/l);
		if(v[0] < 0)
			a = -a;
		double b = acos(v[2]/MeshData::LENGTH(v));

		double Rz[3][3], Rx[3][3], Rt[3][3]; 
		Rz[0][0] = Rz[1][1] = cos(a);
		Rz[0][1] = sin(a);
		Rz[1][0] = -Rz[0][1];
		Rz[2][0] = Rz[2][1] = Rz[0][2] = Rz[1][2] = 0;
		Rz[2][2] = 1;

		Rt[0][0] = Rt[1][1] = cos(angle);
		Rt[0][1] = -sin(angle);
		Rt[1][0] = -Rt[0][1];
		Rt[2][0] = Rt[2][1] = Rt[0][2] = Rt[1][2] = 0;
		Rt[2][2] = 1;

		Rx[0][0] = 1;
		Rx[0][1] = Rx[0][2] = Rx[1][0] = Rx[2][0] = 0;
		Rx[1][1] = Rx[2][2] = cos(b);
		Rx[1][2] = sin(b);
		Rx[2][1] = -Rx[1][2];

		double tmp[3][3];
		MAT_TIMES(tmp, Rz, Rx);
		MAT_TIMES(R, tmp, Rt);
		Rz[1][0] *= -1;
		Rz[0][1] *= -1;
		Rx[1][2] *= -1;
		Rx[2][1] *= -1;
		MAT_TIMES(tmp, R, Rx);
		MAT_TIMES(R, tmp, Rz);
	}

	static inline BOOL INVERSE(double B[10], double A[10]){
		double d = DET(A);
		if(fabs(d) < EPSILON2)
			return false;
		B[0] = (A[3]*A[5] - A[4]*A[4])/d;
		B[1] = (A[2]*A[4] - A[1]*A[5])/d;
		B[2] = (A[1]*A[4] - A[2]*A[3])/d;
		B[3] = (A[0]*A[5] - A[2]*A[2])/d;
		B[4] = (A[1]*A[2] - A[0]*A[4])/d;
		B[5] = (A[0]*A[3] - A[1]*A[1])/d;
		return true;
	}

	static inline double DET(double A[10]){
		return A[0]*A[3]*A[5] + 2.0*A[1]*A[4]*A[2] 
			-A[2]*A[2]*A[3] - A[1]*A[1]*A[5] - A[4]*A[4]*A[0];
	}

	static inline void MAT_BY_VEC(double v[3], double A[10], double b[3]){
		v[0] = A[0]*b[0] + A[1]*b[1] + A[2]*b[2];
		v[1] = A[1]*b[0] + A[3]*b[1] + A[4]*b[2];
		v[2] = A[2]*b[0] + A[4]*b[1] + A[5]*b[2];
	}
};

#endif // !defined(AFX_SUBDIVIDER_H__1168D90D_29F3_455D_8C32_4302007876E6__INCLUDED_)