/**********************************************************************
Algorithm for getting 3-D coordinates 
from 2-D coordinates in mirror symmetric model

Copyright (C) 1993 Shusaku Furushima.

Caution: This source code cannot be executed itself.  Please understand 
the concepts of this algorithm and modify this source code to fit your 
products.

This source code is copyrighted, but you can freely use all or a part of 
source code for your products (commercial use also allowed).
**********************************************************************/

#include <stdio.h>
#include <math.h>

#define real double
	/* max point of qr analysis */
#define RQRMP 10
#define RQRNP 30

/******** Rcoord ********/
typedef struct {
	real x;
	real y;
	real z;
} Rcoord;

/******** 4 x 4 matrix **********/
typedef real Rmat[4][4];

/********** mathematical functions ********/
extern void Rqranalysis(int m, int n, real qa[][RQRNP],
	                real qy[], real qz[]);

/******** general 4x4 matrix *********/
class Rmat44{
	void initmat(Rmat amat);
  public:
	void matmul(Rmat mat2);
	Rmat mat;
	void IRmat44();
	void rotp(real p);
	void rotq(real p);
	void rotr(real p);
	void copyfrom(Rmat mat2);
	void move(real x, real y, real z);
	void bepers();
	void mul(real x, real y, real z);
	void getcoord2(Rcoord *coord, Rcoord *cd, int pers);
};

/********* 4x4 matrix for display objects *********/
class Rmatobj{
	public:
	Rmat44 *devmat;			/* matrix of [3D -> 2D screen] */
	Rmat44 *symm;			/* symmetrical matrix of devmat */
	real pp, pq, pr;		/* rotate angle */
	real dx, dy, dz;		/* move coordinate */
	real bx, by;			/* scale to screen 2d coord */
	int	pers;				/* flag if pers:1 or ortho:0 */

	void IRmatobj();					/* init matrixes */
	void makedev();						/* make devmat */
	void setsymmatx();					/* set symm matrix of devmat */
	Rcoord getdev(Rcoord *pt);			/* get screen coord */
	Rcoord calc3d(Rcoord *pt1, Rcoord *pt2); /* get 3D coord from 2D coord */
};



/********** least square method by QR (without pivoting) ***********/
/********* qranalysis simple version ********/
/*   m    number of unknown param.        */
/*   n    number of data sets.            */
/*   qa   data param matrix.              */
/*   qy   data value vector.              */
/*   qz   returning answer.               */
/*  qrnp  define value of points number.  */
/*  qrmp  define value of param number.   */
/******************************************/
void Rqranalysis(int m,int n,real qa[][RQRNP],
                real qy[],real qz[])
{
 int i,j,k;
 static real a,qq[RQRMP][RQRNP],rq[RQRMP][RQRMP];

 for( k = 1; k <= m; k++)
    {
    a = 0;
      for(i = 1; i <=n; i++)
        a = a + qa[k][i] * qa[k][i];
      rq[k][k] = sqrt(a);
      for( i = 1; i <= n; i++)
        qq[k][i] = qa[k][i] / rq[k][k];
      for( j = k+1; j <= m; j++)
        {
        a = 0;
          for( i = 1; i <= n; i++)
            a = a + qq[k][i] * qa[j][i];
        rq[j][k] = a;
          for( i = 1; i <= n; i++)
            qa[j][i] = qa[j][i] - qq[k][i] * rq[j][k];
        }
    }

  for( i = 1; i <= m; i++)
    {
    a = 0;
      for( j = 1; j <= n; j++)
        a = a + qq[i][j] * qy[j];
    qz[i] = a;
    }
  qz[m] = qz[m] / rq[m][m];
  for( i = m - 1; i >= 1; i--)
    {
      for( j = i + 1; j <= m; j++)
        qz[i] = qz[i] - qz[j] * rq[j][i];
    qz[i] = qz[i] / rq[i][i];
    }
}

/******** mathod for general matrix Rmat44 *********/
/*** initialise matrix ***/
void Rmat44::IRmat44()
{
	initmat(mat);
}

/*** set [I] matrix ***/
void Rmat44::initmat(Rmat amat)
{
	int i, j;
	for (i = 0; i < 4; i++){
		for (j = 0; j < 4; j++){
			amat[i][j] = 0.0;
			if (i == j)
				amat[i][j] = 1.0;
		}
	}
}

