PostGIS  2.2.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 4991 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().

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

Here is the call graph for this function:

Here is the caller graph for this function: