1 /**********************************************************************
2 *
3 * PostGIS - Spatial Types for PostgreSQL
4 * http://postgis.net
5 *
6 * PostGIS is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * PostGIS is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with PostGIS. If not, see <http://www.gnu.org/licenses/>.
18 *
19 **********************************************************************
20 *
21 * Copyright 2020 Paul Ramsey <pramsey@cleverelephant.ca>
22 *
23 **********************************************************************/
24
25 #include "postgres.h"
26 #include "fmgr.h"
27 #include "funcapi.h"
28 #include "access/htup_details.h"
29
30 #include "../postgis_config.h"
31 #include "lwgeom_pg.h"
32 #include "liblwgeom.h"
33
34 #include <math.h>
35 #include <float.h>
36 #include <string.h>
37 #include <stdio.h>
38
39
40 Datum ST_Hexagon(PG_FUNCTION_ARGS);
41 Datum ST_Square(PG_FUNCTION_ARGS);
42 Datum ST_ShapeGrid(PG_FUNCTION_ARGS);
43
44 /* ********* ********* ********* ********* ********* ********* ********* ********* */
45
46 typedef enum
47 {
48 SHAPE_SQUARE,
49 SHAPE_HEXAGON,
50 SHAPE_TRIANGLE
51 } GeometryShape;
52
53 typedef struct GeometryGridState
54 {
55 GeometryShape cell_shape;
56 bool done;
57 GBOX bounds;
58 int32_t srid;
59 double size;
60 int32_t i, j;
61 }
62 GeometryGridState;
63
64 /* ********* ********* ********* ********* ********* ********* ********* ********* */
65
66 /* Ratio of hexagon edge length to edge->center vertical distance = cos(M_PI/6.0) */
67 #define H 0.8660254037844387
68
69 /* Build origin hexagon centered around origin point */
70 static const double hex_x[] = {-1.0, -0.5, 0.5, 1.0, 0.5, -0.5, -1.0};
71 static const double hex_y[] = {0.0, -0.5, -0.5, 0.0, 0.5, 0.5, 0.0};
72
73 static LWGEOM *
hexagon(double origin_x,double origin_y,double size,int cell_i,int cell_j,int32_t srid)74 hexagon(double origin_x, double origin_y, double size, int cell_i, int cell_j, int32_t srid)
75 {
76 POINTARRAY **ppa = lwalloc(sizeof(POINTARRAY*));
77 POINTARRAY *pa = ptarray_construct(0, 0, 7);
78
79 for (uint32_t i = 0; i < 7; ++i)
80 {
81 double height = size * 2 * H;
82 POINT4D pt;
83 pt.x = origin_x + size * (1.5 * cell_i + hex_x[i]);
84 pt.y = origin_y + height * (cell_j + 0.5 * (abs(cell_i) % 2) + hex_y[i]);
85 ptarray_set_point4d(pa, i, &pt);
86 }
87
88 ppa[0] = pa;
89 return lwpoly_as_lwgeom(lwpoly_construct(srid, NULL, 1 /* nrings */, ppa));
90 }
91
92
93 typedef struct HexagonGridState
94 {
95 GeometryShape cell_shape;
96 bool done;
97 GBOX bounds;
98 int32_t srid;
99 double size;
100 int32_t i, j;
101 int32_t column_min, column_max;
102 int32_t row_min_odd, row_max_odd;
103 int32_t row_min_even, row_max_even;
104 }
105 HexagonGridState;
106
107 static HexagonGridState *
hexagon_grid_state(double size,const GBOX * gbox,int32_t srid)108 hexagon_grid_state(double size, const GBOX *gbox, int32_t srid)
109 {
110 HexagonGridState *state = palloc0(sizeof(HexagonGridState));
111 double col_width = 1.5 * size;
112 double row_height = size * 2 * H;
113
114 /* fill in state */
115 state->cell_shape = SHAPE_HEXAGON;
116 state->size = size;
117 state->srid = srid;
118 state->done = false;
119 state->bounds = *gbox;
120
121 /* Column address is just width / column width, with an */
122 /* adjustment to account for partial overlaps */
123 state->column_min = floor(gbox->xmin / col_width);
124 if(gbox->xmin - state->column_min * col_width > size)
125 state->column_min++;
126
127 state->column_max = ceil(gbox->xmax / col_width);
128 if(state->column_max * col_width - gbox->xmax > size)
129 state->column_max--;
130
131 /* Row address range depends on the odd/even column we're on */
132 state->row_min_even = floor(gbox->ymin/row_height + 0.5);
133 state->row_max_even = floor(gbox->ymax/row_height + 0.5);
134 state->row_min_odd = floor(gbox->ymin/row_height);
135 state->row_max_odd = floor(gbox->ymax/row_height);
136
137 /* Set initial state */
138 state->i = state->column_min;
139 state->j = (state->i % 2) ? state->row_min_odd : state->row_min_even;
140
141 return state;
142 }
143
144 static void
hexagon_state_next(HexagonGridState * state)145 hexagon_state_next(HexagonGridState *state)
146 {
147 if (!state || state->done) return;
148 /* Move up one row */
149 state->j++;
150 /* Off the end, increment column counter, reset row counter back to (appropriate) minimum */
151 if (state->j > ((state->i % 2) ? state->row_max_odd : state->row_max_even))
152 {
153 state->i++;
154 state->j = ((state->i % 2) ? state->row_min_odd : state->row_min_even);
155 }
156 /* Column counter over max means we have used up all combinations */
157 if (state->i > state->column_max)
158 {
159 state->done = true;
160 }
161 }
162
163 /* ********* ********* ********* ********* ********* ********* ********* ********* */
164
165 static LWGEOM *
square(double origin_x,double origin_y,double size,int cell_i,int cell_j,int32_t srid)166 square(double origin_x, double origin_y, double size, int cell_i, int cell_j, int32_t srid)
167 {
168 double ll_x = origin_x + (size * cell_i);
169 double ll_y = origin_y + (size * cell_j);
170 double ur_x = origin_x + (size * (cell_i + 1));
171 double ur_y = origin_y + (size * (cell_j + 1));
172 return (LWGEOM*)lwpoly_construct_envelope(srid, ll_x, ll_y, ur_x, ur_y);
173 }
174
175 typedef struct SquareGridState
176 {
177 GeometryShape cell_shape;
178 bool done;
179 GBOX bounds;
180 int32_t srid;
181 double size;
182 int32_t i, j;
183 int32_t column_min, column_max;
184 int32_t row_min, row_max;
185 }
186 SquareGridState;
187
188 static SquareGridState *
square_grid_state(double size,const GBOX * gbox,int32_t srid)189 square_grid_state(double size, const GBOX *gbox, int32_t srid)
190 {
191 SquareGridState *state = palloc0(sizeof(SquareGridState));
192
193 /* fill in state */
194 state->cell_shape = SHAPE_SQUARE;
195 state->size = size;
196 state->srid = srid;
197 state->done = false;
198 state->bounds = *gbox;
199
200 state->column_min = floor(gbox->xmin / size);
201 state->column_max = floor(gbox->xmax / size);
202 state->row_min = floor(gbox->ymin / size);
203 state->row_max = floor(gbox->ymax / size);
204 state->i = state->column_min;
205 state->j = state->row_min;
206 return state;
207 }
208
209 static void
square_state_next(SquareGridState * state)210 square_state_next(SquareGridState *state)
211 {
212 if (!state || state->done) return;
213 /* Move up one row */
214 state->j++;
215 /* Off the end, increment column counter, reset row counter back to (appropriate) minimum */
216 if (state->j > state->row_max)
217 {
218 state->i++;
219 state->j = state->row_min;
220 }
221 /* Column counter over max means we have used up all combinations */
222 if (state->i > state->column_max)
223 {
224 state->done = true;
225 }
226 }
227
228 /**
229 * ST_HexagonGrid(size float8, bounds geometry)
230 * ST_SquareGrid(size float8, bounds geometry)
231 */
232 PG_FUNCTION_INFO_V1(ST_ShapeGrid);
ST_ShapeGrid(PG_FUNCTION_ARGS)233 Datum ST_ShapeGrid(PG_FUNCTION_ARGS)
234 {
235 FuncCallContext *funcctx;
236
237 GSERIALIZED *gbounds;
238 GeometryGridState *state;
239
240 LWGEOM *lwgeom;
241 bool isnull[3] = {0,0,0}; /* needed to say no value is null */
242 Datum tuple_arr[3]; /* used to construct the composite return value */
243 HeapTuple tuple;
244 Datum result; /* the actual composite return value */
245
246 if (SRF_IS_FIRSTCALL())
247 {
248 MemoryContext oldcontext;
249 const char *func_name;
250 char gbounds_is_empty;
251 GBOX bounds;
252 double size;
253 funcctx = SRF_FIRSTCALL_INIT();
254
255 /* Get a local copy of the bounds we are going to fill with shapes */
256 gbounds = PG_GETARG_GSERIALIZED_P(1);
257 size = PG_GETARG_FLOAT8(0);
258
259 gbounds_is_empty = (gserialized_get_gbox_p(gbounds, &bounds) == LW_FAILURE);
260
261 /* quick opt-out if we get nonsensical inputs */
262 if (size <= 0.0 || gbounds_is_empty)
263 {
264 funcctx = SRF_PERCALL_SETUP();
265 SRF_RETURN_DONE(funcctx);
266 }
267
268 oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
269
270 /*
271 * Support both hexagon and square grids with one function,
272 * by checking the calling signature up front.
273 */
274 func_name = get_func_name(fcinfo->flinfo->fn_oid);
275 if (strcmp(func_name, "st_hexagongrid") == 0)
276 {
277 state = (GeometryGridState*)hexagon_grid_state(size, &bounds, gserialized_get_srid(gbounds));
278 }
279 else if (strcmp(func_name, "st_squaregrid") == 0)
280 {
281 state = (GeometryGridState*)square_grid_state(size, &bounds, gserialized_get_srid(gbounds));
282 }
283 else
284 {
285 ereport(ERROR, (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
286 errmsg("%s called from unsupported functional context '%s'", __func__, func_name)));
287 }
288
289 funcctx->user_fctx = state;
290
291 /* get tuple description for return type */
292 if (get_call_result_type(fcinfo, 0, &funcctx->tuple_desc) != TYPEFUNC_COMPOSITE)
293 {
294 ereport(ERROR, (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
295 errmsg("set-valued function called in context that cannot accept a set")));
296 }
297
298 BlessTupleDesc(funcctx->tuple_desc);
299 MemoryContextSwitchTo(oldcontext);
300 }
301
302 /* stuff done on every call of the function */
303 funcctx = SRF_PERCALL_SETUP();
304
305 /* get state */
306 state = funcctx->user_fctx;
307
308 /* Stop when we've used up all the grid squares */
309 if (state->done)
310 {
311 SRF_RETURN_DONE(funcctx);
312 }
313
314 /* Store grid cell coordinates */
315 tuple_arr[1] = Int32GetDatum(state->i);
316 tuple_arr[2] = Int32GetDatum(state->j);
317
318 /* Generate geometry */
319 switch (state->cell_shape)
320 {
321 case SHAPE_HEXAGON:
322 lwgeom = hexagon(0.0, 0.0, state->size, state->i, state->j, state->srid);
323 /* Increment to next cell */
324 hexagon_state_next((HexagonGridState*)state);
325 break;
326 case SHAPE_SQUARE:
327 lwgeom = square(0.0, 0.0, state->size, state->i, state->j, state->srid);
328 /* Increment to next cell */
329 square_state_next((SquareGridState*)state);
330 break;
331 default:
332 ereport(ERROR, (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
333 errmsg("%s called from with unsupported shape '%d'", __func__, state->cell_shape)));
334 }
335
336 tuple_arr[0] = PointerGetDatum(geometry_serialize(lwgeom));
337 lwfree(lwgeom);
338
339 /* Form tuple and return */
340 tuple = heap_form_tuple(funcctx->tuple_desc, tuple_arr, isnull);
341 result = HeapTupleGetDatum(tuple);
342 SRF_RETURN_NEXT(funcctx, result);
343 }
344
345 /**
346 * ST_Hexagon(size double, cell_i integer default 0, cell_j integer default 0, origin geometry default 'POINT(0 0)')
347 */
348 PG_FUNCTION_INFO_V1(ST_Hexagon);
ST_Hexagon(PG_FUNCTION_ARGS)349 Datum ST_Hexagon(PG_FUNCTION_ARGS)
350 {
351 LWPOINT *lwpt;
352 double size = PG_GETARG_FLOAT8(0);
353 GSERIALIZED *ghex;
354 int cell_i = PG_GETARG_INT32(1);
355 int cell_j = PG_GETARG_INT32(2);
356 GSERIALIZED *gorigin = PG_GETARG_GSERIALIZED_P(3);
357 LWGEOM *lwhex;
358 LWGEOM *lworigin = lwgeom_from_gserialized(gorigin);
359
360 if (lwgeom_is_empty(lworigin))
361 {
362 elog(ERROR, "%s: origin point is empty", __func__);
363 PG_RETURN_NULL();
364 }
365
366 lwpt = lwgeom_as_lwpoint(lworigin);
367 if (!lwpt)
368 {
369 elog(ERROR, "%s: origin argument is not a point", __func__);
370 PG_RETURN_NULL();
371 }
372
373 lwhex = hexagon(lwpoint_get_x(lwpt), lwpoint_get_y(lwpt),
374 size, cell_i, cell_j,
375 lwgeom_get_srid(lworigin));
376
377 ghex = geometry_serialize(lwhex);
378 PG_FREE_IF_COPY(gorigin, 3);
379 PG_RETURN_POINTER(ghex);
380 }
381
382
383 /**
384 * ST_Square(size double, cell_i integer, cell_j integer, origin geometry default 'POINT(0 0)')
385 */
386 PG_FUNCTION_INFO_V1(ST_Square);
ST_Square(PG_FUNCTION_ARGS)387 Datum ST_Square(PG_FUNCTION_ARGS)
388 {
389 LWPOINT *lwpt;
390 double size = PG_GETARG_FLOAT8(0);
391 GSERIALIZED *gsqr;
392 int cell_i = PG_GETARG_INT32(1);
393 int cell_j = PG_GETARG_INT32(2);
394 GSERIALIZED *gorigin = PG_GETARG_GSERIALIZED_P(3);
395 LWGEOM *lwsqr;
396 LWGEOM *lworigin = lwgeom_from_gserialized(gorigin);
397
398 if (lwgeom_is_empty(lworigin))
399 {
400 elog(ERROR, "%s: origin point is empty", __func__);
401 PG_RETURN_NULL();
402 }
403
404 lwpt = lwgeom_as_lwpoint(lworigin);
405 if (!lwpt)
406 {
407 elog(ERROR, "%s: origin argument is not a point", __func__);
408 PG_RETURN_NULL();
409 }
410
411 lwsqr = square(lwpoint_get_x(lwpt), lwpoint_get_y(lwpt),
412 size, cell_i, cell_j,
413 lwgeom_get_srid(lworigin));
414
415 gsqr = geometry_serialize(lwsqr);
416 PG_FREE_IF_COPY(gorigin, 3);
417 PG_RETURN_POINTER(gsqr);
418 }
419