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

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;
}

int main( int argc, char ** argv )

{
    SHPHandle	hSHP;
    int			nShapeType, nEntities, i, iPart;
    const char	*pszPlus;
    double		adfMinBound[4], adfMaxBound[4];
    double *xpr,*ypr;
    unsigned long l;
	// 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];
	// Polygon test stuff
	double		latitude,longitude;
	int			name_field=1;
	unsigned long x=0;
    
/* Display a usage message. */
	printf("argc = %d\n",argc);
    if( argc != 5 )
    {
		printf( "%s shp_file name_field latitude longitude\n",argv[0] );
		exit( 1 );
    }

/* Set the latitude and longitude to test */
	latitude = (double) atof(argv[3]);
	longitude = atof(argv[4]);
	name_field = atoi(argv[2]);

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

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

#if 0
/* 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 );

#if 0
    printf( "Shapefile Type: %s   # of Shapes: %d\n\n",
            SHPTypeName( nShapeType ), nEntities );
#endif

		// Check to makesure we are dealing with a polygon shapefile
		if ( nShapeType%5 == 0 )
				{
				printf("Polygon File\n");
				} else {
				printf("Not a polygon File!\n");
				return 1;
				}


    
    printf( "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] );

printf("(Point to test) Latitude: %04.3f Longitude: %04.3f\n",latitude,longitude);

		if ( (longitude < adfMinBound[1]) || (longitude > adfMaxBound[1])
			|| (latitude < adfMinBound[0]) || (latitude > adfMaxBound[0]) )
			{
			fprintf(stdout,"Test point is outside of file bounds.\n");
			exit(10);
			}
    
/* -------------------------------------------------------------------- */
/*		Skim over the list of shapes, printing all the vertices.		*/
/* -------------------------------------------------------------------- */
    for( i = 0; i < nEntities; i++ )
    {
		int				j;
        SHPObject		*psShape;

		psShape = SHPReadObject( hSHP, i );
#if 0
		printf( "\nShape:%d (%s)  nVertices=%d, nParts=%d\n"
                "  Bounds:(%12.3f,%12.3f, %lg, %lg)\n"
                "      to (%12.3f,%12.3f, %lg, %lg)\n",
				i, SHPTypeName(psShape->nSHPType),
                psShape->nVertices, psShape->nParts,
                psShape->dfXMin, psShape->dfYMin,
                psShape->dfZMin, psShape->dfMMin,
                psShape->dfXMax, psShape->dfYMax,
                psShape->dfZMax, psShape->dfMMax );
#endif



		xpr = psShape->padfX;
		ypr = psShape->padfY;

for( x =0 ; x < 15000 ; x++ )
{
		/* check to see if we are in the bounds of this particular object */
		if ( ! (longitude < psShape->dfXMin) || (longitude > psShape->dfXMax)
			|| (latitude < psShape->dfYMin) || (latitude > psShape->dfYMax) )
			{
	
if ( pnpoly(psShape->nVertices, xpr, ypr, latitude, longitude) )
		{
		fprintf(stdout,"Inside polygon:\n ");
		/* print the name of the object from DBF file */
//		fprintf(stderr,"%s\n",DBFReadStringAttribute( hDBF, i, name_field ));
		} else {
		//fprintf(stdout,"Outside polygon... :-(\n");
		}
			}
} // for
}

    SHPClose( hSHP );

#ifdef USE_DBMALLOC
    malloc_dump(2);
#endif

    exit( 0 );
}