#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <GL/glut.h>
#include <windows.h>
#include "glf.h"

#define MAXPOINT 100
#define DIVIDENUMBER 100
#define BARYCENTRIC 1
#define SPLINE 3
#define CLEAR 9
#define CUBICSPLINE 5

double	**dmatrix(int nrl,int nrh,int ncl,int nch);
int		inverse(double *m);
int     gaussj(double **a, int n,double **b, int m);
void	free_dmatrix(double **m,int nrl,int nrh,int ncl,int nch);

typedef struct Lines // 좌표를 담을 struct
{
	double x;
	double y;
}Line;

typedef struct BaryLine // 알파,알파2,1의 값을 넣을 struct
{
	double alpha;
	double beta;
	double gamma;
}BaryLine;

typedef struct spline{
	double a,b,c,d;
}spline;

Line temp[MAXPOINT]; // Hermite 전용 
Line Spline_Line[MAXPOINT]; // Spline 전용
Line bary[3];        // Barycentric 전용
BaryLine bary_line;
spline answer[MAXPOINT];
char print[256];

int point_number = 0 , arrayIndex = -1 , display_mode = -1 , bary_number = 0 , bary_check = 0, Spline_number = 1, SplineIndex = -1;
GLsizei width = 800 , height = 600;
double interpol(int N, Line P[], spline F[], double x);
void   calcSpline(int N, Line P[], spline F[]);
double interpol(int N, Line P[], spline F[], double x);
void   calcSpline(int N, Line P[], spline F[]);
static void  tridia(int n, double M[][4], double S[]); /* tridiagonal matrix */

void myCaption()
{
	float	cscale=14., space=110., sspace=90., lm=80., um=height-30.;
	
	glColor3f(0.3, 0.5, 0.1);
	glPushMatrix();
	glTranslatef(400, um, 0.);
	glScalef(cscale, cscale, 1.);
	glfDrawSolidString(print);
	glPopMatrix();
}

void myBarycentric()
{
	double **M,**I,a1,a2,a3;
	a1 = a2 = a3 = 0;
	M = dmatrix(1,3,1,3); // 3x3 행렬 완성 -> 행렬에 값 넣기
	I = dmatrix(1,3,1,3); // 3x3 행렬 완성 -> 행렬에 값 넣기 (버리는 행렬)
	// 	M[1][1] = bary[0].x;  // 행렬값 추가 (3x3 값)
	// 	M[2][1] = bary[1].x;
	// 	M[3][1] = bary[2].x;
	// 	M[1][2] = bary[0].y;
	// 	M[2][2] = bary[1].y;
	// 	M[3][2] = bary[2].y;
	// 	M[1][3] = M[2][3] = M[3][3] = 1;
	M[1][1] = bary[0].x;  // 행렬값 추가 (3x3 값)
	M[1][2] = bary[1].x;
	M[1][3] = bary[2].x;
	M[2][1] = bary[0].y;
	M[2][2] = bary[1].y;
	M[2][3] = bary[2].y;
	M[3][1] = M[3][2] = M[3][3] = 1;
	
	//bary_line.alpha;  알파  bary_line.beta;  베타 bary_line.gamma;  1
	gaussj(M,3,I,3); // 역행렬을 M에 담아옴
	printf(" 역 행 렬 값 \n");
	for(int i = 1 ; i <= 3 ; i++)
	{
		for(int j = 1 ; j <= 3 ; j++)
			printf("%5.2lf \t",M[i][j]);
		printf("\n");
	}
	// 역행렬과 a1,a2,a3 값을 
	// 	a1 = (M[1][1] * bary_line.alpha) + (M[2][1] * bary_line.beta) + (M[3][1] * bary_line.gamma);
	// 	a2 = (M[1][2] * bary_line.alpha) + (M[2][2] * bary_line.beta) + (M[3][2] * bary_line.gamma);
	// 	a3 = (M[1][3] * bary_line.alpha) + (M[2][3] * bary_line.beta) + (M[3][3] * bary_line.gamma);
	a1 = (M[1][1] * bary_line.alpha) + (M[1][2] * bary_line.beta) + (M[1][3] * bary_line.gamma);
	a2 = (M[2][1] * bary_line.alpha) + (M[2][2] * bary_line.beta) + (M[2][3] * bary_line.gamma);
	a3 = (M[3][1] * bary_line.alpha) + (M[3][2] * bary_line.beta) + (M[3][3] * bary_line.gamma);
	printf("\n알파 : %.2lf , 알파2 : %.2lf , 알파3 : %.2lf \n",bary_line.alpha,bary_line.beta,bary_line.gamma);
	sprintf(print,"Barycentric Coordinate = ( %.2lf , %.2lf , %.2lf )",a1,a2,a3);
	bary_check = 1;
	free_dmatrix(M, 1, 3, 1, 3);
	free_dmatrix(I, 1, 3, 1, 3);
	glutPostRedisplay();
}

void mySpline()
{
	if(Spline_number > 2) // 2개의 점 이상부터 계산 (1구간 이상부터)
		calcSpline(Spline_number,Spline_Line,answer); // 찍은 점에 대한 spline 함수 계산
	
    //Spline_number : 점의 갯수
	//Spline_Line[i] : 찍은 점의 좌표(x,y) 배열 
	//answer : f(x) = ax^3 + bx^2 + cx + d, 각 구간의 계수 (a,b,c,d)
	
	double calc_x, calc_y, t;
	
	glColor3f(.5f,.1f,.3f);
	glBegin(GL_LINE_STRIP);
	if(Spline_number > 2) // 2개의 점 이상부터 계산
	{
		for(int i = 1 ; i < Spline_number ; i++) 
		{
			calc_x = Spline_Line[i].x;
			t = (Spline_Line[i+1].x - calc_x)/DIVIDENUMBER; // 시간 나눔
			
			while(calc_x < Spline_Line[i+1].x) // 다음 점에 도달할때까지 루프
			{
				calc_y = interpol(Spline_number,Spline_Line,answer,calc_x); // y값 계산
				glVertex2d(calc_x, calc_y); // 선 긋기
				calc_x += t; // 시간 추가
			}
			//	glVertex2d(Spline_Line[i+1].x,Spline_Line[i+1].y);
		}
	}
	glEnd();
}

void myHermite()
{
	double left_dx, left_dy, right_dx, right_dy, t, calc_x, calc_y;
	
	glColor3f(.5f,.1f,.3f);
	glBegin(GL_LINE_STRIP);
	for(int i = 0 ; i < point_number - 1 ; i++) // 0부터 전체 입력 좌표
	{
		for(int j = 0 ; j < DIVIDENUMBER ; j++) // 0부터 DIVIDENUMBER의 양에 좌우되는 곡선 표기
		{
			t = (double)j / DIVIDENUMBER;
			
			if(i == 0) // 처음의 변화량
			{
				left_dx = temp[i+1].x - temp[i].x;
				left_dy = temp[i+1].y - temp[i].y;
			}	
			else // 자신의 점에 대한 변화량
			{
				left_dx = temp[i+1].x - temp[i-1].x;
				left_dy = temp[i+1].y - temp[i-1].y;
			}
			
			if(i == point_number - 2)  // 마지막의 변화량을 계산
			{
				right_dx = temp[i+1].x - temp[i].x;
				right_dy = temp[i+1].y - temp[i].y;
			}
			else // 오른쪽 점에 대한 변화량
			{
				right_dx = temp[i+2].x - temp[i].x;
				right_dy = temp[i+2].y - temp[i].y;
			}
			// temp[i].x,temp[i+1].x, temp[i].y,temp[i+1].y
			calc_x = (( 2 * t * t * t)-(3 * t * t) + 1) * temp[i].x +   // x의 첫번째 행 f(x)
				(( t * t * t    )-(2 * t * t) + t) * left_dx +     // x의 두번째 행 f'(x)
				((-2 * t * t * t)+(3 * t * t)    ) * temp[i+1].x + // x의 세번째 행 f(x+1)
				(( t * t * t    )-(t * t)        ) * right_dx;     // x의 네번째 행 f'(x+1)
			
			calc_y = (( 2 * t * t * t)-(3 * t * t) + 1) * temp[i].y +
				(( t * t * t    )-(2 * t * t) + t) * left_dy  +
				((-2 * t * t * t)+(3 * t * t)    ) * temp[i+1].y +
				(( t * t * t    )-(t * t)        ) * right_dy;	
			
			//printf("debug %d 번째 출력 calc_x = %d, calc_y = %d\n",j,calc_x,calc_y);
			glVertex2d(calc_x, calc_y); // 계산 결과를 출력
		}
	}
	glEnd();
}

