/src/postgis/liblwgeom/lwgeom_geos_split.c

Source
/**********************************************************************
 *
 * PostGIS - Spatial Types for PostgreSQL
 * http://postgis.net
 *
 * PostGIS is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 2 of the License, or
 * (at your option) any later version.
 *
 * PostGIS is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with PostGIS.  If not, see <http://www.gnu.org/licenses/>.
 *
 **********************************************************************
 *
 * Copyright 2011-2015 Sandro Santilli <strk@kbt.io>
 *
 **********************************************************************/

#include "../postgis_config.h"
/*#define POSTGIS_DEBUG_LEVEL 4*/
#include "lwgeom_geos.h"
#include "liblwgeom_internal.h"

#include <string.h>
#include <assert.h>

static LWGEOM* lwline_split_by_line(const LWLINE* lwgeom_in, const LWGEOM* blade_in);
static LWGEOM* lwline_split_by_point(const LWLINE* lwgeom_in, const LWPOINT* blade_in);
static LWGEOM* lwline_split_by_mpoint(const LWLINE* lwgeom_in, const LWMPOINT* blade_in);
static LWGEOM* lwline_split(const LWLINE* lwgeom_in, const LWGEOM* blade_in);
static LWGEOM* lwpoly_split_by_line(const LWPOLY* lwgeom_in, const LWGEOM* blade_in);
static LWGEOM* lwcollection_split(const LWCOLLECTION* lwcoll_in, const LWGEOM* blade_in);
static LWGEOM* lwpoly_split(const LWPOLY* lwpoly_in, const LWGEOM* blade_in);

/* Initializes and uses GEOS internally */
static LWGEOM*
lwline_split_by_line(const LWLINE* lwline_in, const LWGEOM* blade_in)
{
  LWGEOM** components;
  LWGEOM* diff;
  LWCOLLECTION* out;
  GEOSGeometry* gdiff; /* difference */
  GEOSGeometry* g1;
  GEOSGeometry* g2;
  int ret;

  /* ASSERT blade_in is LINE or MULTILINE */
  assert (blade_in->type == LINETYPE ||
          blade_in->type == MULTILINETYPE ||
          blade_in->type == POLYGONTYPE ||
          blade_in->type == MULTIPOLYGONTYPE );

  /* Possible outcomes:
   *
   *  1. The lines do not cross or overlap
   *      -> Return a collection with single element
   *  2. The lines cross
   *      -> Return a collection of all elements resulting from the split
   */

  initGEOS(lwgeom_geos_error, lwgeom_geos_error);

  g1 = LWGEOM2GEOS((LWGEOM*)lwline_in, 0);
  if ( ! g1 )
  {
    lwerror("LWGEOM2GEOS: %s", lwgeom_geos_errmsg);
    return NULL;
  }
  g2 = LWGEOM2GEOS(blade_in, 0);
  if ( ! g2 )
  {
    GEOSGeom_destroy(g1);
    lwerror("LWGEOM2GEOS: %s", lwgeom_geos_errmsg);
    return NULL;
  }

  /* If blade is a polygon, pick its boundary */
  if ( blade_in->type == POLYGONTYPE || blade_in->type == MULTIPOLYGONTYPE )
  {
    gdiff = GEOSBoundary(g2);
    GEOSGeom_destroy(g2);
    if ( ! gdiff )
    {
      GEOSGeom_destroy(g1);
      lwerror("GEOSBoundary: %s", lwgeom_geos_errmsg);
      return NULL;
    }
    g2 = gdiff; gdiff = NULL;
  }

  /* If interior intersection is linear we can't split */
  ret = GEOSRelatePattern(g1, g2, "1********");
  if ( 2 == ret )
  {
    lwerror("GEOSRelatePattern: %s", lwgeom_geos_errmsg);
    GEOSGeom_destroy(g1);
    GEOSGeom_destroy(g2);
    return NULL;
  }
  if ( ret )
  {
    GEOSGeom_destroy(g1);
    GEOSGeom_destroy(g2);
    lwerror("Splitter line has linear intersection with input");
    return NULL;
  }


  gdiff = GEOSDifference(g1,g2);
  GEOSGeom_destroy(g1);
  GEOSGeom_destroy(g2);
  if (gdiff == NULL)
  {
    lwerror("GEOSDifference: %s", lwgeom_geos_errmsg);
    return NULL;
  }

  diff = GEOS2LWGEOM(gdiff, FLAGS_GET_Z(lwline_in->flags));
  GEOSGeom_destroy(gdiff);
  if (NULL == diff)
  {
    lwerror("GEOS2LWGEOM: %s", lwgeom_geos_errmsg);
    return NULL;
  }

  out = lwgeom_as_lwcollection(diff);
  if ( ! out )
  {
    components = lwalloc(sizeof(LWGEOM*)*1);
    components[0] = diff;
    out = lwcollection_construct(COLLECTIONTYPE, lwline_in->srid,
                                 NULL, 1, components);
  }
  else
  {
    /* Set SRID */
    lwgeom_set_srid((LWGEOM*)out, lwline_in->srid);
    /* Force collection type */
    out->type = COLLECTIONTYPE;
  }


  return (LWGEOM*)out;
}

static LWGEOM*
lwline_split_by_point(const LWLINE* lwline_in, const LWPOINT* blade_in)
{
  LWMLINE* out;

  out = lwmline_construct_empty(lwline_in->srid,
    FLAGS_GET_Z(lwline_in->flags),
    FLAGS_GET_M(lwline_in->flags));
  if ( lwline_split_by_point_to(lwline_in, blade_in, out) < 2 )
  {
    lwmline_add_lwline(out, lwline_clone_deep(lwline_in));
  }

  /* Turn multiline into collection */
  out->type = COLLECTIONTYPE;

  return (LWGEOM*)out;
}

static LWGEOM*
lwline_split_by_mpoint(const LWLINE* lwline_in, const LWMPOINT* mp)
{
  LWMLINE* out;
  uint32_t i, j;

  out = lwmline_construct_empty(lwline_in->srid,
          FLAGS_GET_Z(lwline_in->flags),
          FLAGS_GET_M(lwline_in->flags));
  lwmline_add_lwline(out, lwline_clone_deep(lwline_in));

  for (i=0; i<mp->ngeoms; ++i)
  {
    for (j=0; j<out->ngeoms; ++j)
    {
      lwline_in = out->geoms[j];
      LWPOINT *blade_in = mp->geoms[i];
      int ret = lwline_split_by_point_to(lwline_in, blade_in, out);
      if ( 2 == ret )
      {
        /* the point splits this line,
         * 2 splits were added to collection.
         * We'll move the latest added into
         * the slot of the current one.
         */
        lwline_free(out->geoms[j]);
        out->geoms[j] = out->geoms[--out->ngeoms];
      }
    }
  }

  /* Turn multiline into collection */
  out->type = COLLECTIONTYPE;

  return (LWGEOM*)out;
}

int
lwline_split_by_point_to(const LWLINE* lwline_in, const LWPOINT* blade_in,
                         LWMLINE* v)
{
  double mindist_sqr = -1;
  POINT4D pt, pt_projected;
  POINT4D p1, p2;
  POINTARRAY *ipa = lwline_in->points;
  POINTARRAY* pa1;
  POINTARRAY* pa2;
  uint32_t i, nsegs, seg = UINT32_MAX;

  /* Possible outcomes:
   *
   *  1. The point is not on the line or on the boundary
   *      -> Leave collection untouched, return 0
   *  2. The point is on the boundary
   *      -> Leave collection untouched, return 1
   *  3. The point is in the line
   *      -> Push 2 elements on the collection:
   *         o start_point - cut_point
   *         o cut_point - last_point
   *      -> Return 2
   */

  getPoint4d_p(blade_in->point, 0, &pt);

  /* Find closest segment */
  if ( ipa->npoints < 1 ) return 0; /* empty input line */
  getPoint4d_p(ipa, 0, &p1);
  nsegs = ipa->npoints - 1;
  for ( i = 0; i < nsegs; i++ )
  {
    getPoint4d_p(ipa, i+1, &p2);
    double dist_sqr = distance2d_sqr_pt_seg((POINT2D *)&pt, (POINT2D *)&p1, (POINT2D *)&p2);
    LWDEBUGF(4, "Distance (squared) of point %.15g %.15g to segment %.15g %.15g, %.15g %.15g: %.15g",
         pt.x, pt.y,
         p1.x, p1.y,
         p2.x, p2.y,
         dist_sqr);
    if (i == 0 || dist_sqr < mindist_sqr)
    {
      mindist_sqr = dist_sqr;
      seg=i;
      if (mindist_sqr == 0.0)
        break; /* can't be closer than ON line */
    }
    p1 = p2;
  }

  LWDEBUGF(3, "Closest segment: %d", seg);
  LWDEBUGF(3, "mindist: %.15g", mindist_sqr);

  /* No intersection */
  if (mindist_sqr > 0)
    return 0;

  /* empty or single-point line, intersection on boundary */
  if ( seg == UINT32_MAX ) return 1;

  /*
   * We need to project the
   * point on the closest segment,
   * to interpolate Z and M if needed
   */
  getPoint4d_p(ipa, seg, &p1);
  getPoint4d_p(ipa, seg+1, &p2);
  closest_point_on_segment(&pt, &p1, &p2, &pt_projected);
  /* But X and Y we want the ones of the input point,
   * as on some architectures the interpolation math moves the
   * coordinates (see #3422)
   */
  pt_projected.x = pt.x;
  pt_projected.y = pt.y;

  LWDEBUGF(3, "Projected point:(%.15g %.15g), seg:%d, p1:(%.15g %.15g), p2:(%.15g %.15g)", pt_projected.x, pt_projected.y, seg, p1.x, p1.y, p2.x, p2.y);

  /* When closest point == an endpoint, this is a boundary intersection */
  if ( ( (seg == nsegs-1) && P4D_SAME_STRICT(&pt_projected, &p2) ) ||
       ( (seg == 0)       && P4D_SAME_STRICT(&pt_projected, &p1) ) )
  {
    return 1;
  }

  /* This is an internal intersection, let's build the two new pointarrays */

  pa1 = ptarray_construct_empty(FLAGS_GET_Z(ipa->flags), FLAGS_GET_M(ipa->flags), seg+2);
  /* TODO: replace with a memcpy ? */
  for (i=0; i<=seg; ++i)
  {
    getPoint4d_p(ipa, i, &p1);
    ptarray_append_point(pa1, &p1, LW_FALSE);
  }
  ptarray_append_point(pa1, &pt_projected, LW_FALSE);

  pa2 = ptarray_construct_empty(FLAGS_GET_Z(ipa->flags), FLAGS_GET_M(ipa->flags), ipa->npoints-seg);
  ptarray_append_point(pa2, &pt_projected, LW_FALSE);
  /* TODO: replace with a memcpy (if so need to check for duplicated point) ? */
  for (i=seg+1; i<ipa->npoints; ++i)
  {
    getPoint4d_p(ipa, i, &p1);
    ptarray_append_point(pa2, &p1, LW_FALSE);
  }

  /* NOTE: I've seen empty pointarrays with loc != 0 and loc != 1 */
  if ( pa1->npoints == 0 || pa2->npoints == 0 ) {
    ptarray_free(pa1);
    ptarray_free(pa2);
    /* Intersection is on the boundary */
    return 1;
  }

  lwmline_add_lwline(v, lwline_construct(SRID_UNKNOWN, NULL, pa1));
  lwmline_add_lwline(v, lwline_construct(SRID_UNKNOWN, NULL, pa2));
  return 2;
}

static LWGEOM*
lwline_split(const LWLINE* lwline_in, const LWGEOM* blade_in)
{
  switch (blade_in->type)
  {
  case POINTTYPE:
    return lwline_split_by_point(lwline_in, (LWPOINT*)blade_in);
  case MULTIPOINTTYPE:
    return lwline_split_by_mpoint(lwline_in, (LWMPOINT*)blade_in);

  case LINETYPE:
  case MULTILINETYPE:
  case POLYGONTYPE:
  case MULTIPOLYGONTYPE:
    return lwline_split_by_line(lwline_in, blade_in);

  default:
    lwerror("Splitting a Line by a %s is unsupported",
            lwtype_name(blade_in->type));
    return NULL;
  }
  return NULL;
}

/* Initializes and uses GEOS internally */
static LWGEOM*
lwpoly_split_by_line(const LWPOLY* lwpoly_in, const LWGEOM* blade_in)
{
  LWCOLLECTION* out;
  GEOSGeometry* g1;
  GEOSGeometry* g2;
  GEOSGeometry* g1_bounds;
  GEOSGeometry* polygons;
  const GEOSGeometry *vgeoms[1];
  int i,n;
  int hasZ = FLAGS_GET_Z(lwpoly_in->flags);


  /* Possible outcomes:
   *
   *  1. The line does not split the polygon
   *      -> Return a collection with single element
   *  2. The line does split the polygon
   *      -> Return a collection of all elements resulting from the split
   */

  initGEOS(lwgeom_geos_error, lwgeom_geos_error);

  g1 = LWGEOM2GEOS((LWGEOM*)lwpoly_in, 0);
  if ( NULL == g1 )
  {
    lwerror("LWGEOM2GEOS: %s", lwgeom_geos_errmsg);
    return NULL;
  }
  g1_bounds = GEOSBoundary(g1);
  if ( NULL == g1_bounds )
  {
    GEOSGeom_destroy(g1);
    lwerror("GEOSBoundary: %s", lwgeom_geos_errmsg);
    return NULL;
  }

  g2 = LWGEOM2GEOS(blade_in, 0);
  if ( NULL == g2 )
  {
    GEOSGeom_destroy(g1);
    GEOSGeom_destroy(g1_bounds);
    lwerror("LWGEOM2GEOS: %s", lwgeom_geos_errmsg);
    return NULL;
  }

  vgeoms[0] = GEOSUnion(g1_bounds, g2);
  if ( NULL == vgeoms[0] )
  {
    GEOSGeom_destroy(g1);
    GEOSGeom_destroy(g2);
    GEOSGeom_destroy(g1_bounds);
    lwerror("GEOSUnion: %s", lwgeom_geos_errmsg);
    return NULL;
  }

  polygons = GEOSPolygonize(vgeoms, 1);
  if ( NULL == polygons )
  {
    GEOSGeom_destroy(g1);
    GEOSGeom_destroy(g2);
    GEOSGeom_destroy(g1_bounds);
    GEOSGeom_destroy((GEOSGeometry*)vgeoms[0]);
    lwerror("GEOSPolygonize: %s", lwgeom_geos_errmsg);
    return NULL;
  }

#ifndef NDEBUG
  if ( GEOSGeomTypeId(polygons) != COLLECTIONTYPE )
  {
    GEOSGeom_destroy(g1);
    GEOSGeom_destroy(g2);
    GEOSGeom_destroy(g1_bounds);
    GEOSGeom_destroy((GEOSGeometry*)vgeoms[0]);
    GEOSGeom_destroy(polygons);
    lwerror("%s [%d] Unexpected return from GEOSpolygonize", __FILE__, __LINE__);
    return 0;
  }
#endif

  /* We should now have all polygons, just skip
   * the ones which are in holes of the original
   * geometries and return the rest in a collection
   */
  n = GEOSGetNumGeometries(polygons);
  out = lwcollection_construct_empty(COLLECTIONTYPE, lwpoly_in->srid,
             hasZ, 0);
  /* Allocate space for all polys */
  out->geoms = lwrealloc(out->geoms, sizeof(LWGEOM*)*n);
  assert(0 == out->ngeoms);
  for (i=0; i<n; ++i)
  {
    GEOSGeometry* pos; /* point on surface */
    const GEOSGeometry* p = GEOSGetGeometryN(polygons, i);
    int contains;

    pos = GEOSPointOnSurface(p);
    if ( ! pos )
    {
      GEOSGeom_destroy(g1);
      GEOSGeom_destroy(g2);
      GEOSGeom_destroy(g1_bounds);
      GEOSGeom_destroy((GEOSGeometry*)vgeoms[0]);
      GEOSGeom_destroy(polygons);
      lwerror("GEOSPointOnSurface: %s", lwgeom_geos_errmsg);
      return NULL;
    }

    contains = GEOSContains(g1, pos);
    if ( 2 == contains )
    {
      GEOSGeom_destroy(g1);
      GEOSGeom_destroy(g2);
      GEOSGeom_destroy(g1_bounds);
      GEOSGeom_destroy((GEOSGeometry*)vgeoms[0]);
      GEOSGeom_destroy(polygons);
      GEOSGeom_destroy(pos);
      lwerror("GEOSContains: %s", lwgeom_geos_errmsg);
      return NULL;
    }

    GEOSGeom_destroy(pos);

    if ( 0 == contains )
    {
      /* Original geometry doesn't contain
       * a point in this ring, must be an hole
       */
      continue;
    }

    out->geoms[out->ngeoms++] = GEOS2LWGEOM(p, hasZ);
  }

  GEOSGeom_destroy(g1);
  GEOSGeom_destroy(g2);
  GEOSGeom_destroy(g1_bounds);
  GEOSGeom_destroy((GEOSGeometry*)vgeoms[0]);
  GEOSGeom_destroy(polygons);

  return (LWGEOM*)out;
}

static LWGEOM*
lwcollection_split(const LWCOLLECTION* lwcoll_in, const LWGEOM* blade_in)
{
  LWGEOM** split_vector=NULL;
  LWCOLLECTION* out;
  size_t split_vector_capacity;
  size_t split_vector_size=0;
  size_t i,j;

  split_vector_capacity=8;
  split_vector = lwalloc(split_vector_capacity * sizeof(LWGEOM*));
  if ( ! split_vector )
  {
    lwerror("Out of virtual memory");
    return NULL;
  }

  for (i=0; i<lwcoll_in->ngeoms; ++i)
  {
    LWCOLLECTION* col;
    LWGEOM* split = lwgeom_split(lwcoll_in->geoms[i], blade_in);
    /* an exception should prevent this from ever returning NULL */
    if ( ! split ) return NULL;

    col = lwgeom_as_lwcollection(split);
    /* Output, if any, will always be a collection */
    assert(col);

    /* Reallocate split_vector if needed */
    if ( split_vector_size + col->ngeoms > split_vector_capacity )
    {
      /* NOTE: we could be smarter on reallocations here */
      split_vector_capacity += col->ngeoms;
      split_vector = lwrealloc(split_vector,
                               split_vector_capacity * sizeof(LWGEOM*));
      if ( ! split_vector )
      {
        lwerror("Out of virtual memory");
        return NULL;
      }
    }

    for (j=0; j<col->ngeoms; ++j)
    {
      col->geoms[j]->srid = SRID_UNKNOWN; /* strip srid */
      split_vector[split_vector_size++] = col->geoms[j];
    }
    lwfree(col->geoms);
    lwfree(col);
  }

  /* Now split_vector has split_vector_size geometries */
  out = lwcollection_construct(COLLECTIONTYPE, lwcoll_in->srid,
                               NULL, split_vector_size, split_vector);

  return (LWGEOM*)out;
}

static LWGEOM*
lwpoly_split(const LWPOLY* lwpoly_in, const LWGEOM* blade_in)
{
  switch (blade_in->type)
  {
  case MULTILINETYPE:
  case LINETYPE:
    return lwpoly_split_by_line(lwpoly_in, blade_in);
  default:
    lwerror("Splitting a Polygon by a %s is unsupported",
            lwtype_name(blade_in->type));
    return NULL;
  }
  return NULL;
}

/* exported */
LWGEOM*
lwgeom_split(const LWGEOM* lwgeom_in, const LWGEOM* blade_in)
{
  switch (lwgeom_in->type)
  {
  case LINETYPE:
    return lwline_split((const LWLINE*)lwgeom_in, blade_in);

  case POLYGONTYPE:
    return lwpoly_split((const LWPOLY*)lwgeom_in, blade_in);

  case MULTIPOLYGONTYPE:
  case MULTILINETYPE:
  case COLLECTIONTYPE:
    return lwcollection_split((const LWCOLLECTION*)lwgeom_in, blade_in);

  default:
    lwerror("Splitting of %s geometries is unsupported",
            lwtype_name(lwgeom_in->type));
    return NULL;
  }

}


Coverage Report

Created: 2026-03-31 06:40

