PostGIS  2.5.0beta1dev-r@@SVN_REVISION@@

◆ _lwt_AddPoint()

static LWT_ELEMID _lwt_AddPoint ( LWT_TOPOLOGY topo,
LWPOINT point,
double  tol,
int  findFace,
int *  moved 
)
static

Definition at line 4886 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, getPoint2d_cp(), getPoint4d_p(), LW_BOUNDARY, LW_SUCCESS, lwalloc(), LWDEBUG, LWDEBUGF, LWDEBUGG, lwerror(), lwfree(), 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_contains_point_partial(), 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().

4888 {
4889  int num, i;
4890  double mindist = FLT_MAX;
4891  LWT_ISO_NODE *nodes, *nodes2;
4892  LWT_ISO_EDGE *edges, *edges2;
4893  LWGEOM *pt = lwpoint_as_lwgeom(point);
4894  int flds;
4895  LWT_ELEMID id = 0;
4896  scored_pointer *sorted;
4897 
4898  /* Get tolerance, if 0 was given */
4899  if ( ! tol ) tol = _LWT_MINTOLERANCE( topo, pt );
4900 
4901  LWDEBUGG(1, pt, "Adding point");
4902 
4903  /*
4904  -- 1. Check if any existing node is closer than the given precision
4905  -- and if so pick the closest
4906  TODO: use WithinBox2D
4907  */
4909  nodes = lwt_be_getNodeWithinDistance2D(topo, point, tol, &num, flds, 0);
4910  if ( num == -1 )
4911  {
4912  lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4913  return -1;
4914  }
4915  if ( num )
4916  {
4917  LWDEBUGF(1, "New point is within %.15g units of %d nodes", tol, num);
4918  /* Order by distance if there are more than a single return */
4919  if ( num > 1 )
4920  {{
4921  sorted= lwalloc(sizeof(scored_pointer)*num);
4922  for (i=0; i<num; ++i)
4923  {
4924  sorted[i].ptr = nodes+i;
4925  sorted[i].score = lwgeom_mindistance2d(lwpoint_as_lwgeom(nodes[i].geom), pt);
4926  LWDEBUGF(1, "Node %" LWTFMT_ELEMID " distance: %.15g",
4927  ((LWT_ISO_NODE*)(sorted[i].ptr))->node_id, sorted[i].score);
4928  }
4929  qsort(sorted, num, sizeof(scored_pointer), compare_scored_pointer);
4930  nodes2 = lwalloc(sizeof(LWT_ISO_NODE)*num);
4931  for (i=0; i<num; ++i)
4932  {
4933  nodes2[i] = *((LWT_ISO_NODE*)sorted[i].ptr);
4934  }
4935  lwfree(sorted);
4936  lwfree(nodes);
4937  nodes = nodes2;
4938  }}
4939 
4940  for ( i=0; i<num; ++i )
4941  {
4942  LWT_ISO_NODE *n = &(nodes[i]);
4943  LWGEOM *g = lwpoint_as_lwgeom(n->geom);
4944  double dist = lwgeom_mindistance2d(g, pt);
4945  /* TODO: move this check in the previous sort scan ... */
4946  /* must be closer than tolerated, unless distance is zero */
4947  if ( dist && dist >= tol ) continue;
4948  if ( ! id || dist < mindist )
4949  {
4950  id = n->node_id;
4951  mindist = dist;
4952  }
4953  }
4954  if ( id )
4955  {
4956  /* found an existing node */
4957  if ( nodes ) _lwt_release_nodes(nodes, num);
4958  if ( moved ) *moved = mindist == 0 ? 0 : 1;
4959  return id;
4960  }
4961  }
4962 
4963  initGEOS(lwnotice, lwgeom_geos_error);
4964 
4965  /*
4966  -- 2. Check if any existing edge falls within tolerance
4967  -- and if so split it by a point projected on it
4968  TODO: use WithinBox2D
4969  */
4971  edges = lwt_be_getEdgeWithinDistance2D(topo, point, tol, &num, flds, 0);
4972  if ( num == -1 )
4973  {
4974  lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
4975  return -1;
4976  }
4977  if ( num )
4978  {
4979  LWDEBUGF(1, "New point is within %.15g units of %d edges", tol, num);
4980 
4981  /* Order by distance if there are more than a single return */
4982  if ( num > 1 )
4983  {{
4984  int j;
4985  sorted = lwalloc(sizeof(scored_pointer)*num);
4986  for (i=0; i<num; ++i)
4987  {
4988  sorted[i].ptr = edges+i;
4989  sorted[i].score = lwgeom_mindistance2d(lwline_as_lwgeom(edges[i].geom), pt);
4990  LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " distance: %.15g",
4991  ((LWT_ISO_EDGE*)(sorted[i].ptr))->edge_id, sorted[i].score);
4992  }
4993  qsort(sorted, num, sizeof(scored_pointer), compare_scored_pointer);
4994  edges2 = lwalloc(sizeof(LWT_ISO_EDGE)*num);
4995  for (j=0, i=0; i<num; ++i)
4996  {
4997  if ( sorted[i].score == sorted[0].score )
4998  {
4999  edges2[j++] = *((LWT_ISO_EDGE*)sorted[i].ptr);
5000  }
5001  else
5002  {
5003  lwline_free(((LWT_ISO_EDGE*)sorted[i].ptr)->geom);
5004  }
5005  }
5006  num = j;
5007  lwfree(sorted);
5008  lwfree(edges);
5009  edges = edges2;
5010  }}
5011 
5012  for (i=0; i<num; ++i)
5013  {
5014  /* The point is on or near an edge, split the edge */
5015  LWT_ISO_EDGE *e = &(edges[i]);
5016  LWGEOM *g = lwline_as_lwgeom(e->geom);
5017  LWGEOM *prj;
5018  int contains;
5019  LWT_ELEMID edge_id = e->edge_id;
5020 
5021  LWDEBUGF(1, "Splitting edge %" LWTFMT_ELEMID, edge_id);
5022 
5023  /* project point to line, split edge by point */
5024  prj = lwgeom_closest_point(g, pt);
5025  if ( moved ) *moved = lwgeom_same(prj,pt) ? 0 : 1;
5026  if ( lwgeom_has_z(pt) )
5027  {{
5028  /*
5029  -- This is a workaround for ClosestPoint lack of Z support:
5030  -- http://trac.osgeo.org/postgis/ticket/2033
5031  */
5032  LWGEOM *tmp;
5033  double z;
5034  POINT4D p4d;
5035  LWPOINT *prjpt;
5036  /* add Z to "prj" */
5037  tmp = lwgeom_force_3dz(prj);
5038  prjpt = lwgeom_as_lwpoint(tmp);
5039  getPoint4d_p(point->point, 0, &p4d);
5040  z = p4d.z;
5041  getPoint4d_p(prjpt->point, 0, &p4d);
5042  p4d.z = z;
5043  ptarray_set_point4d(prjpt->point, 0, &p4d);
5044  lwgeom_free(prj);
5045  prj = tmp;
5046  }}
5047  const POINT2D *pt = getPoint2d_cp(lwgeom_as_lwpoint(prj)->point, 0);
5048  contains = ptarray_contains_point_partial(e->geom->points, pt, 0, NULL) == LW_BOUNDARY;
5049  if ( ! contains )
5050  {{
5051  double snaptol;
5052  LWGEOM *snapedge;
5053  LWLINE *snapline;
5054  POINT4D p1, p2;
5055 
5056  LWDEBUGF(1, "Edge %" LWTFMT_ELEMID
5057  " does not contain projected point to it",
5058  edge_id);
5059 
5060  /* In order to reduce the robustness issues, we'll pick
5061  * an edge that contains the projected point, if possible */
5062  if ( i+1 < num )
5063  {
5064  LWDEBUG(1, "But there's another to check");
5065  lwgeom_free(prj);
5066  continue;
5067  }
5068 
5069  /*
5070  -- The tolerance must be big enough for snapping to happen
5071  -- and small enough to snap only to the projected point.
5072  -- Unfortunately ST_Distance returns 0 because it also uses
5073  -- a projected point internally, so we need another way.
5074  */
5075  snaptol = _lwt_minTolerance(prj);
5076  snapedge = _lwt_toposnap(g, prj, snaptol);
5077  snapline = lwgeom_as_lwline(snapedge);
5078 
5079  LWDEBUGF(1, "Edge snapped with tolerance %g", snaptol);
5080 
5081  /* TODO: check if snapping did anything ? */
5082 #if POSTGIS_DEBUG_LEVEL > 0
5083  {
5084  size_t sz;
5085  char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz);
5086  char *wkt2 = lwgeom_to_wkt(snapedge, WKT_EXTENDED, 15, &sz);
5087  LWDEBUGF(1, "Edge %s snapped became %s", wkt1, wkt2);
5088  lwfree(wkt1);
5089  lwfree(wkt2);
5090  }
5091 #endif
5092 
5093 
5094  /*
5095  -- Snapping currently snaps the first point below tolerance
5096  -- so may possibly move first point. See ticket #1631
5097  */
5098  getPoint4d_p(e->geom->points, 0, &p1);
5099  getPoint4d_p(snapline->points, 0, &p2);
5100  LWDEBUGF(1, "Edge first point is %g %g, "
5101  "snapline first point is %g %g",
5102  p1.x, p1.y, p2.x, p2.y);
5103  if ( p1.x != p2.x || p1.y != p2.y )
5104  {
5105  LWDEBUG(1, "Snapping moved first point, re-adding it");
5106  if ( LW_SUCCESS != ptarray_insert_point(snapline->points, &p1, 0) )
5107  {
5108  lwgeom_free(prj);
5109  lwgeom_free(snapedge);
5110  _lwt_release_edges(edges, num);
5111  lwerror("GEOS exception on Contains: %s", lwgeom_geos_errmsg);
5112  return -1;
5113  }
5114 #if POSTGIS_DEBUG_LEVEL > 0
5115  {
5116  size_t sz;
5117  char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz);
5118  LWDEBUGF(1, "Tweaked snapline became %s", wkt1);
5119  lwfree(wkt1);
5120  }
5121 #endif
5122  }
5123 #if POSTGIS_DEBUG_LEVEL > 0
5124  else {
5125  LWDEBUG(1, "Snapping did not move first point");
5126  }
5127 #endif
5128 
5129  if ( -1 == lwt_ChangeEdgeGeom( topo, edge_id, snapline ) )
5130  {
5131  /* TODO: should have invoked lwerror already, leaking memory */
5132  lwgeom_free(prj);
5133  lwgeom_free(snapedge);
5134  _lwt_release_edges(edges, num);
5135  lwerror("lwt_ChangeEdgeGeom failed");
5136  return -1;
5137  }
5138  lwgeom_free(snapedge);
5139  }}
5140 #if POSTGIS_DEBUG_LEVEL > 0
5141  else
5142  {{
5143  size_t sz;
5144  char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz);
5145  char *wkt2 = lwgeom_to_wkt(prj, WKT_EXTENDED, 15, &sz);
5146  LWDEBUGF(1, "Edge %s contains projected point %s", wkt1, wkt2);
5147  lwfree(wkt1);
5148  lwfree(wkt2);
5149  }}
5150 #endif
5151 
5152  /* TODO: pass 1 as last argument (skipChecks) ? */
5153  id = lwt_ModEdgeSplit( topo, edge_id, lwgeom_as_lwpoint(prj), 0 );
5154  if ( -1 == id )
5155  {
5156  /* TODO: should have invoked lwerror already, leaking memory */
5157  lwgeom_free(prj);
5158  _lwt_release_edges(edges, num);
5159  lwerror("lwt_ModEdgeSplit failed");
5160  return -1;
5161  }
5162 
5163  lwgeom_free(prj);
5164 
5165  /*
5166  * TODO: decimate the two new edges with the given tolerance ?
5167  *
5168  * the edge identifiers to decimate would be: edge_id and "id"
5169  * The problem here is that decimation of existing edges
5170  * may introduce intersections or topological inconsistencies,
5171  * for example:
5172  *
5173  * - A node may end up falling on the other side of the edge
5174  * - The decimated edge might intersect another existing edge
5175  *
5176  */
5177 
5178  break; /* we only want to snap a single edge */
5179  }
5180  _lwt_release_edges(edges, num);
5181  }
5182  else
5183  {
5184  /* The point is isolated, add it as such */
5185  /* TODO: pass 1 as last argument (skipChecks) ? */
5186  id = _lwt_AddIsoNode(topo, -1, point, 0, findFace);
5187  if ( moved ) *moved = 0;
5188  if ( -1 == id )
5189  {
5190  /* should have invoked lwerror already, leaking memory */
5191  lwerror("lwt_AddIsoNode failed");
5192  return -1;
5193  }
5194  }
5195 
5196  return id;
5197 }
double x
Definition: liblwgeom.h:354
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:3171
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:790
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1144
static int compare_scored_pointer(const void *si1, const void *si2)
Definition: lwgeom_topo.c:4866
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:161
POINTARRAY * point
Definition: liblwgeom.h:413
#define _LWT_MINTOLERANCE(topo, geom)
Definition: lwgeom_topo.c:4857
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:930
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:1019
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:330
void lwgeom_geos_error(const char *fmt,...)
const LWT_BE_IFACE * be_iface
int ptarray_contains_point_partial(const POINTARRAY *pa, const POINT2D *pt, int check_closed, int *winding_number)
Definition: ptarray.c:733
static double _lwt_minTolerance(LWGEOM *g)
Definition: lwgeom_topo.c:4839
LWT_ELEMID node_id
LWT_ELEMID edge_id
double z
Definition: liblwgeom.h:354
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition: lwgeom.c:170
#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:582
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:2076
#define LW_BOUNDARY
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 initializing 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:335
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:354
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:424
const POINT2D * getPoint2d_cp(const POINTARRAY *pa, uint32_t n)
Returns a POINT2D pointer into the POINTARRAY serialized_ptlist, suitable for reading from...
Definition: lwgeom_api.c:364
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: