/*!
	\file cuttingPlane.cpp
	\brief Create and draws an arbitrary cutting plane.	
	\author Melvin Quintos
	\date October 2002

	We use a partial implementation of the "Marching Cubes" algorithm to draw an arbitrary
	cutting plane with normal N.	
*/

#include "stdheader.h"
#include <GL/glut.h>
#include <math.h>
#include <stdio.h>

//! Calculates if point Q lies on plane with normal N
/*!
	\sa calculateCuttingPlane()
	\param q1 \a q2 \a q3 Coordinates of point Q, Q=<q1,q2,q3>
	\param n1 \a n2 \a n3 Components of normal vector, N=<n1,n2,n3>

	This function calculates if point Q lies on plane N and returns 0 if it is on the plane,
	<0 if it is less than, >0 if greater than.
*/
inline float planeEq(int q1, int q2, int q3,  // point Q = (q1,q2,q3)
					 float n1, float n2, float n3)  // normalized normal to plane, N = (n1,n2,n3)
{
	// N * < Q - P> = 0
	//
	// N = (n1,n2,n3)
	// Q = (q1,q2,q3)
	// P = g_fPx, g_fPy, g_fPz;
	//
	return n1*(q1-g_fPx) + n2*(q2-g_fPy) + n3*(q3-g_fPz);
}

//! Calculates cutting plane
void calculateCuttingPlane()
{	
	// calculate normalized normal vector to plane
	float n1, n2, n3;

	// calculate any transformations to the plane
	float rotMat[16];
	glMatrixMode(GL_MODELVIEW);
	glPushMatrix();
	glLoadIdentity();
	glRotatef(g_fNormalRotateX, 1,0,0);	
	glRotatef(g_fNormalRotateZ, 0,0,-1);
	glGetFloatv(GL_MODELVIEW_MATRIX,rotMat);
	glPopMatrix();
	g_fNx = rotMat[1]+rotMat[3];
	g_fNy = rotMat[5]+rotMat[9];
	g_fNz = rotMat[9]+rotMat[11];

	
	double val = sqrt(g_fNx*g_fNx + g_fNy*g_fNy + g_fNz*g_fNz);

	if(val == 0)
		val=1;
	n1 = g_fNx/val;
	n2 = g_fNy/val;
	n3 = g_fNz/val;

	for(int i=0; i < g_Header.depth; i++)
		for(int j=0; j < g_Header.height; j++)
			for(int k=0; k < g_Header.width; k++)
				g_bData[i][j][k] = (planeEq(k,j,i, n1, n2, n3) < 0.0000) ? false: true;

	for(i=0; i < g_Header.depth-1; i++)
		for(int j=0; j < g_Header.height-1; j++)
			for(int k=0; k < g_Header.width-1; k++)
			{
				// clear old value
				g_ucCell[i][j][k] = 0;
				// calculate new value
				if(g_bData[i][j][k])		g_ucCell[i][j][k] |= 0x1;
				if(g_bData[i][j][k+1])		g_ucCell[i][j][k] |= 0x2;
				if(g_bData[i][j+1][k+1])	g_ucCell[i][j][k] |= 0x4;
				if(g_bData[i][j+1][k])		g_ucCell[i][j][k] |= 0x8;
				if(g_bData[i+1][j][k])		g_ucCell[i][j][k] |= 0x10;
				if(g_bData[i+1][j][k+1])	g_ucCell[i][j][k] |= 0x20;
				if(g_bData[i+1][j+1][k+1])	g_ucCell[i][j][k] |= 0x40;
				if(g_bData[i+1][j+1][k])	g_ucCell[i][j][k] |= 0x80;				
			}

}