lwline_clip_to_ordinate_range()

static LWCOLLECTION * lwline_clip_to_ordinate_range ( const LWLINE *  line, char  ordinate, double  from, double  to  )

inlinestatic

Take in a LINESTRING and return a MULTILINESTRING of those portions of the LINESTRING between the from/to range for the specified ordinate (XYZM)

Definition at line 528 of file lwlinearreferencing.c.

{
        POINTARRAY *pa_in = NULL;
        LWCOLLECTION *lwgeom_out = NULL;
        POINTARRAY *dp = NULL;
        uint32_t i;
        int added_last_point = 0;
        POINT4D *p = NULL, *q = NULL, *r = NULL;
        double ordinate_value_p = 0.0, ordinate_value_q = 0.0;
        char hasz, hasm;
        char dims;
 
        /* Null input, nothing we can do. */
        assert(line);
        hasz = lwgeom_has_z(lwline_as_lwgeom(line));
        hasm = lwgeom_has_m(lwline_as_lwgeom(line));
        dims = FLAGS_NDIMS(line->flags);
 
        /* Asking for an ordinate we don't have. Error. */
        if ((ordinate == 'Z' && !hasz) || (ordinate == 'M' && !hasm))
        {
                lwerror("Cannot clip on ordinate %d in a %d-d geometry.", ordinate, dims);
                return NULL;
        }
 
        /* Prepare our working point objects. */
        p = lwalloc(sizeof(POINT4D));
        q = lwalloc(sizeof(POINT4D));
        r = lwalloc(sizeof(POINT4D));
 
        /* Construct a collection to hold our outputs. */
        lwgeom_out = lwcollection_construct_empty(MULTILINETYPE, line->srid, hasz, hasm);
 
        /* Get our input point array */
        pa_in = line->points;
 
        for (i = 0; i < pa_in->npoints; i++)
        {
                if (i > 0)
                {
                        *q = *p;
                        ordinate_value_q = ordinate_value_p;
                }
                getPoint4d_p(pa_in, i, p);
                ordinate_value_p = lwpoint_get_ordinate(p, ordinate);
 
                /* Is this point inside the ordinate range? Yes. */
                if (ordinate_value_p >= from && ordinate_value_p <= to)
                {
 
                        if (!added_last_point)
                        {
                                /* We didn't add the previous point, so this is a new segment.
                                 *  Make a new point array. */
                                dp = ptarray_construct_empty(hasz, hasm, 32);
 
                                /* We're transiting into the range so add an interpolated
                                 *  point at the range boundary.
                                 *  If we're on a boundary and crossing from the far side,
                                 *  we also need an interpolated point. */
                                if (i > 0 &&
                                    (/* Don't try to interpolate if this is the first point */
                                     (ordinate_value_p > from && ordinate_value_p < to) ||  /* Inside */
                                     (ordinate_value_p == from && ordinate_value_q > to) || /* Hopping from above */
                                     (ordinate_value_p == to && ordinate_value_q < from)))  /* Hopping from below */
                                {
                                        double interpolation_value;
                                        (ordinate_value_q > to) ? (interpolation_value = to)
                                                                : (interpolation_value = from);
                                        point_interpolate(q, p, r, hasz, hasm, ordinate, interpolation_value);
                                        ptarray_append_point(dp, r, LW_FALSE);
                                }
                        }
                        /* Add the current vertex to the point array. */
                        ptarray_append_point(dp, p, LW_FALSE);
                        if (ordinate_value_p == from || ordinate_value_p == to)
                                added_last_point = 2; /* Added on boundary. */
                        else
                                added_last_point = 1; /* Added inside range. */
                }
                /* Is this point inside the ordinate range? No. */
                else
                {
                        if (added_last_point == 1)
                        {
                                /* We're transiting out of the range, so add an interpolated point
                                 *  to the point array at the range boundary. */
                                double interpolation_value;
                                (ordinate_value_p > to) ? (interpolation_value = to) : (interpolation_value = from);
                                point_interpolate(q, p, r, hasz, hasm, ordinate, interpolation_value);
                                ptarray_append_point(dp, r, LW_FALSE);
                        }
                        else if (added_last_point == 2)
                        {
                                /* We're out and the last point was on the boundary.
                                 *  If the last point was the near boundary, nothing to do.
                                 *  If it was the far boundary, we need an interpolated point. */
                                if (from != to && ((ordinate_value_q == from && ordinate_value_p > from) ||
                                                   (ordinate_value_q == to && ordinate_value_p < to)))
                                {
                                        double interpolation_value;
                                        (ordinate_value_p > to) ? (interpolation_value = to)
                                                                : (interpolation_value = from);
                                        point_interpolate(q, p, r, hasz, hasm, ordinate, interpolation_value);
                                        ptarray_append_point(dp, r, LW_FALSE);
                                }
                        }
                        else if (i && ordinate_value_q < from && ordinate_value_p > to)
                        {
                                /* We just hopped over the whole range, from bottom to top,
                                 *  so we need to add *two* interpolated points! */
                                dp = ptarray_construct(hasz, hasm, 2);
                                /* Interpolate lower point. */
                                point_interpolate(p, q, r, hasz, hasm, ordinate, from);
                                ptarray_set_point4d(dp, 0, r);
                                /* Interpolate upper point. */
                                point_interpolate(p, q, r, hasz, hasm, ordinate, to);
                                ptarray_set_point4d(dp, 1, r);
                        }
                        else if (i && ordinate_value_q > to && ordinate_value_p < from)
                        {
                                /* We just hopped over the whole range, from top to bottom,
                                 *  so we need to add *two* interpolated points! */
                                dp = ptarray_construct(hasz, hasm, 2);
                                /* Interpolate upper point. */
                                point_interpolate(p, q, r, hasz, hasm, ordinate, to);
                                ptarray_set_point4d(dp, 0, r);
                                /* Interpolate lower point. */
                                point_interpolate(p, q, r, hasz, hasm, ordinate, from);
                                ptarray_set_point4d(dp, 1, r);
                        }
                        /* We have an extant point-array, save it out to a multi-line. */
                        if (dp)
                        {
                                /* Only one point, so we have to make an lwpoint to hold this
                                 *  and set the overall output type to a generic collection. */
                                if (dp->npoints == 1)
                                {
                                        LWPOINT *opoint = lwpoint_construct(line->srid, NULL, dp);
                                        lwgeom_out->type = COLLECTIONTYPE;
                                        lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(opoint));
                                }
                                else
                                {
                                        LWLINE *oline = lwline_construct(line->srid, NULL, dp);
                                        lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwline_as_lwgeom(oline));
                                }
 
                                /* Pointarray is now owned by lwgeom_out, so drop reference to it */
                                dp = NULL;
                        }
                        added_last_point = 0;
                }
        }
 
        /* Still some points left to be saved out. */
        if (dp)
        {
                if (dp->npoints == 1)
                {
                        LWPOINT *opoint = lwpoint_construct(line->srid, NULL, dp);
                        lwgeom_out->type = COLLECTIONTYPE;
                        lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(opoint));
                }
                else if (dp->npoints > 1)
                {
                        LWLINE *oline = lwline_construct(line->srid, NULL, dp);
                        lwgeom_out = lwcollection_add_lwgeom(lwgeom_out, lwline_as_lwgeom(oline));
                }
                else
                        ptarray_free(dp);
        }
 
        lwfree(p);
        lwfree(q);
        lwfree(r);
 
        if (line->bbox && lwgeom_out->ngeoms > 0)
                lwgeom_refresh_bbox((LWGEOM *)lwgeom_out);
 
        return lwgeom_out;
}

References LWLINE::bbox, COLLECTIONTYPE, LWLINE::flags, FLAGS_NDIMS, getPoint4d_p(), LW_FALSE, lwalloc(), lwcollection_add_lwgeom(), lwcollection_construct_empty(), lwerror(), lwfree(), lwgeom_has_m(), lwgeom_has_z(), lwgeom_refresh_bbox(), lwline_as_lwgeom(), lwline_construct(), lwpoint_as_lwgeom(), lwpoint_construct(), lwpoint_get_ordinate(), MULTILINETYPE, LWCOLLECTION::ngeoms, POINTARRAY::npoints, point_interpolate(), LWLINE::points, ptarray_append_point(), ptarray_construct(), ptarray_construct_empty(), ptarray_free(), ptarray_set_point4d(), r, LWLINE::srid, and LWCOLLECTION::type.

Referenced by lwgeom_clip_to_ordinate_range().

Here is the call graph for this function:

Here is the caller graph for this function: