/****************************************************************************/
/*                                                                          */
/*  VolVis is a volume visualization system for investigating, manipulating */
/*  and rendering geometric and volumetric data.                            */
/*                                                                          */
/*  Copyright (C) 1993 by the Research Foundation of the State University   */
/*                            of New York                                   */
/*                                                                          */
/*  This program is free software; you can redistribute it and/or modify    */
/*  it under the terms of the GNU General Public License as published by    */
/*  the Free Software Foundation; either version 1, or (at your option)     */
/*  any later version.                                                      */
/*                                                                          */
/*  This program is distributed in the hope that it will be useful,         */
/*  but WITHOUT ANY WARRANTY; without even the implied warranty of          */
/*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the           */
/*  GNU General Public License for more details.                            */
/*                                                                          */
/*  You should have received a copy of the GNU General Public License       */
/*  along with this program; if not, write to the Free Software             */
/*  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.               */
/*                                                                          */
/*  For information on VolVis, contact us at:                               */
/*                                                                          */
/*                volvis@cs.sunysb.edu                         (email)      */
/*                                                                          */
/*                Lisa Sobierajski & Ricardo Avila             (US Mail)    */
/*                Department of Computer Science                            */
/*                State University of New York at Stony Brook               */
/*                Stony Brook, New York  11794-4400                         */
/*                                                                          */
/****************************************************************************/




/*
 *			File: C_solve_cubic.c
 *		      Author: Rick Avila
 *			Date: 06/10/93
 *		 Description: Solution To A Cubic Equation When Coeffs
 *			      Are Real
 *	Modification History:
 *
 *		Who?		When?		Why?
 *	--------------------------------------------------------------------
 *
 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "C_volvis.h"
#include "C_raytrace.h"

test_intersect()
{
	C_FPosition		start;
	C_FVector		vector;
	C_LineIntersectInfo	intersect_info;
	int			i;

	double	c0, c1, c2, c3;	/* Coefficients Of Cubic Equation */
	double	r1, r2, r3;	/* Roots Of Equation */
	double	temp;		/* Swap Variable */
	int	num_roots;	/* Number Of Unique Roots To Equation */
	int	root;		/* Loops Through Roots */

	void			C_trilin_line_intersection();

	/* Set Ray Origin */
	start.x =  1.0;
	start.y =  1.0;
	start.z =  1.0;

	/* Set Ray Direction */
	vector.x = -0.57735;
	vector.y = -0.57735;
	vector.z = -0.57735;

	C_trilin_line_intersection( 
			&start, &vector,
			0.0, 0.0, 0.0, 0.0,	/* A, B, C, D */
			0.0, 0.0, 0.0, 1.0,	/* E, F, G, H */
			0.5, 			/* Treshold */
			&intersect_info );	/* Solutions */

printf("Solutions: %d\n", intersect_info.num_intersections );

	for( i=0; i < intersect_info.num_intersections; i++ )
	{
		printf("Solution %d: %f %f %f %f\n", i+1,
			intersect_info.local_position[i].x,
			intersect_info.local_position[i].y,
			intersect_info.local_position[i].z,
			intersect_info.local_distance[i] );
	}
}

