PostGIS  3.0.0dev-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 4885 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().

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