/*
	trimat.c -- Recalculates a topography on a regular grid

	Date:  Pisa,16-APR-1991

Authors: G. Macedonio and M.T. Pareschi

Modified by R. Ambroziak, USGS/OEMG, Reston, VA 703-648-6168

Usage:
	trimat ptsfile trifile nrows ncols xinf yinf stepx stepy

where:
	ptsfile     -- input file of vertices
	trifile     -- input file of triangles  (built by trimake)
	nrows       -- number of rows    in the output matrix
   ncols       -- number of columns in the output matrix
	xinf,yinf   -- coordinates of the lower left corner
	stepx,stepy -- pixel dimensions
*/
#include <stdio.h>

#define MAXPTS  6000
#define MAXTRI 12000
#define MAXCOL  1024

int Jnd[3]={1,2,0};						/* Auxilary vector */
int huge Vt[3][MAXTRI];					/* List of triangle vertices */
int huge Nt[3][MAXTRI];					/* List of neighbor triangles */
float Px[MAXPTS],Py[MAXPTS];			/* Vertices coordinates */
float Pz[MAXPTS];							/* Heights */
float Out[MAXCOL];						/* Output row */
float huge Par[3][MAXTRI];				/* Interp. plane parameters */

int ptfind(double xa,double ya,int tr);/* Finds external triangle */
int interp(int tr,float *al,float *bl,float *cl);/* Compts plane params */
/* Comuputes matrix determinate */
double dmat2(float pv1[],float pv2[],int p1,int p2); 
double dmat3(float pv1[],float pv2[],float pv[3],int p1,int p2,int p3);
double dmat3r(float pv1[],float pv2[],int p1,int p2,int p3);
void exit();

main(int argc,char *argv[])
{
	int npts;							/* Number of verticies				*/
	int ntri;							/* End of triangle list				*/
	int i,j,it,itt;					/* Auxiliary variables				*/
	int    nr,nc;						/* number of rows and columns			*/
	double x,y;							/* Auxiliary variables				*/
	double xinf,yinf;					/* Coord. of left lower corner			*/
	double xmin,ymin,ymax;			/* Coord. of left lower pixel			*/
	double stepx,stepy;				/* Pixel dimensions					*/
	FILE *ptstrm;						/* Stream of input vertices file		*/
	FILE *trstrm;						/* Stream of output triang. file		*/

	if(argc<9)
	{
		printf("trimat:  Imvalid number of arguments\n\n");
		printf("Command Line:\n");
		printf("\n");
		printf(
"trimat [ptsfile] [trifile] [nrows] [ncols] [xinf] [yinf] [stepx] [stepy]\n");
		printf("\n\n");
		printf("where:\n");
		printf("   ptsfile     -- input file of vertices\n");
		printf("   trifile     -- input file of triangles  (built by trimake)\n");
		printf("   nrows       -- number of rows    in the output matrix\n");
		printf("   ncols       -- number of columns in the output matrix\n");
		printf("   xinf,yinf   -- coordinates of the lower left corner\n");
		printf("   stepx,stepy -- pixel dimensions\n");
		exit(1);
	}
	/* Open input file */
	if((ptstrm=fopen(argv[1],"r"))==NULL)
	{
		printf("Can't open input file '%s'\n",argv[1]);
		exit(1);
	}
	if((trstrm=fopen(argv[2],"r"))==NULL)
	{
		printf("Can't open input file '%s'\n",argv[2]);
		exit(1);
	}
	/* Determine output matrix limits */
	if(sscanf(argv[3],"%d",&nr)<1)
	{
		printf("Invalid 'nrows' -- '%s'\n",argv[3]);
		exit(1);
	}
	if(sscanf(argv[4],"%d",&nc)<1)
	{
		printf("Invalid 'ncols' -- '%s'\n",argv[4]);
		exit(1);
	}
	if(sscanf(argv[5],"%lf",&xinf)<1)
	{
		printf("Invalid 'xinf' -- '%s'\n",argv[5]);
		exit(1);
	}
	if(sscanf(argv[6],"%lf",&yinf)<1)
	{
		printf("Invalid 'yinf' -- '%s'\n",argv[6]);
		exit(1);
	}
	if(sscanf(argv[7],"%lf",&stepx)<1)
	{
		printf("Invalid 'stepx' -- '%s'\n",argv[7]);
		exit(1);
	}
	if(sscanf(argv[8],"%lf",&stepy)<1)
	{
		printf("Invalid 'stepy' -- '%s'\n",argv[8]);
		exit(1);
	}
	xmin=xinf;
	ymin=yinf;
	ymax=(nr-1)*stepy+ymin;

	/* Reads points and quotes */
	printf("Reading data ...\n");
	npts=0;
	while(fscanf(ptstrm,"%f%f%f",&Px[npts],&Py[npts],&Pz[npts])!=EOF)
	{
		npts++;
		if(npts>MAXPTS)
		{
			printf("Too many points\n");
			exit(1);
		}
	}
	printf("There are %d vertices.\n",npts);
	/* Read triangle vertices and neighbor triangles */
	ntri=0;
	while(fscanf(trstrm,"%d%d%d%d%d%d",
		&Vt[0][ntri],&Vt[1][ntri],&Vt[2][ntri],
		&Nt[0][ntri],&Nt[1][ntri],&Nt[2][ntri])!=EOF)
	{
		ntri++;
		if(ntri>MAXTRI)
		{
			printf("Too many triangels.\n");
			exit(1);
		}
	}
	printf("There are %d triangles.\n",ntri);
	/* Begin recalcutation */
	printf("Interpolation parameters calculation ...\n");
	for(it=0;it<ntri;it++)
		interp(it,&Par[0][it],&Par[1][it],&Par[2][it]);
	printf("Beginning matrix calculation ...\n");
	it=0;
	x=xmin;
	y=ymax;
	for(i=0;i<nr;i++)
	{
		for(j=0;j<nc;j++)
		{
			if((itt=ptfind(x,y,it))!=-1)
			{
				it=itt;
				Out[j]=x*Par[0][it]+y*Par[1][it]+Par[2][it];
			}
			else
			{
				Out[j]=0.0;
			}
			x+=stepx;
		}
/* an example of output format */
		if(stepx>0.0)
		{
			for(j=0;j<nc;j++)
				printf(" %7.3f",Out[j]);
		}
		else
		{
			for(j=nc-1;j>=0;j--)
				printf(" %7.3f",Out[j]);
		}
		putchar('\n');
		stepx=-stepx;
		x+=stepx;
		y-=stepy;
	}
}

