/*
        Equation solving for Triangle Transformations
*/
#include <float.h>
#include "fdestria.h"
#include "fdesign.h"
#include "fdesfile.h"
#include "fdesmenu.h"
#include "fdesmous.h"
#include "fdesplot.h"


/************************************************************************
        SOLVING EQUATIONS FOR TRIANGLE TRANSFORMATION TO IFS CODE
************************************************************************/
/*
	Compute the IFS array for the given triangles
*/
float det(float a, float b, float c, float d)
{
	return(a*d-b*c);
}
float solve3(float x1,float x2,float x1h,float y1,float y2,float y1h,
	     float z1,float z2,float z1h, float *a,float *b,float *e)	/* solve linear system */
{
	/*  Format for the equations is this:
		x1*a + x2*b + e = x1h
		y1*a + y2*b + e = y1h
		z1*a + z2*b + e = z1h
	*/
float det1;
	det1 = x1 * det(y2,1.0,z2,1.0) - x2 * det(y1,1.0,z1,1.0) + 1*det(y1,y2,z1,z2);
	if (det1 == 0.0) return (det1);
	*a = (x1h * det(y2,1.0,z2,1.0) - x2 * det(y1h,1.0,z1h,1.0) + 1*det(y1h,y2,z1h,z2))/det1;
	*b = (x1 * det(y1h,1.0,z1h,1.0) - x1h * det(y1,1.0,z1,1.0) + 1*det(y1,y1h,z1,z1h))/det1;
	*e = (x1 * det(y2,y1h,z2,z1h) - x2 * det(y1,y1h,z1,z1h) + x1h*det(y1,y2,z1,z2))/det1;
	return(det1);
}

/*
	compute IFS of triangle and place codes at IFS pointer
	t1 is the transform triangle, t0 is the reference triangle
*/
void IFS_compute(float *IFS, triangle *t1, triangle *t0)
{
        solve3((*t0).col[0],(*t0).row[0],(*t1).col[0],
               (*t0).col[1],(*t0).row[1],(*t1).col[1],
               (*t0).col[2],(*t0).row[2],(*t1).col[2],
	       &IFS[0],&IFS[1],&IFS[4]);
        solve3((*t0).col[0],(*t0).row[0],(*t1).row[0],
               (*t0).col[1],(*t0).row[1],(*t1).row[1],
               (*t0).col[2],(*t0).row[2],(*t1).row[2],
	       &IFS[2],&IFS[3],&IFS[5]);
	}
/*
	compute IFS codes for all of the triangles
*/
void IFS_compute_all(int howmany, triangle *t, triangle *t0)
{
int i;
float total_area;
	IFS[0] = howmany;
	total_area = 0.0;
	for (i=0; i<howmany; i++)
	{
		IFS_compute(&IFS[1+i*7],&t[i],t0);
		/* pseudo-probability measure */
		total_area += (IFS[1+i*7+6] = triangle_area(&t[i]));
	}
	/* normalize probabilities to 1.0 */
	for (i=0; i<howmany; i++)
	{
		IFS[1+i*7+6] /= total_area;
	}
}