/*
 * TriMesh: out of a domain outline file ( Argus format ), 
 * use the Triangle package to create a triangular mesh 
 *
 */

#include "./TriMesh.h"


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


	/*Matlab arrays: */
	mxArray* pmxa_array=NULL;
	int i,j;
	int counter,counter2,backcounter;
	int prhs_counter;
	
	/* returned quantities: */

	double* index=NULL;
	double* x=NULL;
	double* y=NULL;
	double* segments=NULL;
	double*    segmentmarkerlist=NULL;

	/* input: */
	char*  domainname=NULL;
	char*  riftname=NULL;
	double area;
	char*  order=NULL;
	
	/*Domain outline variables: */
	int      nprof;
	int*     profngrids=NULL;
	double** pprofx=NULL;
	double** pprofy=NULL;
	double*  xprof=NULL;
	double*  yprof=NULL;
	int      numberofpoints;

	/*Rift outline variables: */
	int      numrifts;
	int*     riftsnumgrids=NULL;
	double** riftsgridsx=NULL;
	double** riftsgridsy=NULL;

	/* error handling: */
	int		noerr=1;	
  
	/* Triangle structures: */
	struct triangulateio in,out;
	char   options[256];

	/* verify correct usage: */
	if (nlhs==0 && nrhs==0) {
		/* special case: */
		TriMeshUsage();
		return;
	}

	if (!(  (nlhs==5) ||(nrhs==2) || (nrhs==3)  || (nrhs==4) )){
		mexPrintf("   %s format error.\n", __FUNCT__);
		TriMeshUsage();
		printf("   ");
		mexErrMsgTxt(" ");
	}

	/*Fetch data needed by Triangle: */

	prhs_counter=0;
	/*First recover the domain outline file name: */
	if (!mxIsChar(prhs[prhs_counter])){
		mexPrintf("%s%s\n",__FUNCT__," error message; first argument should be the domain outline file name!");
		mexErrMsgTxt(" ");
	}
	domainname = (char *) mxMalloc((mxGetN(prhs[prhs_counter])+1)*sizeof(char));
	mxGetString(prhs[prhs_counter],domainname,mxGetN(prhs[prhs_counter])+1);

	/*Look for optional rifts file name: */
	prhs_counter++;
	if (mxIsChar(prhs[prhs_counter])){
		riftname = (char *) mxMalloc((mxGetN(prhs[prhs_counter])+1)*sizeof(char));
		mxGetString(prhs[prhs_counter],riftname,mxGetN(prhs[prhs_counter])+1);
		prhs_counter++;
	}

	/*Recover the mesh density desired:*/
	area=mxGetScalar(prhs[prhs_counter]);

	/*Optionaly, recover desired order: */
	prhs_counter++;
	if (mxIsChar(prhs[prhs_counter])){
		order = (char *) mxMalloc((mxGetN(prhs[prhs_counter])+1)*sizeof(char));
		mxGetString(prhs[prhs_counter],order,mxGetN(prhs[prhs_counter])+1);
	}
	
	#ifdef _ISSM_DEBUG_
	printf("Domain name: %s\n",domainname);
	printf("Rift name: %s\n",riftname);
	printf("Maximum area: %16.16lf\n",area);
	printf("Order: %s\n",order);
	#endif
	
	
	/*Start reading the domain outline file: */
	noerr=DomainOutlineRead(&nprof,&profngrids,&pprofx,&pprofy,domainname);
	if(!noerr){
		printf("%s%s%s\n",__FUNCT__," error message reading domain outline ",domainname);
		mexErrMsgTxt(" ");
	}

	/*Read rifts file if present: */
	if(riftname){
		noerr=DomainOutlineRead(&numrifts,&riftsnumgrids,&riftsgridsx,&riftsgridsy,riftname);
		if(!noerr){
			printf("%s%s%s\n",__FUNCT__," error message reading rifts outline ",riftname);
			mexErrMsgTxt(" ");
		}
	}

	#ifdef _ISSM_DEBUG_
	printf("# points in domain outline: %i\n",nprof); 
	#endif
	#ifdef _DEBUG2_
	for (i=0;i<nprof;i++){
		printf("Profile #%i\n",i);
		x=pprofx[i];
		y=pprofy[i];
		for (j=0;j<profngrids[i];j++){
			printf("   %lf %lf\n",x[j],y[j]);
		}
	}
	#endif


	#ifdef _ISSM_DEBUG_
	printf("# points in rift outline: %i\n",numrifts); 
	#endif
	#ifdef _DEBUG2_
	for (i=0;i<numrifts;i++){
		printf("Profile #%i\n",i);
		x=riftsgridsx[i];
		y=riftsgridsy[i];
		for (j=0;j<riftsnumgrids[i];j++){
			printf("   %lf %lf\n",x[j],y[j]);
		}
	}
	#endif

	/*Create initial triangulation to call triangulate():*/
	numberofpoints=0;
	for (i=0;i<nprof;i++){
		numberofpoints+=profngrids[i];
	}
	if (riftname){
		for (i=0;i<numrifts;i++){
			numberofpoints+=riftsnumgrids[i];
		}
	}
	in.numberofpoints=numberofpoints;

	in.numberofpointattributes=1;
	in.pointlist = (REAL *) mxMalloc(in.numberofpoints * 2 * sizeof(REAL));

	counter=0;
	for (i=0;i<nprof;i++){
		xprof=pprofx[i];
		yprof=pprofy[i];
		for (j=0;j<profngrids[i];j++){
			in.pointlist[2*counter+0]=xprof[j];
			in.pointlist[2*counter+1]=yprof[j];
			counter++;
		}
	}
	if(riftname){
		for (i=0;i<numrifts;i++){
			xprof=riftsgridsx[i];
			yprof=riftsgridsy[i];
			for (j=0;j<riftsnumgrids[i];j++){
				in.pointlist[2*counter+0]=xprof[j];
				in.pointlist[2*counter+1]=yprof[j];
				counter++;
			}
		}
	}
	
	in.pointattributelist = (REAL *) mxMalloc(in.numberofpoints *
										  in.numberofpointattributes *
										  sizeof(REAL));
	for (i=0;i<in.numberofpoints;i++){
		in.pointattributelist[i] = 0.0;
	}
	in.pointmarkerlist = (int *) mxMalloc(in.numberofpoints * sizeof(int));
	for(i=0;i<in.numberofpoints;i++){
		in.pointmarkerlist[i] = 0;
	}
	

	/*Build segments: */
	/*Figure out number of segments: holes and closed outlines have as many segments as grids, 
	 *for rifts, we have one less segment as we have grids*/
	in.numberofsegments=0;
	for (i=0;i<nprof;i++){
		in.numberofsegments+=profngrids[i];
	}
	if (riftname){
		for (i=0;i<numrifts;i++){
			in.numberofsegments+=riftsnumgrids[i]-1;
		}
	}
	
	in.segmentlist = (int *) mxMalloc(in.numberofsegments * 2 * sizeof(int));
	in.segmentmarkerlist = (int *) mxCalloc(in.numberofsegments,sizeof(int));
	counter=0;
	backcounter=0;
	for (i=0;i<nprof;i++){
		for (j=0;j<(profngrids[i]-1);j++){
			in.segmentlist[2*counter+0]=counter;
			in.segmentlist[2*counter+1]=counter+1;
			in.segmentmarkerlist[counter]=0;
			counter++;
		}
		/*Close this profile: */
		 in.segmentlist[2*counter+0]=counter;
		 in.segmentlist[2*counter+1]=backcounter;
		 in.segmentmarkerlist[counter]=0;
		 counter++;
		 backcounter=counter;
	}
	counter2=counter;
	if(riftname){
		for (i=0;i<numrifts;i++){
			for (j=0;j<(riftsnumgrids[i]-1);j++){
				in.segmentlist[2*counter2+0]=counter;
				in.segmentlist[2*counter2+1]=counter+1;
				in.segmentmarkerlist[counter2]=2+i;
				counter2++;
				counter++;
			}
			counter++;
		}
	}

	
	/*Build regions: */
	in.numberofregions = 0;

	/*Build holes: */
	in.numberofholes = nprof-1; /*everything is a hole, but for the first profile.*/
	in.holelist = (REAL *) mxMalloc(in.numberofholes * 2 * sizeof(REAL));
	for (i=0;i<nprof-1;i++){
		/*We are looking for a grid that lies inside the hole: */
		GridInsideHole(&in.holelist[2*i+0],&in.holelist[2*i+1],profngrids[i+1],pprofx[i+1],pprofy[i+1]);
	}

	/* Make necessary initializations so that Triangle can return a */
	/*   triangulation in `out': */

	out.pointlist = (REAL *) NULL;            
	out.pointattributelist = (REAL *) NULL;
	out.pointmarkerlist = (int *) NULL; 
	out.trianglelist = (int *) NULL;          
	out.triangleattributelist = (REAL *) NULL;
	out.neighborlist = (int *) NULL;         
	out.segmentlist = (int *) NULL;
	out.segmentmarkerlist = (int *) NULL;
	out.edgelist = (int *) NULL;             
	out.edgemarkerlist = (int *) NULL;   

	/* Triangulate the points:.  Switches are chosen to read and write a  */
	/*   PSLG (p), preserve the convex hull (c), number everything from  */
	/*   zero (z), assign a regional attribute to each element (A), and  */
	/*   produce an edge list (e), a Voronoi diagram (v), and a triangle */
	/*   neighbor list (n).                                              */

	sprintf(options,"%s%lf","pQzDq30ia",area); /*replace V by Q to quiet down the logging*/
  
	triangulate(options, &in, &out, NULL);
	/*report(&out, 0, 1, 1, 1, 1, 0);*/

	/*Allocate index, x and y: */
	index=(double*)mxMalloc(3*out.numberoftriangles*sizeof(double));
	x=(double*)mxMalloc(out.numberofpoints*sizeof(double));
	y=(double*)mxMalloc(out.numberofpoints*sizeof(double));
	segments=(double*)mxMalloc(3*out.numberofsegments*sizeof(double));
	segmentmarkerlist=(double*)mxMalloc(out.numberofsegments*sizeof(double));

	for (i = 0; i < out.numberoftriangles; i++) {
		for (j = 0; j < out.numberofcorners; j++) {
			*(index+3*i+j)=(double)out.trianglelist[i * out.numberofcorners + j]+1;
		}
	}
	for (i = 0; i < out.numberofpoints; i++) {
		x[i]=out.pointlist[i * 2 + 0];
		y[i]=out.pointlist[i * 2 + 1];
	}
	
	for (i = 0; i < out.numberofsegments; i++) {
		segments[3*i+0]=(double)out.segmentlist[i*2+0]+1;
		segments[3*i+1]=(double)out.segmentlist[i*2+1]+1;
		segmentmarkerlist[i]=(double)out.segmentmarkerlist[i];
	}

	/*Associate elements with segments: */
	noerr=AssociateSegmentToElement(&segments,out.numberofsegments,index,out.numberoftriangles);

	/*Order segments so that their normals point outside the domain: */
	if(!strcmp(order,"yes")){
		OrderSegments(&segments,out.numberofsegments, index,out.numberoftriangles);
	}

	/*Output : */
	pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
	mxSetM(pmxa_array,3);
	mxSetN(pmxa_array,out.numberoftriangles);
	mxSetPr(pmxa_array,index);
	mexCallMATLAB( 1, &plhs[0], 1, &pmxa_array, "'");

	pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
	mxSetM(pmxa_array,1);
	mxSetN(pmxa_array,out.numberofpoints);
	mxSetPr(pmxa_array,x);
	mexCallMATLAB( 1, &plhs[1], 1, &pmxa_array, "'");

	pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
	mxSetM(pmxa_array,1);
	mxSetN(pmxa_array,out.numberofpoints);
	mxSetPr(pmxa_array,y);
	mexCallMATLAB( 1, &plhs[2], 1, &pmxa_array, "'");

	pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
	mxSetM(pmxa_array,3);
	mxSetN(pmxa_array,out.numberofsegments);
	mxSetPr(pmxa_array,segments);
	mexCallMATLAB( 1, &plhs[3], 1, &pmxa_array, "'");

	pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
	mxSetM(pmxa_array,1);
	mxSetN(pmxa_array,out.numberofsegments);
	mxSetPr(pmxa_array,segmentmarkerlist);
	mexCallMATLAB( 1, &plhs[4], 1, &pmxa_array, "'");
	
	return;
}

void TriMeshUsage(void)
{
	printf("\n");
	printf("   usage: [index,x,y,segments,segmentmarkers]=TriMesh(domainoutlinefilename,riftsoutlinename,area,ordered) \n");
	printf("      where: index,x,y defines a triangulation, segments is an array made \n");
	printf("      of exterior segments to the mesh domain outline, segmentmarkers is an array flagging each segment \n");
	printf("      (if rifts are present, markers >=2 flag them ), outlinefilename an Argus domain outline file.\n");
	printf("      riftsoutlinename is an Argus domain file, defining rifts (ie: open profiles), \n");
	printf("      area is the maximum area desired for any element of the resulting mesh. \n");
	printf("      and ordered is a string ('yes' or 'no') that determines whether segments are output in the \n");
	printf("      order they are made by Triangle (ie none), or ordered counter clockwise around the domain outline.\n");
	printf("      riftsoutlinename and ordered are optional arguments.\n");
	printf("\n");
}