/***********************************************************************
**
**		ptfind() -- finds the triangle containing point xa,ya
**						tr = Starting triangle
**
************************************************************************/

int ptfind(double xa,double ya, int tr)

{
	int i,aux,newtr;

	newtr=tr;
	do
	{
		aux=0;
		for(i=0;i<=2;i++)
		{
			double diff;
			int i1,ind1,ind2;

			i1=Jnd[i];
			ind1=Vt[i][newtr];
			ind2=Vt[i1][newtr];
			diff=(ya-Py[ind1])*(Px[ind2]-Px[ind1])-
					(Py[ind2]-Py[ind1])*(xa-Px[ind1]);
			if(diff<0.0)
			{
				newtr=Nt[i1][newtr];
				if(newtr==-1)
					return(newtr);
				aux=1;
				break;
			}
		}
	}while(aux==1);
	return(newtr);
}

/***********************************************************************
**
**		interp() -- Computes interpolating plane parameters
**
**
************************************************************************/

int interp(int tr,float *a1,float *b1,float *c1)

{
	int p1,p2,p3;
	double det;

	p1=Vt[0][tr];
	p2=Vt[1][tr];
	p3=Vt[2][tr];
	det=dmat3r(Px,Py,p1,p2,p3);
	*a1=dmat3r(Pz,Py,p1,p2,p3)/det;
	*b1=dmat3r(Px,Pz,p1,p2,p3)/det;
	*c1=dmat3(Px,Py,Pz,p1,p2,p3)/det;

	return(0);
}

/***********************************************************************
**
**		dmat2() -- Computes determinants of a 2*2 matrix
**
**
************************************************************************/

double dmat2(float pv1[],float pv2[],int p1,int p2)

{
	double det;

	det=pv1[p1]*pv2[p2]-pv1[p2]*pv2[p1];
	return(det);
}

/***********************************************************************
**
**		dmat3() -- Computes determinants of a 3*3 matrix
**
**
************************************************************************/

double dmat3(float pv1[],float pv2[],float pv3[],int p1,int p2,int p3)

{
	double det;

	det=pv3[p3]*dmat2(pv1,pv2,p1,p2)-pv3[p2]*dmat2(pv1,pv2,p1,p3)+
		 pv3[p1]*dmat2(pv1,pv2,p2,p3);
	return(det);
}

/***********************************************************************
**
**		dmat3r() -- Computes determinate of a 3*3 matrix with a 
**                constatnt column
**
************************************************************************/

double dmat3r(float pv1[],float pv2[],int p1,int p2,int p3)

{
	double det;

	det=dmat2(pv1,pv2,p1,p2)-dmat2(pv1,pv2,p1,p3)+dmat2(pv1,pv2,p2,p3);
	return(det);
}






        