/* polymx.c: blanking polygon mask generation for zgrid. */

/* This is inspired by the "inside" function that is part of
Roger Lukas CONTOUR package; but it is coded entirely afresh.
It is heavily specialized for use with zgrid.  It would not
be too hard to modify it for more general applications.

The function is written on the assumption that the x and y grid
vectors are monotonically increasing.  zblank.m is a wrapper that
enforces this requirement before passing the arguments on to
polymx.

Eric Firing
97/05/12   (working; no cleanup or refinement yet)
*/

#include "mex.h"
#include "matrix.h"


#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* Defines that substitute mex versions of C library functions. */
#define calloc mxCalloc
#define realloc mxRealloc
#define printf mexPrintf
#define free   mxFree


/* The following are set so that multiplying by minus 1 switches
   from in to out and vice versa
*/
#define IS_IN     1
#define IS_ON     0
#define IS_OUT    -1

/* Similarly, the code depends on DIR_HOR being 0 here so it
   can be added to another direction without changing it.
*/
#define DIR_UP    1
#define DIR_DOWN -1
#define DIR_HOR   0

typedef struct
{
   int n;
   int *right;
   double *x;
} REGION_TYPE;

typedef struct
{
   int n;
   int *direction;
   double *x, *y;
} POLYGON_TYPE;

typedef struct
{
   int n;
   double *v;
} VECTOR_TYPE;

typedef struct
{
   int m,n;
   double *z;
} XY_ARRAY_TYPE;

typedef struct
{
   int is_vertex;    /* 1 if on vertex, 0 if interpolated */
   int side;         /* Index into polygon of first vertex */
   double x;
} CROSSING_TYPE;

typedef struct
{
   int n;
   int size;
   double y;  /* redundant, but maybe convenient */
   CROSSING_TYPE *cross;
}   CROSSING_SET_TYPE;


void set_mask(int i, XY_ARRAY_TYPE *mask, VECTOR_TYPE *x,
               REGION_TYPE *r)
{
   double *z;
   int ii, ir;

   z = mask->z + i * mask->m;

   ii = x->n -1;
   for (ir = r->n-1; ir >= 0; ir--)
   {
      while ( ii >= 0 && x->v[ii] > r->x[ir])
      {
         if (r->right[ir] == IS_OUT) z[ii] = 0.0;
         else                        z[ii] = 1.0;
         ii--;
      }
      if ( ii >= 0 && x->v[ii] == r->x[ir])
      {
         z[ii] = 1.0; /* on the boundary = in */
         ii--;
      }
   }

   while (ii >=0)
   {
      z[ii] = 0.0;
      ii--;
   }
}

void sort_crossings(CROSSING_SET_TYPE *c)
{
   CROSSING_TYPE cr_temp;
   int i, j;

   for (i = 1; i < c->n; i++)
   {
      cr_temp = c->cross[i];
      for (j = i-1; j >= 0; j--)
      {
         if (c->cross[j].x > cr_temp.x)
         {
            c->cross[j+1] = c->cross[j];
         }
         else break;
      }
      c->cross[j+1] = cr_temp;
   }
}

void find_regions(CROSSING_SET_TYPE *c, POLYGON_TYPE *p, REGION_TYPE *r)
{
   int i, ir;
   int prev_direction, present_direction, iv, iv_prev;
   int start_direction, end_direction;
   int on_edge = 0;
   int prev_region = IS_OUT;

   r->n = 0;

   for (i = 0; i < c->n; i++)
   {
      iv = c->cross[i].side;
      iv_prev = iv - 1;
      if (iv_prev < 0) iv_prev = p->n -2; /* 2 because of duplicated vertex */
      prev_direction = p->direction[iv_prev];
      present_direction = p->direction[iv];

      /* A crossing is a region boundary if it is not a vertex in
         the middle of a polygon side lying on the y line.
      */
      if (prev_direction != DIR_HOR || present_direction != DIR_HOR)
      {
         ir = r->n;
         r->n++;
         r->x[ir] = c->cross[i].x;
         if (!c->cross[i].is_vertex)
         {
            r->right[ir] = - prev_region;
            prev_region = - prev_region;
         }
         else if (prev_direction != DIR_HOR && present_direction != DIR_HOR)
         {
            if (prev_direction == present_direction)
               prev_region = - prev_region;
            r->right[ir] = prev_region;
         }
         else if (!on_edge)
         {
            start_direction = prev_direction + present_direction;
            on_edge = 1;
            r->right[ir] = IS_ON;
         }
         else
         {
            end_direction = prev_direction + present_direction;
            on_edge = 0;
            if (start_direction == end_direction)
               prev_region = - prev_region;
            r->right[ir] =  prev_region;
         }
      } /* creation of new region boundary */
   } /* loop through the crossings */
   if (r->right[ir] != IS_OUT)
      printf("Error in find_regions.");
}




void add_crossing(double xv[], double yv[], double y, int iv,
                     CROSSING_SET_TYPE *c)
{
   int i;
   CROSSING_TYPE *cr;

   if (c->n == c->size)
   {
      c->cross = realloc(c->cross, (c->size+10)*sizeof(CROSSING_TYPE));
      if (c->cross == NULL)
         mexErrMsgTxt("Allocation failure.");
      c->size += 10;
   }

   i = c->n;
   c->n++;
   cr = &(c->cross[i]); /* pointer to the new entry */
   cr->side = iv;       /* vertex number */
   if (y == yv[0])
   {
      cr->is_vertex = 1;
      cr->x = xv[0];
   }
   else if (y == yv[1])
   {
      cr->is_vertex = 1;
      cr->x = xv[1];
   }
   else
   {
      cr->is_vertex = 0;
      cr->x = xv[0] + (y - yv[0])*(xv[1] - xv[0])/(yv[1] - yv[0]);
   }

}


void find_crossings( POLYGON_TYPE *p, VECTOR_TYPE *x,
                     VECTOR_TYPE *y, CROSSING_SET_TYPE c[])
{
   int iv, iv1;            /* vertex */
   int iy, iy0, iy1;       /* for y */
   int i0, i1;
   double xv[2], yv[2];
   iy = y->n;
   /* Loop through the vertices, finding all crossings of each vertex pair. */
   /* The last vertex is a duplicate; not in the loop. */
   for (iv = 0; iv < p->n-1; iv++)
   {
      iv1 = iv+1;
      if (p->y[iv] == p->y[iv1])
         p->direction[iv] = DIR_HOR;
      else if (p->y[iv] > p->y[iv1])
         p->direction[iv] = DIR_DOWN;
      else
         p->direction[iv] = DIR_UP;
      if (p->direction[iv] == DIR_DOWN)
      {
         i0 = iv1;
         i1 = iv;
      }
      else
      {
         i0 = iv;
         i1 = iv1;
      }
      xv[0] = p->x[i0];
      yv[0] = p->y[i0];
      xv[1] = p->x[i1];
      yv[1] = p->y[i1];

      iy0 = iy1 = -1; /* Flag value */
      /* Find the first array y that is >= yv[0]. */
      while (iy > 0 && y->v[iy-1] >= yv[0])
         iy--;
      while (iy < y->n && y->v[iy] < yv[0])
         iy++;
      if (p->direction[iv] == -1 && y->v[iy] == yv[0])
         iy++; /* eliminate duplication */
      if (iy < y->n && y->v[iy] >= yv[0] && y->v[iy] <= yv[1])
      {
         iy0 = iy;
         while (iy < y->n && y->v[iy] <= yv[1])
            iy++;
         /* Overshot: back up. */
         iy--;
         if (p->direction[iv] >= 0 && y->v[iy] == yv[1])
            iy--; /* eliminate duplication */
         if (iy > iy0 && y->v[iy] <= yv[1] )
         {
            iy1 = iy;
         }
      }

      if (iy0 > -1)
      {
         if (iy1 == -1) iy1 = iy0;
         for (iy = iy0; iy <= iy1; iy++)
         {
            add_crossing(xv, yv, y->v[iy], iv, &(c[iy]));
         }
      }
   }

}
/* Inputs: grid: xg, yg; polygon: xp, yp. */