/*** rotate around Y axis ***/
void Rmat44::rotp(real p)
{
	Rmat tmp;
	
	initmat(tmp);
	tmp[0][0] = cos(p);
	tmp[2][0] = sin(p);
	tmp[0][2] = - sin(p);
	tmp[2][2] = cos(p);
	matmul(tmp);
}

/*** rotate around X axis ***/
void Rmat44::rotq(real p)
{
	Rmat tmp;
	
	initmat(tmp);
	tmp[1][1] = cos(p);
	tmp[2][1] = -sin(p);
	tmp[1][2] = sin(p);
	tmp[2][2] = cos(p);
	matmul(tmp);
}

/*** rotate around Z axis ***/
void Rmat44::rotr(real p)
{
	Rmat tmp;
	
	initmat(tmp);
	tmp[0][0] = cos(p);
	tmp[1][0] = -sin(p);
	tmp[0][1] = sin(p);
	tmp[1][1] = cos(p);
	matmul(tmp);
}

/*** shift coord ***/
void Rmat44::move(real x, real y, real z)
{
	Rmat tmp;
	
	initmat(tmp);
	tmp[0][3] = x;
	tmp[1][3] = y;
	tmp[2][3] = z;
	matmul(tmp);
}

/*** set to perspective matrix ***/
void Rmat44::bepers()
{
	Rmat tmp;
	
	initmat(tmp);
		/* project plane is always z=1 and eye point is x=y=z=0 */
	tmp[3][2] = 1.0;
	matmul(tmp);
}

/*** magnify value ***/
void Rmat44::mul(real x, real y, real z)
{
	Rmat tmp;
	
	initmat(tmp);
	tmp[0][0] = x;
	tmp[1][1] = y;
	tmp[2][2] = z;
	matmul(tmp);
}

/*** multiple two matrices ***/
void Rmat44::matmul(Rmat mat2)
{
	int i, j, k;
	Rmat tmp;
	
	for (i = 0; i < 4; i++){
		for (j = 0; j < 4; j++){
			tmp[i][j] = mat[i][j];
		}
	}
	for (i = 0; i < 4; i++){
		for (j = 0; j < 4; j++){
			mat[i][j] = 0.0;
			for (k = 0; k < 4; k++){
				mat[i][j] += tmp[k][j] * mat2[i][k];
			}
		}
	}
}

/*** copy matrix ***/
void Rmat44::copyfrom(Rmat mat2)
{
	int i,j;
	for (i = 0; i < 4; i++){
		for (j = 0; j < 4; j++){
			mat[i][j] = mat2[i][j];
		}
	}
}

/*** get 2D coordinates from 3D coordinates from the matrix ***/
void Rmat44::getcoord2(Rcoord *coord, Rcoord *cd, int pers)
{
	int i,j;
	real w;
	
	coord->x = mat[0][0] * cd->x + mat[0][1] * cd->y +
				mat[0][2] * cd->z + mat[0][3];
	coord->y = mat[1][0] * cd->x + mat[1][1] * cd->y +
				mat[1][2] * cd->z + mat[1][3];
	coord->z = 0.0;
	if (pers){
		w	 = mat[3][0] * cd->x + mat[3][1] * cd->y +
				mat[3][2] * cd->z + mat[3][3];
		w = 1.0 / w;
		coord->x *= w;
		coord->y *= w;
	}
}

/************* method for Roku matrix Rmatobj **************/
/*** initialize matrices ***/
void Rmatobj::IRmatobj()
{
	devmat = new(Rmat44);
	symm = new(Rmat44);
	devmat->IRmat44();
	symm->IRmat44();
						/*** You can change these values in your application ***/
	pp = -0.9;				/* default angle */
	pq = -0.25;
	pr = 0.0;
	
	dx = 0.0;				/* default position */		
	dy = 0.0;
	dz = 1.5;
	
	bx = by = 100.0;		/* default magnify value */
	
	pers = 1;				/* default is perspective (not orthogonal) */
}

/*** make matrix of [3D coordinates -> 2D screen coordinates] ***/
void Rmatobj::makedev()
{
	devmat->IRmat44();
	devmat->rotp(pp);			/* rotate around Y axis */
	devmat->rotq(pq);			/* rotate around X axis */
	devmat->rotr(pr);			/* rotate around Z axis */
	devmat->move(dx, dy, dz);	/* shift position */
	if (pers)
		devmat->bepers();		/* set to perspective matrix */
	dsmat->mul(bx, by, 1.0); 	/* magnify to be screen scale */
}

/*** get 2D screen coordinates from 3D coordinates ***/
Rcoord Rmatobj::getdev(Rcoord *pt)
{
	Rcoord tmp;

	devmat->getcoord2	(&tmp, pt, pers);
	return(tmp);
}

