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