#include "shapefil.h"
#include <stdio.h>
#include <stdlib.h>

struct PT {
	double x;
	double y;
};
typedef struct PT POINT;

double atof();

//
//  Input parameters:
//     npol               : number of vertices
//     xp[npol], yp[npol] : x,y coord of vertices
//     x,y                : coord of point to test 
//
//  Return Value:
//     0 : test point is outside polygon
//     1 : test point is inside polygon
//
//  Notes:
//     if test point is on the border, 0 or 1 is returned.
//     If there exists an adiacent polygon, the point is
//     _in_ only in one of the two.
//

int pnpoly(int npol, double *xp, double *yp, float x, float y) {
	int i, j, c = 0;
	for (i = 0, j = npol-1; i < npol; j = i++) { 
		if (	( ((yp[i]<=y) && (y<yp[j])) || ((yp[j]<=y) && (y<yp[i]))) &&
				(x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i])) {
          c = !c;
		}
	}
return c;
}

/* Returns true if structure is filled with a valid lat / lon pair */
int getTestPoint(POINT *search) {
	char buffer[256];

	if ( (0 != fgets(buffer,256,stdin)) ) {
		//sscanf(buffer,"%lf %lf",&search->x,&search->y);
		sscanf(buffer,"%lf %lf",&search->y,&search->x); // lat / lon
		return 1;
	}
return 0;
}

typedef struct {
	SHPObject	*o;
	int	value;
	int	hitCount;
	int	record;
}	SHPstr;

int SHPstrcmp( SHPstr *a, SHPstr *b)
{
return	b->value - a->value;	
	
}

int main( int argc, char **argv ) {
	SHPHandle hSHP;
	SHPstr *s;
	int nShapeType, nEntities, i, iPart;
	const char *pszPlus;
	double adfMinBound[4], adfMaxBound[4];
	
	// DBF Stuff
	DBFHandle hDBF;
	int *panWidth, iRecord;
	char szFormat[32], szField[1024], *pszFilename = NULL;
	int nWidth, nDecimals;
	int bHeader = 0;
	int bRaw = 0;
	int bMultiLine = 0;
	char szTitle[12];

	// PnPoly stuff
	int totalVerticies=0;
	POINT search;
	int found;
	int j;
	

	
	/* Open the passed shapefile. */
	hSHP = SHPOpen( argv[1], "rb" );

	if( hSHP == NULL ) {
		printf( "Unable to open:%s\n", argv[1] );
		exit( 1 );
	}
#if 1

	/* Open DBF file associated with shapefile. */
	pszFilename = argv[1];	// DBF
	hDBF = DBFOpen( pszFilename, "rb" );
	if( hDBF == NULL ) {
		printf( "DBFOpen(%s,\"r\") failed.\n", argv[1] );
		exit( 2 );
	}
	/* Warn user that DBF file doesn't contain anything */
	if( DBFGetFieldCount(hDBF) == 0 ) {
		printf( "There are no fields in this table!\n" );
		exit( 3 );
	}

#endif
	/* Print out the file bounds. */
	SHPGetInfo( hSHP, &nEntities, &nShapeType, adfMinBound, adfMaxBound );


	// Check to makesure we are dealing with a polygon shapefile
	if ( nShapeType != SHPT_POLYGON && nShapeType != SHPT_POLYGONZ && nShapeType != SHPT_POLYGONM) {
		printf("Not a polygon File!\n");
		return 1;
	}

	fprintf(stderr, "File Bounds: (%12.3f,%12.3f,%lg,%lg)\n"
			"         to  (%12.3f,%12.3f,%lg,%lg)\n", 
				adfMinBound[0], 
            adfMinBound[1], 
            adfMinBound[2], 
            adfMinBound[3], 
            adfMaxBound[0], 
            adfMaxBound[1], 
            adfMaxBound[2], 
            adfMaxBound[3] );
	
	/*	Load an array with pointers to all the shapes.	*/
s = calloc(nEntities,sizeof(SHPstr));
	for ( i = 0; i < nEntities; i++ ) {
		s[i].record = i;
		s[i].o = SHPReadObject( hSHP, i );
		s[i].value = atoi(DBFReadStringAttribute( hDBF, i, 7 ));
		totalVerticies += s[i].o->nVertices;
	} 

qsort(s,nEntities,sizeof(SHPstr),SHPstrcmp);

	for ( i=0 ; i<nEntities ; i++ ) {
		printf("[%d] %d=population %d\n",i,s[i].value,s[i].record);	
	}

	fprintf(stderr,"%d Verticies and %d Polygons loaded into memory\n",totalVerticies,nEntities);

	/* Reading to start checking if input falls within the polygon */
	while ( getTestPoint(&search) ) {
//		printf("Search point: (%lf, %lf)\n",search.x,search.y);

		for ( i=0,found=0 ; ( !found && i<nEntities ) ; i++ ) {
			// printf("Searching shape %d\n",i);

			// BUG - probably only works for north west hemisphere
			if (	(search.x < s[i].o->dfXMax) && (search.x > s[i].o->dfXMin) &&  
			 		(search.y < s[i].o->dfYMax) && (search.y > s[i].o->dfYMin) ) {
				if ( pnpoly(s[i].o->nVertices, s[i].o->padfX, s[i].o->padfY, search.x, search.y) ) {
					/* Each time we get a hit we move the
					 * state up in priority so it gets 
					 * searched sooner */
					SHPstr x;
					printf("Found it in shape %d\n",i);
					/* Mark the hit */
					s[i].hitCount++;
					found=1;
#if 1
					if ( 0 < i ) {
						x=s[i-1];
						s[i-1]=s[i];
						s[i]=x;
					}
#endif

				}
			}
		}
#if 0
		if ( !found ) {
			printf("Shape not found\n");
		}
#endif
		
	}
#if 1
	for ( i=0 ; i < nEntities ; i++ ) {
		//j=atoi(DBFReadStringAttribute( hDBF, searchNext(i), 4 ));
		j=atoi(DBFReadStringAttribute( hDBF,i, 2 ));
		printf("%d %d\n",j,s[i].hitCount);
	}
#endif
	
	fprintf(stderr,"done\n");
	

	free(s);
	SHPClose( hSHP );
	DBFClose( hDBF );

#ifdef USE_DBMALLOC
    malloc_dump(2);
#endif

	exit( 0 );
}