void myDisplay() // 전체 display
{
	glClear(GL_COLOR_BUFFER_BIT);
	glClearColor (0.85, 0.85, 0.85, 1.0);
	glPointSize(10.0f);
	//myCaption();
	if(display_mode == -1)
	{
		if(bary_number != 0)
		{
			for(int i = 0 ; i < bary_number ; i++)
				bary_line.alpha = bary_line.beta = bary_line.gamma = 0;
			bary_number = 0;
			bary_check = 0;
		}
		glBegin(GL_POINTS);        // 버텍스 
		for(int i = 0 ; i < point_number ; i++) 
		{
			if(i == arrayIndex)    // 선택, 마지막 선택 좌표 색 표기
				glColor3f(1.0f,0.0f,0.0f);
			else
				glColor3f(.3f,.5f,.1f);
			glVertex2i(temp[i].x,temp[i].y);
		}
		glEnd();
		
		glColor3f(0.1f,0.5f,0.3f); // 엣지 그리기
		glBegin(GL_LINE_STRIP);  
		for(i = 0 ; i < point_number ; i++)
			glVertex2d(temp[i].x,temp[i].y);
		glEnd();
	}
	else if(display_mode == 1) // barycentric coordinates
	{
		if(point_number != 0)
		{
			for(int i = 0 ; i < point_number ; i++)
				temp[i].x = temp[i].y = 0;
			point_number = 0;		
		}
		if(Spline_number != 1)
		{
			for(int i = 1 ; i < Spline_number ; i++)
				Spline_Line[i].x = Spline_Line[i].y = 0;
			Spline_number = 1;
		}

		glBegin(GL_POINTS);        // 버텍스 
		for(int i = 0 ; i < bary_number ; i++) 
		{
			glColor3f(.3f,.5f,.1f);
			glVertex2i(bary[i].x,bary[i].y);
		}
		glEnd();
		if(bary_number == 3)
		{
			glColor3f(0.1f,0.5f,0.3f); 
			glBegin(GL_TRIANGLES);  // 삼각형 그리기
			for(int i = 0 ; i < bary_number; i++)
				glVertex2d(bary[i].x,bary[i].y);
			glEnd();
		}
	}
	if(display_mode == 3) // Piecewise cubic hermite interpolation
	{
		if(bary_number != 0)
		{
			for(int i = 0 ; i < bary_number ; i++)
				bary_line.alpha = bary_line.beta = bary_line.gamma = 0;
			bary_number = 0;
			bary_check = 0;
		}
		if(Spline_number != 1)
		{
			for(int i = 1 ; i < Spline_number ; i++)
				Spline_Line[i].x = Spline_Line[i].y = 0;
			Spline_number = 1;
		}
		glClear(GL_COLOR_BUFFER_BIT);
		glBegin(GL_POINTS);        // 버텍스 
		for(int i = 0 ; i < point_number ; i++) 
		{
			if(i == arrayIndex)    // 선택, 마지막 선택 좌표 색 표기
				glColor3f(1.0f,0.0f,0.0f);
			else
				glColor3f(.3f,.5f,.1f);
			glVertex2i(temp[i].x,temp[i].y);
		}
		glEnd();
		
		glColor3f(0.1f,0.5f,0.3f); // 엣지 그리기
		glBegin(GL_LINE_STRIP);  
		for(i = 0 ; i < point_number ; i++)
			glVertex2d(temp[i].x,temp[i].y);
		glEnd();
		
		myHermite();
	}
	if(display_mode == 5) // Spline
	{
		if(bary_number != 0)
		{
			for(int i = 0 ; i < bary_number ; i++)
				bary_line.alpha = bary_line.beta = bary_line.gamma = 0;
			bary_number = 0;
			bary_check = 0;
		}
		if(point_number != 0)
		{
			for(int i = 0 ; i < point_number ; i++)
				temp[i].x = temp[i].y = 0;
			point_number = 0;		
		}

		glClear(GL_COLOR_BUFFER_BIT);
		glBegin(GL_POINTS);  // 버텍스 
		for(int i = 1 ; i < Spline_number ; i++) 
		{
			if(i == SplineIndex)    // 선택, 마지막 선택 좌표 색 표기
				glColor3f(1.0f,0.0f,0.0f);
			else
				glColor3f(.3f,.5f,.1f);
					
			glVertex2i(Spline_Line[i].x,Spline_Line[i].y);
		}
		glEnd();

		glColor3f(0.1f,0.5f,0.3f); // 엣지 그리기
		glBegin(GL_LINE_STRIP);  
		for(i = 1 ; i < Spline_number ; i++)
			glVertex2d(Spline_Line[i].x,Spline_Line[i].y);
		glEnd();
	
		mySpline();
	}
	if(display_mode == 9) // 초기화
	{
		glClear(GL_COLOR_BUFFER_BIT);
		for(int i = 0 ; i < point_number ; i++)
			temp[i].x = temp[i].y = 0;
		point_number = 0;
		display_mode = -1;
		for(i = 0 ; i < bary_number ; i++)
			bary_line.alpha = bary_line.beta = bary_line.gamma = 0;
		bary_number = 0;
		bary_check = 0;
		for(i = 1 ; i < Spline_number ; i++)
			Spline_Line[i].x = Spline_Line[i].y = 0;
		Spline_number = 1;
	}
	if(bary_check == 1)
		myCaption();
	glFlush();
	glutSwapBuffers();
}