/************************************************************************/
/* 	This routine computes the intersection(s) of a vector and an 	*/
/* isosurface within the trilinear interpolation function. The starting	*/
/* position of the vector is given in variable "start" and the direction*/
/* of the vector is given in the variable "vec". The scalar values at 	*/
/* the vertices of the [0.0 <-> 1.0] cube are supplied within variables */
/* "A"-"H" (See macro C_Trilin()). */
/*									*/
/*	Scalar Field:							*/
/*									*/
/*			C_Trilin( x, y, z, A, B, C, D, E, F, G, H )	*/
/*									*/
/* 	Parametric Line Equation:					*/
/*									*/
/* 			x = x0 + at					*/
/* 			y = y0 + bt					*/
/* 			z = z0 + ct					*/
/*									*/
/*	Isosurface Threshold Value:					*/
/*									*/
/*			iso						*/
/*									*/
/*	Intermediate Calculations:					*/
/*									*/
/*			P =  A - B - C + D				*/
/*			Q =  A - C - E + G				*/
/*			R =  A - B - E + F				*/
/*			S = -A + B + C - D + E - F - G + H		*/
/*			T = a * b * c * S				*/
/*									*/
/*	Trilinear Interpolation With Parametric Substitutions:		*/
/*									*/
/*			c0*t^3 + c1*t^2 + c2*t + c3 = 0			*/
/*									*/
/*	Where:								*/
/*									*/
/*	c0 = a*b*c*S							*/
/*									*/
/*      c1 = a*b*P + b*c*Q + a*c*R + (x0*b*c + a*(y0*c + z0*b))*S	*/
/*									*/
/*	c2 = (x0*b + y0*a)*P + (y0*c + z0*b)*Q + (x0*c + z0*a)*R +	*/
/*	     (a*y0*z0 + x0*(y0*c + z0*b))*S + 				*/
/*	     (B - A)*a + (C - A)*b + (E - A)*c				*/
/*									*/
/*	c3 = (1.0-x0-y0-z0)*A + B*x0 + C*y0 + E*z0 + 			*/
/*	      x0*y0*P + y0*z0*Q + x0*z0*R + x0*y0*z0*S - iso		*/
/*									*/
/************************************************************************/
void C_trilin_line_intersection( start, vec, A,B,C,D,E,F,G,H, iso, solution )
C_FPosition		*start;		/* Starting Point Of Parametric Line */
C_FVector		*vec; 		/* Ray Direction */
double			A, B, C, D, E, F, G, H;	/* Vertex Values */
double			iso;		/* Isosurface Threshold Value */
C_LineIntersectInfo	*solution;	/* Solution To Problem */
{
	double	c0, c1, c2, c3;	/* Coefficients Of Cubic Equation */
	double	r1, r2, r3;	/* Roots Of Equation */
	double	temp;		/* Swap Variable */
	int	num_roots;	/* Number Of Unique Roots To Equation */
	int	root;		/* Loops Through Roots */
	int	pos_dist_num;	/* Number Of Positive Distance Roots */

	double	x0, y0, z0;
	double	a, b, c;
	double	P, Q, R, S, T, U;
	double	x, y, z;
	double	dist;

	x0 = start->x;
	y0 = start->y;
	z0 = start->z;

	a  = vec->x;
	b  = vec->y;
	c  = vec->z;

	P =  A - B - C + D;
	Q =  A - C - E + G;
	R =  A - B - E + F;
	S = -A + B + C - D + E - F - G + H;
	T =  a * b * c * S;

	/* Initialize the Number Of Intersections To Zero */
	solution->num_intersections = 0;

	c0 = T;

	c1 = (a*b*P + b*c*Q + a*c*R + (x0*b*c + a*(y0*c + z0*b))*S);

	c2 = ( (x0*b + y0*a)*P + (y0*c + z0*b)*Q + (x0*c + z0*a)*R +
	       (a*y0*z0 + x0*(y0*c + z0*b))*S + 
	       (B - A)*a + (C - A)*b + (E - A)*c
	     );

	c3 = ( (1.0-x0-y0-z0)*A + B*x0 + C*y0 + E*z0 + 
		x0*y0*P + y0*z0*Q + x0*z0*R + x0*y0*z0*S - iso
	     );

	C_solve_cubic( c0, c1, c2, c3, &r1, &r2, &r3, &num_roots );

/****
printf("Equation: (%10.10f)x^3 + (%10.10f)x^2 + (%10.10f)x + (%10.10f) = 0 :\n",
	c0,c1,c2,c3);

	switch( num_roots )
	{
		case 1:
			printf("Root: %f\n", r1 );
			break;
		case 2:
			printf("Roots: %f %f\n", r1, r2 );
			break;
		case 3:
			printf("Roots: %f %f %f\n", r1, r2, r3 );
	}
****/

	/* Remove Negative Solutions And Store In Distance Array */
	pos_dist_num = 0;
	for( root=0; root < num_roots; root++ )
	{
		switch( root )
		{
			case 0:
				dist = r1;
				break;
			case 1:
				dist = r2;
				break;
			case 2:
				dist = r3;
		}

		if( dist >= 0.0 )
		{
			solution->local_distance[pos_dist_num] = dist;
			pos_dist_num += 1;
		}
	}

	solution->num_intersections = pos_dist_num;

	/* Sort The Solutions Based On Distance */
	if( pos_dist_num == 2 )
	{
		if( solution->local_distance[0] > solution->local_distance[1] )
		{
		   temp = solution->local_distance[0];
		   solution->local_distance[0] = solution->local_distance[1];
		   solution->local_distance[1] = temp;
		}
	}
	else if( pos_dist_num == 3 )
	{
		if( solution->local_distance[0] > solution->local_distance[1] )
		{
		   temp = solution->local_distance[0];
		   solution->local_distance[0] = solution->local_distance[1];
		   solution->local_distance[1] = temp;
		}
		if( solution->local_distance[1] > solution->local_distance[2] )
		{
		   temp = solution->local_distance[1];
		   solution->local_distance[1] = solution->local_distance[2];
		   solution->local_distance[2] = temp;
		}
		if( solution->local_distance[0] > solution->local_distance[1] )
		{
		   temp = solution->local_distance[0];
		   solution->local_distance[0] = solution->local_distance[1];
		   solution->local_distance[1] = temp;
		}
	}

	for( root=0; root < solution->num_intersections; root++ )
	{
		/**********************************************/
		/* Determine The (x,y,z) Position Of Solution */
		/**********************************************/
		x = x0 + a * solution->local_distance[root];
		y = y0 + b * solution->local_distance[root];
		z = z0 + c * solution->local_distance[root];

		solution->local_position[root].x = x;
		solution->local_position[root].y = y;
		solution->local_position[root].z = z;
	}
}