void mexFunction(
   int nlhs, mxArray *plhs[],
   int nrhs, const mxArray* prhs[])
{
   double tmp;
   POLYGON_TYPE p;
   VECTOR_TYPE x, y;

   CROSSING_SET_TYPE *c;
   REGION_TYPE r;
   XY_ARRAY_TYPE mask;

   int i, j, ii, icr, inout;
   double yy;


   if (nlhs != 1)
      mexErrMsgTxt("Need 1 output: mask array.");
   if (nrhs != 4)
      mexErrMsgTxt("Need 4 inputs to poly mex.");

   x.n = mxGetM(prhs[0]) * mxGetN(prhs[0]);
   x.v = mxGetPr(prhs[0]);
   y.n = mxGetM(prhs[1]) * mxGetN(prhs[1]);
   y.v = mxGetPr(prhs[1]);

   p.n = mxGetM(prhs[2]) * mxGetN(prhs[2]);
   tmp = mxGetM(prhs[3]) * mxGetN(prhs[3]);
   if (tmp != p.n)
      mexErrMsgTxt("xp and yp lengths must match");
   p.x = mxGetPr(prhs[2]);
   p.y = mxGetPr(prhs[3]);

   p.direction = (int *)calloc(p.n, sizeof(int));

   mask.m = x.n;
   mask.n = y.n;
   mask.z = (double *)calloc(mask.m*mask.n, sizeof(double));

   /* room for array of CROSSING_SET_TYPE, one for each y. */
   c = calloc(y.n, sizeof(CROSSING_SET_TYPE));
   for (i=0; i<y.n; i++)
   {
      c[i].n = 0;
      c[i].size = 0;
      c[i].cross = NULL;
      c[i].y = y.v[i];  /* Maybe eliminate; redundant. */
   }

   /* Allocate max possible space for REGION_TYPE. */
   r.right = calloc(p.n, sizeof(int));
   r.x     = calloc(p.n, sizeof(double));
   find_crossings(&p, &x, &y, c);

   for (i=0; i<y.n; i++)
   {
      int offset = i*x.n;
      if (c[i].n == 0)
      {
         for (ii = 0; ii < x.n; ii++)
            mask.z[offset+ii] = 0.0;
      }
      else if (c[i].n == 1)
      {
         for (ii = 0; ii < x.n; ii++)
            if (x.v[ii] == c[i].cross[0].x)
               mask.z[offset+ii] = 1.0;
            else
               mask.z[offset+ii] = 0.0;
      }
      else
      {
         sort_crossings(&(c[i]));
         find_regions(&(c[i]), &p, &r);
         set_mask(i, &mask, &x, &r);
      }
   }

   plhs[0] = mxCreateDoubleMatrix(x.n, y.n, mxREAL);
   mxSetPr(plhs[0], mask.z);

}  /* mexFunction */



#if 0

/* Test program, for use prior to conversion to mex file: */
void main(int argc, char **argv)
{
   POLYGON_TYPE p;
   VECTOR_TYPE x, y;
   CROSSING_SET_TYPE *c;
   REGION_TYPE r;
   XY_ARRAY_TYPE mask;

   int i, j, ii, icr, inout;
   double yy;

   p.n = 5;
   p.x = (double *)calloc(p.n, sizeof(double));
   p.y = (double *)calloc(p.n, sizeof(double));
   p.direction = (int *)calloc(p.n, sizeof(int));

   p.x[0] = 0.; p.x[1] = 1.; p.x[2] = 0.75; p.x[3] = 0.25;
         p.x[4] = 0.;
   p.y[0] = 0.; p.y[1] = 0.; p.y[2] = 0.5; p.y[3] = 1.;
         p.y[4] = 0.;

   y.n = 20;
   y.v = (double *)calloc(y.n, sizeof(double));
   yy = -0.2;
   for (i = 0; i < y.n; i++)
   {
      y.v[i] = yy;
      yy += 0.1;
   }

   x.n = 22;
   x.v = (double *)calloc(x.n, sizeof(double));
   yy = -0.2;
   for (i = 0; i < x.n; i++)
   {
      x.v[i] = yy;
      yy += 0.09;
   }

   mask.m = x.n;
   mask.n = y.n;
   mask.z = (double *)calloc(mask.m*mask.n, sizeof(double));

   /* room for array of CROSSING_SET_TYPE, one for each y. */
   c = calloc(y.n, sizeof(CROSSING_SET_TYPE));
   for (i=0; i<y.n; i++)
   {
      c[i].n = 0;
      c[i].size = 0;
      c[i].cross = NULL;
      c[i].y = y.v[i];  /* Maybe eliminate; redundant. */
   }

   /* Allocate max possible space for REGION_TYPE. */
   r.right = calloc(p.n, sizeof(int));
   r.x     = calloc(p.n, sizeof(double));

   printf("\n\n");
   find_crossings(&p, &x, &y, c);
   for (i=0; i<y.n; i++)
   {
      int offset = i*x.n;
      if (c[i].n == 0)
      {
         for (ii = 0; ii < x.n; ii++)
            mask.z[offset+ii] = 0.0;
      }
      else if (c[i].n == 1)
      {
         for (ii = 0; ii < x.n; ii++)
            if (x.v[ii] == c[i].cross[0].x)
               mask.z[offset+ii] = 1.0;
            else
               mask.z[offset+ii] = 0.0;
      }
      else
      {
         sort_crossings(&(c[i]));
         find_regions(&(c[i]), &p, &r);
         set_mask(i, &mask, &x, &r);
      }

   }

   for (j = 0; j<mask.n; j++)
   {
      int offset = j*mask.m;

      for (i = 0; i<mask.m; i++)
      {
         printf(" %1d", (int)mask.z[offset+i]);
      }
      printf("\n");
   }
}

#endif