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

double atof();

double latitude[15000];
double longitude[15000];
int test_size;


void load_test_points(void)
{
char buffer[512];
int i=0;
float x,y;

while ( (i < 15000) && (0 != fgets(buffer,512,stdin)) )
	{
	sscanf(buffer,"%f%f",&x,&y);
//	printf("load: %d %f %f\n",i,y,x);
	longitude[i] = (double) x;
	latitude[i] = (double) y;
	i++;
	}
test_size = i;
printf("test size: %d\n",test_size);
}

//
//  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 != 3 )
    {
		printf( "%s shp_file name_field\n",argv[0] );
		exit( 1 );
    }

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

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

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

#if 0
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);
			}
#endif
    
/* -------------------------------------------------------------------- */
/*		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 < test_size ; x++ )
//for ( x=0 ; x < 100 ; x++ )
{
//		printf("%d %f %f\n",x,latitude[x],longitude[x]);
		/* check to see if we are in the bounds of this particular object */
		if ( latitude[x] != 500 )
		if ( ! (longitude[x] < psShape->dfXMin) || 
			(longitude[x] > psShape->dfXMax) || 
			(latitude[x] < psShape->dfYMin) || 
			(latitude[x] > psShape->dfYMax) )
			{
	
if ( pnpoly(psShape->nVertices, xpr, ypr, latitude[x], longitude[x]) )
		{
		fprintf(stdout,"%0d: ",x);
		/* print the name of the object from DBF file */
		fprintf(stdout,"%s\n",DBFReadStringAttribute( hDBF, i, name_field ));
		latitude[x] = 500;
		} else {
//		fprintf(stdout,"Outside polygon... :-(\n");
		}
			}
} // for
}

    SHPClose( hSHP );

#ifdef USE_DBMALLOC
    malloc_dump(2);
#endif

    exit( 0 );
}