Suspect that attaching large amounts of code is a breach of etiquette. So apologies in advance.

Needed to write a contains(polygon,point) function that gave a level of accuracy greater than the current 'poly@>point' operator. The new function works just fine. While at it, thought it would be nice to add in all the 'box op point' and 'point op box' routines that are not currently defined.

The code in contains.c implements these functions. In contains.sql you can find how these functions are defined, tided to operators, and then placed in an operator class. (Permission granted to copy any code in here you like, it's all trivial)

The perl script test1 tests all the operators. It compares the output of "box op point" with "box op box(point,point)". It is called with one parameter, the name of the database to use. It creates a table called points within that database. It also creates two files in the current directory a.out and b.out. tkdiff is called to show the differences. Running this shows that all the operators work fine.

The perl script test2 tests 1000 points against 1000 boxes. It is called with one parameter, the name of the database to use. It creates two tables called points and boxes within that database. It also creates two files in the current directory a.out and b.out. tkdiff is called to show the differences. Running this shows that the indexing works fine. (Explain says it uses the index)

Problem:

The following code works in production with 20,000,000 points against 10,000 polygons.
SELECT
      W.geocode,
      F.state_code,
      F.area_code,
      F.area_name
FROM
      work     as W,
      boundary as B,
      features as F
WHERE
      B.boundbox @> box(W.geocode,W.geocode)
  AND contains(B.boundary,W.geocode)
  AND B.boundout = 'T'
  AND (B.feature_id) NOT IN (
      SELECT feature_id
        FROM boundary
       WHERE boundbox @> box(W.geocode,W.geocode)
         AND contains(boundary,W.geocode)
         AND boundout = 'F' )
  AND B.feature_id = F.feature_id;
However the following does not work. It returns an empty set. Note that all that has changed is that "box @> box(point,point)" has been changed to "box @> point". From test1 we know that the operators work, and from test2 we know the indexing works. Confused and bewildered at this point.
SELECT
      W.geocode,
      F.state_code,
      F.area_code,
      F.area_name
FROM
      work     as W,
      boundary as B,
      features as F
WHERE
      B.boundbox @> W.geocode
  AND contains(B.boundary,W.geocode)
  AND B.boundout = 'T'
  AND (B.feature_id) NOT IN (
      SELECT feature_id
        FROM boundary
       WHERE boundbox @> W.geocode
         AND contains(boundary,W.geocode)
         AND boundout = 'F' )
  AND B.feature_id = F.feature_id;
/*=========================================================================
 * 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 <w...@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;
}
---------------------------------------------------------------------------
-- Declare Box+Point functions
---------------------------------------------------------------------------
CREATE OR REPLACE FUNCTION leftof(box,point) RETURNS boolean
LANGUAGE C IMMUTABLE STRICT
AS 'contains.so', 'box_point_leftof';

CREATE OR REPLACE FUNCTION notright(box,point) RETURNS boolean
LANGUAGE C IMMUTABLE STRICT
AS 'contains.so', 'box_point_notright';

CREATE OR REPLACE FUNCTION notleft(box,point) RETURNS boolean
LANGUAGE C IMMUTABLE STRICT
AS 'contains.so', 'box_point_notleft';

CREATE OR REPLACE FUNCTION rightof(box,point) RETURNS boolean
LANGUAGE C IMMUTABLE STRICT
AS 'contains.so', 'box_point_rightof';

CREATE OR REPLACE FUNCTION contains(box,point) RETURNS boolean
LANGUAGE C IMMUTABLE STRICT
AS 'contains.so', 'box_point_contains';

CREATE OR REPLACE FUNCTION notabove(box,point) RETURNS boolean
LANGUAGE C IMMUTABLE STRICT
AS 'contains.so', 'box_point_notabove';

CREATE OR REPLACE FUNCTION below(box,point) RETURNS boolean
LANGUAGE C IMMUTABLE STRICT
AS 'contains.so', 'box_point_below';

CREATE OR REPLACE FUNCTION above(box,point) RETURNS boolean
LANGUAGE C IMMUTABLE STRICT
AS 'contains.so', 'box_point_above';

CREATE OR REPLACE FUNCTION notbelow(box,point) RETURNS boolean
LANGUAGE C IMMUTABLE STRICT
AS 'contains.so', 'box_point_notbelow';

---------------------------------------------------------------------------
-- Grant Execute
---------------------------------------------------------------------------
GRANT EXECUTE ON FUNCTION leftof   (box,point) TO PUBLIC;
GRANT EXECUTE ON FUNCTION notright (box,point) TO PUBLIC;
GRANT EXECUTE ON FUNCTION notleft  (box,point) TO PUBLIC;
GRANT EXECUTE ON FUNCTION rightof  (box,point) TO PUBLIC;
GRANT EXECUTE ON FUNCTION contains (box,point) TO PUBLIC;
GRANT EXECUTE ON FUNCTION notabove (box,point) TO PUBLIC;
GRANT EXECUTE ON FUNCTION below    (box,point) TO PUBLIC;
GRANT EXECUTE ON FUNCTION above    (box,point) TO PUBLIC;
GRANT EXECUTE ON FUNCTION notbelow (box,point) TO PUBLIC;

---------------------------------------------------------------------------
-- Create Operators
---------------------------------------------------------------------------
DROP OPERATOR IF EXISTS <<(box,point);
CREATE OPERATOR << (
  LEFTARG    = box,
  RIGHTARG   = point,
  PROCEDURE  = leftof,
  COMMUTATOR = >>,
  RESTRICT   = positionsel,
  JOIN       = positionjoinsel
);

DROP OPERATOR IF EXISTS &<(box,point);
CREATE OPERATOR &< (
  LEFTARG    = box,
  RIGHTARG   = point,
  PROCEDURE  = notright,
  RESTRICT   = positionsel,
  JOIN       = positionjoinsel
);

DROP OPERATOR IF EXISTS &&(box,point);
CREATE OPERATOR && (
  LEFTARG    = box,
  RIGHTARG   = point,
  PROCEDURE  = contains,
  RESTRICT   = areasel,
  JOIN       = areajoinsel
);

DROP OPERATOR IF EXISTS &>(box,point);
CREATE OPERATOR &> (
  LEFTARG    = box,
  RIGHTARG   = point,
  PROCEDURE  = notleft,
  RESTRICT   = positionsel,
  JOIN       = positionjoinsel
);

DROP OPERATOR IF EXISTS >>(box,point);
CREATE OPERATOR >> (
  LEFTARG    = box,
  RIGHTARG   = point,
  PROCEDURE  = rightof,
  COMMUTATOR = <<,
  RESTRICT   = positionsel,
  JOIN       = positionjoinsel
);

DROP OPERATOR IF EXISTS @>(box,point);
CREATE OPERATOR @> (
  LEFTARG    = box,
  RIGHTARG   = point,
  PROCEDURE  = contains,
  COMMUTATOR = <@,
  RESTRICT   = contsel,
  JOIN       = contjoinsel
);

DROP OPERATOR IF EXISTS &<|(box,point);
CREATE OPERATOR &<| (
  LEFTARG    = box,
  RIGHTARG   = point,
  PROCEDURE  = notabove,
  RESTRICT   = positionsel,
  JOIN       = positionjoinsel
);

DROP OPERATOR IF EXISTS <<|(box,point);
CREATE OPERATOR <<| (
  LEFTARG    = box,
  RIGHTARG   = point,
  PROCEDURE  = below,
  COMMUTATOR = |>>,  
  RESTRICT   = positionsel,
  JOIN       = positionjoinsel
);

DROP OPERATOR IF EXISTS |>>(box,point);
CREATE OPERATOR |>> (
  LEFTARG    = box,
  RIGHTARG   = point,
  PROCEDURE  = above,
  COMMUTATOR = <<|,  
  RESTRICT   = positionsel,
  JOIN       = positionjoinsel
);

DROP OPERATOR IF EXISTS |&>(box,point);
CREATE OPERATOR |&> (
  LEFTARG    = box,
  RIGHTARG   = point,
  PROCEDURE  = notbelow,
  RESTRICT   = positionsel,
  JOIN       = positionjoinsel
);



---------------------------------------------------------------------------
-- Declare Point+Box Functions
---------------------------------------------------------------------------
CREATE OR REPLACE FUNCTION leftof(point,box) RETURNS boolean
LANGUAGE C IMMUTABLE STRICT
AS 'contains.so', 'point_box_leftof';

CREATE OR REPLACE FUNCTION notright(point,box) RETURNS boolean
LANGUAGE C IMMUTABLE STRICT
AS 'contains.so', 'point_box_notright';

CREATE OR REPLACE FUNCTION overlaps(point,box) RETURNS boolean
LANGUAGE C IMMUTABLE STRICT
AS 'contains.so', 'point_box_overlaps';

CREATE OR REPLACE FUNCTION notleft(point,box) RETURNS boolean
LANGUAGE C IMMUTABLE STRICT
AS 'contains.so', 'point_box_notleft';

CREATE OR REPLACE FUNCTION rightof(point,box) RETURNS boolean
LANGUAGE C IMMUTABLE STRICT
AS 'contains.so', 'point_box_rightof';

--CREATE OR REPLACE FUNCTION contained(point,box) RETURNS boolean
--LANGUAGE C IMMUTABLE STRICT
--AS 'contains.so', 'point_box_contained';

CREATE OR REPLACE FUNCTION notabove(point,box) RETURNS boolean
LANGUAGE C IMMUTABLE STRICT
AS 'contains.so', 'point_box_notabove';

CREATE OR REPLACE FUNCTION below(point,box) RETURNS boolean
LANGUAGE C IMMUTABLE STRICT
AS 'contains.so', 'point_box_below';

CREATE OR REPLACE FUNCTION above(point,box) RETURNS boolean
LANGUAGE C IMMUTABLE STRICT
AS 'contains.so', 'point_box_above';

CREATE OR REPLACE FUNCTION notbelow(point,box) RETURNS boolean
LANGUAGE C IMMUTABLE STRICT
AS 'contains.so', 'point_box_notbelow';

---------------------------------------------------------------------------
-- Grant Execute
---------------------------------------------------------------------------
GRANT EXECUTE ON FUNCTION leftof   (point,box) TO PUBLIC;
GRANT EXECUTE ON FUNCTION notright (point,box) TO PUBLIC;
GRANT EXECUTE ON FUNCTION overlaps (point,box) TO PUBLIC;
GRANT EXECUTE ON FUNCTION notleft  (point,box) TO PUBLIC;
GRANT EXECUTE ON FUNCTION rightof  (point,box) TO PUBLIC;
--GRANT EXECUTE ON FUNCTION contained(point,box) TO PUBLIC;
GRANT EXECUTE ON FUNCTION notabove (point,box) TO PUBLIC;
GRANT EXECUTE ON FUNCTION below    (point,box) TO PUBLIC;
GRANT EXECUTE ON FUNCTION above    (point,box) TO PUBLIC;
GRANT EXECUTE ON FUNCTION notbelow (point,box) TO PUBLIC;

---------------------------------------------------------------------------
-- Create Point+Box operators
---------------------------------------------------------------------------
DROP OPERATOR IF EXISTS <<(point,box);
CREATE OPERATOR << (
  LEFTARG    = point,
  RIGHTARG   = box,
  PROCEDURE  = leftof,
  COMMUTATOR = >>,  
  RESTRICT   = positionsel,
  JOIN       = positionjoinsel
);

DROP OPERATOR IF EXISTS &<(point,box);
CREATE OPERATOR &< (
  LEFTARG    = point,
  RIGHTARG   = box,
  PROCEDURE  = notright,
  RESTRICT   = positionsel,
  JOIN       = positionjoinsel
);

DROP OPERATOR IF EXISTS &&(point,box);
CREATE OPERATOR && (
  LEFTARG    = point,
  RIGHTARG   = box,
  PROCEDURE  = overlaps,
  RESTRICT   = areasel,
  JOIN       = areajoinsel
);

DROP OPERATOR IF EXISTS &>(point,box);
CREATE OPERATOR &> (
  LEFTARG    = point,
  RIGHTARG   = box,
  PROCEDURE  = notleft,
  RESTRICT   = positionsel,
  JOIN       = positionjoinsel
);

DROP OPERATOR IF EXISTS >>(point,box);
CREATE OPERATOR >> (
  LEFTARG    = point,
  RIGHTARG   = box,
  PROCEDURE  = rightof,
  COMMUTATOR = <<,  
  RESTRICT   = positionsel,
  JOIN       = positionjoinsel
);

--DROP OPERATOR IF EXISTS <@(point,box);
--CREATE OPERATOR <@ (
--  LEFTARG    = point,
--  RIGHTARG   = box,
--  PROCEDURE  = contained,
--  RESTRICT   = contsel,
--  JOIN       = contjoinsel
--);

DROP OPERATOR IF EXISTS &<|(point,box);
CREATE OPERATOR &<| (
  LEFTARG    = point,
  RIGHTARG   = box,
  PROCEDURE  = notabove,
  RESTRICT   = positionsel,
  JOIN       = positionjoinsel
);

DROP OPERATOR IF EXISTS <<|(point,box);
CREATE OPERATOR <<| (
  LEFTARG    = point,
  RIGHTARG   = box,
  PROCEDURE  = below,
  COMMUTATOR = |>>,  
  RESTRICT   = positionsel,
  JOIN       = positionjoinsel
);

DROP OPERATOR IF EXISTS |>>(point,box);
CREATE OPERATOR |>> (
  LEFTARG    = point,
  RIGHTARG   = box,
  PROCEDURE  = above,
  COMMUTATOR = <<|,  
  RESTRICT   = positionsel,
  JOIN       = positionjoinsel
);

DROP OPERATOR IF EXISTS |&>(point,box);
CREATE OPERATOR |&> (
  LEFTARG    = point,
  RIGHTARG   = box,
  PROCEDURE  = notbelow,
  RESTRICT   = positionsel,
  JOIN       = positionjoinsel
);


-----------------------------------------------------------------------------
---- Alter the existing box_ops family
-----------------------------------------------------------------------------
ALTER OPERATOR FAMILY box_ops USING GiST ADD
  OPERATOR 1  <<  (box,point),
  OPERATOR 2  &<  (box,point),
  OPERATOR 3  &&  (box,point),
  OPERATOR 4  &>  (box,point),
  OPERATOR 5  >>  (box,point),  
  OPERATOR 7  @>  (box,point), 
  OPERATOR 9  &<| (box,point), 
  OPERATOR 10 <<| (box,point), 
  OPERATOR 11 |>> (box,point),   
  OPERATOR 12 |&> (box,point),
  OPERATOR 1  <<  (point,box),
  OPERATOR 2  &<  (point,box),
  OPERATOR 3  &&  (point,box),
  OPERATOR 4  &>  (point,box),
  OPERATOR 5  >>  (point,box),  
--OPERATOR 8  <@  (point,box),
  OPERATOR 9  &<| (point,box), 
  OPERATOR 10 <<| (point,box), 
  OPERATOR 11 |>> (point,box),   
  OPERATOR 12 |&> (point,box);   


---------------------------------------------------------------------------
-- Polygon Functions
---------------------------------------------------------------------------

-- Polygon contains Point
CREATE OR REPLACE FUNCTION contains(polygon,point) RETURNS boolean
LANGUAGE C IMMUTABLE STRICT
AS 'contains.so', 'contains_polygon_point';

GRANT EXECUTE ON FUNCTION contains(polygon,point) TO PUBLIC;


-- Point contained in Polygon
CREATE OR REPLACE FUNCTION contained(point,polygon) RETURNS boolean
LANGUAGE C IMMUTABLE STRICT
AS 'contains.so', 'contained_point_polygon';

GRANT EXECUTE ON FUNCTION contained(point,polygon) TO PUBLIC;
#!/usr/bin/perl
use DBI;
use warnings;
use strict;

my $dbname = shift or die;
my ($dbh, $stm);

$dbh = DBI->connect("dbi:Pg:dbname=$dbname", "", "") or die $!;
$dbh->do(<<EOSQL) or die $!;
DROP TABLE IF EXISTS points;
CREATE TABLE points( 
  pointpoint point,
  pointbox   box
);
EOSQL

for( my $i=0; $i<1000; $i++ ) {
  my $x = rand(2.0)-1.0;
  my $y = rand(2.0)-1.0;
  $dbh->do(<<EOSQL) or die $!;
INSERT INTO points VALUES ( 
  Point('$x,$y'),
  Box(Point('$x,$y'),Point('$x,$y'))
);
EOSQL
}

check_box_point('<<');
check_box_point('&<');
check_box_point('&&');
check_box_point('&>');
check_box_point('>>');
check_box_point('@>');
check_box_point('&<|');
check_box_point('<<|');
check_box_point('|>>');
check_box_point('|&>');

check_point_box('<<');
check_point_box('&<');
check_point_box('&&');
check_point_box('&>');
check_point_box('>>');
check_point_box('<@');
check_point_box('&<|');
check_point_box('<<|');
check_point_box('|>>');
check_point_box('|&>');


#--------------------------------------------------------------------------
# Check 'box OP point' against 'box OP box(point,point)'
#--------------------------------------------------------------------------
sub check_box_point {
  my $op = shift;
  print "box $op point\n";
  
  $stm = $dbh->prepare(<<EOSQL) or die $!;
SELECT
  pointpoint[0] as x,
  pointpoint[1] as y
FROM
  points                                   AS p,
  box(point('-0.5,-0.5'),point('0.5,0.5')) AS b
WHERE 
  b $op pointbox
ORDER BY
  x, y;
EOSQL

  $stm->execute() or die $!;
  open FH, ">a.out";
  while( my $ref = $stm->fetchrow_hashref() ) {
    my $x = $ref->{x};
    my $y = $ref->{y};
    print FH "$x,$y\n";
  }
  close FH;

  $stm = $dbh->prepare(<<EOSQL) or die $!;
SELECT
  pointpoint[0] as x,
  pointpoint[1] as y
FROM
  points                                   AS p,
  box(point('-0.5,-0.5'),point('0.5,0.5')) AS b
WHERE
  b $op pointpoint
ORDER BY
  x, y;
EOSQL

  $stm->execute() or die $!;
  open FH, ">b.out";
  while( my $ref = $stm->fetchrow_hashref() ) {
    my $x = $ref->{x};
    my $y = $ref->{y};
    print FH "$x,$y\n";
  }
  close FH;

  system( "tkdiff a.out b.out >/dev/null 2>&1" );
}


#--------------------------------------------------------------------------
# Check 'point OP box' against 'box(point,point) OP box'
#--------------------------------------------------------------------------
sub check_point_box {
  my $op = shift;
  print "point $op box\n";
  
  $stm = $dbh->prepare(<<EOSQL) or die $!;
SELECT
  pointpoint[0] as x,
  pointpoint[1] as y
FROM
  points                                     AS p,
  box(point('0.25,0.25'),point('0.75,0.75')) AS b
WHERE 
  pointbox $op b
ORDER BY
  x, y;
EOSQL

  $stm->execute() or die $!;
  open FH, ">a.out";
  while( my $ref = $stm->fetchrow_hashref() ) {
    my $x = $ref->{x};
    my $y = $ref->{y};
    print FH "$x,$y\n";
  }
  close FH;

  $stm = $dbh->prepare(<<EOSQL) or die $!;
SELECT
  pointpoint[0] as x,
  pointpoint[1] as y
FROM
  points                                     AS p,
  box(point('0.25,0.25'),point('0.75,0.75')) AS b
WHERE
  pointpoint $op b
ORDER BY
  x, y;
EOSQL

  $stm->execute() or die $!;
  open FH, ">b.out";
  while( my $ref = $stm->fetchrow_hashref() ) {
    my $x = $ref->{x};
    my $y = $ref->{y};
    print FH "$x,$y\n";
  }
  close FH;

  system( "tkdiff a.out b.out >/dev/null 2>&1" );
}
#!/usr/bin/perl
use DBI;
use warnings;
use strict;

my $dbname = shift or die;
my ($dbh, $stm);

print STDERR "Creating Tables ...\n";
$dbh = DBI->connect("dbi:Pg:dbname=$dbname", "", "") or die $!;
$dbh->do(<<EOSQL) or die $!;
DROP TABLE IF EXISTS points;
DROP TABLE IF EXISTS boxes;
CREATE TABLE points(
  pointpoint point,
  pointbox   box
);
CREATE TABLE boxes(
  boxid int4,
  bound box
);
EOSQL

print STDERR "Creating Points ...\n";
for( my $i=0; $i<1000; $i++ ) {
  my $x = rand(1.0);
  my $y = rand(1.0);
  $dbh->do(<<EOSQL) or die $!;
INSERT INTO points VALUES ( 
  Point('$x,$y'),
  Box(Point('$x,$y'),Point('$x,$y'))
);
EOSQL
}

print STDERR "Creating Boxes ...\n";
for( my $i=0; $i<1000; $i++ ) {
  my $w = rand(0.90);
  my $h = rand(0.90);
  my $x = rand(0.05)+0.5;
  my $y = rand(0.05)+0.5;
  my $xmin = $x;
  my $xmax = $x+$w;
  my $ymin = $y;
  my $ymax = $y+$h;
  $dbh->do(<<EOSQL) or die $!;
INSERT INTO boxes VALUES ( 
  $i,
  Box(Point('$xmin,$ymin'),Point('$xmax,$ymax'))
);
EOSQL
}

print STDERR "Creating Index ...\n";
$dbh->do(<<EOSQL) or die $!;
CREATE INDEX boxes_bound ON boxes USING GiST (bound);
EOSQL

print STDERR "Analyse Boxes ...\n";
$dbh->do(<<EOSQL) or die $!;
VACUUM boxes;
EOSQL

print STDERR "Analyse Boxes ...\n";
$dbh->do(<<EOSQL) or die $!;
ANALYSE boxes;
EOSQL

print STDERR "Query 1 ...\n";
$stm = $dbh->prepare(<<EOSQL) or die $!;
set enable_seqscan = off;
SELECT
  pointpoint[0] as x,
  pointpoint[1] as y,
  boxid
FROM
  points AS p,
  boxes  AS b
WHERE 
  bound @> pointbox
EOSQL

$stm->execute() or die $!;
open FH, "| sort > a.out";
while( my $ref = $stm->fetchrow_hashref() ) {
  my $x  = $ref->{x};
  my $y  = $ref->{y};
  my $id = $ref->{boxid};
  print FH "$x,$y $id\n";
}
close FH;

print STDERR "Query 2 ...\n";
$stm = $dbh->prepare(<<EOSQL) or die $!;
set enable_seqscan = off;  
SELECT
  pointpoint[0] as x,
  pointpoint[1] as y,
  boxid
FROM
  points AS p,
  boxes  AS b
WHERE 
  bound @> pointpoint
EOSQL

$stm->execute() or die $!;
open FH, "| sort > b.out";
while( my $ref = $stm->fetchrow_hashref() ) {
  my $x = $ref->{x};
  my $y = $ref->{y};
  my $id = $ref->{boxid};
  print FH "$x,$y $id\n";
}
close FH;

system( "tkdiff a.out b.out >/dev/null 2>&1" );
-- 
Sent via pgsql-hackers mailing list (pgsql-hackers@postgresql.org)
To make changes to your subscription:
http://www.postgresql.org/mailpref/pgsql-hackers

Reply via email to