PostGIS  2.2.7dev-r@@SVN_REVISION@@
geography_measurement_trees.c
Go to the documentation of this file.
2 
3 
4 /*
5 * Specific tree types include all the generic slots and
6 * their own slots for their trees. We put the implementation
7 * for the CircTreeGeomCache here because we can't shove
8 * the PgSQL specific bits of the code (fcinfo) back into
9 * liblwgeom, where most of the circtree logic lives.
10 */
11 typedef struct {
12  int type; // <GeomCache>
15  size_t geom1_size; //
16  size_t geom2_size; //
17  int32 argnum; // </GeomCache>
20 
21 
22 
26 static int
27 CircTreeBuilder(const LWGEOM* lwgeom, GeomCache* cache)
28 {
29  CircTreeGeomCache* circ_cache = (CircTreeGeomCache*)cache;
30  CIRC_NODE* tree = lwgeom_calculate_circ_tree(lwgeom);
31 
32  if ( circ_cache->index )
33  {
34  circ_tree_free(circ_cache->index);
35  circ_cache->index = 0;
36  }
37  if ( ! tree )
38  return LW_FAILURE;
39 
40  circ_cache->index = tree;
41  return LW_SUCCESS;
42 }
43 
44 static int
45 CircTreeFreer(GeomCache* cache)
46 {
47  CircTreeGeomCache* circ_cache = (CircTreeGeomCache*)cache;
48  if ( circ_cache->index )
49  {
50  circ_tree_free(circ_cache->index);
51  circ_cache->index = 0;
52  circ_cache->argnum = 0;
53  }
54  return LW_SUCCESS;
55 }
56 
57 static GeomCache*
59 {
60  CircTreeGeomCache* cache = palloc(sizeof(CircTreeGeomCache));
61  memset(cache, 0, sizeof(CircTreeGeomCache));
62  return (GeomCache*)cache;
63 }
64 
65 static GeomCacheMethods CircTreeCacheMethods =
66 {
67  CIRC_CACHE_ENTRY,
71 };
72 
73 static CircTreeGeomCache*
74 GetCircTreeGeomCache(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, const GSERIALIZED* g2)
75 {
76  return (CircTreeGeomCache*)GetGeomCache(fcinfo, &CircTreeCacheMethods, g1, g2);
77 }
78 
79 
80 static int
81 CircTreePIP(const CIRC_NODE* tree1, const GSERIALIZED* g1, const POINT4D* in_point)
82 {
83  int tree1_type = gserialized_get_type(g1);
84  GBOX gbox1;
85  GEOGRAPHIC_POINT in_gpoint;
86  POINT3D in_point3d;
87 
88  POSTGIS_DEBUGF(3, "tree1_type=%d", tree1_type);
89 
90  /* If the tree'ed argument is a polygon, do the P-i-P using the tree-based P-i-P */
91  if ( tree1_type == POLYGONTYPE || tree1_type == MULTIPOLYGONTYPE )
92  {
93  POSTGIS_DEBUG(3, "tree is a polygon, using tree PiP");
94  /* Need a gbox to calculate an outside point */
95  if ( LW_FAILURE == gserialized_get_gbox_p(g1, &gbox1) )
96  {
97  LWGEOM* lwgeom1 = lwgeom_from_gserialized(g1);
98  POSTGIS_DEBUG(3, "unable to read gbox from gserialized, calculating from scratch");
99  lwgeom_calculate_gbox_geodetic(lwgeom1, &gbox1);
100  lwgeom_free(lwgeom1);
101  }
102 
103  /* Flip the candidate point into geographics */
104  geographic_point_init(in_point->x, in_point->y, &in_gpoint);
105  geog2cart(&in_gpoint, &in_point3d);
106 
107  /* If the candidate isn't in the tree box, it's not in the tree area */
108  if ( ! gbox_contains_point3d(&gbox1, &in_point3d) )
109  {
110  POSTGIS_DEBUG(3, "in_point3d is not inside the tree gbox, CircTreePIP returning FALSE");
111  return LW_FALSE;
112  }
113  /* The candidate point is in the box, so it *might* be inside the tree */
114  else
115  {
116  POINT2D pt2d_outside; /* latlon */
117  POINT2D pt2d_inside;
118  pt2d_inside.x = in_point->x;
119  pt2d_inside.y = in_point->y;
120  /* Calculate a definitive outside point */
121  gbox_pt_outside(&gbox1, &pt2d_outside);
122  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);
123  /* Test the candidate point for strict containment */
124  POSTGIS_DEBUG(3, "calling circ_tree_contains_point for PiP test");
125  return circ_tree_contains_point(tree1, &pt2d_inside, &pt2d_outside, NULL);
126  }
127  }
128  else
129  {
130  POSTGIS_DEBUG(3, "tree1 not polygonal, so CircTreePIP returning FALSE");
131  return LW_FALSE;
132  }
133 }
134 
135 
136 static int
137 geography_distance_cache_tolerance(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, const GSERIALIZED* g2, const SPHEROID* s, double tolerance, double* distance)
138 {
139  CircTreeGeomCache* tree_cache = NULL;
140 
141  int type1 = gserialized_get_type(g1);
142  int type2 = gserialized_get_type(g2);
143 
144  Assert(distance);
145 
146  /* Two points? Get outa here... */
147  if ( type1 == POINTTYPE && type2 == POINTTYPE )
148  return LW_FAILURE;
149 
150  /* Fetch/build our cache, if appropriate, etc... */
151  tree_cache = GetCircTreeGeomCache(fcinfo, g1, g2);
152 
153  /* OK, we have an index at the ready! Use it for the one tree argument and */
154  /* fill in the other tree argument */
155  if ( tree_cache && tree_cache->argnum && tree_cache->index )
156  {
157  CIRC_NODE* circtree_cached = tree_cache->index;
158  CIRC_NODE* circtree = NULL;
159  const GSERIALIZED* g_cached;
160  const GSERIALIZED* g;
161  LWGEOM* lwgeom = NULL;
162  int geomtype_cached;
163  int geomtype;
164  POINT4D p4d;
165 
166  /* We need to dynamically build a tree for the uncached side of the function call */
167  if ( tree_cache->argnum == 1 )
168  {
169  g_cached = g1;
170  g = g2;
171  geomtype_cached = type1;
172  geomtype = type2;
173  }
174  else if ( tree_cache->argnum == 2 )
175  {
176  g_cached = g2;
177  g = g1;
178  geomtype_cached = type2;
179  geomtype = type1;
180  }
181  else
182  {
183  lwpgerror("geography_distance_cache this cannot happen!");
184  return LW_FAILURE;
185  }
186 
187  lwgeom = lwgeom_from_gserialized(g);
188  if ( geomtype_cached == POLYGONTYPE || geomtype_cached == MULTIPOLYGONTYPE )
189  {
190  lwgeom_startpoint(lwgeom, &p4d);
191  if ( CircTreePIP(circtree_cached, g_cached, &p4d) )
192  {
193  *distance = 0.0;
194  lwgeom_free(lwgeom);
195  return LW_SUCCESS;
196  }
197  }
198 
199  circtree = lwgeom_calculate_circ_tree(lwgeom);
200  if ( geomtype == POLYGONTYPE || geomtype == MULTIPOLYGONTYPE )
201  {
202  POINT2D p2d;
203  circ_tree_get_point(circtree_cached, &p2d);
204  p4d.x = p2d.x;
205  p4d.y = p2d.y;
206  if ( CircTreePIP(circtree, g, &p4d) )
207  {
208  *distance = 0.0;
209  circ_tree_free(circtree);
210  lwgeom_free(lwgeom);
211  return LW_SUCCESS;
212  }
213  }
214 
215  *distance = circ_tree_distance_tree(circtree_cached, circtree, s, tolerance);
216  circ_tree_free(circtree);
217  lwgeom_free(lwgeom);
218  return LW_SUCCESS;
219  }
220  else
221  {
222  return LW_FAILURE;
223  }
224 }
225 
226 
227 int
228 geography_distance_cache(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, const GSERIALIZED* g2, const SPHEROID* s, double* distance)
229 {
230  return geography_distance_cache_tolerance(fcinfo, g1, g2, s, FP_TOLERANCE, distance);
231 }
232 
233 int
234 geography_dwithin_cache(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, const GSERIALIZED* g2, const SPHEROID* s, double tolerance, int* dwithin)
235 {
236  double distance;
237  /* Ticket #2422, difference between sphere and spheroid distance can trip up the */
238  /* threshold shortcircuit (stopping a calculation before the spheroid distance is actually */
239  /* below the threshold. Lower in the code line, we actually reduce the threshold a little to */
240  /* avoid this. */
241  /* Correct fix: propogate the spheroid information all the way to the bottom of the calculation */
242  /* so the "right thing" can be done in all cases. */
243  if ( LW_SUCCESS == geography_distance_cache_tolerance(fcinfo, g1, g2, s, tolerance, &distance) )
244  {
245  *dwithin = (distance <= (tolerance + FP_TOLERANCE) ? LW_TRUE : LW_FALSE);
246  return LW_SUCCESS;
247  }
248  return LW_FAILURE;
249 }
250 
251 int
252 geography_tree_distance(const GSERIALIZED* g1, const GSERIALIZED* g2, const SPHEROID* s, double tolerance, double* distance)
253 {
254  CIRC_NODE* circ_tree1 = NULL;
255  CIRC_NODE* circ_tree2 = NULL;
256  LWGEOM* lwgeom1 = NULL;
257  LWGEOM* lwgeom2 = NULL;
258  POINT4D pt1, pt2;
259 
260  lwgeom1 = lwgeom_from_gserialized(g1);
261  lwgeom2 = lwgeom_from_gserialized(g2);
262  circ_tree1 = lwgeom_calculate_circ_tree(lwgeom1);
263  circ_tree2 = lwgeom_calculate_circ_tree(lwgeom2);
264  lwgeom_startpoint(lwgeom1, &pt1);
265  lwgeom_startpoint(lwgeom2, &pt2);
266 
267  if ( CircTreePIP(circ_tree1, g1, &pt2) || CircTreePIP(circ_tree2, g2, &pt1) )
268  {
269  *distance = 0.0;
270  }
271  else
272  {
273  /* Calculate tree/tree distance */
274  *distance = circ_tree_distance_tree(circ_tree1, circ_tree2, s, tolerance);
275  }
276 
277  circ_tree_free(circ_tree1);
278  circ_tree_free(circ_tree2);
279  lwgeom_free(lwgeom1);
280  lwgeom_free(lwgeom2);
281  return LW_SUCCESS;
282 }
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:371
double x
Definition: liblwgeom.h:336
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:55
int lwgeom_calculate_gbox_geodetic(const LWGEOM *geom, GBOX *gbox)
Calculate the geodetic bounding box for an LWGEOM.
Definition: lwgeodetic.c:2615
static int CircTreeFreer(GeomCache *cache)
unsigned int int32
Definition: shpopen.c:273
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:1837
#define POLYGONTYPE
Definition: liblwgeom.h:72
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1050
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:224
#define LW_SUCCESS
Definition: liblwgeom.h:65
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:32
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:64
static int CircTreePIP(const CIRC_NODE *tree1, const GSERIALIZED *g1, const POINT4D *in_point)
double x
Definition: liblwgeom.h:312
#define LW_FALSE
Definition: liblwgeom.h:62
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:61
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:312
Datum distance(PG_FUNCTION_ARGS)
void geog2cart(const GEOGRAPHIC_POINT *g, POINT3D *p)
Convert spherical coordinates to cartesion coordinates on unit sphere.
Definition: lwgeodetic.c:354
#define MULTIPOLYGONTYPE
Definition: liblwgeom.h:75
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:1444
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:156
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:70
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:336