/*=========================================================================
 * Contains PostgreSQL extention
 *
 * PostgreSQL (poly@>point) operator is currently broken beyond a certain 
 * precision. This module provides a contains(poly,point) to use as a 
 * replacement.
 *
 * It addition it adds the currently missing box to point operators, and
 * conversly the point to box operators as well.
 *
 * (c) 2009 Paul Matthews. LGPL ?? .. (to be advised)
 *=======================================================================*/
#include "postgres.h"
#include "utils/geo_decls.h"
#include "fmgr.h"

#ifdef PG_MODULE_MAGIC
PG_MODULE_MAGIC;
#endif

/*-------------------------------------------------------------------------
 * Box op Point routines
 *
 * strictly left of           1    <<    box_point_leftof
 * does not extend right of   2    &<    box_point_notright
 * overlaps                   3    &&    box_point_contains
 * does not extend left of    4    &>    box_point_notleft
 * strictly right of          5    >>    box_point_rightof
 * contains                   7    @>    box_point_contains
 * does not extend above      9    &<|   box_point_notabove
 * strictly below             10   <<|   box_point_below
 * strictly above             11   |>>   box_point_above
 * does not extend below      12   |&>   box_point_notbelow
 *-----------------------------------------------------------------------*/
PG_FUNCTION_INFO_V1(box_point_leftof);
PG_FUNCTION_INFO_V1(box_point_notright);
PG_FUNCTION_INFO_V1(box_point_notleft);
PG_FUNCTION_INFO_V1(box_point_rightof);
PG_FUNCTION_INFO_V1(box_point_contains);
PG_FUNCTION_INFO_V1(box_point_notabove);
PG_FUNCTION_INFO_V1(box_point_below);
PG_FUNCTION_INFO_V1(box_point_above);
PG_FUNCTION_INFO_V1(box_point_notbelow);

/*
 * Box strictly left of point. box << point
 */
Datum box_point_leftof(PG_FUNCTION_ARGS)
{
  BOX   *box   = PG_GETARG_BOX_P(0);
  Point *point = PG_GETARG_POINT_P(1);

  PG_RETURN_BOOL( box->high.x < point->x );
}

/*
 * Box does not extend right of point. box &< point
 */
Datum box_point_notright(PG_FUNCTION_ARGS)
{
  BOX   *box   = PG_GETARG_BOX_P(0);
  Point *point = PG_GETARG_POINT_P(1);

  PG_RETURN_BOOL( box->high.x <= point->x );
}

/*
 * Box does not extend left of point. box &> point
 */
Datum box_point_notleft(PG_FUNCTION_ARGS)
{
  BOX   *box   = PG_GETARG_BOX_P(0);
  Point *point = PG_GETARG_POINT_P(1);

  PG_RETURN_BOOL( box->low.x >= point->x );
}

/*
 * Box strictly right of point. box >> point
 */
Datum box_point_rightof(PG_FUNCTION_ARGS)
{
  BOX   *box   = PG_GETARG_BOX_P(0);
  Point *point = PG_GETARG_POINT_P(1);

  PG_RETURN_BOOL( box->low.x > point->x );
}

/*
 * Box contains point. box @> point.
 */
Datum box_point_contains(PG_FUNCTION_ARGS)
{
  BOX   *box   = PG_GETARG_BOX_P(0);
  Point *point = PG_GETARG_POINT_P(1);
  int    isin  = point->x >= box->low.x  &&
                 point->x <= box->high.x &&
                 point->y >= box->low.y  &&
                 point->y <= box->high.y;
  PG_RETURN_BOOL(isin);
}

/*
 * Box does not extend above point. box &<| point
 */
Datum box_point_notabove(PG_FUNCTION_ARGS)
{
  BOX   *box   = PG_GETARG_BOX_P(0);
  Point *point = PG_GETARG_POINT_P(1);

  PG_RETURN_BOOL( box->high.y <= point->y );
}

/*
 * Box strictly below point. box <<| point
 */
Datum box_point_below(PG_FUNCTION_ARGS)
{
  BOX   *box   = PG_GETARG_BOX_P(0);
  Point *point = PG_GETARG_POINT_P(1);

  PG_RETURN_BOOL( box->high.y < point->y );
}

/*
 * Box strictly above point. box |>> point
 */
Datum box_point_above(PG_FUNCTION_ARGS)
{
  BOX   *box   = PG_GETARG_BOX_P(0);
  Point *point = PG_GETARG_POINT_P(1);

  PG_RETURN_BOOL( box->low.y > point->y );
}

/* 
 * Box does not extend below point. box |&> point
 */
Datum box_point_notbelow(PG_FUNCTION_ARGS)
{
  BOX   *box   = PG_GETARG_BOX_P(0);
  Point *point = PG_GETARG_POINT_P(1);

  PG_RETURN_BOOL( box->low.y >= point->y );
}

/*-------------------------------------------------------------------------
 * Point op Box routines
 *
 * strictly left of           1    <<    point_box_leftof
 * does not extend right of   2    &<    point_box_notright
 * overlaps                   3    &&    point_box_contains
 * does not extend left of    4    &>    point_box_notleft
 * strictly right of          5    >>    point_box_rightof
 * contained by               8    <@    point_box_contained
 * does not extend above      9    &<|   point_box_notabove
 * strictly below             10   <<|   point_box_below
 * strictly above             11   |>>   point_box_above
 * does not extend below      12   |&>   point_box_notbelow
 *
 *-----------------------------------------------------------------------*/
PG_FUNCTION_INFO_V1(point_box_leftof);
PG_FUNCTION_INFO_V1(point_box_notright);
PG_FUNCTION_INFO_V1(point_box_overlaps);
PG_FUNCTION_INFO_V1(point_box_notleft);
PG_FUNCTION_INFO_V1(point_box_rightof);
PG_FUNCTION_INFO_V1(point_box_contained);
PG_FUNCTION_INFO_V1(point_box_notabove);
PG_FUNCTION_INFO_V1(point_box_below);
PG_FUNCTION_INFO_V1(point_box_above);
PG_FUNCTION_INFO_V1(point_box_notbelow);

/*
 * Point left of box. point << box.
 */
Datum point_box_leftof(PG_FUNCTION_ARGS)
{
  Point *point = PG_GETARG_POINT_P(0);
  BOX   *box   = PG_GETARG_BOX_P(1);
  PG_RETURN_BOOL( point->x < box->low.x );
}

/*
 * Point does not extend right of box. point &< box.
 */
Datum point_box_notright(PG_FUNCTION_ARGS)
{
  Point *point = PG_GETARG_POINT_P(0);
  BOX   *box   = PG_GETARG_BOX_P(1);
  PG_RETURN_BOOL( point->x <= box->high.x );
}

/*
 * Point overlaps with box. point && box.
 */
Datum point_box_overlaps(PG_FUNCTION_ARGS)
{
  Point *point = PG_GETARG_POINT_P(0);
  BOX   *box   = PG_GETARG_BOX_P(1);
  int    isin  = point->x >= box->low.x  &&
                 point->x <= box->high.x &&
                 point->y >= box->low.y  &&
                 point->y <= box->high.y;
  PG_RETURN_BOOL(isin);
}

/*
 * Point does not extend left of box. point &> box.
 */
Datum point_box_notleft(PG_FUNCTION_ARGS)
{
  Point *point = PG_GETARG_POINT_P(0);
  BOX   *box   = PG_GETARG_BOX_P(1);
  PG_RETURN_BOOL( point->x >= box->low.x );
}

/*
 * Point right of box. point >> box.
 */
Datum point_box_rightof(PG_FUNCTION_ARGS)
{
  Point *point = PG_GETARG_POINT_P(0);
  BOX   *box   = PG_GETARG_BOX_P(1);
  PG_RETURN_BOOL( point->x > box->high.x );
}

/*
 * Point containd by box. point <@ box.
 */
Datum point_box_contained(PG_FUNCTION_ARGS)
{
  Point *point = PG_GETARG_POINT_P(0);
  BOX   *box   = PG_GETARG_BOX_P(1);
  int    isin  = point->x >= box->low.x  &&
                 point->x <= box->high.x &&
                 point->y >= box->low.y  &&
                 point->y <= box->high.y;
  PG_RETURN_BOOL(isin);
}

/*
 * Point does not extend above box. point &<| box
 */
Datum point_box_notabove(PG_FUNCTION_ARGS)
{
  Point *point = PG_GETARG_POINT_P(0);
  BOX   *box   = PG_GETARG_BOX_P(1);
  PG_RETURN_BOOL( point->y <= box->high.y );
}

/*
 * Point strictly below box. point <<| box
 */
Datum point_box_below(PG_FUNCTION_ARGS)
{
  Point *point = PG_GETARG_POINT_P(0);
  BOX   *box   = PG_GETARG_BOX_P(1);
  PG_RETURN_BOOL( point->y < box->low.y );
}

/*
 * Point strictly above box. point |>> box
 */
Datum point_box_above(PG_FUNCTION_ARGS)
{
  Point *point = PG_GETARG_POINT_P(0);
  BOX   *box   = PG_GETARG_BOX_P(1);
  PG_RETURN_BOOL( point->y > box->high.y );
}

/* 
 * Point does not extend below point. point |&> box
 */
Datum point_box_notbelow(PG_FUNCTION_ARGS)
{
  Point *point = PG_GETARG_POINT_P(0);
  BOX   *box   = PG_GETARG_BOX_P(1);
  PG_RETURN_BOOL( point->y >= box->low.y );
}

/*-------------------------------------------------------------------------
 * contains_polygon_point
 *
 * Returns TRUE if the given point is contained within the given polygon.
 *
 * Below quoted from the comp.graphics.algorithms FAQ.
 *
 * The essence of the ray-crossing method is as follows.
 * Think of standing inside a field with a fence representing the polygon.
 * Then walk north. If you have to jump the fence you know you are now
 * outside the poly. If you have to cross again you know you are now
 * inside again; i.e., if you were inside the field to start with, the
 * total number of fence jumps you would make will be odd, whereas if you
 * were ouside the jumps will be even.
 *
 * The code below is from Wm. Randolph Franklin <wrf@ecse.rpi.edu>
 * (see URL below) with some minor modifications for speed.  It returns
 * 1 for strictly interior points, 0 for strictly exterior, and 0 or 1
 * for points on the boundary.  The boundary behavior is complex but
 * determined; in particular, for a partition of a region into polygons,
 * each point is "in" exactly one polygon. (See p.243 of [O'Rourke (C)]
 * for a discussion of boundary behavior.)
 *
 * A longer explantion of the algorithm can be found at :
 * http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
 *-----------------------------------------------------------------------*/
PG_FUNCTION_INFO_V1(contains_polygon_point);
PG_FUNCTION_INFO_V1(contained_point_polygon);

Datum contains_polygon_point(PG_FUNCTION_ARGS)
{
  POLYGON* polygon;
  Point*   point;
  int      isin;

  polygon = PG_GETARG_POLYGON_P(0);
  point   = PG_GETARG_POINT_P(1);
  isin    = contains_polypt( polygon->npts, polygon->p, point );

  PG_RETURN_BOOL(isin);
}

Datum contained_point_polygon(PG_FUNCTION_ARGS)
{
  POLYGON* polygon;
  Point*   point;
  int      isin;

  point   = PG_GETARG_POINT_P(0);
  polygon = PG_GETARG_POLYGON_P(1);
  isin    = contains_polypt( polygon->npts, polygon->p, point );

  PG_RETURN_BOOL(isin);
}

/*
 * Point in Polygon?
 */
int contains_polypt( int nvert, Point* vertex, Point* test )
{
  int i, j, c = 0;
  for( i=0, j=nvert-1; i<nvert; j=i++ ) {
    if( ((vertex[i].y>test->y) != (vertex[j].y>test->y)) &&
         (test->x < (vertex[j].x-vertex[i].x) * (test->y-vertex[i].y) /
         (vertex[j].y-vertex[i].y) + vertex[i].x) )
      c = !c;
  }
  return c;
}