void myMouse(int btn, int state, int x, int y) // 마우스 이벤트
{
	y = height-y;
	if(btn == GLUT_LEFT_BUTTON && state == GLUT_DOWN && arrayIndex == -1 && display_mode != 1 && display_mode != 5) // 클릭할때마다 좌표 저장
	{  
		temp[point_number].x = x;
		temp[point_number].y = y;
		point_number++;
		printf("%d 번째 좌표 입력 x = %d , y = %d\n",point_number,x,y); // 좌표 입력에 대한 안내
	}
	if(btn == GLUT_LEFT_BUTTON && state == GLUT_DOWN && display_mode == 1 && bary_number != 3) // 클릭할때마다 좌표 저장
	{
		bary[bary_number].x = x;
		bary[bary_number].y = y;
		bary_number++;
		printf("barycentric %d 번째 좌표 입력 x = %d , y = %d\n",bary_number,x,y); // 좌표 입력에 대한 안내
	}
	if(btn == GLUT_LEFT_BUTTON && state == GLUT_DOWN && SplineIndex == -1 && display_mode == 5) // 클릭할때마다 좌표 저장
	{
		Spline_Line[Spline_number].x = x;
		Spline_Line[Spline_number].y = y;
		printf("Spline %d 번째 좌표 입력 x = %d , y = %d\n",Spline_number,x,y); // 좌표 입력에 대한 안내
		Spline_number++;
	}
	myDisplay(); // 갱신 호출 
}

void myMotion(int x, int y) // 좌표 이동에 대한 갱신
{
	y = height-y;
	if(display_mode == 3)
	{
		temp[arrayIndex].x = x;
		temp[arrayIndex].y = y;
	//printf("옮기는 인덱스 %d , 좌표 : x = %d , y = %d \n",arrayIndex,x,y);
	}
	if(display_mode == 5)
	{
		Spline_Line[SplineIndex].x = x;
		Spline_Line[SplineIndex].y = y;
	//printf("옮기는 인덱스 %d , 좌표 : x = %d , y = %d \n",arrayIndex,x,y);
	}
	glutPostRedisplay();

}

