Updated TODO and refactored ST_Intersects(geometry, raster) to use

ST_BandSurface()

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

diff --git a/raster/TODO b/raster/TODO
index 1ca9c99198efb17e8dfeecfa56bbc03b37aee164..ea3eba5260e0001607e577d10710e527a0ebdf90 100644 (file)
--- a/raster/TODO
+++ b/raster/TODO
@@ -6,11 +6,6 @@ yet to be ticketed in PostGIS trac and lack specifics.
 
 == Simple Projects ==
 
--- C version of ST_Intersects() in vector-space --
-
-A C version of this function should be faster than the current plpgsql
-version.
-
 == Larger Projects ==
 
 -- Spatial Relationships functions --

diff --git a/raster/rt_pg/rtpostgis.sql.in.c b/raster/rt_pg/rtpostgis.sql.in.c
index a025010c191cdfa10c183b65803af94d9be165c1..4fb31eb594ec84c31df085efe76ee15edb55ccf1 100644 (file)
--- a/raster/rt_pg/rtpostgis.sql.in.c
+++ b/raster/rt_pg/rtpostgis.sql.in.c
@@ -3190,32 +3190,9 @@ CREATE OR REPLACE FUNCTION st_intersects(rast raster, nband integer, geom geomet
 CREATE OR REPLACE FUNCTION _st_intersects(geom geometry, rast raster, nband integer DEFAULT NULL)
        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;
+               hasnodata boolean := TRUE;
+               surface geometry;
        BEGIN
                convexhull := ST_ConvexHull(rast);
                IF nband IS NOT NULL THEN
@@ -3228,88 +3205,13 @@ CREATE OR REPLACE FUNCTION _st_intersects(geom geometry, rast raster, nband inte
                        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;
+               -- get band surface
+               surface := ST_BandSurface(rast, nband);
 
-               -- 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);
+               IF surface IS NOT NULL THEN
+                       RETURN ST_Intersects(geom, surface);
                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_world2rastercoordx(rast, x1w, y1w);
-               y1 := st_world2rastercoordy(rast, x1w, y1w);
-               x2 := st_world2rastercoordx(rast, x2w, y1w);
-               y2 := st_world2rastercoordy(rast, x2w, y1w);
-               x3 := st_world2rastercoordx(rast, x1w, y2w);
-               y3 := st_world2rastercoordy(rast, x1w, y2w);
-               x4 := st_world2rastercoordx(rast, x2w, y2w);
-               y4 := st_world2rastercoordy(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;
        END;
        $$ LANGUAGE 'plpgsql' IMMUTABLE