PostGIS  2.5.0dev-r@@SVN_REVISION@@
geography_measurement_trees.c
Go to the documentation of this file.
1 /**********************************************************************
2  *
3  * PostGIS - Spatial Types for PostgreSQL
4  * http://postgis.net
5  *
6  * PostGIS is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation, either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * PostGIS is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with PostGIS. If not, see <http://www.gnu.org/licenses/>.
18  *
19  **********************************************************************
20  *
21  * ^copyright^
22  *
23  **********************************************************************/
24 
26 
27 
28 /*
29 * Specific tree types include all the generic slots and
30 * their own slots for their trees. We put the implementation
31 * for the CircTreeGeomCache here because we can't shove
32 * the PgSQL specific bits of the code (fcinfo) back into
33 * liblwgeom, where most of the circtree logic lives.
34 */
35 typedef struct {
36  GeomCache gcache;
39 
40 
41 
45 static int
46 CircTreeBuilder(const LWGEOM* lwgeom, GeomCache* cache)
47 {
48  CircTreeGeomCache* circ_cache = (CircTreeGeomCache*)cache;
49  CIRC_NODE* tree = lwgeom_calculate_circ_tree(lwgeom);
50 
51  if ( circ_cache->index )
52  {
53  circ_tree_free(circ_cache->index);
54  circ_cache->index = 0;
55  }
56  if ( ! tree )
57  return LW_FAILURE;
58 
59  circ_cache->index = tree;
60  return LW_SUCCESS;
61 }
62 
63 static int
64 CircTreeFreer(GeomCache* cache)
65 {
66  CircTreeGeomCache* circ_cache = (CircTreeGeomCache*)cache;
67  if ( circ_cache->index )
68  {
69  circ_tree_free(circ_cache->index);
70  circ_cache->index = 0;
71  circ_cache->gcache.argnum = 0;
72  }
73  return LW_SUCCESS;
74 }
75 
76 static GeomCache*
78 {
79  CircTreeGeomCache* cache = palloc(sizeof(CircTreeGeomCache));
80  memset(cache, 0, sizeof(CircTreeGeomCache));
81  return (GeomCache*)cache;
82 }
83 
84 static GeomCacheMethods CircTreeCacheMethods =
85 {
86  CIRC_CACHE_ENTRY,
90 };
91 
92 static CircTreeGeomCache*
93 GetCircTreeGeomCache(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, const GSERIALIZED* g2)
94 {
95  return (CircTreeGeomCache*)GetGeomCache(fcinfo, &CircTreeCacheMethods, g1, g2);
96 }
97 
98 
99 static int
100 CircTreePIP(const CIRC_NODE* tree1, const GSERIALIZED* g1, const POINT4D* in_point)
101 {
102  int tree1_type = gserialized_get_type(g1);
103  GBOX gbox1;
104  GEOGRAPHIC_POINT in_gpoint;
105  POINT3D in_point3d;
106 
107  POSTGIS_DEBUGF(3, "tree1_type=%d", tree1_type);
108 
109  /* If the tree'ed argument is a polygon, do the P-i-P using the tree-based P-i-P */
110  if ( tree1_type == POLYGONTYPE || tree1_type == MULTIPOLYGONTYPE )
111  {
112  POSTGIS_DEBUG(3, "tree is a polygon, using tree PiP");
113  /* Need a gbox to calculate an outside point */
114  if ( LW_FAILURE == gserialized_get_gbox_p(g1, &gbox1) )
115  {
116  LWGEOM* lwgeom1 = lwgeom_from_gserialized(g1);
117  POSTGIS_DEBUG(3, "unable to read gbox from gserialized, calculating from scratch");
118  lwgeom_calculate_gbox_geodetic(lwgeom1, &gbox1);
119  lwgeom_free(lwgeom1);
120  }
121 
122  /* Flip the candidate point into geographics */
123  geographic_point_init(in_point->x, in_point->y, &in_gpoint);
124  geog2cart(&in_gpoint, &in_point3d);
125 
126  /* If the candidate isn't in the tree box, it's not in the tree area */
127  if ( ! gbox_contains_point3d(&gbox1, &in_point3d) )
128  {
129  POSTGIS_DEBUG(3, "in_point3d is not inside the tree gbox, CircTreePIP returning FALSE");
130  return LW_FALSE;
131  }
132  /* The candidate point is in the box, so it *might* be inside the tree */
133  else
134  {
135  POINT2D pt2d_outside; /* latlon */
136  POINT2D pt2d_inside;
137  pt2d_inside.x = in_point->x;
138  pt2d_inside.y = in_point->y;
139  /* Calculate a definitive outside point */
140  gbox_pt_outside(&gbox1, &pt2d_outside);
141  POSTGIS_DEBUGF(3, "p2d_inside=POINT(%g %g) p2d_outside=POINT(%g %g)", pt2d_inside.x, pt2d_inside.y, pt2d_outside.x, pt2d_outside.y);
142  /* Test the candidate point for strict containment */
143  POSTGIS_DEBUG(3, "calling circ_tree_contains_point for PiP test");
144  return circ_tree_contains_point(tree1, &pt2d_inside, &pt2d_outside, NULL);
145  }
146  }
147  else
148  {
149  POSTGIS_DEBUG(3, "tree1 not polygonal, so CircTreePIP returning FALSE");
150  return LW_FALSE;
151  }
152 }
153 
154 
155 static int
156 geography_distance_cache_tolerance(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, const GSERIALIZED* g2, const SPHEROID* s, double tolerance, double* distance)
157 {
158  CircTreeGeomCache* tree_cache = NULL;
159 
160  int type1 = gserialized_get_type(g1);
161  int type2 = gserialized_get_type(g2);
162 
163  Assert(distance);
164 
165  /* Two points? Get outa here... */
166  if ( type1 == POINTTYPE && type2 == POINTTYPE )
167  return LW_FAILURE;
168 
169  /* Fetch/build our cache, if appropriate, etc... */
170  tree_cache = GetCircTreeGeomCache(fcinfo, g1, g2);
171 
172  /* OK, we have an index at the ready! Use it for the one tree argument and */
173  /* fill in the other tree argument */
174  if ( tree_cache && tree_cache->gcache.argnum && tree_cache->index )
175  {
176  CIRC_NODE* circtree_cached = tree_cache->index;
177  CIRC_NODE* circtree = NULL;
178  const GSERIALIZED* g_cached;
179  const GSERIALIZED* g;
180  LWGEOM* lwgeom = NULL;
181  int geomtype_cached;
182  int geomtype;
183  POINT4D p4d;
184 
185  /* We need to dynamically build a tree for the uncached side of the function call */
186  if ( tree_cache->gcache.argnum == 1 )
187  {
188  g_cached = g1;
189  g = g2;
190  geomtype_cached = type1;
191  geomtype = type2;
192  }
193  else if ( tree_cache->gcache.argnum == 2 )
194  {
195  g_cached = g2;
196  g = g1;
197  geomtype_cached = type2;
198  geomtype = type1;
199  }
200  else
201  {
202  lwpgerror("geography_distance_cache this cannot happen!");
203  return LW_FAILURE;
204  }
205 
206  lwgeom = lwgeom_from_gserialized(g);
207  if ( geomtype_cached == POLYGONTYPE || geomtype_cached == MULTIPOLYGONTYPE )
208  {
209  lwgeom_startpoint(lwgeom, &p4d);
210  if ( CircTreePIP(circtree_cached, g_cached, &p4d) )
211  {
212  *distance = 0.0;
213  lwgeom_free(lwgeom);
214  return LW_SUCCESS;
215  }
216  }
217 
218  circtree = lwgeom_calculate_circ_tree(lwgeom);
219  if ( geomtype == POLYGONTYPE || geomtype == MULTIPOLYGONTYPE )
220  {
221  POINT2D p2d;
222  circ_tree_get_point(circtree_cached, &p2d);
223  p4d.x = p2d.x;
224  p4d.y = p2d.y;
225  if ( CircTreePIP(circtree, g, &p4d) )
226  {
227  *distance = 0.0;
228  circ_tree_free(circtree);
229  lwgeom_free(lwgeom);
230  return LW_SUCCESS;
231  }
232  }
233 
234  *distance = circ_tree_distance_tree(circtree_cached, circtree, s, tolerance);
235  circ_tree_free(circtree);
236  lwgeom_free(lwgeom);
237  return LW_SUCCESS;
238  }
239  else
240  {
241  return LW_FAILURE;
242  }
243 }
244 
245 
246 int
247 geography_distance_cache(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, const GSERIALIZED* g2, const SPHEROID* s, double* distance)
248 {
249  return geography_distance_cache_tolerance(fcinfo, g1, g2, s, FP_TOLERANCE, distance);
250 }
251 
252 int
253 geography_dwithin_cache(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, const GSERIALIZED* g2, const SPHEROID* s, double tolerance, int* dwithin)
254 {
255  double distance;
256  /* Ticket #2422, difference between sphere and spheroid distance can trip up the */
257  /* threshold shortcircuit (stopping a calculation before the spheroid distance is actually */
258  /* below the threshold. Lower in the code line, we actually reduce the threshold a little to */
259  /* avoid this. */
260  /* Correct fix: propogate the spheroid information all the way to the bottom of the calculation */
261  /* so the "right thing" can be done in all cases. */
262  if ( LW_SUCCESS == geography_distance_cache_tolerance(fcinfo, g1, g2, s, tolerance, &distance) )
263  {
264  *dwithin = (distance <= (tolerance + FP_TOLERANCE) ? LW_TRUE : LW_FALSE);
265  return LW_SUCCESS;
266  }
267  return LW_FAILURE;
268 }
269 
270 int
271 geography_tree_distance(const GSERIALIZED* g1, const GSERIALIZED* g2, const SPHEROID* s, double tolerance, double* distance)
272 {
273  CIRC_NODE* circ_tree1 = NULL;
274  CIRC_NODE* circ_tree2 = NULL;
275  LWGEOM* lwgeom1 = NULL;
276  LWGEOM* lwgeom2 = NULL;
277  POINT4D pt1, pt2;
278 
279  lwgeom1 = lwgeom_from_gserialized(g1);
280  lwgeom2 = lwgeom_from_gserialized(g2);
281  circ_tree1 = lwgeom_calculate_circ_tree(lwgeom1);
282  circ_tree2 = lwgeom_calculate_circ_tree(lwgeom2);
283  lwgeom_startpoint(lwgeom1, &pt1);
284  lwgeom_startpoint(lwgeom2, &pt2);
285 
286  if ( CircTreePIP(circ_tree1, g1, &pt2) || CircTreePIP(circ_tree2, g2, &pt1) )
287  {
288  *distance = 0.0;
289  }
290  else
291  {
292  /* Calculate tree/tree distance */
293  *distance = circ_tree_distance_tree(circ_tree1, circ_tree2, s, tolerance);
294  }
295 
296  circ_tree_free(circ_tree1);
297  circ_tree_free(circ_tree2);
298  lwgeom_free(lwgeom1);
299  lwgeom_free(lwgeom2);
300  return LW_SUCCESS;
301 }
int gserialized_get_gbox_p(const GSERIALIZED *g, GBOX *box)
Read the bounding box off a serialization and calculate one if it is not already there.
Definition: g_serialized.c:639
double x
Definition: liblwgeom.h:354
uint32_t gserialized_get_type(const GSERIALIZED *s)
Extract the geometry type from the serialized form (it hides in the anonymous data area...
Definition: g_serialized.c:86
int lwgeom_calculate_gbox_geodetic(const LWGEOM *geom, GBOX *gbox)
Calculate the geodetic bounding box for an LWGEOM.
Definition: lwgeodetic.c:2940
static int CircTreeFreer(GeomCache *cache)
Note that p1 and p2 are pointers into an independent POINTARRAY, do not free them.
LWGEOM * lwgeom_from_gserialized(const GSERIALIZED *g)
Allocate a new LWGEOM from a GSERIALIZED.
int circ_tree_get_point(const CIRC_NODE *node, POINT2D *pt)
Returns a POINT2D that is a vertex of the input shape.
int lwgeom_startpoint(const LWGEOM *lwgeom, POINT4D *pt)
Definition: lwgeom.c:2115
#define POLYGONTYPE
Definition: liblwgeom.h:86
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1144
static int geography_distance_cache_tolerance(FunctionCallInfoData *fcinfo, const GSERIALIZED *g1, const GSERIALIZED *g2, const SPHEROID *s, double tolerance, double *distance)
int gbox_contains_point3d(const GBOX *gbox, const POINT3D *pt)
Return true if the point is inside the gbox.
Definition: g_box.c:254
#define LW_SUCCESS
Definition: liblwgeom.h:79
static int CircTreeBuilder(const LWGEOM *lwgeom, GeomCache *cache)
Builder, freeer and public accessor for cached CIRC_NODE trees.
static GeomCache * CircTreeAllocator(void)
double circ_tree_distance_tree(const CIRC_NODE *n1, const CIRC_NODE *n2, const SPHEROID *spheroid, double threshold)
int geography_tree_distance(const GSERIALIZED *g1, const GSERIALIZED *g2, const SPHEROID *s, double tolerance, double *distance)
Point in spherical coordinates on the world.
Definition: lwgeodetic.h:47
CIRC_NODE * lwgeom_calculate_circ_tree(const LWGEOM *lwgeom)
static CircTreeGeomCache * GetCircTreeGeomCache(FunctionCallInfoData *fcinfo, const GSERIALIZED *g1, const GSERIALIZED *g2)
#define LW_FAILURE
Definition: liblwgeom.h:78
static int CircTreePIP(const CIRC_NODE *tree1, const GSERIALIZED *g1, const POINT4D *in_point)
double x
Definition: liblwgeom.h:330
#define LW_FALSE
Definition: liblwgeom.h:76
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
int circ_tree_contains_point(const CIRC_NODE *node, const POINT2D *pt, const POINT2D *pt_outside, int *on_boundary)
Walk the tree and count intersections between the stab line and the edges.
#define FP_TOLERANCE
Floating point comparators.
char * s
Definition: cu_in_wkt.c:23
double y
Definition: liblwgeom.h:330
Datum distance(PG_FUNCTION_ARGS)
void geog2cart(const GEOGRAPHIC_POINT *g, POINT3D *p)
Convert spherical coordinates to cartesian coordinates on unit sphere.
Definition: lwgeodetic.c:373
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:89
void circ_tree_free(CIRC_NODE *node)
Recurse from top of node tree and free all children.
void gbox_pt_outside(const GBOX *gbox, POINT2D *pt_outside)
Calculate a spherical point that falls outside the geocentric gbox.
Definition: lwgeodetic.c:1465
int geography_dwithin_cache(FunctionCallInfoData *fcinfo, const GSERIALIZED *g1, const GSERIALIZED *g2, const SPHEROID *s, double tolerance, int *dwithin)
void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g)
Initialize a geographic point.
Definition: lwgeodetic.c:171
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:84
int geography_distance_cache(FunctionCallInfoData *fcinfo, const GSERIALIZED *g1, const GSERIALIZED *g2, const SPHEROID *s, double *distance)
static GeomCacheMethods CircTreeCacheMethods
double y
Definition: liblwgeom.h:354