[postgis-devel] PostGIS release

[Michael Fuhr]

On Thu, Jul 12, 2007 at 06:37:44PM -0600, Michael Fuhr wrote:
> On Thu, Jul 12, 2007 at 06:49:02PM +0100, Mark Cave-Ayland wrote:
> > As for the regress_ogc results, Mike posted them on the list here:
> > http://postgis.refractions.net/pipermail/postgis-devel/2007-July/002621.html.
> 
> I posted additional test cases here:
> 
> http://postgis.refractions.net/pipermail/postgis-devel/2007-July/002633.html

I think I've tracked down what's happening.  If I define PGIS_DEBUG
I get the following:

findLineSegments 0x83f7cb8: not contained

This comes from lines 236ff in lwgeom/lwgeom_rtree.c:

        if(!isContained(root->interval, value))
        {
#ifdef PGIS_DEBUG
                lwnotice("findLineSegments %p: not contained.", root);
#endif
                return NULL;
        }

Later in the same file we find isContained():

uint32 isContained(INTERVAL *interval, double value)
{
        return FP_CONTAINS_INCL(interval->min, value, interval->max) ? 1 : 0;
}

Here's a test program using the relevant macros from lwgeom/liblwgeom.h:

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

#define PGIS_EPSILON 1e-12
#define FP_LTEQ(A, B) ((A - PGIS_EPSILON) < B)
#define FP_CONTAINS_INCL(A, X, B) (FP_LTEQ(A, X) && FP_LTEQ(X, B))

int
main(int argc, char *argv[])
{
    double  a, x, b;

    if (argc != 4) {
        fprintf(stderr, "Syntax: %s a x b\n", argv[0]);
        return 1;
    }

    a = atof(argv[1]);
    x = atof(argv[2]);
    b = atof(argv[3]);

    printf("FP_CONTAINS_INCL(%f, %f, %f) = %d\n", a, x, b, FP_CONTAINS_INCL(a, x, b));

    return 0;
}

Here are a couple of test runs:

./pgis-test 16384 16384 16385
FP_CONTAINS_INCL(16384.000000, 16384.000000, 16385.000000) = 1

./pgis-test 16384.0001 16384.0001 16385
FP_CONTAINS_INCL(16384.000100, 16384.000100, 16385.000000) = 0

When A gets large enough, A - PGIS_EPSILON requires more precision
than available to be distinguishable from A, so A - PGIS_EPSILON
equals A and is no longer less than B.

-- 
Michael Fuhr