Skip to content

Commit

Permalink
Distance calculation support for arc features (#2018). Commit adds lo…
Browse files Browse the repository at this point in the history 
…west level primitive support for distance calculations on single arcs.

git-svn-id: http://svn.osgeo.org/postgis/trunk@10335 b70326c6-7e19-0410-871a-916f4a2858ee
  • Loading branch information
@pramsey
pramsey committed Sep 27, 2012
1 parent 7a2faca commit 503d866
Show file tree
Hide file tree
Showing 7 changed files with 654 additions and 34 deletions.
200 changes: 200 additions & 0 deletions liblwgeom/cunit/cu_measures.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -345,6 +345,203 @@ test_lwgeom_locate_along(void)
lwfree(str);
}

static void
test_lw_dist2d_pt_arc(void)
{
/* int lw_dist2d_pt_arc(const POINT2D* P, const POINT2D* A1, const POINT2D* A2, const POINT2D* A3, DISTPTS* dl) */
DISTPTS dl;
POINT2D P, A1, A2, A3;
int rv;


/* Point within unit semicircle, 0.5 units from arc */
A1.x = -1; A1.y = 0;
A2.x = 0 ; A2.y = 1;
A3.x = 1 ; A3.y = 0;
P.x = 0 ; P.y = 0.5;

lw_dist2d_distpts_init(&dl, DIST_MIN);
rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);

/* Point outside unit semicircle, 0.5 units from arc */
P.x = 0 ; P.y = 1.5;
lw_dist2d_distpts_init(&dl, DIST_MIN);
rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);

/* Point outside unit semicircle, sqrt(2) units from arc end point*/
P.x = 0 ; P.y = -1;
lw_dist2d_distpts_init(&dl, DIST_MIN);
rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0), 0.000001);

/* Point outside unit semicircle, sqrt(2)-1 units from arc end point*/
P.x = 1 ; P.y = 1;
lw_dist2d_distpts_init(&dl, DIST_MIN);
rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1, 0.000001);

/* Point on unit semicircle midpoint */
P.x = 0 ; P.y = 1;
lw_dist2d_distpts_init(&dl, DIST_MIN);
rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);

/* Point on unit semicircle endpoint */
P.x = 1 ; P.y = 0;
lw_dist2d_distpts_init(&dl, DIST_MIN);
rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);

/* Point inside closed circle */
P.x = 0 ; P.y = 0.5;
A2.x = 1; A2.y = 0;
A3 = A1;
lw_dist2d_distpts_init(&dl, DIST_MIN);
rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
//printf("distance %g\n", dl.distance);
CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
}

static void
test_lw_dist2d_seg_arc(void)
{
/* int lw_dist2d_seg_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *B1, const POINT2D *B2, const POINT2D *B3, DISTPTS *dl) */

DISTPTS dl;
POINT2D A1, A2, B1, B2, B3;
int rv;

/* Unit semicircle */
B1.x = -1; B1.y = 0;
B2.x = 0 ; B2.y = 1;
B3.x = 1 ; B3.y = 0;

/* Edge above the unit semicircle */
lw_dist2d_distpts_init(&dl, DIST_MIN);
A1.x = -2; A1.y = 2;
A2.x = 2 ; A2.y = 2;
rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);

/* Edge to the right of the unit semicircle */
lw_dist2d_distpts_init(&dl, DIST_MIN);
A1.x = 2; A1.y = -2;
A2.x = 2; A2.y = 2;
rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);

/* Edge to the left of the unit semicircle */
lw_dist2d_distpts_init(&dl, DIST_MIN);
A1.x = -2; A1.y = -2;
A2.x = -2; A2.y = 2;
rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);

/* Edge within the unit semicircle */
lw_dist2d_distpts_init(&dl, DIST_MIN);
A1.x = 0; A1.y = 0;
A2.x = 0; A2.y = 0.5;
rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);

/* Edge grazing the unit semicircle */
lw_dist2d_distpts_init(&dl, DIST_MIN);
A1.x = -2; A1.y = 1;
A2.x = 2; A2.y = 1;
rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0., 0.000001);

/* Line grazing the unit semicircle, but edge not */
lw_dist2d_distpts_init(&dl, DIST_MIN);
A1.x = 1; A1.y = 1;
A2.x = 2; A2.y = 1;
rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1, 0.000001);

/* Edge intersecting the unit semicircle */
lw_dist2d_distpts_init(&dl, DIST_MIN);
A1.x = 0; A1.y = 0;
A2.x = 2; A2.y = 2;
rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);

/* Line intersecting the unit semicircle, but edge not */
lw_dist2d_distpts_init(&dl, DIST_MIN);
A1.x = -1; A1.y = 1;
A2.x = -2; A2.y = 2;
rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
//printf("distance %g\n", dl.distance);
CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1, 0.000001);
}

static void
test_lw_dist2d_arc_arc(void)
{
/* lw_dist2d_arc_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3,
const POINT2D *B1, const POINT2D *B2, const POINT2D *B3,
DISTPTS *dl) */
DISTPTS dl;
POINT2D A1, A2, A3, B1, B2, B3;
int rv;

/* Unit semicircle at 0,0 */
B1.x = -1; B1.y = 0;
B2.x = 0 ; B2.y = 1;
B3.x = 1 ; B3.y = 0;

/* Arc above the unit semicircle */
lw_dist2d_distpts_init(&dl, DIST_MIN);
A1.x = -1; A1.y = 3;
A2.x = 0 ; A2.y = 2;
A3.x = 1 ; A3.y = 3;
rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);

/* Arc grazes the unit semicircle */
lw_dist2d_distpts_init(&dl, DIST_MIN);
A1.x = -1; A1.y = 2;
A2.x = 0 ; A2.y = 1;
A3.x = 1 ; A3.y = 2;
rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);

/* Circles intersect, but arcs do not */
lw_dist2d_distpts_init(&dl, DIST_MIN);
A1.x = -1; A1.y = 1;
A2.x = 0; A2.y = 2;
A3.x = 1; A3.y = 1;
rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2)-1, 0.000001);

/* Circles and arcs intersect */
lw_dist2d_distpts_init(&dl, DIST_MIN);
A1.x = -1; A1.y = 1;
A2.x = 0; A2.y = 0;
A3.x = 1; A3.y = 1;
rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);

/* inscribed and closest on arcs */
lw_dist2d_distpts_init(&dl, DIST_MIN);
A1.x = -0.5; A1.y = 0.0;
A2.x = 0.0; A2.y = 0.5;
A3.x = 0.5; A3.y = 0.0;
rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
//printf("distance %g\n", dl.distance);
CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);

/* inscribed and closest not on arcs */
lw_dist2d_distpts_init(&dl, DIST_MIN);
A1.x = -0.5; A1.y = 0.0;
A2.x = 0.0; A2.y = -0.5;
A3.x = 0.5; A3.y = 0.0;
rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
//printf("distance %g\n", dl.distance);
CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
}

/*
** Used by test harness to register the tests in this file.
*/
Expand All @@ -355,6 +552,9 @@ CU_TestInfo measures_tests[] =
PG_TEST(test_rect_tree_intersects_tree),
PG_TEST(test_lwgeom_segmentize2d),
PG_TEST(test_lwgeom_locate_along),
PG_TEST(test_lw_dist2d_pt_arc),
PG_TEST(test_lw_dist2d_seg_arc),
PG_TEST(test_lw_dist2d_arc_arc),
CU_TEST_INFO_NULL
};
CU_SuiteInfo measures_suite = {"PostGIS Measures Suite", NULL, NULL, measures_tests};
4 changes: 2 additions & 2 deletions liblwgeom/cunit/cu_tree.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -191,9 +191,9 @@ static void test_tree_circ_distance(void)
d2 = lwgeom_distance_spheroid(lwg1, lwg2, &s, threshold);
// printf("d1 = %g d2 = %g\n", d1 * WGS84_RADIUS, d2 * WGS84_RADIUS);
// printf("line\n");
circ_tree_print(c1, 0);
// circ_tree_print(c1, 0);
// printf("poly\n");
circ_tree_print(c2, 0);
// circ_tree_print(c2, 0);
circ_tree_free(c1);
circ_tree_free(c2);
lwgeom_free(lwg1);
Expand Down
7 changes: 4 additions & 3 deletions liblwgeom/g_box.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -213,7 +213,8 @@ int gbox_overlaps(const GBOX *g1, const GBOX *g2)
return LW_TRUE;
}

