ptarray_locate_point_spheroid()

double ptarray_locate_point_spheroid	(	const POINTARRAY *	pa,
		const POINT4D *	p4d,
		const SPHEROID *	s,
		double	tolerance,
		double *	mindistout,
		POINT4D *	proj4d
	)

Locate a point along the point array defining a geographic line.

Definition at line 361 of file lwgeodetic_measures.c.

{
        GEOGRAPHIC_EDGE e;
        GEOGRAPHIC_POINT a, b, nearest = {0}; /* make compiler quiet */
        POINT4D p1, p2;
        const POINT2D *p;
        POINT2D proj;
        uint32_t i, seg = 0;
        int use_sphere = (s->a == s->b ? 1 : 0);
        int hasz;
        double za = 0.0, zb = 0.0;
        double distance,
                length,   /* Used for computing lengths */
                seglength = 0.0,  /* length of the segment where the closest point is located */
                partlength = 0.0, /* length from the beginning of the point array to the closest point */
                totlength = 0.0;  /* length of the point array */
 
        /* Initialize our point */
        geographic_point_init(p4d->x, p4d->y, &a);
 
        /* Handle point/point case here */
        if ( pa->npoints <= 1)
        {
                double mindist = 0.0;
                if ( pa->npoints == 1 )
                {
                        p = getPoint2d_cp(pa, 0);
                        geographic_point_init(p->x, p->y, &b);
                        /* Sphere special case, axes equal */
                        mindist = s->radius * sphere_distance(&a, &b);
                        /* If close or greater than tolerance, get the real answer to be sure */
                        if ( ! use_sphere || mindist > 0.95 * tolerance )
                        {
                                mindist = spheroid_distance(&a, &b, s);
                        }
                }
                if ( mindistout ) *mindistout = mindist;
                return 0.0;
        }
 
        /* Make result really big, so that everything will be smaller than it */
        distance = FLT_MAX;
 
        /* Initialize start of line */
        p = getPoint2d_cp(pa, 0);
        geographic_point_init(p->x, p->y, &(e.start));
 
        /* Iterate through the edges in our line */
        for ( i = 1; i < pa->npoints; i++ )
        {
                double d;
                p = getPoint2d_cp(pa, i);
                geographic_point_init(p->x, p->y, &(e.end));
                /* Get the spherical distance between point and edge */
                d = s->radius * edge_distance_to_point(&e, &a, &b);
                /* New shortest distance! Record this distance / location / segment */
                if ( d < distance )
                {
                        distance = d;
                        nearest = b;
                        seg = i - 1;
                }
            /* We've gotten closer than the tolerance... */
                if ( d < tolerance )
                {
                        /* Working on a sphere? The answer is correct, return */
                        if ( use_sphere )
                        {
                                break;
                        }
                        /* Far enough past the tolerance that the spheroid calculation won't change things */
                        else if ( d < tolerance * 0.95 )
                        {
                                break;
                        }
                        /* On a spheroid and near the tolerance? Confirm that we are *actually* closer than tolerance */
                        else
                        {
                                d = spheroid_distance(&a, &nearest, s);
                                /* Yes, closer than tolerance, return! */
                                if ( d < tolerance )
                                        break;
                        }
                }
                e.start = e.end;
        }
 
        if ( mindistout ) *mindistout = distance;
 
        /* See if we have a third dimension */
        hasz = (bool) FLAGS_GET_Z(pa->flags);
 
        /* Initialize first point of array */
        getPoint4d_p(pa, 0, &p1);
        geographic_point_init(p1.x, p1.y, &a);
        if ( hasz )
                za = p1.z;
 
        /* Loop and sum the length for each segment */
        for ( i = 1; i < pa->npoints; i++ )
        {
                getPoint4d_p(pa, i, &p1);
                geographic_point_init(p1.x, p1.y, &b);
                if ( hasz )
                        zb = p1.z;
 
                /* Special sphere case */
                if ( s->a == s->b )
                        length = s->radius * sphere_distance(&a, &b);
                /* Spheroid case */
                else
                        length = spheroid_distance(&a, &b, s);
 
                /* Add in the vertical displacement if we're in 3D */
                if ( hasz )
                        length = sqrt( (zb-za)*(zb-za) + length*length );
 
                /* Add this segment length to the total length */
                totlength += length;
 
                /* Add this segment length to the partial length */
                if (i - 1 < seg)
                        partlength += length;
                else if (i - 1 == seg)
                        /* Save segment length for computing the final value of partlength */
                        seglength = length;
 
                /* B gets incremented in the next loop, so we save the value here */
                a = b;
                za = zb;
        }
 
        /* Copy nearest into 2D/4D holder */
        proj4d->x = proj.x = rad2deg(nearest.lon);
        proj4d->y = proj.y = rad2deg(nearest.lat);
 
        /* Compute distance from beginning of the segment to closest point */
 
        /* Start of the segment */
        getPoint4d_p(pa, seg, &p1);
        geographic_point_init(p1.x, p1.y, &a);
 
        /* Closest point */
        geographic_point_init(proj4d->x, proj4d->y, &b);
 
        /* Special sphere case */
        if ( s->a == s->b )
                length = s->radius * sphere_distance(&a, &b);
        /* Spheroid case */
        else
                length = spheroid_distance(&a, &b, s);
 
        if ( hasz )
        {
                /* Compute Z and M values for closest point */
                double f = length / seglength;
                getPoint4d_p(pa, seg + 1, &p2);
                proj4d->z = p1.z + ((p2.z - p1.z) * f);
                proj4d->m = p1.m + ((p2.m - p1.m) * f);
                /* Add in the vertical displacement if we're in 3D */
                za = p1.z;
                zb = proj4d->z;
                length = sqrt( (zb-za)*(zb-za) + length*length );
        }
 
        /* Add this segment length to the total */
        partlength += length;
 
        /* Location of any point on a zero-length line is 0 */
        /* See http://trac.osgeo.org/postgis/ticket/1772#comment:2 */
        if ( partlength == 0 || totlength == 0 )
                return 0.0;
 
        /* For robustness, force 0 when closest point == startpoint */
        p = getPoint2d_cp(pa, 0);
        if ( seg == 0 && p2d_same(&proj, p) )
                return 0.0;
 
        /* For robustness, force 1 when closest point == endpoint */
        p = getPoint2d_cp(pa, pa->npoints - 1);
        if ( (seg >= (pa->npoints-2)) && p2d_same(&proj, p) )
                return 1.0;
 
        return partlength / totlength;
}

References distance(), edge_distance_to_point(), GEOGRAPHIC_EDGE::end, POINTARRAY::flags, FLAGS_GET_Z, geographic_point_init(), getPoint2d_cp(), getPoint4d_p(), GEOGRAPHIC_POINT::lat, GEOGRAPHIC_POINT::lon, POINT4D::m, POINTARRAY::npoints, p2d_same(), rad2deg, s, sphere_distance(), spheroid_distance(), GEOGRAPHIC_EDGE::start, POINT2D::x, POINT4D::x, POINT2D::y, POINT4D::y, and POINT4D::z.

Referenced by geography_line_locate_point().

Here is the call graph for this function:

Here is the caller graph for this function: