lwgeom_tcpa()

double lwgeom_tcpa	(	const LWGEOM *	g1,
		const LWGEOM *	g2,
		double *	mindist
	)

Find the time of closest point of approach.

Parameters

mindist

if not null will be set to the minimum distance between the trajectories at the closest point of approach.

Returns

the time value in which the minimum distance was reached, -1 if inputs are invalid (lwerror is called in that case), -2 if the trajectories do not share any point in time.

Definition at line 1071 of file lwlinearreferencing.c.

References compare_double(), distance3d_pt_pt(), FP_MAX, FP_MIN, lwalloc(), LWDEBUG, LWDEBUGF, lwerror(), lwfree(), lwgeom_as_lwline(), lwgeom_calculate_gbox(), lwgeom_has_m(), GBOX::mmax, GBOX::mmin, POINTARRAY::npoints, LWLINE::points, ptarray_collect_mvals(), ptarray_locate_along_linear(), segments_tcpa(), uniq(), POINT4D::x, POINT4D::y, and POINT4D::z.

Referenced by ST_ClosestPointOfApproach(), ST_DistanceCPA(), and test_lwgeom_tcpa().

{
        LWLINE *l1, *l2;
        int i;
        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;
        }

        l1 = lwgeom_as_lwline(g1);
        l2 = lwgeom_as_lwline(g2);

        if ( ! l1 || ! l2 )
        {
                lwerror("Both input geometries must be linestrings");
                return -1;
        }

        if ( l1->points->npoints < 2 || l2->points->npoints < 2 )
        {
                lwerror("Both input lines must have at least 2 points");
                return -1;
        }

        /* We use lwgeom_calculate_gbox() instead of lwgeom_get_gbox() */
        /* because we cannot afford the float rounding inaccuracy when */
        /* we compare the ranges for overlap below */
        lwgeom_calculate_gbox(g1, &gbox1);
        lwgeom_calculate_gbox(g2, &gbox2);

        /*
         * 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);
                */

                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;
}

Here is the call graph for this function:

Here is the caller graph for this function: