From 03cbe9584514c680b3e226b51c1695e03cd4101d Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Tue, 15 Sep 2009 19:50:24 +0000 Subject: [PATCH] Geocentric bounding box roughed in and compiles. git-svn-id: http://svn.osgeo.org/postgis/trunk@4499 b70326c6-7e19-0410-871a-916f4a2858ee --- liblwgeom/cunit/Makefile.in | 6 +- liblwgeom/cunit/cu_algorithm.h | 6 +- liblwgeom/cunit/cu_geodetic.c | 206 ++++++++++++++++ liblwgeom/cunit/cu_geodetic.h | 29 +++ liblwgeom/cunit/cu_libgeom.c | 145 +---------- liblwgeom/cunit/cu_libgeom.h | 9 +- liblwgeom/cunit/cu_measures.h | 6 +- liblwgeom/cunit/cu_tester.c | 13 +- liblwgeom/cunit/cu_tester.h | 16 ++ liblwgeom/g_box.c | 23 +- liblwgeom/g_serialized.c | 22 +- liblwgeom/libgeom.h | 11 +- liblwgeom/lwgeodetic.c | 435 +++++++++++++++++++++++++++++++-- liblwgeom/lwgeodetic.h | 26 +- postgis/geography_gist.c | 2 +- postgis/geography_inout.c | 4 +- 16 files changed, 749 insertions(+), 210 deletions(-) create mode 100644 liblwgeom/cunit/cu_geodetic.c create mode 100644 liblwgeom/cunit/cu_geodetic.h create mode 100644 liblwgeom/cunit/cu_tester.h diff --git a/liblwgeom/cunit/Makefile.in b/liblwgeom/cunit/Makefile.in index f1c262dad..0a0520191 100644 --- a/liblwgeom/cunit/Makefile.in +++ b/liblwgeom/cunit/Makefile.in @@ -16,10 +16,12 @@ CFLAGS=@CFLAGS@ @WARNFLAGS@ CUNIT_LDFLAGS=@CUNIT_LDFLAGS@ CUNIT_CPPFLAGS=@CUNIT_CPPFLAGS@ -I.. -OBJS= cu_tester.o \ +OBJS= \ cu_algorithm.o \ + cu_geodetic.o \ cu_measures.o \ - cu_libgeom.o + cu_libgeom.o \ + cu_tester.o # If we couldn't find the cunit library then display a helpful message ifeq ($(CUNIT_LDFLAGS),) diff --git a/liblwgeom/cunit/cu_algorithm.h b/liblwgeom/cunit/cu_algorithm.h index 4538212fb..637e121e4 100644 --- a/liblwgeom/cunit/cu_algorithm.h +++ b/liblwgeom/cunit/cu_algorithm.h @@ -16,16 +16,12 @@ #include "CUnit/Basic.h" #include "lwalgorithm.h" +#include "cu_tester.h" /*********************************************************************** ** for Computational Geometry Suite */ -/* Set-up / clean-up functions */ -CU_pSuite register_cg_suite(void); -int init_cg_suite(void); -int clean_cg_suite(void); - /* Test functions */ void test_lw_segment_side(void); void test_lw_segment_intersects(void); diff --git a/liblwgeom/cunit/cu_geodetic.c b/liblwgeom/cunit/cu_geodetic.c new file mode 100644 index 000000000..1e8373061 --- /dev/null +++ b/liblwgeom/cunit/cu_geodetic.c @@ -0,0 +1,206 @@ +/********************************************************************** + * $Id$ + * + * PostGIS - Spatial Types for PostgreSQL + * http://postgis.refractions.net + * Copyright 2009 Paul Ramsey + * + * This is free software; you can redistribute and/or modify it under + * the terms of the GNU General Public Licence. See the COPYING file. + * + **********************************************************************/ + +#include "cu_geodetic.h" + +/* +** Called from test harness to register the tests in this file. +*/ +CU_pSuite register_geodetic_suite(void) +{ + CU_pSuite pSuite; + pSuite = CU_add_suite("Geodetic Suite", init_geodetic_suite, clean_geodetic_suite); + if (NULL == pSuite) + { + CU_cleanup_registry(); + return NULL; + } + + if ( + (NULL == CU_add_test(pSuite, "test_gbox_from_spherical_coordinates()", test_gbox_from_spherical_coordinates)) || + (NULL == CU_add_test(pSuite, "test_gserialized_get_gbox_geocentric()", test_gserialized_get_gbox_geocentric)) || + (NULL == CU_add_test(pSuite, "test_gbox_calculation()", test_gbox_calculation)) + ) + { + CU_cleanup_registry(); + return NULL; + } + return pSuite; +} + +/* +** The suite initialization function. +** Create any re-used objects. +*/ +int init_geodetic_suite(void) +{ + return 0; +} + +/* +** The suite cleanup function. +** Frees any global objects. +*/ +int clean_geodetic_suite(void) +{ + return 0; +} + +void test_gbox_from_spherical_coordinates(void) +{ + double ll[64]; + GBOX *box = gbox_new(gflags(0, 0, 1)); + int rv; + POINTARRAY *pa; + + ll[0] = 0.0; + ll[1] = 0.0; + ll[2] = 45.0; + ll[3] = 45.0; + + pa = pointArray_construct((uchar*)ll, 0, 0, 2); + rv = ptarray_calculate_gbox_geodetic(pa, box); + //printf("%s\n", gbox_to_string(box, gflags(0, 0, 1))); + CU_ASSERT_DOUBLE_EQUAL(box->xmin, 0.707106781187, 0.001); + lwfree(pa); + + ll[0] = 45.0; + ll[1] = 45.0; + ll[2] = 45.5; + ll[3] = 45.5; + ll[4] = 46.0; + ll[5] = 46.0; + + pa = pointArray_construct((uchar*)ll, 0, 0, 3); + rv = ptarray_calculate_gbox_geodetic(pa, box); + //printf("%s\n", gbox_to_string(box, gflags(0, 0, 1))); + CU_ASSERT_DOUBLE_EQUAL(box->xmin, 0.694658370459, 0.001); + lwfree(pa); + + + ll[0] = 45.0; + ll[1] = 45.0; + + pa = pointArray_construct((uchar*)ll, 0, 0, 1); + rv = ptarray_calculate_gbox_geodetic(pa, box); + //printf("%s\n", gbox_to_string(box, gflags(0, 0, 1))); + CU_ASSERT_DOUBLE_EQUAL(box->xmin, 0.707106781187, 0.001); + lwfree(pa); + + lwfree(box); + +} + +void test_gserialized_get_gbox_geocentric(void) +{ + LWGEOM *lwg; + GSERIALIZED *g; + GBOX *gbox; + int i; + + char wkt[][512] = { + "POINT(45 45)", + "MULTILINESTRING((45 45, 45 45),(45 45, 45 45))", + "MULTILINESTRING((45 45, 45 45),(45 45, 45.5 45.5),(45 45, 46 46))", + "GEOMETRYCOLLECTION(MULTILINESTRING((45 45, 45 45),(45 45, 45.5 45.5),(45 45, 46 46)))", + "MULTIPOLYGON(((45 45, 45 45, 45 45, 45 45),(45 45, 45 45, 45 45, 45 45)))", + "LINESTRING(45 45, 45.1 45.1)", + "MULTIPOINT(45 45, 45.05 45.05, 45.1 45.1)" + }; + + double xmin[] = { + 0.707106781187, + 0.707106781187, + 0.694658370459, + 0.694658370459, + 0.707106781187, + 0.705871570679, + 0.705871570679 + }; + + for ( i = 0; i < 7; i++ ) + { + //printf("%s\n", wkt[i]); + lwg = lwgeom_from_ewkt(wkt[i], PARSER_CHECK_NONE); + g = gserialized_from_lwgeom(lwg, 1, 0); + g->flags = FLAGS_SET_GEODETIC(g->flags, 1); + lwgeom_free(lwg); + gbox = gserialized_calculate_gbox_geocentric(g); + //printf("%s\n", gbox_to_string(box, g->flags)); + CU_ASSERT_DOUBLE_EQUAL(gbox->xmin, xmin[i], 0.00000001); + lwfree(g); + lwfree(gbox); + } + +} + +/* +* Build LWGEOM on top of *aligned* structure so we can use the read-only +* point access methods on them. +*/ +static LWGEOM* lwgeom_over_gserialized(char *wkt) +{ + LWGEOM *lwg; + GSERIALIZED *g; + + lwg = lwgeom_from_ewkt(wkt, PARSER_CHECK_NONE); + g = gserialized_from_lwgeom(lwg, 1, 0); + lwgeom_free(lwg); + return lwgeom_from_gserialized(g); +} + +void test_gbox_calculation(void) +{ + + LWGEOM *geom; + int i = 0; + GBOX *gbox = gbox_new(gflags(0,0,0)); + BOX3D *box3d; + + char ewkt[][512] = { + "POINT(0 0.2)", + "LINESTRING(-1 -1,-1 2.5,2 2,2 -1)", + "SRID=1;MULTILINESTRING((-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1))", + "POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0))", + "SRID=4326;MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)),((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)))", + "SRID=4326;GEOMETRYCOLLECTION(POINT(0 1),POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0)),MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5))))", + "POINT(0 220.2)", + "LINESTRING(-1 -1,-1231 2.5,2 2,2 -1)", + "SRID=1;MULTILINESTRING((-1 -131,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1))", + "POLYGON((-1 -1,-1 2.5,2 2,2 -133,-1 -1),(0 0,0 1,1 1,1 0,0 0))", + "SRID=4326;MULTIPOLYGON(((-1 -1,-1 2.5,211 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)),((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)))", + "SRID=4326;GEOMETRYCOLLECTION(POINT(0 1),POLYGON((-1 -1,-1111 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0)),MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5))))", + }; + + for( i = 0; i < 6; i++ ) + { + geom = lwgeom_over_gserialized(ewkt[i]); + lwgeom_calculate_gbox(geom, gbox); + box3d = lwgeom_compute_box3d(geom); + //printf("%g %g\n", gbox->xmin, box3d->xmin); + CU_ASSERT_EQUAL(gbox->xmin, box3d->xmin); + CU_ASSERT_EQUAL(gbox->xmax, box3d->xmax); + CU_ASSERT_EQUAL(gbox->ymin, box3d->ymin); + CU_ASSERT_EQUAL(gbox->ymax, box3d->ymax); + lwfree(box3d); + } + + lwfree(gbox); + + +} + + + + + + diff --git a/liblwgeom/cunit/cu_geodetic.h b/liblwgeom/cunit/cu_geodetic.h new file mode 100644 index 000000000..a8f6266c1 --- /dev/null +++ b/liblwgeom/cunit/cu_geodetic.h @@ -0,0 +1,29 @@ +/********************************************************************** + * $Id$ + * + * PostGIS - Spatial Types for PostgreSQL + * http://postgis.refractions.net + * Copyright 2009 Paul Ramsey + * + * This is free software; you can redistribute and/or modify it under + * the terms of the GNU General Public Licence. See the COPYING file. + * + **********************************************************************/ + +#include +#include +#include +#include "CUnit/Basic.h" + +#include "lwgeodetic.h" +#include "cu_tester.h" + +/*********************************************************************** +** For new geography library lwgeodetic.h +*/ + + +/* Test functions */ +void test_gbox_from_spherical_coordinates(void); +void test_gserialized_get_gbox_geocentric(void); +void test_gbox_calculation(void); diff --git a/liblwgeom/cunit/cu_libgeom.c b/liblwgeom/cunit/cu_libgeom.c index c4d659143..64fcf38d1 100644 --- a/liblwgeom/cunit/cu_libgeom.c +++ b/liblwgeom/cunit/cu_libgeom.c @@ -35,10 +35,7 @@ CU_pSuite register_libgeom_suite(void) (NULL == CU_add_test(pSuite, "test_geometry_type_from_string()", test_geometry_type_from_string)) || (NULL == CU_add_test(pSuite, "test_lwgeom_check_geodetic()", test_lwgeom_check_geodetic)) || (NULL == CU_add_test(pSuite, "test_lwgeom_count_vertices()", test_lwgeom_count_vertices)) || - (NULL == CU_add_test(pSuite, "test_on_gser_lwgeom_count_vertices()", test_on_gser_lwgeom_count_vertices)) || - (NULL == CU_add_test(pSuite, "test_gbox_from_spherical_coordinates()", test_gbox_from_spherical_coordinates)) || - (NULL == CU_add_test(pSuite, "test_gserialized_get_gbox_geocentric()", test_gserialized_get_gbox_geocentric)) || - (NULL == CU_add_test(pSuite, "test_gbox_calculation()", test_gbox_calculation)) + (NULL == CU_add_test(pSuite, "test_on_gser_lwgeom_count_vertices()", test_on_gser_lwgeom_count_vertices)) ) { CU_cleanup_registry(); @@ -345,7 +342,7 @@ void test_lwgeom_check_geodetic(void) for( i = 0; i < 6; i++ ) { geom = lwgeom_over_gserialized(ewkt[i]); - CU_ASSERT_EQUAL(lwgeom_check_geodetic(geom), G_TRUE); + CU_ASSERT_EQUAL(lwgeom_check_geodetic(geom), LW_TRUE); lwgeom_free(geom); } @@ -353,7 +350,7 @@ void test_lwgeom_check_geodetic(void) { //char *out_ewkt; geom = lwgeom_over_gserialized(ewkt[i]); - CU_ASSERT_EQUAL(lwgeom_check_geodetic(geom), G_FALSE); + CU_ASSERT_EQUAL(lwgeom_check_geodetic(geom), LW_FALSE); //out_ewkt = lwgeom_to_ewkt(geom, PARSER_CHECK_NONE); //printf("%s\n", out_ewkt); lwgeom_free(geom); @@ -404,139 +401,3 @@ void test_on_gser_lwgeom_count_vertices(void) lwfree(g_ser1); } - - -void test_gbox_from_spherical_coordinates(void) -{ - double ll[64]; - GBOX *box = gbox_new(gflags(0, 0, 1)); - int rv; - POINTARRAY *pa; - - ll[0] = 0.0; - ll[1] = 0.0; - ll[2] = 45.0; - ll[3] = 45.0; - - pa = pointArray_construct((uchar*)ll, 0, 0, 2); - rv = ptarray_calculate_gbox_geodetic(pa, box); - //printf("%s\n", gbox_to_string(box, gflags(0, 0, 1))); - CU_ASSERT_DOUBLE_EQUAL(box->xmin, 0.707106781187, 0.001); - lwfree(pa); - - ll[0] = 45.0; - ll[1] = 45.0; - ll[2] = 45.5; - ll[3] = 45.5; - ll[4] = 46.0; - ll[5] = 46.0; - - pa = pointArray_construct((uchar*)ll, 0, 0, 3); - rv = ptarray_calculate_gbox_geodetic(pa, box); - //printf("%s\n", gbox_to_string(box, gflags(0, 0, 1))); - CU_ASSERT_DOUBLE_EQUAL(box->xmin, 0.694658370459, 0.001); - lwfree(pa); - - - ll[0] = 45.0; - ll[1] = 45.0; - - pa = pointArray_construct((uchar*)ll, 0, 0, 1); - rv = ptarray_calculate_gbox_geodetic(pa, box); - //printf("%s\n", gbox_to_string(box, gflags(0, 0, 1))); - CU_ASSERT_DOUBLE_EQUAL(box->xmin, 0.707106781187, 0.001); - lwfree(pa); - - lwfree(box); - -} - -void test_gserialized_get_gbox_geocentric(void) -{ - LWGEOM *lwg; - GSERIALIZED *g; - GBOX *gbox; - int i; - - char wkt[][512] = { - "POINT(45 45)", - "MULTILINESTRING((45 45, 45 45),(45 45, 45 45))", - "MULTILINESTRING((45 45, 45 45),(45 45, 45.5 45.5),(45 45, 46 46))", - "GEOMETRYCOLLECTION(MULTILINESTRING((45 45, 45 45),(45 45, 45.5 45.5),(45 45, 46 46)))", - "MULTIPOLYGON(((45 45, 45 45, 45 45, 45 45),(45 45, 45 45, 45 45, 45 45)))", - "LINESTRING(45 45, 45.1 45.1)", - "MULTIPOINT(45 45, 45.05 45.05, 45.1 45.1)" - }; - - double xmin[] = { - 0.707106781187, - 0.707106781187, - 0.694658370459, - 0.694658370459, - 0.707106781187, - 0.705871570679, - 0.705871570679 - }; - - for ( i = 0; i < 7; i++ ) - { - //printf("%s\n", wkt[i]); - lwg = lwgeom_from_ewkt(wkt[i], PARSER_CHECK_NONE); - g = gserialized_from_lwgeom(lwg, 1, 0); - g->flags = FLAGS_SET_GEODETIC(g->flags, 1); - lwgeom_free(lwg); - gbox = gserialized_calculate_gbox_geocentric(g); - //printf("%s\n", gbox_to_string(box, g->flags)); - CU_ASSERT_DOUBLE_EQUAL(gbox->xmin, xmin[i], 0.00000001); - lwfree(g); - lwfree(gbox); - } - -} - -void test_gbox_calculation(void) -{ - - LWGEOM *geom; - int i = 0; - GBOX *gbox = gbox_new(gflags(0,0,0)); - BOX3D *box3d; - - char ewkt[][512] = { - "POINT(0 0.2)", - "LINESTRING(-1 -1,-1 2.5,2 2,2 -1)", - "SRID=1;MULTILINESTRING((-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1))", - "POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0))", - "SRID=4326;MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)),((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)))", - "SRID=4326;GEOMETRYCOLLECTION(POINT(0 1),POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0)),MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5))))", - "POINT(0 220.2)", - "LINESTRING(-1 -1,-1231 2.5,2 2,2 -1)", - "SRID=1;MULTILINESTRING((-1 -131,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1))", - "POLYGON((-1 -1,-1 2.5,2 2,2 -133,-1 -1),(0 0,0 1,1 1,1 0,0 0))", - "SRID=4326;MULTIPOLYGON(((-1 -1,-1 2.5,211 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)),((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)))", - "SRID=4326;GEOMETRYCOLLECTION(POINT(0 1),POLYGON((-1 -1,-1111 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0)),MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5))))", - }; - - for( i = 0; i < 6; i++ ) - { - geom = lwgeom_over_gserialized(ewkt[i]); - lwgeom_calculate_gbox(geom, gbox); - box3d = lwgeom_compute_box3d(geom); - //printf("%g %g\n", gbox->xmin, box3d->xmin); - CU_ASSERT_EQUAL(gbox->xmin, box3d->xmin); - CU_ASSERT_EQUAL(gbox->xmax, box3d->xmax); - CU_ASSERT_EQUAL(gbox->ymin, box3d->ymin); - CU_ASSERT_EQUAL(gbox->ymax, box3d->ymax); - lwfree(box3d); - } - - lwfree(gbox); - - -} - - - - - - diff --git a/liblwgeom/cunit/cu_libgeom.h b/liblwgeom/cunit/cu_libgeom.h index 717b638fb..c8bfc14f9 100644 --- a/liblwgeom/cunit/cu_libgeom.h +++ b/liblwgeom/cunit/cu_libgeom.h @@ -16,15 +16,13 @@ #include "CUnit/Basic.h" #include "libgeom.h" +#include "cu_tester.h" /*********************************************************************** ** For new geometry library libgeom.h */ -/* Set-up / clean-up functions */ -CU_pSuite register_libgeom_suite(void); -int init_libgeom_suite(void); -int clean_libgeom_suite(void); + /* Test functions */ void test_typmod_macros(void); @@ -37,6 +35,3 @@ void test_geometry_type_from_string(void); void test_lwgeom_check_geodetic(void); void test_lwgeom_count_vertices(void); void test_on_gser_lwgeom_count_vertices(void); -void test_gbox_from_spherical_coordinates(void); -void test_gserialized_get_gbox_geocentric(void); -void test_gbox_calculation(void); diff --git a/liblwgeom/cunit/cu_measures.h b/liblwgeom/cunit/cu_measures.h index b0b94c2e8..3590b9a64 100644 --- a/liblwgeom/cunit/cu_measures.h +++ b/liblwgeom/cunit/cu_measures.h @@ -16,15 +16,13 @@ #include "CUnit/Basic.h" #include "liblwgeom.h" +#include "cu_tester.h" /*********************************************************************** ** for Computational Geometry Suite */ -/* Set-up / clean-up functions */ -CU_pSuite register_measures_suite(void); -int init_measures_suite(void); -int clean_measures_suite(void); + /* Test functions */ void test_mindistance2d_recursive_tolerance(void); diff --git a/liblwgeom/cunit/cu_tester.c b/liblwgeom/cunit/cu_tester.c index 15faf549d..0d1a428c6 100644 --- a/liblwgeom/cunit/cu_tester.c +++ b/liblwgeom/cunit/cu_tester.c @@ -13,10 +13,8 @@ #include #include #include "CUnit/Basic.h" - -#include "cu_algorithm.h" -#include "cu_measures.h" -#include "cu_libgeom.h" +#include "liblwgeom.h" +#include "cu_tester.h" /* @@ -62,6 +60,13 @@ int main() return CU_get_error(); } + /* Add the geodetic suite to the registry */ + if (NULL == register_geodetic_suite()) + { + CU_cleanup_registry(); + return CU_get_error(); + } + /* Run all tests using the CUnit Basic interface */ CU_basic_set_mode(CU_BRM_VERBOSE); CU_basic_run_tests(); diff --git a/liblwgeom/cunit/cu_tester.h b/liblwgeom/cunit/cu_tester.h new file mode 100644 index 000000000..60def6fb4 --- /dev/null +++ b/liblwgeom/cunit/cu_tester.h @@ -0,0 +1,16 @@ + +CU_pSuite register_measures_suite(void); +CU_pSuite register_geodetic_suite(void); +CU_pSuite register_libgeom_suite(void); +CU_pSuite register_cg_suite(void); + +int init_measures_suite(void); +int init_geodetic_suite(void); +int init_libgeom_suite(void); +int init_cg_suite(void); + +int clean_measures_suite(void); +int clean_geodetic_suite(void); +int clean_libgeom_suite(void); +int clean_cg_suite(void); + diff --git a/liblwgeom/g_box.c b/liblwgeom/g_box.c index 3f093bcd4..64a1c143b 100644 --- a/liblwgeom/g_box.c +++ b/liblwgeom/g_box.c @@ -19,6 +19,17 @@ GBOX* gbox_new(uchar flags) return g; } +int gbox_merge_point3d(GBOX *gbox, POINT3D *p) +{ + if( gbox->xmin > p->x ) gbox->xmin = p->x; + if( gbox->ymin > p->y ) gbox->xmin = p->y; + if( gbox->zmin > p->z ) gbox->xmin = p->z; + if( gbox->xmax < p->x ) gbox->xmax = p->x; + if( gbox->ymax < p->y ) gbox->xmax = p->y; + if( gbox->zmax < p->z ) gbox->xmax = p->z; + return G_SUCCESS; +} + int gbox_merge(GBOX *master_box, GBOX *new_box) { assert(master_box); @@ -53,18 +64,18 @@ int gbox_overlaps(GBOX *g1, GBOX *g2) if( g1->xmax < g2->xmin || g1->ymax < g2->ymin || g1->xmin > g2->xmax || g1->ymin > g2->ymax ) - return G_FALSE; + return LW_FALSE; if( FLAGS_GET_Z(g1->flags) || FLAGS_GET_GEODETIC(g1->flags) ) { if( g1->zmax < g2->zmin || g1->zmin > g2->zmax ) - return G_FALSE; + return LW_FALSE; } if( FLAGS_GET_M(g1->flags) ) { if( g1->mmax < g2->mmin || g1->mmin > g2->mmax ) - return G_FALSE; + return LW_FALSE; } - return G_TRUE; + return LW_TRUE; } char* gbox_to_string(GBOX *gbox) @@ -429,7 +440,7 @@ static int lwcollection_calculate_gbox(LWCOLLECTION *coll, GBOX *gbox) GBOX subbox; int i; int result = G_FAILURE; - int first = G_TRUE; + int first = LW_TRUE; assert(coll); if( coll->ngeoms == 0 ) return G_FAILURE; @@ -447,7 +458,7 @@ static int lwcollection_calculate_gbox(LWCOLLECTION *coll, GBOX *gbox) if( first ) { gbox_duplicate(gbox, &subbox); - first = G_FALSE; + first = LW_FALSE; } else { diff --git a/liblwgeom/g_serialized.c b/liblwgeom/g_serialized.c index 756814bd1..cc33d67d4 100644 --- a/liblwgeom/g_serialized.c +++ b/liblwgeom/g_serialized.c @@ -747,31 +747,31 @@ static LWCIRCSTRING* lwcircstring_from_gserialized_buffer(uchar *data_ptr, uchar static int lwcollection_from_gserialized_allowed_types(int collectiontype, int subtype) { if( collectiontype == COLLECTIONTYPE ) - return G_TRUE; + return LW_TRUE; if( collectiontype == MULTIPOINTTYPE && subtype == POINTTYPE ) - return G_TRUE; + return LW_TRUE; if( collectiontype == MULTILINETYPE && subtype == LINETYPE ) - return G_TRUE; + return LW_TRUE; if( collectiontype == MULTIPOLYGONTYPE && subtype == POLYGONTYPE ) - return G_TRUE; + return LW_TRUE; if( collectiontype == COMPOUNDTYPE && (subtype == LINETYPE || subtype == CIRCSTRINGTYPE) ) - return G_TRUE; + return LW_TRUE; if( collectiontype == CURVEPOLYTYPE && (subtype == CIRCSTRINGTYPE || subtype == LINETYPE || subtype == COMPOUNDTYPE) ) - return G_TRUE; + return LW_TRUE; if( collectiontype == MULTICURVETYPE && (subtype == CIRCSTRINGTYPE || subtype == LINETYPE || subtype == COMPOUNDTYPE) ) - return G_TRUE; + return LW_TRUE; if( collectiontype == MULTISURFACETYPE && (subtype == POLYGONTYPE || subtype == CURVEPOLYTYPE) ) - return G_TRUE; + return LW_TRUE; /* Must be a bad combination! */ - return G_FALSE; + return LW_FALSE; } static LWCOLLECTION* lwcollection_from_gserialized_buffer(uchar *data_ptr, uchar g_flags, size_t *g_size) @@ -1044,7 +1044,7 @@ static int gserialized_calculate_gbox_geocentric_from_collection(uchar *data_ptr uchar *start_ptr = data_ptr; int ngeoms = 0; int i; - int first = G_TRUE; + int first = LW_TRUE; int result = G_FAILURE; assert(data_ptr); @@ -1065,7 +1065,7 @@ static int gserialized_calculate_gbox_geocentric_from_collection(uchar *data_ptr if ( first ) { gbox_duplicate(&subbox, gbox); - first = G_FALSE; + first = LW_FALSE; } else { diff --git a/liblwgeom/libgeom.h b/liblwgeom/libgeom.h index 6305be938..ff8e66bbe 100644 --- a/liblwgeom/libgeom.h +++ b/liblwgeom/libgeom.h @@ -28,8 +28,8 @@ #define G_FAILURE 0 #define G_SUCCESS 1 -#define G_TRUE 1 -#define G_FALSE 0 +#define LW_TRUE 1 +#define LW_FALSE 0 /** * Maximum allowed SRID value. @@ -395,6 +395,11 @@ extern GBOX* gbox_new(uchar flags); */ extern int gbox_merge(GBOX *merged_box, GBOX *new_box); +/** +* Update the #GBOX to be large enough to include itself and the new point. +*/ +extern int gbox_merge_point3d(GBOX *gbox, POINT3D *p); + /** * Allocate a string representation of the #GBOX, based on dimensionality of flags. */ @@ -406,7 +411,7 @@ extern char* gbox_to_string(GBOX *gbox); extern GBOX* gbox_copy(GBOX *gbox); /** -* Return #G_TRUE if the #GBOX overlaps, #G_FALSE otherwise. +* Return #LW_TRUE if the #GBOX overlaps, #LW_FALSE otherwise. */ extern int gbox_overlaps(GBOX *g1, GBOX *g2); diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c index 482bc6c8f..002bf3179 100644 --- a/liblwgeom/lwgeodetic.c +++ b/liblwgeom/lwgeodetic.c @@ -11,6 +11,132 @@ #include "lwgeodetic.h" +/** +* Convert spherical coordinates to cartesion coordinates on unit sphere +*/ +void inline geog2cart(GEOGRAPHIC_POINT g, POINT3D *p) +{ + p->x = cos(g.lat) * cos(g.lon); + p->y = cos(g.lat) * sin(g.lon); + p->z = sin(g.lat); +} + +/** +* Convert cartesion coordinates to spherical coordinates on unit sphere +*/ +void inline cart2geog(POINT3D p, GEOGRAPHIC_POINT *g) +{ + g->lon = atan2(p.x, p.y); + g->lat = asin(p.z); +} + +/** +* Calculate the dot product of two unit vectors +*/ +double inline dot_product(POINT3D p1, POINT3D p2) +{ + return (p1.x*p2.x) + (p1.y*p2.y) + (p1.z*p2.z); +} + +void inline unit_normal(POINT3D a, POINT3D b, POINT3D *n) +{ + double d; + double x = a.y * b.z - a.z * b.y; + double y = a.z * b.x - a.x * b.z; + double z = a.x * b.y - a.y * b.x; + d = sqrt(x*x + y*y + z*z); + n->x = x/d; + n->y = y/d; + n->z = z/d; + return; +} + +/** +* Normalize to a unit vector. +*/ +void normalize(POINT3D *p) +{ + double d = sqrt(p->x*p->x + p->y*p->y + p->z*p->z); + if(FP_IS_ZERO(d)) + { + p->x = p->y = p->z = 0.0; + return; + } + p->x = p->x / d; + p->y = p->y / d; + p->z = p->z / d; +} + +/** +* Computes the cross product of two unit vectors using their lat, lng representations. +* Good for small distances between p and q. +*/ +void robust_cross_product(GEOGRAPHIC_POINT p, GEOGRAPHIC_POINT q, POINT3D *a) +{ + double lon_qpp = (q.lon + p.lon) / -2.0; + double lon_qmp = (q.lon - p.lon) / 2.0; + a->x = sin(p.lat-q.lat) * sin(lon_qpp) * cos(lon_qmp) - + sin(p.lat+q.lat) * cos(lon_qpp) * sin(lon_qmp); + a->y = sin(p.lat-q.lat) * cos(lon_qpp) * cos(lon_qmp) + + sin(p.lat+q.lat) * sin(lon_qpp) * sin(lon_qmp); + a->z = cos(p.lat) * cos(q.lat) * sin(q.lon-p.lon); + normalize(a); +} + +void x_to_z(POINT3D *p) +{ + double tmp = p->z; + p->z = p->x; + p->x = tmp; +} + +void y_to_z(POINT3D *p) +{ + double tmp = p->z; + p->z = p->y; + p->y = tmp; +} + +/** +* Returns true if the point p is on the great circle plane. +* Forms the scalar triple product of A,B,p and if the volume of the +* resulting parallelepiped is near zero the point p is on the +* great circle plane. +*/ +int edge_point_on_plane(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p) +{ + POINT3D normal, pt; + double w; + robust_cross_product(e.start, e.end, &normal); + geog2cart(p, &pt); + w = dot_product(normal, pt); + if( FP_IS_ZERO(w) ) + return LW_FALSE; + return LW_TRUE; +} + +/** +* True if the longitude of p is within the range of the longitude of the ends of e +*/ +int edge_contains_longitude(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p, int flipped_longitude) +{ + if( flipped_longitude ) + { + if( p.lon > e.start.lon && p.lon > e.end.lon ) + return LW_TRUE; + if( p.lon < e.start.lon && p.lon < e.end.lon ) + return LW_TRUE; + } + else + { + if( e.start.lon < p.lon && p.lon < e.end.lon ) + return LW_TRUE; + if( e.end.lon < p.lon && p.lon < e.start.lon ) + return LW_TRUE; + } + return LW_FALSE; +} + /** * Given two points on a unit sphere, calculate their distance apart. @@ -41,23 +167,298 @@ double sphere_direction(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e) return heading; } +/** +* Returns true if the point p is on the minor edge defined by the +* end points of e. +*/ +int edge_contains_point(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p, int flipped_longitude) +{ + if( edge_point_on_plane(e, p) && edge_contains_longitude(e, p, flipped_longitude) ) + return LW_TRUE; + return LW_FALSE; +} + +/** +* Used in great circle to compute the pole of the great circle. +*/ +double z_to_latitude(double z) +{ + double sign = signum(z); + double tlat = acos(z); + if(fabs(tlat) > (PI/2.0) ) + { + tlat = sign * (PI - fabs(tlat)); + } + else + { + tlat = sign * tlat; + } + return tlat; +} + +/** +* Computes the pole of the great circle disk which is the intersection of +* the great circle with the line of maximum/minimum gradiant that lies on +* the great circle plane. +*/ +int clairaut_cartesian(POINT3D start, POINT3D end, int top, GEOGRAPHIC_POINT *g) +{ + POINT3D t1, t2; + GEOGRAPHIC_POINT vN1, vN2; + unit_normal(start, end, &t1); + unit_normal(end, start, &t2); + cart2geog(t1, &vN1); + cart2geog(t2, &vN2); + if( top ) + { + g->lat = z_to_latitude(t1.z); + g->lon = vN2.lon; + } + else + { + g->lat = z_to_latitude(t2.z); + g->lon = vN1.lon; + } + return G_SUCCESS; +} + +/** +* Computes the pole of the great circle disk which is the intersection of +* the great circle with the line of maximum/minimum gradiant that lies on +* the great circle plane. +*/ +int clairaut_geographic(GEOGRAPHIC_POINT start, GEOGRAPHIC_POINT end, int top, GEOGRAPHIC_POINT *g) +{ + POINT3D t1, t2; + GEOGRAPHIC_POINT vN1, vN2; + robust_cross_product(start, end, &t1); + robust_cross_product(end, start, &t2); + cart2geog(t1, &vN1); + cart2geog(t2, &vN2); + if( top ) + { + g->lat = z_to_latitude(t1.z); + g->lon = vN2.lon; + } + else + { + g->lat = z_to_latitude(t2.z); + g->lon = vN1.lon; + } + return G_SUCCESS; +} + + /** * Given a starting location r, a distance and an azimuth * to the new point, compute the location of the projected point on the unit sphere. */ -void sphere_project(GEOGRAPHIC_POINT r, double distance, double azimuth, GEOGRAPHIC_POINT *n) +int sphere_project(GEOGRAPHIC_POINT r, double distance, double azimuth, GEOGRAPHIC_POINT *n) { double d = distance; double lat1 = r.lat; - n->lat = asin(sin(lat1) * cos(d) + - cos(lat1) * sin(d) * cos(azimuth)); double a = cos(lat1) * cos(d) - sin(lat1) * sin(d) * cos(azimuth); double b = signum(d) * sin(azimuth); + n->lat = asin(sin(lat1) * cos(d) + + cos(lat1) * sin(d) * cos(azimuth)); n->lon = atan(b/a) + r.lon; + return G_SUCCESS; } -void sphere_gbox(GEOGRAPHIC_EDGE e, GBOX *gbox) + + +int edge_calculate_gbox(GEOGRAPHIC_EDGE e, GBOX *gbox) { + double deltaLongitude; + double distance = sphere_distance(e.start, e.end); + int flipped_longitude = LW_FALSE; + int gimbal_lock = LW_FALSE; + POINT3D p, start, end, startXZ, endXZ, startYZ, endYZ, nT1B, nT2B; + GEOGRAPHIC_EDGE g; + GEOGRAPHIC_POINT vT1, vT2; + + /* Initialize our working copy of the edge */ + g = e; + + /* Initialize box with the start and end points of the edge. */ + geog2cart(g.start, &start); + geog2cart(g.end, &end); + gbox->xmin = FP_MIN(start.x, end.x); + gbox->ymin = FP_MIN(start.y, end.y); + gbox->zmin = FP_MIN(start.z, end.z); + gbox->xmax = FP_MAX(start.x, end.x); + gbox->ymax = FP_MAX(start.y, end.y); + gbox->zmax = FP_MAX(start.z, end.z); + + /* Edge is zero length, just return the naive box */ + if( FP_IS_ZERO(distance) ) + return G_SUCCESS; + + /* Edge is antipodal (one point on each side of the globe), + set the box to contain the whole world and return */ + if( FP_EQUALS(distance, PI) ) + { + gbox->xmin = gbox->ymin = gbox->zmin = -1.0; + gbox->xmax = gbox->ymax = gbox->zmax = 1.0; + return G_SUCCESS; + } + + /* Calculate the difference in longitude between the two points. */ + if( signum(g.start.lon) == signum(g.end.lon) ) + { + deltaLongitude = fabs(g.start.lon - g.end.lon); + } + else + { + double dl = fabs(g.start.lon) + fabs(g.end.lon); + /* Less then a hemisphere apart */ + if( dl < PI ) + { + deltaLongitude = dl; + } + /* Exactly a hemisphere apart */ + else if ( FP_EQUALS( dl, PI ) ) + { + deltaLongitude = PI; + } + /* More than a hemisphere apart, return the other half of the sphere + and note that we are crossing the dateline */ + else + { + flipped_longitude = LW_TRUE; + deltaLongitude = dl - PI; + } + } + + /* If we are crossing the dateline, flip the calculation to the other + side of the globe. We'll flip our output box back at the end of the + calculation. */ + if ( flipped_longitude ) + { + if ( g.start.lon > 0.0 ) + g.start.lon -= PI; + else + g.start.lon += PI; + if ( g.end.lon > 0.0 ) + g.end.lon -= PI; + else + g.end.lon += PI; + } + + /* Check for pole crossings. */ + if( FP_EQUALS(deltaLongitude, PI) ) + { + /* Crosses the north pole, adjust box to contain pole */ + if( (g.start.lat + g.end.lat) > 0.0 ) + { + gbox->zmax = 1.0; + } + /* Crosses the south pole, adjust box to contain pole */ + else + { + gbox->zmin = -1.0; + } + } + /* How about maximal latitudes in this great circle. Are any + of them contained within this arc? */ + else + { + clairaut_cartesian(start, end, LW_TRUE, &vT1); + clairaut_cartesian(start, end, LW_FALSE, &vT2); + if( edge_contains_point(g, vT1, flipped_longitude) ) + { + geog2cart(vT1, &p); + gbox_merge_point3d(gbox, &p); + } + else if ( edge_contains_point(g, vT2, flipped_longitude) ) + { + geog2cart(vT2, &p); + gbox_merge_point3d(gbox, &p); + } + } + + /* Flip the X axis to Z and check for maximal latitudes again. */ + startXZ = start; + x_to_z(&startXZ); + endXZ = end; + x_to_z(&endXZ); + clairaut_cartesian(startXZ, endXZ, LW_TRUE, &vT1); + clairaut_cartesian(startXZ, endXZ, LW_FALSE, &vT2); + gimbal_lock = LW_FALSE; + if ( FP_IS_ZERO(vT1.lat) ) + { + gimbal_lock = LW_TRUE; + } + geog2cart(vT1, &nT1B); + geog2cart(vT2, &nT2B); + x_to_z(&nT1B); + x_to_z(&nT2B); + cart2geog(nT1B, &vT1); + cart2geog(nT2B, &vT2); + if( gimbal_lock ) + { + vT1.lon = PI; + vT2.lon = -1.0 * PI; + } + if( edge_contains_point(g, vT1, flipped_longitude) ) + { + geog2cart(vT1, &p); + gbox_merge_point3d(gbox, &p); + } + else if ( edge_contains_point(g, vT2, flipped_longitude) ) + { + geog2cart(vT2, &p); + gbox_merge_point3d(gbox, &p); + } + + /* Flip the Y axis to Z and check for maximal latitudes again. */ + startYZ = start; + y_to_z(&startYZ); + endYZ = end; + y_to_z(&endYZ); + clairaut_cartesian(startYZ, endYZ, LW_TRUE, &vT1); + clairaut_cartesian(startYZ, endYZ, LW_FALSE, &vT2); + gimbal_lock = LW_FALSE; + if ( FP_IS_ZERO(vT1.lat) ) + { + gimbal_lock = LW_TRUE; + } + geog2cart(vT1, &nT1B); + geog2cart(vT2, &nT2B); + y_to_z(&nT1B); + y_to_z(&nT2B); + cart2geog(nT1B, &vT1); + cart2geog(nT2B, &vT2); + if( gimbal_lock ) + { + vT1.lon = PI; + vT2.lon = -1.0 * PI; + } + if( edge_contains_point(g, vT1, flipped_longitude) ) + { + geog2cart(vT1, &p); + gbox_merge_point3d(gbox, &p); + } + else if ( edge_contains_point(g, vT2, flipped_longitude) ) + { + geog2cart(vT2, &p); + gbox_merge_point3d(gbox, &p); + } + + /* Our cartesian gbox is complete! + If we flipped our longitudes, now we have to flip our cartesian box. */ + if (flipped_longitude) + { + double tmp; + tmp = gbox->xmax; + gbox->xmax = -1.0 * gbox->xmin; + gbox->xmin = -1.0 * tmp; + tmp = gbox->ymax; + gbox->ymax = -1.0 * gbox->ymin; + gbox->ymin = -1.0 * tmp; + } + + return G_SUCCESS; } @@ -152,7 +553,7 @@ static int lwpolygon_calculate_gbox_geodetic(LWPOLY *poly, GBOX *gbox) { GBOX ringbox; int i; - int first = G_TRUE; + int first = LW_TRUE; assert(poly); if( poly->nrings == 0 ) return G_FAILURE; @@ -164,7 +565,7 @@ static int lwpolygon_calculate_gbox_geodetic(LWPOLY *poly, GBOX *gbox) if( first ) { gbox_duplicate(gbox, &ringbox); - first = G_FALSE; + first = LW_FALSE; } else { @@ -179,7 +580,7 @@ static int lwcollection_calculate_gbox_geodetic(LWCOLLECTION *coll, GBOX *gbox) GBOX subbox; int i; int result = G_FAILURE; - int first = G_TRUE; + int first = LW_TRUE; assert(coll); if( coll->ngeoms == 0 ) return G_FAILURE; @@ -197,7 +598,7 @@ static int lwcollection_calculate_gbox_geodetic(LWCOLLECTION *coll, GBOX *gbox) if( first ) { gbox_duplicate(gbox, &subbox); - first = G_FALSE; + first = LW_FALSE; } else { @@ -255,10 +656,10 @@ static int ptarray_check_geodetic(POINTARRAY *pa) getPoint2d_p_ro(pa, t, &pt); //printf( "%d (%g, %g)\n", t, pt->x, pt->y); if ( pt->x < -180.0 || pt->y < -90.0 || pt->x > 180.0 || pt->y > 90.0 ) - return G_FALSE; + return LW_FALSE; } - return G_TRUE; + return LW_TRUE; } static int lwpoint_check_geodetic(LWPOINT *point) @@ -280,10 +681,10 @@ static int lwpoly_check_geodetic(LWPOLY *poly) for( i = 0; i < poly->nrings; i++ ) { - if( ptarray_check_geodetic(poly->rings[i]) == G_FALSE ) - return G_FALSE; + if( ptarray_check_geodetic(poly->rings[i]) == LW_FALSE ) + return LW_FALSE; } - return G_TRUE; + return LW_TRUE; } static int lwcollection_check_geodetic(LWCOLLECTION *col) @@ -293,10 +694,10 @@ static int lwcollection_check_geodetic(LWCOLLECTION *col) for( i = 0; i < col->ngeoms; i++ ) { - if( lwgeom_check_geodetic(col->geoms[i]) == G_FALSE ) - return G_FALSE; + if( lwgeom_check_geodetic(col->geoms[i]) == LW_FALSE ) + return LW_FALSE; } - return G_TRUE; + return LW_TRUE; } int lwgeom_check_geodetic(const LWGEOM *geom) @@ -317,7 +718,7 @@ int lwgeom_check_geodetic(const LWGEOM *geom) default: lwerror("unsupported input geometry type: %d", TYPE_GETTYPE(geom->type)); } - return G_FALSE; + return LW_FALSE; } diff --git a/liblwgeom/lwgeodetic.h b/liblwgeom/lwgeodetic.h index bbdf586c8..09c6bf837 100644 --- a/liblwgeom/lwgeodetic.h +++ b/liblwgeom/lwgeodetic.h @@ -24,8 +24,8 @@ typedef struct { * Two-point great circle segment from a to b. */ typedef struct { - GEOGRAPHIC_POINT a; - GEOGRAPHIC_POINT b; + GEOGRAPHIC_POINT start; + GEOGRAPHIC_POINT end; } GEOGRAPHIC_EDGE; /** @@ -42,9 +42,23 @@ typedef struct { /* ** Prototypes for internal functions. */ -double sphere_distance(GEOGRAPHIC_POINT a, GEOGRAPHIC_POINT b); -double sphere_direction(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e); -void sphere_project(GEOGRAPHIC_POINT r, double distance, double azimuth, GEOGRAPHIC_POINT *n); -void sphere_gbox(GEOGRAPHIC_EDGE e, GBOX *gbox); +void inline geog2cart(GEOGRAPHIC_POINT g, POINT3D *p); +void inline cart2geog(POINT3D p, GEOGRAPHIC_POINT *g); +double inline dot_product(POINT3D p1, POINT3D p2); +void inline unit_normal(POINT3D a, POINT3D b, POINT3D *n); +void normalize(POINT3D *p); +void robust_cross_product(GEOGRAPHIC_POINT p, GEOGRAPHIC_POINT q, POINT3D *a); +void x_to_z(POINT3D *p); +void y_to_z(POINT3D *p); +int edge_point_on_plane(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p); +int edge_contains_longitude(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p, int flipped_longitude); +double sphere_distance(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e); +double sphere_direction(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e); +int edge_contains_point(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p, int flipped_longitude); +double z_to_latitude(double z); +int clairaut_cartesian(POINT3D start, POINT3D end, int top, GEOGRAPHIC_POINT *g); +int clairaut_geographic(GEOGRAPHIC_POINT start, GEOGRAPHIC_POINT end, int top, GEOGRAPHIC_POINT *g); +int sphere_project(GEOGRAPHIC_POINT r, double distance, double azimuth, GEOGRAPHIC_POINT *n); +int edge_calculate_gbox(GEOGRAPHIC_EDGE e, GBOX *gbox); diff --git a/postgis/geography_gist.c b/postgis/geography_gist.c index a39ab6348..b6bd9c121 100644 --- a/postgis/geography_gist.c +++ b/postgis/geography_gist.c @@ -576,7 +576,7 @@ Datum geography_overlaps(PG_FUNCTION_ARGS) geography_datum_gidx(PG_GETARG_DATUM(0), gbox1); geography_datum_gidx(PG_GETARG_DATUM(1), gbox2); - if ( gidx_overlaps(gbox1, gbox2) == G_TRUE ) + if ( gidx_overlaps(gbox1, gbox2) == LW_TRUE ) PG_RETURN_BOOL(TRUE); PG_RETURN_BOOL(FALSE); diff --git a/postgis/geography_inout.c b/postgis/geography_inout.c index 089b57608..988f08b7c 100644 --- a/postgis/geography_inout.c +++ b/postgis/geography_inout.c @@ -223,7 +223,7 @@ Datum geography_in(PG_FUNCTION_ARGS) lwgeom = lwgeom_from_gserialized(g_ser); /* Check if the geography has valid coordinate range. */ - if( lwgeom_check_geodetic(lwgeom) == G_FALSE ) + if( lwgeom_check_geodetic(lwgeom) == LW_FALSE ) { ereport(ERROR, ( errcode(ERRCODE_INVALID_PARAMETER_VALUE), @@ -659,7 +659,7 @@ Datum geography_from_geometry(PG_FUNCTION_ARGS) lwgeom = lwgeom_from_gserialized(g_ser); /* Check if the geography has valid coordinate range. */ - if( lwgeom_check_geodetic(lwgeom) == G_FALSE ) + if( lwgeom_check_geodetic(lwgeom) == LW_FALSE ) { ereport(ERROR, ( errcode(ERRCODE_INVALID_PARAMETER_VALUE), -- 2.50.0