From: Bborie Park <bkpark at ucdavis.edu>
Date: Sun, 22 Sep 2013 02:40:44 +0000 (+0000)
Subject: clean up code for ST_Intersects(geometry, raster)
X-Git-Tag: 2.2.0rc1~1366
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=96a070d8483ce4b2325763ffbe9ccfbc386c9331;p=postgis

clean up code for ST_Intersects(geometry, raster)

git-svn-id: http://svn.osgeo.org/postgis/trunk@11966 b70326c6-7e19-0410-871a-916f4a2858ee
---

diff --git a/raster/rt_pg/rtpostgis.sql.in b/raster/rt_pg/rtpostgis.sql.in
index 43dcad6b2..c67b25dc6 100644
--- a/raster/rt_pg/rtpostgis.sql.in
+++ b/raster/rt_pg/rtpostgis.sql.in
@@ -5870,130 +5870,25 @@ CREATE OR REPLACE FUNCTION _st_intersects(geom geometry, rast raster, nband inte
 	RETURNS boolean AS $$
 	DECLARE
 		hasnodata boolean := TRUE;
-		nodata float8 := 0.0;
-		convexhull geometry;
-		geomintersect geometry;
-		x1w double precision := 0.0;
-		x2w double precision := 0.0;
-		y1w double precision := 0.0;
-		y2w double precision := 0.0;
-		x1 integer := 0;
-		x2 integer := 0;
-		x3 integer := 0;
-		x4 integer := 0;
-		y1 integer := 0;
-		y2 integer := 0;
-		y3 integer := 0;
-		y4 integer := 0;
-		x integer := 0;
-		y integer := 0;
-		xinc integer := 0;
-		yinc integer := 0;
-		pixelval double precision;
-		bintersect boolean := FALSE;
-		gtype text;
-		scale float8;
-		w int;
-		h int;
+		_geom geometry;
 	BEGIN
 		IF ST_SRID(rast) != ST_SRID(geom) THEN
 			RAISE EXCEPTION 'Raster and geometry do not have the same SRID';
 		END IF;
 
-		convexhull := ST_ConvexHull(rast);
+		_geom := ST_ConvexHull(rast);
 		IF nband IS NOT NULL THEN
 			SELECT CASE WHEN bmd.nodatavalue IS NULL THEN FALSE ELSE NULL END INTO hasnodata FROM ST_BandMetaData(rast, nband) AS bmd;
 		END IF;
 
-		IF ST_Intersects(geom, convexhull) IS NOT TRUE THEN
+		IF ST_Intersects(geom, _geom) IS NOT TRUE THEN
 			RETURN FALSE;
 		ELSEIF nband IS NULL OR hasnodata IS FALSE THEN
 			RETURN TRUE;
 		END IF;
 
-		-- Get the intersection between with the geometry.
-		-- We will search for withvalue pixel only in this area.
-		geomintersect := st_intersection(geom, convexhull);
-
---RAISE NOTICE 'geomintersect=%', st_astext(geomintersect);
-
-		-- If the intersection is empty, return false
-		IF st_isempty(geomintersect) THEN
-			RETURN FALSE;
-		END IF;
-
-		-- We create a minimalistic buffer around the intersection in order to scan every pixels
-		-- that would touch the edge or intersect with the geometry
-		SELECT sqrt(scalex * scalex + skewy * skewy), width, height INTO scale, w, h FROM ST_Metadata(rast);
-		IF scale != 0 THEN
-			geomintersect := st_buffer(geomintersect, scale / 1000000);
-		END IF;
-
---RAISE NOTICE 'geomintersect2=%', st_astext(geomintersect);
-
-		-- Find the world coordinates of the bounding box of the intersecting area
-		x1w := st_xmin(geomintersect);
-		y1w := st_ymin(geomintersect);
-		x2w := st_xmax(geomintersect);
-		y2w := st_ymax(geomintersect);
-		nodata := st_bandnodatavalue(rast, nband);
-
---RAISE NOTICE 'x1w=%, y1w=%, x2w=%, y2w=%', x1w, y1w, x2w, y2w;
-
-		-- Convert world coordinates to raster coordinates
-		x1 := st_worldtorastercoordx(rast, x1w, y1w);
-		y1 := st_worldtorastercoordy(rast, x1w, y1w);
-		x2 := st_worldtorastercoordx(rast, x2w, y1w);
-		y2 := st_worldtorastercoordy(rast, x2w, y1w);
-		x3 := st_worldtorastercoordx(rast, x1w, y2w);
-		y3 := st_worldtorastercoordy(rast, x1w, y2w);
-		x4 := st_worldtorastercoordx(rast, x2w, y2w);
-		y4 := st_worldtorastercoordy(rast, x2w, y2w);
-
---RAISE NOTICE 'x1=%, y1=%, x2=%, y2=%, x3=%, y3=%, x4=%, y4=%', x1, y1, x2, y2, x3, y3, x4, y4;
-
-		-- Order the raster coordinates for the upcoming FOR loop.
-		x1 := int4smaller(int4smaller(int4smaller(x1, x2), x3), x4);
-		y1 := int4smaller(int4smaller(int4smaller(y1, y2), y3), y4);
-		x2 := int4larger(int4larger(int4larger(x1, x2), x3), x4);
-		y2 := int4larger(int4larger(int4larger(y1, y2), y3), y4);
-
-		-- Make sure the range is not lower than 1.
-		-- This can happen when world coordinate are exactly on the left border
-		-- of the raster and that they do not span on more than one pixel.
-		x1 := int4smaller(int4larger(x1, 1), w);
-		y1 := int4smaller(int4larger(y1, 1), h);
-
-		-- Also make sure the range does not exceed the width and height of the raster.
-		-- This can happen when world coordinate are exactly on the lower right border
-		-- of the raster.
-		x2 := int4smaller(x2, w);
-		y2 := int4smaller(y2, h);
-
---RAISE NOTICE 'x1=%, y1=%, x2=%, y2=%', x1, y1, x2, y2;
-
-		-- Search exhaustively for withvalue pixel on a moving 3x3 grid
-		-- (very often more efficient than searching on a mere 1x1 grid)
-		FOR xinc in 0..2 LOOP
-			FOR yinc in 0..2 LOOP
-				FOR x IN x1+xinc..x2 BY 3 LOOP
-					FOR y IN y1+yinc..y2 BY 3 LOOP
-						-- Check first if the pixel intersects with the geometry. Often many won't.
-						bintersect := NOT st_isempty(st_intersection(st_pixelaspolygon(rast, x, y), geom));
-
-						IF bintersect THEN
-							-- If the pixel really intersects, check its value. Return TRUE if with value.
-							pixelval := st_value(rast, nband, x, y);
-							IF pixelval != nodata THEN
-								RETURN TRUE;
-							END IF;
-						END IF;
-					END LOOP;
-				END LOOP;
-			END LOOP;
-		END LOOP;
-
-		RETURN FALSE;
+		SELECT ST_Collect(t.geom) INTO _geom FROM ST_PixelAsPolygons(rast, nband) AS t;
+		RETURN ST_Intersects(geom, _geom);
 	END;
 	$$ LANGUAGE 'plpgsql' IMMUTABLE
 	COST 1000;