_lwt_AddPoint()

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

static

Definition at line 4952 of file lwgeom_topo.c.

{
        uint64_t num, i;
        double mindist = FLT_MAX;
        LWT_ISO_NODE *nodes, *nodes2;
        LWT_ISO_EDGE *edges, *edges2;
        LWGEOM *pt = lwpoint_as_lwgeom(point);
        int flds;
        LWT_ELEMID id = 0;
        scored_pointer *sorted;
 
        /* Get tolerance, if 0 was given */
        if (!tol)
                tol = _LWT_MINTOLERANCE(topo, pt);
 
        LWDEBUGG(1, pt, "Adding point");
 
        /*
        -- 1. Check if any existing node is closer than the given precision
        --    and if so pick the closest
        TODO: use WithinBox2D
        */
        flds = LWT_COL_NODE_NODE_ID | LWT_COL_NODE_GEOM;
        nodes = lwt_be_getNodeWithinDistance2D(topo, point, tol, &num, flds, 0);
        if (num == UINT64_MAX)
        {
                lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
                return -1;
        }
  if ( num )
  {
    LWDEBUGF(1, "New point is within %.15g units of %d nodes", tol, num);
    /* Order by distance if there are more than a single return */
    if ( num > 1 )
    {{
      sorted= lwalloc(sizeof(scored_pointer)*num);
      for (i=0; i<num; ++i)
      {
        sorted[i].ptr = nodes+i;
        sorted[i].score = lwgeom_mindistance2d(lwpoint_as_lwgeom(nodes[i].geom), pt);
        LWDEBUGF(1, "Node %" LWTFMT_ELEMID " distance: %.15g",
          ((LWT_ISO_NODE*)(sorted[i].ptr))->node_id, sorted[i].score);
      }
      qsort(sorted, num, sizeof(scored_pointer), compare_scored_pointer);
      nodes2 = lwalloc(sizeof(LWT_ISO_NODE)*num);
      for (i=0; i<num; ++i)
      {
        nodes2[i] = *((LWT_ISO_NODE*)sorted[i].ptr);
      }
      lwfree(sorted);
      lwfree(nodes);
      nodes = nodes2;
    }}
 
    for ( i=0; i<num; ++i )
    {
      LWT_ISO_NODE *n = &(nodes[i]);
      LWGEOM *g = lwpoint_as_lwgeom(n->geom);
      double dist = lwgeom_mindistance2d(g, pt);
      /* TODO: move this check in the previous sort scan ... */
      /* must be closer than tolerated, unless distance is zero */
      if ( dist && dist >= tol ) continue;
      if ( ! id || dist < mindist )
      {
        id = n->node_id;
        mindist = dist;
      }
    }
    if ( id )
    {
      /* found an existing node */
      if ( nodes ) _lwt_release_nodes(nodes, num);
      if ( moved ) *moved = mindist == 0 ? 0 : 1;
      return id;
    }
  }
 
  initGEOS(lwnotice, lwgeom_geos_error);
 
  /*
  -- 2. Check if any existing edge falls within tolerance
  --    and if so split it by a point projected on it
  TODO: use WithinBox2D
  */
  flds = LWT_COL_EDGE_EDGE_ID|LWT_COL_EDGE_GEOM;
  edges = lwt_be_getEdgeWithinDistance2D(topo, point, tol, &num, flds, 0);
  if (num == UINT64_MAX)
  {
    lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
    return -1;
  }
  if ( num )
  {
  LWDEBUGF(1, "New point is within %.15g units of %d edges", tol, num);
 
  /* Order by distance if there are more than a single return */
  if ( num > 1 )
  {{
    int j;
    sorted = lwalloc(sizeof(scored_pointer)*num);
    for (i=0; i<num; ++i)
    {
      sorted[i].ptr = edges+i;
      sorted[i].score = lwgeom_mindistance2d(lwline_as_lwgeom(edges[i].geom), pt);
      LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " distance: %.15g",
        ((LWT_ISO_EDGE*)(sorted[i].ptr))->edge_id, sorted[i].score);
    }
    qsort(sorted, num, sizeof(scored_pointer), compare_scored_pointer);
    edges2 = lwalloc(sizeof(LWT_ISO_EDGE)*num);
    for (j=0, i=0; i<num; ++i)
    {
      if ( sorted[i].score == sorted[0].score )
      {
        edges2[j++] = *((LWT_ISO_EDGE*)sorted[i].ptr);
      }
      else
      {
        lwline_free(((LWT_ISO_EDGE*)sorted[i].ptr)->geom);
      }
    }
    num = j;
    lwfree(sorted);
    lwfree(edges);
    edges = edges2;
  }}
 
  for (i=0; i<num; ++i)
  {
    /* The point is on or near an edge, split the edge */
    LWT_ISO_EDGE *e = &(edges[i]);
    LWGEOM *g = lwline_as_lwgeom(e->geom);
    LWGEOM *prj;
    int contains;
    LWT_ELEMID edge_id = e->edge_id;
 
    LWDEBUGF(1, "Splitting edge %" LWTFMT_ELEMID, edge_id);
 
    /* project point to line, split edge by point */
    prj = lwgeom_closest_point(g, pt);
    if ( moved ) *moved = lwgeom_same(prj,pt) ? 0 : 1;
    if ( lwgeom_has_z(pt) )
    {{
      /*
      -- This is a workaround for ClosestPoint lack of Z support:
      -- http://trac.osgeo.org/postgis/ticket/2033
      */
      LWGEOM *tmp;
      double z;
      POINT4D p4d;
      LWPOINT *prjpt;
      /* add Z to "prj" */
      tmp = lwgeom_force_3dz(prj, 0);
      prjpt = lwgeom_as_lwpoint(tmp);
      getPoint4d_p(point->point, 0, &p4d);
      z = p4d.z;
      getPoint4d_p(prjpt->point, 0, &p4d);
      p4d.z = z;
      ptarray_set_point4d(prjpt->point, 0, &p4d);
      lwgeom_free(prj);
      prj = tmp;
    }}
    const POINT2D *pt = getPoint2d_cp(lwgeom_as_lwpoint(prj)->point, 0);
    contains = ptarray_contains_point_partial(e->geom->points, pt, 0, NULL) == LW_BOUNDARY;
    if ( ! contains )
    {{
      double snaptol;
      LWGEOM *snapedge;
      LWLINE *snapline;
      POINT4D p1, p2;
 
      LWDEBUGF(1, "Edge %" LWTFMT_ELEMID
                  " does not contain projected point to it",
                  edge_id);
 
      /* In order to reduce the robustness issues, we'll pick
       * an edge that contains the projected point, if possible */
      if ( i+1 < num )
      {
        LWDEBUG(1, "But there's another to check");
        lwgeom_free(prj);
        continue;
      }
 
      /*
      -- The tolerance must be big enough for snapping to happen
      -- and small enough to snap only to the projected point.
      -- Unfortunately ST_Distance returns 0 because it also uses
      -- a projected point internally, so we need another way.
      */
      snaptol = _lwt_minTolerance(prj);
      snapedge = _lwt_toposnap(g, prj, snaptol);
      snapline = lwgeom_as_lwline(snapedge);
 
      LWDEBUGF(1, "Edge snapped with tolerance %g", snaptol);
 
      /* TODO: check if snapping did anything ? */
#if POSTGIS_DEBUG_LEVEL > 0
      {
      size_t sz;
      char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz);
      char *wkt2 = lwgeom_to_wkt(snapedge, WKT_EXTENDED, 15, &sz);
      LWDEBUGF(1, "Edge %s snapped became %s", wkt1, wkt2);
      lwfree(wkt1);
      lwfree(wkt2);
      }
#endif
 
 
      /*
      -- Snapping currently snaps the first point below tolerance
      -- so may possibly move first point. See ticket #1631
      */
      getPoint4d_p(e->geom->points, 0, &p1);
      getPoint4d_p(snapline->points, 0, &p2);
      LWDEBUGF(1, "Edge first point is %g %g, "
                  "snapline first point is %g %g",
                  p1.x, p1.y, p2.x, p2.y);
      if ( p1.x != p2.x || p1.y != p2.y )
      {
        LWDEBUG(1, "Snapping moved first point, re-adding it");
        if ( LW_SUCCESS != ptarray_insert_point(snapline->points, &p1, 0) )
        {
          lwgeom_free(prj);
          lwgeom_free(snapedge);
          _lwt_release_edges(edges, num);
          lwerror("GEOS exception on Contains: %s", lwgeom_geos_errmsg);
          return -1;
        }
#if POSTGIS_DEBUG_LEVEL > 0
        {
        size_t sz;
        char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz);
        LWDEBUGF(1, "Tweaked snapline became %s", wkt1);
        lwfree(wkt1);
        }
#endif
      }
#if POSTGIS_DEBUG_LEVEL > 0
      else {
        LWDEBUG(1, "Snapping did not move first point");
      }
#endif
 
      if ( -1 == lwt_ChangeEdgeGeom( topo, edge_id, snapline ) )
      {
        /* TODO: should have invoked lwerror already, leaking memory */
        lwgeom_free(prj);
        lwgeom_free(snapedge);
        _lwt_release_edges(edges, num);
        lwerror("lwt_ChangeEdgeGeom failed");
        return -1;
      }
      lwgeom_free(snapedge);
    }}
#if POSTGIS_DEBUG_LEVEL > 0
    else
    {{
      size_t sz;
      char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz);
      char *wkt2 = lwgeom_to_wkt(prj, WKT_EXTENDED, 15, &sz);
      LWDEBUGF(1, "Edge %s contains projected point %s", wkt1, wkt2);
      lwfree(wkt1);
      lwfree(wkt2);
    }}
#endif
 
    /* TODO: pass 1 as last argument (skipChecks) ? */
    id = lwt_ModEdgeSplit( topo, edge_id, lwgeom_as_lwpoint(prj), 0 );
    if ( -1 == id )
    {
      /* TODO: should have invoked lwerror already, leaking memory */
      lwgeom_free(prj);
      _lwt_release_edges(edges, num);
      lwerror("lwt_ModEdgeSplit failed");
      return -1;
    }
 
    lwgeom_free(prj);
 
    /*
     * TODO: decimate the two new edges with the given tolerance ?
     *
     * the edge identifiers to decimate would be: edge_id and "id"
     * The problem here is that decimation of existing edges
     * may introduce intersections or topological inconsistencies,
     * for example:
     *
     *  - A node may end up falling on the other side of the edge
     *  - The decimated edge might intersect another existing edge
     *
     */
 
    break; /* we only want to snap a single edge */
  }
  _lwt_release_edges(edges, num);
  }
  else
  {
    /* The point is isolated, add it as such */
    /* TODO: pass 1 as last argument (skipChecks) ? */
    id = _lwt_AddIsoNode(topo, -1, point, 0, findFace);
    if ( moved ) *moved = 0;
    if ( -1 == id )
    {
      /* should have invoked lwerror already, leaking memory */
      lwerror("lwt_AddIsoNode failed");
      return -1;
    }
  }
 
  return id;
}

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().

Here is the call graph for this function:

Here is the caller graph for this function: