1 : /******************************************************************************
2 : * $Id: gdalcutline.cpp 17041 2009-05-17 22:35:15Z warmerdam $
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 : *
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 "gdalwarper.h"
31 : #include "gdal_alg.h"
32 : #include "ogr_api.h"
33 : #include "ogr_geos.h"
34 : #include "ogr_geometry.h"
35 : #include "cpl_string.h"
36 :
37 : CPL_CVSID("$Id: gdalcutline.cpp 17041 2009-05-17 22:35:15Z warmerdam $");
38 :
39 : /************************************************************************/
40 : /* BlendMaskGenerator() */
41 : /************************************************************************/
42 :
43 : static CPLErr
44 1 : BlendMaskGenerator( int nXOff, int nYOff, int nXSize, int nYSize,
45 : GByte *pabyPolyMask, float *pafValidityMask,
46 : OGRGeometryH hPolygon, double dfBlendDist )
47 :
48 : {
49 : #ifndef HAVE_GEOS
50 : CPLError( CE_Failure, CPLE_AppDefined,
51 : "Blend distance support not available without the GEOS library.");
52 : return CE_Failure;
53 :
54 : #else /* HAVE_GEOS */
55 :
56 : /* -------------------------------------------------------------------- */
57 : /* Convert the polygon into a collection of lines so that we */
58 : /* measure distance from the edge even on the inside. */
59 : /* -------------------------------------------------------------------- */
60 : OGRGeometry *poLines
61 : = OGRGeometryFactory::forceToMultiLineString(
62 1 : ((OGRGeometry *) hPolygon)->clone() );
63 :
64 : /* -------------------------------------------------------------------- */
65 : /* Prepare a clipping polygon a bit bigger than the area of */
66 : /* interest in the hopes of simplifying the cutline down to */
67 : /* stuff that will be relavent for this area of interest. */
68 : /* -------------------------------------------------------------------- */
69 1 : CPLString osClipRectWKT;
70 :
71 : osClipRectWKT.Printf( "POLYGON((%g %g,%g %g,%g %g,%g %g,%g %g))",
72 : nXOff - (dfBlendDist+1),
73 : nYOff - (dfBlendDist+1),
74 : nXOff + nXSize + (dfBlendDist+1),
75 : nYOff - (dfBlendDist+1),
76 : nXOff + nXSize + (dfBlendDist+1),
77 : nYOff + nYSize + (dfBlendDist+1),
78 : nXOff - (dfBlendDist+1),
79 : nYOff + nYSize + (dfBlendDist+1),
80 : nXOff - (dfBlendDist+1),
81 1 : nYOff - (dfBlendDist+1) );
82 :
83 1 : OGRPolygon *poClipRect = NULL;
84 1 : char *pszWKT = (char *) osClipRectWKT.c_str();
85 :
86 : OGRGeometryFactory::createFromWkt( &pszWKT, NULL,
87 1 : (OGRGeometry**) (&poClipRect) );
88 :
89 1 : if( poClipRect )
90 : {
91 : OGRGeometry *poClippedLines =
92 1 : poLines->Intersection( poClipRect );
93 1 : delete poLines;
94 1 : poLines = poClippedLines;
95 1 : delete poClipRect;
96 : }
97 :
98 : /* -------------------------------------------------------------------- */
99 : /* Convert our polygon into GEOS format, and compute an */
100 : /* envelope to accelerate later distance operations. */
101 : /* -------------------------------------------------------------------- */
102 1 : OGREnvelope sEnvelope;
103 : int iXMin, iYMin, iXMax, iYMax;
104 : GEOSGeom poGEOSPoly;
105 :
106 1 : poGEOSPoly = poLines->exportToGEOS();
107 1 : OGR_G_GetEnvelope( hPolygon, &sEnvelope );
108 :
109 1 : delete poLines;
110 :
111 1 : if( sEnvelope.MinY - dfBlendDist > nYOff+nYSize
112 : || sEnvelope.MaxY + dfBlendDist < nYOff
113 : || sEnvelope.MinX - dfBlendDist > nXOff+nXSize
114 : || sEnvelope.MaxX + dfBlendDist < nXOff )
115 0 : return CE_None;
116 :
117 1 : iXMin = MAX(0,(int) floor(sEnvelope.MinX - dfBlendDist - nXOff));
118 1 : iXMax = MIN(nXSize, (int) ceil(sEnvelope.MaxX + dfBlendDist - nXOff));
119 1 : iYMin = MAX(0,(int) floor(sEnvelope.MinY - dfBlendDist - nYOff));
120 1 : iYMax = MIN(nYSize, (int) ceil(sEnvelope.MaxY + dfBlendDist - nYOff));
121 :
122 : /* -------------------------------------------------------------------- */
123 : /* Loop over potential area within blend line distance, */
124 : /* processing each pixel. */
125 : /* -------------------------------------------------------------------- */
126 : int iY, iX;
127 : double dfLastDist;
128 :
129 101 : for( iY = 0; iY < nYSize; iY++ )
130 : {
131 100 : dfLastDist = 0.0;
132 :
133 100 : for( iX = 0; iX < nXSize; iX++ )
134 : {
135 10000 : if( iX < iXMin || iX >= iXMax
136 : || iY < iYMin || iY > iYMax
137 : || dfLastDist > dfBlendDist + 1.5 )
138 : {
139 7944 : if( pabyPolyMask[iX + iY * nXSize] == 0 )
140 7784 : pafValidityMask[iX + iY * nXSize] = 0.0;
141 :
142 7944 : dfLastDist -= 1.0;
143 7944 : continue;
144 : }
145 :
146 : double dfDist, dfRatio;
147 2056 : CPLString osPointWKT;
148 : GEOSGeom poGEOSPoint;
149 :
150 2056 : osPointWKT.Printf( "POINT(%d.5 %d.5)", iX + nXOff, iY + nYOff );
151 2056 : poGEOSPoint = GEOSGeomFromWKT( osPointWKT );
152 :
153 2056 : GEOSDistance( poGEOSPoly, poGEOSPoint, &dfDist );
154 2056 : GEOSGeom_destroy( poGEOSPoint );
155 :
156 2056 : dfLastDist = dfDist;
157 :
158 2056 : if( dfDist > dfBlendDist )
159 : {
160 584 : if( pabyPolyMask[iX + iY * nXSize] == 0 )
161 366 : pafValidityMask[iX + iY * nXSize] = 0.0;
162 :
163 584 : continue;
164 : }
165 :
166 1472 : if( pabyPolyMask[iX + iY * nXSize] == 0 )
167 : {
168 : /* outside */
169 850 : dfRatio = 0.5 - (dfDist / dfBlendDist) * 0.5;
170 : }
171 : else
172 : {
173 : /* inside */
174 622 : dfRatio = 0.5 + (dfDist / dfBlendDist) * 0.5;
175 : }
176 :
177 1472 : pafValidityMask[iX + iY * nXSize] *= dfRatio;
178 : }
179 : }
180 :
181 : /* -------------------------------------------------------------------- */
182 : /* Cleanup */
183 : /* -------------------------------------------------------------------- */
184 1 : GEOSGeom_destroy( poGEOSPoly );
185 :
186 1 : return CE_None;
187 :
188 : #endif /* HAVE_GEOS */
189 : }
190 :
191 : /************************************************************************/
192 : /* CutlineTransformer() */
193 : /* */
194 : /* A simple transformer for the cutline that just offsets */
195 : /* relative to the current chunk. */
196 : /************************************************************************/
197 :
198 6 : static int CutlineTransformer( void *pTransformArg, int bDstToSrc,
199 : int nPointCount,
200 : double *x, double *y, double *z,
201 : int *panSuccess )
202 :
203 : {
204 6 : int nXOff = ((int *) pTransformArg)[0];
205 6 : int nYOff = ((int *) pTransformArg)[1];
206 : int i;
207 :
208 6 : if( bDstToSrc )
209 : {
210 0 : nXOff *= -1;
211 0 : nYOff *= -1;
212 : }
213 :
214 43 : for( i = 0; i < nPointCount; i++ )
215 : {
216 37 : x[i] -= nXOff;
217 37 : y[i] -= nYOff;
218 : }
219 :
220 6 : return TRUE;
221 : }
222 :
223 :
224 : /************************************************************************/
225 : /* GDALWarpCutlineMasker() */
226 : /* */
227 : /* This function will generate a source mask based on a */
228 : /* provided cutline, and optional blend distance. */
229 : /************************************************************************/
230 :
231 : CPLErr
232 6 : GDALWarpCutlineMasker( void *pMaskFuncArg, int nBandCount, GDALDataType eType,
233 : int nXOff, int nYOff, int nXSize, int nYSize,
234 : GByte ** /*ppImageData */,
235 : int bMaskIsFloat, void *pValidityMask )
236 :
237 : {
238 6 : GDALWarpOptions *psWO = (GDALWarpOptions *) pMaskFuncArg;
239 6 : float *pafMask = (float *) pValidityMask;
240 : CPLErr eErr;
241 : GDALDriverH hMemDriver;
242 :
243 6 : if( nXSize < 1 || nYSize < 1 )
244 0 : return CE_None;
245 :
246 : /* -------------------------------------------------------------------- */
247 : /* Do some minimal checking. */
248 : /* -------------------------------------------------------------------- */
249 6 : if( !bMaskIsFloat )
250 : {
251 : CPLAssert( FALSE );
252 0 : return CE_Failure;
253 : }
254 :
255 6 : if( psWO == NULL || psWO->hCutline == NULL )
256 : {
257 : CPLAssert( FALSE );
258 0 : return CE_Failure;
259 : }
260 :
261 6 : hMemDriver = GDALGetDriverByName("MEM");
262 6 : if (hMemDriver == NULL)
263 : {
264 0 : CPLError(CE_Failure, CPLE_AppDefined, "GDALWarpCutlineMasker needs MEM driver");
265 0 : return CE_Failure;
266 : }
267 :
268 : /* -------------------------------------------------------------------- */
269 : /* Check the polygon. */
270 : /* -------------------------------------------------------------------- */
271 6 : OGRGeometryH hPolygon = (OGRGeometryH) psWO->hCutline;
272 6 : OGREnvelope sEnvelope;
273 :
274 6 : if( wkbFlatten(OGR_G_GetGeometryType(hPolygon)) != wkbPolygon
275 : && wkbFlatten(OGR_G_GetGeometryType(hPolygon)) != wkbMultiPolygon )
276 : {
277 : CPLAssert( FALSE );
278 0 : return CE_Failure;
279 : }
280 :
281 6 : OGR_G_GetEnvelope( hPolygon, &sEnvelope );
282 6 : if( sEnvelope.MaxX + psWO->dfCutlineBlendDist < nXOff
283 : || sEnvelope.MinX - psWO->dfCutlineBlendDist > nXOff + nXSize
284 : || sEnvelope.MaxY + psWO->dfCutlineBlendDist < nYOff
285 : || sEnvelope.MinY - psWO->dfCutlineBlendDist > nYOff + nYSize )
286 : {
287 : // We are far from the blend line - everything is masked to zero.
288 : // It would be nice to realize no work is required for this whole
289 : // chunk!
290 0 : memset( pafMask, 0, sizeof(float) * nXSize * nYSize );
291 0 : return CE_None;
292 : }
293 :
294 : /* -------------------------------------------------------------------- */
295 : /* Create a byte buffer into which we can burn the */
296 : /* mask polygon and wrap it up as a memory dataset. */
297 : /* -------------------------------------------------------------------- */
298 6 : GByte *pabyPolyMask = (GByte *) CPLCalloc( nXSize, nYSize );
299 : GDALDatasetH hMemDS;
300 6 : double adfGeoTransform[6] = { 0.0, 1.0, 0.0, 0.0, 0.0, 1.0 };
301 :
302 : char szDataPointer[100];
303 6 : char *apszOptions[] = { szDataPointer, NULL };
304 :
305 6 : memset( szDataPointer, 0, sizeof(szDataPointer) );
306 6 : sprintf( szDataPointer, "DATAPOINTER=" );
307 : CPLPrintPointer( szDataPointer+strlen(szDataPointer),
308 : pabyPolyMask,
309 6 : sizeof(szDataPointer) - strlen(szDataPointer) );
310 :
311 : hMemDS = GDALCreate( hMemDriver, "warp_temp",
312 6 : nXSize, nYSize, 0, GDT_Byte, NULL );
313 6 : GDALAddBand( hMemDS, GDT_Byte, apszOptions );
314 6 : GDALSetGeoTransform( hMemDS, adfGeoTransform );
315 :
316 : /* -------------------------------------------------------------------- */
317 : /* Burn the polygon into the mask with 1.0 values. */
318 : /* -------------------------------------------------------------------- */
319 6 : int nTargetBand = 1;
320 6 : double dfBurnValue = 255.0;
321 : int anXYOff[2];
322 6 : char **papszRasterizeOptions = NULL;
323 :
324 :
325 6 : if( CSLFetchBoolean( psWO->papszWarpOptions, "CUTLINE_ALL_TOUCHED", FALSE ))
326 : papszRasterizeOptions =
327 1 : CSLSetNameValue( papszRasterizeOptions, "ALL_TOUCHED", "TRUE" );
328 :
329 6 : anXYOff[0] = nXOff;
330 6 : anXYOff[1] = nYOff;
331 :
332 : eErr =
333 : GDALRasterizeGeometries( hMemDS, 1, &nTargetBand,
334 : 1, &hPolygon,
335 : CutlineTransformer, anXYOff,
336 : &dfBurnValue, papszRasterizeOptions,
337 6 : NULL, NULL );
338 :
339 6 : CSLDestroy( papszRasterizeOptions );
340 :
341 : // Close and ensure data flushed to underlying array.
342 6 : GDALClose( hMemDS );
343 :
344 : /* -------------------------------------------------------------------- */
345 : /* In the case with no blend distance, we just apply this as a */
346 : /* mask, zeroing out everything outside the polygon. */
347 : /* -------------------------------------------------------------------- */
348 6 : if( psWO->dfCutlineBlendDist == 0.0 )
349 : {
350 : int i;
351 :
352 50005 : for( i = nXSize * nYSize - 1; i >= 0; i-- )
353 : {
354 50000 : if( pabyPolyMask[i] == 0 )
355 41931 : ((float *) pValidityMask)[i] = 0.0;
356 : }
357 : }
358 : else
359 : {
360 : eErr = BlendMaskGenerator( nXOff, nYOff, nXSize, nYSize,
361 : pabyPolyMask, (float *) pValidityMask,
362 1 : hPolygon, psWO->dfCutlineBlendDist );
363 : }
364 :
365 : /* -------------------------------------------------------------------- */
366 : /* Clean up. */
367 : /* -------------------------------------------------------------------- */
368 6 : CPLFree( pabyPolyMask );
369 :
370 6 : return eErr;
371 : }
372 :
|