int gbox_overlaps_2d(const GBOX *g1, const GBOX *g2)
int
gbox_overlaps_2d(const GBOX *g1, const GBOX *g2)
{

/* Make sure our boxes are consistent */
Expand Down Expand Up @@ -323,13 +324,13 @@ size_t gbox_serialized_size(uint8_t flags)
static int lwcircle_calculate_gbox_cartesian(const POINT4D *p1, const POINT4D *p2, const POINT4D *p3, GBOX *gbox)
{
POINT2D xmin, ymin, xmax, ymax;
POINT4D center;
POINT2D center;
int p2_side;
double radius;

LWDEBUG(2, "lwcircle_calculate_gbox called.");

radius = lwcircle_center(p1, p2, p3, &center);
radius = lwcircle_center((POINT2D*)p1, (POINT2D*)p2, (POINT2D*)p3, &center);

/* Negative radius signals straight line, p1/p2/p3 are colinear */
if (radius < 0.0)
Expand Down
2 changes: 1 addition & 1 deletion liblwgeom/liblwgeom.h.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1802,7 +1802,7 @@ extern char *lwmessage_truncate(char *str, int startpos, int endpos, int maxleng
******************************************************************************/

int lwgeom_has_arc(const LWGEOM *geom);
double lwcircle_center(const POINT4D *p1, const POINT4D *p2, const POINT4D *p3, POINT4D *result);
double lwcircle_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, POINT2D *result);
LWGEOM *lwgeom_segmentize(LWGEOM *geom, uint32_t perQuad);
LWGEOM *lwgeom_desegmentize(LWGEOM *geom);

Expand Down
17 changes: 8 additions & 9 deletions liblwgeom/lwsegmentize.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -82,13 +82,13 @@ lwgeom_has_arc(const LWGEOM *geom)
* point is coincident with either end point, they are taken as colinear.
*/
double
lwcircle_center(const POINT4D *p1, const POINT4D *p2, const POINT4D *p3, POINT4D *result)
lwcircle_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, POINT2D *result)
{
POINT4D c;
POINT2D c;
double cx, cy, cr;
double temp, bc, cd, det;

c.x = c.y = c.z = c.m = 0.0;
c.x = c.y = 0.0;

LWDEBUGF(2, "lwcircle_center called (%.16f,%.16f), (%.16f,%.16f), (%.16f,%.16f).", p1->x, p1->y, p2->x, p2->y, p3->x, p3->y);

Expand Down Expand Up @@ -157,7 +157,7 @@ static double interpolate_arc(double angle, double a1, double a2, double a3, dou
static POINTARRAY *
lwcircle_segmentize(POINT4D *p1, POINT4D *p2, POINT4D *p3, uint32_t perQuad)
{
POINT4D center;
POINT2D center;
POINT4D pt;
int p2_side = 0;
int clockwise = LW_TRUE;
Expand All @@ -170,7 +170,7 @@ lwcircle_segmentize(POINT4D *p1, POINT4D *p2, POINT4D *p3, uint32_t perQuad)

LWDEBUG(2, "lwcircle_calculate_gbox called.");

radius = lwcircle_center(p1, p2, p3, &center);
radius = lwcircle_center((POINT2D*)p1, (POINT2D*)p2, (POINT2D*)p3, &center);
p2_side = signum(lw_segment_side((POINT2D*)p1, (POINT2D*)p3, (POINT2D*)p2));

/* Matched start/end points imply circle */
Expand Down Expand Up @@ -527,16 +527,15 @@ lwgeom_segmentize(LWGEOM *geom, uint32_t perQuad)
*/
static int pt_continues_arc(const POINT4D *a1, const POINT4D *a2, const POINT4D *a3, const POINT4D *b)
{
POINT4D center;
POINT4D *centerptr=&center;
double radius = lwcircle_center(a1, a2, a3, &center);
POINT2D center;
double radius = lwcircle_center((POINT2D*)a1, (POINT2D*)a2, (POINT2D*)a3, &center);
double b_distance, diff;

/* Co-linear a1/a2/a3 */
if ( radius < 0.0 )
return LW_FALSE;

b_distance = distance2d_pt_pt((POINT2D*)b, (POINT2D*)centerptr);
b_distance = distance2d_pt_pt((POINT2D*)b, &center);
diff = fabs(radius - b_distance);
LWDEBUGF(4, "circle_radius=%g, b_distance=%g, diff=%g, percentage=%g", radius, b_distance, diff, diff/radius);

Expand Down
Loading

0 comments on commit 503d866

Please sign in to comment.