geography_measurement_trees.c 8.79 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
/**********************************************************************
 *
 * PostGIS - Spatial Types for PostgreSQL
 * http://postgis.net
 *
 * PostGIS is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 2 of the License, or
 * (at your option) any later version.
 *
 * PostGIS is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with PostGIS.  If not, see <http://www.gnu.org/licenses/>.
 *
 **********************************************************************
 *
 * ^copyright^
 *
 **********************************************************************/

25 26 27 28
#include "geography_measurement_trees.h"


/*
29
* Specific tree types include all the generic slots and
30 31 32 33 34 35
* their own slots for their trees. We put the implementation
* for the CircTreeGeomCache here because we can't shove
* the PgSQL specific bits of the code (fcinfo) back into
* liblwgeom, where most of the circtree logic lives.
*/
typedef struct {
36
	int                     type;       // <GeomCache>
37 38 39 40
	GSERIALIZED*                geom1;      //
	GSERIALIZED*                geom2;      //
	size_t                      geom1_size; //
	size_t                      geom2_size; //
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
	int32                       argnum;     // </GeomCache>
	CIRC_NODE*                  index;
} CircTreeGeomCache;



/**
* Builder, freeer and public accessor for cached CIRC_NODE trees
*/
static int
CircTreeBuilder(const LWGEOM* lwgeom, GeomCache* cache)
{
	CircTreeGeomCache* circ_cache = (CircTreeGeomCache*)cache;
	CIRC_NODE* tree = lwgeom_calculate_circ_tree(lwgeom);

	if ( circ_cache->index )
	{
		circ_tree_free(circ_cache->index);
		circ_cache->index = 0;
	}
	if ( ! tree )
		return LW_FAILURE;
	
	circ_cache->index = tree;
	return LW_SUCCESS;
}

static int
CircTreeFreer(GeomCache* cache)
{
	CircTreeGeomCache* circ_cache = (CircTreeGeomCache*)cache;
72
	if ( circ_cache->index )
73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102
	{
		circ_tree_free(circ_cache->index);
		circ_cache->index = 0;
		circ_cache->argnum = 0;
	}
	return LW_SUCCESS;
}

static GeomCache*
CircTreeAllocator(void)
{
	CircTreeGeomCache* cache = palloc(sizeof(CircTreeGeomCache));
	memset(cache, 0, sizeof(CircTreeGeomCache));
	return (GeomCache*)cache;
}

static GeomCacheMethods CircTreeCacheMethods =
{
	CIRC_CACHE_ENTRY,
	CircTreeBuilder,
	CircTreeFreer,
	CircTreeAllocator
};

static CircTreeGeomCache*
GetCircTreeGeomCache(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, const GSERIALIZED* g2)
{
	return (CircTreeGeomCache*)GetGeomCache(fcinfo, &CircTreeCacheMethods, g1, g2);
}

103

104
static int
105
CircTreePIP(const CIRC_NODE* tree1, const GSERIALIZED* g1, const POINT4D* in_point)
106 107 108 109 110 111
{
	int tree1_type = gserialized_get_type(g1);
	GBOX gbox1;
	GEOGRAPHIC_POINT in_gpoint;
	POINT3D in_point3d;

112
	POSTGIS_DEBUGF(3, "tree1_type=%d", tree1_type);
113 114 115 116 117 118 119 120 121 122 123 124 125

	/* If the tree'ed argument is a polygon, do the P-i-P using the tree-based P-i-P */
	if ( tree1_type == POLYGONTYPE || tree1_type == MULTIPOLYGONTYPE )
	{
		POSTGIS_DEBUG(3, "tree is a polygon, using tree PiP");
		/* Need a gbox to calculate an outside point */
		if ( LW_FAILURE == gserialized_get_gbox_p(g1, &gbox1) )
		{
			LWGEOM* lwgeom1 = lwgeom_from_gserialized(g1);
			POSTGIS_DEBUG(3, "unable to read gbox from gserialized, calculating from scratch");
			lwgeom_calculate_gbox_geodetic(lwgeom1, &gbox1);
			lwgeom_free(lwgeom1);
		}
126
			
127
		/* Flip the candidate point into geographics */
128
		geographic_point_init(in_point->x, in_point->y, &in_gpoint);
129 130 131 132 133 134 135 136 137 138 139 140 141
		geog2cart(&in_gpoint, &in_point3d);
		
		/* If the candidate isn't in the tree box, it's not in the tree area */
		if ( ! gbox_contains_point3d(&gbox1, &in_point3d) )
		{
			POSTGIS_DEBUG(3, "in_point3d is not inside the tree gbox, CircTreePIP returning FALSE");
			return LW_FALSE;
		}
		/* The candidate point is in the box, so it *might* be inside the tree */
		else
		{
			POINT2D pt2d_outside; /* latlon */
			POINT2D pt2d_inside;
142
			pt2d_inside.x = in_point->x;
143
			pt2d_inside.y = in_point->y;
144 145 146 147 148 149 150 151 152 153
			/* Calculate a definitive outside point */
			gbox_pt_outside(&gbox1, &pt2d_outside);
			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);
			/* Test the candidate point for strict containment */
			POSTGIS_DEBUG(3, "calling circ_tree_contains_point for PiP test");
			return circ_tree_contains_point(tree1, &pt2d_inside, &pt2d_outside, NULL);
		}
	}
	else
	{
154
		POSTGIS_DEBUG(3, "tree1 not polygonal, so CircTreePIP returning FALSE");
155 156 157 158
		return LW_FALSE;
	}		
}

159 160 161

static int
geography_distance_cache_tolerance(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, const GSERIALIZED* g2, const SPHEROID* s, double tolerance, double* distance)
162 163
{
	CircTreeGeomCache* tree_cache = NULL;
164 165 166

	int type1 = gserialized_get_type(g1);
	int type2 = gserialized_get_type(g2);
167 168 169 170
	
	Assert(distance);
	
	/* Two points? Get outa here... */
171
	if ( type1 == POINTTYPE && type2 == POINTTYPE )
172 173 174 175 176
		return LW_FAILURE;

	/* Fetch/build our cache, if appropriate, etc... */
	tree_cache = GetCircTreeGeomCache(fcinfo, g1, g2);
	
177 178 179
	/* OK, we have an index at the ready! Use it for the one tree argument and */
	/* fill in the other tree argument */
	if ( tree_cache && tree_cache->argnum && tree_cache->index )
180
	{
181 182 183 184
		CIRC_NODE* circtree_cached = tree_cache->index;
		CIRC_NODE* circtree = NULL;
		const GSERIALIZED* g_cached;
		const GSERIALIZED* g;
185
		LWGEOM* lwgeom = NULL;
186 187 188 189
		int geomtype_cached;
		int geomtype;
		POINT4D p4d;
		
190 191 192
		/* We need to dynamically build a tree for the uncached side of the function call */
		if ( tree_cache->argnum == 1 )
		{
193 194 195 196
			g_cached = g1;
			g = g2;
			geomtype_cached = type1;
			geomtype = type2;
197 198 199
		}
		else if ( tree_cache->argnum == 2 )
		{
200 201 202 203
			g_cached = g2;
			g = g1;
			geomtype_cached = type2;
			geomtype = type1;
204 205 206
		}
		else
		{
207
			lwpgerror("geography_distance_cache this cannot happen!");
208
			return LW_FAILURE;
209 210
		}
		
211 212 213 214 215 216 217 218 219 220 221 222 223
		lwgeom = lwgeom_from_gserialized(g);
		if ( geomtype_cached == POLYGONTYPE || geomtype_cached == MULTIPOLYGONTYPE )
		{
			lwgeom_startpoint(lwgeom, &p4d);
			if ( CircTreePIP(circtree_cached, g_cached, &p4d) )
			{
				*distance = 0.0;
				lwgeom_free(lwgeom);
				return LW_SUCCESS;
			}
		}
		
		circtree = lwgeom_calculate_circ_tree(lwgeom);
224
		if ( geomtype == POLYGONTYPE || geomtype == MULTIPOLYGONTYPE )
225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240
		{
			POINT2D p2d;
			circ_tree_get_point(circtree_cached, &p2d);
			p4d.x = p2d.x;
			p4d.y = p2d.y;
			if ( CircTreePIP(circtree, g, &p4d) )
			{
				*distance = 0.0;
				circ_tree_free(circtree);
				lwgeom_free(lwgeom);
				return LW_SUCCESS;
			}
		}

		*distance = circ_tree_distance_tree(circtree_cached, circtree, s, tolerance);
		circ_tree_free(circtree);
241 242 243 244 245 246 247 248 249
		lwgeom_free(lwgeom);	
		return LW_SUCCESS;
	}
	else
	{
		return LW_FAILURE;
	}
}

250 251 252 253 254 255 256

int
geography_distance_cache(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, const GSERIALIZED* g2, const SPHEROID* s, double* distance)
{
	return geography_distance_cache_tolerance(fcinfo, g1, g2, s, FP_TOLERANCE, distance);
}

257 258 259 260
int
geography_dwithin_cache(FunctionCallInfoData* fcinfo, const GSERIALIZED* g1, const GSERIALIZED* g2, const SPHEROID* s, double tolerance, int* dwithin)
{
	double distance;
261 262 263 264 265 266 267
	/* Ticket #2422, difference between sphere and spheroid distance can trip up the */
	/* threshold shortcircuit (stopping a calculation before the spheroid distance is actually */
	/* below the threshold. Lower in the code line, we actually reduce the threshold a little to */
	/* avoid this. */
	/* Correct fix: propogate the spheroid information all the way to the bottom of the calculation */
	/* so the "right thing" can be done in all cases. */
	if ( LW_SUCCESS == geography_distance_cache_tolerance(fcinfo, g1, g2, s, tolerance, &distance) )
268
	{
269
		*dwithin = (distance <= (tolerance + FP_TOLERANCE) ? LW_TRUE : LW_FALSE);
270 271
		return LW_SUCCESS;
	}
272
	return LW_FAILURE;
273
}
274
	
275 276 277 278 279 280 281
int
geography_tree_distance(const GSERIALIZED* g1, const GSERIALIZED* g2, const SPHEROID* s, double tolerance, double* distance)
{
	CIRC_NODE* circ_tree1 = NULL;
	CIRC_NODE* circ_tree2 = NULL;
	LWGEOM* lwgeom1 = NULL;
	LWGEOM* lwgeom2 = NULL;
282
	POINT4D pt1, pt2;
283 284 285 286 287
	
	lwgeom1 = lwgeom_from_gserialized(g1);
	lwgeom2 = lwgeom_from_gserialized(g2);
	circ_tree1 = lwgeom_calculate_circ_tree(lwgeom1);
	circ_tree2 = lwgeom_calculate_circ_tree(lwgeom2);
288 289
	lwgeom_startpoint(lwgeom1, &pt1);
	lwgeom_startpoint(lwgeom2, &pt2);
290
	
291
	if ( CircTreePIP(circ_tree1, g1, &pt2) || CircTreePIP(circ_tree2, g2, &pt1) )
292 293 294
	{
		*distance = 0.0;
	}
295
	else
296 297 298 299 300 301 302 303 304 305 306
	{
		/* Calculate tree/tree distance */
		*distance = circ_tree_distance_tree(circ_tree1, circ_tree2, s, tolerance);
	}
	
	circ_tree_free(circ_tree1);
	circ_tree_free(circ_tree2);
	lwgeom_free(lwgeom1);
	lwgeom_free(lwgeom2);
	return LW_SUCCESS;
}