From f350207f28e473c84b7b9824725988279776adff Mon Sep 17 00:00:00 2001 From: Bborie Park Date: Fri, 13 Jul 2012 23:44:29 +0000 Subject: [PATCH] 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 --- raster/TODO | 5 -- raster/rt_pg/rtpostgis.sql.in.c | 110 ++------------------------------ 2 files changed, 6 insertions(+), 109 deletions(-) diff --git a/raster/TODO b/raster/TODO index 1ca9c9919..ea3eba526 100644 --- 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 a025010c1..4fb31eb59 100644 --- 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 -- 2.40.0