Line	Count	Source
1		/**********************************************************************
2		*
3		* PostGIS - Spatial Types for PostgreSQL
4		* http://postgis.net
5		*
6		* PostGIS is free software: you can redistribute it and/or modify
7		* it under the terms of the GNU General Public License as published by
8		* the Free Software Foundation, either version 2 of the License, or
9		* (at your option) any later version.
10		*
11		* PostGIS is distributed in the hope that it will be useful,
12		* but WITHOUT ANY WARRANTY; without even the implied warranty of
13		* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14		* GNU General Public License for more details.
15		*
16		* You should have received a copy of the GNU General Public License
17		* along with PostGIS. If not, see <http://www.gnu.org/licenses/>.
18		*
19		**********************************************************************
20		*
21		* Copyright 2011-2015 Sandro Santilli <strk@kbt.io>
22		*
23		**********************************************************************/
24
25		#include "../postgis_config.h"
26		/#define POSTGIS_DEBUG_LEVEL 4/
27		#include "lwgeom_geos.h"
28		#include "liblwgeom_internal.h"
29
30		#include <string.h>
31		#include <assert.h>
32
33		static LWGEOM* lwline_split_by_line(const LWLINE* lwgeom_in, const LWGEOM* blade_in);
34		static LWGEOM* lwline_split_by_point(const LWLINE* lwgeom_in, const LWPOINT* blade_in);
35		static LWGEOM* lwline_split_by_mpoint(const LWLINE* lwgeom_in, const LWMPOINT* blade_in);
36		static LWGEOM* lwline_split(const LWLINE* lwgeom_in, const LWGEOM* blade_in);
37		static LWGEOM* lwpoly_split_by_line(const LWPOLY* lwgeom_in, const LWGEOM* blade_in);
38		static LWGEOM* lwcollection_split(const LWCOLLECTION* lwcoll_in, const LWGEOM* blade_in);
39		static LWGEOM* lwpoly_split(const LWPOLY* lwpoly_in, const LWGEOM* blade_in);
40
41		/* Initializes and uses GEOS internally */
42		static LWGEOM*
43		lwline_split_by_line(const LWLINE* lwline_in, const LWGEOM* blade_in)
44	0	{
45	0	LWGEOM** components;
46	0	LWGEOM* diff;
47	0	LWCOLLECTION* out;
48	0	GEOSGeometry* gdiff; /* difference */
49	0	GEOSGeometry* g1;
50	0	GEOSGeometry* g2;
51	0	int ret;
52
53		/* ASSERT blade_in is LINE or MULTILINE */
54	0	assert (blade_in->type == LINETYPE \|\|
55	0	blade_in->type == MULTILINETYPE \|\|
56	0	blade_in->type == POLYGONTYPE \|\|
57	0	blade_in->type == MULTIPOLYGONTYPE );
58
59		/* Possible outcomes:
60		*
61		* 1. The lines do not cross or overlap
62		* -> Return a collection with single element
63		* 2. The lines cross
64		* -> Return a collection of all elements resulting from the split
65		*/
66
67	0	initGEOS(lwgeom_geos_error, lwgeom_geos_error);
68
69	0	g1 = LWGEOM2GEOS((LWGEOM*)lwline_in, 0);
70	0	if ( ! g1 )
71	0	{
72	0	lwerror("LWGEOM2GEOS: %s", lwgeom_geos_errmsg);
73	0	return NULL;
74	0	}
75	0	g2 = LWGEOM2GEOS(blade_in, 0);
76	0	if ( ! g2 )
77	0	{
78	0	GEOSGeom_destroy(g1);
79	0	lwerror("LWGEOM2GEOS: %s", lwgeom_geos_errmsg);
80	0	return NULL;
81	0	}
82
83		/* If blade is a polygon, pick its boundary */
84	0	if ( blade_in->type == POLYGONTYPE \|\| blade_in->type == MULTIPOLYGONTYPE )
85	0	{
86	0	gdiff = GEOSBoundary(g2);
87	0	GEOSGeom_destroy(g2);
88	0	if ( ! gdiff )
89	0	{
90	0	GEOSGeom_destroy(g1);
91	0	lwerror("GEOSBoundary: %s", lwgeom_geos_errmsg);
92	0	return NULL;
93	0	}
94	0	g2 = gdiff; gdiff = NULL;
95	0	}
96
97		/* If interior intersection is linear we can't split */
98	0	ret = GEOSRelatePattern(g1, g2, "1********");
99	0	if ( 2 == ret )
100	0	{
101	0	lwerror("GEOSRelatePattern: %s", lwgeom_geos_errmsg);
102	0	GEOSGeom_destroy(g1);
103	0	GEOSGeom_destroy(g2);
104	0	return NULL;
105	0	}
106	0	if ( ret )
107	0	{
108	0	GEOSGeom_destroy(g1);
109	0	GEOSGeom_destroy(g2);
110	0	lwerror("Splitter line has linear intersection with input");
111	0	return NULL;
112	0	}
113
114
115	0	gdiff = GEOSDifference(g1,g2);
116	0	GEOSGeom_destroy(g1);
117	0	GEOSGeom_destroy(g2);
118	0	if (gdiff == NULL)
119	0	{
120	0	lwerror("GEOSDifference: %s", lwgeom_geos_errmsg);
121	0	return NULL;
122	0	}
123
124	0	diff = GEOS2LWGEOM(gdiff, FLAGS_GET_Z(lwline_in->flags));
125	0	GEOSGeom_destroy(gdiff);
126	0	if (NULL == diff)
127	0	{
128	0	lwerror("GEOS2LWGEOM: %s", lwgeom_geos_errmsg);
129	0	return NULL;
130	0	}
131
132	0	out = lwgeom_as_lwcollection(diff);
133	0	if ( ! out )
134	0	{
135	0	components = lwalloc(sizeof(LWGEOM)1);
136	0	components[0] = diff;
137	0	out = lwcollection_construct(COLLECTIONTYPE, lwline_in->srid,
138	0	NULL, 1, components);
139	0	}
140	0	else
141	0	{
142		/* Set SRID */
143	0	lwgeom_set_srid((LWGEOM*)out, lwline_in->srid);
144		/* Force collection type */
145	0	out->type = COLLECTIONTYPE;
146	0	}
147
148
149	0	return (LWGEOM*)out;
150	0	}
151
152		static LWGEOM*
153		lwline_split_by_point(const LWLINE* lwline_in, const LWPOINT* blade_in)
154	0	{
155	0	LWMLINE* out;
156
157	0	out = lwmline_construct_empty(lwline_in->srid,
158	0	FLAGS_GET_Z(lwline_in->flags),
159	0	FLAGS_GET_M(lwline_in->flags));
160	0	if ( lwline_split_by_point_to(lwline_in, blade_in, out) < 2 )
161	0	{
162	0	lwmline_add_lwline(out, lwline_clone_deep(lwline_in));
163	0	}
164
165		/* Turn multiline into collection */
166	0	out->type = COLLECTIONTYPE;
167
168	0	return (LWGEOM*)out;
169	0	}
170
171		static LWGEOM*
172		lwline_split_by_mpoint(const LWLINE* lwline_in, const LWMPOINT* mp)
173	0	{
174	0	LWMLINE* out;
175	0	uint32_t i, j;
176
177	0	out = lwmline_construct_empty(lwline_in->srid,
178	0	FLAGS_GET_Z(lwline_in->flags),
179	0	FLAGS_GET_M(lwline_in->flags));
180	0	lwmline_add_lwline(out, lwline_clone_deep(lwline_in));
181
182	0	for (i=0; i<mp->ngeoms; ++i)
183	0	{
184	0	for (j=0; j<out->ngeoms; ++j)
185	0	{
186	0	lwline_in = out->geoms[j];
187	0	LWPOINT *blade_in = mp->geoms[i];
188	0	int ret = lwline_split_by_point_to(lwline_in, blade_in, out);
189	0	if ( 2 == ret )
190	0	{
191		/* the point splits this line,
192		* 2 splits were added to collection.
193		* We'll move the latest added into
194		* the slot of the current one.
195		*/
196	0	lwline_free(out->geoms[j]);
197	0	out->geoms[j] = out->geoms[--out->ngeoms];
198	0	}
199	0	}
200	0	}
201
202		/* Turn multiline into collection */
203	0	out->type = COLLECTIONTYPE;
204
205	0	return (LWGEOM*)out;
206	0	}
207
208		int
209		lwline_split_by_point_to(const LWLINE* lwline_in, const LWPOINT* blade_in,
210		LWMLINE* v)
211	0	{
212	0	double mindist_sqr = -1;
213	0	POINT4D pt, pt_projected;
214	0	POINT4D p1, p2;
215	0	POINTARRAY *ipa = lwline_in->points;
216	0	POINTARRAY* pa1;
217	0	POINTARRAY* pa2;
218	0	uint32_t i, nsegs, seg = UINT32_MAX;
219
220		/* Possible outcomes:
221		*
222		* 1. The point is not on the line or on the boundary
223		* -> Leave collection untouched, return 0
224		* 2. The point is on the boundary
225		* -> Leave collection untouched, return 1
226		* 3. The point is in the line
227		* -> Push 2 elements on the collection:
228		* o start_point - cut_point
229		* o cut_point - last_point
230		* -> Return 2
231		*/
232
233	0	getPoint4d_p(blade_in->point, 0, &pt);
234
235		/* Find closest segment */
236	0	if ( ipa->npoints < 1 ) return 0; /* empty input line */
237	0	getPoint4d_p(ipa, 0, &p1);
238	0	nsegs = ipa->npoints - 1;
239	0	for ( i = 0; i < nsegs; i++ )
240	0	{
241	0	getPoint4d_p(ipa, i+1, &p2);
242	0	double dist_sqr = distance2d_sqr_pt_seg((POINT2D )&pt, (POINT2D )&p1, (POINT2D *)&p2);
243	0	LWDEBUGF(4, "Distance (squared) of point %.15g %.15g to segment %.15g %.15g, %.15g %.15g: %.15g",
244	0	pt.x, pt.y,
245	0	p1.x, p1.y,
246	0	p2.x, p2.y,
247	0	dist_sqr);
248	0	if (i == 0 \|\| dist_sqr < mindist_sqr)
249	0	{
250	0	mindist_sqr = dist_sqr;
251	0	seg=i;
252	0	if (mindist_sqr == 0.0)
253	0	break; /* can't be closer than ON line */
254	0	}
255	0	p1 = p2;
256	0	}
257
258	0	LWDEBUGF(3, "Closest segment: %d", seg);
259	0	LWDEBUGF(3, "mindist: %.15g", mindist_sqr);
260
261		/* No intersection */
262	0	if (mindist_sqr > 0)
263	0	return 0;
264
265		/* empty or single-point line, intersection on boundary */
266	0	if ( seg == UINT32_MAX ) return 1;
267
268		/*
269		* We need to project the
270		* point on the closest segment,
271		* to interpolate Z and M if needed
272		*/
273	0	getPoint4d_p(ipa, seg, &p1);
274	0	getPoint4d_p(ipa, seg+1, &p2);
275	0	closest_point_on_segment(&pt, &p1, &p2, &pt_projected);
276		/* But X and Y we want the ones of the input point,
277		* as on some architectures the interpolation math moves the
278		* coordinates (see #3422)
279		*/
280	0	pt_projected.x = pt.x;
281	0	pt_projected.y = pt.y;
282
283	0	LWDEBUGF(3, "Projected point:(%.15g %.15g), seg:%d, p1:(%.15g %.15g), p2:(%.15g %.15g)", pt_projected.x, pt_projected.y, seg, p1.x, p1.y, p2.x, p2.y);
284
285		/* When closest point == an endpoint, this is a boundary intersection */
286	0	if ( ( (seg == nsegs-1) && P4D_SAME_STRICT(&pt_projected, &p2) ) \|\|
287	0	( (seg == 0) && P4D_SAME_STRICT(&pt_projected, &p1) ) )
288	0	{
289	0	return 1;
290	0	}
291
292		/* This is an internal intersection, let's build the two new pointarrays */
293
294	0	pa1 = ptarray_construct_empty(FLAGS_GET_Z(ipa->flags), FLAGS_GET_M(ipa->flags), seg+2);
295		/* TODO: replace with a memcpy ? */
296	0	for (i=0; i<=seg; ++i)
297	0	{
298	0	getPoint4d_p(ipa, i, &p1);
299	0	ptarray_append_point(pa1, &p1, LW_FALSE);
300	0	}
301	0	ptarray_append_point(pa1, &pt_projected, LW_FALSE);
302
303	0	pa2 = ptarray_construct_empty(FLAGS_GET_Z(ipa->flags), FLAGS_GET_M(ipa->flags), ipa->npoints-seg);
304	0	ptarray_append_point(pa2, &pt_projected, LW_FALSE);
305		/* TODO: replace with a memcpy (if so need to check for duplicated point) ? */
306	0	for (i=seg+1; i<ipa->npoints; ++i)
307	0	{
308	0	getPoint4d_p(ipa, i, &p1);
309	0	ptarray_append_point(pa2, &p1, LW_FALSE);
310	0	}
311
312		/* NOTE: I've seen empty pointarrays with loc != 0 and loc != 1 */
313	0	if ( pa1->npoints == 0 \|\| pa2->npoints == 0 ) {
314	0	ptarray_free(pa1);
315	0	ptarray_free(pa2);
316		/* Intersection is on the boundary */
317	0	return 1;
318	0	}
319
320	0	lwmline_add_lwline(v, lwline_construct(SRID_UNKNOWN, NULL, pa1));
321	0	lwmline_add_lwline(v, lwline_construct(SRID_UNKNOWN, NULL, pa2));
322	0	return 2;
323	0	}
324
325		static LWGEOM*
326		lwline_split(const LWLINE* lwline_in, const LWGEOM* blade_in)
327	0	{
328	0	switch (blade_in->type)
329	0	{
330	0	case POINTTYPE:
331	0	return lwline_split_by_point(lwline_in, (LWPOINT*)blade_in);
332	0	case MULTIPOINTTYPE:
333	0	return lwline_split_by_mpoint(lwline_in, (LWMPOINT*)blade_in);
334
335	0	case LINETYPE:
336	0	case MULTILINETYPE:
337	0	case POLYGONTYPE:
338	0	case MULTIPOLYGONTYPE:
339	0	return lwline_split_by_line(lwline_in, blade_in);
340
341	0	default:
342	0	lwerror("Splitting a Line by a %s is unsupported",
343	0	lwtype_name(blade_in->type));
344	0	return NULL;
345	0	}
346	0	return NULL;
347	0	}
348
349		/* Initializes and uses GEOS internally */
350		static LWGEOM*
351		lwpoly_split_by_line(const LWPOLY* lwpoly_in, const LWGEOM* blade_in)
352	0	{
353	0	LWCOLLECTION* out;
354	0	GEOSGeometry* g1;
355	0	GEOSGeometry* g2;
356	0	GEOSGeometry* g1_bounds;
357	0	GEOSGeometry* polygons;
358	0	const GEOSGeometry *vgeoms[1];
359	0	int i,n;
360	0	int hasZ = FLAGS_GET_Z(lwpoly_in->flags);
361
362
363		/* Possible outcomes:
364		*
365		* 1. The line does not split the polygon
366		* -> Return a collection with single element
367		* 2. The line does split the polygon
368		* -> Return a collection of all elements resulting from the split
369		*/
370
371	0	initGEOS(lwgeom_geos_error, lwgeom_geos_error);
372
373	0	g1 = LWGEOM2GEOS((LWGEOM*)lwpoly_in, 0);
374	0	if ( NULL == g1 )
375	0	{
376	0	lwerror("LWGEOM2GEOS: %s", lwgeom_geos_errmsg);
377	0	return NULL;
378	0	}
379	0	g1_bounds = GEOSBoundary(g1);
380	0	if ( NULL == g1_bounds )
381	0	{
382	0	GEOSGeom_destroy(g1);
383	0	lwerror("GEOSBoundary: %s", lwgeom_geos_errmsg);
384	0	return NULL;
385	0	}
386
387	0	g2 = LWGEOM2GEOS(blade_in, 0);
388	0	if ( NULL == g2 )
389	0	{
390	0	GEOSGeom_destroy(g1);
391	0	GEOSGeom_destroy(g1_bounds);
392	0	lwerror("LWGEOM2GEOS: %s", lwgeom_geos_errmsg);
393	0	return NULL;
394	0	}
395
396	0	vgeoms[0] = GEOSUnion(g1_bounds, g2);
397	0	if ( NULL == vgeoms[0] )
398	0	{
399	0	GEOSGeom_destroy(g1);
400	0	GEOSGeom_destroy(g2);
401	0	GEOSGeom_destroy(g1_bounds);
402	0	lwerror("GEOSUnion: %s", lwgeom_geos_errmsg);
403	0	return NULL;
404	0	}
405
406	0	polygons = GEOSPolygonize(vgeoms, 1);
407	0	if ( NULL == polygons )
408	0	{
409	0	GEOSGeom_destroy(g1);
410	0	GEOSGeom_destroy(g2);
411	0	GEOSGeom_destroy(g1_bounds);
412	0	GEOSGeom_destroy((GEOSGeometry*)vgeoms[0]);
413	0	lwerror("GEOSPolygonize: %s", lwgeom_geos_errmsg);
414	0	return NULL;
415	0	}
416
417		#ifndef NDEBUG
418		if ( GEOSGeomTypeId(polygons) != COLLECTIONTYPE )
419		{
420		GEOSGeom_destroy(g1);
421		GEOSGeom_destroy(g2);
422		GEOSGeom_destroy(g1_bounds);
423		GEOSGeom_destroy((GEOSGeometry*)vgeoms[0]);
424		GEOSGeom_destroy(polygons);
425		lwerror("%s [%d] Unexpected return from GEOSpolygonize", __FILE__, __LINE__);
426		return 0;
427		}
428		#endif
429
430		/* We should now have all polygons, just skip
431		* the ones which are in holes of the original
432		* geometries and return the rest in a collection
433		*/
434	0	n = GEOSGetNumGeometries(polygons);
435	0	out = lwcollection_construct_empty(COLLECTIONTYPE, lwpoly_in->srid,
436	0	hasZ, 0);
437		/* Allocate space for all polys */
438	0	out->geoms = lwrealloc(out->geoms, sizeof(LWGEOM)n);
439	0	assert(0 == out->ngeoms);
440	0	for (i=0; i<n; ++i)
441	0	{
442	0	GEOSGeometry* pos; /* point on surface */
443	0	const GEOSGeometry* p = GEOSGetGeometryN(polygons, i);
444	0	int contains;
445
446	0	pos = GEOSPointOnSurface(p);
447	0	if ( ! pos )
448	0	{
449	0	GEOSGeom_destroy(g1);
450	0	GEOSGeom_destroy(g2);
451	0	GEOSGeom_destroy(g1_bounds);
452	0	GEOSGeom_destroy((GEOSGeometry*)vgeoms[0]);
453	0	GEOSGeom_destroy(polygons);
454	0	lwerror("GEOSPointOnSurface: %s", lwgeom_geos_errmsg);
455	0	return NULL;
456	0	}
457
458	0	contains = GEOSContains(g1, pos);
459	0	if ( 2 == contains )
460	0	{
461	0	GEOSGeom_destroy(g1);
462	0	GEOSGeom_destroy(g2);
463	0	GEOSGeom_destroy(g1_bounds);
464	0	GEOSGeom_destroy((GEOSGeometry*)vgeoms[0]);
465	0	GEOSGeom_destroy(polygons);
466	0	GEOSGeom_destroy(pos);
467	0	lwerror("GEOSContains: %s", lwgeom_geos_errmsg);
468	0	return NULL;
469	0	}
470
471	0	GEOSGeom_destroy(pos);
472
473	0	if ( 0 == contains )
474	0	{
475		/* Original geometry doesn't contain
476		* a point in this ring, must be an hole
477		*/
478	0	continue;
479	0	}
480
481	0	out->geoms[out->ngeoms++] = GEOS2LWGEOM(p, hasZ);
482	0	}
483
484	0	GEOSGeom_destroy(g1);
485	0	GEOSGeom_destroy(g2);
486	0	GEOSGeom_destroy(g1_bounds);
487	0	GEOSGeom_destroy((GEOSGeometry*)vgeoms[0]);
488	0	GEOSGeom_destroy(polygons);
489
490	0	return (LWGEOM*)out;
491	0	}
492
493		static LWGEOM*
494		lwcollection_split(const LWCOLLECTION* lwcoll_in, const LWGEOM* blade_in)
495	0	{
496	0	LWGEOM** split_vector=NULL;
497	0	LWCOLLECTION* out;
498	0	size_t split_vector_capacity;
499	0	size_t split_vector_size=0;
500	0	size_t i,j;
501
502	0	split_vector_capacity=8;
503	0	split_vector = lwalloc(split_vector_capacity * sizeof(LWGEOM*));
504	0	if ( ! split_vector )
505	0	{
506	0	lwerror("Out of virtual memory");
507	0	return NULL;
508	0	}
509
510	0	for (i=0; i<lwcoll_in->ngeoms; ++i)
511	0	{
512	0	LWCOLLECTION* col;
513	0	LWGEOM* split = lwgeom_split(lwcoll_in->geoms[i], blade_in);
514		/* an exception should prevent this from ever returning NULL */
515	0	if ( ! split ) return NULL;
516
517	0	col = lwgeom_as_lwcollection(split);
518		/* Output, if any, will always be a collection */
519	0	assert(col);
520
521		/* Reallocate split_vector if needed */
522	0	if ( split_vector_size + col->ngeoms > split_vector_capacity )
523	0	{
524		/* NOTE: we could be smarter on reallocations here */
525	0	split_vector_capacity += col->ngeoms;
526	0	split_vector = lwrealloc(split_vector,
527	0	split_vector_capacity * sizeof(LWGEOM*));
528	0	if ( ! split_vector )
529	0	{
530	0	lwerror("Out of virtual memory");
531	0	return NULL;
532	0	}
533	0	}
534
535	0	for (j=0; j<col->ngeoms; ++j)
536	0	{
537	0	col->geoms[j]->srid = SRID_UNKNOWN; /* strip srid */
538	0	split_vector[split_vector_size++] = col->geoms[j];
539	0	}
540	0	lwfree(col->geoms);
541	0	lwfree(col);
542	0	}
543
544		/* Now split_vector has split_vector_size geometries */
545	0	out = lwcollection_construct(COLLECTIONTYPE, lwcoll_in->srid,
546	0	NULL, split_vector_size, split_vector);
547
548	0	return (LWGEOM*)out;
549	0	}
550
551		static LWGEOM*
552		lwpoly_split(const LWPOLY* lwpoly_in, const LWGEOM* blade_in)
553	0	{
554	0	switch (blade_in->type)
555	0	{
556	0	case MULTILINETYPE:
557	0	case LINETYPE:
558	0	return lwpoly_split_by_line(lwpoly_in, blade_in);
559	0	default:
560	0	lwerror("Splitting a Polygon by a %s is unsupported",
561	0	lwtype_name(blade_in->type));
562	0	return NULL;
563	0	}
564	0	return NULL;
565	0	}
566
567		/* exported */
568		LWGEOM*
569		lwgeom_split(const LWGEOM* lwgeom_in, const LWGEOM* blade_in)
570	0	{
571	0	switch (lwgeom_in->type)
572	0	{
573	0	case LINETYPE:
574	0	return lwline_split((const LWLINE*)lwgeom_in, blade_in);
575
576	0	case POLYGONTYPE:
577	0	return lwpoly_split((const LWPOLY*)lwgeom_in, blade_in);
578
579	0	case MULTIPOLYGONTYPE:
580	0	case MULTILINETYPE:
581	0	case COLLECTIONTYPE:
582	0	return lwcollection_split((const LWCOLLECTION*)lwgeom_in, blade_in);
583
584	0	default:
585	0	lwerror("Splitting of %s geometries is unsupported",
586	0	lwtype_name(lwgeom_in->type));
587	0	return NULL;
588	0	}
589
590	0	}
591