point_set(-45.0, 60.0, &gs);
point_set(135.0, 60.0, &ge);
- geog2cart(gs, &vs);
- geog2cart(ge, &ve);
+ geog2cart(&gs, &vs);
+ geog2cart(&ge, &ve);
- clairaut_cartesian(vs, ve, &v_out_top, &v_out_bottom);
- clairaut_geographic(gs, ge, &g_out_top, &g_out_bottom);
+ clairaut_cartesian(&vs, &ve, &v_out_top, &v_out_bottom);
+ clairaut_geographic(&gs, &ge, &g_out_top, &g_out_bottom);
CU_ASSERT_DOUBLE_EQUAL(v_out_top.lat, g_out_top.lat, 0.000001);
CU_ASSERT_DOUBLE_EQUAL(v_out_top.lon, g_out_top.lon, 0.000001);
gs.lon = -1.3387392931438733;
ge.lon = 1.80285336044592;
- geog2cart(gs, &vs);
- geog2cart(ge, &ve);
+ geog2cart(&gs, &vs);
+ geog2cart(&ge, &ve);
- clairaut_cartesian(vs, ve, &v_out_top, &v_out_bottom);
- clairaut_geographic(gs, ge, &g_out_top, &g_out_bottom);
+ clairaut_cartesian(&vs, &ve, &v_out_top, &v_out_bottom);
+ clairaut_geographic(&gs, &ge, &g_out_top, &g_out_bottom);
CU_ASSERT_DOUBLE_EQUAL(v_out_top.lat, g_out_top.lat, 0.000001);
CU_ASSERT_DOUBLE_EQUAL(v_out_top.lon, g_out_top.lon, 0.000001);
/* Covers case, end-to-end intersection */
edge_set(50, -10.999999999999998224, -10.0, 50.0, &e1);
edge_set(-10.0, 50.0, -10.272779983831613393, -16.937003313332997578, &e2);
- rv = edge_intersection(e1, e2, &g);
+ rv = edge_intersection(&e1, &e2, &g);
CU_ASSERT_EQUAL(rv, LW_TRUE);
/* Medford case, very short segment vs very long one */
e2.start.lon = -2.1461493501950630325;
e2.end.lat = 0.70971354024834598651;
e2.end.lon = 2.1082194552519770703;
- rv = edge_intersection(e1, e2, &g);
+ rv = edge_intersection(&e1, &e2, &g);
CU_ASSERT_EQUAL(rv, LW_FALSE);
/* Again, this time with a less exact input edge. */
edge_set(-123.165031277506, 42.4696787216231, -123.165031605021, 42.4697127292275, &e1);
- rv = edge_intersection(e1, e2, &g);
+ rv = edge_intersection(&e1, &e2, &g);
CU_ASSERT_EQUAL(rv, LW_FALSE);
/* Second Medford case, very short segment vs very long one
/* Intersection at (0 0) */
edge_set(-1.0, 0.0, 1.0, 0.0, &e1);
edge_set(0.0, -1.0, 0.0, 1.0, &e2);
- rv = edge_intersection(e1, e2, &g);
+ rv = edge_intersection(&e1, &e2, &g);
point_rad2deg(&g);
CU_ASSERT_DOUBLE_EQUAL(g.lat, 0.0, 0.00001);
CU_ASSERT_DOUBLE_EQUAL(g.lon, 0.0, 0.00001);
/* No intersection at (0 0) */
edge_set(-1.0, 0.0, 1.0, 0.0, &e1);
edge_set(0.0, -1.0, 0.0, -2.0, &e2);
- rv = edge_intersection(e1, e2, &g);
+ rv = edge_intersection(&e1, &e2, &g);
CU_ASSERT_EQUAL(rv, LW_FALSE);
/* End touches middle of segment at (0 0) */
edge_set(-1.0, 0.0, 1.0, 0.0, &e1);
edge_set(0.0, -1.0, 0.0, 0.0, &e2);
- rv = edge_intersection(e1, e2, &g);
+ rv = edge_intersection(&e1, &e2, &g);
point_rad2deg(&g);
#if 0
printf("\n");
/* End touches end of segment at (0 0) */
edge_set(0.0, 0.0, 1.0, 0.0, &e1);
edge_set(0.0, -1.0, 0.0, 0.0, &e2);
- rv = edge_intersection(e1, e2, &g);
+ rv = edge_intersection(&e1, &e2, &g);
point_rad2deg(&g);
#if 0
printf("\n");
/* Intersection at (180 0) */
edge_set(-179.0, -1.0, 179.0, 1.0, &e1);
edge_set(-179.0, 1.0, 179.0, -1.0, &e2);
- rv = edge_intersection(e1, e2, &g);
+ rv = edge_intersection(&e1, &e2, &g);
point_rad2deg(&g);
CU_ASSERT_DOUBLE_EQUAL(g.lat, 0.0, 0.00001);
CU_ASSERT_DOUBLE_EQUAL(fabs(g.lon), 180.0, 0.00001);
/* Intersection at (180 0) */
edge_set(-170.0, 0.0, 170.0, 0.0, &e1);
edge_set(180.0, -10.0, 180.0, 10.0, &e2);
- rv = edge_intersection(e1, e2, &g);
+ rv = edge_intersection(&e1, &e2, &g);
point_rad2deg(&g);
CU_ASSERT_DOUBLE_EQUAL(g.lat, 0.0, 0.00001);
CU_ASSERT_DOUBLE_EQUAL(fabs(g.lon), 180.0, 0.00001);
/* Intersection at north pole */
edge_set(-180.0, 80.0, 0.0, 80.0, &e1);
edge_set(90.0, 80.0, -90.0, 80.0, &e2);
- rv = edge_intersection(e1, e2, &g);
+ rv = edge_intersection(&e1, &e2, &g);
point_rad2deg(&g);
CU_ASSERT_DOUBLE_EQUAL(g.lat, 90.0, 0.00001);
CU_ASSERT_EQUAL(rv, LW_TRUE);
/* Equal edges return true */
edge_set(45.0, 10.0, 50.0, 20.0, &e1);
edge_set(45.0, 10.0, 50.0, 20.0, &e2);
- rv = edge_intersection(e1, e2, &g);
+ rv = edge_intersection(&e1, &e2, &g);
point_rad2deg(&g);
CU_ASSERT_EQUAL(rv, LW_TRUE);
/* Parallel edges (same great circle, different end points) return true */
edge_set(40.0, 0.0, 70.0, 0.0, &e1);
edge_set(60.0, 0.0, 50.0, 0.0, &e2);
- rv = edge_intersection(e1, e2, &g);
+ rv = edge_intersection(&e1, &e2, &g);
point_rad2deg(&g);
CU_ASSERT_EQUAL(rv, 2); /* Hack, returning 2 as the 'co-linear' value */
/* End touches arc at north pole */
edge_set(-180.0, 80.0, 0.0, 80.0, &e1);
edge_set(90.0, 80.0, -90.0, 90.0, &e2);
- rv = edge_intersection(e1, e2, &g);
+ rv = edge_intersection(&e1, &e2, &g);
point_rad2deg(&g);
#if 0
printf("\n");
/* End touches end at north pole */
edge_set(-180.0, 80.0, 0.0, 90.0, &e1);
edge_set(90.0, 80.0, -90.0, 90.0, &e2);
- rv = edge_intersection(e1, e2, &g);
+ rv = edge_intersection(&e1, &e2, &g);
point_rad2deg(&g);
#if 0
printf("\n");
/* closest point at origin, one degree away */
edge_set(-50.0, 0.0, 50.0, 0.0, &e);
point_set(0.0, 1.0, &g);
- d = edge_distance_to_point(e, g, 0);
+ d = edge_distance_to_point(&e, &g, 0);
CU_ASSERT_DOUBLE_EQUAL(d, M_PI / 180.0, 0.00001);
/* closest point at origin, one degree away */
edge_set(-50.0, 0.0, 50.0, 0.0, &e);
point_set(0.0, 2.0, &g);
- d = edge_distance_to_point(e, g, &closest);
+ d = edge_distance_to_point(&e, &g, &closest);
#if 0
printf("LINESTRING(%.8g %.8g, %.8g %.8g)\n", e.start.lon, e.start.lat, e.end.lon, e.end.lat);
printf("POINT(%.9g %.9g)\n", g.lon, g.lat);
/* closest point at origin, one degree away */
edge_set(-50.0, 0.0, 50.0, 0.0, &e1);
edge_set(-5.0, 20.0, 0.0, 1.0, &e2);
- d = edge_distance_to_edge(e1, e2, &c1, &c2);
+ d = edge_distance_to_edge(&e1, &e2, &c1, &c2);
#if 0
printf("LINESTRING(%.8g %.8g, %.8g %.8g)\n", e1.start.lon, e1.start.lat, e1.end.lon, e1.end.lat);
printf("LINESTRING(%.8g %.8g, %.8g %.8g)\n", e2.start.lon, e2.start.lat, e2.end.lon, e2.end.lat);
pt_to_test.y = 1.05;
pt_outside.x = 1.2;
pt_outside.y = 1.15;
- result = ptarray_point_in_ring(poly->rings[0], pt_outside, pt_to_test);
+ result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
CU_ASSERT_EQUAL(result, LW_TRUE);
lwgeom_free(lwg);
pt_to_test.y = 1.15;
pt_outside.x = 1.2;
pt_outside.y = 1.2;
- result = ptarray_point_in_ring(poly->rings[0], pt_outside, pt_to_test);
+ result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
CU_ASSERT_EQUAL(result, LW_FALSE);
lwgeom_free(lwg);
pt_to_test.y = 0.9;
pt_outside.x = 1.2;
pt_outside.y = 1.05;
- result = ptarray_point_in_ring(poly->rings[0], pt_outside, pt_to_test);
+ result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
CU_ASSERT_EQUAL(result, LW_FALSE);
lwgeom_free(lwg);
pt_to_test.y = 1.0;
pt_outside.x = 1.0;
pt_outside.y = 10.0;
- result = ptarray_point_in_ring(poly->rings[0], pt_outside, pt_to_test);
+ result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
CU_ASSERT_EQUAL(result, LW_TRUE);
lwgeom_free(lwg);
pt_to_test.y = 1.05;
pt_outside.x = 1.2;
pt_outside.y = 1.05;
- result = ptarray_point_in_ring(poly->rings[0], pt_outside, pt_to_test);
+ result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
CU_ASSERT_EQUAL(result, LW_TRUE);
lwgeom_free(lwg);
pt_to_test.y = 1.0;
pt_outside.x = 1.2;
pt_outside.y = 1.05;
- result = ptarray_point_in_ring(poly->rings[0], pt_outside, pt_to_test);
+ result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
CU_ASSERT_EQUAL(result, LW_TRUE);
lwgeom_free(lwg);
pt_to_test.y = 1.1;
pt_outside.x = 1.2;
pt_outside.y = 1.05;
- result = ptarray_point_in_ring(poly->rings[0], pt_outside, pt_to_test);
+ result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
CU_ASSERT_EQUAL(result, LW_TRUE);
lwgeom_free(lwg);
pt_to_test.y = 1.05;
pt_outside.x = 1.1;
pt_outside.y = 1.3;
- result = ptarray_point_in_ring(poly->rings[0], pt_outside, pt_to_test);
+ result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
CU_ASSERT_EQUAL(result, LW_TRUE);
lwgeom_free(lwg);
pt_to_test.y = 1.0;
pt_outside.x = 1.5;
pt_outside.y = 2.0;
- result = ptarray_point_in_ring(poly->rings[0], pt_outside, pt_to_test);
+ result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
CU_ASSERT_EQUAL(result, LW_FALSE);
lwgeom_free(lwg);
pt_to_test.y = 0.0;
pt_outside.x = 1.0;
pt_outside.y = 2.0;
- result = ptarray_point_in_ring(poly->rings[0], pt_outside, pt_to_test);
+ result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
CU_ASSERT_EQUAL(result, LW_FALSE);
lwgeom_free(lwg);
pt_to_test.y = 50.0;
pt_outside.x = -10.2727799838316134;
pt_outside.y = -16.9370033133329976;
- result = ptarray_point_in_ring(poly->rings[0], pt_outside, pt_to_test);
+ result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
CU_ASSERT_EQUAL(result, LW_TRUE);
lwgeom_free(lwg);
pt_to_test.y = 42.2702301207017328;
pt_outside.x = 120.695136159150778;
pt_outside.y = 40.6920926049588516;
- result = ptarray_point_in_ring(poly->rings[0], pt_outside, pt_to_test);
+ result = ptarray_point_in_ring(poly->rings[0], &pt_outside, &pt_to_test);
CU_ASSERT_EQUAL(result, LW_FALSE);
lwgeom_free(lwg);
#endif
poly = (LWPOLY*)lwg;
pt_to_test.x = -10.0;
pt_to_test.y = 50.0;
- result = lwpoly_covers_point2d(poly, gbox, pt_to_test);
+ result = lwpoly_covers_point2d(poly, &gbox, &pt_to_test);
CU_ASSERT_EQUAL(result, LW_TRUE);
lwgeom_free(lwg);
}
lwg2 = lwgeom_from_ewkt("LINESTRING(-10 -5, -5 0, 5 0, 10 -5)", PARSER_CHECK_NONE);
lwgeom_calculate_gbox_geodetic(lwg1, &gbox1);
lwgeom_calculate_gbox_geodetic(lwg2, &gbox2);
- d = lwgeom_distance_spheroid(lwg1, lwg2, gbox1, gbox2, s, 0.0);
+ d = lwgeom_distance_spheroid(lwg1, lwg2, &gbox1, &gbox2, &s, 0.0);
CU_ASSERT_DOUBLE_EQUAL(d, s.radius * M_PI / 180.0, 0.00001);
lwgeom_free(lwg1);
lwgeom_free(lwg2);
lwg2 = lwgeom_from_ewkt("LINESTRING(-10 -5, -5 20, 5 0, 10 -5)", PARSER_CHECK_NONE);
lwgeom_calculate_gbox_geodetic(lwg1, &gbox1);
lwgeom_calculate_gbox_geodetic(lwg2, &gbox2);
- d = lwgeom_distance_spheroid(lwg1, lwg2, gbox1, gbox2, s, 0.0);
+ d = lwgeom_distance_spheroid(lwg1, lwg2, &gbox1, &gbox2, &s, 0.0);
CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001);
lwgeom_free(lwg1);
lwgeom_free(lwg2);
lwg2 = lwgeom_from_ewkt("LINESTRING(-10 -5, -5 0, 5 0, 10 -5)", PARSER_CHECK_NONE);
lwgeom_calculate_gbox_geodetic(lwg1, &gbox1);
lwgeom_calculate_gbox_geodetic(lwg2, &gbox2);
- d = lwgeom_distance_spheroid(lwg1, lwg2, gbox1, gbox2, s, 0.0);
+ d = lwgeom_distance_spheroid(lwg1, lwg2, &gbox1, &gbox2, &s, 0.0);
CU_ASSERT_DOUBLE_EQUAL(d, s.radius * M_PI / 180.0, 0.00001);
lwgeom_free(lwg1);
lwgeom_free(lwg2);
lwg2 = lwgeom_from_ewkt("POINT(-4 -1)", PARSER_CHECK_NONE);
lwgeom_calculate_gbox_geodetic(lwg1, &gbox1);
lwgeom_calculate_gbox_geodetic(lwg2, &gbox2);
- d = lwgeom_distance_spheroid(lwg1, lwg2, gbox1, gbox2, s, 0.0);
+ d = lwgeom_distance_spheroid(lwg1, lwg2, &gbox1, &gbox2, &s, 0.0);
CU_ASSERT_DOUBLE_EQUAL(d, s.radius * M_PI / 90.0, 0.00001);
lwgeom_free(lwg1);
lwgeom_free(lwg2);
lwg2 = lwgeom_from_ewkt("POINT(-1 -1)", PARSER_CHECK_NONE);
lwgeom_calculate_gbox_geodetic(lwg1, &gbox1);
lwgeom_calculate_gbox_geodetic(lwg2, &gbox2);
- d = lwgeom_distance_spheroid(lwg1, lwg2, gbox1, gbox2, s, 0.0);
+ d = lwgeom_distance_spheroid(lwg1, lwg2, &gbox1, &gbox2, &s, 0.0);
CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001);
lwgeom_free(lwg1);
lwgeom_free(lwg2);
lwg2 = lwgeom_from_ewkt("POINT(-1 -1)", PARSER_CHECK_NONE);
lwgeom_calculate_gbox_geodetic(lwg1, &gbox1);
lwgeom_calculate_gbox_geodetic(lwg2, &gbox2);
- d = lwgeom_distance_spheroid(lwg1, lwg2, gbox1, gbox2, s, 0.0);
+ d = lwgeom_distance_spheroid(lwg1, lwg2, &gbox1, &gbox2, &s, 0.0);
CU_ASSERT_DOUBLE_EQUAL(d, 111178.142466, 0.1);
lwgeom_free(lwg1);
lwgeom_free(lwg2);
lwg2 = lwgeom_from_ewkt("POINT(2 2)", PARSER_CHECK_NONE);
lwgeom_calculate_gbox_geodetic(lwg1, &gbox1);
lwgeom_calculate_gbox_geodetic(lwg2, &gbox2);
- d = lwgeom_distance_spheroid(lwg1, lwg2, gbox1, gbox2, s, 0.0);
+ d = lwgeom_distance_spheroid(lwg1, lwg2, &gbox1, &gbox2, &s, 0.0);
CU_ASSERT_DOUBLE_EQUAL(d, 0.0, 0.00001);
lwgeom_free(lwg1);
lwgeom_free(lwg2);
lwg2 = lwgeom_from_ewkt("0106000020E61000000100000001030000000100000007000000280EC3FB8CCA5EC0A5CDC747233C45402787C8F58CCA5EC0659EA2761E3C45400CED58DF8FCA5EC0C37FAE6E1E3C4540AE97B8E08FCA5EC00346F58B1F3C4540250359FD8ECA5EC05460628E1F3C45403738F4018FCA5EC05DC84042233C4540280EC3FB8CCA5EC0A5CDC747233C4540", PARSER_CHECK_NONE);
lwgeom_calculate_gbox_geodetic(lwg1, &gbox1);
lwgeom_calculate_gbox_geodetic(lwg2, &gbox2);
- d = lwgeom_distance_spheroid(lwg1, lwg2, gbox1, gbox2, s, 0.0);
+ d = lwgeom_distance_spheroid(lwg1, lwg2, &gbox1, &gbox2, &s, 0.0);
CU_ASSERT_DOUBLE_EQUAL(d, 23630.8003, 0.1);
lwgeom_free(lwg1);
lwgeom_free(lwg2);
/* One vertical degree */
point_set(0.0, 0.0, &g1);
point_set(0.0, 1.0, &g2);
- d = spheroid_distance(g1, g2, s);
+ d = spheroid_distance(&g1, &g2, &s);
CU_ASSERT_DOUBLE_EQUAL(d, 110574.388615329, 0.001);
/* Ten horizontal degrees */
point_set(-10.0, 0.0, &g1);
point_set(0.0, 0.0, &g2);
- d = spheroid_distance(g1, g2, s);
+ d = spheroid_distance(&g1, &g2, &s);
CU_ASSERT_DOUBLE_EQUAL(d, 1113194.90793274, 0.001);
/* One horizonal degree */
point_set(-1.0, 0.0, &g1);
point_set(0.0, 0.0, &g2);
- d = spheroid_distance(g1, g2, s);
+ d = spheroid_distance(&g1, &g2, &s);
CU_ASSERT_DOUBLE_EQUAL(d, 111319.490779, 0.001);
/* Around world w/ slight bend */
point_set(-180.0, 0.0, &g1);
point_set(0.0, 1.0, &g2);
- d = spheroid_distance(g1, g2, s);
+ d = spheroid_distance(&g1, &g2, &s);
CU_ASSERT_DOUBLE_EQUAL(d, 19893357.0704483, 0.001);
/* Up to pole */
point_set(-180.0, 0.0, &g1);
point_set(0.0, 90.0, &g2);
- d = spheroid_distance(g1, g2, s);
+ d = spheroid_distance(&g1, &g2, &s);
CU_ASSERT_DOUBLE_EQUAL(d, 10001965.7295318, 0.001);
}
/* Medford lot test polygon */
lwg = lwgeom_from_ewkt("POLYGON((-122.848227067007 42.5007249610493,-122.848309475585 42.5007179884263,-122.848327688675 42.500835880696,-122.848245279942 42.5008428533324,-122.848227067007 42.5007249610493))", PARSER_CHECK_NONE);
lwgeom_calculate_gbox_geodetic(lwg, &gbox);
- a1 = lwgeom_area_sphere(lwg, gbox, s);
- a2 = lwgeom_area_spheroid(lwg, gbox, s);
+ a1 = lwgeom_area_sphere(lwg, &gbox, &s);
+ a2 = lwgeom_area_spheroid(lwg, &gbox, &s);
//printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2);
CU_ASSERT_DOUBLE_EQUAL(a1, a2, 0.2);
/* Big-ass polygon */
lwg = lwgeom_from_ewkt("POLYGON((-2 3, -2 4, -1 4, -1 3, -2 3))", PARSER_CHECK_NONE);
lwgeom_calculate_gbox_geodetic(lwg, &gbox);
- a1 = lwgeom_area_sphere(lwg, gbox, s);
- a2 = lwgeom_area_spheroid(lwg, gbox, s);
+ a1 = lwgeom_area_sphere(lwg, &gbox, &s);
+ a2 = lwgeom_area_spheroid(lwg, &gbox, &s);
//printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2);
CU_ASSERT_DOUBLE_EQUAL(a1, a2, 100000000.0);
return g;
}
-int gbox_merge_point3d(POINT3D p, GBOX *gbox)
+int gbox_merge_point3d(const POINT3D *p, GBOX *gbox)
{
- if( gbox->xmin > p.x ) gbox->xmin = p.x;
- if( gbox->ymin > p.y ) gbox->ymin = p.y;
- if( gbox->zmin > p.z ) gbox->zmin = p.z;
- if( gbox->xmax < p.x ) gbox->xmax = p.x;
- if( gbox->ymax < p.y ) gbox->ymax = p.y;
- if( gbox->zmax < p.z ) gbox->zmax = p.z;
+ if( gbox->xmin > p->x ) gbox->xmin = p->x;
+ if( gbox->ymin > p->y ) gbox->ymin = p->y;
+ if( gbox->zmin > p->z ) gbox->zmin = p->z;
+ if( gbox->xmax < p->x ) gbox->xmax = p->x;
+ if( gbox->ymax < p->y ) gbox->ymax = p->y;
+ if( gbox->zmax < p->z ) gbox->zmax = p->z;
return G_SUCCESS;
}
-int gbox_contains_point3d(GBOX gbox, POINT3D pt)
+int gbox_contains_point3d(const GBOX *gbox, const POINT3D *pt)
{
- if( gbox.xmin > pt.x || gbox.ymin > pt.y || gbox.zmin > pt.z ||
- gbox.xmax < pt.x || gbox.ymax < pt.y || gbox.zmax < pt.z )
+ if( gbox->xmin > pt->x || gbox->ymin > pt->y || gbox->zmin > pt->z ||
+ gbox->xmax < pt->x || gbox->ymax < pt->y || gbox->zmax < pt->z )
{
return LW_FALSE;
}
return LW_TRUE;
}
-int gbox_merge(GBOX new_box, GBOX *merge_box)
+int gbox_merge(const GBOX *new_box, GBOX *merge_box)
{
assert(merge_box);
- if( merge_box->flags != new_box.flags )
+ if( merge_box->flags != new_box->flags )
return G_FAILURE;
- if( new_box.xmin < merge_box->xmin) merge_box->xmin = new_box.xmin;
- if( new_box.ymin < merge_box->ymin) merge_box->ymin = new_box.ymin;
- if( new_box.xmax > merge_box->xmax) merge_box->xmax = new_box.xmax;
- if( new_box.ymax > merge_box->ymax) merge_box->ymax = new_box.ymax;
+ if( new_box->xmin < merge_box->xmin) merge_box->xmin = new_box->xmin;
+ if( new_box->ymin < merge_box->ymin) merge_box->ymin = new_box->ymin;
+ if( new_box->xmax > merge_box->xmax) merge_box->xmax = new_box->xmax;
+ if( new_box->ymax > merge_box->ymax) merge_box->ymax = new_box->ymax;
if( FLAGS_GET_Z(merge_box->flags) || FLAGS_GET_GEODETIC(merge_box->flags) )
{
- if( new_box.zmin < merge_box->zmin) merge_box->zmin = new_box.zmin;
- if( new_box.zmax > merge_box->zmax) merge_box->zmax = new_box.zmax;
+ if( new_box->zmin < merge_box->zmin) merge_box->zmin = new_box->zmin;
+ if( new_box->zmax > merge_box->zmax) merge_box->zmax = new_box->zmax;
}
if( FLAGS_GET_M(merge_box->flags) )
{
- if( new_box.mmin < merge_box->mmin) merge_box->mmin = new_box.mmin;
- if( new_box.mmax > merge_box->mmax) merge_box->mmax = new_box.mmax;
+ if( new_box->mmin < merge_box->mmin) merge_box->mmin = new_box->mmin;
+ if( new_box->mmax > merge_box->mmax) merge_box->mmax = new_box->mmax;
}
return G_SUCCESS;
}
-int gbox_overlaps(GBOX g1, GBOX g2)
+int gbox_overlaps(const GBOX *g1, const GBOX *g2)
{
/* Make sure our boxes have the same dimensionality */
- if( ! (FLAGS_GET_Z(g1.flags) == FLAGS_GET_Z(g2.flags) &&
- FLAGS_GET_M(g1.flags) == FLAGS_GET_M(g2.flags) &&
- FLAGS_GET_GEODETIC(g1.flags) == FLAGS_GET_GEODETIC(g2.flags) ) )
+ if( ! (FLAGS_GET_Z(g1->flags) == FLAGS_GET_Z(g2->flags) &&
+ FLAGS_GET_M(g1->flags) == FLAGS_GET_M(g2->flags) &&
+ FLAGS_GET_GEODETIC(g1->flags) == FLAGS_GET_GEODETIC(g2->flags) ) )
{
lwerror("gbox_overlaps: geometries have mismatched dimensionality");
}
- if( g1.xmax < g2.xmin || g1.ymax < g2.ymin ||
- g1.xmin > g2.xmax || g1.ymin > g2.ymax )
+ if( g1->xmax < g2->xmin || g1->ymax < g2->ymin ||
+ g1->xmin > g2->xmax || g1->ymin > g2->ymax )
return LW_FALSE;
- if( FLAGS_GET_Z(g1.flags) || FLAGS_GET_GEODETIC(g1.flags) )
+ if( FLAGS_GET_Z(g1->flags) || FLAGS_GET_GEODETIC(g1->flags) )
{
- if( g1.zmax < g2.zmin || g1.zmin > g2.zmax )
+ if( g1->zmax < g2->zmin || g1->zmin > g2->zmax )
return LW_FALSE;
}
- if( FLAGS_GET_M(g1.flags) )
+ if( FLAGS_GET_M(g1->flags) )
{
- if( g1.mmax < g2.mmin || g1.mmin > g2.mmax )
+ if( g1->mmax < g2->mmin || g1->mmin > g2->mmax )
return LW_FALSE;
}
return LW_TRUE;
* Warning, this function is only good for x/y/z boxes, used
* in unit testing of geodetic box generation.
*/
-GBOX* gbox_from_string(char *str)
+GBOX* gbox_from_string(const char *str)
{
- char *ptr = str;
+ const char *ptr = str;
char *nextptr;
char *gbox_start = strstr(str, "GBOX((");
GBOX *gbox = gbox_new(gflags(0,0,1));
return gbox;
}
-char* gbox_to_string(GBOX *gbox)
+char* gbox_to_string(const GBOX *gbox)
{
static int sz = 128;
char *str = NULL;
return str;
}
-GBOX* gbox_copy(GBOX *box)
+GBOX* gbox_copy(const GBOX *box)
{
GBOX *copy = (GBOX*)lwalloc(sizeof(GBOX));
memcpy(copy, box, sizeof(GBOX));
return copy;
}
-void gbox_duplicate(GBOX original, GBOX *duplicate)
+void gbox_duplicate(const GBOX *original, GBOX *duplicate)
{
assert(duplicate);
- if( original.flags != duplicate->flags )
+ if( original->flags != duplicate->flags )
lwerror("gbox_duplicate: geometries have inconsistent dimensionality");
- duplicate->xmin = original.xmin;
- duplicate->ymin = original.ymin;
- duplicate->xmax = original.xmax;
- duplicate->ymax = original.ymax;
+ duplicate->xmin = original->xmin;
+ duplicate->ymin = original->ymin;
+ duplicate->xmax = original->xmax;
+ duplicate->ymax = original->ymax;
- if( FLAGS_GET_GEODETIC(original.flags) || FLAGS_GET_Z(original.flags) )
+ if( FLAGS_GET_GEODETIC(original->flags) || FLAGS_GET_Z(original->flags) )
{
- duplicate->zmin = original.zmin;
- duplicate->zmax = original.zmax;
+ duplicate->zmin = original->zmin;
+ duplicate->zmax = original->zmax;
}
- if( FLAGS_GET_M(original.flags) )
+ if( FLAGS_GET_M(original->flags) )
{
- duplicate->mmin = original.mmin;
- duplicate->mmax = original.mmax;
+ duplicate->mmin = original->mmin;
+ duplicate->mmax = original->mmax;
}
return;
}
return 2 * FLAGS_NDIMS(flags) * sizeof(float);
}
-int gbox_from_gserialized(GSERIALIZED *g, GBOX *gbox)
+int gbox_from_gserialized(const GSERIALIZED *g, GBOX *gbox)
{
/* Null input! */
return G_SUCCESS;
}
-int ptarray_calculate_gbox( POINTARRAY *pa, GBOX *gbox )
+int ptarray_calculate_gbox(const POINTARRAY *pa, GBOX *gbox )
{
int i;
POINT4D p;
if (lwcircle_calculate_gbox(p1, p2, p3, &tmp) == G_FAILURE)
continue;
- gbox_merge(tmp, gbox);
+ gbox_merge(&tmp, gbox);
}
return G_SUCCESS;
{
if( first )
{
- gbox_duplicate(subbox, gbox);
+ gbox_duplicate(&subbox, gbox);
first = LW_FALSE;
}
else
{
- gbox_merge(subbox, gbox);
+ gbox_merge(&subbox, gbox);
}
result = G_SUCCESS;
}
* GSERIALIZED metadata utility functions.
*/
-uint32 gserialized_get_type(GSERIALIZED *s)
+uint32 gserialized_get_type(const GSERIALIZED *s)
{
uint32 *ptr;
assert(s);
return *ptr;
}
-uint32 gserialized_get_srid(GSERIALIZED *s)
+uint32 gserialized_get_srid(const GSERIALIZED *s)
{
uint32 srid = 0;
srid = srid | (s->srid[0] << 16);
/* Private functions */
-static size_t gserialized_from_any_size(LWGEOM *geom); /* Local prototype */
+static size_t gserialized_from_any_size(const LWGEOM *geom); /* Local prototype */
-static size_t gserialized_from_lwpoint_size(LWPOINT *point)
+static size_t gserialized_from_lwpoint_size(const LWPOINT *point)
{
size_t size = 4; /* Type number. */
return size;
}
-static size_t gserialized_from_lwline_size(LWLINE *line)
+static size_t gserialized_from_lwline_size(const LWLINE *line)
{
size_t size = 4; /* Type number. */
return size;
}
-static size_t gserialized_from_lwpoly_size(LWPOLY *poly)
+static size_t gserialized_from_lwpoly_size(const LWPOLY *poly)
{
size_t size = 4; /* Type number. */
int i = 0;
return size;
}
-static size_t gserialized_from_lwcircstring_size(LWCIRCSTRING *curve)
+static size_t gserialized_from_lwcircstring_size(const LWCIRCSTRING *curve)
{
size_t size = 4; /* Type number. */
return size;
}
-static size_t gserialized_from_lwcollection_size(LWCOLLECTION *col)
+static size_t gserialized_from_lwcollection_size(const LWCOLLECTION *col)
{
size_t size = 4; /* Type number. */
int i = 0;
return size;
}
-static size_t gserialized_from_any_size(LWGEOM *geom)
+static size_t gserialized_from_any_size(const LWGEOM *geom)
{
int type = TYPE_GETTYPE(geom->type);
/* Public function */
-size_t gserialized_from_lwgeom_size(LWGEOM *geom)
+size_t gserialized_from_lwgeom_size(const LWGEOM *geom)
{
size_t size = 8; /* Header overhead. */
assert(geom);
/* Private functions */
-static size_t gserialized_from_lwgeom_any(LWGEOM *geom, uchar *buf);
+static size_t gserialized_from_lwgeom_any(const LWGEOM *geom, uchar *buf);
-static size_t gserialized_from_lwpoint(LWPOINT *point, uchar *buf)
+static size_t gserialized_from_lwpoint(const LWPOINT *point, uchar *buf)
{
uchar *loc;
int ptsize = pointArray_ptsize(point->point);
return (size_t)(loc - buf);
}
-static size_t gserialized_from_lwline(LWLINE *line, uchar *buf)
+static size_t gserialized_from_lwline(const LWLINE *line, uchar *buf)
{
uchar *loc;
int ptsize;
return (size_t)(loc - buf);
}
-static size_t gserialized_from_lwpoly(LWPOLY *poly, uchar *buf)
+static size_t gserialized_from_lwpoly(const LWPOLY *poly, uchar *buf)
{
int i;
uchar *loc;
return (size_t)(loc - buf);
}
-static size_t gserialized_from_lwcircstring(LWCIRCSTRING *curve, uchar *buf)
+static size_t gserialized_from_lwcircstring(const LWCIRCSTRING *curve, uchar *buf)
{
uchar *loc;
int ptsize;
return (size_t)(loc - buf);
}
-static size_t gserialized_from_lwcollection(LWCOLLECTION *coll, uchar *buf)
+static size_t gserialized_from_lwcollection(const LWCOLLECTION *coll, uchar *buf)
{
size_t subsize = 0;
uchar *loc;
return (size_t)(loc - buf);
}
-static size_t gserialized_from_lwgeom_any(LWGEOM *geom, uchar *buf)
+static size_t gserialized_from_lwgeom_any(const LWGEOM *geom, uchar *buf)
{
int type = 0;
return 0;
}
-static size_t gserialized_from_gbox(GBOX gbox, uchar *buf)
+static size_t gserialized_from_gbox(const GBOX *gbox, uchar *buf)
{
uchar *loc;
float f;
loc = buf;
- f = nextDown_f(gbox.xmin);
+ f = nextDown_f(gbox->xmin);
memcpy(loc, &f, sizeof(float));
loc += sizeof(float);
- f = nextUp_f(gbox.xmax);
+ f = nextUp_f(gbox->xmax);
memcpy(loc, &f, sizeof(float));
loc += sizeof(float);
- f = nextDown_f(gbox.ymin);
+ f = nextDown_f(gbox->ymin);
memcpy(loc, &f, sizeof(float));
loc += sizeof(float);
- f = nextUp_f(gbox.ymax);
+ f = nextUp_f(gbox->ymax);
memcpy(loc, &f, sizeof(float));
loc += sizeof(float);
- if( FLAGS_GET_GEODETIC(gbox.flags) )
+ if( FLAGS_GET_GEODETIC(gbox->flags) )
{
- f = nextDown_f(gbox.zmin);
+ f = nextDown_f(gbox->zmin);
memcpy(loc, &f, sizeof(float));
loc += sizeof(float);
- f = nextUp_f(gbox.zmax);
+ f = nextUp_f(gbox->zmax);
memcpy(loc, &f, sizeof(float));
loc += sizeof(float);
return return_size;
}
- if( FLAGS_GET_Z(gbox.flags) )
+ if( FLAGS_GET_Z(gbox->flags) )
{
- f = nextDown_f(gbox.zmin);
+ f = nextDown_f(gbox->zmin);
memcpy(loc, &f, sizeof(float));
loc += sizeof(float);
- f = nextUp_f(gbox.zmax);
+ f = nextUp_f(gbox->zmax);
memcpy(loc, &f, sizeof(float));
loc += sizeof(float);
}
- if( FLAGS_GET_M(gbox.flags) )
+ if( FLAGS_GET_M(gbox->flags) )
{
- f = nextDown_f(gbox.mmin);
+ f = nextDown_f(gbox->mmin);
memcpy(loc, &f, sizeof(float));
loc += sizeof(float);
- f = nextUp_f(gbox.mmax);
+ f = nextUp_f(gbox->mmax);
memcpy(loc, &f, sizeof(float));
loc += sizeof(float);
}
/* Public function */
-GSERIALIZED* gserialized_from_lwgeom(LWGEOM *geom, int is_geodetic, size_t *size)
+GSERIALIZED* gserialized_from_lwgeom(const LWGEOM *geom, int is_geodetic, size_t *size)
{
size_t expected_box_size = 0;
size_t expected_size = 0;
/* Write in the serialized form of the gbox, if necessary. */
if( FLAGS_GET_BBOX(gbox.flags) )
- ptr += gserialized_from_gbox(gbox, ptr);
+ ptr += gserialized_from_gbox(&gbox, ptr);
/* Write in the serialized form of the geometry. */
ptr += gserialized_from_lwgeom_any(geom, ptr);
}
-LWGEOM* lwgeom_from_gserialized(GSERIALIZED *g)
+LWGEOM* lwgeom_from_gserialized(const GSERIALIZED *g)
{
uchar g_flags = 0;
uchar has_srid = 0;
{
if ( first )
{
- gbox_duplicate(subbox, gbox);
+ gbox_duplicate(&subbox, gbox);
first = LW_FALSE;
}
else
{
- gbox_merge(subbox, gbox);
+ gbox_merge(&subbox, gbox);
}
result = G_SUCCESS;
}
}
}
-int gserialized_calculate_gbox_geocentric_p(GSERIALIZED *g, GBOX *g_box)
+int gserialized_calculate_gbox_geocentric_p(const GSERIALIZED *g, GBOX *g_box)
{
uchar *data_ptr = NULL;
size_t g_size = 0;
return result;
}
-GBOX* gserialized_calculate_gbox_geocentric(GSERIALIZED *g)
+GBOX* gserialized_calculate_gbox_geocentric(const GSERIALIZED *g)
{
GBOX g_box;
int result = G_SUCCESS;
}
-GSERIALIZED* gserialized_copy(GSERIALIZED *g)
+GSERIALIZED* gserialized_copy(const GSERIALIZED *g)
{
GSERIALIZED *g_out = NULL;
assert(g);
return g_out;
}
-char* gserialized_to_string(GSERIALIZED *g)
+char* gserialized_to_string(const GSERIALIZED *g)
{
LWGEOM_UNPARSER_RESULT lwg_unparser_result;
LWGEOM *lwgeom = lwgeom_from_gserialized(g);
* Type == 0 in the case of the string "GEOMETRY" or "GEOGRAPHY".
* Return G_SUCCESS for success.
*/
-int geometry_type_from_string(char *str, int *type, int *z, int *m)
+int geometry_type_from_string(const char *str, int *type, int *z, int *m)
{
char *tmpstr;
int tmpstartpos, tmpendpos;
* Extract the geometry type from the serialized form (it hides in
* the anonymous data area, so this is a handy function).
*/
-extern uint32 gserialized_get_type(GSERIALIZED *g);
+extern uint32 gserialized_get_type(const GSERIALIZED *g);
/**
* Extract the SRID from the serialized form (it is packed into
* three bytes so this is a handy function).
*/
-extern uint32 gserialized_get_srid(GSERIALIZED *g);
+extern uint32 gserialized_get_srid(const GSERIALIZED *g);
/**
* Write the SRID into the serialized form (it is packed into
* form of the geodetic coordinates. Only accepts serialized geographies
* flagged as geodetic. Caller is responsible for disposing of the GBOX.
*/
-extern GBOX* gserialized_calculate_gbox_geocentric(GSERIALIZED *g);
+extern GBOX* gserialized_calculate_gbox_geocentric(const GSERIALIZED *g);
/**
* Calculate the geocentric bounding box directly from the serialized
* form of the geodetic coordinates. Only accepts serialized geographies
* flagged as geodetic.
*/
-int gserialized_calculate_gbox_geocentric_p(GSERIALIZED *g, GBOX *g_box);
+int gserialized_calculate_gbox_geocentric_p(const GSERIALIZED *g, GBOX *g_box);
/**
* Return a WKT representation of the gserialized geometry.
* Caller is responsible for disposing of the char*.
*/
-extern char* gserialized_to_string(GSERIALIZED *g);
+extern char* gserialized_to_string(const GSERIALIZED *g);
/**
* Return a copy of the input serialized geometry.
*/
-extern GSERIALIZED* gserialized_copy(GSERIALIZED *g);
+extern GSERIALIZED* gserialized_copy(const GSERIALIZED *g);
/**
* Check that coordinates of LWGEOM are all within the geodetic range.
* A spheroid with major axis == minor axis will be treated as a sphere.
* Pass in a tolerance in spheroid units.
*/
-extern double lwgeom_distance_spheroid(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox1, GBOX gbox2, SPHEROID spheroid, double tolerance);
+extern double lwgeom_distance_spheroid(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2, const GBOX *gbox1, const GBOX *gbox2, const SPHEROID *spheroid, double tolerance);
/**
* Calculate the geodetic area of a lwgeom on the sphere. The result
* will be multiplied by the average radius of the supplied spheroid.
*/
-extern double lwgeom_area_sphere(LWGEOM *lwgeom, GBOX gbox, SPHEROID spheroid);
+extern double lwgeom_area_sphere(const LWGEOM *lwgeom, const GBOX *gbox, const SPHEROID *spheroid);
/**
* Calculate the geodetic area of a lwgeom on the spheroid. The result
* will have the squared units of the spheroid axes.
*/
-extern double lwgeom_area_spheroid(LWGEOM *lwgeom, GBOX gbox, SPHEROID spheroid);
+extern double lwgeom_area_spheroid(const LWGEOM *lwgeom, const GBOX *gbox, const SPHEROID *spheroid);
/**
* Calculate the geodetic length of a lwgeom on the unit sphere. The result
* will have to by multiplied by the real radius to get the real length.
*/
-extern double lwgeom_length_spheroid(LWGEOM *geom, SPHEROID s);
+extern double lwgeom_length_spheroid(const LWGEOM *geom, const SPHEROID *s);
/**
* Calculate covers predicate for two lwgeoms on the sphere. Currently
* only handles point-in-polygon.
*/
-extern int lwgeom_covers_lwgeom_sphere(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2, GBOX gbox1, GBOX gbox2);
+extern int lwgeom_covers_lwgeom_sphere(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2, const GBOX *gbox1, const GBOX *gbox2);
/**
* New function to read doubles directly from the double* coordinate array
/**
* Calculate geodetic (x/y/z) box and add values to gbox. Return #G_SUCCESS on success.
*/
-extern int ptarray_calculate_gbox_geodetic(POINTARRAY *pa, GBOX *gbox);
+extern int ptarray_calculate_gbox_geodetic(const POINTARRAY *pa, GBOX *gbox);
/**
* Calculate box (x/y) and add values to gbox. Return #G_SUCCESS on success.
*/
-extern int ptarray_calculate_gbox( POINTARRAY *pa, GBOX *gbox );
+extern int ptarray_calculate_gbox(const POINTARRAY *pa, GBOX *gbox );
/**
* Calculate a spherical point that falls outside the geocentric gbox
*/
-void gbox_pt_outside(GBOX gbox, POINT2D *pt_outside);
+void gbox_pt_outside(const GBOX *gbox, POINT2D *pt_outside);
/**
* Create a new gbox with the dimensionality indicated by the flags. Caller
/**
* Update the merged #GBOX to be large enough to include itself and the new box.
*/
-extern int gbox_merge(GBOX new_box, GBOX *merged_box);
+extern int gbox_merge(const GBOX *new_box, GBOX *merged_box);
/**
* Update the #GBOX to be large enough to include itself and the new point.
*/
-extern int gbox_merge_point3d(POINT3D p, GBOX *gbox);
+extern int gbox_merge_point3d(const POINT3D *p, GBOX *gbox);
/**
* Return true if the point is inside the gbox
*/
-extern int gbox_contains_point3d(GBOX gbox, POINT3D pt);
+extern int gbox_contains_point3d(const GBOX *gbox, const POINT3D *pt);
/**
* Allocate a string representation of the #GBOX, based on dimensionality of flags.
*/
-extern char* gbox_to_string(GBOX *gbox);
+extern char* gbox_to_string(const GBOX *gbox);
/**
* Return a copy of the #GBOX, based on dimensionality of flags.
*/
-extern GBOX* gbox_copy(GBOX *gbox);
+extern GBOX* gbox_copy(const GBOX *gbox);
/**
* Warning, do not use this function, it is very particular about inputs.
*/
-extern GBOX* gbox_from_string(char *str);
+extern GBOX* gbox_from_string(const char *str);
/**
* Given a serialized form, extract the box if it exists, calculate it if it does not.
*/
-extern int gbox_from_gserialized(GSERIALIZED *g, GBOX *gbox);
+extern int gbox_from_gserialized(const GSERIALIZED *g, GBOX *gbox);
/**
* Return #LW_TRUE if the #GBOX overlaps, #LW_FALSE otherwise.
*/
-extern int gbox_overlaps(GBOX g1, GBOX g2);
+extern int gbox_overlaps(const GBOX *g1, const GBOX *g2);
/**
* Copy the values of original #GBOX into duplicate.
*/
-extern void gbox_duplicate(GBOX original, GBOX *duplicate);
+extern void gbox_duplicate(const GBOX *original, GBOX *duplicate);
/**
* Return the number of bytes necessary to hold a #GBOX of this dimension in
* Utility function to get type number from string. For example, a string 'POINTZ'
* would return type of 1 and z of 1 and m of 0. Valid
*/
-extern int geometry_type_from_string(char *str, int *type, int *z, int *m);
+extern int geometry_type_from_string(const char *str, int *type, int *z, int *m);
/**
* Calculate required memory segment to contain a serialized form of the LWGEOM.
* Primarily used internally by the serialization code. Exposed to allow the cunit
* tests to exercise it.
*/
-extern size_t gserialized_from_lwgeom_size(LWGEOM *geom);
+extern size_t gserialized_from_lwgeom_size(const LWGEOM *geom);
/**
* Allocate a new #GSERIALIZED from an #LWGEOM. For all non-point types, a bounding
* will contain the size of the final output, which is useful for setting the PgSQL
* VARSIZE information.
*/
-extern GSERIALIZED* gserialized_from_lwgeom(LWGEOM *geom, int is_geodetic, size_t *size);
+extern GSERIALIZED* gserialized_from_lwgeom(const LWGEOM *geom, int is_geodetic, size_t *size);
/**
* Allocate a new #LWGEOM from a #GSERIALIZED. The resulting #LWGEOM will have coordinates
* that are double aligned and suitable for direct reading using getPoint2d_p_ro
*/
-extern LWGEOM* lwgeom_from_gserialized(GSERIALIZED *g);
+extern LWGEOM* lwgeom_from_gserialized(const GSERIALIZED *g);
/**
* Serialize/deserialize from/to #G_GEOMETRY into #GSERIALIZED
* index operations can generate a large number of box-retrieval operations
* when scanning keys.
*/
-extern int lwgeom_needs_bbox(LWGEOM *geom);
+extern int lwgeom_needs_bbox(const LWGEOM *geom);
/**
* Count the total number of vertices in any #LWGEOM.
* Return true of false depending on whether a geometry is an "empty"
* geometry (no vertices members)
*/
-extern int lwgeom_is_empty(LWGEOM *geom);
+extern int lwgeom_is_empty(const LWGEOM *geom);
/**
* Return the dimensionality (relating to point/line/poly) of an lwgeom
}
-int geographic_point_equals(GEOGRAPHIC_POINT g1, GEOGRAPHIC_POINT g2)
+int geographic_point_equals(const GEOGRAPHIC_POINT *g1, const GEOGRAPHIC_POINT *g2)
{
- return FP_EQUALS(g1.lat, g2.lat) && FP_EQUALS(g1.lon, g2.lon);
+ return FP_EQUALS(g1->lat, g2->lat) && FP_EQUALS(g1->lon, g2->lon);
}
void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g)
/**
* Convert spherical coordinates to cartesion coordinates on unit sphere
*/
-void geog2cart(GEOGRAPHIC_POINT g, POINT3D *p)
+void geog2cart(const 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);
+ 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 cart2geog(POINT3D p, GEOGRAPHIC_POINT *g)
+void cart2geog(const POINT3D *p, GEOGRAPHIC_POINT *g)
{
- g->lon = atan2(p.y, p.x);
- g->lat = asin(p.z);
+ g->lon = atan2(p->y, p->x);
+ g->lat = asin(p->z);
}
/**
* Calculate the dot product of two unit vectors
* (-1 == opposite, 0 == orthogonal, 1 == identical)
*/
-static double dot_product(POINT3D p1, POINT3D p2)
+static double dot_product(const POINT3D *p1, const POINT3D *p2)
{
- return (p1.x*p2.x) + (p1.y*p2.y) + (p1.z*p2.z);
+ return (p1->x*p2->x) + (p1->y*p2->y) + (p1->z*p2->z);
}
/**
* Calculate the cross product of two vectors
*/
-static void cross_product(POINT3D a, POINT3D b, POINT3D *n)
+static void cross_product(const POINT3D *a, const POINT3D *b, POINT3D *n)
{
- n->x = a.y * b.z - a.z * b.y;
- n->y = a.z * b.x - a.x * b.z;
- n->z = a.x * b.y - a.y * b.x;
+ n->x = a->y * b->z - a->z * b->y;
+ n->y = a->z * b->x - a->x * b->z;
+ n->z = a->x * b->y - a->y * b->x;
return;
}
/**
* Calculate the sum of two vectors
*/
-static void vector_sum(POINT3D a, POINT3D b, POINT3D *n)
+static void vector_sum(const POINT3D *a, const POINT3D *b, POINT3D *n)
{
- n->x = a.x + b.x;
- n->y = a.y + b.y;
- n->z = a.z + b.z;
+ n->x = a->x + b->x;
+ n->y = a->y + b->y;
+ n->z = a->z + b->z;
return;
}
/**
* Calculate the difference of two vectors
*/
-static void vector_difference(POINT3D a, POINT3D b, POINT3D *n)
+static void vector_difference(const POINT3D *a, const POINT3D *b, POINT3D *n)
{
- n->x = a.x - b.x;
- n->y = a.y - b.y;
- n->z = a.z - b.z;
+ n->x = a->x - b->x;
+ n->y = a->y - b->y;
+ n->z = a->z - b->z;
return;
}
return;
}
-static void unit_normal(POINT3D a, POINT3D b, POINT3D *n)
+static void unit_normal(const POINT3D *a, const POINT3D *b, POINT3D *n)
{
cross_product(a, b, n);
normalize(n);
* Computes the cross product of two vectors using their lat, lng representations.
* Good even for small distances between p and q.
*/
-void robust_cross_product(GEOGRAPHIC_POINT p, GEOGRAPHIC_POINT q, POINT3D *a)
+void robust_cross_product(const GEOGRAPHIC_POINT *p, const GEOGRAPHIC_POINT *q, POINT3D *a)
{
- double lon_qpp = (q.lon + p.lon) / -2.0;
- double lon_qmp = (q.lon - p.lon) / 2.0;
- double sin_p_lat_minus_q_lat = sin(p.lat-q.lat);
- double sin_p_lat_plus_q_lat = sin(p.lat+q.lat);
+ double lon_qpp = (q->lon + p->lon) / -2.0;
+ double lon_qmp = (q->lon - p->lon) / 2.0;
+ double sin_p_lat_minus_q_lat = sin(p->lat-q->lat);
+ double sin_p_lat_plus_q_lat = sin(p->lat+q->lat);
double sin_lon_qpp = sin(lon_qpp);
double sin_lon_qmp = sin(lon_qmp);
double cos_lon_qpp = cos(lon_qpp);
sin_p_lat_plus_q_lat * cos_lon_qpp * sin_lon_qmp;
a->y = sin_p_lat_minus_q_lat * cos_lon_qpp * cos_lon_qmp +
sin_p_lat_plus_q_lat * sin_lon_qpp * sin_lon_qmp;
- a->z = cos(p.lat) * cos(q.lat) * sin(q.lon-p.lon);
+ a->z = cos(p->lat) * cos(q->lat) * sin(q->lon-p->lon);
}
void x_to_z(POINT3D *p)
}
-int crosses_dateline(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e)
+int crosses_dateline(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e)
{
- double sign_s = signum(s.lon);
- double sign_e = signum(e.lon);
- double ss = fabs(s.lon);
- double ee = fabs(e.lon);
+ double sign_s = signum(s->lon);
+ double sign_e = signum(e->lon);
+ double ss = fabs(s->lon);
+ double ee = fabs(e->lon);
if( sign_s == sign_e )
{
return LW_FALSE;
* 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)
+int edge_point_on_plane(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p)
{
POINT3D normal, pt;
double w;
/* Normal to the plane defined by e */
- robust_cross_product(e.start, e.end, &normal);
+ robust_cross_product(&(e->start), &(e->end), &normal);
normalize(&normal);
geog2cart(p, &pt);
/* We expect the dot product of with normal with any vector in the plane to be zero */
- w = dot_product(normal, pt);
+ w = dot_product(&normal, &pt);
LWDEBUGF(4,"dot product %.9g",w);
if ( FP_IS_ZERO(w) )
{
* Returns true if the point p is inside the cone defined by the
* two ends of the edge e.
*/
-int edge_point_in_cone(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p)
+int edge_point_in_cone(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p)
{
POINT3D vcp, vs, ve, vp;
double vs_dot_vcp, vp_dot_vcp;
- geog2cart(e.start, &vs);
- geog2cart(e.end, &ve);
+ geog2cart(&(e->start), &vs);
+ geog2cart(&(e->end), &ve);
/* Antipodal case, everything is inside. */
if( vs.x == -1.0 * ve.x && vs.y == -1.0 * ve.y && vs.z == -1.0 * ve.z )
return LW_TRUE;
geog2cart(p, &vp);
/* The normalized sum bisects the angle between start and end. */
- vector_sum(vs, ve, &vcp);
+ vector_sum(&vs, &ve, &vcp);
normalize(&vcp);
/* The projection of start onto the center defines the minimum similarity */
- vs_dot_vcp = dot_product(vs, vcp);
+ vs_dot_vcp = dot_product(&vs, &vcp);
LWDEBUGF(4,"vs_dot_vcp %.19g",vs_dot_vcp);
/* The projection of candidate p onto the center */
- vp_dot_vcp = dot_product(vp, vcp);
+ vp_dot_vcp = dot_product(&vp, &vcp);
LWDEBUGF(4,"vp_dot_vcp %.19g",vp_dot_vcp);
/* If p is more similar than start then p is inside the cone */
LWDEBUGF(4,"fabs(vp_dot_vcp - vs_dot_vcp) %.39g",fabs(vp_dot_vcp - vs_dot_vcp));
/**
* True if the longitude of p is within the range of the longitude of the ends of e
*/
-int edge_contains_coplanar_point(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p)
+int edge_contains_coplanar_point(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p)
{
GEOGRAPHIC_EDGE g;
GEOGRAPHIC_POINT q;
- double slon = fabs(e.start.lon) + fabs(e.end.lon);
- double dlon = fabs(fabs(e.start.lon) - fabs(e.end.lon));
- double slat = e.start.lat + e.end.lat;
+ double slon = fabs((e->start).lon) + fabs((e->end).lon);
+ double dlon = fabs(fabs((e->start).lon) - fabs((e->end).lon));
+ double slat = (e->start).lat + (e->end).lat;
- LWDEBUGF(4, "e.start == GPOINT(%.6g %.6g) ", e.start.lat, e.start.lon);
- LWDEBUGF(4, "e.end == GPOINT(%.6g %.6g) ", e.end.lat, e.end.lon);
- LWDEBUGF(4, "p == GPOINT(%.6g %.6g) ", p.lat, p.lon);
+ LWDEBUGF(4, "e.start == GPOINT(%.6g %.6g) ", (e->start).lat, (e->start).lon);
+ LWDEBUGF(4, "e.end == GPOINT(%.6g %.6g) ", (e->end).lat, (e->end).lon);
+ LWDEBUGF(4, "p == GPOINT(%.6g %.6g) ", p->lat, p->lon);
/* Copy values into working registers */
- g = e;
- q = p;
+ g = *e;
+ q = *p;
/* Vertical plane, we need to do this calculation in latitude */
if( FP_EQUALS( g.start.lon, g.end.lon ) )
/**
* Given two points on a unit sphere, calculate their distance apart in radians.
*/
-double sphere_distance(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e)
+double sphere_distance(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e)
{
- double d_lon = e.lon - s.lon;
+ double d_lon = e->lon - s->lon;
double cos_d_lon = cos(d_lon);
- double cos_lat_e = cos(e.lat);
- double sin_lat_e = sin(e.lat);
- double cos_lat_s = cos(s.lat);
- double sin_lat_s = sin(s.lat);
+ double cos_lat_e = cos(e->lat);
+ double sin_lat_e = sin(e->lat);
+ double cos_lat_s = cos(s->lat);
+ double sin_lat_s = sin(s->lat);
double a1 = pow(cos_lat_e * sin(d_lon), 2.0);
double a2 = pow(cos_lat_s * sin_lat_e - sin_lat_s * cos_lat_e * cos_d_lon, 2.0);
/**
* Given two unit vectors, calculate their distance apart in radians.
*/
-double sphere_distance_cartesian(POINT3D s, POINT3D e)
+double sphere_distance_cartesian(const POINT3D *s, const POINT3D *e)
{
return acos(dot_product(s, e));
}
/**
* Given two points on a unit sphere, calculate the direction from s to e.
*/
-static double sphere_direction(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e, double d)
+static double sphere_direction(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e, double d)
{
double heading = 0.0;
- if( FP_IS_ZERO(cos(s.lat)) )
- return (s.lat > 0.0) ? M_PI : 0.0;
+ if( FP_IS_ZERO(cos(s->lat)) )
+ return (s->lat > 0.0) ? M_PI : 0.0;
- if( sin(e.lon - s.lon) < 0.0 )
+ if( sin(e->lon - s->lon) < 0.0 )
{
- heading = acos((sin(e.lat) - sin(s.lat) * cos(d)) / (sin(d) * cos(s.lat)));
+ heading = acos((sin(e->lat) - sin(s->lat) * cos(d)) / (sin(d) * cos(s->lat)));
}
else
{
- heading = -1.0 * acos((sin(e.lat) - sin(s.lat) * cos(d)) / (sin(d) * cos(s.lat)));
+ heading = -1.0 * acos((sin(e->lat) - sin(s->lat) * cos(d)) / (sin(d) * cos(s->lat)));
}
return heading;
}
* @param c The last triangle vertex.
* @return the signed spherical excess.
*/
-static double sphere_excess(GEOGRAPHIC_POINT a, GEOGRAPHIC_POINT b, GEOGRAPHIC_POINT c)
+static double sphere_excess(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const GEOGRAPHIC_POINT *c)
{
double a_dist = sphere_distance(b, c);
double b_dist = sphere_distance(c, a);
* 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 edge_contains_point(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p)
{
if ( edge_point_in_cone(e, p) && edge_point_on_plane(e, p) )
/* if ( edge_contains_coplanar_point(e, p) && edge_point_on_plane(e, p) ) */
* the great circle with the line of maximum/minimum gradiant that lies on
* the great circle plane.
*/
-int clairaut_cartesian(POINT3D start, POINT3D end, GEOGRAPHIC_POINT *g_top, GEOGRAPHIC_POINT *g_bottom)
+int clairaut_cartesian(const POINT3D *start, const POINT3D *end, GEOGRAPHIC_POINT *g_top, GEOGRAPHIC_POINT *g_bottom)
{
POINT3D t1, t2;
GEOGRAPHIC_POINT vN1, vN2;
unit_normal(end, start, &t2);
LWDEBUGF(4, "unit normal t1 == POINT(%.8g %.8g %.8g)", t1.x, t1.y, t1.z);
LWDEBUGF(4, "unit normal t2 == POINT(%.8g %.8g %.8g)", t2.x, t2.y, t2.z);
- cart2geog(t1, &vN1);
- cart2geog(t2, &vN2);
+ cart2geog(&t1, &vN1);
+ cart2geog(&t2, &vN2);
g_top->lat = z_to_latitude(t1.z,LW_TRUE);
g_top->lon = vN2.lon;
g_bottom->lat = z_to_latitude(t2.z,LW_FALSE);
* 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, GEOGRAPHIC_POINT *g_top, GEOGRAPHIC_POINT *g_bottom)
+int clairaut_geographic(const GEOGRAPHIC_POINT *start, const GEOGRAPHIC_POINT *end, GEOGRAPHIC_POINT *g_top, GEOGRAPHIC_POINT *g_bottom)
{
POINT3D t1, t2;
GEOGRAPHIC_POINT vN1, vN2;
normalize(&t2);
LWDEBUGF(4, "unit normal t1 == POINT(%.8g %.8g %.8g)", t1.x, t1.y, t1.z);
LWDEBUGF(4, "unit normal t2 == POINT(%.8g %.8g %.8g)", t2.x, t2.y, t2.z);
- cart2geog(t1, &vN1);
- cart2geog(t2, &vN2);
+ cart2geog(&t1, &vN1);
+ cart2geog(&t2, &vN2);
g_top->lat = z_to_latitude(t1.z,LW_TRUE);
g_top->lon = vN2.lon;
g_bottom->lat = z_to_latitude(t2.z,LW_FALSE);
* Returns true if an intersection can be calculated, and places it in *g.
* Returns false otherwise.
*/
-int edge_intersection(GEOGRAPHIC_EDGE e1, GEOGRAPHIC_EDGE e2, GEOGRAPHIC_POINT *g)
+int edge_intersection(const GEOGRAPHIC_EDGE *e1, const GEOGRAPHIC_EDGE *e2, GEOGRAPHIC_POINT *g)
{
POINT3D ea, eb, v;
- LWDEBUGF(4, "e1 start(%.20g %.20g) end(%.20g %.20g)", e1.start.lat, e1.start.lon, e1.end.lat, e1.end.lon);
- LWDEBUGF(4, "e2 start(%.20g %.20g) end(%.20g %.20g)", e2.start.lat, e2.start.lon, e2.end.lat, e2.end.lon);
+ LWDEBUGF(4, "e1 start(%.20g %.20g) end(%.20g %.20g)", e1->start.lat, e1->start.lon, e1->end.lat, e1->end.lon);
+ LWDEBUGF(4, "e2 start(%.20g %.20g) end(%.20g %.20g)", e2->start.lat, e2->start.lon, e2->end.lat, e2->end.lon);
- LWDEBUGF(4, "e1 start(%.20g %.20g) end(%.20g %.20g)", rad2deg(e1.start.lon), rad2deg(e1.start.lat), rad2deg(e1.end.lon), rad2deg(e1.end.lat));
- LWDEBUGF(4, "e2 start(%.20g %.20g) end(%.20g %.20g)", rad2deg(e2.start.lon), rad2deg(e2.start.lat), rad2deg(e2.end.lon), rad2deg(e2.end.lat));
+ LWDEBUGF(4, "e1 start(%.20g %.20g) end(%.20g %.20g)", rad2deg(e1->start.lon), rad2deg(e1->start.lat), rad2deg(e1->end.lon), rad2deg(e1->end.lat));
+ LWDEBUGF(4, "e2 start(%.20g %.20g) end(%.20g %.20g)", rad2deg(e2->start.lon), rad2deg(e2->start.lat), rad2deg(e2->end.lon), rad2deg(e2->end.lat));
- if( geographic_point_equals(e1.start, e2.start) )
+ if( geographic_point_equals(&(e1->start), &(e2->start)) )
{
- *g = e1.start;
+ *g = e1->start;
return LW_TRUE;
}
- if( geographic_point_equals(e1.end, e2.end) )
+ if( geographic_point_equals(&(e1->end), &(e2->end)) )
{
- *g = e1.end;
+ *g = e1->end;
return LW_TRUE;
}
- if( geographic_point_equals(e1.end, e2.start) )
+ if( geographic_point_equals(&(e1->end), &(e2->start)) )
{
- *g = e1.end;
+ *g = e1->end;
return LW_TRUE;
}
- if( geographic_point_equals(e1.start, e2.end) )
+ if( geographic_point_equals(&(e1->start), &(e2->end)) )
{
- *g = e1.start;
+ *g = e1->start;
return LW_TRUE;
}
- robust_cross_product(e1.start, e1.end, &ea);
+ robust_cross_product(&(e1->start), &(e1->end), &ea);
normalize(&ea);
- robust_cross_product(e2.start, e2.end, &eb);
+ robust_cross_product(&(e2->start), &(e2->end), &eb);
normalize(&eb);
LWDEBUGF(4, "e1 cross product == POINT(%.12g %.12g %.12g)", ea.x, ea.y, ea.z);
LWDEBUGF(4, "e2 cross product == POINT(%.12g %.12g %.12g)", eb.x, eb.y, eb.z);
- LWDEBUGF(4, "fabs(dot_product(ea, eb)) == %.14g", fabs(dot_product(ea, eb)));
- if( FP_EQUALS(fabs(dot_product(ea, eb)), 1.0) )
+ LWDEBUGF(4, "fabs(dot_product(ea, eb)) == %.14g", fabs(dot_product(&ea, &eb)));
+ if( FP_EQUALS(fabs(dot_product(&ea, &eb)), 1.0) )
{
- LWDEBUGF(4, "parallel edges found! dot_product = %.12g", dot_product(ea, eb));
+ LWDEBUGF(4, "parallel edges found! dot_product = %.12g", dot_product(&ea, &eb));
/* Parallel (maybe equal) edges! */
/* Hack alert, only returning ONE end of the edge right now, most do better later. */
/* Hack alart #2, returning a value of 2 to indicate a co-linear crossing event. */
- if ( edge_contains_point(e1, e2.start) )
+ if ( edge_contains_point(e1, &(e2->start)) )
{
- *g = e2.start;
+ *g = e2->start;
return 2;
}
- if ( edge_contains_point(e1, e2.end) )
+ if ( edge_contains_point(e1, &(e2->end)) )
{
- *g = e2.end;
+ *g = e2->end;
return 2;
}
- if ( edge_contains_point(e2, e1.start) )
+ if ( edge_contains_point(e2, &(e1->start)) )
{
- *g = e1.start;
+ *g = e1->start;
return 2;
}
- if ( edge_contains_point(e2, e1.end) )
+ if ( edge_contains_point(e2, &(e1->end)) )
{
- *g = e1.end;
+ *g = e1->end;
return 2;
}
}
- unit_normal(ea, eb, &v);
+ unit_normal(&ea, &eb, &v);
LWDEBUGF(4, "v == POINT(%.12g %.12g %.12g)", v.x, v.y, v.z);
g->lat = atan2(v.z, sqrt(v.x * v.x + v.y * v.y));
g->lon = atan2(v.y, v.x);
LWDEBUGF(4, "g == GPOINT(%.12g %.12g)", g->lat, g->lon);
LWDEBUGF(4, "g == POINT(%.12g %.12g)", rad2deg(g->lon), rad2deg(g->lat));
- if ( edge_contains_point(e1, *g) && edge_contains_point(e2, *g) )
+ if ( edge_contains_point(e1, g) && edge_contains_point(e2, g) )
{
return LW_TRUE;
}
{
g->lon = -1.0 * (2.0 * M_PI - g->lon);
}
- if ( edge_contains_point(e1, *g) && edge_contains_point(e2, *g) )
+ if ( edge_contains_point(e1, g) && edge_contains_point(e2, g) )
{
return LW_TRUE;
}
return LW_FALSE;
}
-double edge_distance_to_point(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT gp, GEOGRAPHIC_POINT *closest)
+double edge_distance_to_point(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *gp, GEOGRAPHIC_POINT *closest)
{
double d1 = 1000000000.0, d2, d3, d_nearest;
POINT3D n, p, k;
GEOGRAPHIC_POINT gk, g_nearest;
/* Zero length edge, */
- if( geographic_point_equals(e.start,e.end) )
- return sphere_distance(e.start, gp);
+ if( geographic_point_equals(&(e->start), &(e->end)) )
+ return sphere_distance(&(e->start), gp);
- robust_cross_product(e.start, e.end, &n);
+ robust_cross_product(&(e->start), &(e->end), &n);
normalize(&n);
geog2cart(gp, &p);
- vector_scale(&n, dot_product(p, n));
- vector_difference(p, n, &k);
+ vector_scale(&n, dot_product(&p, &n));
+ vector_difference(&p, &n, &k);
normalize(&k);
- cart2geog(k, &gk);
- if( edge_contains_point(e, gk) )
+ cart2geog(&k, &gk);
+ if( edge_contains_point(e, &gk) )
{
- d1 = sphere_distance(gp, gk);
+ d1 = sphere_distance(gp, &gk);
}
- d2 = sphere_distance(gp, e.start);
- d3 = sphere_distance(gp, e.end);
+ d2 = sphere_distance(gp, &(e->start));
+ d3 = sphere_distance(gp, &(e->end));
d_nearest = d1;
g_nearest = gk;
if( d2 < d_nearest )
{
d_nearest = d2;
- g_nearest = e.start;
+ g_nearest = e->start;
}
if( d3 < d_nearest )
{
d_nearest = d3;
- g_nearest = e.end;
+ g_nearest = e->end;
}
if(closest)
*closest = g_nearest;
return d_nearest;
}
-double edge_distance_to_edge(GEOGRAPHIC_EDGE e1, GEOGRAPHIC_EDGE e2, GEOGRAPHIC_POINT *closest1, GEOGRAPHIC_POINT *closest2)
+double edge_distance_to_edge(const GEOGRAPHIC_EDGE *e1, const GEOGRAPHIC_EDGE *e2, GEOGRAPHIC_POINT *closest1, GEOGRAPHIC_POINT *closest2)
{
double d;
GEOGRAPHIC_POINT gcp1s, gcp1e, gcp2s, gcp2e, c1, c2;
- double d1s = edge_distance_to_point(e1, e2.start, &gcp1s);
- double d1e = edge_distance_to_point(e1, e2.end, &gcp1e);
- double d2s = edge_distance_to_point(e2, e1.start, &gcp2s);
- double d2e = edge_distance_to_point(e2, e1.end, &gcp2e);
+ double d1s = edge_distance_to_point(e1, &(e2->start), &gcp1s);
+ double d1e = edge_distance_to_point(e1, &(e2->end), &gcp1e);
+ double d2s = edge_distance_to_point(e2, &(e1->start), &gcp2s);
+ double d2e = edge_distance_to_point(e2, &(e1->end), &gcp2e);
d = d1s;
c1 = gcp1s;
- c2 = e2.start;
+ c2 = e2->start;
if( d1e < d )
{
d = d1e;
c1 = gcp1e;
- c2 = e2.end;
+ c2 = e2->end;
}
if( d2s < d )
{
d = d2s;
- c1 = e1.start;
+ c1 = e1->start;
c2 = gcp2s;
}
if( d2e < d )
{
d = d2e;
- c1 = e1.end;
+ c1 = e1->end;
c2 = gcp2e;
}
* 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.
*/
-int sphere_project(GEOGRAPHIC_POINT r, double distance, double azimuth, GEOGRAPHIC_POINT *n)
+int sphere_project(const GEOGRAPHIC_POINT *r, double distance, double azimuth, GEOGRAPHIC_POINT *n)
{
double d = distance;
- double lat1 = r.lat;
+ double lat1 = r->lat;
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;
+ n->lon = atan(b/a) + r->lon;
return G_SUCCESS;
}
-int edge_calculate_gbox_slow(GEOGRAPHIC_EDGE e, GBOX *gbox)
+int edge_calculate_gbox_slow(const GEOGRAPHIC_EDGE *e, GBOX *gbox)
{
int steps = 1000000;
int i;
double dx, dy, dz;
- double distance = sphere_distance(e.start, e.end);
+ double distance = sphere_distance(&(e->start), &(e->end));
POINT3D pn, p, start, end;
/* Edge is zero length, just return the naive box */
if ( FP_IS_ZERO(distance) )
{
LWDEBUG(4, "edge is zero length. returning");
- geog2cart(e.start, &start);
- geog2cart(e.end, &end);
+ geog2cart(&(e->start), &start);
+ geog2cart(&(e->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);
/* Walk along the chord between start and end incrementally,
normalizing at each step. */
- geog2cart(e.start, &start);
- geog2cart(e.end, &end);
+ geog2cart(&(e->start), &start);
+ geog2cart(&(e->end), &end);
dx = (end.x - start.x)/steps;
dy = (end.y - start.y)/steps;
dz = (end.z - start.z)/steps;
p.z += dz;
pn = p;
normalize(&pn);
- gbox_merge_point3d(pn, gbox);
+ gbox_merge_point3d(&pn, gbox);
}
return G_SUCCESS;
}
* for edges over poles and fully around the globe. An edge is assumed
* to follow the shortest great circle route between the end points.
*/
-int edge_calculate_gbox(GEOGRAPHIC_EDGE e, GBOX *gbox)
+int edge_calculate_gbox(const GEOGRAPHIC_EDGE *e, GBOX *gbox)
{
double deltaLongitude;
- double distance = sphere_distance(e.start, e.end);
+ 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, nT1, nT2;
}
/* Initialize our working copy of the edge */
- g = e;
+ g = *e;
LWDEBUG(4, "entered function");
LWDEBUGF(4, "edge length: %.8g", distance);
if ( FP_IS_ZERO(distance) )
{
LWDEBUG(4, "edge is zero length. returning");
- geog2cart(g.start, &start);
- geog2cart(g.end, &end);
+ 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);
LWDEBUGF(4, "edge values: (%.6g %.6g, %.6g %.6g)", g.start.lon, g.start.lat, g.end.lon, g.end.lat);
/* Initialize box with the start and end points of the edge. */
- geog2cart(g.start, &start);
- geog2cart(g.end, &end);
+ 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);
else
{
LWDEBUG(4, "not a pole crossing, calculating clairaut points");
- clairaut_cartesian(start, end, &vT1, &vT2);
+ clairaut_cartesian(&start, &end, &vT1, &vT2);
LWDEBUGF(4, "vT1 == GPOINT(%.6g %.6g) ", vT1.lat, vT1.lon);
LWDEBUGF(4, "vT2 == GPOINT(%.6g %.6g) ", vT2.lat, vT2.lon);
- if ( edge_contains_point(g, vT1) )
+ if ( edge_contains_point(&g, &vT1) )
{
- geog2cart(vT1, &p);
+ geog2cart(&vT1, &p);
LWDEBUGF(4, "p == POINT(%.8g %.8g %.8g)", p.x, p.y, p.z);
- gbox_merge_point3d(p, gbox);
+ gbox_merge_point3d(&p, gbox);
LWDEBUG(4, "edge contained vT1");
LWDEBUGF(4, "gbox: %s", gbox_to_string(gbox));
}
- else if ( edge_contains_point(g, vT2) )
+ else if ( edge_contains_point(&g, &vT2) )
{
- geog2cart(vT2, &p);
+ geog2cart(&vT2, &p);
LWDEBUGF(4, "p == POINT(%.8g %.8g %.8g)", p.x, p.y, p.z);
- gbox_merge_point3d(p, gbox);
+ gbox_merge_point3d(&p, gbox);
LWDEBUG(4, "edge contained vT2");
LWDEBUGF(4, "gbox: %s", gbox_to_string(gbox));
}
endXZ = end;
x_to_z(&startXZ);
x_to_z(&endXZ);
- clairaut_cartesian(startXZ, endXZ, &vT1, &vT2);
+ clairaut_cartesian(&startXZ, &endXZ, &vT1, &vT2);
gimbal_lock = LW_FALSE;
LWDEBUG(4, "vT1/vT2 before flipping back z to x");
LWDEBUGF(4, "vT1 == GPOINT(%.6g %.6g) ", vT1.lat, vT1.lon);
{
gimbal_lock = LW_TRUE;
}
- geog2cart(vT1, &nT1);
- geog2cart(vT2, &nT2);
+ geog2cart(&vT1, &nT1);
+ geog2cart(&vT2, &nT2);
x_to_z(&nT1);
x_to_z(&nT2);
- cart2geog(nT1, &vT1);
- cart2geog(nT2, &vT2);
+ cart2geog(&nT1, &vT1);
+ cart2geog(&nT2, &vT2);
LWDEBUG(4, "vT1/vT2 after flipping back z to x");
LWDEBUGF(4, "vT1 == GPOINT(%.6g %.6g) ", vT1.lat, vT1.lon);
LWDEBUGF(4, "vT2 == GPOINT(%.6g %.6g) ", vT2.lat, vT2.lon);
LWDEBUGF(4, "p1 == POINT(%.8g %.8g %.8g)", nT1.x, nT1.y, nT1.z);
LWDEBUGF(4, "p2 == POINT(%.8g %.8g %.8g)", nT2.x, nT2.y, nT2.z);
*/
- if ( edge_contains_point(g, vT1) )
+ if ( edge_contains_point(&g, &vT1) )
{
- geog2cart(vT1, &p);
+ geog2cart(&vT1, &p);
LWDEBUGF(4, "p == POINT(%.8g %.8g %.8g)", p.x, p.y, p.z);
- gbox_merge_point3d(p, gbox);
+ gbox_merge_point3d(&p, gbox);
LWDEBUG(4, "edge contained vT1");
LWDEBUGF(4, "gbox: %s", gbox_to_string(gbox));
}
- else if ( edge_contains_point(g, vT2) )
+ else if ( edge_contains_point(&g, &vT2) )
{
- geog2cart(vT2, &p);
+ geog2cart(&vT2, &p);
LWDEBUGF(4, "p == POINT(%.8g %.8g %.8g)", p.x, p.y, p.z);
- gbox_merge_point3d(p, gbox);
+ gbox_merge_point3d(&p, gbox);
LWDEBUG(4, "edge contained vT2");
LWDEBUGF(4, "gbox: %s", gbox_to_string(gbox));
}
endYZ = end;
y_to_z(&startYZ);
y_to_z(&endYZ);
- clairaut_cartesian(startYZ, endYZ, &vT1, &vT2);
+ clairaut_cartesian(&startYZ, &endYZ, &vT1, &vT2);
gimbal_lock = LW_FALSE;
LWDEBUG(4, "vT1/vT2 before flipping back z to y");
LWDEBUGF(4, "vT1 == GPOINT(%.6g %.6g) ", vT1.lat, vT1.lon);
{
gimbal_lock = LW_TRUE;
}
- geog2cart(vT1, &nT1);
- geog2cart(vT2, &nT2);
+ geog2cart(&vT1, &nT1);
+ geog2cart(&vT2, &nT2);
y_to_z(&nT1);
y_to_z(&nT2);
- cart2geog(nT1, &vT1);
- cart2geog(nT2, &vT2);
+ cart2geog(&nT1, &vT1);
+ cart2geog(&nT2, &vT2);
LWDEBUG(4, "vT1/vT2 after flipping back z to y");
LWDEBUGF(4, "vT1 == GPOINT(%.6g %.6g) ", vT1.lat, vT1.lon);
LWDEBUGF(4, "vT2 == GPOINT(%.6g %.6g) ", vT2.lat, vT2.lon);
LWDEBUGF(4, "p1 == POINT(%.8g %.8g %.8g)", nT1.x, nT1.y, nT1.z);
LWDEBUGF(4, "p2 == POINT(%.8g %.8g %.8g)", nT2.x, nT2.y, nT2.z);
*/
- if ( edge_contains_point(g, vT1) )
+ if ( edge_contains_point(&g, &vT1) )
{
- geog2cart(vT1, &p);
+ geog2cart(&vT1, &p);
LWDEBUGF(4, "p == POINT(%.8g %.8g %.8g)", p.x, p.y, p.z);
- gbox_merge_point3d(p, gbox);
+ gbox_merge_point3d(&p, gbox);
LWDEBUG(4, "edge contained vT1");
LWDEBUGF(4, "gbox: %s", gbox_to_string(gbox));
}
- else if ( edge_contains_point(g, vT2) )
+ else if ( edge_contains_point(&g, &vT2) )
{
- geog2cart(vT2, &p);
+ geog2cart(&vT2, &p);
LWDEBUGF(4, "p == POINT(%.8g %.8g %.8g)", p.x, p.y, p.z);
- gbox_merge_point3d(p, gbox);
+ gbox_merge_point3d(&p, gbox);
LWDEBUG(4, "edge contained vT2");
LWDEBUGF(4, "gbox: %s", gbox_to_string(gbox));
}
* Given a unit geocentric gbox, return a lon/lat (degrees) coordinate point point that is
* guaranteed to be outside the box (and therefore anything it contains).
*/
-void gbox_pt_outside(GBOX gbox, POINT2D *pt_outside)
+void gbox_pt_outside(const GBOX *gbox, POINT2D *pt_outside)
{
static double grow = M_PI / 180.0 / 60.0; /* one arc-minute */
double d;
GEOGRAPHIC_POINT g;
/* Assign our box and expand it slightly. */
- ge = gbox;
+ ge = *gbox;
ge.xmin -= grow;
ge.ymin -= grow;
ge.zmin -= grow;
for( i = 0; i < 8; i++ )
{
normalize(&(corners[i]));
- if( ! gbox_contains_point3d(gbox, corners[i]) )
+ if( ! gbox_contains_point3d(gbox, &(corners[i])) )
{
pt = corners[i];
normalize(&pt);
- cart2geog(pt, &g);
+ cart2geog(&pt, &g);
pt_outside->x = rad2deg(g.lon);
pt_outside->y = rad2deg(g.lat);
return;
pt.y = 0.0;
pt.z = 0.0;
- if((1.0 - gbox.xmax) > 0.1)
+ if((1.0 - gbox->xmax) > 0.1)
{
- pt.x = gbox.xmax + (1.0 - gbox.xmax) * 0.01;
+ pt.x = gbox->xmax + (1.0 - gbox->xmax) * 0.01;
d = sqrt((1.0 - pow(pt.x, 2.0))/2.0);
pt.y = d;
pt.z = d;
}
- else if((1.0 - gbox.ymax) > 0.1)
+ else if((1.0 - gbox->ymax) > 0.1)
{
- pt.y = gbox.ymax + (1.0 - gbox.ymax) * 0.01;
+ pt.y = gbox->ymax + (1.0 - gbox->ymax) * 0.01;
d = sqrt((1.0 - pow(pt.y, 2.0))/2.0);
pt.x = d;
pt.z = d;
}
- else if((1.0 - gbox.zmax) > 0.1)
+ else if((1.0 - gbox->zmax) > 0.1)
{
- pt.z = gbox.zmax + (1.0 - gbox.zmax) * 0.01;
+ pt.z = gbox->zmax + (1.0 - gbox->zmax) * 0.01;
d = sqrt((1.0 - pow(pt.z, 2.0))/2.0);
pt.x = d;
pt.y = d;
}
normalize(&pt);
- cart2geog(pt, &g);
+ cart2geog(&pt, &g);
pt_outside->x = rad2deg(g.lon);
pt_outside->y = rad2deg(g.lat);
return;
* Returns the area of the ring (ring must be closed) in square radians (surface of
* the sphere is 4*PI).
*/
-double ptarray_area_sphere(POINTARRAY *pa, POINT2D pt_outside)
+double ptarray_area_sphere(const POINTARRAY *pa, const POINT2D *pt_outside)
{
GEOGRAPHIC_POINT a, b, c;
POINT2D p;
if( ! pa || pa->npoints < 4 )
return 0.0;
- geographic_point_init(pt_outside.x, pt_outside.y, &c);
+ geographic_point_init(pt_outside->x, pt_outside->y, &c);
/* Initialize first point */
getPoint2d_p(pa, 0, &p);
getPoint2d_p(pa, i, &p);
geographic_point_init(p.x, p.y, &b);
- if( crosses_dateline(a, b) )
+ if( crosses_dateline(&a, &b) )
{
GEOGRAPHIC_POINT a1 = a, b1 = b, c1 = c;
double shift;
point_shift(&a1, shift);
point_shift(&b1, shift);
point_shift(&c1, shift);
- excess = sphere_excess(a1, b1, c1);
+ excess = sphere_excess(&a1, &b1, &c1);
}
else
{
- excess = sphere_excess(a, b, c);
+ excess = sphere_excess(&a, &b, &c);
}
area += excess;
* to derive one in postgis, or the gbox_pt_outside() function if you don't mind burning CPU cycles
* building a gbox first).
*/
-int ptarray_point_in_ring(POINTARRAY *pa, POINT2D pt_outside, POINT2D pt_to_test)
+int ptarray_point_in_ring(const POINTARRAY *pa, const POINT2D *pt_outside, const POINT2D *pt_to_test)
{
GEOGRAPHIC_EDGE crossing_edge, edge;
POINT2D p;
return LW_FALSE;
/* Set up our stab line */
- geographic_point_init(pt_to_test.x, pt_to_test.y, &(crossing_edge.start));
- geographic_point_init(pt_outside.x, pt_outside.y, &(crossing_edge.end));
+ geographic_point_init(pt_to_test->x, pt_to_test->y, &(crossing_edge.start));
+ geographic_point_init(pt_outside->x, pt_outside->y, &(crossing_edge.end));
/* Initialize first point */
getPoint2d_p(pa, first_point, &p);
geographic_point_init(p.x, p.y, &(edge.start));
/* If the start point is on the stab line, back up until it isn't */
- while(edge_contains_point(crossing_edge, edge.start) && ! geographic_point_equals(crossing_edge.start, edge.start) )
+ while(edge_contains_point(&crossing_edge, &(edge.start)) && ! geographic_point_equals(&(crossing_edge.start), &(edge.start)) )
{
first_point--;
LWDEBUGF(4,"first point was on stab line, reversing %d points", first_point);
LWDEBUGF(4,"testing edge (%d)", i);
/* Our test point is on an edge! Point is "in ring" by our definition */
- if( geographic_point_equals(crossing_edge.start, edge.start) ||
- geographic_point_equals(crossing_edge.start, edge.end) ||
- edge_contains_point(edge, crossing_edge.start) )
+ if( geographic_point_equals(&(crossing_edge.start), &(edge.start)) ||
+ geographic_point_equals(&(crossing_edge.start), &(edge.end)) ||
+ edge_contains_point(&edge, &(crossing_edge.start)) )
{
LWDEBUGF(4,"edge (%d) contains the test point, returning true", i);
return LW_TRUE;
/* If the end of our edge is on the stab line, extend the edge to the
next end point, by skipping the start->end assignment step at the
end of this loop */
- if(edge_contains_point(crossing_edge, edge.end))
+ if(edge_contains_point(&crossing_edge, &(edge.end)))
{
LWDEBUGF(4,"edge (%d) end point is on the stab line, continuing", i);
continue;
LWDEBUG(4,"testing edge crossing");
- if( edge_intersection(edge, crossing_edge, &g) )
+ if( edge_intersection(&edge, &crossing_edge, &g) )
{
count++;
LWDEBUGF(4,"edge (%d) crossed, count == %d", i, count);
}
-static double ptarray_distance_spheroid(POINTARRAY *pa1, POINTARRAY *pa2, SPHEROID s, double tolerance, int check_intersection)
+static double ptarray_distance_spheroid(const POINTARRAY *pa1, const POINTARRAY *pa2, const SPHEROID *s, double tolerance, int check_intersection)
{
GEOGRAPHIC_EDGE e1, e2;
GEOGRAPHIC_POINT g1, g2;
POINT2D p;
double distance;
int i, j;
+ int use_sphere = (s->a == s->b ? 1 : 0);
/* Make result really big, so that everything will be smaller than it */
distance = MAXFLOAT;
getPoint2d_p(pa2, 0, &p);
geographic_point_init(p.x, p.y, &g2);
/* Sphere special case, axes equal */
- if( s.a == s.b )
- distance = s.radius * sphere_distance(g1, g2);
+ if( use_sphere )
+ distance = s->radius * sphere_distance(&g1, &g2);
else
- distance = spheroid_distance(g1, g2, s);
+ distance = spheroid_distance(&g1, &g2, s);
return distance;
}
{
/* Handle one/many case here */
int i;
- POINTARRAY *pa_one, *pa_many;
+ const POINTARRAY *pa_one;
+ const POINTARRAY *pa_many;
if( pa1->npoints == 1 )
{
getPoint2d_p(pa_many, i, &p);
geographic_point_init(p.x, p.y, &(e1.end));
/* Get the spherical distance between point and edge */
- d = s.radius * edge_distance_to_point(e1, g1, &g2);
+ d = s->radius * edge_distance_to_point(&e1, &g1, &g2);
/* New shortest distance! Record this distance / location */
if( d < distance )
{
if( d < tolerance )
{
/* Working on a sphere? The answer is correct, return */
- if( s.a == s.b )
+ if( use_sphere )
{
return distance;
}
/* On a spheroid? Confirm that we are *actually* closer than tolerance */
else
{
- d = spheroid_distance(g1, nearest2, s);
+ d = spheroid_distance(&g1, &nearest2, s);
/* Yes, closer than tolerance, return! */
if( d < tolerance )
return d;
e1.start = e1.end;
}
/* On sphere, return answer */
- if( s.a == s.b )
+ if( use_sphere )
return distance;
/* On spheroid, calculate final answer based on closest approach */
else
- return spheroid_distance(g1, nearest2, s);
+ return spheroid_distance(&g1, &nearest2, s);
}
/* Initialize start of line 1 */
LWDEBUGF(4, "e2.start == GPOINT(%.6g %.6g) ", e2.start.lat, e2.start.lon);
LWDEBUGF(4, "e2.end == GPOINT(%.6g %.6g) ", e2.end.lat, e2.end.lon);
- if ( check_intersection && edge_intersection(e1, e2, &g1) )
+ if ( check_intersection && edge_intersection(&e1, &e2, &g1) )
{
LWDEBUG(4,"edge intersection! returning 0.0");
return 0.0;
}
- d = s.radius * edge_distance_to_edge(e1, e2, &g1, &g2);
+ d = s->radius * edge_distance_to_edge(&e1, &e2, &g1, &g2);
LWDEBUGF(4,"got edge_distance_to_edge %.8g", d);
if( d < distance )
}
if( d < tolerance )
{
- if( s.a == s.b )
+ if( use_sphere )
{
return d;
}
else
{
- d = spheroid_distance(nearest1, nearest2, s);
+ d = spheroid_distance(&nearest1, &nearest2, s);
if( d < tolerance )
return d;
}
}
LWDEBUGF(4,"finished all loops, returning %.8g", distance);
- if( s.a == s.b )
+ if( use_sphere )
return distance;
else
- return spheroid_distance(nearest1, nearest2, s);
+ return spheroid_distance(&nearest1, &nearest2, s);
}
* calculate external ring area and subtract internal ring area. A GBOX is
* required to calculate an outside point.
*/
-double lwgeom_area_sphere(LWGEOM *lwgeom, GBOX gbox, SPHEROID spheroid)
+double lwgeom_area_sphere(const LWGEOM *lwgeom, const GBOX *gbox, const SPHEROID *spheroid)
{
int type;
POINT2D pt_outside;
- double radius2 = spheroid.radius * spheroid.radius;
+ double radius2 = spheroid->radius * spheroid->radius;
assert(lwgeom);
return 0.0;
/* First, the area of the outer ring */
- area += radius2 * ptarray_area_sphere(poly->rings[0], pt_outside);
+ area += radius2 * ptarray_area_sphere(poly->rings[0], &pt_outside);
/* Subtract areas of inner rings */
for( i = 1; i < poly->nrings; i++ )
{
- area -= radius2 * ptarray_area_sphere(poly->rings[i], pt_outside);
+ area -= radius2 * ptarray_area_sphere(poly->rings[i], &pt_outside);
}
return area;
}
* below the tolerance (useful for dwithin calculations).
* Return a negative distance for incalculable cases.
*/
-double lwgeom_distance_spheroid(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox1, GBOX gbox2, SPHEROID spheroid, double tolerance)
+double lwgeom_distance_spheroid(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2, const GBOX *gbox1, const GBOX *gbox2, const SPHEROID *spheroid, double tolerance)
{
int type1, type2;
int check_intersection = LW_FALSE;
POINT2D p;
LWPOLY *lwpoly;
LWPOINT *lwpt;
- GBOX gbox;
+ const GBOX *gbox;
double distance = MAXFLOAT;
int i;
getPoint2d_p(lwpt->point, 0, &p);
/* Point in polygon implies zero distance */
- if( lwpoly_covers_point2d(lwpoly, gbox, p) )
+ if( lwpoly_covers_point2d(lwpoly, gbox, &p) )
return 0.0;
/* Not inside, so what's the actual distance? */
POINT2D p;
LWPOLY *lwpoly;
LWLINE *lwline;
- GBOX gbox;
+ const GBOX *gbox;
double distance = MAXFLOAT;
int i;
LWDEBUG(4, "checking if a point of line is in polygon");
/* Point in polygon implies zero distance */
- if( lwpoly_covers_point2d(lwpoly, gbox, p) )
+ if( lwpoly_covers_point2d(lwpoly, gbox, &p) )
return 0.0;
LWDEBUG(4, "checking ring distances");
/* Point of 2 in polygon 1 implies zero distance */
getPoint2d_p(lwpoly1->rings[0], 0, &p);
- if( lwpoly_covers_point2d(lwpoly2, gbox2, p) )
+ if( lwpoly_covers_point2d(lwpoly2, gbox2, &p) )
return 0.0;
/* Point of 1 in polygon 2 implies zero distance */
getPoint2d_p(lwpoly2->rings[0], 0, &p);
- if( lwpoly_covers_point2d(lwpoly1, gbox1, p) )
+ if( lwpoly_covers_point2d(lwpoly1, gbox1, &p) )
return 0.0;
/* Not contained, so what's the actual distance? */
}
-int lwgeom_covers_lwgeom_sphere(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2, GBOX gbox1, GBOX gbox2)
+int lwgeom_covers_lwgeom_sphere(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2, const GBOX *gbox1, const GBOX *gbox2)
{
int type1, type2;
{
POINT2D pt_to_test;
getPoint2d_p(((LWPOINT*)lwgeom2)->point, 0, &pt_to_test);
- return lwpoly_covers_point2d((LWPOLY*)lwgeom1, gbox1, pt_to_test);
+ return lwpoly_covers_point2d((LWPOLY*)lwgeom1, gbox1, &pt_to_test);
}
/* If any of the first argument parts covers the second argument, it's true */
* a guaranteed outside point (lon/lat decimal degrees) (calculate with gbox_pt_outside())
* return LW_TRUE if point is inside or on edge of polygon.
*/
-int lwpoly_covers_point2d(const LWPOLY *poly, GBOX gbox, POINT2D pt_to_test)
+int lwpoly_covers_point2d(const LWPOLY *poly, const GBOX *gbox, const POINT2D *pt_to_test)
{
int i;
int in_hole_count = 0;
}
/* Point not in box? Done! */
- geographic_point_init(pt_to_test.x, pt_to_test.y, &gpt_to_test);
- geog2cart(gpt_to_test, &p);
- if( ! gbox_contains_point3d(gbox, p) )
+ geographic_point_init(pt_to_test->x, pt_to_test->y, &gpt_to_test);
+ geog2cart(&gpt_to_test, &p);
+ if( ! gbox_contains_point3d(gbox, &p) )
return LW_FALSE;
/* Calculate our outside point from the gbox */
gbox_pt_outside(gbox, &pt_outside);
LWDEBUGF(4, "pt_outside POINT(%.18g %.18g)", pt_outside.x, pt_outside.y);
- LWDEBUGF(4, "pt_to_test POINT(%.18g %.18g)", pt_to_test.x, pt_to_test.y);
+ LWDEBUGF(4, "pt_to_test POINT(%.18g %.18g)", pt_to_test->x, pt_to_test->y);
LWDEBUGF(4, "polygon %s", lwgeom_to_ewkt((LWGEOM*)poly, PARSER_CHECK_NONE));
- LWDEBUGF(4, "gbox %s", gbox_to_string(&gbox));
+ LWDEBUGF(4, "gbox %s", gbox_to_string(gbox));
/* Not in outer ring? We're done! */
- if( ! ptarray_point_in_ring(poly->rings[0], pt_outside, pt_to_test) )
+ if( ! ptarray_point_in_ring(poly->rings[0], &pt_outside, pt_to_test) )
{
LWDEBUG(4,"returning false, point is outside ring");
return LW_FALSE;
{
LWDEBUGF(4, "ring test loop %d", i);
/* Count up hole containment. Odd => outside boundary. */
- if( ptarray_point_in_ring(poly->rings[i], pt_outside, pt_to_test) )
+ if( ptarray_point_in_ring(poly->rings[i], &pt_outside, pt_to_test) )
in_hole_count++;
}
return G_SUCCESS;
}
-int ptarray_calculate_gbox_geodetic(POINTARRAY *pa, GBOX *gbox)
+int ptarray_calculate_gbox_geodetic(const POINTARRAY *pa, GBOX *gbox)
{
int i;
int first = LW_TRUE;
GEOGRAPHIC_POINT gp;
getPoint2d_p(pa, 0, &in_pt);
geographic_point_init(in_pt.x, in_pt.y, &gp);
- geog2cart(gp, &out_pt);
+ geog2cart(&gp, &out_pt);
gbox->xmin = gbox->xmax = out_pt.x;
gbox->ymin = gbox->ymax = out_pt.y;
gbox->zmin = gbox->zmax = out_pt.z;
getPoint2d_p(pa, i, &end_pt);
geographic_point_init(end_pt.x, end_pt.y, &(edge.end));
- edge_calculate_gbox(edge, &edge_gbox);
+ edge_calculate_gbox(&edge, &edge_gbox);
LWDEBUGF(4, "edge_gbox: %s", gbox_to_string(&edge_gbox));
/* Initialize the box */
if ( first )
{
- gbox_duplicate(edge_gbox,gbox);
+ gbox_duplicate(&edge_gbox, gbox);
LWDEBUGF(4, "gbox_duplicate: %s", gbox_to_string(gbox));
first = LW_FALSE;
}
/* Expand the box where necessary */
else
{
- gbox_merge(edge_gbox, gbox);
+ gbox_merge(&edge_gbox, gbox);
LWDEBUGF(4, "gbox_merge: %s", gbox_to_string(gbox));
}
}
-static int lwpoint_calculate_gbox_geodetic(LWPOINT *point, GBOX *gbox)
+static int lwpoint_calculate_gbox_geodetic(const LWPOINT *point, GBOX *gbox)
{
assert(point);
if ( ptarray_calculate_gbox_geodetic(point->point, gbox) == G_FAILURE )
return G_SUCCESS;
}
-static int lwline_calculate_gbox_geodetic(LWLINE *line, GBOX *gbox)
+static int lwline_calculate_gbox_geodetic(const LWLINE *line, GBOX *gbox)
{
assert(line);
if ( ptarray_calculate_gbox_geodetic(line->points, gbox) == G_FAILURE )
return G_SUCCESS;
}
-static int lwpolygon_calculate_gbox_geodetic(LWPOLY *poly, GBOX *gbox)
+static int lwpolygon_calculate_gbox_geodetic(const LWPOLY *poly, GBOX *gbox)
{
GBOX ringbox;
int i;
return G_FAILURE;
if ( first )
{
- gbox_duplicate(ringbox, gbox);
+ gbox_duplicate(&ringbox, gbox);
first = LW_FALSE;
}
else
{
- gbox_merge(ringbox, gbox);
+ gbox_merge(&ringbox, gbox);
}
}
return G_SUCCESS;
}
-static int lwcollection_calculate_gbox_geodetic(LWCOLLECTION *coll, GBOX *gbox)
+static int lwcollection_calculate_gbox_geodetic(const LWCOLLECTION *coll, GBOX *gbox)
{
GBOX subbox;
int i;
{
if ( first )
{
- gbox_duplicate(subbox, gbox);
+ gbox_duplicate(&subbox, gbox);
first = LW_FALSE;
}
else
{
- gbox_merge(subbox, gbox);
+ gbox_merge(&subbox, gbox);
}
result = G_SUCCESS;
}
-static int ptarray_check_geodetic(POINTARRAY *pa)
+static int ptarray_check_geodetic(const POINTARRAY *pa)
{
int t;
POINT2D pt;
return LW_TRUE;
}
-static int lwpoint_check_geodetic(LWPOINT *point)
+static int lwpoint_check_geodetic(const LWPOINT *point)
{
assert(point);
return ptarray_check_geodetic(point->point);
}
-static int lwline_check_geodetic(LWLINE *line)
+static int lwline_check_geodetic(const LWLINE *line)
{
assert(line);
return ptarray_check_geodetic(line->points);
}
-static int lwpoly_check_geodetic(LWPOLY *poly)
+static int lwpoly_check_geodetic(const LWPOLY *poly)
{
int i = 0;
assert(poly);
return LW_TRUE;
}
-static int lwcollection_check_geodetic(LWCOLLECTION *col)
+static int lwcollection_check_geodetic(const LWCOLLECTION *col)
{
int i = 0;
assert(col);
return LW_FALSE;
}
-double ptarray_length_spheroid(POINTARRAY *pa, SPHEROID s)
+double ptarray_length_spheroid(const POINTARRAY *pa, const SPHEROID *s)
{
GEOGRAPHIC_POINT a, b;
POINT2D p;
geographic_point_init(p.x, p.y, &b);
/* Special sphere case */
- if( s.a == s.b )
- length += s.radius * sphere_distance(a, b);
+ if( s->a == s->b )
+ length += s->radius * sphere_distance(&a, &b);
/* Spheroid case */
else
- length += spheroid_distance(a, b, s);
+ length += spheroid_distance(&a, &b, s);
/* B gets incremented in the next loop, so we save the value here */
a = b;
return length;
}
-double lwgeom_length_spheroid(LWGEOM *geom, SPHEROID s)
+double lwgeom_length_spheroid(const LWGEOM *geom, const SPHEROID *s)
{
int type;
int i = 0;
/*
** Prototypes for internal functions.
*/
-void geog2cart(GEOGRAPHIC_POINT g, POINT3D *p);
-void cart2geog(POINT3D p, GEOGRAPHIC_POINT *g);
-void robust_cross_product(GEOGRAPHIC_POINT p, GEOGRAPHIC_POINT q, POINT3D *a);
+
+void geog2cart(const GEOGRAPHIC_POINT *g, POINT3D *p);
+void cart2geog(const POINT3D *p, GEOGRAPHIC_POINT *g);
+void robust_cross_product(const GEOGRAPHIC_POINT *p, const 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_point_in_cone(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p);
-int edge_contains_coplanar_point(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p);
-int edge_contains_point(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p);
+int edge_point_on_plane(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p);
+int edge_point_in_cone(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p);
+int edge_contains_coplanar_point(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p);
+int edge_contains_point(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p);
double z_to_latitude(double z, int top);
-int clairaut_cartesian(POINT3D start, POINT3D end, GEOGRAPHIC_POINT *g_top, GEOGRAPHIC_POINT *g_bottom);
-int clairaut_geographic(GEOGRAPHIC_POINT start, GEOGRAPHIC_POINT end, GEOGRAPHIC_POINT *g_top, GEOGRAPHIC_POINT *g_bottom);
-double sphere_distance(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e);
-double sphere_distance_cartesian(POINT3D s, POINT3D e);
-int sphere_project(GEOGRAPHIC_POINT r, double distance, double azimuth, GEOGRAPHIC_POINT *n);
-int edge_calculate_gbox(GEOGRAPHIC_EDGE e, GBOX *gbox);
-int edge_calculate_gbox_slow(GEOGRAPHIC_EDGE e, GBOX *gbox);
-int edge_intersection(GEOGRAPHIC_EDGE e1, GEOGRAPHIC_EDGE e2, GEOGRAPHIC_POINT *g);
-double edge_distance_to_point(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT gp, GEOGRAPHIC_POINT *closest);
-double edge_distance_to_edge(GEOGRAPHIC_EDGE e1, GEOGRAPHIC_EDGE e2, GEOGRAPHIC_POINT *closest1, GEOGRAPHIC_POINT *closest2);
+int clairaut_cartesian(const POINT3D *start, const POINT3D *end, GEOGRAPHIC_POINT *g_top, GEOGRAPHIC_POINT *g_bottom);
+int clairaut_geographic(const GEOGRAPHIC_POINT *start, const GEOGRAPHIC_POINT *end, GEOGRAPHIC_POINT *g_top, GEOGRAPHIC_POINT *g_bottom);
+double sphere_distance(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e);
+double sphere_distance_cartesian(const POINT3D *s, const POINT3D *e);
+int sphere_project(const GEOGRAPHIC_POINT *r, double distance, double azimuth, GEOGRAPHIC_POINT *n);
+int edge_calculate_gbox(const GEOGRAPHIC_EDGE *e, GBOX *gbox);
+int edge_calculate_gbox_slow(const GEOGRAPHIC_EDGE *e, GBOX *gbox);
+int edge_intersection(const GEOGRAPHIC_EDGE *e1, const GEOGRAPHIC_EDGE *e2, GEOGRAPHIC_POINT *g);
+double edge_distance_to_point(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *gp, GEOGRAPHIC_POINT *closest);
+double edge_distance_to_edge(const GEOGRAPHIC_EDGE *e1, const GEOGRAPHIC_EDGE *e2, GEOGRAPHIC_POINT *closest1, GEOGRAPHIC_POINT *closest2);
void edge_deg2rad(GEOGRAPHIC_EDGE *e);
void edge_rad2deg(GEOGRAPHIC_EDGE *e);
void point_deg2rad(GEOGRAPHIC_POINT *p);
void point_rad2deg(GEOGRAPHIC_POINT *p);
void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g);
-int ptarray_point_in_ring_winding(POINTARRAY *pa, POINT2D pt_to_test);
-int lwpoly_covers_point2d(const LWPOLY *poly, GBOX gbox, POINT2D pt_to_test);
-int ptarray_point_in_ring(POINTARRAY *pa, POINT2D pt_outside, POINT2D pt_to_test);
-double ptarray_area_sphere(POINTARRAY *pa, POINT2D pt_outside);
+int ptarray_point_in_ring_winding(const POINTARRAY *pa, const POINT2D *pt_to_test);
+int lwpoly_covers_point2d(const LWPOLY *poly, const GBOX *gbox, const POINT2D *pt_to_test);
+int ptarray_point_in_ring(const POINTARRAY *pa, const POINT2D *pt_outside, const POINT2D *pt_to_test);
+double ptarray_area_sphere(const POINTARRAY *pa, const POINT2D *pt_outside);
double latitude_degrees_normalize(double lat);
double longitude_degrees_normalize(double lon);
-double ptarray_length_spheroid(POINTARRAY *pa, SPHEROID s);
-int geographic_point_equals(GEOGRAPHIC_POINT g1, GEOGRAPHIC_POINT g2);
-int crosses_dateline(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e);
+double ptarray_length_spheroid(const POINTARRAY *pa, const SPHEROID *s);
+int geographic_point_equals(const GEOGRAPHIC_POINT *g1, const GEOGRAPHIC_POINT *g2);
+int crosses_dateline(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e);
void point_shift(GEOGRAPHIC_POINT *p, double shift);
/*
** Prototypes for spheroid functions.
*/
-double spheroid_distance(GEOGRAPHIC_POINT a, GEOGRAPHIC_POINT b, SPHEROID spheroid);
-double spheroid_direction(GEOGRAPHIC_POINT r, GEOGRAPHIC_POINT s, SPHEROID spheroid);
-int spheroid_project(GEOGRAPHIC_POINT r, SPHEROID spheroid, double distance, double azimuth, GEOGRAPHIC_POINT *g);
+double spheroid_distance(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const SPHEROID *spheroid);
+double spheroid_direction(const GEOGRAPHIC_POINT *r, const GEOGRAPHIC_POINT *s, const SPHEROID *spheroid);
+int spheroid_project(const GEOGRAPHIC_POINT *r, const SPHEROID *spheroid, double distance, double azimuth, GEOGRAPHIC_POINT *g);
};
-int lwgeom_needs_bbox(LWGEOM *geom)
+int lwgeom_needs_bbox(const LWGEOM *geom)
{
assert(geom);
if( TYPE_GETTYPE(geom->type) == POINTTYPE )
return result;
}
-static int lwpoint_is_empty(LWPOINT *point)
+static int lwpoint_is_empty(const LWPOINT *point)
{
if( ! point->point || point->point->npoints == 0 )
return LW_TRUE;
return LW_FALSE;
}
-static int lwline_is_empty(LWLINE *line)
+static int lwline_is_empty(const LWLINE *line)
{
if( !line->points || line->points->npoints == 0 )
return LW_TRUE;
return LW_FALSE;
}
-static int lwpoly_is_empty(LWPOLY *poly)
+static int lwpoly_is_empty(const LWPOLY *poly)
{
if( !poly->rings || poly->nrings == 0 )
return LW_TRUE;
return LW_FALSE;
}
-static int lwcircstring_is_empty(LWCIRCSTRING *circ)
+static int lwcircstring_is_empty(const LWCIRCSTRING *circ)
{
if( !circ->points || circ->points->npoints == 0 )
return LW_TRUE;
return LW_FALSE;
}
-static int lwcollection_is_empty(LWCOLLECTION *col)
+static int lwcollection_is_empty(const LWCOLLECTION *col)
{
if( !col->geoms || col->ngeoms == 0 )
return LW_TRUE;
return LW_FALSE;
}
-int lwgeom_is_empty(LWGEOM *geom)
+int lwgeom_is_empty(const LWGEOM *geom)
{
int result = LW_FALSE;
LWDEBUGF(4, "got type %d", TYPE_GETTYPE(geom->type));
s->radius = (2.0 * a + b ) / 3.0;
}
-static double spheroid_mu2(double alpha, SPHEROID s)
+static double spheroid_mu2(double alpha, const SPHEROID *s)
{
- double b2 = POW2(s.b);
- return POW2(cos(alpha)) * (POW2(s.a) - b2) / b2;
+ double b2 = POW2(s->b);
+ return POW2(cos(alpha)) * (POW2(s->a) - b2) / b2;
}
static double spheroid_big_a(double u2)
* @param s - spheroid to calculate on
* @return spheroidal distance between a and b in spheroid units.
*/
-double spheroid_distance(GEOGRAPHIC_POINT a, GEOGRAPHIC_POINT b, SPHEROID spheroid)
+double spheroid_distance(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const SPHEROID *spheroid)
{
- double lambda = (b.lon - a.lon);
- double f = spheroid.f;
- double omf = 1 - spheroid.f;
+ double lambda = (b->lon - a->lon);
+ double f = spheroid->f;
+ double omf = 1 - spheroid->f;
double u1, u2;
double cos_u1, cos_u2;
double sin_u1, sin_u2;
return 0.0;
}
- u1 = atan(omf * tan(a.lat));
+ u1 = atan(omf * tan(a->lat));
cos_u1 = cos(u1);
sin_u1 = sin(u1);
- u2 = atan(omf * tan(b.lat));
+ u2 = atan(omf * tan(b->lat));
cos_u2 = cos(u2);
sin_u2 = sin(u2);
delta_sigma = big_b * sin_sigma * (cos2_sigma_m + (big_b / 4.0) * (cos_sigma * (-1.0 + 2.0 * POW2(cos2_sigma_m)) -
(big_b / 6.0) * cos2_sigma_m * (-3.0 + 4.0 * sqrsin_sigma) * (-3.0 + 4.0 * POW2(cos2_sigma_m))));
- distance = spheroid.b * big_a * (sigma - delta_sigma);
+ distance = spheroid->b * big_a * (sigma - delta_sigma);
/* Algorithm failure, distance == NaN, fallback to sphere */
if( distance != distance )
{
- lwerror("spheroid_distance returned NaN: (%.20g %.20g) (%.20g %.20g) a = %.20g b = %.20g",a.lat, a.lon, b.lat, b.lon, spheroid.a, spheroid.b);
- return spheroid.radius * sphere_distance(a, b);
+ lwerror("spheroid_distance returned NaN: (%.20g %.20g) (%.20g %.20g) a = %.20g b = %.20g",a->lat, a->lon, b->lat, b->lon, spheroid->a, spheroid->b);
+ return spheroid->radius * sphere_distance(a, b);
}
return distance;
* @param s - location of second point
* @return azimuth of line joining r and s
*/
-double spheroid_direction(GEOGRAPHIC_POINT r, GEOGRAPHIC_POINT s, SPHEROID spheroid)
+double spheroid_direction(const GEOGRAPHIC_POINT *r, const GEOGRAPHIC_POINT *s, const SPHEROID *spheroid)
{
int i = 0;
- double lambda = s.lon - r.lon;
- double omf = 1 - spheroid.f;
- double u1 = atan(omf * tan(r.lat));
+ double lambda = s->lon - r->lon;
+ double omf = 1 - spheroid->f;
+ double u1 = atan(omf * tan(r->lat));
double cos_u1 = cos(u1);
double sin_u1 = sin(u1);
- double u2 = atan(omf * tan(s.lat));
+ double u2 = atan(omf * tan(s->lat));
double cos_u2 = cos(u2);
double sin_u2 = sin(u2);
if( cos2_sigma_m < -1.0 )
cos2_sigma_m = -1.0;
- C = (spheroid.f / 16.0) * cos_alphasq * (4.0 + spheroid.f * (4.0 - 3.0 * cos_alphasq));
+ C = (spheroid->f / 16.0) * cos_alphasq * (4.0 + spheroid->f * (4.0 - 3.0 * cos_alphasq));
last_lambda = lambda;
- lambda = omega + (1.0 - C) * spheroid.f * sin(alpha) * (sigma + C * sin(sigma) *
+ lambda = omega + (1.0 - C) * spheroid->f * sin(alpha) * (sigma + C * sin(sigma) *
(cos2_sigma_m + C * cos(sigma) * (-1.0 + 2.0 * POW2(cos2_sigma_m))));
i++;
}
* @param azimuth - azimuth in radians.
* @return s - location of projected point.
*/
-int spheroid_project(GEOGRAPHIC_POINT r, SPHEROID spheroid, double distance, double azimuth, GEOGRAPHIC_POINT *g)
+int spheroid_project(const GEOGRAPHIC_POINT *r, const SPHEROID *spheroid, double distance, double azimuth, GEOGRAPHIC_POINT *g)
{
- double omf = 1 - spheroid.f;
- double tan_u1 = omf * tan(r.lat);
+ double omf = 1 - spheroid->f;
+ double tan_u1 = omf * tan(r->lat);
double u1 = atan(tan_u1);
double sigma, last_sigma, delta_sigma, two_sigma_m;
double sigma1, sin_alpha, alpha, cos_alphasq;
A = spheroid_big_a(u2);
B = spheroid_big_b(u2);
- sigma = (distance / (spheroid.b * A));
+ sigma = (distance / (spheroid->b * A));
do
{
two_sigma_m = 2.0 * sigma1 + sigma;
delta_sigma = B * sin(sigma) * (cos(two_sigma_m) + (B / 4.0) * (cos(sigma) * (-1.0 + 2.0 * POW2(cos(two_sigma_m)) - (B / 6.0) * cos(two_sigma_m) * (-3.0 + 4.0 * POW2(sin(sigma))) * (-3.0 + 4.0 * POW2(cos(two_sigma_m))))));
last_sigma = sigma;
- sigma = (distance / (spheroid.b * A)) + delta_sigma;
+ sigma = (distance / (spheroid->b * A)) + delta_sigma;
i++;
}
while (i < 999 && fabs((last_sigma - sigma) / sigma) > 1.0e-9);
cos(azimuth)))));
lambda = atan2((sin(sigma) * sin(azimuth)), (cos(u1) * cos(sigma) -
sin(u1) * sin(sigma) * cos(azimuth)));
- C = (spheroid.f / 16.0) * cos_alphasq * (4.0 + spheroid.f * (4.0 - 3.0 * cos_alphasq));
- omega = lambda - (1.0 - C) * spheroid.f * sin_alpha * (sigma + C * sin(sigma) *
+ C = (spheroid->f / 16.0) * cos_alphasq * (4.0 + spheroid->f * (4.0 - 3.0 * cos_alphasq));
+ omega = lambda - (1.0 - C) * spheroid->f * sin_alpha * (sigma + C * sin(sigma) *
(cos(two_sigma_m) + C * cos(sigma) * (-1.0 + 2.0 * POW2(cos(two_sigma_m)))));
- lambda2 = r.lon + omega;
+ lambda2 = r->lon + omega;
g->lat = lat2;
g->lon = lambda2;
return G_SUCCESS;
}
-static inline double spheroid_prime_vertical_radius_of_curvature(double latitude, SPHEROID spheroid)
+static inline double spheroid_prime_vertical_radius_of_curvature(double latitude, const SPHEROID *spheroid)
{
- return spheroid.a / (sqrt(1.0 - spheroid.e_sq * POW2(sin(latitude))));
+ return spheroid->a / (sqrt(1.0 - spheroid->e_sq * POW2(sin(latitude))));
}
-static inline double spheroid_parallel_arc_length(double latitude, double deltaLongitude, SPHEROID spheroid)
+static inline double spheroid_parallel_arc_length(double latitude, double deltaLongitude, const SPHEROID *spheroid)
{
return spheroid_prime_vertical_radius_of_curvature(latitude, spheroid)
* cos(latitude)
* @param northEastCorner - upper right corner of bounding box.
* @return area in square meters.
*/
-static double spheroid_boundingbox_area(GEOGRAPHIC_POINT southWestCorner, GEOGRAPHIC_POINT northEastCorner, SPHEROID spheroid)
+static double spheroid_boundingbox_area(const GEOGRAPHIC_POINT *southWestCorner, const GEOGRAPHIC_POINT *northEastCorner, const SPHEROID *spheroid)
{
- double z0 = (northEastCorner.lon - southWestCorner.lon) * POW2(spheroid.b) / 2.0;
- double e = sqrt(spheroid.e_sq);
- double sinPhi1 = sin(southWestCorner.lat);
- double sinPhi2 = sin(northEastCorner.lat);
- double t1p1 = sinPhi1 / (1.0 - spheroid.e_sq * sinPhi1 * sinPhi1);
- double t1p2 = sinPhi2 / (1.0 - spheroid.e_sq * sinPhi2 * sinPhi2);
+ double z0 = (northEastCorner->lon - southWestCorner->lon) * POW2(spheroid->b) / 2.0;
+ double e = sqrt(spheroid->e_sq);
+ double sinPhi1 = sin(southWestCorner->lat);
+ double sinPhi2 = sin(northEastCorner->lat);
+ double t1p1 = sinPhi1 / (1.0 - spheroid->e_sq * sinPhi1 * sinPhi1);
+ double t1p2 = sinPhi2 / (1.0 - spheroid->e_sq * sinPhi2 * sinPhi2);
double oneOver2e = 1.0 / (2.0 * e);
double t2p1 = oneOver2e * log((1.0 + e * sinPhi1) / (1.0 - e * sinPhi1));
double t2p2 = oneOver2e * log((1.0 + e * sinPhi2) / (1.0 - e * sinPhi2));
* This function doesn't work for edges crossing the dateline or in the southern
* hemisphere. Points are pre-conditioned in ptarray_area_spheroid.
*/
-static double spheroid_striparea(GEOGRAPHIC_POINT a, GEOGRAPHIC_POINT b, double latitude_min, SPHEROID spheroid)
+static double spheroid_striparea(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, double latitude_min, const SPHEROID *spheroid)
{
- GEOGRAPHIC_POINT A = a;
- GEOGRAPHIC_POINT B = b;
+ GEOGRAPHIC_POINT A = *a;
+ GEOGRAPHIC_POINT B = *b;
GEOGRAPHIC_POINT mL, nR;
double deltaLng, baseArea, topArea;
double bE, tE, ratio, sign;
mL.lon = FP_MIN(A.lon, B.lon);
nR.lat = FP_MIN(A.lat, B.lat);
nR.lon = FP_MAX(A.lon, B.lon);
- baseArea = spheroid_boundingbox_area(mL, nR, spheroid);
+ baseArea = spheroid_boundingbox_area(&mL, &nR, spheroid);
mL.lat = FP_MIN(A.lat, B.lat);
mL.lon = FP_MIN(A.lon, B.lon);
nR.lat = FP_MAX(A.lat, B.lat);
nR.lon = FP_MAX(A.lon, B.lon);
- topArea = spheroid_boundingbox_area(mL, nR, spheroid);
+ topArea = spheroid_boundingbox_area(&mL, &nR, spheroid);
deltaLng = B.lon - A.lon;
bE = spheroid_parallel_arc_length(A.lat, deltaLng, spheroid);
return (baseArea + topArea / ratio) * sign;
}
-static double ptarray_area_spheroid(POINTARRAY *pa, SPHEROID spheroid)
+static double ptarray_area_spheroid(const POINTARRAY *pa, const SPHEROID *spheroid)
{
GEOGRAPHIC_POINT a, b;
POINT2D p;
LWDEBUGF(4, "in_south %d", in_south);
- if( crosses_dateline(a, b) )
+ if( crosses_dateline(&a, &b) )
{
double shift;
point_shift(&b1, shift);
}
- LWDEBUGF(4, "crosses_dateline(a, b) %d", crosses_dateline(a, b) );
+ LWDEBUGF(4, "crosses_dateline(a, b) %d", crosses_dateline(&a, &b) );
delta_lon = fabs(b1.lon - a1.lon);
{
if( delta_lon < delta_lon_tolerance )
{
- strip_area = spheroid_striparea(a1, b1, latitude_min, spheroid);
+ strip_area = spheroid_striparea(&a1, &b1, latitude_min, spheroid);
LWDEBUGF(4, "strip_area %.12g", strip_area);
area += strip_area;
}
{
GEOGRAPHIC_POINT p, q;
double step = floor(delta_lon / delta_lon_tolerance);
- double distance = spheroid_distance(a1, b1, spheroid);
+ double distance = spheroid_distance(&a1, &b1, spheroid);
double pDistance = 0.0;
int j = 0;
LWDEBUGF(4, "step %.18g", step);
p = a1;
while (pDistance < (distance - step * 1.01))
{
- double azimuth = spheroid_direction(p, b1, spheroid);
+ double azimuth = spheroid_direction(&p, &b1, spheroid);
j++;
LWDEBUGF(4, " iteration %d", j);
LWDEBUGF(4, " azimuth %.12g", azimuth);
pDistance = pDistance + step;
LWDEBUGF(4, " pDistance %.12g", pDistance);
- spheroid_project(p, spheroid, step, azimuth, &q);
- strip_area = spheroid_striparea(p, q, latitude_min, spheroid);
+ spheroid_project(&p, spheroid, step, azimuth, &q);
+ strip_area = spheroid_striparea(&p, &q, latitude_min, spheroid);
LWDEBUGF(4, " strip_area %.12g", strip_area);
area += strip_area;
LWDEBUGF(4, " area %.12g", area);
p.lat = q.lat;
p.lon = q.lon;
}
- strip_area = spheroid_striparea(p, b1, latitude_min, spheroid);
+ strip_area = spheroid_striparea(&p, &b1, latitude_min, spheroid);
area += strip_area;
}
}
* required to check relationship to equator an outside point.
* WARNING: Does NOT WORK for polygons over equator or pole.
*/
-double lwgeom_area_spheroid(LWGEOM *lwgeom, GBOX gbox, SPHEROID spheroid)
+double lwgeom_area_spheroid(const LWGEOM *lwgeom, const GBOX *gbox, const SPHEROID *spheroid)
{
int type;
PG_RETURN_NULL();
}
- distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, gbox1, gbox2, s, 0.0);
+ distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, &gbox1, &gbox2, &s, 0.0);
/* Something went wrong, negative return... should already be eloged, return NULL */
if( distance < 0.0 )
PG_RETURN_BOOL(FALSE);
}
- distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, gbox1, gbox2, s, tolerance);
+ distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, &gbox1, &gbox2, &s, tolerance);
/* Something went wrong... should already be eloged, return FALSE */
if( distance < 0.0 )
/* Calculate the area */
if( use_spheroid )
- area = lwgeom_area_spheroid(lwgeom, gbox, s);
+ area = lwgeom_area_spheroid(lwgeom, &gbox, &s);
else
- area = lwgeom_area_sphere(lwgeom, gbox, s);
+ area = lwgeom_area_sphere(lwgeom, &gbox, &s);
/* Something went wrong... */
if( area < 0.0 )
s.a = s.b = s.radius;
/* Calculate the length */
- length = lwgeom_length_spheroid(lwgeom, s);
+ length = lwgeom_length_spheroid(lwgeom, &s);
/* Something went wrong... */
if( length < 0.0 )
}
/* Get an exterior point, based on this gbox */
- gbox_pt_outside(gbox, &pt);
+ gbox_pt_outside(&gbox, &pt);
lwpoint = make_lwpoint2d(4326, pt.x, pt.y);
}
/* Calculate answer */
- result = lwgeom_covers_lwgeom_sphere(lwgeom1, lwgeom2, gbox1, gbox2);
+ result = lwgeom_covers_lwgeom_sphere(lwgeom1, lwgeom2, &gbox1, &gbox2);
/* Clean up, but not all the way to the point arrays */
lwgeom_release(lwgeom1);
PG_LWGEOM *geom = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
SPHEROID *sphere = (SPHEROID *) PG_GETARG_POINTER(1);
LWGEOM *lwgeom = lwgeom_deserialize(SERIALIZED_FORM(geom));
- double dist = lwgeom_length_spheroid(lwgeom, *sphere);
+ double dist = lwgeom_length_spheroid(lwgeom, sphere);
lwgeom_release(lwgeom);
PG_RETURN_FLOAT8(dist);
}
PG_RETURN_NULL();
};
- distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, gbox1, gbox2, *sphere, 0.0);
+ distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, &gbox1, &gbox2, sphere, 0.0);
PG_RETURN_FLOAT8(distance);
projects / postgis / commitdiff