/*! \file  MeshProfileIntersection.cpp
    \brief: takes a  .exp file (made of several profiles), and figures out its intersection 
	with a mesh.

	usage:
	[segments]=MeshProfileIntersection(index,x,y,filename);

	where:
	input:
		index,x,y is a triangulation
		filename: name of Argus style .exp file containing the segments (can be groups of disconnected segments)
	output:
		segments: array made of x1,y1,x2,y2,element_id lines (x1,y1) and (x2,y2) are segment extremitis for a segment 
		belonging to the elemnt_id element. there are as many lines in segments as there are segments intersecting the 
		mesh.
*/
	
#include "./MeshProfileIntersection.h"

void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {

	int i,j;

	/* required input: */
	//mesh
	double* double_index=NULL;
	int*    index=NULL;
	int     nel;
	double* x=NULL;
	double* y=NULL;
	int     nods;
	int     dummy;
	
	//contours
	mxArray*  matlabstructure=NULL;
	Contour** contours=NULL;
	int       numcontours;
	Contour*  contouri=NULL;

	/* output: */
	double* segments=NULL;
	int     numsegs;

	/*Boot module: */
	MODULEBOOT();

	/*checks on arguments on the matlab side: */
	CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&MeshProfileIntersectionUsage);

	/*First, call expread on filename to build a contour array in the matlab workspace: */
	mexCallMATLAB( 1, &matlabstructure, 1, (mxArray**)&FILENAME, "expread");

	/*Fetch inputs: */
	//index
	FetchData(&double_index,&nel,&dummy,INDEX);
	if(dummy!=3 && dummy!=6)_error_(" element triangulation should be of 3 or 6 column width!");
	index=(int*)xmalloc(nel*3*sizeof(int));
	for(i=0;i<nel;i++){
		for(j=0;j<3;j++){
			*(index+3*i+j)=(int)*(double_index+dummy*i+j)-1; //"C" style indexing
		}
	}
	//x and y
	FetchData(&x,&nods,X);
	FetchData(&y,&dummy,Y);

	//contours
	numcontours=mxGetNumberOfElements(matlabstructure);
	contours=(Contour**)xmalloc(numcontours*sizeof(Contour*));
	for(i=0;i<numcontours;i++){
		//allocate
		contouri=(Contour*)xmalloc(sizeof(Contour));
		//retrieve dimension of this contour.
		contouri->nods=(int)mxGetScalar(mxGetField(matlabstructure,i,"nods"));
		//set pointers.
		contouri->x=mxGetPr(mxGetField(matlabstructure,i,"x"));
		contouri->y=mxGetPr(mxGetField(matlabstructure,i,"y"));
		*(contours+i)=contouri;
	}

	/* Debugging of contours :{{{1*/
	/*for(i=0;i<numcontours;i++){
		printf("\nContour echo: contour number  %i / %i\n",i+1,numcontours);
		contouri=*(contours+i);
		printf("   Number of vertices %i\n",contouri->nods);
		for (j=0;j<contouri->nods;j++){
			printf("   %lf %lf\n",*(contouri->x+j),*(contouri->y+j));
		}
	}*/
	/*}}}*/

	/*Run interpolation routine: */
	MeshProfileIntersectionx(&segments,&numsegs,index,x,y,nel,nods,contours,numcontours);

	/* output: */
	WriteData(SEGMENTS,segments,numsegs,5);

	/*end module: */
	MODULEEND();

}

void MeshProfileIntersectionUsage(void){
	printf("   usage:\n");
	printf("   [segments]=MeshProfileIntersection(index,x,y,filename);\n");
	printf("   where:\n");
	printf("   input:\n");
	printf("        index,x,y is a triangulation\n");
	printf("        filename: name of Argus style .exp file containing the segments (can be groups of disconnected segments)\n");
	printf("   output:\n");
	printf("        segments: array made of x1,y1,x2,y2,element_id lines (x1,y1) and (x2,y2) are segment extremitis for a segment \n");
	printf("        belonging to the elemnt_id element. there are as many lines in segments as there are segments intersecting the \n");
	printf("        mesh.\n");
}
