← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/postgis/oneiric-proposed

« back to all changes in this revision

Viewing changes to liblwgeom/lwgeodetic.c

  • Committer: Bazaar Package Importer
  • Author(s): Alan Boudreault
  • Date: 2010-09-29 09:16:10 UTC
  • mfrom: (1.1.12 upstream)
  • Revision ID: james.westby@ubuntu.com-20100929091610-vj4efw8woq34hdn7
Tags: 1.5.2-1
* New upstream release, with a few bug fixes.
* Added shp2pgsql-gui binary.
* Removed patches, applied upstream: getopt.    

Show diffs side-by-side

added added

removed removed

Lines of Context:
liblwgeom/lwgeodetic.c
1
1
/**********************************************************************
2
 
 * $Id: lwgeodetic.c 5344 2010-02-25 15:13:39Z pramsey $
 
2
 * $Id: lwgeodetic.c 5851 2010-08-19 19:48:52Z pramsey $
3
3
 *
4
4
 * PostGIS - Spatial Types for PostgreSQL
5
5
 * Copyright 2009 Paul Ramsey <pramsey@cleverelephant.ca>
158
158
*/
159
159
static int gbox_check_poles(GBOX *gbox)
160
160
{
 
161
        int rv = LW_FALSE;
 
162
        LWDEBUG(4, "checking poles");
 
163
        LWDEBUGF(4, "gbox %s", gbox_to_string(gbox));
161
164
        /* Z axis */
162
165
        if ( gbox->xmin < 0.0 && gbox->xmax > 0.0 &&
163
 
                gbox->ymin < 0.0 && gbox->ymax > 0.0 )
 
166
             gbox->ymin < 0.0 && gbox->ymax > 0.0 )
164
167
        {
165
 
                if ( (gbox->zmin+gbox->zmin)/2 > 0.0 )
 
168
                if ( (gbox->zmin + gbox->zmax) > 0.0 )
 
169
                {
 
170
                        LWDEBUG(4, "enclosed positive z axis");
166
171
                        gbox->zmax = 1.0;
 
172
                }
167
173
                else
 
174
                {
 
175
                        LWDEBUG(4, "enclosed negative z axis");
168
176
                        gbox->zmin = -1.0;
169
 
                return LW_TRUE;
 
177
                }
 
178
                rv = LW_TRUE;
170
179
        }
171
180
 
172
181
        /* Y axis */
173
182
        if ( gbox->xmin < 0.0 && gbox->xmax > 0.0 &&
174
 
                gbox->zmin < 0.0 && gbox->zmax > 0.0 )
 
183
             gbox->zmin < 0.0 && gbox->zmax > 0.0 )
175
184
        {
176
 
                if ( (gbox->ymin+gbox->ymin)/2 > 0.0 )
 
185
                if ( (gbox->ymin + gbox->ymax) > 0.0 )
 
186
                {
 
187
                        LWDEBUG(4, "enclosed positive y axis");
177
188
                        gbox->ymax = 1.0;
 
189
                }
178
190
                else
 
191
                {
 
192
                        LWDEBUG(4, "enclosed negative y axis");
179
193
                        gbox->ymin = -1.0;
180
 
                return LW_TRUE;
 
194
                }
 
195
                rv = LW_TRUE;
181
196
        }
182
197
 
183
198
        /* X axis */
184
199
        if ( gbox->ymin < 0.0 && gbox->ymax > 0.0 &&
185
 
                gbox->zmin < 0.0 && gbox->zmax > 0.0 )
 
200
             gbox->zmin < 0.0 && gbox->zmax > 0.0 )
186
201
        {
187
 
                if ( (gbox->xmin+gbox->xmin)/2 > 0.0 )
 
202
                if ( (gbox->xmin + gbox->xmax) > 0.0 )
 
203
                {
 
204
                        LWDEBUG(4, "enclosed positive x axis");
188
205
                        gbox->xmax = 1.0;
 
206
                }
189
207
                else
 
208
                {
 
209
                        LWDEBUG(4, "enclosed negative x axis");
190
210
                        gbox->xmin = -1.0;
191
 
                return LW_TRUE;
 
211
                }
 
212
                rv = LW_TRUE;
192
213
        }
193
214
 
194
 
        return LW_FALSE;
 
215
        return rv;
195
216
}
196
217
 
197
218
 
786
807
                LWDEBUG(4, "flipping point to other side of sphere");
787
808
                g->lat = -1.0 * g->lat;
788
809
                g->lon = g->lon + M_PI;
 
810
                LWDEBUGF(4, "g == GPOINT(%.12g %.12g)", g->lat, g->lon);
 
811
                LWDEBUGF(4, "g == POINT(%.12g %.12g)", rad2deg(g->lon), rad2deg(g->lat));
789
812
                if ( g->lon > M_PI )
790
813
                {
791
814
                        g->lon = -1.0 * (2.0 * M_PI - g->lon);
1248
1271
*/
1249
1272
void gbox_pt_outside(const GBOX *gbox, POINT2D *pt_outside)
1250
1273
{
1251
 
        static double grow = M_PI / 180.0 / 60.0; /* one arc-minute */
1252
 
        double d;
 
1274
        double grow = M_PI / 180.0 / 60.0; /* one arc-minute */
1253
1275
        int i;
1254
1276
        GBOX ge;
1255
1277
        POINT3D corners[8];
1256
1278
        POINT3D pt;
1257
1279
        GEOGRAPHIC_POINT g;
1258
1280
 
1259
 
        /* Assign our box and expand it slightly. */
1260
 
        ge = *gbox;
1261
 
        ge.xmin -= grow;
1262
 
        ge.ymin -= grow;
1263
 
        ge.zmin -= grow;
1264
 
        ge.xmax += grow;
1265
 
        ge.ymax += grow;
1266
 
        ge.zmax += grow;
1267
 
 
1268
 
        /* Build our eight corner points */
1269
 
        corners[0].x = ge.xmin;
1270
 
        corners[0].y = ge.ymin;
1271
 
        corners[0].z = ge.zmin;
1272
 
 
1273
 
        corners[1].x = ge.xmin;
1274
 
        corners[1].y = ge.ymax;
1275
 
        corners[1].z = ge.zmin;
1276
 
 
1277
 
        corners[2].x = ge.xmin;
1278
 
        corners[2].y = ge.ymin;
1279
 
        corners[2].z = ge.zmax;
1280
 
 
1281
 
        corners[3].x = ge.xmax;
1282
 
        corners[3].y = ge.ymin;
1283
 
        corners[3].z = ge.zmin;
1284
 
 
1285
 
        corners[4].x = ge.xmax;
1286
 
        corners[4].y = ge.ymax;
1287
 
        corners[4].z = ge.zmin;
1288
 
 
1289
 
        corners[5].x = ge.xmax;
1290
 
        corners[5].y = ge.ymin;
1291
 
        corners[5].z = ge.zmax;
1292
 
 
1293
 
        corners[6].x = ge.xmin;
1294
 
        corners[6].y = ge.ymax;
1295
 
        corners[6].z = ge.zmax;
1296
 
 
1297
 
        corners[7].x = ge.xmax;
1298
 
        corners[7].y = ge.ymax;
1299
 
        corners[7].z = ge.zmax;
1300
 
 
1301
 
        for ( i = 0; i < 8; i++ )
 
1281
        while( grow < M_PI ) 
1302
1282
        {
1303
 
                normalize(&(corners[i]));
1304
 
                if ( ! gbox_contains_point3d(gbox, &(corners[i])) )
 
1283
                /* Assign our box and expand it slightly. */
 
1284
                ge = *gbox;
 
1285
                if ( ge.xmin > -1 ) ge.xmin -= grow;
 
1286
                if ( ge.ymin > -1 ) ge.ymin -= grow;
 
1287
                if ( ge.zmin > -1 ) ge.zmin -= grow;
 
1288
                if ( ge.xmax < 1 )  ge.xmax += grow;
 
1289
                if ( ge.ymax < 1 )  ge.ymax += grow;
 
1290
                if ( ge.zmax < 1 )  ge.zmax += grow;
 
1291
 
 
1292
                /* Build our eight corner points */
 
1293
                corners[0].x = ge.xmin;
 
1294
                corners[0].y = ge.ymin;
 
1295
                corners[0].z = ge.zmin;
 
1296
 
 
1297
                corners[1].x = ge.xmin;
 
1298
                corners[1].y = ge.ymax;
 
1299
                corners[1].z = ge.zmin;
 
1300
 
 
1301
                corners[2].x = ge.xmin;
 
1302
                corners[2].y = ge.ymin;
 
1303
                corners[2].z = ge.zmax;
 
1304
 
 
1305
                corners[3].x = ge.xmax;
 
1306
                corners[3].y = ge.ymin;
 
1307
                corners[3].z = ge.zmin;
 
1308
 
 
1309
                corners[4].x = ge.xmax;
 
1310
                corners[4].y = ge.ymax;
 
1311
                corners[4].z = ge.zmin;
 
1312
 
 
1313
                corners[5].x = ge.xmax;
 
1314
                corners[5].y = ge.ymin;
 
1315
                corners[5].z = ge.zmax;
 
1316
 
 
1317
                corners[6].x = ge.xmin;
 
1318
                corners[6].y = ge.ymax;
 
1319
                corners[6].z = ge.zmax;
 
1320
 
 
1321
                corners[7].x = ge.xmax;
 
1322
                corners[7].y = ge.ymax;
 
1323
                corners[7].z = ge.zmax;
 
1324
 
 
1325
                LWDEBUG(4, "trying to use a box corner point...");
 
1326
                for ( i = 0; i < 8; i++ )
1305
1327
                {
1306
 
                        pt = corners[i];
1307
 
                        normalize(&pt);
1308
 
                        cart2geog(&pt, &g);
1309
 
                        pt_outside->x = rad2deg(g.lon);
1310
 
                        pt_outside->y = rad2deg(g.lat);
1311
 
                        return;
 
1328
                        normalize(&(corners[i]));
 
1329
                        LWDEBUGF(4, "testing corner %d: POINT(%.8g %.8g %.8g)", i, corners[i].x, corners[i].y, corners[i].z);
 
1330
                        if ( ! gbox_contains_point3d(gbox, &(corners[i])) )
 
1331
                        {
 
1332
                                LWDEBUGF(4, "corner %d is outside our gbox", i);
 
1333
                                pt = corners[i];
 
1334
                                normalize(&pt);
 
1335
                                cart2geog(&pt, &g);
 
1336
                                pt_outside->x = rad2deg(g.lon);
 
1337
                                pt_outside->y = rad2deg(g.lat);
 
1338
                                LWDEBUGF(4, "returning POINT(%.8g %.8g) as outside point", pt_outside->x, pt_outside->y);
 
1339
                                return;
 
1340
                        }
1312
1341
                }
1313
 
        }
1314
 
 
1315
 
        pt.x = 1.0;
1316
 
        pt.y = 0.0;
1317
 
        pt.z = 0.0;
1318
 
 
1319
 
        if ((1.0 - gbox->xmax) > 0.1)
1320
 
        {
1321
 
                pt.x = gbox->xmax + (1.0 - gbox->xmax) * 0.01;
1322
 
                d = sqrt((1.0 - POW2(pt.x))/2.0);
1323
 
                pt.y = d;
1324
 
                pt.z = d;
1325
 
        }
1326
 
        else if ((1.0 - gbox->ymax) > 0.1)
1327
 
        {
1328
 
                pt.y = gbox->ymax + (1.0 - gbox->ymax) * 0.01;
1329
 
                d = sqrt((1.0 - POW2(pt.y))/2.0);
1330
 
                pt.x = d;
1331
 
                pt.z = d;
1332
 
        }
1333
 
        else if ((1.0 - gbox->zmax) > 0.1)
1334
 
        {
1335
 
                pt.z = gbox->zmax + (1.0 - gbox->zmax) * 0.01;
1336
 
                d = sqrt((1.0 - POW2(pt.z))/2.0);
1337
 
                pt.x = d;
1338
 
                pt.y = d;
1339
 
        }
1340
 
        normalize(&pt);
1341
 
        cart2geog(&pt, &g);
1342
 
        pt_outside->x = rad2deg(g.lon);
1343
 
        pt_outside->y = rad2deg(g.lat);
 
1342
                
 
1343
                /* Try a wider growth to push the corners outside the original box. */
 
1344
                grow *= 2.0;
 
1345
        }
 
1346
 
 
1347
        /* This should never happen! */
 
1348
        lwerror("BOOM! Could not generate outside point!");
1344
1349
        return;
1345
1350
}
1346
1351
 
1421
1426
 
1422
1427
 
1423
1428
/**
1424
 
* This routine returns LW_TRUE if the point is inside the ring or on the boundary, LW_FALSE otherwise.
 
1429
* This routine returns LW_TRUE if the stabline joining the pt_outside and pt_to_test 
 
1430
* crosses the ring an odd number of times, or if the pt_to_test is on the ring boundary itself, 
 
1431
* returning LW_FALSE otherwise.
1425
1432
* The pt_outside must be guaranteed to be outside the ring (use the geography_pt_outside() function
1426
1433
* to derive one in postgis, or the gbox_pt_outside() function if you don't mind burning CPU cycles
1427
1434
* building a gbox first).
2067
2074
                return LW_FALSE;
2068
2075
        }
2069
2076
 
 
2077
        LWDEBUGF(4, "pt_to_test POINT(%.18g %.18g)", pt_to_test->x, pt_to_test->y);
 
2078
        LWDEBUGF(4, "gbox %s", gbox_to_string(gbox));
 
2079
 
2070
2080
        /* Point not in box? Done! */
2071
2081
        geographic_point_init(pt_to_test->x, pt_to_test->y, &gpt_to_test);
2072
2082
        geog2cart(&gpt_to_test, &p);
2073
 
        if ( ! gbox_contains_point3d(gbox, &p) )
 
2083
        LWDEBUGF(4, "p POINT(%.18g %.18g %.18g)", p.x, p.y, p.z);       
 
2084
        if ( ! gbox_contains_point3d(gbox, &p) ) 
 
2085
        {
 
2086
                LWDEBUG(4, "the point is not in the box!");
2074
2087
                return LW_FALSE;
 
2088
        }
2075
2089
 
2076
2090
        /* Calculate our outside point from the gbox */
2077
2091
        gbox_pt_outside(gbox, &pt_outside);
2078
2092
 
2079
2093
        LWDEBUGF(4, "pt_outside POINT(%.18g %.18g)", pt_outside.x, pt_outside.y);
2080
 
        LWDEBUGF(4, "pt_to_test POINT(%.18g %.18g)", pt_to_test->x, pt_to_test->y);
2081
2094
        LWDEBUGF(4, "polygon %s", lwgeom_to_ewkt((LWGEOM*)poly, PARSER_CHECK_NONE));
2082
 
        LWDEBUGF(4, "gbox %s", gbox_to_string(gbox));
2083
2095
 
2084
2096
        /* Not in outer ring? We're done! */
2085
2097
        if ( ! ptarray_point_in_ring(poly->rings[0], &pt_outside, pt_to_test) )