#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;
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. */
pn->y += cos(theta) * offset;
}
- return LW_TRUE;
+ return LW_TRUE;
}
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++ )
{
/* Add our new point to the array */
ptarray_append_point(dpa, &pn, 0);
- }
+ }
return dpa;
}
/* Return degenerates upwards */
if ( ! lwline ) return NULL;
-
+
/* Create empty return shell */
srid = lwgeom_get_srid(lwg);
hasz = lwgeom_has_z(lwg);
/* 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);
}
/* 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);
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;
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));
if ( FP_EQUALS(m, point_m) )
{
lwmpoint_add_lwpoint(lwout, lwpoint_clone(lwin->geoms[i]));
- }
+ }
}
return lwout;
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:
}
-/**
+/**
* Given a point, ordinate number and value, set that ordinate on the
* point.
*/
}
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;
/* 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 )
{
/* 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);
{
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 )
{
lwgeom_drop_bbox((LWGEOM*)lwgeom_out);
lwgeom_add_bbox((LWGEOM*)lwgeom_out);
}
-
+
if ( ! homogeneous )
{
lwgeom_out->type = COLLECTIONTYPE;
LWPOINT *opoint = lwpoint_construct(line->srid, NULL, dp);
lwgeom_out->type = COLLECTIONTYPE;
lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(opoint));
-
+
}
else
{
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++ )
{
}
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));
}
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
*/
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; i<pa->npoints; ++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; i<pa->npoints; ++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;
}
static int
uniq(double *vals, int nvals)
{
- int i, last=0;
- for (i=1; i<nvals; ++i)
- {
- // lwnotice("(I%d):%g", i, vals[i]);
- if ( vals[i] != vals[last] )
- {
- vals[++last] = vals[i];
- // lwnotice("(O%d):%g", last, vals[last]);
- }
- }
- return last+1;
+ int i, last=0;
+ for (i=1; i<nvals; ++i)
+ {
+ // lwnotice("(I%d):%g", i, vals[i]);
+ if ( vals[i] != vals[last] )
+ {
+ vals[++last] = vals[i];
+ // lwnotice("(O%d):%g", last, vals[last]);
+ }
+ }
+ return last+1;
}
/*
{
int i = from;
POINT4D p1, p2;
-
+
/* Walk through each segment in the point array */
- getPoint4d_p(pa, i, &p1);
+ getPoint4d_p(pa, i, &p1);
for ( i = from+1; i < pa->npoints; i++ )
{
getPoint4d_p(pa, i, &p2);
if ( segment_locate_along(&p1, &p2, m, 0, p) == LW_TRUE )
return i-1; /* found */
- p1 = p2;
- }
+ p1 = p2;
+ }
return -1; /* not found */
}
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; i<nmvals; ++i)
- {
- double t0 = mvals[i-1];
- double t1 = mvals[i];
- double t;
- POINT4D p0, p1, q0, q1;
- int seg;
- double dist2;
-
- // lwnotice("T %g-%g", t0, t1);
-
- seg = ptarray_locate_along_linear(l1->points, 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; i<nmvals; ++i)
+ {
+ double t0 = mvals[i-1];
+ double t1 = mvals[i];
+ double t;
+ POINT4D p0, p1, q0, q1;
+ int seg;
+ double dist2;
+
+ // lwnotice("T %g-%g", t0, t1);
+
+ seg = ptarray_locate_along_linear(l1->points, 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; i<nmvals; ++i)
- {
- double t0 = mvals[i-1];
- double t1 = mvals[i];
- POINT4D p0, p1, q0, q1;
- int seg;
- double dist2;
-
- // lwnotice("T %g-%g", t0, t1);
-
- seg = ptarray_locate_along_linear(l1->points, 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; i<nmvals; ++i)
+ {
+ double t0 = mvals[i-1];
+ double t1 = mvals[i];
+ POINT4D p0, p1, q0, q1;
+ int seg;
+ double dist2;
+
+ // lwnotice("T %g-%g", t0, t1);
+
+ seg = ptarray_locate_along_linear(l1->points, 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;
}
projects / postgis / commitdiff