PostGIS  2.5.0dev-r@@SVN_REVISION@@
static LWT_ELEMID _lwt_AddPoint ( LWT_TOPOLOGY topo,
LWPOINT point,
double  tol,
int  findFace,
int *  moved 
)
static

Definition at line 5092 of file lwgeom_topo.c.

References _lwt_AddIsoNode(), _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_same(), lwgeom_to_wkt(), lwline_as_lwgeom(), lwline_free(), lwnotice(), lwpoint_as_lwgeom(), 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(), and lwt_AddPoint().

5094 {
5095  int num, i;
5096  double mindist = FLT_MAX;
5097  LWT_ISO_NODE *nodes, *nodes2;
5098  LWT_ISO_EDGE *edges, *edges2;
5099  LWGEOM *pt = lwpoint_as_lwgeom(point);
5100  int flds;
5101  LWT_ELEMID id = 0;
5102  scored_pointer *sorted;
5103 
5104  /* Get tolerance, if 0 was given */
5105  if ( ! tol ) tol = _LWT_MINTOLERANCE( topo, pt );
5106 
5107  LWDEBUGG(1, pt, "Adding point");
5108 
5109  /*
5110  -- 1. Check if any existing node is closer than the given precision
5111  -- and if so pick the closest
5112  TODO: use WithinBox2D
5113  */
5115  nodes = lwt_be_getNodeWithinDistance2D(topo, point, tol, &num, flds, 0);
5116  if ( num == -1 )
5117  {
5118  lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
5119  return -1;
5120  }
5121  if ( num )
5122  {
5123  LWDEBUGF(1, "New point is within %.15g units of %d nodes", tol, num);
5124  /* Order by distance if there are more than a single return */
5125  if ( num > 1 )
5126  {{
5127  sorted= lwalloc(sizeof(scored_pointer)*num);
5128  for (i=0; i<num; ++i)
5129  {
5130  sorted[i].ptr = nodes+i;
5131  sorted[i].score = lwgeom_mindistance2d(lwpoint_as_lwgeom(nodes[i].geom), pt);
5132  LWDEBUGF(1, "Node %" LWTFMT_ELEMID " distance: %.15g",
5133  ((LWT_ISO_NODE*)(sorted[i].ptr))->node_id, sorted[i].score);
5134  }
5135  qsort(sorted, num, sizeof(scored_pointer), compare_scored_pointer);
5136  nodes2 = lwalloc(sizeof(LWT_ISO_NODE)*num);
5137  for (i=0; i<num; ++i)
5138  {
5139  nodes2[i] = *((LWT_ISO_NODE*)sorted[i].ptr);
5140  }
5141  lwfree(sorted);
5142  lwfree(nodes);
5143  nodes = nodes2;
5144  }}
5145 
5146  for ( i=0; i<num; ++i )
5147  {
5148  LWT_ISO_NODE *n = &(nodes[i]);
5149  LWGEOM *g = lwpoint_as_lwgeom(n->geom);
5150  double dist = lwgeom_mindistance2d(g, pt);
5151  /* TODO: move this check in the previous sort scan ... */
5152  /* must be closer than tolerated, unless distance is zero */
5153  if ( dist && dist >= tol ) continue;
5154  if ( ! id || dist < mindist )
5155  {
5156  id = n->node_id;
5157  mindist = dist;
5158  }
5159  }
5160  if ( id )
5161  {
5162  /* found an existing node */
5163  if ( nodes ) _lwt_release_nodes(nodes, num);
5164  if ( moved ) *moved = mindist == 0 ? 0 : 1;
5165  return id;
5166  }
5167  }
5168 
5169  initGEOS(lwnotice, lwgeom_geos_error);
5170 
5171  /*
5172  -- 2. Check if any existing edge falls within tolerance
5173  -- and if so split it by a point projected on it
5174  TODO: use WithinBox2D
5175  */
5177  edges = lwt_be_getEdgeWithinDistance2D(topo, point, tol, &num, flds, 0);
5178  if ( num == -1 )
5179  {
5180  lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
5181  return -1;
5182  }
5183  if ( num )
5184  {
5185  LWDEBUGF(1, "New point is within %.15g units of %d edges", tol, num);
5186 
5187  /* Order by distance if there are more than a single return */
5188  if ( num > 1 )
5189  {{
5190  int j;
5191  sorted = lwalloc(sizeof(scored_pointer)*num);
5192  for (i=0; i<num; ++i)
5193  {
5194  sorted[i].ptr = edges+i;
5195  sorted[i].score = lwgeom_mindistance2d(lwline_as_lwgeom(edges[i].geom), pt);
5196  LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " distance: %.15g",
5197  ((LWT_ISO_EDGE*)(sorted[i].ptr))->edge_id, sorted[i].score);
5198  }
5199  qsort(sorted, num, sizeof(scored_pointer), compare_scored_pointer);
5200  edges2 = lwalloc(sizeof(LWT_ISO_EDGE)*num);
5201  for (j=0, i=0; i<num; ++i)
5202  {
5203  if ( sorted[i].score == sorted[0].score )
5204  {
5205  edges2[j++] = *((LWT_ISO_EDGE*)sorted[i].ptr);
5206  }
5207  else
5208  {
5209  lwline_free(((LWT_ISO_EDGE*)sorted[i].ptr)->geom);
5210  }
5211  }
5212  num = j;
5213  lwfree(sorted);
5214  lwfree(edges);
5215  edges = edges2;
5216  }}
5217 
5218  for (i=0; i<num; ++i)
5219  {
5220  /* The point is on or near an edge, split the edge */
5221  LWT_ISO_EDGE *e = &(edges[i]);
5222  LWGEOM *g = lwline_as_lwgeom(e->geom);
5223  LWGEOM *prj;
5224  int contains;
5225  GEOSGeometry *prjg, *gg;
5226  LWT_ELEMID edge_id = e->edge_id;
5227 
5228  LWDEBUGF(1, "Splitting edge %" LWTFMT_ELEMID, edge_id);
5229 
5230  /* project point to line, split edge by point */
5231  prj = lwgeom_closest_point(g, pt);
5232  if ( moved ) *moved = lwgeom_same(prj,pt) ? 0 : 1;
5233  if ( lwgeom_has_z(pt) )
5234  {{
5235  /*
5236  -- This is a workaround for ClosestPoint lack of Z support:
5237  -- http://trac.osgeo.org/postgis/ticket/2033
5238  */
5239  LWGEOM *tmp;
5240  double z;
5241  POINT4D p4d;
5242  LWPOINT *prjpt;
5243  /* add Z to "prj" */
5244  tmp = lwgeom_force_3dz(prj);
5245  prjpt = lwgeom_as_lwpoint(tmp);
5246  getPoint4d_p(point->point, 0, &p4d);
5247  z = p4d.z;
5248  getPoint4d_p(prjpt->point, 0, &p4d);
5249  p4d.z = z;
5250  ptarray_set_point4d(prjpt->point, 0, &p4d);
5251  lwgeom_free(prj);
5252  prj = tmp;
5253  }}
5254  prjg = LWGEOM2GEOS(prj, 0);
5255  if ( ! prjg ) {
5256  lwgeom_free(prj);
5257  _lwt_release_edges(edges, num);
5258  lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg);
5259  return -1;
5260  }
5261  gg = LWGEOM2GEOS(g, 0);
5262  if ( ! gg ) {
5263  lwgeom_free(prj);
5264  _lwt_release_edges(edges, num);
5265  GEOSGeom_destroy(prjg);
5266  lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg);
5267  return -1;
5268  }
5269  contains = GEOSContains(gg, prjg);
5270  GEOSGeom_destroy(prjg);
5271  GEOSGeom_destroy(gg);
5272  if ( contains == 2 )
5273  {
5274  lwgeom_free(prj);
5275  _lwt_release_edges(edges, num);
5276  lwerror("GEOS exception on Contains: %s", lwgeom_geos_errmsg);
5277  return -1;
5278  }
5279  if ( ! contains )
5280  {{
5281  double snaptol;
5282  LWGEOM *snapedge;
5283  LWLINE *snapline;
5284  POINT4D p1, p2;
5285 
5286  LWDEBUGF(1, "Edge %" LWTFMT_ELEMID
5287  " does not contain projected point to it",
5288  edge_id);
5289 
5290  /* In order to reduce the robustness issues, we'll pick
5291  * an edge that contains the projected point, if possible */
5292  if ( i+1 < num )
5293  {
5294  LWDEBUG(1, "But there's another to check");
5295  lwgeom_free(prj);
5296  continue;
5297  }
5298 
5299  /*
5300  -- The tolerance must be big enough for snapping to happen
5301  -- and small enough to snap only to the projected point.
5302  -- Unfortunately ST_Distance returns 0 because it also uses
5303  -- a projected point internally, so we need another way.
5304  */
5305  snaptol = _lwt_minTolerance(prj);
5306  snapedge = _lwt_toposnap(g, prj, snaptol);
5307  snapline = lwgeom_as_lwline(snapedge);
5308 
5309  LWDEBUGF(1, "Edge snapped with tolerance %g", snaptol);
5310 
5311  /* TODO: check if snapping did anything ? */
5312 #if POSTGIS_DEBUG_LEVEL > 0
5313  {
5314  size_t sz;
5315  char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz);
5316  char *wkt2 = lwgeom_to_wkt(snapedge, WKT_EXTENDED, 15, &sz);
5317  LWDEBUGF(1, "Edge %s snapped became %s", wkt1, wkt2);
5318  lwfree(wkt1);
5319  lwfree(wkt2);
5320  }
5321 #endif
5322 
5323 
5324  /*
5325  -- Snapping currently snaps the first point below tolerance
5326  -- so may possibly move first point. See ticket #1631
5327  */
5328  getPoint4d_p(e->geom->points, 0, &p1);
5329  getPoint4d_p(snapline->points, 0, &p2);
5330  LWDEBUGF(1, "Edge first point is %g %g, "
5331  "snapline first point is %g %g",
5332  p1.x, p1.y, p2.x, p2.y);
5333  if ( p1.x != p2.x || p1.y != p2.y )
5334  {
5335  LWDEBUG(1, "Snapping moved first point, re-adding it");
5336  if ( LW_SUCCESS != ptarray_insert_point(snapline->points, &p1, 0) )
5337  {
5338  lwgeom_free(prj);
5339  lwgeom_free(snapedge);
5340  _lwt_release_edges(edges, num);
5341  lwerror("GEOS exception on Contains: %s", lwgeom_geos_errmsg);
5342  return -1;
5343  }
5344 #if POSTGIS_DEBUG_LEVEL > 0
5345  {
5346  size_t sz;
5347  char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz);
5348  LWDEBUGF(1, "Tweaked snapline became %s", wkt1);
5349  lwfree(wkt1);
5350  }
5351 #endif
5352  }
5353 #if POSTGIS_DEBUG_LEVEL > 0
5354  else {
5355  LWDEBUG(1, "Snapping did not move first point");
5356  }
5357 #endif
5358 
5359  if ( -1 == lwt_ChangeEdgeGeom( topo, edge_id, snapline ) )
5360  {
5361  /* TODO: should have invoked lwerror already, leaking memory */
5362  lwgeom_free(prj);
5363  lwgeom_free(snapedge);
5364  _lwt_release_edges(edges, num);
5365  lwerror("lwt_ChangeEdgeGeom failed");
5366  return -1;
5367  }
5368  lwgeom_free(snapedge);
5369  }}
5370 #if POSTGIS_DEBUG_LEVEL > 0
5371  else
5372  {{
5373  size_t sz;
5374  char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz);
5375  char *wkt2 = lwgeom_to_wkt(prj, WKT_EXTENDED, 15, &sz);
5376  LWDEBUGF(1, "Edge %s contains projected point %s", wkt1, wkt2);
5377  lwfree(wkt1);
5378  lwfree(wkt2);
5379  }}
5380 #endif
5381 
5382  /* TODO: pass 1 as last argument (skipChecks) ? */
5383  id = lwt_ModEdgeSplit( topo, edge_id, lwgeom_as_lwpoint(prj), 0 );
5384  if ( -1 == id )
5385  {
5386  /* TODO: should have invoked lwerror already, leaking memory */
5387  lwgeom_free(prj);
5388  _lwt_release_edges(edges, num);
5389  lwerror("lwt_ModEdgeSplit failed");
5390  return -1;
5391  }
5392 
5393  lwgeom_free(prj);
5394 
5395  /*
5396  * TODO: decimate the two new edges with the given tolerance ?
5397  *
5398  * the edge identifiers to decimate would be: edge_id and "id"
5399  * The problem here is that decimation of existing edges
5400  * may introduce intersections or topological inconsistencies,
5401  * for example:
5402  *
5403  * - A node may end up falling on the other side of the edge
5404  * - The decimated edge might intersect another existing edge
5405  *
5406  */
5407 
5408  break; /* we only want to snap a single edge */
5409  }
5410  _lwt_release_edges(edges, num);
5411  }
5412  else
5413  {
5414  /* The point is isolated, add it as such */
5415  /* TODO: pass 1 as last argument (skipChecks) ? */
5416  id = _lwt_AddIsoNode(topo, -1, point, 0, findFace);
5417  if ( moved ) *moved = 0;
5418  if ( -1 == id )
5419  {
5420  /* should have invoked lwerror already, leaking memory */
5421  lwerror("lwt_AddIsoNode failed");
5422  return -1;
5423  }
5424  }
5425 
5426  return id;
5427 }
double x
Definition: liblwgeom.h:351
static void _lwt_release_nodes(LWT_ISO_NODE *nodes, int num_nodes)
Definition: lwgeom_topo.c:467
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:3331
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:675
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:177
void lwfree(void *mem)
Definition: lwutil.c:244
LWGEOM * lwgeom_force_3dz(const LWGEOM *geom)
Definition: lwgeom.c:783
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1137
static int compare_scored_pointer(const void *si1, const void *si2)
Definition: lwgeom_topo.c:5072
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]
void ptarray_set_point4d(POINTARRAY *pa, uint32_t n, const POINT4D *p4d)
Definition: lwgeom_api.c:425
LWLINE * geom
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:83
LWPOINT * lwgeom_as_lwpoint(const LWGEOM *lwgeom)
Definition: lwgeom.c:160
POINTARRAY * point
Definition: liblwgeom.h:410
#define _LWT_MINTOLERANCE(topo, geom)
Definition: lwgeom_topo.c:5063
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:923
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:1046
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:329
void lwgeom_geos_error(const char *fmt,...)
const LWT_BE_IFACE * be_iface
static double _lwt_minTolerance(LWGEOM *g)
Definition: lwgeom_topo.c:5045
LWT_ELEMID node_id
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:169
#define LWDEBUGG(level, geom, msg)
Definition: lwgeom_log.h:93
char lwgeom_same(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2)
geom1 same as geom2 iff
Definition: lwgeom.c:575
int getPoint4d_p(const POINTARRAY *pa, uint32_t n, POINT4D *point)
Definition: lwgeom_api.c:113
#define LWT_COL_EDGE_EDGE_ID
Edge fields.
#define LWT_COL_NODE_NODE_ID
Node fields.
#define WKT_EXTENDED
Definition: liblwgeom.h:2070
static void _lwt_release_edges(LWT_ISO_EDGE *edges, int num_edges)
Definition: lwgeom_topo.c:457
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
int ptarray_insert_point(POINTARRAY *pa, const POINT4D *p, uint32_t where)
Insert a point into an existing POINTARRAY.
Definition: ptarray.c:96
#define LWT_COL_NODE_GEOM
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:334
LWT_INT64 LWT_ELEMID
Identifier of topology element.
#define LWT_COL_EDGE_GEOM
void * lwalloc(size_t size)
Definition: lwutil.c:229
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
uint32_t e
Definition: geobuf.h:57
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:190
#define LWTFMT_ELEMID
Definition: lwgeom_topo.c:44
POINTARRAY * points
Definition: liblwgeom.h:421
static LWT_ELEMID _lwt_AddIsoNode(LWT_TOPOLOGY *topo, LWT_ELEMID face, LWPOINT *pt, int skipISOChecks, int checkFace)
Definition: lwgeom_topo.c:523

Here is the call graph for this function:

Here is the caller graph for this function: