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;