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 2016 Sandro Santilli <strk@kbt.io>
22 *
23 **********************************************************************/
24
25 #include "../postgis_config.h"
26 /*#define POSTGIS_DEBUG_LEVEL 4*/
27 #include "lwgeom_geos.h"
28 #include "liblwgeom_internal.h"
29
30 #include <string.h>
31 #include <assert.h>
32
33 LWGEOM* lwgeom_wrapx(const LWGEOM* lwgeom_in, double cutx, double amount);
34 static LWCOLLECTION* lwcollection_wrapx(const LWCOLLECTION* lwcoll_in, double cutx, double amount);
35
36 static LWGEOM*
lwgeom_split_wrapx(const LWGEOM * geom_in,double cutx,double amount)37 lwgeom_split_wrapx(const LWGEOM* geom_in, double cutx, double amount)
38 {
39 LWGEOM *blade, *split;
40 POINTARRAY *bladepa;
41 POINT4D pt;
42 const GBOX *box_in;
43 AFFINE affine = {
44 1, 0, 0,
45 0, 1, 0,
46 0, 0, 1,
47 amount, 0, 0,
48 };
49
50 /* Extract box */
51 /* TODO: check if the bbox should be force-recomputed */
52 box_in = lwgeom_get_bbox(geom_in);
53 if ( ! box_in ) {
54 /* must be empty */
55 return lwgeom_clone_deep(geom_in);
56 }
57
58 LWDEBUGF(2, "BOX X range is %g..%g, cutx:%g, amount:%g", box_in->xmin, box_in->xmax, cutx, amount);
59
60 /* Check if geometry is fully on the side needing shift */
61 if ( ( amount < 0 && box_in->xmin >= cutx ) || ( amount > 0 && box_in->xmax <= cutx ) )
62 {
63 split = lwgeom_clone_deep(geom_in);
64 lwgeom_affine(split, &affine);
65 LWDEBUGG(2, split, "returning the translated geometry");
66 return split;
67 }
68
69 /* Check if geometry is fully on the side needing no shift */
70 if ( ( amount < 0 && box_in->xmax <= cutx ) || ( amount > 0 && box_in->xmin >= cutx ) )
71 {
72 split = lwgeom_clone_deep(geom_in);
73 LWDEBUGG(2, split, "returning the cloned geometry");
74 return split;
75 }
76
77 /* We need splitting here */
78
79 /* construct blade */
80 bladepa = ptarray_construct(0, 0, 2);
81 pt.x = cutx;
82 pt.y = box_in->ymin - 1;
83 ptarray_set_point4d(bladepa, 0, &pt);
84 pt.y = box_in->ymax + 1;
85 ptarray_set_point4d(bladepa, 1, &pt);
86 blade = lwline_as_lwgeom(lwline_construct(geom_in->srid, NULL, bladepa));
87
88 LWDEBUG(2, "splitting the geometry");
89
90 /* split by blade */
91 split = lwgeom_split(geom_in, blade);
92 lwgeom_free(blade);
93 if ( ! split ) {
94 lwerror("%s:%d - lwgeom_split_wrapx: %s", __FILE__, __LINE__, lwgeom_geos_errmsg);
95 return NULL;
96 }
97 LWDEBUGG(2, split, "split geometry");
98
99
100 /* iterate over components, translate if needed */
101 const LWCOLLECTION *col = lwgeom_as_lwcollection(split);
102 if ( ! col ) {
103 /* not split, this is unexpected */
104 lwnotice("WARNING: unexpected lack of split in lwgeom_split_wrapx");
105 return lwgeom_clone_deep(geom_in);
106 }
107 LWCOLLECTION *col_out = lwcollection_wrapx(col, cutx, amount);
108 lwgeom_free(split);
109
110 /* unary-union the result (homogenize too ?) */
111 LWGEOM* out = lwgeom_unaryunion(lwcollection_as_lwgeom(col_out));
112 LWDEBUGF(2, "col_out:%p, unaryunion_out:%p", col_out, out);
113 LWDEBUGG(2, out, "unary-unioned");
114
115 lwcollection_free(col_out);
116
117 return out;
118 }
119
120 static LWCOLLECTION*
lwcollection_wrapx(const LWCOLLECTION * lwcoll_in,double cutx,double amount)121 lwcollection_wrapx(const LWCOLLECTION* lwcoll_in, double cutx, double amount)
122 {
123 LWGEOM** wrap_geoms;
124 LWCOLLECTION* out;
125 uint32_t i;
126 int outtype = lwcoll_in->type;
127
128 wrap_geoms = lwalloc(lwcoll_in->ngeoms * sizeof(LWGEOM*));
129 if ( ! wrap_geoms )
130 {
131 lwerror("Out of virtual memory");
132 return NULL;
133 }
134
135 for (i=0; i<lwcoll_in->ngeoms; ++i)
136 {
137 LWDEBUGF(3, "Wrapping collection element %d", i);
138 wrap_geoms[i] = lwgeom_wrapx(lwcoll_in->geoms[i], cutx, amount);
139 /* an exception should prevent this from ever returning NULL */
140 if ( ! wrap_geoms[i] ) {
141 uint32_t j;
142 lwnotice("Error wrapping geometry, cleaning up");
143 for (j = 0; j < i; j++)
144 {
145 lwnotice("cleaning geometry %d (%p)", j, wrap_geoms[j]);
146 lwgeom_free(wrap_geoms[j]);
147 }
148 lwfree(wrap_geoms);
149 lwnotice("cleanup complete");
150 return NULL;
151 }
152 if ( outtype != COLLECTIONTYPE ) {
153 if ( MULTITYPE[wrap_geoms[i]->type] != outtype )
154 {
155 outtype = COLLECTIONTYPE;
156 }
157 }
158 }
159
160 /* Now wrap_geoms has wrap_geoms_size geometries */
161 out = lwcollection_construct(outtype, lwcoll_in->srid, NULL,
162 lwcoll_in->ngeoms, wrap_geoms);
163
164 return out;
165 }
166
167 /* exported */
168 LWGEOM*
lwgeom_wrapx(const LWGEOM * lwgeom_in,double cutx,double amount)169 lwgeom_wrapx(const LWGEOM* lwgeom_in, double cutx, double amount)
170 {
171 /* Nothing to wrap in an empty geom */
172 if ( lwgeom_is_empty(lwgeom_in) )
173 {
174 LWDEBUG(2, "geom is empty, cloning");
175 return lwgeom_clone_deep(lwgeom_in);
176 }
177
178 /* Nothing to wrap if shift amount is zero */
179 if ( amount == 0 )
180 {
181 LWDEBUG(2, "amount is zero, cloning");
182 return lwgeom_clone_deep(lwgeom_in);
183 }
184
185 switch (lwgeom_in->type)
186 {
187 case LINETYPE:
188 case POLYGONTYPE:
189 LWDEBUG(2, "split-wrapping line or polygon");
190 return lwgeom_split_wrapx(lwgeom_in, cutx, amount);
191
192 case POINTTYPE:
193 {
194 const LWPOINT *pt = lwgeom_as_lwpoint(lwgeom_clone_deep(lwgeom_in));
195 POINT4D pt4d;
196 getPoint4d_p(pt->point, 0, &pt4d);
197
198 LWDEBUGF(2, "POINT X is %g, cutx:%g, amount:%g", pt4d.x, cutx, amount);
199
200 if ( ( amount < 0 && pt4d.x > cutx ) || ( amount > 0 && pt4d.x < cutx ) )
201 {
202 pt4d.x += amount;
203 ptarray_set_point4d(pt->point, 0, &pt4d);
204 }
205 return lwpoint_as_lwgeom(pt);
206 }
207
208 case MULTIPOINTTYPE:
209 case MULTIPOLYGONTYPE:
210 case MULTILINETYPE:
211 case COLLECTIONTYPE:
212 LWDEBUG(2, "collection-wrapping multi");
213 return lwcollection_as_lwgeom(
214 lwcollection_wrapx((const LWCOLLECTION*)lwgeom_in, cutx, amount)
215 );
216
217 default:
218 lwerror("Wrapping of %s geometries is unsupported",
219 lwtype_name(lwgeom_in->type));
220 return NULL;
221 }
222
223 }
224