void myPassiveMotion(int x, int y) // 찍은 좌표의 범위에 포함되는지 확인
{
	y = height-y;
	//printf("이동좌표 : %d \t %d\n",x,y);
	if(display_mode != 1)
	{
		if(display_mode == 3)
		{
			for(int i = 0 ; i < point_number ; i++)
			{
				if(temp[i].x-10 <= x && temp[i].x+10 >= x && temp[i].y-10 < y && temp[i].y+10 >= y) // 범위 포함
				{
					arrayIndex = i; // 그 점의 인덱스 설정 (선택 위치)
					myDisplay();
					return;
				}
			}
			arrayIndex = -1; // 없다면 원래대로 
		}
		if(display_mode == 5)
		{
			for(int i = 1 ; i < Spline_number ; i++)
			{
				if(Spline_Line[i].x-10 <= x && Spline_Line[i].x+10 >= x && Spline_Line[i].y-10 < y && Spline_Line[i].y+10 >= y) // 범위 포함
				{
					SplineIndex = i; // 그 점의 인덱스 설정 (선택 위치)
					myDisplay();
					return;
				}
			}
			SplineIndex = -1; // 없다면 원래대로 
		}
	}

	if(display_mode == 1 && bary_number == 3) // barycentric mode
	{
		bary_line.alpha = x;
		bary_line.beta = y;
		bary_line.gamma = 1;
		myBarycentric();
		
	}
}

void right_menu(int id) // 우클릭 메뉴
{
	if(id == 1)
	{
		display_mode = BARYCENTRIC;
		arrayIndex = -1;
		glutPostRedisplay();
	}
	else if(id == 3)
	{
		display_mode = SPLINE;
		arrayIndex = -1;
		glutPostRedisplay();
	}
	else if(id == 5)
	{
		display_mode = CUBICSPLINE;
		arrayIndex = -1;
		glutPostRedisplay();
	}
	
	else if(id == 9)
	{
		display_mode = CLEAR;
		arrayIndex = -1;
		glutPostRedisplay();
	}
}

void myReshape(GLsizei w, GLsizei h) 
{
	glClearColor (0.85, 0.85, 0.85, 1.0);
	glViewport(0,0,w,h);
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity(); 
	glOrtho(0.0, (GLdouble)w, 0.0, (GLdouble)h, -1.0, 1.0);
	glMatrixMode(GL_MODELVIEW); 
	glLoadIdentity(); 
	glLineWidth(1.5);
	width = w;
	height = h; 
	
}

void calcSpline(int n, Line P[], spline F[]) // n개의 점이 주어질 때, 각 구간의 spline함수 answer를 구함
{ 
	int   i;
	double (*M)[4];
	double *S;
	double *H;
	double *D;
	
	M = (double (*)[4])calloc(n,sizeof(double[4])); // 1 ~ n-1
	D = (double *)calloc(n,sizeof(double));         // 0 ~ n-1
	S = (double *)calloc(n+1,sizeof(double));       // 0 ~ n
	H = (double *)calloc(n,sizeof(double));         // 0 ~ n
	
	// H[i] = 구간 수 n개
	for(i=0; i < n; i++)
	{
		H[i] = P[i+1].x - P[i].x;
	}
	
	// D[i] = (P[i+1].y-P[i].y)/H[i], i = 0 ~ n-1
	for(i=0; i < n; i++)
	{
		D[i] = (P[i+1].y - P[i].y)/H[i];
	}
	
	// M[i][0 ~ 2] = { H[i-1], 2 * (H[i-1]+H[i]), H[i] }, i = 2 ~ n-2
	for(i=2; i < n-1; i++) 
	{ 
		M[i][0] = H[i-1];
		M[i][1] = 2.*(H[i-1]+H[i]);
		M[i][2] = H[i];
	}
	
	// M[i][3] = 6*(D[i+1] - D[i]), i=1,...,n-1
	for(i=1;i < n;i++)
	{
		M[i][3] = 6.0*(D[i] - D[i-1]);
	}
	
	S[0] = S[n]  = 0.0;  // 양 끝점의 2차미분 값 
	
	// M[1][0~2] = {0, 2*(H[0]+H[1]), H[1]}, i=1
	M[1][0] = 0.;                   
	M[1][1] = 2.0*(H[0]+H[1]);      
	M[1][2] = H[1];                 
	
	// M[n-1][0~2] = { H[n-2]-H[n-1]*H[n-1]/H[n-2], 2*H[n-2]+3*H[n-1]+H[n-1]*H[n-1]/H[n-2], 0 } i=n-1
	M[n-1][0] = H[n-2]-H[n-1]*H[n-1]/H[n-2];               
	M[n-1][1] = 2.0*H[n-2]+3.*H[n-1]+H[n-1]*H[n-1]/H[n-2]; 
	M[n-1][2] = 0.; 
	
	// 2차 미분 값을 구함 S[1],...,S[n-1]을 구함.
	tridia(n-1, M+1, S+1);
	
	// n개 구간의 함수의 계수 계산
	for(i = 0 ; i < n ; i++)
	{
		F[i].a = (S[i+1]-S[i])/(6.0*H[i]);
		F[i].b = S[i]/2.0;
		F[i].c = (P[i+1].y-P[i].y)/H[i] - H[i]*(2.0*S[i]+S[i+1])/6.0;
		F[i].d = P[i].y;
	}
	
	free(D);	free(H);	free(S);	free(M);
}

void tridia(int n, double M[][4], double S[]) // 2차미분 함수
{
	int i;
	
	//Modify the coefficients.
	M[0][2] = M[0][2]/M[0][1];                      //Division by zero risk.
	M[0][3] = M[0][3]/M[0][1];
	
	double id;
	for(i = 1; i != n; i++)
	{
		id = 1.0/(M[i][1] - M[i-1][2]*M[i][0]);      //Division by zero risk.
		M[i][2] = M[i][2] * id;                         //Last value calculated is redundant.
		M[i][3] = (M[i][3] - M[i][0]*M[i-1][3])*id;
	}
	
	//Now back substitute.
	S[n-1] = M[n-1][3];
	for(i = n - 2; i != -1; i--)
	{
		S[i] = M[i][3] - M[i][2]*S[i+1];
	}
}

double interpol(int n, Line P[], spline F[], double x)
{ 
	int   i;
	double fx, xi;
	
	for(i = 0 ; i < n ; i++)
	{ 
		if((x >= P[i].x) && (x < P[i+1].x)) 
		{ 
			xi = x - P[i].x;
			fx = xi * (xi * (xi * F[i].a + F[i].b) + F[i].c) + F[i].d;
			return fx;
		}
	}
	return 0;
}

void main(int argc, char** argv)
{
	glutInit(&argc,argv);
	glutInitDisplayMode(GLUT_DOUBLE|GLUT_RGB);  
	glutInitWindowPosition(300,200);
	glutInitWindowSize(width,height); 
	glutCreateWindow("Homework 1 - 20021546 Yeon Jae-Hyuk"); 
	glfInit();
	glfLoadFont("fonts/garamond1.glf");
	glutReshapeFunc(myReshape);
	glutMotionFunc(myMotion);
	glutPassiveMotionFunc(myPassiveMotion);
	glutDisplayFunc(myDisplay);
	glutMouseFunc(myMouse);
	glutCreateMenu(right_menu);
	glutAddMenuEntry("Barycentric Coodinates",1);
	glutAddMenuEntry("Piecewise cubic hermite interpolation",3);
	glutAddMenuEntry("Cubic Spline",5);
	glutAddMenuEntry("Clear",9);
	glutAttachMenu(GLUT_RIGHT_BUTTON);
	glutMainLoop(); 
}