From e21896ff317e3c28d5c44b6ad00d425953bfc0bf Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Fri, 27 Jan 2012 20:55:22 +0000 Subject: [PATCH] Tighten up on-arc test a bit more. git-svn-id: http://svn.osgeo.org/postgis/trunk@8948 b70326c6-7e19-0410-871a-916f4a2858ee --- liblwgeom/cunit/cu_ptarray.c | 13 ++ liblwgeom/lwsegmentize.c | 397 +---------------------------------- 2 files changed, 20 insertions(+), 390 deletions(-) diff --git a/liblwgeom/cunit/cu_ptarray.c b/liblwgeom/cunit/cu_ptarray.c index 10d974d9d..e25484b53 100644 --- a/liblwgeom/cunit/cu_ptarray.c +++ b/liblwgeom/cunit/cu_ptarray.c @@ -262,6 +262,19 @@ static void test_ptarray_desegmentize() { LWGEOM *in, *out; char *str; + + /* It would be nice if this example returned two arcs (it's the intersection of two circles) + but it looks like the intersection itself is too sloppy in generating the derived point + to accurately reconstruct the circles. + in = lwgeom_from_text("POLYGON((0.5 0,0.471177920604846 -0.292635483024192,0.38581929876693 -0.574025148547634,0.247204418453818 -0.833355349529403,0.0606601717798223 -1.06066017177982,-5.44089437167602e-17 -1.11044268820754,-0.0606601717798188 -1.06066017177982,-0.247204418453816 -0.833355349529406,-0.385819298766929 -0.574025148547639,-0.471177920604845 -0.292635483024197,-0.5 -4.84663372329776e-15,-0.471177920604847 0.292635483024187,-0.385819298766932 0.57402514854763,-0.247204418453821 0.833355349529398,-0.0606601717798256 1.06066017177982,3.45538806345173e-16 1.11044268820754,0.0606601717798183 1.06066017177982,0.247204418453816 0.833355349529407,0.385819298766929 0.574025148547638,0.471177920604845 0.292635483024196,0.5 0))"); + out = lwgeom_desegmentize(in); + str = lwgeom_to_wkt(out, WKT_ISO, 8, NULL); + printf("%s\n", str); + CU_ASSERT_STRING_EQUAL(str, "CIRCULARSTRING(-1 0,0 1,0 -1)"); + lwgeom_free(in); + lwgeom_free(out); + lwfree(str); + */ in = lwgeom_from_text("LINESTRING(-1 0, 0 1, 1 0, 0 -1)"); out = lwgeom_desegmentize(in); diff --git a/liblwgeom/lwsegmentize.c b/liblwgeom/lwsegmentize.c index 08cd45a7b..f8a70cbf8 100644 --- a/liblwgeom/lwsegmentize.c +++ b/liblwgeom/lwsegmentize.c @@ -22,7 +22,7 @@ LWMLINE *lwmcurve_segmentize(LWMCURVE *mcurve, uint32_t perQuad); LWMPOLY *lwmsurface_segmentize(LWMSURFACE *msurface, uint32_t perQuad); LWCOLLECTION *lwcollection_segmentize(LWCOLLECTION *collection, uint32_t perQuad); -LWGEOM *append_segment(LWGEOM *geom, POINTARRAY *pts, int type, int srid); + LWGEOM *pta_desegmentize(POINTARRAY *points, int type, int srid); LWGEOM *lwline_desegmentize(LWLINE *line); LWGEOM *lwpolygon_desegmentize(LWPOLY *poly); @@ -31,7 +31,6 @@ LWGEOM *lwmpolygon_desegmentize(LWMPOLY *mpoly); LWGEOM *lwgeom_desegmentize(LWGEOM *geom); - /* * Tolerance used to determine equality. */ @@ -521,181 +520,6 @@ lwgeom_segmentize(LWGEOM *geom, uint32_t perQuad) return ogeom; } -/******************************************************************************* - * End curve segmentize functions - ******************************************************************************/ -LWGEOM * -append_segment(LWGEOM *geom, POINTARRAY *pts, int type, int srid) -{ - LWGEOM *result; - int currentType, i; - - LWDEBUGF(2, "append_segment called %p, %p, %d, %d", geom, pts, type, srid); - - if (geom == NULL) - { - if (type == LINETYPE) - { - LWDEBUG(3, "append_segment: line to NULL"); - - return (LWGEOM *)lwline_construct(srid, NULL, pts); - } - else if (type == CIRCSTRINGTYPE) - { -#if POSTGIS_DEBUG_LEVEL >= 4 - POINT4D tmp; - - LWDEBUGF(4, "append_segment: curve to NULL %d", pts->npoints); - - for (i=0; inpoints; i++) - { - getPoint4d_p(pts, i, &tmp); - LWDEBUGF(4, "new point: (%.16f,%.16f)",tmp.x,tmp.y); - } -#endif - return (LWGEOM *)lwcircstring_construct(srid, NULL, pts); - } - else - { - lwerror("Invalid segment type %d - %s.", type, lwtype_name(type)); - } - } - - currentType = geom->type; - if (currentType == LINETYPE && type == LINETYPE) - { - POINTARRAY *newPoints; - POINT4D pt; - LWLINE *line = (LWLINE *)geom; - - LWDEBUG(3, "append_segment: line to line"); - - newPoints = ptarray_construct(FLAGS_GET_Z(pts->flags), FLAGS_GET_M(pts->flags), pts->npoints + line->points->npoints - 1); - for (i=0; ipoints->npoints; i++) - { - getPoint4d_p(line->points, i, &pt); - LWDEBUGF(5, "copying to %p [%d]", newPoints, i); - ptarray_set_point4d(newPoints, i, &pt); - } - for (i=1; inpoints; i++) - { - getPoint4d_p(pts, i, &pt); - ptarray_set_point4d(newPoints, i + line->points->npoints - 1, &pt); - } - result = (LWGEOM *)lwline_construct(srid, NULL, newPoints); - lwgeom_release(geom); - return result; - } - else if (currentType == CIRCSTRINGTYPE && type == CIRCSTRINGTYPE) - { - POINTARRAY *newPoints; - POINT4D pt; - LWCIRCSTRING *curve = (LWCIRCSTRING *)geom; - - LWDEBUG(3, "append_segment: circularstring to circularstring"); - - newPoints = ptarray_construct(FLAGS_GET_Z(pts->flags), FLAGS_GET_M(pts->flags), pts->npoints + curve->points->npoints - 1); - - LWDEBUGF(3, "New array length: %d", pts->npoints + curve->points->npoints - 1); - - for (i=0; ipoints->npoints; i++) - { - getPoint4d_p(curve->points, i, &pt); - - LWDEBUGF(3, "orig point %d: (%.16f,%.16f)", i, pt.x, pt.y); - - ptarray_set_point4d(newPoints, i, &pt); - } - for (i=1; inpoints; i++) - { - getPoint4d_p(pts, i, &pt); - - LWDEBUGF(3, "new point %d: (%.16f,%.16f)", i + curve->points->npoints - 1, pt.x, pt.y); - - ptarray_set_point4d(newPoints, i + curve->points->npoints - 1, &pt); - } - result = (LWGEOM *)lwcircstring_construct(srid, NULL, newPoints); - lwgeom_release(geom); - return result; - } - else if (currentType == CIRCSTRINGTYPE && type == LINETYPE) - { - LWLINE *line; - LWGEOM **geomArray; - - LWDEBUG(3, "append_segment: line to curve"); - - geomArray = lwalloc(sizeof(LWGEOM *)*2); - geomArray[0] = lwgeom_clone(geom); - - line = lwline_construct(srid, NULL, pts); - geomArray[1] = lwgeom_clone((LWGEOM *)line); - - result = (LWGEOM *)lwcollection_construct(COMPOUNDTYPE, srid, NULL, 2, geomArray); - lwfree((LWGEOM *)line); - lwgeom_release(geom); - return result; - } - else if (currentType == LINETYPE && type == CIRCSTRINGTYPE) - { - LWCIRCSTRING *curve; - LWGEOM **geomArray; - - LWDEBUG(3, "append_segment: circularstring to line"); - - geomArray = lwalloc(sizeof(LWGEOM *)*2); - geomArray[0] = lwgeom_clone(geom); - - curve = lwcircstring_construct(srid, NULL, pts); - geomArray[1] = lwgeom_clone((LWGEOM *)curve); - - result = (LWGEOM *)lwcollection_construct(COMPOUNDTYPE, srid, NULL, 2, geomArray); - lwfree((LWGEOM *)curve); - lwgeom_release(geom); - return result; - } - else if (currentType == COMPOUNDTYPE) - { - LWGEOM *newGeom; - LWCOMPOUND *compound; - int count; - LWGEOM **geomArray; - - compound = (LWCOMPOUND *)geom; - count = compound->ngeoms + 1; - geomArray = lwalloc(sizeof(LWGEOM *)*count); - for (i=0; ingeoms; i++) - { - geomArray[i] = lwgeom_clone(compound->geoms[i]); - } - if (type == LINETYPE) - { - LWDEBUG(3, "append_segment: line to compound"); - - newGeom = (LWGEOM *)lwline_construct(srid, NULL, pts); - } - else if (type == CIRCSTRINGTYPE) - { - LWDEBUG(3, "append_segment: circularstring to compound"); - - newGeom = (LWGEOM *)lwcircstring_construct(srid, NULL, pts); - } - else - { - lwerror("Invalid segment type %d - %s.", type, lwtype_name(type)); - return NULL; - } - geomArray[compound->ngeoms] = lwgeom_clone(newGeom); - - result = (LWGEOM *)lwcollection_construct(COMPOUNDTYPE, srid, NULL, count, geomArray); - lwfree(newGeom); - lwgeom_release(geom); - return result; - } - lwerror("Invalid state [%d] %s - [%d] %s", - currentType, lwtype_name(currentType), type, lwtype_name(type)); - return NULL; -} /** * Returns LW_TRUE if b is on the arc formed by a1/a2/a3, but not within @@ -706,14 +530,18 @@ static int pt_continues_arc(const POINT4D *a1, const POINT4D *a2, const POINT4D POINT4D center; POINT4D *centerptr=¢er; double radius = lwcircle_center(a1, a2, a3, ¢er); + 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); + diff = fabs(radius - b_distance); + LWDEBUGF(4, "circle_radius=%g, b_distance=%g, diff=%g, percentage=%g", radius, b_distance, diff, diff/radius); + /* Is the point b on the circle? */ - if ( fabs(radius - distance2d_pt_pt((POINT2D*)b, (POINT2D*)centerptr)) < EPSILON_SQLMM ) -// if ( FP_EQUALS(radius, distance2d_pt_pt((POINT2D*)b, (POINT2D*)centerptr)) ) + if ( diff < EPSILON_SQLMM ) { int a2_side = signum(lw_segment_side((POINT2D*)a1, (POINT2D*)a3, (POINT2D*)a2)); int b_side = signum(lw_segment_side((POINT2D*)a1, (POINT2D*)a3, (POINT2D*)b)); @@ -726,8 +554,6 @@ static int pt_continues_arc(const POINT4D *a1, const POINT4D *a2, const POINT4D return LW_FALSE; } -LWGEOM* pta_desegmentize2(POINTARRAY *points, int type, int srid); - static LWGEOM* linestring_from_pa(const POINTARRAY *pa, int srid, int start, int end) { @@ -888,215 +714,6 @@ pta_desegmentize(POINTARRAY *points, int type, int srid) return lwcollection_as_lwgeom(outcol); } -LWGEOM * -pta_desegmentize2(POINTARRAY *points, int type, int srid) -{ - int i, j, commit, isline, count; - double last_angle, last_length; - double dxab, dyab, dxbc, dybc, theta, length; - POINT4D a, b, c, tmp; - POINTARRAY *pts; - LWGEOM *geom = NULL; - - LWDEBUG(2, "pta_desegmentize called."); - - getPoint4d_p(points, 0, &a); - getPoint4d_p(points, 1, &b); - getPoint4d_p(points, 2, &c); - - dxab = b.x - a.x; - dyab = b.y - a.y; - dxbc = c.x - b.x; - dybc = c.y - b.y; - - theta = atan2(dyab, dxab); - last_angle = theta - atan2(dybc, dxbc); - last_length = sqrt(dxbc*dxbc+dybc*dybc); - length = sqrt(dxab*dxab+dyab*dyab); - if ((last_length - length) < EPSILON_SQLMM) - { - isline = -1; - LWDEBUG(3, "Starting as unknown."); - } - else - { - isline = 1; - LWDEBUG(3, "Starting as line."); - } - - commit = 0; - for (i=3; inpoints; i++) - { - getPoint4d_p(points, i-2, &a); - getPoint4d_p(points, i-1, &b); - getPoint4d_p(points, i, &c); - - dxab = b.x - a.x; - dyab = b.y - a.y; - dxbc = c.x - b.x; - dybc = c.y - b.y; - - LWDEBUGF(3, "(dxab, dyab, dxbc, dybc) (%.16f, %.16f, %.16f, %.16f)", dxab, dyab, dxbc, dybc); - - theta = atan2(dyab, dxab); - theta = theta - atan2(dybc, dxbc); - length = sqrt(dxbc*dxbc+dybc*dybc); - - LWDEBUGF(3, "Last/current length and angle %.16f/%.16f, %.16f/%.16f", last_angle, theta, last_length, length); - - /* Found a line segment */ - if (fabs(length - last_length) > EPSILON_SQLMM || - fabs(theta - last_angle) > EPSILON_SQLMM) - { - last_length = length; - last_angle = theta; - /* We are tracking a line, keep going */ - if (isline > 0) - {} - /* We were tracking a circularstring, commit it and start line*/ - else if (isline == 0) - { - LWDEBUGF(3, "Building circularstring, %d - %d", commit, i); - - count = i - commit; - pts = ptarray_construct( - FLAGS_GET_Z(points->flags), - FLAGS_GET_M(points->flags), - 3); - getPoint4d_p(points, commit, &tmp); - ptarray_set_point4d(pts, 0, &tmp); - getPoint4d_p(points, - commit + count/2, &tmp); - ptarray_set_point4d(pts, 1, &tmp); - getPoint4d_p(points, i - 1, &tmp); - ptarray_set_point4d(pts, 2, &tmp); - - commit = i-1; - geom = append_segment(geom, pts, CIRCSTRINGTYPE, srid); - isline = -1; - - /* - * We now need to move ahead one point to - * determine if it's a potential new curve, - * since the last_angle value is corrupt. - * - * Note we can only look ahead one point if - * we are not already at the end of our - * set of points. - */ - if (i < points->npoints - 1) - { - i++; - - getPoint4d_p(points, i-2, &a); - getPoint4d_p(points, i-1, &b); - getPoint4d_p(points, i, &c); - - dxab = b.x - a.x; - dyab = b.y - a.y; - dxbc = c.x - b.x; - dybc = c.y - b.y; - - theta = atan2(dyab, dxab); - last_angle = theta - atan2(dybc, dxbc); - last_length = sqrt(dxbc*dxbc+dybc*dybc); - length = sqrt(dxab*dxab+dyab*dyab); - if ((last_length - length) < EPSILON_SQLMM) - { - isline = -1; - LWDEBUG(3, "Restarting as unknown."); - } - else - { - isline = 1; - LWDEBUG(3, "Restarting as line."); - } - } - else - { - /* Indicate that we have tracked a CIRCSTRING */ - isline = 0; - } - } - /* We didn't know what we were tracking, now we do. */ - else - { - LWDEBUG(3, "It's a line"); - isline = 1; - } - } - /* Found a circularstring segment */ - else - { - /* We were tracking a circularstring, commit it and start line */ - if (isline > 0) - { - LWDEBUGF(3, "Building line, %d - %d", commit, i-2); - - count = i - commit - 2; - - pts = ptarray_construct( - FLAGS_GET_Z(points->flags), - FLAGS_GET_M(points->flags), - count); - for (j=commit; j 2) - { - LWDEBUGF(3, "Finishing circularstring %d,%d.", commit, i); - - pts = ptarray_construct( - FLAGS_GET_Z(points->flags), - FLAGS_GET_M(points->flags), - 3); - getPoint4d_p(points, commit, &tmp); - ptarray_set_point4d(pts, 0, &tmp); - getPoint4d_p(points, commit + count/2, &tmp); - ptarray_set_point4d(pts, 1, &tmp); - getPoint4d_p(points, i - 1, &tmp); - ptarray_set_point4d(pts, 2, &tmp); - - geom = append_segment(geom, pts, CIRCSTRINGTYPE, srid); - } - else - { - LWDEBUGF(3, "Finishing line %d,%d.", commit, i); - - pts = ptarray_construct( - FLAGS_GET_Z(points->flags), - FLAGS_GET_M(points->flags), - count); - for (j=commit; j