1---------------------------------------------------------------------- 2-- 3-- 4-- Copyright (c) 2009-2010 Pierre Racine <pierre.racine@sbf.ulaval.ca> 5-- 6---------------------------------------------------------------------- 7 8-- Note: The functions found in this file are for exclusive usage of ST_MapAlgebra2 9 10CREATE OR REPLACE FUNCTION max(a int, b int) 11 RETURNS int 12 AS 'SELECT CASE WHEN $1 < $2 THEN $2 ELSE $1 END' 13 LANGUAGE 'SQL' IMMUTABLE STRICT; 14 15CREATE OR REPLACE FUNCTION min(a int, b int) 16 RETURNS int 17 AS 'SELECT CASE WHEN $1 < $2 THEN $1 ELSE $2 END' 18 LANGUAGE 'SQL' IMMUTABLE STRICT; 19 20--DROP FUNCTION _MapAlgebraParts(r1x int, r1y int, r1w int, r1h int, r2x int, r2y int, r2w int, r2h int); 21CREATE OR REPLACE FUNCTION _MapAlgebraParts(r1x int, 22 r1y int, 23 r1w int, 24 r1h int, 25 r2x int, 26 r2y int, 27 r2w int, 28 r2h int) 29 RETURNS int[] 30 AS $$ 31 DECLARE 32 z11x int; 33 z11y int; 34 z11w int; 35 z11h int; 36 z12x int; 37 z12y int; 38 z12w int; 39 z12h int; 40 z13x int; 41 z13y int; 42 z13w int; 43 z13h int; 44 z14x int; 45 z14y int; 46 z14w int; 47 z14h int; 48 z21x int; 49 z21y int; 50 z21w int; 51 z21h int; 52 z22x int; 53 z22y int; 54 z22w int; 55 z22h int; 56 z23x int; 57 z23y int; 58 z23w int; 59 z23h int; 60 z24x int; 61 z24y int; 62 z24w int; 63 z24h int; 64 zcx int; 65 zcy int; 66 zcw int; 67 zch int; 68 BEGIN 69 z11x := r1x; 70 z11y := r1y; 71 z11w := r1w; 72 z11h := min(max(0, r2y - r1y), r1h); 73 74 z12x := r1x; 75 z12y := z11y + z11h; 76 z12w := min(max(0, r2x - r1x), r1w); 77 z12h := max(0, min(max(0, r2y + r2h - (r1y + z11h)), z11y + r1h - z12y)); 78 79 z13x := max(min(r1x + r1w, r2x + r2w), r1x); 80 z13y := z12y; 81 z13w := min(max(0, r1x + r1w - (r2x + r2w)), r1w); 82 z13h := z12h; 83 84 z14x := r1x; 85 z14y := z12y + z12h; 86 z14w := r1w; 87 z14h := min(max(0, r1y + r1h - (r2y + r2h)), r1h); 88 89 z21x := r2x; 90 z21y := r2y; 91 z21w := r2w; 92 z21h := min(max(0, r1y - r2y), r2h); 93 94 z22x := r2x; 95 z22y := z21y + z21h; 96 z22w := min(max(0, r1x - r2x), r2w); 97 z22h := max(0, min(max(0, r1y + r1h - (r2y + z21h)), z21y + r2h - z22y)); 98 99 z23x := max(min(r2x + r2w, r1x + r1w), r2x); 100 z23y := z22y; 101 z23w := min(max(0, r2x + r2w - (r1x + r1w)), r2w); 102 z23h := z22h; 103 104 z24x := r2x; 105 z24y := z22y + z22h; 106 z24w := r2w; 107 z24h := min(max(0, r2y + r2h - (r1y + r1h)), r2h); 108 109 zcx := z12x + z12w; 110 zcy := z12y; 111 zcw := z13x - (z12x + z12w); 112 zch := z14y - z12y; 113 114 -- Merge z11 with z12 if they are continuous parts of the same vertical bar 115 IF z12h > 0 AND z12x = z11x AND z12w = z11w THEN 116 z12y := z11y; 117 z12h := z11h + z12h; 118 z11h := 0; 119 END IF; 120 121 -- Merge z11 with z13 if they are continuous parts of the same vertical bar 122 IF z13h > 0 AND z13x = z11x AND z13w = z11w THEN 123 z13y := z11y; 124 z13h := z11h + z13h; 125 z11h := 0; 126 END IF; 127 128 -- Merge z12 with z14 if they are continuous parts of the same vertical bar 129 IF z14h > 0 AND z14x = z12x AND z14w = z12w THEN 130 z14y := z12y; 131 z14h := z12h + z14h; 132 z12h := 0; 133 END IF; 134 135 -- Merge z13 with z14 if they are continuous parts of the same vertical bar 136 IF z14h > 0 AND z14x = z13x AND z14w = z13w THEN 137 z14y := z13y; 138 z14h := z13h + z14h; 139 z13h := 0; 140 END IF; 141 142 -- Merge z21 with z22 if they are continuous parts of the same vertical bar 143 IF z22h > 0 AND z22x = z21x AND z22w = z21w THEN 144 z22y := z21y; 145 z22h := z21h + z22h; 146 z21h := 0; 147 END IF; 148 149 -- Merge z21 with z23 if they are continuous parts of the same vertical bar 150 IF z23h > 0 AND z23x = z21x AND z23w = z21w THEN 151 z23y := z21y; 152 z23h := z21h + z23h; 153 z21h := 0; 154 END IF; 155 156 -- Merge z22 with z24 if they are continuous parts of the same vertical bar 157 IF z24h > 0 AND z24x = z22x AND z24w = z22w THEN 158 z24y := z22y; 159 z24h := z22h + z24h; 160 z22h := 0; 161 END IF; 162 163 -- Merge z23 with z24 if they are continuous parts of the same vertical bar 164 IF z24h > 0 AND z24x = z23x AND z24w = z23w THEN 165 z24y := z23y; 166 z24h := z23h + z24h; 167 z23h := 0; 168 END IF; 169 RETURN ARRAY[z11x, z11y, z11w, z11h, z12x, z12y, z12w, z12h, z13x, z13y, z13w, z13h, z14x, z14y, z14w, z14h, z21x, z21y, z21w, z21h, z22x, z22y, z22w, z22h, z23x, z23y, z23w, z23h, z24x, z24y, z24w, z24h, zcx, zcy, zcw, zch]; 170 END; 171 $$ 172 LANGUAGE 'plpgsql'; 173 174--DROP FUNCTION _MapAlgebraPartsGeom(r1x int, r1y int, r1w int, r1h int, r2x int, r2y int, r2w int, r2h int); 175CREATE OR REPLACE FUNCTION _MapAlgebraPartsGeom(nx int, 176 ny int, 177 r1x int, 178 r1y int, 179 r1w int, 180 r1h int, 181 r2x int, 182 r2y int, 183 r2w int, 184 r2h int) 185 RETURNS SETOF geometry AS 186 $$ 187 DECLARE 188 BEGIN 189 RETURN NEXT ST_MakeBox2D(ST_Point(10 * ny + r1x, -10 * nx + 5 - r1y), ST_Point(10 * ny + r1x + r1w, -10 * nx + 5 - (r1y + r1h)))::geometry; 190 RETURN NEXT ST_MakeBox2D(ST_Point(10 * ny + r2x, -10 * nx + 5 - r2y), ST_Point(10 * ny + r2x + r2w, -10 * nx + 5 - (r2y + r2h)))::geometry; 191 RETURN; 192 END; 193 $$ 194 LANGUAGE 'plpgsql'; 195 196CREATE OR REPLACE FUNCTION _MapAlgebraAllPartsGeom(r1x int, 197 r1y int, 198 r1w int, 199 r1h int, 200 r2x int, 201 r2y int, 202 r2w int, 203 r2h int) 204 RETURNS SETOF geometry AS 205 $$ 206 DECLARE 207 z int[]; 208 BEGIN 209 z := _MapAlgebraParts(r1x, r1y, r1w, r1h, r2x, r2y, r2w, r2h); 210 RETURN NEXT ST_MakeBox2D(ST_Point(z[1], z[2]), ST_Point(z[1] + z[3], z[2] + z[4]))::geometry; 211 RETURN NEXT ST_MakeBox2D(ST_Point(z[5], z[6]), ST_Point(z[5] + z[7], z[6] + z[8]))::geometry; 212 RETURN NEXT ST_MakeBox2D(ST_Point(z[9], z[10]), ST_Point(z[9] + z[11], z[10] + z[12]))::geometry; 213 RETURN NEXT ST_MakeBox2D(ST_Point(z[13], z[14]), ST_Point(z[13] + z[15], z[14] + z[16]))::geometry; 214 215 RETURN NEXT ST_MakeBox2D(ST_Point(z[17], z[18]), ST_Point(z[17] + z[19], z[18] + z[20]))::geometry; 216 RETURN NEXT ST_MakeBox2D(ST_Point(z[21], z[22]), ST_Point(z[21] + z[23], z[22] + z[24]))::geometry; 217 RETURN NEXT ST_MakeBox2D(ST_Point(z[25], z[26]), ST_Point(z[25] + z[27], z[26] + z[28]))::geometry; 218 RETURN NEXT ST_MakeBox2D(ST_Point(z[29], z[30]), ST_Point(z[29] + z[31], z[30] + z[32]))::geometry; 219 220 RETURN NEXT ST_MakeBox2D(ST_Point(z[33], z[34]), ST_Point(z[33] + z[35], z[34] + z[36]))::geometry; 221 RETURN; 222 END; 223 $$ 224 LANGUAGE 'plpgsql'; 225 226SELECT asbinary(_MapAlgebraAllPartsGeom(0, 0, 1, 1, 1, 0, 1, 1)) 227 228CREATE OR REPLACE FUNCTION X1W1X2W2() 229 RETURNS SETOF record AS 230 $$ 231 DECLARE 232 x1w1x2w2 record; 233 BEGIN 234 x1w1x2w2 := (0, 3-4, 2, 0-4, 2); 235 RETURN NEXT x1w1x2w2; 236 x1w1x2w2 := (1, 2-4, 3, 0-4, 2); 237 RETURN NEXT x1w1x2w2; 238 x1w1x2w2 := (2, 2-4, 3, 0-4, 3); 239 RETURN NEXT x1w1x2w2; 240 x1w1x2w2 := (3, 2-4, 3, 0-4, 5); 241 RETURN NEXT x1w1x2w2; 242 x1w1x2w2 := (4, 1-4, 3, 0-4, 5); 243 RETURN NEXT x1w1x2w2; 244 x1w1x2w2 := (5, 1-4, 4, 1-4, 2); 245 RETURN NEXT x1w1x2w2; 246 x1w1x2w2 := (6, 1-4, 3, 1-4, 3); 247 RETURN NEXT x1w1x2w2; 248 x1w1x2w2 := (7, 1-4, 2, 1-4, 4); 249 RETURN NEXT x1w1x2w2; 250 x1w1x2w2 := (8, 0-4, 5, 1-4, 3); 251 RETURN NEXT x1w1x2w2; 252 x1w1x2w2 := (9, 0-4, 5, 1-4, 4); 253 RETURN NEXT x1w1x2w2; 254 x1w1x2w2 := (10, 0-4, 3, 2-4, 3); 255 RETURN NEXT x1w1x2w2; 256 x1w1x2w2 := (11, 0-4, 3, 3-4, 2); 257 RETURN NEXT x1w1x2w2; 258 x1w1x2w2 := (12, 0-4, 2, 3-4, 2); 259 RETURN NEXT x1w1x2w2; 260 RETURN; 261 END; 262 $$ 263 LANGUAGE 'plpgsql'; 264 265CREATE OR REPLACE FUNCTION Y1H1Y2H2() 266 RETURNS SETOF record AS 267 $$ 268 DECLARE 269 y1h1y2h2 record; 270 BEGIN 271 y1h1y2h2 := (0, 3-4, 2, 0-4, 2); 272 RETURN NEXT y1h1y2h2; 273 y1h1y2h2 := (1, 2-4, 3, 0-4, 2); 274 RETURN NEXT y1h1y2h2; 275 y1h1y2h2 := (2, 2-4, 3, 0-4, 3); 276 RETURN NEXT y1h1y2h2; 277 y1h1y2h2 := (3, 2-4, 3, 0-4, 5); 278 RETURN NEXT y1h1y2h2; 279 y1h1y2h2 := (4, 1-4, 3, 0-4, 5); 280 RETURN NEXT y1h1y2h2; 281 y1h1y2h2 := (5, 1-4, 4, 1-4, 2); 282 RETURN NEXT y1h1y2h2; 283 y1h1y2h2 := (6, 1-4, 3, 1-4, 3); 284 RETURN NEXT y1h1y2h2; 285 y1h1y2h2 := (7, 1-4, 2, 1-4, 4); 286 RETURN NEXT y1h1y2h2; 287 y1h1y2h2 := (8, 0-4, 5, 1-4, 3); 288 RETURN NEXT y1h1y2h2; 289 y1h1y2h2 := (9, 0-4, 5, 1-4, 4); 290 RETURN NEXT y1h1y2h2; 291 y1h1y2h2 := (10, 0-4, 3, 2-4, 3); 292 RETURN NEXT y1h1y2h2; 293 y1h1y2h2 := (11, 0-4, 3, 3-4, 2); 294 RETURN NEXT y1h1y2h2; 295 y1h1y2h2 := (12, 0-4, 2, 3-4, 2); 296 RETURN NEXT y1h1y2h2; 297 RETURN; 298 END; 299 $$ 300 LANGUAGE 'plpgsql'; 301 302_MapAlgebraParts(r1x, r1y, r1w, r1h, r2x, r2y, r2w, r2h) 303 304SELECT nx, x1, w1, x2, w2 FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int); 305 306SELECT nx, ny, x1, w1, x2, w2, y1, h1, y2, h2 307FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int); 308 309SELECT nx, ny, x1, w1, x2, w2, y1, h1, y2, h2, asbinary(_MapAlgebraPartsGeom(nx, ny, x1, y1, w1, h1, x2, y2, w2, h2)) 310FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int); 311 312SELECT asbinary(_MapAlgebraPartsGeom(nx, ny, x1, y1, w1, h1, x2, y2, w2, h2)) 313FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int); 314 315-- First series of zones covering raster 1 316SELECT nx, ny, map[1], map[2], map[3], map[4], 317 asbinary(ST_MakeBox2D(ST_Point(10 * ny + map[1], -10 * nx + 5 - map[2]), ST_Point(10 * ny + map[1] + map[3], -10 * nx + 5 - (map[2] + map[4])))::geometry) 318FROM ( 319 SELECT nx, ny, x1, y1, w1, h1, x2, y2, w2, h2, _MapAlgebraParts(x1, y1, w1, h1, x2, y2, w2, h2) as map 320 FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int) 321 ) as foo; 322 323-- Second series of zones covering raster 1 324SELECT nx, ny, map[5], map[6], map[7], map[8], 325 asbinary(ST_Point(10 * ny + map[5], -10 * nx + 5 - map[6])::geometry) 326FROM ( 327 SELECT nx, ny, x1, y1, w1, h1, x2, y2, w2, h2, _MapAlgebraParts(x1, y1, w1, h1, x2, y2, w2, h2) as map 328 FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int) 329 ) as foo; 330 331SELECT nx, ny, map[5], map[6], map[7], map[8], 332 asbinary(ST_MakeBox2D(ST_Point(10 * ny + map[5], -10 * nx + 5 - map[6]), ST_Point(10 * ny + map[5] + map[7], -10 * nx + 5 - (map[6] + map[8])))::geometry) 333FROM ( 334 SELECT nx, ny, x1, y1, w1, h1, x2, y2, w2, h2, _MapAlgebraParts(x1, y1, w1, h1, x2, y2, w2, h2) as map 335 FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int) 336 ) as foo; 337 338-- Third series of zones covering raster 1 339SELECT nx, ny, map[9], map[10], map[11], map[12], 340 asbinary(ST_MakeBox2D(ST_Point(10 * ny + map[9], -10 * nx + 5 - map[10]), ST_Point(10 * ny + map[9] + map[11], -10 * nx + 5 - (map[10] + map[12])))::geometry) 341FROM ( 342 SELECT nx, ny, x1, y1, w1, h1, x2, y2, w2, h2, _MapAlgebraParts(x1, y1, w1, h1, x2, y2, w2, h2) as map 343 FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int) 344 ) as foo; 345 346-- Fourth series of zones covering raster 1 347SELECT nx, ny, map[13], map[14], map[15], map[16], 348 asbinary(ST_MakeBox2D(ST_Point(10 * ny + map[13], -10 * nx + 5 - map[14]), ST_Point(10 * ny + map[13] + map[15], -10 * nx + 5 - (map[14] + map[16])))::geometry) 349FROM ( 350 SELECT nx, ny, x1, y1, w1, h1, x2, y2, w2, h2, _MapAlgebraParts(x1, y1, w1, h1, x2, y2, w2, h2) as map 351 FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int) 352 ) as foo; 353 354-- First series of zones covering raster 2 355SELECT nx, ny, map[17], map[18], map[19], map[20], 356 asbinary(ST_MakeBox2D(ST_Point(10 * ny + map[17], -10 * nx + 5 - map[18]), ST_Point(10 * ny + map[17] + map[19], -10 * nx + 5 - (map[18] + map[20])))::geometry) 357FROM ( 358 SELECT nx, ny, x1, y1, w1, h1, x2, y2, w2, h2, _MapAlgebraParts(x1, y1, w1, h1, x2, y2, w2, h2) as map 359 FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int) 360 ) as foo; 361 362-- Second series of zones covering raster 2 363SELECT nx, ny, map[21], map[22], map[23], map[24], 364 asbinary(ST_Point(10 * ny + map[21], -10 * nx + 5 - map[22])::geometry) 365FROM ( 366 SELECT nx, ny, x1, y1, w1, h1, x2, y2, w2, h2, _MapAlgebraParts(x1, y1, w1, h1, x2, y2, w2, h2) as map 367 FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int) 368 ) as foo; 369 370SELECT nx, ny, map[21], map[22], map[23], map[24], 371 asbinary(ST_MakeBox2D(ST_Point(10 * ny + map[21], -10 * nx + 5 - map[22]), ST_Point(10 * ny + map[21] + map[23], -10 * nx + 5 - (map[22] + map[24])))::geometry) 372FROM ( 373 SELECT nx, ny, x1, y1, w1, h1, x2, y2, w2, h2, _MapAlgebraParts(x1, y1, w1, h1, x2, y2, w2, h2) as map 374 FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int) 375 ) as foo; 376 377-- Third series of zones covering raster 2 378SELECT nx, ny, map[25], map[26], map[27], map[28], 379 asbinary(ST_MakeBox2D(ST_Point(10 * ny + map[25], -10 * nx + 5 - map[26]), ST_Point(10 * ny + map[25] + map[27], -10 * nx + 5 - (map[26] + map[28])))::geometry) 380FROM ( 381 SELECT nx, ny, x1, y1, w1, h1, x2, y2, w2, h2, _MapAlgebraParts(x1, y1, w1, h1, x2, y2, w2, h2) as map 382 FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int) 383 ) as foo; 384 385-- Fourth series of zones covering raster 2 386SELECT nx, ny, map[29], map[30], map[31], map[32], 387 asbinary(ST_MakeBox2D(ST_Point(10 * ny + map[28], -10 * nx + 5 - map[29]), ST_Point(10 * ny + map[28] + map[30], -10 * nx + 5 - (map[29] + map[31])))::geometry) 388FROM ( 389 SELECT nx, ny, x1, y1, w1, h1, x2, y2, w2, h2, _MapAlgebraParts(x1, y1, w1, h1, x2, y2, w2, h2) as map 390 FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int) 391 ) as foo; 392 393-- Common zone 394SELECT nx, ny, map[33], map[34], map[35], map[36], 395 asbinary(ST_MakeBox2D(ST_Point(10 * ny + map[33], -10 * nx + 5 - map[34]), ST_Point(10 * ny + map[33] + map[35], -10 * nx + 5 - (map[34] + map[36])))::geometry) 396FROM ( 397 SELECT nx, ny, x1, y1, w1, h1, x2, y2, w2, h2, _MapAlgebraParts(x1, y1, w1, h1, x2, y2, w2, h2) as map 398 FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int) 399 ) as foo;