/************************************************************/
/* Solves A Cubic Equation When c0, c1, c2, And c3 Are REAL */
/* Solution Taken From Numerical Recipes In C 2nd Ed        */
/************************************************************/
C_solve_cubic( c0, c1, c2, c3, r1, r2, r3, num_roots )
double	c0, c1, c2, c3;	/* Coefficients To Cubic Equation */
double	*r1, *r2, *r3;	/* Ptrs To Roots */
int	*num_roots;	/* Number Of Roots */
{
	double	Q, R;
	double	R_squared;	/* R*R */
	double	Q_cubed;	/* Q*Q*Q */
	double	theta;
	double	A, B;

	void	C_solve_quadratic();

	/*************************************************/
	/*                     3       2                 */
	/* Cubic Equation: c0*t  + c1*t  + c2*t + c3 = 0 */
	/*                                               */
	/*************************************************/

	/* Make Sure This Is A Bonafide Cubic Equation */
	if( c0 != 0.0 )
	{
	   /* Put Coefficients In Right Form */
	   c1 = c1/c0;
	   c2 = c2/c0;
	   c3 = c3/c0;

	   Q = ((c1*c1) - 3*c2)/9.0;

	   R = (2.0*(c1*c1*c1) - 9.0*(c1*c2) + 27.0*c3)/54.0;

	   R_squared = R*R;
	   Q_cubed   = Q*Q*Q;

	   if( R_squared <= Q_cubed )
	   {
		if( Q_cubed == 0.0 )
		{
			*r1 = -c1/3.0;
			*r2 = *r1;
			*r3 = *r1;
			*num_roots = 1;
		} 
		else
		{
		   theta = acos( R / (sqrt(Q_cubed) ) );

		   *r1 = -2.0*sqrt(Q)*cos( theta/3.0 ) - c1/3.0;
		   *r2 = -2.0*sqrt(Q)*cos( (theta + 2.0*M_PI)/3.0 ) - c1/3.0;
		   *r3 = -2.0*sqrt(Q)*cos( (theta - 2.0*M_PI)/3.0 ) - c1/3.0;

		   *num_roots = 3;

		   /*********************************/
		   /* Reduce Number Of Roots To Two */
		   /*********************************/
		   if( *r1 == *r2 )
		   {
			*num_roots = 2;
			*r2 = *r3;
		   }
		   else if( *r1 == *r3 )
		   {
			*num_roots = 2;
		   }

		   if( (*r2 == *r3) && (*num_roots == 3) )
		   {
			*num_roots = 2;
		   }

		   /*********************************/
		   /* Reduce Number Of Roots To One */
		   /*********************************/
		   if( (*r1 == *r2) )
		   {
			*num_roots = 1;
		   }
		}
	   }
	   else 
	   {
		A = -C_Sign(R) * C_cbrt(fabs(R) + sqrt(R_squared - Q_cubed));

		if( A == 0.0 )
			B = 0.0;
		else
			B = Q/A;

		*r1 =  (A + B) - c1/3.0;
		*r2 = *r1;
		*r3 = *r1;

		*num_roots = 1;
	   }
	}
	else
	{
		/*********************************************/
		/*                         2                 */
		/* Quadratic Equation: c1*t  + c2*t + c3 = 0 */
		/*                                           */
		/*********************************************/

		/* Okay This Was Not A Cubic - Lets Try Quadratic?? */
		C_solve_quadratic( c1, c2, c3, r1, r2, num_roots );
	}
}

/************************************************************/
/* Solves A Quadratic Equation When c1, c2, And c3 Are REAL */
/* Solution Taken From Numerical Recipes In C 2nd Ed        */
/************************************************************/
void C_solve_quadratic( c1, c2, c3, r1, r2, num_roots )
double	c1, c2, c3;	/* Coefficients To Cubic Equation */
double	*r1, *r2;	/* Ptrs To Roots */
int	*num_roots;	/* Number Of Roots */
{
	double	Q;
	double	determinant;

	/*********************************************/
	/*                         2                 */
	/* Quadratic Equation: c1*t  + c2*t + c3 = 0 */
	/*                                           */
	/*********************************************/

	/* Make Sure This Is A Bonafide Cubic Equation */
	if( c1 != 0.0 )
	{
		determinant = c2*c2 - 4*c1*c3;

		if( determinant >= 0.0 )
		{
			Q = -0.5 * (c2 + C_Sign(c2)*sqrt(determinant));

			*r1 = Q / c1;

			if( Q == 0.0 )
				*r2 = 0.0;
			else
				*r2 = c3 / Q;

			*num_roots = 2;

			/*********************************/
			/* Reduce Number Of Roots To One */
			/*********************************/
			if( *r1 == *r2 )
			{
				*num_roots = 1;
			}
		}
		else	/* Equation Does Not Have Real Roots */
		{
			*num_roots = 0;
		}
	}
	else
	{
		/**********************************/
		/*                                */
		/* Linear Equation: c2*t + c3 = 0 */
		/*                                */
		/**********************************/

		/* Now This Had Better Be Linear */
		if( c2 != 0.0 )
		{
			*r1 = -c3 / c2;

			*num_roots = 1;
		}
		else
		{
			*num_roots = 0;
		}
	}
}
