/* \file exp.h
 * \brief: header file for contour (argus type, files in .exp extension) operations
 */

#ifndef _EXP_H_
#define _EXP_H_

#include <cstring>
#include "../../classes/objects/Contour.h"
#include "../../shared/Numerics/recast.h"
#include "../../Container/Container.h"

int IsInPolySerial(double* in,double* xc,double* yc,int numvertices,double* x,double* y,int nods, int edgevalue);
int DomainOutlineWrite(int nprof,int* profnvertices,double** pprofx,double** pprofy,bool* closed,char* domainname);
int pnpoly(int npol, double *xp, double *yp, double x, double y, int edgevalue);

/*IsInPoly {{{*/
template <class doubletype>
int IsInPoly(SeqVec<doubletype>* in,double* xc,double* yc,int numvertices,double* x,double* y,int i0,int i1, int edgevalue){

	int i;
	double x0,y0;
	doubletype value;
	double xmin=xc[0];
	double xmax=xc[0];
	double ymin=yc[0];
	double ymax=yc[0];

	/*Get extrema*/
	for (i=1;i<numvertices;i++){
		if(xc[i]<xmin) xmin=xc[i];
		if(xc[i]>xmax) xmax=xc[i];
		if(yc[i]<ymin) ymin=yc[i];
		if(yc[i]>ymax) ymax=yc[i];
	}

	/*Go through all vertices of the mesh:*/
	for (i=i0;i<i1;i++){

		//Get current value of value[i] -> do not change it if != 0
		in->GetValue(&value,i);
		if (reCast<bool,doubletype>(value)){
			/*this vertex already is inside one of the contours, continue*/
			continue;
		}

		/*pick up vertex (x[i],y[i]) and figure out if located inside contour (xc,yc)*/
		x0=x[i]; y0=y[i];
		if(x0<xmin || x0>xmax || y0<ymin || y0>ymax){
			value=0;
		}
		else{
			value=pnpoly(numvertices,xc,yc,x0,y0,edgevalue);
		}
		in->SetValue(i,value,INS_VAL);
	}
	 return 1;
}/*}}}*/
/*DomainOutlineRead{{{*/
template <class doubletype>
int DomainOutlineRead(int* pnprof,int** pprofnvertices,doubletype*** ppprofx,doubletype*** ppprofy,bool** pclosed,char* domainname){

	/*indexing: */
	int i,counter;

	/*I/O: */
	FILE   *fid = NULL;
	char    chardummy[256];
	double  ddummy;

	/*output: */
	int          nprof;                //number of profiles in the domainname file
	int         *profnvertices = NULL; //array holding the number of vertices for the nprof profiles
	doubletype **pprofx        = NULL; //array of profiles x coordinates
	doubletype **pprofy        = NULL; //array of profiles y coordinates
	bool        *closed        = NULL; //array holding closed flags for the nprof profiles

	/*For each profile: */
	int         n;
	doubletype *x  = NULL;
	doubletype *y  = NULL;
	bool        cl;

	/*open domain outline file for reading: */
	if ((fid=fopen(domainname,"r"))==NULL){
		_error_("could not find domain file " << domainname); 
	}

	/*Do a first pass through the domainname file, to figure out how many profiles we need to read: */
	nprof=1;
	for(;;){
		fscanf(fid,"%255s %255s\n",chardummy,chardummy);
		fscanf(fid,"%255s %255s\n",chardummy,chardummy);
		fscanf(fid,"%255s %255s %255s %255s\n",chardummy,chardummy,chardummy,chardummy);
		fscanf(fid,"%20u %255s\n",&n,chardummy);
		fscanf(fid,"%255s %255s %255s %255s %255s\n",chardummy,chardummy,chardummy,chardummy,chardummy);
		for (i=0;i<n;i++) fscanf(fid,"%20lf %20lf\n",&ddummy,&ddummy);
		/*check whether we are at the end of the file, otherwise, keep reading next profile:*/
		if (feof(fid)) break;
		nprof++;
	}

	/*Allocate and initialize all the profiles: */
	profnvertices = xNew<int>(nprof);
	pprofx        = xNew<doubletype*>(nprof);
	pprofy        = xNew<doubletype*>(nprof);
	for (i=0;i<nprof;i++){
		pprofx[i] = NULL;
		pprofy[i] = NULL;
	}
	closed=xNew<bool>(nprof);

	/*Reaset file pointer to beginning of file: */
	fseek(fid,0,SEEK_SET);

	/*Start reading profiles: */
	for(counter=0;counter<nprof;counter++){

		/*Skip header: */
		fscanf(fid,"%255s %255s\n",chardummy,chardummy);
		fscanf(fid,"%255s %255s\n",chardummy,chardummy);
		fscanf(fid,"%255s %255s %255s %255s\n",chardummy,chardummy,chardummy,chardummy);

		/*Get number of profile vertices: */
		fscanf(fid,"%20u %255s\n",&n,chardummy);

		/*Skip next line: */
		fscanf(fid,"%255s %255s %255s %255s %255s\n",chardummy,chardummy,chardummy,chardummy,chardummy);

		/*Allocate vertices: */
		x=xNew<doubletype>(n);
		y=xNew<doubletype>(n);

		/*Read vertices: */
		for (i=0;i<n;i++){
			fscanf(fid,"%20lf %20lf\n",&x[i],&y[i]);
		}

		/*Now check that we are dealing with open contours: */
		cl=false;
		if((x[0]==x[n-1]) && (y[0]==y[n-1])){
			cl=true;
		}

		/*Assign pointers: */
		profnvertices[counter]=n;
		pprofx[counter]=x;
		pprofy[counter]=y;
		closed[counter]=cl;
	}

	/*close domain outline file: */
	fclose(fid);

	/*Assign output pointers: */
	*pnprof=nprof;
	*pprofnvertices=profnvertices;
	*ppprofx=pprofx;
	*ppprofy=pprofy;
	if(pclosed) *pclosed=closed;
	else         xDelete<bool>(closed);
} /*}}}*/
/*DataSet* DomainOutlineRead(char* domainname){{{*/
template <class doubletype>
DataSet* DomainOutlineRead(char* domainname){

	/*intermediary: */
	int                  nprof;
	int                 *profnvertices = NULL;
	doubletype         **pprofx        = NULL;
	doubletype         **pprofy        = NULL;

	/*output: */
	DataSet *domain = NULL;

	/*If domainname is an empty string, return empty dataset*/
	if (strcmp(domainname,"")==0){
		nprof=0;
	}
	else{
		DomainOutlineRead(&nprof,&profnvertices,&pprofx, &pprofy, NULL,domainname);
	}

	/*now create dataset of contours: */
	domain=new DataSet(0);

	for(int i=0;i<nprof;i++){
		domain->AddObject(new Contour<doubletype>(i,profnvertices[i],pprofx[i],pprofy[i],1));
	}
	return domain;
} /*}}}*/

#endif