/*** set plane symmetry matrix ***/
void Rmatobj::setsymmatx()
{
	symm->IRmat44();
	symm->mat[0][0] = - 1.0;
	symm->matmul(devmat->mat);
}

/*** calc 3D coordinates from two 2D coordinates ***/
Rcoord Rmatobj::calc3d(Rcoord *pt1, Rcoord *pt2)
{
	Rcoord p;
	real qa[RQRMP][RQRNP], qy[RQRNP], qz[RQRMP];
	
	qa[1][1] = devmat->mat[0][0] - devmat->mat[3][0] * pt1->x;
	qa[2][1] = devmat->mat[0][1] - devmat->mat[3][1] * pt1->x;
	qa[3][1] = devmat->mat[0][2] - devmat->mat[3][2] * pt1->x;
	qy[1] = - devmat->mat[0][3] + devmat->mat[3][3] * pt1->x;

	qa[1][2] = devmat->mat[1][0] - devmat->mat[3][0] * pt1->y;
	qa[2][2] = devmat->mat[1][1] - devmat->mat[3][1] * pt1->y;
	qa[3][2] = devmat->mat[1][2] - devmat->mat[3][2] * pt1->y;
	qy[2] = - devmat->mat[1][3] + devmat->mat[3][3] * pt1->y;

	qa[1][3] = symm->mat[0][0] - symm->mat[3][0] * pt2->x;
	qa[2][3] = symm->mat[0][1] - symm->mat[3][1] * pt2->x;
	qa[3][3] = symm->mat[0][2] - symm->mat[3][2] * pt2->x;
	qy[3] = - symm->mat[0][3] + symm->mat[3][3] * pt2->x;

	qa[1][4] = symm->mat[1][0] - symm->mat[3][0] * pt2->y;
	qa[2][4] = symm->mat[1][1] - symm->mat[3][1] * pt2->y;
	qa[3][4] = symm->mat[1][2] - symm->mat[3][2] * pt2->y;
	qy[4] = - symm->mat[1][3] + symm->mat[3][3] * pt2->y;
	
	Rqranalysis(3, 4, qa, qy, qz);
	p.x = qz[1];
	p.y = qz[2];
	p.z = qz[3];
	return(p);
}

/************** summary of 2D-3D conversion routine in HexaGreat-3D *****************/
void main()
{
	Rmatobj *rmat;
	
	Rcoord gridscpt[MAX], gridpt[MAX];			/* Grid data */
	int ngpoint, ngline, gridline[MAX][2];
	
	Rcoord inputscpt[MAX][2], inputpt[MAX][2];	/* Drawing data */
	int ninputpt, ninputln, inputline[MAX][2];
	
	short x, y;
	int i;
	
			:
	// set grid coordinates and line informations.
			:
	rmat = new(Rmatobj);
	rmat->IRmatobj();					/* initialize matrix */
			:
	// changing display parameters (such as rmat->pp)
			:
	rmat->makedev();					/* get conversion matrix */

	for (i = 0; i < ngpoint; i++){		/* calculate screen coordinates of grid */
		gridscpt[i] = rmat->getdev(&gridpt[i]);
	}
	for (i = 0; i < ngline; i++){		/* display grid to screen */
		x = (short)(gridscpt[ gridline[i][0] ].x);
		y = (short)(gridscpt[ gridline[i][0] ].y);
		MoveTo(x, y);
		x = (short)(gridscpt[ gridline[i][1] ].x);
		y = (short)(gridscpt[ gridline[i][1] ].y);
		LineTo(x, y);
	}
			:
	// input 2D drawings.
	// at this moment, inputscpt[i][0] contains 2D screen coordinate of 
	// left side point, and inputscpt[i][1] contains 2D screen coordinate 
	// of right side point.
	// (only x and y values are used. z value is not used.
			:
	rmat->setsymmatx();					/* set symmetric matrix */
	for (i = 0; i < ninputpt; i++){
										/* calculate 3D coordinate */
		inputpt[i][0] = rmat->calc3d(&inputscpt[i][0], &inputscpt[i][1]);
										/* invert X axis sign */
		inputpt[i][1] = inputpt[i][0];
		inputpt[i][1].x = - inputpt[i][1].x;
	}
			:
	// consequently, inputpt[i][0] contains 3D coordinate of left side point, 
	// and inputpt[i][1] contains 3D coordinate of right side point.
	
}



