// test_point_in_poly.c count number of edge crossing horizontal or verticle
//                      line from test point
// c==1, is inside  pnpoly uses horizontal line from test point
//                  pvpoly uses vertical line from test point

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "datread.h" // to get data structure
#include "point_in_poly.h"
#define abs(x) ((x)<0.0?(-(x)):(x))
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))

int main()
{
  float testx = 3.0;
  float testy = 2.0;
  float vertx[5] = {1.0, 3.0, 4.0, 2.0, 2.0}; // no repeat last vert
  float verty[5] = {3.0, 4.0, 2.0, 1.0, 2.0};
  int c, pts, pt;
  int ix, iy, i, j, k;
  int cls;
  char symb[] = {' ', '*'};
  char plot[50][60];
  float stestx = 0.0;
  float stesty = 0.0;
  float starx[10] = { 0.00000,  1.12257,  4.75528,  1.81636,  2.93893,
                     -0.00000, -2.93893, -1.81636, -4.75528, -1.12257};
  float stary[10] = { 5.00000,  1.54508,  1.54508, -0.59017, -4.04508,
		     -1.90983, -4.04508, -0.59017,  1.54508,  1.54508};
  dpts *data_points;
  int num_points;
  int *data_polys;
  int num_polys; 
  float size;
  int nvert;
  // float vertx[3];
  // float verty[3];
  float vertz[3];
  float testz;
  float x, y, z;
  double xmin, xmax, ymin, ymax, zmin, zmax, hx, hy, hz;
  char spikec[11][11][11];

  printf("test_point_in_poly running \n");
  c = pnpoly(5, vertx, verty, testx, testy);
  printf("pnpoly c=%d \n", c);
  for(iy=0; iy<25; iy++)
    for(ix=0; ix<25; ix++)
      plot[iy][ix] = symb[pnpoly(5, vertx, verty,
                                 0.2*(double)ix, 0.2*(double)iy)];
  for(iy=24; iy>=0; iy--)
  {
    for(ix=0; ix<25; ix++) printf("%c", plot[iy][ix]);
    printf(" \n");
  }

  c = pvpoly(5, vertx, verty, testx, testy);
  printf("\npvpoly c=%d \n", c);
  for(iy=0; iy<25; iy++)
    for(ix=0; ix<25; ix++)
      plot[iy][ix] = symb[pvpoly(5, vertx, verty,
                                 0.2*(double)ix, 0.2*(double)iy)];
  for(iy=24; iy>=0; iy--)
  {
    for(ix=0; ix<25; ix++) printf("%c", plot[iy][ix]);
    printf(" \n");
  }
  c = pnpoly(10, starx, stary, stestx, stesty);
  printf("pnpoly star c=%d \n", c);
  for(iy=0; iy<50; iy++)
    for(ix=0; ix<60; ix++)
      plot[iy][ix] = symb[pvpoly(10, starx, stary,
                                 -5.0+0.2*(double)ix, -4.0+0.2*(double)iy)];
  for(iy=49; iy>=0; iy--)
  {
    for(ix=0; ix<60; ix++) printf("%c", plot[iy][ix]);
    printf(" \n");
  }

  for(i=0; i<11; i++)
    for(j=0; j<11; j++)
      for(k=0; k<11; k++) spikec[i][j][k] = ' ';

  datread("spike.dat", &data_points, &num_points,
          &data_polys, &num_polys, &size);
  printf("read %s, num_points=%d, num_polys=%d, size=%f\n",
         "spike.dat", num_points, num_polys, size);
  datclean(            &data_points, &num_points,
          &data_polys, &num_polys, &size);
  printf("read %s, num_points=%d, num_polys=%d, size=%f\n",
         "spike.dat", num_points, num_polys, size);
  cls = datclosed(&data_points, &num_points, &data_polys, &num_polys);
  if(cls!=1) printf("polyhedra not closed \n");

  c = 0;
  nvert = 1;
  testx = 0.0;
  testy = 0.0;
  testz = 0.0;

  c = pitri(nvert, vertx, verty, vertz, testx, testy, testz);
  printf("pitri x=%f, y=%f, z=%f, c=%d\n", testx, testy, testz, c);

  // find points inside volume
  // first get limits
  xmin = data_points[0].x;
  ymin = data_points[0].y;
  zmin = data_points[0].z;
  xmax = data_points[0].x;
  ymax = data_points[0].y;
  zmax = data_points[0].z;
  for(i=1; i<num_points; i++)
  {
    xmin = min(xmin, data_points[i].x);
    ymin = min(ymin, data_points[i].y);
    zmin = min(zmin, data_points[i].z);
    xmax = max(xmax, data_points[i].x);
    ymax = max(ymax, data_points[i].y);
    zmax = max(zmax, data_points[i].z);
  }
  hx = (xmax-xmin)/(double)9;
  hy = (ymax-ymin)/(double)9;
  hz = (zmax-zmin)/(double)9;
  printf("xmin=%f, xmax=%f, hx=%f \n", xmin, xmax, hx);
  printf("ymin=%f, xmax=%f, hy=%f \n", ymin, ymax, hy);
  printf("zmin=%f, xmax=%f, hz=%f \n", zmin, zmax, hz);

  for(x=xmin+hx/2.0; x<xmax-hx/2.0; x+=hx)
  {
    for(y=ymin+hy/2.0; y<ymax-hy/2.0; y+=hy)
    {
      for(z=zmin+hz/2.0; z<zmax-hz/2.0; z+=hz)
      {
        testx = x;
        testy = y;
        testz = z;
        c = ptdat(num_points, num_polys, data_points, data_polys,
	          testx, testy, testz);
        printf("ptdat x=%f, y=%f, z=%f, c=%d\n", testx, testy, testz, c);
        if(c%2==1)
	{
          i = (1.0+x)*4.5;
          j = (1.0+y)*4.5;
          k = (1.0+z)*4.5;
          spikec[i][j][k] = '*';
	}
      } // end z
    } // end y
  } // end x

  for(i=0; i<11; i++)
  {
    printf("x=%f \n", xmin+hx/2.0 + (double)i* hx);
    for(j=0; j<11; j++)
    {
      printf("  ");
      for(k=0; k<11; k++) 
      {
        printf("%c", spikec[i][j][k]);
      }
      printf(" \n");
    }
    printf(" \n \n");
  }
  printf("test_point_in_poly finished \n");
  return 0;
} // end test_point_in_poly.c