From: Paul Ramsey Date: Tue, 30 Jun 2015 13:17:34 +0000 (+0000) Subject: reformat to style guide X-Git-Tag: 2.2.0rc1~280 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=b98f69a5381ec817932c006f3e5e068ace771d0d;p=postgis reformat to style guide git-svn-id: http://svn.osgeo.org/postgis/trunk@13758 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/liblwgeom/lwlinearreferencing.c b/liblwgeom/lwlinearreferencing.c index 3ff32246d..0476ed4df 100644 --- a/liblwgeom/lwlinearreferencing.c +++ b/liblwgeom/lwlinearreferencing.c @@ -15,7 +15,7 @@ #include "lwgeom_log.h" #include "measures3d.h" -static int +static int segment_locate_along(const POINT4D *p1, const POINT4D *p2, double m, double offset, POINT4D *pn) { double m1 = p1->m; @@ -28,7 +28,7 @@ segment_locate_along(const POINT4D *p1, const POINT4D *p2, double m, double offs return LW_FALSE; } - if( m1 == m2 ) + if( m1 == m2 ) { /* Degenerate case: same M on both points. If they are the same point we just return one of them. */ @@ -58,7 +58,7 @@ segment_locate_along(const POINT4D *p1, const POINT4D *p2, double m, double offs pn->y += cos(theta) * offset; } - return LW_TRUE; + return LW_TRUE; } @@ -68,10 +68,10 @@ ptarray_locate_along(const POINTARRAY *pa, double m, double offset) int i; POINT4D p1, p2, pn; POINTARRAY *dpa = NULL; - + /* Can't do anything with degenerate point arrays */ if ( ! pa || pa->npoints < 2 ) return NULL; - + /* Walk through each segment in the point array */ for ( i = 1; i < pa->npoints; i++ ) { @@ -88,7 +88,7 @@ ptarray_locate_along(const POINTARRAY *pa, double m, double offset) /* Add our new point to the array */ ptarray_append_point(dpa, &pn, 0); - } + } return dpa; } @@ -103,7 +103,7 @@ lwline_locate_along(const LWLINE *lwline, double m, double offset) /* Return degenerates upwards */ if ( ! lwline ) return NULL; - + /* Create empty return shell */ srid = lwgeom_get_srid(lwg); hasz = lwgeom_has_z(lwg); @@ -114,7 +114,7 @@ lwline_locate_along(const LWLINE *lwline, double m, double offset) /* Find points along */ opa = ptarray_locate_along(lwline->points, m, offset); } - else + else { LWLINE *lwline_measured = lwline_measured_from_lwline(lwline, 0.0, 1.0); opa = ptarray_locate_along(lwline_measured->points, m, offset); @@ -122,9 +122,9 @@ lwline_locate_along(const LWLINE *lwline, double m, double offset) } /* Return NULL as EMPTY */ - if ( ! opa ) + if ( ! opa ) return lwmpoint_construct_empty(srid, hasz, hasm); - + /* Convert pointarray into a multipoint */ mp = lwmpoint_construct(srid, opa); ptarray_free(opa); @@ -150,15 +150,16 @@ lwmline_locate_along(const LWMLINE *lwmline, double m, double offset) LWMPOINT *along = lwline_locate_along(lwmline->geoms[i], m, offset); if ( along ) { - if ( ! lwgeom_is_empty((LWGEOM*)along) ) + if ( ! lwgeom_is_empty((LWGEOM*)along) ) { - for ( j = 0; j < along->ngeoms; j++ ) + for ( j = 0; j < along->ngeoms; j++ ) { lwmpoint_add_lwpoint(lwmpoint, along->geoms[j]); } } /* Free the containing geometry, but leave the sub-geometries around */ - along->ngeoms = 0; lwmpoint_free(along); + along->ngeoms = 0; + lwmpoint_free(along); } } return lwmpoint; @@ -183,7 +184,7 @@ lwmpoint_locate_along(const LWMPOINT *lwin, double m, double offset) LWGEOM *lwg = lwmpoint_as_lwgeom(lwin); LWMPOINT *lwout = NULL; int i; - + /* Construct return */ lwout = lwmpoint_construct_empty(lwgeom_get_srid(lwg), lwgeom_has_z(lwg), lwgeom_has_m(lwg)); @@ -193,7 +194,7 @@ lwmpoint_locate_along(const LWMPOINT *lwin, double m, double offset) if ( FP_EQUALS(m, point_m) ) { lwmpoint_add_lwpoint(lwout, lwpoint_clone(lwin->geoms[i])); - } + } } return lwout; @@ -205,8 +206,8 @@ lwgeom_locate_along(const LWGEOM *lwin, double m, double offset) if ( ! lwin ) return NULL; if ( ! lwgeom_has_m(lwin) ) - lwerror("Input geometry does not have a measure dimension"); - + lwerror("Input geometry does not have a measure dimension"); + switch (lwin->type) { case POINTTYPE: @@ -262,7 +263,7 @@ double lwpoint_get_ordinate(const POINT4D *p, char ordinate) } -/** +/** * Given a point, ordinate number and value, set that ordinate on the * point. */ @@ -319,7 +320,7 @@ int point_interpolate(const POINT4D *p1, const POINT4D *p2, POINT4D *p, int hasz } if ( FP_MIN(p1_value, p2_value) > interpolation_value || - FP_MAX(p1_value, p2_value) < interpolation_value ) + FP_MAX(p1_value, p2_value) < interpolation_value ) { lwerror("Cannot interpolate to a value (%g) not between the input points (%g, %g).", interpolation_value, p1_value, p2_value); return 0; @@ -369,19 +370,19 @@ lwpoint_clip_to_ordinate_range(const LWPOINT *point, char ordinate, double from, /* Read Z/M info */ hasz = lwgeom_has_z(lwpoint_as_lwgeom(point)); hasm = lwgeom_has_m(lwpoint_as_lwgeom(point)); - + /* Prepare return object */ lwgeom_out = lwcollection_construct_empty(MULTIPOINTTYPE, point->srid, hasz, hasm); /* Test if ordinate is in range */ - lwpoint_getPoint4d_p(point, &p4d); + lwpoint_getPoint4d_p(point, &p4d); ordinate_value = lwpoint_get_ordinate(&p4d, ordinate); if ( from <= ordinate_value && to >= ordinate_value ) { LWPOINT *lwp = lwpoint_clone(point); lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(lwp)); } - + /* Set the bbox, if necessary */ if ( lwgeom_out->bbox ) { @@ -419,7 +420,7 @@ lwmpoint_clip_to_ordinate_range(const LWMPOINT *mpoint, char ordinate, double fr /* Read Z/M info */ hasz = lwgeom_has_z(lwmpoint_as_lwgeom(mpoint)); hasm = lwgeom_has_m(lwmpoint_as_lwgeom(mpoint)); - + /* Prepare return object */ lwgeom_out = lwcollection_construct_empty(MULTIPOINTTYPE, mpoint->srid, hasz, hasm); @@ -428,17 +429,17 @@ lwmpoint_clip_to_ordinate_range(const LWMPOINT *mpoint, char ordinate, double fr { POINT4D p4d; double ordinate_value; - + lwpoint_getPoint4d_p(mpoint->geoms[i], &p4d); ordinate_value = lwpoint_get_ordinate(&p4d, ordinate); - + if ( from <= ordinate_value && to >= ordinate_value ) { LWPOINT *lwp = lwpoint_clone(mpoint->geoms[i]); lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(lwp)); } } - + /* Set the bbox, if necessary */ if ( lwgeom_out->bbox ) { @@ -516,7 +517,7 @@ lwmline_clip_to_ordinate_range(const LWMLINE *mline, char ordinate, double from, lwgeom_drop_bbox((LWGEOM*)lwgeom_out); lwgeom_add_bbox((LWGEOM*)lwgeom_out); } - + if ( ! homogeneous ) { lwgeom_out->type = COLLECTIONTYPE; @@ -707,7 +708,7 @@ lwline_clip_to_ordinate_range(const LWLINE *line, char ordinate, double from, do LWPOINT *opoint = lwpoint_construct(line->srid, NULL, dp); lwgeom_out->type = COLLECTIONTYPE; lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(opoint)); - + } else { @@ -766,41 +767,41 @@ lwgeom_clip_to_ordinate_range(const LWGEOM *lwin, char ordinate, double from, do LWCOLLECTION *out_col; LWCOLLECTION *out_offset; int i; - + if ( ! lwin ) lwerror("lwgeom_clip_to_ordinate_range: null input geometry!"); - + switch ( lwin->type ) { - case LINETYPE: - out_col = lwline_clip_to_ordinate_range((LWLINE*)lwin, ordinate, from, to); - break; - case MULTILINETYPE: - out_col = lwmline_clip_to_ordinate_range((LWMLINE*)lwin, ordinate, from, to); - break; - case MULTIPOINTTYPE: - out_col = lwmpoint_clip_to_ordinate_range((LWMPOINT*)lwin, ordinate, from, to); - break; - case POINTTYPE: - out_col = lwpoint_clip_to_ordinate_range((LWPOINT*)lwin, ordinate, from, to); - break; - default: - lwerror("This function does not accept %s geometries.", lwtype_name(lwin->type)); - return NULL;; + case LINETYPE: + out_col = lwline_clip_to_ordinate_range((LWLINE*)lwin, ordinate, from, to); + break; + case MULTILINETYPE: + out_col = lwmline_clip_to_ordinate_range((LWMLINE*)lwin, ordinate, from, to); + break; + case MULTIPOINTTYPE: + out_col = lwmpoint_clip_to_ordinate_range((LWMPOINT*)lwin, ordinate, from, to); + break; + case POINTTYPE: + out_col = lwpoint_clip_to_ordinate_range((LWPOINT*)lwin, ordinate, from, to); + break; + default: + lwerror("This function does not accept %s geometries.", lwtype_name(lwin->type)); + return NULL;; } - + /* Stop if result is NULL */ if ( out_col == NULL ) lwerror("lwgeom_clip_to_ordinate_range clipping routine returned NULL"); - + /* Return if we aren't going to offset the result */ if ( FP_EQUALS(offset, 0.0) || lwgeom_is_empty(lwcollection_as_lwgeom(out_col)) ) return out_col; - + /* Construct a collection to hold our outputs. */ /* Things get ugly: GEOS offset drops Z's and M's so we have to drop ours */ out_offset = lwcollection_construct_empty(MULTILINETYPE, lwin->srid, 0, 0); - + /* Try and offset the linear portions of the return value */ for ( i = 0; i < out_col->ngeoms; i++ ) { @@ -820,7 +821,7 @@ lwgeom_clip_to_ordinate_range(const LWGEOM *lwin, char ordinate, double from, do } lwcollection_add_lwgeom(out_offset, lwoff); } - else + else { lwerror("lwgeom_clip_to_ordinate_range found an unexpected type (%s) in the offset routine",lwtype_name(type)); } @@ -833,45 +834,45 @@ LWCOLLECTION* lwgeom_locate_between(const LWGEOM *lwin, double from, double to, double offset) { if ( ! lwgeom_has_m(lwin) ) - lwerror("Input geometry does not have a measure dimension"); + lwerror("Input geometry does not have a measure dimension"); return lwgeom_clip_to_ordinate_range(lwin, 'M', from, to, offset); -} +} double lwgeom_interpolate_point(const LWGEOM *lwin, const LWPOINT *lwpt) { POINT4D p, p_proj; double ret = 0.0; - + if ( ! lwin ) lwerror("lwgeom_interpolate_point: null input geometry!"); if ( ! lwgeom_has_m(lwin) ) - lwerror("Input geometry does not have a measure dimension"); - + lwerror("Input geometry does not have a measure dimension"); + if ( lwgeom_is_empty(lwin) || lwpoint_is_empty(lwpt) ) - lwerror("Input geometry is empty"); - + lwerror("Input geometry is empty"); + switch ( lwin->type ) { - case LINETYPE: - { - LWLINE *lwline = lwgeom_as_lwline(lwin); - lwpoint_getPoint4d_p(lwpt, &p); - ret = ptarray_locate_point(lwline->points, &p, NULL, &p_proj); - ret = p_proj.m; - break; - } - default: - lwerror("This function does not accept %s geometries.", lwtype_name(lwin->type)); + case LINETYPE: + { + LWLINE *lwline = lwgeom_as_lwline(lwin); + lwpoint_getPoint4d_p(lwpt, &p); + ret = ptarray_locate_point(lwline->points, &p, NULL, &p_proj); + ret = p_proj.m; + break; + } + default: + lwerror("This function does not accept %s geometries.", lwtype_name(lwin->type)); } return ret; } /* * Time of closest point of approach - * + * * Given two vectors (p1-p2 and q1-q2) and * a time range (t1-t2) return the time in which * a point p is closest to a point q on their @@ -895,111 +896,114 @@ lwgeom_interpolate_point(const LWGEOM *lwin, const LWPOINT *lwpt) */ static double segments_tcpa(POINT4D* p0, const POINT4D* p1, - POINT4D* q0, const POINT4D* q1, - double t0, double t1) + POINT4D* q0, const POINT4D* q1, + double t0, double t1) { - POINT3DZ pv; /* velocity of p, aka u */ - POINT3DZ qv; /* velocity of q, aka v */ - POINT3DZ dv; /* velocity difference */ - POINT3DZ w0; /* vector between first points */ - -/* - lwnotice("FROM %g,%g,%g,%g -- %g,%g,%g,%g", - p0->x, p0->y, p0->z, p0->m, - p1->x, p1->y, p1->z, p1->m); - lwnotice(" TO %g,%g,%g,%g -- %g,%g,%g,%g", - q0->x, q0->y, q0->z, q0->m, - q1->x, q1->y, q1->z, q1->m); -*/ - - /* PV aka U */ - pv.x = ( p1->x - p0->x ); - pv.y = ( p1->y - p0->y ); - pv.z = ( p1->z - p0->z ); - /*lwnotice("PV: %g, %g, %g", pv.x, pv.y, pv.z);*/ - - /* QV aka V */ - qv.x = ( q1->x - q0->x ); - qv.y = ( q1->y - q0->y ); - qv.z = ( q1->z - q0->z ); - /*lwnotice("QV: %g, %g, %g", qv.x, qv.y, qv.z);*/ - - dv.x = pv.x - qv.x; - dv.y = pv.y - qv.y; - dv.z = pv.z - qv.z; - /*lwnotice("DV: %g, %g, %g", dv.x, dv.y, dv.z);*/ - - double dv2 = DOT(dv,dv); - /*lwnotice("DOT: %g", dv2);*/ - - if ( dv2 == 0.0 ) - { - /* Distance is the same at any time, we pick the earliest */ - return t0; - } + POINT3DZ pv; /* velocity of p, aka u */ + POINT3DZ qv; /* velocity of q, aka v */ + POINT3DZ dv; /* velocity difference */ + POINT3DZ w0; /* vector between first points */ + + /* + lwnotice("FROM %g,%g,%g,%g -- %g,%g,%g,%g", + p0->x, p0->y, p0->z, p0->m, + p1->x, p1->y, p1->z, p1->m); + lwnotice(" TO %g,%g,%g,%g -- %g,%g,%g,%g", + q0->x, q0->y, q0->z, q0->m, + q1->x, q1->y, q1->z, q1->m); + */ + + /* PV aka U */ + pv.x = ( p1->x - p0->x ); + pv.y = ( p1->y - p0->y ); + pv.z = ( p1->z - p0->z ); + /*lwnotice("PV: %g, %g, %g", pv.x, pv.y, pv.z);*/ + + /* QV aka V */ + qv.x = ( q1->x - q0->x ); + qv.y = ( q1->y - q0->y ); + qv.z = ( q1->z - q0->z ); + /*lwnotice("QV: %g, %g, %g", qv.x, qv.y, qv.z);*/ + + dv.x = pv.x - qv.x; + dv.y = pv.y - qv.y; + dv.z = pv.z - qv.z; + /*lwnotice("DV: %g, %g, %g", dv.x, dv.y, dv.z);*/ + + double dv2 = DOT(dv,dv); + /*lwnotice("DOT: %g", dv2);*/ + + if ( dv2 == 0.0 ) + { + /* Distance is the same at any time, we pick the earliest */ + return t0; + } - /* Distance at any given time, with t0 */ - w0.x = ( p0->x - q0->x ); - w0.y = ( p0->y - q0->y ); - w0.z = ( p0->z - q0->z ); + /* Distance at any given time, with t0 */ + w0.x = ( p0->x - q0->x ); + w0.y = ( p0->y - q0->y ); + w0.z = ( p0->z - q0->z ); - /*lwnotice("W0: %g, %g, %g", w0.x, w0.y, w0.z);*/ + /*lwnotice("W0: %g, %g, %g", w0.x, w0.y, w0.z);*/ - /* Check that at distance dt w0 is distance */ + /* Check that at distance dt w0 is distance */ - /* This is the fraction of measure difference */ - double t = -DOT(w0,dv) / dv2; - /*lwnotice("CLOSEST TIME (fraction): %g", t);*/ + /* This is the fraction of measure difference */ + double t = -DOT(w0,dv) / dv2; + /*lwnotice("CLOSEST TIME (fraction): %g", t);*/ - if ( t > 1.0 ) { - /* Getting closer as we move to the end */ - /*lwnotice("Converging");*/ - t = 1; - } else if ( t < 0.0 ) { - /*lwnotice("Diverging");*/ - t = 0; - } + if ( t > 1.0 ) + { + /* Getting closer as we move to the end */ + /*lwnotice("Converging");*/ + t = 1; + } + else if ( t < 0.0 ) + { + /*lwnotice("Diverging");*/ + t = 0; + } - /* Interpolate the actual points now */ + /* Interpolate the actual points now */ - p0->x += pv.x * t; - p0->y += pv.y * t; - p0->z += pv.z * t; + p0->x += pv.x * t; + p0->y += pv.y * t; + p0->z += pv.z * t; - q0->x += qv.x * t; - q0->y += qv.y * t; - q0->z += qv.z * t; + q0->x += qv.x * t; + q0->y += qv.y * t; + q0->z += qv.z * t; - t = t0 + (t1 - t0) * t; - /*lwnotice("CLOSEST TIME (real): %g", t);*/ + t = t0 + (t1 - t0) * t; + /*lwnotice("CLOSEST TIME (real): %g", t);*/ - return t; + return t; } static int ptarray_collect_mvals(const POINTARRAY *pa, double tmin, double tmax, double *mvals) { - POINT4D pbuf; - int i, n=0; - for (i=0; inpoints; ++i) - { - getPoint4d_p(pa, i, &pbuf); /* could be optimized */ - if ( pbuf.m >= tmin && pbuf.m <= tmax ) - mvals[n++] = pbuf.m; - } - return n; + POINT4D pbuf; + int i, n=0; + for (i=0; inpoints; ++i) + { + getPoint4d_p(pa, i, &pbuf); /* could be optimized */ + if ( pbuf.m >= tmin && pbuf.m <= tmax ) + mvals[n++] = pbuf.m; + } + return n; } static int compare_double(const void *pa, const void *pb) { - double a = *((double *)pa); - double b = *((double *)pb); - if ( a < b ) + double a = *((double *)pa); + double b = *((double *)pb); + if ( a < b ) return -1; - else if ( a > b ) + else if ( a > b ) return 1; - else + else return 0; } @@ -1007,17 +1011,17 @@ compare_double(const void *pa, const void *pb) static int uniq(double *vals, int nvals) { - int i, last=0; - for (i=1; inpoints; i++ ) { getPoint4d_p(pa, i, &p2); @@ -1049,8 +1053,8 @@ ptarray_locate_along_linear(const POINTARRAY *pa, double m, POINT4D *p, int from if ( segment_locate_along(&p1, &p2, m, 0, p) == LW_TRUE ) return i-1; /* found */ - p1 = p2; - } + p1 = p2; + } return -1; /* not found */ } @@ -1058,309 +1062,322 @@ ptarray_locate_along_linear(const POINTARRAY *pa, double m, POINT4D *p, int from double lwgeom_tcpa(const LWGEOM *g1, const LWGEOM *g2, double *mindist) { - LWLINE *l1, *l2; - int i; - const GBOX *gbox1, *gbox2; - double tmin, tmax; - double *mvals; - int nmvals = 0; - double mintime; - double mindist2 = FLT_MAX; /* minimum distance, squared */ - - if ( ! lwgeom_has_m(g1) || ! lwgeom_has_m(g2) ) { + LWLINE *l1, *l2; + int i; + const GBOX *gbox1, *gbox2; + double tmin, tmax; + double *mvals; + int nmvals = 0; + double mintime; + double mindist2 = FLT_MAX; /* minimum distance, squared */ + + if ( ! lwgeom_has_m(g1) || ! lwgeom_has_m(g2) ) + { lwerror("Both input geometries must have a measure dimension"); - return -1; - } + return -1; + } - l1 = lwgeom_as_lwline(g1); - l2 = lwgeom_as_lwline(g2); + l1 = lwgeom_as_lwline(g1); + l2 = lwgeom_as_lwline(g2); - if ( ! l1 || ! l2 ) { + if ( ! l1 || ! l2 ) + { lwerror("Both input geometries must be linestrings"); - return -1; - } + return -1; + } - if ( l1->points->npoints < 2 || l2->points->npoints < 2 ) { + if ( l1->points->npoints < 2 || l2->points->npoints < 2 ) + { lwerror("Both input lines must have at least 2 points"); - return -1; - } - - /* WARNING: these ranges may be wider than real ones */ - gbox1 = lwgeom_get_bbox(g1); - gbox2 = lwgeom_get_bbox(g2); - - assert(gbox1); /* or the npoints check above would have failed */ - assert(gbox2); /* or the npoints check above would have failed */ - - /* - * Find overlapping M range - * WARNING: may be larger than the real one - */ - - tmin = FP_MAX(gbox1->mmin, gbox2->mmin); - tmax = FP_MIN(gbox1->mmax, gbox2->mmax); - - if ( tmax < tmin ) { - LWDEBUG(1, "Inputs never exist at the same time"); - return -2; - } - - // lwnotice("Min:%g, Max:%g", tmin, tmax); - - /* - * Collect M values in common time range from inputs - */ - - mvals = lwalloc( sizeof(double) * - ( l1->points->npoints + l2->points->npoints ) ); - - /* TODO: also clip the lines ? */ - nmvals = ptarray_collect_mvals(l1->points, tmin, tmax, mvals); - nmvals += ptarray_collect_mvals(l2->points, tmin, tmax, mvals + nmvals); - - /* Sort values in ascending order */ - qsort(mvals, nmvals, sizeof(double), compare_double); - - /* Remove duplicated values */ - nmvals = uniq(mvals, nmvals); - - if ( nmvals < 2 ) - {{ - /* there's a single time, must be that one... */ - double t0 = mvals[0]; - POINT4D p0, p1; - LWDEBUGF(1, "Inputs only exist both at a single time (%g)", t0); - if ( mindist ) - { - if ( -1 == ptarray_locate_along_linear(l1->points, t0, &p0, 0) ) - { - lwfree(mvals); - lwerror("Could not find point with M=%g on first geom", t0); - return -1; - } - if ( -1 == ptarray_locate_along_linear(l2->points, t0, &p1, 0) ) - { - lwfree(mvals); - lwerror("Could not find point with M=%g on second geom", t0); - return -1; - } - *mindist = distance3d_pt_pt((POINT3D*)&p0, (POINT3D*)&p1); - } - lwfree(mvals); - return t0; - }} - - /* - * For each consecutive pair of measures, compute time of closest point - * approach and actual distance between points at that time - */ - mintime = tmin; - for (i=1; ipoints, t0, &p0, 0); - if ( -1 == seg ) continue; /* possible, if GBOX is approximated */ - // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t0, seg, p0.x, p0.y, p0.z); - - seg = ptarray_locate_along_linear(l1->points, t1, &p1, seg); - if ( -1 == seg ) continue; /* possible, if GBOX is approximated */ - // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t1, seg, p1.x, p1.y, p1.z); - - seg = ptarray_locate_along_linear(l2->points, t0, &q0, 0); - if ( -1 == seg ) continue; /* possible, if GBOX is approximated */ - // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t0, seg, q0.x, q0.y, q0.z); - - seg = ptarray_locate_along_linear(l2->points, t1, &q1, seg); - if ( -1 == seg ) continue; /* possible, if GBOX is approximated */ - // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t1, seg, q1.x, q1.y, q1.z); - - t = segments_tcpa(&p0, &p1, &q0, &q1, t0, t1); + return -1; + } + + /* WARNING: these ranges may be wider than real ones */ + gbox1 = lwgeom_get_bbox(g1); + gbox2 = lwgeom_get_bbox(g2); + + assert(gbox1); /* or the npoints check above would have failed */ + assert(gbox2); /* or the npoints check above would have failed */ + + /* + * Find overlapping M range + * WARNING: may be larger than the real one + */ + + tmin = FP_MAX(gbox1->mmin, gbox2->mmin); + tmax = FP_MIN(gbox1->mmax, gbox2->mmax); + + if ( tmax < tmin ) + { + LWDEBUG(1, "Inputs never exist at the same time"); + return -2; + } + + // lwnotice("Min:%g, Max:%g", tmin, tmax); + + /* + * Collect M values in common time range from inputs + */ + + mvals = lwalloc( sizeof(double) * + ( l1->points->npoints + l2->points->npoints ) ); + + /* TODO: also clip the lines ? */ + nmvals = ptarray_collect_mvals(l1->points, tmin, tmax, mvals); + nmvals += ptarray_collect_mvals(l2->points, tmin, tmax, mvals + nmvals); + + /* Sort values in ascending order */ + qsort(mvals, nmvals, sizeof(double), compare_double); + + /* Remove duplicated values */ + nmvals = uniq(mvals, nmvals); + + if ( nmvals < 2 ) + { + { + /* there's a single time, must be that one... */ + double t0 = mvals[0]; + POINT4D p0, p1; + LWDEBUGF(1, "Inputs only exist both at a single time (%g)", t0); + if ( mindist ) + { + if ( -1 == ptarray_locate_along_linear(l1->points, t0, &p0, 0) ) + { + lwfree(mvals); + lwerror("Could not find point with M=%g on first geom", t0); + return -1; + } + if ( -1 == ptarray_locate_along_linear(l2->points, t0, &p1, 0) ) + { + lwfree(mvals); + lwerror("Could not find point with M=%g on second geom", t0); + return -1; + } + *mindist = distance3d_pt_pt((POINT3D*)&p0, (POINT3D*)&p1); + } + lwfree(mvals); + return t0; + } + } + + /* + * For each consecutive pair of measures, compute time of closest point + * approach and actual distance between points at that time + */ + mintime = tmin; + for (i=1; ipoints, t0, &p0, 0); + if ( -1 == seg ) continue; /* possible, if GBOX is approximated */ + // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t0, seg, p0.x, p0.y, p0.z); + + seg = ptarray_locate_along_linear(l1->points, t1, &p1, seg); + if ( -1 == seg ) continue; /* possible, if GBOX is approximated */ + // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t1, seg, p1.x, p1.y, p1.z); + + seg = ptarray_locate_along_linear(l2->points, t0, &q0, 0); + if ( -1 == seg ) continue; /* possible, if GBOX is approximated */ + // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t0, seg, q0.x, q0.y, q0.z); + + seg = ptarray_locate_along_linear(l2->points, t1, &q1, seg); + if ( -1 == seg ) continue; /* possible, if GBOX is approximated */ + // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t1, seg, q1.x, q1.y, q1.z); + + t = segments_tcpa(&p0, &p1, &q0, &q1, t0, t1); /* - lwnotice("Closest points: %g,%g,%g and %g,%g,%g at time %g", - p0.x, p0.y, p0.z, - q0.x, q0.y, q0.z, t); + lwnotice("Closest points: %g,%g,%g and %g,%g,%g at time %g", + p0.x, p0.y, p0.z, + q0.x, q0.y, q0.z, t); */ - dist2 = ( q0.x - p0.x ) * ( q0.x - p0.x ) + - ( q0.y - p0.y ) * ( q0.y - p0.y ) + - ( q0.z - p0.z ) * ( q0.z - p0.z ); - if ( dist2 < mindist2 ) - { - mindist2 = dist2; - mintime = t; - // lwnotice("MINTIME: %g", mintime); - } - } - - /* - * Release memory - */ - - lwfree(mvals); - - if ( mindist ) { - *mindist = sqrt(mindist2); - } - /*lwnotice("MINDIST: %g", sqrt(mindist2));*/ - - return mintime; + dist2 = ( q0.x - p0.x ) * ( q0.x - p0.x ) + + ( q0.y - p0.y ) * ( q0.y - p0.y ) + + ( q0.z - p0.z ) * ( q0.z - p0.z ); + if ( dist2 < mindist2 ) + { + mindist2 = dist2; + mintime = t; + // lwnotice("MINTIME: %g", mintime); + } + } + + /* + * Release memory + */ + + lwfree(mvals); + + if ( mindist ) + { + *mindist = sqrt(mindist2); + } + /*lwnotice("MINDIST: %g", sqrt(mindist2));*/ + + return mintime; } int lwgeom_cpa_within(const LWGEOM *g1, const LWGEOM *g2, double maxdist) { - LWLINE *l1, *l2; - int i; - const GBOX *gbox1, *gbox2; - double tmin, tmax; - double *mvals; - int nmvals = 0; - double maxdist2 = maxdist * maxdist; - int within = LW_FALSE; - - if ( ! lwgeom_has_m(g1) || ! lwgeom_has_m(g2) ) { + LWLINE *l1, *l2; + int i; + const GBOX *gbox1, *gbox2; + double tmin, tmax; + double *mvals; + int nmvals = 0; + double maxdist2 = maxdist * maxdist; + int within = LW_FALSE; + + if ( ! lwgeom_has_m(g1) || ! lwgeom_has_m(g2) ) + { lwerror("Both input geometries must have a measure dimension"); - return LW_FALSE; - } + return LW_FALSE; + } - l1 = lwgeom_as_lwline(g1); - l2 = lwgeom_as_lwline(g2); + l1 = lwgeom_as_lwline(g1); + l2 = lwgeom_as_lwline(g2); - if ( ! l1 || ! l2 ) { + if ( ! l1 || ! l2 ) + { lwerror("Both input geometries must be linestrings"); - return LW_FALSE; - } + return LW_FALSE; + } - if ( l1->points->npoints < 2 || l2->points->npoints < 2 ) { - /* TODO: return distance between these two points */ + if ( l1->points->npoints < 2 || l2->points->npoints < 2 ) + { + /* TODO: return distance between these two points */ lwerror("Both input lines must have at least 2 points"); - return LW_FALSE; - } - - /* WARNING: these ranges may be wider than real ones */ - gbox1 = lwgeom_get_bbox(g1); - gbox2 = lwgeom_get_bbox(g2); - - assert(gbox1); /* or the npoints check above would have failed */ - assert(gbox2); /* or the npoints check above would have failed */ - - /* - * Find overlapping M range - * WARNING: may be larger than the real one - */ - - tmin = FP_MAX(gbox1->mmin, gbox2->mmin); - tmax = FP_MIN(gbox1->mmax, gbox2->mmax); - - if ( tmax < tmin ) { - LWDEBUG(1, "Inputs never exist at the same time"); - return LW_FALSE; - } - - // lwnotice("Min:%g, Max:%g", tmin, tmax); - - /* - * Collect M values in common time range from inputs - */ - - mvals = lwalloc( sizeof(double) * - ( l1->points->npoints + l2->points->npoints ) ); - - /* TODO: also clip the lines ? */ - nmvals = ptarray_collect_mvals(l1->points, tmin, tmax, mvals); - nmvals += ptarray_collect_mvals(l2->points, tmin, tmax, mvals + nmvals); - - /* Sort values in ascending order */ - qsort(mvals, nmvals, sizeof(double), compare_double); - - /* Remove duplicated values */ - nmvals = uniq(mvals, nmvals); - - if ( nmvals < 2 ) - {{ - /* there's a single time, must be that one... */ - double t0 = mvals[0]; - POINT4D p0, p1; - LWDEBUGF(1, "Inputs only exist both at a single time (%g)", t0); - if ( -1 == ptarray_locate_along_linear(l1->points, t0, &p0, 0) ) - { - lwerror("Could not find point with M=%g on first geom", t0); - return -1; - } - if ( -1 == ptarray_locate_along_linear(l2->points, t0, &p1, 0) ) - { - lwerror("Could not find point with M=%g on second geom", t0); - return -1; - } - if ( distance3d_pt_pt((POINT3D*)&p0, (POINT3D*)&p1) <= maxdist ) - within = LW_TRUE; - lwfree(mvals); - return within; - }} - - /* - * For each consecutive pair of measures, compute time of closest point - * approach and actual distance between points at that time - */ - for (i=1; ipoints, t0, &p0, 0); - if ( -1 == seg ) continue; /* possible, if GBOX is approximated */ - // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t0, seg, p0.x, p0.y, p0.z); - - seg = ptarray_locate_along_linear(l1->points, t1, &p1, seg); - if ( -1 == seg ) continue; /* possible, if GBOX is approximated */ - // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t1, seg, p1.x, p1.y, p1.z); - - seg = ptarray_locate_along_linear(l2->points, t0, &q0, 0); - if ( -1 == seg ) continue; /* possible, if GBOX is approximated */ - // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t0, seg, q0.x, q0.y, q0.z); - - seg = ptarray_locate_along_linear(l2->points, t1, &q1, seg); - if ( -1 == seg ) continue; /* possible, if GBOX is approximated */ - // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t1, seg, q1.x, q1.y, q1.z); - - segments_tcpa(&p0, &p1, &q0, &q1, t0, t1); + return LW_FALSE; + } + + /* WARNING: these ranges may be wider than real ones */ + gbox1 = lwgeom_get_bbox(g1); + gbox2 = lwgeom_get_bbox(g2); + + assert(gbox1); /* or the npoints check above would have failed */ + assert(gbox2); /* or the npoints check above would have failed */ + + /* + * Find overlapping M range + * WARNING: may be larger than the real one + */ + + tmin = FP_MAX(gbox1->mmin, gbox2->mmin); + tmax = FP_MIN(gbox1->mmax, gbox2->mmax); + + if ( tmax < tmin ) + { + LWDEBUG(1, "Inputs never exist at the same time"); + return LW_FALSE; + } + + // lwnotice("Min:%g, Max:%g", tmin, tmax); + + /* + * Collect M values in common time range from inputs + */ + + mvals = lwalloc( sizeof(double) * + ( l1->points->npoints + l2->points->npoints ) ); + + /* TODO: also clip the lines ? */ + nmvals = ptarray_collect_mvals(l1->points, tmin, tmax, mvals); + nmvals += ptarray_collect_mvals(l2->points, tmin, tmax, mvals + nmvals); + + /* Sort values in ascending order */ + qsort(mvals, nmvals, sizeof(double), compare_double); + + /* Remove duplicated values */ + nmvals = uniq(mvals, nmvals); + + if ( nmvals < 2 ) + { + { + /* there's a single time, must be that one... */ + double t0 = mvals[0]; + POINT4D p0, p1; + LWDEBUGF(1, "Inputs only exist both at a single time (%g)", t0); + if ( -1 == ptarray_locate_along_linear(l1->points, t0, &p0, 0) ) + { + lwerror("Could not find point with M=%g on first geom", t0); + return -1; + } + if ( -1 == ptarray_locate_along_linear(l2->points, t0, &p1, 0) ) + { + lwerror("Could not find point with M=%g on second geom", t0); + return -1; + } + if ( distance3d_pt_pt((POINT3D*)&p0, (POINT3D*)&p1) <= maxdist ) + within = LW_TRUE; + lwfree(mvals); + return within; + } + } + + /* + * For each consecutive pair of measures, compute time of closest point + * approach and actual distance between points at that time + */ + for (i=1; ipoints, t0, &p0, 0); + if ( -1 == seg ) continue; /* possible, if GBOX is approximated */ + // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t0, seg, p0.x, p0.y, p0.z); + + seg = ptarray_locate_along_linear(l1->points, t1, &p1, seg); + if ( -1 == seg ) continue; /* possible, if GBOX is approximated */ + // lwnotice("Measure %g on segment %d of line 1: %g, %g, %g", t1, seg, p1.x, p1.y, p1.z); + + seg = ptarray_locate_along_linear(l2->points, t0, &q0, 0); + if ( -1 == seg ) continue; /* possible, if GBOX is approximated */ + // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t0, seg, q0.x, q0.y, q0.z); + + seg = ptarray_locate_along_linear(l2->points, t1, &q1, seg); + if ( -1 == seg ) continue; /* possible, if GBOX is approximated */ + // lwnotice("Measure %g on segment %d of line 2: %g, %g, %g", t1, seg, q1.x, q1.y, q1.z); + + segments_tcpa(&p0, &p1, &q0, &q1, t0, t1); /* - lwnotice("Closest points: %g,%g,%g and %g,%g,%g at time %g", - p0.x, p0.y, p0.z, - q0.x, q0.y, q0.z, t); + lwnotice("Closest points: %g,%g,%g and %g,%g,%g at time %g", + p0.x, p0.y, p0.z, + q0.x, q0.y, q0.z, t); */ - dist2 = ( q0.x - p0.x ) * ( q0.x - p0.x ) + - ( q0.y - p0.y ) * ( q0.y - p0.y ) + - ( q0.z - p0.z ) * ( q0.z - p0.z ); - if ( dist2 <= maxdist2 ) - { - LWDEBUGF(1, "Within distance %g at time %g, breaking", *mindist, mintime); - within = LW_TRUE; - break; - } - } + dist2 = ( q0.x - p0.x ) * ( q0.x - p0.x ) + + ( q0.y - p0.y ) * ( q0.y - p0.y ) + + ( q0.z - p0.z ) * ( q0.z - p0.z ); + if ( dist2 <= maxdist2 ) + { + LWDEBUGF(1, "Within distance %g at time %g, breaking", *mindist, mintime); + within = LW_TRUE; + break; + } + } - /* - * Release memory - */ + /* + * Release memory + */ - lwfree(mvals); + lwfree(mvals); - return within; + return within; }