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 43dcad6b220f3267e4a75771ade4ada4559d07d6..c67b25dc6d40c437bcb109d42a9a008a3670f3d6 100644 (file)
--- 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;