PostGIS  2.3.7dev-r@@SVN_REVISION@@
LWT_ELEMID lwt_AddPoint ( LWT_TOPOLOGY topo,
LWPOINT point,
double  tol 
)

Adds a point to the topology.

The given point will snap to existing nodes or edges within given tolerance. An existing edge may be split by the point.

Parameters
topothe topology to operate on
pointthe point to add
tolsnap tolerance, the topology tolerance will be used if 0
Returns
identifier of added (or pre-existing) node or -1 on error (liblwgeom error handler will be invoked with error message)

Definition at line 5021 of file lwgeom_topo.c.

References _lwt_minTolerance(), _LWT_MINTOLERANCE, _lwt_release_edges(), _lwt_release_nodes(), _lwt_toposnap(), LWT_TOPOLOGY_T::be_iface, compare_scored_pointer(), contains(), LWT_ISO_EDGE::edge_id, LWT_ISO_NODE::geom, LWT_ISO_EDGE::geom, getPoint4d_p(), LW_SUCCESS, lwalloc(), LWDEBUG, LWDEBUGF, LWDEBUGG, lwerror(), lwfree(), LWGEOM2GEOS(), lwgeom_as_lwline(), lwgeom_as_lwpoint(), lwgeom_closest_point(), lwgeom_force_3dz(), lwgeom_free(), lwgeom_geos_errmsg, lwgeom_geos_error(), lwgeom_has_z(), lwgeom_mindistance2d(), lwgeom_to_wkt(), lwline_as_lwgeom(), lwline_free(), lwnotice(), lwpoint_as_lwgeom(), lwt_AddIsoNode(), lwt_be_getEdgeWithinDistance2D(), lwt_be_getNodeWithinDistance2D(), lwt_be_lastErrorMessage(), lwt_ChangeEdgeGeom(), LWT_COL_EDGE_EDGE_ID, LWT_COL_EDGE_GEOM, LWT_COL_NODE_GEOM, LWT_COL_NODE_NODE_ID, lwt_ModEdgeSplit(), LWTFMT_ELEMID, LWT_ISO_NODE::node_id, LWPOINT::point, LWLINE::points, ptarray_insert_point(), ptarray_set_point4d(), scored_pointer_t::ptr, scored_pointer_t::score, WKT_EXTENDED, POINT4D::x, POINT4D::y, and POINT4D::z.

Referenced by _lwt_AddLineEdge().

5022 {
5023  int num, i;
5024  double mindist = FLT_MAX;
5025  LWT_ISO_NODE *nodes, *nodes2;
5026  LWT_ISO_EDGE *edges, *edges2;
5027  LWGEOM *pt = lwpoint_as_lwgeom(point);
5028  int flds;
5029  LWT_ELEMID id = 0;
5030  scored_pointer *sorted;
5031 
5032  /* Get tolerance, if 0 was given */
5033  if ( ! tol ) tol = _LWT_MINTOLERANCE( topo, pt );
5034 
5035  LWDEBUGG(1, pt, "Adding point");
5036 
5037  /*
5038  -- 1. Check if any existing node is closer than the given precision
5039  -- and if so pick the closest
5040  TODO: use WithinBox2D
5041  */
5043  nodes = lwt_be_getNodeWithinDistance2D(topo, point, tol, &num, flds, 0);
5044  if ( num == -1 )
5045  {
5046  lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
5047  return -1;
5048  }
5049  if ( num )
5050  {
5051  LWDEBUGF(1, "New point is within %.15g units of %d nodes", tol, num);
5052  /* Order by distance if there are more than a single return */
5053  if ( num > 1 )
5054  {{
5055  sorted= lwalloc(sizeof(scored_pointer)*num);
5056  for (i=0; i<num; ++i)
5057  {
5058  sorted[i].ptr = nodes+i;
5059  sorted[i].score = lwgeom_mindistance2d(lwpoint_as_lwgeom(nodes[i].geom), pt);
5060  LWDEBUGF(1, "Node %" LWTFMT_ELEMID " distance: %.15g",
5061  ((LWT_ISO_NODE*)(sorted[i].ptr))->node_id, sorted[i].score);
5062  }
5063  qsort(sorted, num, sizeof(scored_pointer), compare_scored_pointer);
5064  nodes2 = lwalloc(sizeof(LWT_ISO_NODE)*num);
5065  for (i=0; i<num; ++i)
5066  {
5067  nodes2[i] = *((LWT_ISO_NODE*)sorted[i].ptr);
5068  }
5069  lwfree(sorted);
5070  lwfree(nodes);
5071  nodes = nodes2;
5072  }}
5073 
5074  for ( i=0; i<num; ++i )
5075  {
5076  LWT_ISO_NODE *n = &(nodes[i]);
5077  LWGEOM *g = lwpoint_as_lwgeom(n->geom);
5078  double dist = lwgeom_mindistance2d(g, pt);
5079  /* TODO: move this check in the previous sort scan ... */
5080  if ( dist >= tol ) continue; /* must be closer than tolerated */
5081  if ( ! id || dist < mindist )
5082  {
5083  id = n->node_id;
5084  mindist = dist;
5085  }
5086  }
5087  if ( id )
5088  {
5089  /* found an existing node */
5090  if ( nodes ) _lwt_release_nodes(nodes, num);
5091  return id;
5092  }
5093  }
5094 
5095  initGEOS(lwnotice, lwgeom_geos_error);
5096 
5097  /*
5098  -- 2. Check if any existing edge falls within tolerance
5099  -- and if so split it by a point projected on it
5100  TODO: use WithinBox2D
5101  */
5103  edges = lwt_be_getEdgeWithinDistance2D(topo, point, tol, &num, flds, 0);
5104  if ( num == -1 )
5105  {
5106  lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
5107  return -1;
5108  }
5109  if ( num )
5110  {
5111  LWDEBUGF(1, "New point is within %.15g units of %d edges", tol, num);
5112 
5113  /* Order by distance if there are more than a single return */
5114  if ( num > 1 )
5115  {{
5116  int j;
5117  sorted = lwalloc(sizeof(scored_pointer)*num);
5118  for (i=0; i<num; ++i)
5119  {
5120  sorted[i].ptr = edges+i;
5121  sorted[i].score = lwgeom_mindistance2d(lwline_as_lwgeom(edges[i].geom), pt);
5122  LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " distance: %.15g",
5123  ((LWT_ISO_EDGE*)(sorted[i].ptr))->edge_id, sorted[i].score);
5124  }
5125  qsort(sorted, num, sizeof(scored_pointer), compare_scored_pointer);
5126  edges2 = lwalloc(sizeof(LWT_ISO_EDGE)*num);
5127  for (j=0, i=0; i<num; ++i)
5128  {
5129  if ( sorted[i].score == sorted[0].score )
5130  {
5131  edges2[j++] = *((LWT_ISO_EDGE*)sorted[i].ptr);
5132  }
5133  else
5134  {
5135  lwline_free(((LWT_ISO_EDGE*)sorted[i].ptr)->geom);
5136  }
5137  }
5138  num = j;
5139  lwfree(sorted);
5140  lwfree(edges);
5141  edges = edges2;
5142  }}
5143 
5144  for (i=0; i<num; ++i)
5145  {
5146  /* The point is on or near an edge, split the edge */
5147  LWT_ISO_EDGE *e = &(edges[i]);
5148  LWGEOM *g = lwline_as_lwgeom(e->geom);
5149  LWGEOM *prj;
5150  int contains;
5151  GEOSGeometry *prjg, *gg;
5152  LWT_ELEMID edge_id = e->edge_id;
5153 
5154  LWDEBUGF(1, "Splitting edge %" LWTFMT_ELEMID, edge_id);
5155 
5156  /* project point to line, split edge by point */
5157  prj = lwgeom_closest_point(g, pt);
5158  if ( lwgeom_has_z(pt) )
5159  {{
5160  /*
5161  -- This is a workaround for ClosestPoint lack of Z support:
5162  -- http://trac.osgeo.org/postgis/ticket/2033
5163  */
5164  LWGEOM *tmp;
5165  double z;
5166  POINT4D p4d;
5167  LWPOINT *prjpt;
5168  /* add Z to "prj" */
5169  tmp = lwgeom_force_3dz(prj);
5170  prjpt = lwgeom_as_lwpoint(tmp);
5171  getPoint4d_p(point->point, 0, &p4d);
5172  z = p4d.z;
5173  getPoint4d_p(prjpt->point, 0, &p4d);
5174  p4d.z = z;
5175  ptarray_set_point4d(prjpt->point, 0, &p4d);
5176  lwgeom_free(prj);
5177  prj = tmp;
5178  }}
5179  prjg = LWGEOM2GEOS(prj, 0);
5180  if ( ! prjg ) {
5181  lwgeom_free(prj);
5182  _lwt_release_edges(edges, num);
5183  lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg);
5184  return -1;
5185  }
5186  gg = LWGEOM2GEOS(g, 0);
5187  if ( ! gg ) {
5188  lwgeom_free(prj);
5189  _lwt_release_edges(edges, num);
5190  GEOSGeom_destroy(prjg);
5191  lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg);
5192  return -1;
5193  }
5194  contains = GEOSContains(gg, prjg);
5195  GEOSGeom_destroy(prjg);
5196  GEOSGeom_destroy(gg);
5197  if ( contains == 2 )
5198  {
5199  lwgeom_free(prj);
5200  _lwt_release_edges(edges, num);
5201  lwerror("GEOS exception on Contains: %s", lwgeom_geos_errmsg);
5202  return -1;
5203  }
5204  if ( ! contains )
5205  {{
5206  double snaptol;
5207  LWGEOM *snapedge;
5208  LWLINE *snapline;
5209  POINT4D p1, p2;
5210 
5211  LWDEBUGF(1, "Edge %" LWTFMT_ELEMID
5212  " does not contain projected point to it",
5213  edge_id);
5214 
5215  /* In order to reduce the robustness issues, we'll pick
5216  * an edge that contains the projected point, if possible */
5217  if ( i+1 < num )
5218  {
5219  LWDEBUG(1, "But there's another to check");
5220  lwgeom_free(prj);
5221  continue;
5222  }
5223 
5224  /*
5225  -- The tolerance must be big enough for snapping to happen
5226  -- and small enough to snap only to the projected point.
5227  -- Unfortunately ST_Distance returns 0 because it also uses
5228  -- a projected point internally, so we need another way.
5229  */
5230  snaptol = _lwt_minTolerance(prj);
5231  snapedge = _lwt_toposnap(g, prj, snaptol);
5232  snapline = lwgeom_as_lwline(snapedge);
5233 
5234  LWDEBUGF(1, "Edge snapped with tolerance %g", snaptol);
5235 
5236  /* TODO: check if snapping did anything ? */
5237 #if POSTGIS_DEBUG_LEVEL > 0
5238  {
5239  size_t sz;
5240  char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz);
5241  char *wkt2 = lwgeom_to_wkt(snapedge, WKT_EXTENDED, 15, &sz);
5242  LWDEBUGF(1, "Edge %s snapped became %s", wkt1, wkt2);
5243  lwfree(wkt1);
5244  lwfree(wkt2);
5245  }
5246 #endif
5247 
5248 
5249  /*
5250  -- Snapping currently snaps the first point below tolerance
5251  -- so may possibly move first point. See ticket #1631
5252  */
5253  getPoint4d_p(e->geom->points, 0, &p1);
5254  getPoint4d_p(snapline->points, 0, &p2);
5255  LWDEBUGF(1, "Edge first point is %g %g, "
5256  "snapline first point is %g %g",
5257  p1.x, p1.y, p2.x, p2.y);
5258  if ( p1.x != p2.x || p1.y != p2.y )
5259  {
5260  LWDEBUG(1, "Snapping moved first point, re-adding it");
5261  if ( LW_SUCCESS != ptarray_insert_point(snapline->points, &p1, 0) )
5262  {
5263  lwgeom_free(prj);
5264  lwgeom_free(snapedge);
5265  _lwt_release_edges(edges, num);
5266  lwerror("GEOS exception on Contains: %s", lwgeom_geos_errmsg);
5267  return -1;
5268  }
5269 #if POSTGIS_DEBUG_LEVEL > 0
5270  {
5271  size_t sz;
5272  char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz);
5273  LWDEBUGF(1, "Tweaked snapline became %s", wkt1);
5274  lwfree(wkt1);
5275  }
5276 #endif
5277  }
5278 #if POSTGIS_DEBUG_LEVEL > 0
5279  else {
5280  LWDEBUG(1, "Snapping did not move first point");
5281  }
5282 #endif
5283 
5284  if ( -1 == lwt_ChangeEdgeGeom( topo, edge_id, snapline ) )
5285  {
5286  /* TODO: should have invoked lwerror already, leaking memory */
5287  lwgeom_free(prj);
5288  lwgeom_free(snapedge);
5289  _lwt_release_edges(edges, num);
5290  lwerror("lwt_ChangeEdgeGeom failed");
5291  return -1;
5292  }
5293  lwgeom_free(snapedge);
5294  }}
5295 #if POSTGIS_DEBUG_LEVEL > 0
5296  else
5297  {{
5298  size_t sz;
5299  char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz);
5300  char *wkt2 = lwgeom_to_wkt(prj, WKT_EXTENDED, 15, &sz);
5301  LWDEBUGF(1, "Edge %s contains projected point %s", wkt1, wkt2);
5302  lwfree(wkt1);
5303  lwfree(wkt2);
5304  }}
5305 #endif
5306 
5307  /* TODO: pass 1 as last argument (skipChecks) ? */
5308  id = lwt_ModEdgeSplit( topo, edge_id, lwgeom_as_lwpoint(prj), 0 );
5309  if ( -1 == id )
5310  {
5311  /* TODO: should have invoked lwerror already, leaking memory */
5312  lwgeom_free(prj);
5313  _lwt_release_edges(edges, num);
5314  lwerror("lwt_ModEdgeSplit failed");
5315  return -1;
5316  }
5317 
5318  lwgeom_free(prj);
5319 
5320  /*
5321  * TODO: decimate the two new edges with the given tolerance ?
5322  *
5323  * the edge identifiers to decimate would be: edge_id and "id"
5324  * The problem here is that decimation of existing edges
5325  * may introduce intersections or topological inconsistencies,
5326  * for example:
5327  *
5328  * - A node may end up falling on the other side of the edge
5329  * - The decimated edge might intersect another existing edge
5330  *
5331  */
5332 
5333  break; /* we only want to snap a single edge */
5334  }
5335  _lwt_release_edges(edges, num);
5336  }
5337  else
5338  {
5339  /* The point is isolated, add it as such */
5340  /* TODO: pass 1 as last argument (skipChecks) ? */
5341  id = lwt_AddIsoNode(topo, -1, point, 0);
5342  if ( -1 == id )
5343  {
5344  /* should have invoked lwerror already, leaking memory */
5345  lwerror("lwt_AddIsoNode failed");
5346  return -1;
5347  }
5348  }
5349 
5350  return id;
5351 }
void ptarray_set_point4d(POINTARRAY *pa, int n, const POINT4D *p4d)
Definition: lwgeom_api.c:549
double x
Definition: liblwgeom.h:351
static void _lwt_release_nodes(LWT_ISO_NODE *nodes, int num_nodes)
Definition: lwgeom_topo.c:477
int lwt_ChangeEdgeGeom(LWT_TOPOLOGY *topo, LWT_ELEMID edge_id, LWLINE *geom)
Changes the shape of an edge without affecting the topology structure.
Definition: lwgeom_topo.c:3298
LWGEOM * lwgeom_closest_point(const LWGEOM *lw1, const LWGEOM *lw2)
Definition: measures.c:54
LWPOINT * geom
char * lwgeom_to_wkt(const LWGEOM *geom, uint8_t variant, int precision, size_t *size_out)
WKT emitter function.
Definition: lwout_wkt.c:669
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:89
void lwfree(void *mem)
Definition: lwutil.c:242
LWGEOM * lwgeom_force_3dz(const LWGEOM *geom)
Definition: lwgeom.c:709
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1063
static int compare_scored_pointer(const void *si1, const void *si2)
Definition: lwgeom_topo.c:5008
LWT_ISO_NODE * lwt_be_getNodeWithinDistance2D(LWT_TOPOLOGY *topo, LWPOINT *pt, double dist, int *numelems, int fields, int limit)
Definition: lwgeom_topo.c:163
#define LW_SUCCESS
Definition: liblwgeom.h:79
void lwline_free(LWLINE *line)
Definition: lwline.c:76
Datum contains(PG_FUNCTION_ARGS)
char lwgeom_geos_errmsg[LWGEOM_GEOS_ERRMSG_MAXSIZE]
LWLINE * geom
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
LWPOINT * lwgeom_as_lwpoint(const LWGEOM *lwgeom)
Definition: lwgeom.c:93
POINTARRAY * point
Definition: liblwgeom.h:410
#define _LWT_MINTOLERANCE(topo, geom)
Definition: lwgeom_topo.c:4999
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:849
LWT_ELEMID lwt_ModEdgeSplit(LWT_TOPOLOGY *topo, LWT_ELEMID edge, LWPOINT *pt, int skipISOChecks)
Split an edge by a node, modifying the original edge and adding a new one.
Definition: lwgeom_topo.c:1034
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:262
void lwgeom_geos_error(const char *fmt,...)
const LWT_BE_IFACE * be_iface
static double _lwt_minTolerance(LWGEOM *g)
Definition: lwgeom_topo.c:4981
LWT_ELEMID node_id
int ptarray_insert_point(POINTARRAY *pa, const POINT4D *p, int where)
Insert a point into an existing POINTARRAY.
Definition: ptarray.c:96
LWT_ELEMID edge_id
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, int autofix)
double z
Definition: liblwgeom.h:351
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition: lwgeom.c:102
LWT_ELEMID lwt_AddIsoNode(LWT_TOPOLOGY *topo, LWT_ELEMID face, LWPOINT *pt, int skipISOChecks)
Add an isolated node.
Definition: lwgeom_topo.c:525
#define LWDEBUGG(level, geom, msg)
Definition: lwgeom_log.h:93
#define LWT_COL_EDGE_EDGE_ID
Edge fields.
#define LWT_COL_NODE_NODE_ID
Node fields.
#define WKT_EXTENDED
Definition: liblwgeom.h:2057
static void _lwt_release_edges(LWT_ISO_EDGE *edges, int num_edges)
Definition: lwgeom_topo.c:467
double lwgeom_mindistance2d(const LWGEOM *lw1, const LWGEOM *lw2)
Function initialazing min distance calculation.
Definition: measures.c:202
LWT_ISO_EDGE * lwt_be_getEdgeWithinDistance2D(LWT_TOPOLOGY *topo, LWPOINT *pt, double dist, int *numelems, int fields, int limit)
Definition: lwgeom_topo.c:260
#define LWT_COL_NODE_GEOM
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:267
LWT_INT64 LWT_ELEMID
Identifier of topology element.
#define LWT_COL_EDGE_GEOM
void * lwalloc(size_t size)
Definition: lwutil.c:227
double y
Definition: liblwgeom.h:351
static LWGEOM * _lwt_toposnap(LWGEOM *src, LWGEOM *tgt, double tol)
Definition: lwgeom_topo.c:420
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
const char * lwt_be_lastErrorMessage(const LWT_BE_IFACE *be)
Definition: lwgeom_topo.c:120
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:102
int getPoint4d_p(const POINTARRAY *pa, int n, POINT4D *point)
Definition: lwgeom_api.c:244
#define LWTFMT_ELEMID
Definition: lwgeom_topo.c:44
POINTARRAY * points
Definition: liblwgeom.h:421

Here is the call graph for this function:

Here is the caller graph for this function: