1 /******************************************************************************
2 *
3 * Project: High Performance Image Reprojector
4 * Purpose: Implement cutline/blend mask generator.
5 * Author: Frank Warmerdam, warmerdam@pobox.com
6 *
7 ******************************************************************************
8 * Copyright (c) 2008, Frank Warmerdam <warmerdam@pobox.com>
9 * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
10 *
11 * Permission is hereby granted, free of charge, to any person obtaining a
12 * copy of this software and associated documentation files (the "Software"),
13 * to deal in the Software without restriction, including without limitation
14 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 * and/or sell copies of the Software, and to permit persons to whom the
16 * Software is furnished to do so, subject to the following conditions:
17 *
18 * The above copyright notice and this permission notice shall be included
19 * in all copies or substantial portions of the Software.
20 *
21 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 * DEALINGS IN THE SOFTWARE.
28 ****************************************************************************/
29
30 #include "cpl_port.h"
31 #include "gdalwarper.h"
32
33 #include <cmath>
34 #include <cstdio>
35 #include <cstring>
36 #include <algorithm>
37
38 #include "cpl_conv.h"
39 #include "cpl_error.h"
40 #include "cpl_string.h"
41 #include "gdal.h"
42 #include "gdal_alg.h"
43 #include "ogr_api.h"
44 #include "ogr_core.h"
45 #include "ogr_geometry.h"
46 #include "ogr_geos.h"
47
48 CPL_CVSID("$Id: gdalcutline.cpp ab04d3bd6e63f3824c0d4a3bd40e3e3f3d84740a 2020-10-09 01:31:35 +0200 Even Rouault $")
49
50 /************************************************************************/
51 /* BlendMaskGenerator() */
52 /************************************************************************/
53
54 #ifndef HAVE_GEOS
55
56 static CPLErr
BlendMaskGenerator(int,int,int,int,GByte *,float *,OGRGeometryH,double)57 BlendMaskGenerator( int /* nXOff */, int /* nYOff */,
58 int /* nXSize */, int /* nYSize */,
59 GByte * /* pabyPolyMask */,
60 float * /* pafValidityMask */,
61 OGRGeometryH /* hPolygon */,
62 double /* dfBlendDist */ )
63 {
64 CPLError(CE_Failure, CPLE_AppDefined,
65 "Blend distance support not available without the GEOS library.");
66 return CE_Failure;
67 }
68 #else
69 static CPLErr
70 BlendMaskGenerator( int nXOff, int nYOff, int nXSize, int nYSize,
71 GByte *pabyPolyMask, float *pafValidityMask,
72 OGRGeometryH hPolygon, double dfBlendDist )
73 {
74
75 /* -------------------------------------------------------------------- */
76 /* Convert the polygon into a collection of lines so that we */
77 /* measure distance from the edge even on the inside. */
78 /* -------------------------------------------------------------------- */
79 OGRGeometry *poLines =
80 OGRGeometryFactory::forceToMultiLineString(
81 reinterpret_cast<OGRGeometry *>(hPolygon)->clone() );
82
83 /* -------------------------------------------------------------------- */
84 /* Prepare a clipping polygon a bit bigger than the area of */
85 /* interest in the hopes of simplifying the cutline down to */
86 /* stuff that will be relevant for this area of interest. */
87 /* -------------------------------------------------------------------- */
88 CPLString osClipRectWKT;
89
90 osClipRectWKT.Printf( "POLYGON((%g %g,%g %g,%g %g,%g %g,%g %g))",
91 nXOff - (dfBlendDist + 1),
92 nYOff - (dfBlendDist + 1),
93 nXOff + nXSize + (dfBlendDist + 1),
94 nYOff - (dfBlendDist + 1),
95 nXOff + nXSize + (dfBlendDist + 1),
96 nYOff + nYSize + (dfBlendDist + 1),
97 nXOff - (dfBlendDist + 1),
98 nYOff + nYSize + (dfBlendDist + 1),
99 nXOff - (dfBlendDist + 1),
100 nYOff - (dfBlendDist + 1) );
101
102 OGRPolygon *poClipRect = nullptr;
103 OGRGeometryFactory::createFromWkt( osClipRectWKT.c_str(), nullptr,
104 reinterpret_cast<OGRGeometry**>(&poClipRect) );
105
106 if( poClipRect )
107 {
108 // If it does not intersect the polym, zero the mask and return.
109 if( !reinterpret_cast<OGRGeometry *>(hPolygon)->Intersects(poClipRect) )
110 {
111 memset( pafValidityMask, 0, sizeof(float) * nXSize * nYSize );
112
113 delete poLines;
114 delete poClipRect;
115
116 return CE_None;
117 }
118
119 // If it does not intersect the line at all, just return.
120 else if( !static_cast<OGRGeometry *>(poLines)->Intersects(poClipRect) )
121 {
122 delete poLines;
123 delete poClipRect;
124
125 return CE_None;
126 }
127
128 OGRGeometry *poClippedLines = poLines->Intersection(poClipRect);
129 delete poLines;
130 poLines = poClippedLines;
131 delete poClipRect;
132 }
133
134 /* -------------------------------------------------------------------- */
135 /* Convert our polygon into GEOS format, and compute an */
136 /* envelope to accelerate later distance operations. */
137 /* -------------------------------------------------------------------- */
138 OGREnvelope sEnvelope;
139 GEOSContextHandle_t hGEOSCtxt = OGRGeometry::createGEOSContext();
140 GEOSGeom poGEOSPoly = poLines->exportToGEOS(hGEOSCtxt);
141 OGR_G_GetEnvelope( hPolygon, &sEnvelope );
142
143 delete poLines;
144
145 // This check was already done in the calling
146 // function and should never be true.
147
148 // if( sEnvelope.MinY - dfBlendDist > nYOff+nYSize
149 // || sEnvelope.MaxY + dfBlendDist < nYOff
150 // || sEnvelope.MinX - dfBlendDist > nXOff+nXSize
151 // || sEnvelope.MaxX + dfBlendDist < nXOff )
152 // return CE_None;
153
154 const int iXMin =
155 std::max(0,
156 static_cast<int>(floor(sEnvelope.MinX - dfBlendDist - nXOff)));
157 const int iXMax =
158 std::min(nXSize,
159 static_cast<int>(ceil(sEnvelope.MaxX + dfBlendDist - nXOff)));
160 const int iYMin =
161 std::max(0,
162 static_cast<int>(floor(sEnvelope.MinY - dfBlendDist - nYOff)));
163 const int iYMax =
164 std::min(nYSize,
165 static_cast<int>(ceil(sEnvelope.MaxY + dfBlendDist - nYOff)));
166
167 /* -------------------------------------------------------------------- */
168 /* Loop over potential area within blend line distance, */
169 /* processing each pixel. */
170 /* -------------------------------------------------------------------- */
171 for( int iY = 0; iY < nYSize; iY++ )
172 {
173 double dfLastDist = 0.0;
174
175 for( int iX = 0; iX < nXSize; iX++ )
176 {
177 if( iX < iXMin || iX >= iXMax
178 || iY < iYMin || iY > iYMax
179 || dfLastDist > dfBlendDist + 1.5 )
180 {
181 if( pabyPolyMask[iX + iY * nXSize] == 0 )
182 pafValidityMask[iX + iY * nXSize] = 0.0;
183
184 dfLastDist -= 1.0;
185 continue;
186 }
187
188 CPLString osPointWKT;
189 osPointWKT.Printf( "POINT(%d.5 %d.5)", iX + nXOff, iY + nYOff );
190
191 GEOSGeom poGEOSPoint = GEOSGeomFromWKT_r( hGEOSCtxt, osPointWKT );
192
193 double dfDist = 0.0;
194 GEOSDistance_r( hGEOSCtxt, poGEOSPoly, poGEOSPoint, &dfDist );
195 GEOSGeom_destroy_r( hGEOSCtxt, poGEOSPoint );
196
197 dfLastDist = dfDist;
198
199 if( dfDist > dfBlendDist )
200 {
201 if( pabyPolyMask[iX + iY * nXSize] == 0 )
202 pafValidityMask[iX + iY * nXSize] = 0.0;
203
204 continue;
205 }
206
207 const double dfRatio =
208 pabyPolyMask[iX + iY * nXSize] == 0
209 ? 0.5 - (dfDist / dfBlendDist) * 0.5 // Outside.
210 : 0.5 + (dfDist / dfBlendDist) * 0.5; // Inside.
211
212 pafValidityMask[iX + iY * nXSize] *= static_cast<float>(dfRatio);
213 }
214 }
215
216 /* -------------------------------------------------------------------- */
217 /* Cleanup */
218 /* -------------------------------------------------------------------- */
219 GEOSGeom_destroy_r( hGEOSCtxt, poGEOSPoly );
220 OGRGeometry::freeGEOSContext( hGEOSCtxt );
221
222 return CE_None;
223
224 }
225 #endif // HAVE_GEOS
226
227 /************************************************************************/
228 /* CutlineTransformer() */
229 /* */
230 /* A simple transformer for the cutline that just offsets */
231 /* relative to the current chunk. */
232 /************************************************************************/
233
CutlineTransformer(void * pTransformArg,int bDstToSrc,int nPointCount,double * x,double * y,double *,int *)234 static int CutlineTransformer( void *pTransformArg,
235 int bDstToSrc,
236 int nPointCount,
237 double *x,
238 double *y,
239 double * /* z */,
240 int * /* panSuccess */ )
241 {
242 int nXOff = static_cast<int *>(pTransformArg)[0];
243 int nYOff = static_cast<int *>(pTransformArg)[1];
244
245 if( bDstToSrc )
246 {
247 nXOff *= -1;
248 nYOff *= -1;
249 }
250
251 for( int i = 0; i < nPointCount; i++ )
252 {
253 x[i] -= nXOff;
254 y[i] -= nYOff;
255 }
256
257 return TRUE;
258 }
259
260 /************************************************************************/
261 /* GDALWarpCutlineMasker() */
262 /* */
263 /* This function will generate a source mask based on a */
264 /* provided cutline, and optional blend distance. */
265 /************************************************************************/
266
267 CPLErr
GDALWarpCutlineMasker(void * pMaskFuncArg,int,GDALDataType,int nXOff,int nYOff,int nXSize,int nYSize,GByte **,int bMaskIsFloat,void * pValidityMask)268 GDALWarpCutlineMasker( void *pMaskFuncArg,
269 int /* nBandCount */,
270 GDALDataType /* eType */,
271 int nXOff, int nYOff, int nXSize, int nYSize,
272 GByte ** /*ppImageData */,
273 int bMaskIsFloat, void *pValidityMask )
274
275 {
276 if( nXSize < 1 || nYSize < 1 )
277 return CE_None;
278
279 /* -------------------------------------------------------------------- */
280 /* Do some minimal checking. */
281 /* -------------------------------------------------------------------- */
282 if( !bMaskIsFloat )
283 {
284 CPLAssert( false );
285 return CE_Failure;
286 }
287
288 GDALWarpOptions *psWO = static_cast<GDALWarpOptions *>(pMaskFuncArg);
289
290 if( psWO == nullptr || psWO->hCutline == nullptr )
291 {
292 CPLAssert( false );
293 return CE_Failure;
294 }
295
296 GDALDriverH hMemDriver = GDALGetDriverByName("MEM");
297 if( hMemDriver == nullptr )
298 {
299 CPLError(CE_Failure, CPLE_AppDefined,
300 "GDALWarpCutlineMasker needs MEM driver");
301 return CE_Failure;
302 }
303
304 /* -------------------------------------------------------------------- */
305 /* Check the polygon. */
306 /* -------------------------------------------------------------------- */
307 OGRGeometryH hPolygon = static_cast<OGRGeometryH>(psWO->hCutline);
308
309 if( wkbFlatten(OGR_G_GetGeometryType(hPolygon)) != wkbPolygon
310 && wkbFlatten(OGR_G_GetGeometryType(hPolygon)) != wkbMultiPolygon )
311 {
312 CPLError(CE_Failure, CPLE_NotSupported,
313 "Cutline should be a polygon or a multipolygon");
314 return CE_Failure;
315 }
316
317 OGREnvelope sEnvelope;
318 OGR_G_GetEnvelope( hPolygon, &sEnvelope );
319
320 float *pafMask = static_cast<float *>(pValidityMask);
321
322 if( sEnvelope.MaxX + psWO->dfCutlineBlendDist < nXOff
323 || sEnvelope.MinX - psWO->dfCutlineBlendDist > nXOff + nXSize
324 || sEnvelope.MaxY + psWO->dfCutlineBlendDist < nYOff
325 || sEnvelope.MinY - psWO->dfCutlineBlendDist > nYOff + nYSize )
326 {
327 // We are far from the blend line - everything is masked to zero.
328 // It would be nice to realize no work is required for this whole
329 // chunk!
330 memset( pafMask, 0, sizeof(float) * nXSize * nYSize );
331 return CE_None;
332 }
333
334 /* -------------------------------------------------------------------- */
335 /* Create a byte buffer into which we can burn the */
336 /* mask polygon and wrap it up as a memory dataset. */
337 /* -------------------------------------------------------------------- */
338 GByte *pabyPolyMask = static_cast<GByte *>(CPLCalloc(nXSize, nYSize));
339
340 char szDataPointer[100] = {};
341
342 // cppcheck-suppress redundantCopy
343 snprintf( szDataPointer, sizeof(szDataPointer), "DATAPOINTER=" );
344 CPLPrintPointer(
345 szDataPointer+strlen(szDataPointer),
346 pabyPolyMask,
347 static_cast<int>(sizeof(szDataPointer) - strlen(szDataPointer)) );
348
349 GDALDatasetH hMemDS = GDALCreate( hMemDriver, "warp_temp",
350 nXSize, nYSize, 0, GDT_Byte, nullptr );
351 char *apszOptions[] = { szDataPointer, nullptr };
352 GDALAddBand( hMemDS, GDT_Byte, apszOptions );
353
354 double adfGeoTransform[6] = { 0.0, 1.0, 0.0, 0.0, 0.0, 1.0 };
355 GDALSetGeoTransform( hMemDS, adfGeoTransform );
356
357 /* -------------------------------------------------------------------- */
358 /* Burn the polygon into the mask with 1.0 values. */
359 /* -------------------------------------------------------------------- */
360 int nTargetBand = 1;
361 double dfBurnValue = 255.0;
362 char **papszRasterizeOptions = nullptr;
363
364 if( CPLFetchBool( psWO->papszWarpOptions, "CUTLINE_ALL_TOUCHED", false ))
365 papszRasterizeOptions =
366 CSLSetNameValue( papszRasterizeOptions, "ALL_TOUCHED", "TRUE" );
367
368 int anXYOff[2] = { nXOff, nYOff };
369
370 CPLErr eErr =
371 GDALRasterizeGeometries( hMemDS, 1, &nTargetBand,
372 1, &hPolygon,
373 CutlineTransformer, anXYOff,
374 &dfBurnValue, papszRasterizeOptions,
375 nullptr, nullptr );
376
377 CSLDestroy( papszRasterizeOptions );
378
379 // Close and ensure data flushed to underlying array.
380 GDALClose( hMemDS );
381
382 /* -------------------------------------------------------------------- */
383 /* In the case with no blend distance, we just apply this as a */
384 /* mask, zeroing out everything outside the polygon. */
385 /* -------------------------------------------------------------------- */
386 if( psWO->dfCutlineBlendDist == 0.0 )
387 {
388 for( int i = nXSize * nYSize - 1; i >= 0; i-- )
389 {
390 if( pabyPolyMask[i] == 0 )
391 static_cast<float *>(pValidityMask)[i] = 0.0;
392 }
393 }
394 else
395 {
396 eErr = BlendMaskGenerator( nXOff, nYOff, nXSize, nYSize,
397 pabyPolyMask,
398 static_cast<float *>(pValidityMask),
399 hPolygon, psWO->dfCutlineBlendDist );
400 }
401
402 /* -------------------------------------------------------------------- */
403 /* Clean up. */
404 /* -------------------------------------------------------------------- */
405 CPLFree( pabyPolyMask );
406
407 return eErr;
408 }
409