1 : /******************************************************************************
2 : * $Id: gdalwarpoperation.cpp 18277 2009-12-12 17:33:37Z rouault $
3 : *
4 : * Project: High Performance Image Reprojector
5 : * Purpose: Implementation of the GDALWarpOperation class.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2003, 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 "cpl_string.h"
32 : #include "cpl_multiproc.h"
33 : #include "ogr_api.h"
34 :
35 : CPL_CVSID("$Id: gdalwarpoperation.cpp 18277 2009-12-12 17:33:37Z rouault $");
36 :
37 : /* Defined in gdalwarpkernel.cpp */
38 : int GWKGetFilterRadius(GDALResampleAlg eResampleAlg);
39 :
40 :
41 : /************************************************************************/
42 : /* ==================================================================== */
43 : /* GDALWarpOperation */
44 : /* ==================================================================== */
45 : /************************************************************************/
46 :
47 : /**
48 : * \class GDALWarpOperation "gdalwarper.h"
49 : *
50 : * High level image warping class.
51 :
52 : <h2>Warper Design</h2>
53 :
54 : The overall GDAL high performance image warper is split into a few components.
55 :
56 : - The transformation between input and output file coordinates is handled
57 : via GDALTransformerFunc() implementations such as the one returned by
58 : GDALCreateGenImgProjTransformer(). The transformers are ultimately responsible
59 : for translating pixel/line locations on the destination image to pixel/line
60 : locations on the source image.
61 :
62 : - In order to handle images too large to hold in RAM, the warper needs to
63 : segment large images. This is the responsibility of the GDALWarpOperation
64 : class. The GDALWarpOperation::ChunkAndWarpImage() invokes
65 : GDALWarpOperation::WarpRegion() on chunks of output and input image that
66 : are small enough to hold in the amount of memory allowed by the application.
67 : This process is described in greater detail in the <b>Image Chunking</b>
68 : section.
69 :
70 : - The GDALWarpOperation::WarpRegion() function creates and loads an output
71 : image buffer, and then calls WarpRegionToBuffer().
72 :
73 : - GDALWarpOperation::WarpRegionToBuffer() is responsible for loading the
74 : source imagery corresponding to a particular output region, and generating
75 : masks and density masks from the source and destination imagery using
76 : the generator functions found in the GDALWarpOptions structure. Binds this
77 : all into an instance of GDALWarpKernel on which the
78 : GDALWarpKernel::PerformWarp() method is called.
79 :
80 : - GDALWarpKernel does the actual image warping, but is given an input image
81 : and an output image to operate on. The GDALWarpKernel does no IO, and in
82 : fact knows nothing about GDAL. It invokes the transformation function to
83 : get sample locations, builds output values based on the resampling algorithm
84 : in use. It also takes any validity and density masks into account during
85 : this operation.
86 :
87 : <h3>Chunk Size Selection</h3>
88 :
89 : The GDALWarpOptions ChunkAndWarpImage() method is responsible for invoking
90 : the WarpRegion() method on appropriate sized output chunks such that the
91 : memory required for the output image buffer, input image buffer and any
92 : required density and validity buffers is less than or equal to the application
93 : defined maximum memory available for use.
94 :
95 : It checks the memory requrired by walking the edges of the output region,
96 : transforming the locations back into source pixel/line coordinates and
97 : establishing a bounding rectangle of source imagery that would be required
98 : for the output area. This is actually accomplished by the private
99 : GDALWarpOperation::ComputeSourceWindow() method.
100 :
101 : Then memory requirements are used by totaling the memory required for all
102 : output bands, input bands, validity masks and density masks. If this is
103 : greater than the GDALWarpOptions::dfWarpMemoryLimit then the destination
104 : region is divided in two (splitting the longest dimension), and
105 : ChunkAndWarpImage() recursively invoked on each destination subregion.
106 :
107 : <h3>Validity and Density Masks Generation</h3>
108 :
109 : Fill in ways in which the validity and density masks may be generated here.
110 : Note that detailed semantics of the masks should be found in
111 : GDALWarpKernel.
112 :
113 : */
114 :
115 : /************************************************************************/
116 : /* GDALWarpOperation() */
117 : /************************************************************************/
118 :
119 304 : GDALWarpOperation::GDALWarpOperation()
120 :
121 : {
122 304 : psOptions = NULL;
123 :
124 304 : dfProgressBase = 0.0;
125 304 : dfProgressScale = 1.0;
126 :
127 304 : hThread1Mutex = NULL;
128 304 : hThread2Mutex = NULL;
129 304 : hIOMutex = NULL;
130 304 : hWarpMutex = NULL;
131 :
132 304 : nChunkListCount = 0;
133 304 : nChunkListMax = 0;
134 304 : panChunkList = NULL;
135 :
136 304 : bReportTimings = FALSE;
137 304 : nLastTimeReported = 0;
138 304 : }
139 :
140 : /************************************************************************/
141 : /* ~GDALWarpOperation() */
142 : /************************************************************************/
143 :
144 343 : GDALWarpOperation::~GDALWarpOperation()
145 :
146 : {
147 304 : WipeOptions();
148 :
149 304 : if( hThread1Mutex != NULL )
150 : {
151 1 : CPLDestroyMutex( hThread1Mutex );
152 1 : CPLDestroyMutex( hThread2Mutex );
153 1 : CPLDestroyMutex( hIOMutex );
154 1 : CPLDestroyMutex( hWarpMutex );
155 : }
156 :
157 304 : WipeChunkList();
158 343 : }
159 :
160 : /************************************************************************/
161 : /* GetOptions() */
162 : /************************************************************************/
163 :
164 118 : const GDALWarpOptions *GDALWarpOperation::GetOptions()
165 :
166 : {
167 118 : return psOptions;
168 : }
169 :
170 : /************************************************************************/
171 : /* WipeOptions() */
172 : /************************************************************************/
173 :
174 304 : void GDALWarpOperation::WipeOptions()
175 :
176 : {
177 304 : if( psOptions != NULL )
178 : {
179 304 : GDALDestroyWarpOptions( psOptions );
180 304 : psOptions = NULL;
181 : }
182 304 : }
183 :
184 : /************************************************************************/
185 : /* ValidateOptions() */
186 : /************************************************************************/
187 :
188 304 : int GDALWarpOperation::ValidateOptions()
189 :
190 : {
191 304 : if( psOptions == NULL )
192 : {
193 : CPLError( CE_Failure, CPLE_IllegalArg,
194 : "GDALWarpOptions.Validate()\n"
195 0 : " no options currently initialized." );
196 0 : return FALSE;
197 : }
198 :
199 304 : if( psOptions->dfWarpMemoryLimit < 100000.0 )
200 : {
201 : CPLError( CE_Failure, CPLE_IllegalArg,
202 : "GDALWarpOptions.Validate()\n"
203 : " dfWarpMemoryLimit=%g is unreasonably small.",
204 0 : psOptions->dfWarpMemoryLimit );
205 0 : return FALSE;
206 : }
207 :
208 304 : if( psOptions->eResampleAlg != GRA_NearestNeighbour
209 : && psOptions->eResampleAlg != GRA_Bilinear
210 : && psOptions->eResampleAlg != GRA_Cubic
211 : && psOptions->eResampleAlg != GRA_CubicSpline
212 : && psOptions->eResampleAlg != GRA_Lanczos )
213 : {
214 : CPLError( CE_Failure, CPLE_IllegalArg,
215 : "GDALWarpOptions.Validate()\n"
216 : " eResampleArg=%d is not a supported value.",
217 0 : psOptions->eResampleAlg );
218 0 : return FALSE;
219 : }
220 :
221 304 : if( (int) psOptions->eWorkingDataType < 1
222 : && (int) psOptions->eWorkingDataType >= GDT_TypeCount )
223 : {
224 : CPLError( CE_Failure, CPLE_IllegalArg,
225 : "GDALWarpOptions.Validate()\n"
226 : " eWorkingDataType=%d is not a supported value.",
227 0 : psOptions->eWorkingDataType );
228 0 : return FALSE;
229 : }
230 :
231 304 : if( psOptions->hSrcDS == NULL )
232 : {
233 : CPLError( CE_Failure, CPLE_IllegalArg,
234 : "GDALWarpOptions.Validate()\n"
235 0 : " hSrcDS is not set." );
236 0 : return FALSE;
237 : }
238 :
239 304 : if( psOptions->nBandCount == 0 )
240 : {
241 : CPLError( CE_Failure, CPLE_IllegalArg,
242 : "GDALWarpOptions.Validate()\n"
243 0 : " nBandCount=0, no bands configured!" );
244 0 : return FALSE;
245 : }
246 :
247 304 : if( psOptions->panSrcBands == NULL )
248 : {
249 : CPLError( CE_Failure, CPLE_IllegalArg,
250 : "GDALWarpOptions.Validate()\n"
251 0 : " panSrcBands is NULL." );
252 0 : return FALSE;
253 : }
254 :
255 304 : if( psOptions->hDstDS != NULL && psOptions->panDstBands == NULL )
256 : {
257 : CPLError( CE_Failure, CPLE_IllegalArg,
258 : "GDALWarpOptions.Validate()\n"
259 0 : " panDstBands is NULL." );
260 0 : return FALSE;
261 : }
262 :
263 642 : for( int iBand = 0; iBand < psOptions->nBandCount; iBand++ )
264 : {
265 676 : if( psOptions->panSrcBands[iBand] < 1
266 338 : || psOptions->panSrcBands[iBand]
267 : > GDALGetRasterCount( psOptions->hSrcDS ) )
268 : {
269 : CPLError( CE_Failure, CPLE_IllegalArg,
270 : "panSrcBands[%d] = %d ... out of range for dataset.",
271 0 : iBand, psOptions->panSrcBands[iBand] );
272 0 : return FALSE;
273 : }
274 1014 : if( psOptions->hDstDS != NULL
275 338 : && (psOptions->panDstBands[iBand] < 1
276 338 : || psOptions->panDstBands[iBand]
277 : > GDALGetRasterCount( psOptions->hDstDS ) ) )
278 : {
279 : CPLError( CE_Failure, CPLE_IllegalArg,
280 : "panDstBands[%d] = %d ... out of range for dataset.",
281 0 : iBand, psOptions->panDstBands[iBand] );
282 0 : return FALSE;
283 : }
284 :
285 338 : if( psOptions->hDstDS != NULL
286 : && GDALGetRasterAccess(
287 : GDALGetRasterBand(psOptions->hDstDS,
288 : psOptions->panDstBands[iBand]))
289 : == GA_ReadOnly )
290 : {
291 : CPLError( CE_Failure, CPLE_IllegalArg,
292 : "Destination band %d appears to be read-only.",
293 0 : psOptions->panDstBands[iBand] );
294 0 : return FALSE;
295 : }
296 : }
297 :
298 304 : if( psOptions->nBandCount == 0 )
299 : {
300 : CPLError( CE_Failure, CPLE_IllegalArg,
301 : "GDALWarpOptions.Validate()\n"
302 0 : " nBandCount=0, no bands configured!" );
303 0 : return FALSE;
304 : }
305 :
306 304 : if( psOptions->padfSrcNoDataReal != NULL
307 : && psOptions->padfSrcNoDataImag == NULL )
308 : {
309 : CPLError( CE_Failure, CPLE_IllegalArg,
310 : "GDALWarpOptions.Validate()\n"
311 0 : " padfSrcNoDataReal set, but padfSrcNoDataImag not set." );
312 0 : return FALSE;
313 : }
314 :
315 304 : if( psOptions->pfnProgress == NULL )
316 : {
317 : CPLError( CE_Failure, CPLE_IllegalArg,
318 : "GDALWarpOptions.Validate()\n"
319 0 : " pfnProgress is NULL." );
320 0 : return FALSE;
321 : }
322 :
323 304 : if( psOptions->pfnTransformer == NULL )
324 : {
325 : CPLError( CE_Failure, CPLE_IllegalArg,
326 : "GDALWarpOptions.Validate()\n"
327 0 : " pfnTransformer is NULL." );
328 0 : return FALSE;
329 : }
330 :
331 304 : if( CSLFetchNameValue( psOptions->papszWarpOptions,
332 : "SAMPLE_STEPS" ) != NULL )
333 : {
334 0 : if( atoi(CSLFetchNameValue( psOptions->papszWarpOptions,
335 : "SAMPLE_STEPS" )) < 2 )
336 : {
337 : CPLError( CE_Failure, CPLE_IllegalArg,
338 : "GDALWarpOptions.Validate()\n"
339 0 : " SAMPLE_STEPS warp option has illegal value." );
340 0 : return FALSE;
341 : }
342 : }
343 :
344 304 : if( psOptions->nSrcAlphaBand > 0
345 : && psOptions->pfnSrcDensityMaskFunc != NULL )
346 : {
347 : CPLError( CE_Failure, CPLE_IllegalArg,
348 : "GDALWarpOptions.Validate()\n"
349 0 : " pfnSrcDensityMaskFunc provided as well as a SrcAlphaBand." );
350 0 : return FALSE;
351 : }
352 :
353 304 : if( psOptions->nDstAlphaBand > 0
354 : && psOptions->pfnDstDensityMaskFunc != NULL )
355 : {
356 : CPLError( CE_Failure, CPLE_IllegalArg,
357 : "GDALWarpOptions.Validate()\n"
358 0 : " pfnDstDensityMaskFunc provided as well as a DstAlphaBand." );
359 0 : return FALSE;
360 : }
361 :
362 304 : return TRUE;
363 : }
364 :
365 : /************************************************************************/
366 : /* Initialize() */
367 : /************************************************************************/
368 :
369 : /**
370 : * \fn CPLErr GDALWarpOperation::Initialize( const GDALWarpOptions * );
371 : *
372 : * This method initializes the GDALWarpOperation's concept of the warp
373 : * options in effect. It creates an internal copy of the GDALWarpOptions
374 : * structure and defaults a variety of additional fields in the internal
375 : * copy if not set in the provides warp options.
376 : *
377 : * Defaulting operations include:
378 : * - If the nBandCount is 0, it will be set to the number of bands in the
379 : * source image (which must match the output image) and the panSrcBands
380 : * and panDstBands will be populated.
381 : *
382 : * @param psNewOptions input set of warp options. These are copied and may
383 : * be destroyed after this call by the application.
384 : *
385 : * @return CE_None on success or CE_Failure if an error occurs.
386 : */
387 :
388 304 : CPLErr GDALWarpOperation::Initialize( const GDALWarpOptions *psNewOptions )
389 :
390 : {
391 304 : CPLErr eErr = CE_None;
392 :
393 : /* -------------------------------------------------------------------- */
394 : /* Copy the passed in options. */
395 : /* -------------------------------------------------------------------- */
396 304 : if( psOptions != NULL )
397 0 : WipeOptions();
398 :
399 304 : psOptions = GDALCloneWarpOptions( psNewOptions );
400 :
401 : /* -------------------------------------------------------------------- */
402 : /* Default band mapping if missing. */
403 : /* -------------------------------------------------------------------- */
404 304 : if( psOptions->nBandCount == 0
405 : && psOptions->hSrcDS != NULL
406 : && psOptions->hDstDS != NULL
407 : && GDALGetRasterCount( psOptions->hSrcDS )
408 : == GDALGetRasterCount( psOptions->hDstDS ) )
409 : {
410 : int i;
411 :
412 0 : psOptions->nBandCount = GDALGetRasterCount( psOptions->hSrcDS );
413 :
414 : psOptions->panSrcBands = (int *)
415 0 : CPLMalloc(sizeof(int) * psOptions->nBandCount );
416 : psOptions->panDstBands = (int *)
417 0 : CPLMalloc(sizeof(int) * psOptions->nBandCount );
418 :
419 0 : for( i = 0; i < psOptions->nBandCount; i++ )
420 : {
421 0 : psOptions->panSrcBands[i] = i+1;
422 0 : psOptions->panDstBands[i] = i+1;
423 : }
424 : }
425 :
426 : /* -------------------------------------------------------------------- */
427 : /* If no working data type was provided, set one now. */
428 : /* */
429 : /* Default to the highest resolution output band. But if the */
430 : /* input band is higher resolution and has a nodata value "out */
431 : /* of band" with the output type we may need to use the higher */
432 : /* resolution input type to ensure we can identify nodata values. */
433 : /* -------------------------------------------------------------------- */
434 304 : if( psOptions->eWorkingDataType == GDT_Unknown
435 : && psOptions->hSrcDS != NULL
436 : && psOptions->hDstDS != NULL
437 : && psOptions->nBandCount >= 1 )
438 : {
439 : int iBand;
440 268 : psOptions->eWorkingDataType = GDT_Byte;
441 :
442 550 : for( iBand = 0; iBand < psOptions->nBandCount; iBand++ )
443 : {
444 : GDALRasterBandH hDstBand = GDALGetRasterBand(
445 282 : psOptions->hDstDS, psOptions->panDstBands[iBand] );
446 : GDALRasterBandH hSrcBand = GDALGetRasterBand(
447 282 : psOptions->hSrcDS, psOptions->panSrcBands[iBand] );
448 :
449 282 : if( hDstBand != NULL )
450 : psOptions->eWorkingDataType =
451 : GDALDataTypeUnion( psOptions->eWorkingDataType,
452 282 : GDALGetRasterDataType( hDstBand ) );
453 :
454 282 : if( hSrcBand != NULL
455 : && psOptions->padfSrcNoDataReal != NULL )
456 : {
457 5 : int bMergeSource = FALSE;
458 :
459 10 : if( psOptions->padfSrcNoDataImag != NULL
460 5 : && psOptions->padfSrcNoDataImag[iBand] != 0.0
461 : && !GDALDataTypeIsComplex( psOptions->eWorkingDataType ) )
462 0 : bMergeSource = TRUE;
463 5 : else if( psOptions->padfSrcNoDataReal[iBand] < 0.0
464 : && (psOptions->eWorkingDataType == GDT_Byte
465 : || psOptions->eWorkingDataType == GDT_UInt16
466 : || psOptions->eWorkingDataType == GDT_UInt32) )
467 0 : bMergeSource = TRUE;
468 5 : else if( psOptions->padfSrcNoDataReal[iBand] < -32768.0
469 : && psOptions->eWorkingDataType == GDT_Int16 )
470 0 : bMergeSource = TRUE;
471 5 : else if( psOptions->padfSrcNoDataReal[iBand] < -2147483648.0
472 : && psOptions->eWorkingDataType == GDT_Int32 )
473 0 : bMergeSource = TRUE;
474 5 : else if( psOptions->padfSrcNoDataReal[iBand] > 256
475 : && psOptions->eWorkingDataType == GDT_Byte )
476 0 : bMergeSource = TRUE;
477 5 : else if( psOptions->padfSrcNoDataReal[iBand] > 32767
478 : && psOptions->eWorkingDataType == GDT_Int16 )
479 0 : bMergeSource = TRUE;
480 5 : else if( psOptions->padfSrcNoDataReal[iBand] > 65535
481 : && psOptions->eWorkingDataType == GDT_UInt16 )
482 0 : bMergeSource = TRUE;
483 5 : else if( psOptions->padfSrcNoDataReal[iBand] > 2147483648.0
484 : && psOptions->eWorkingDataType == GDT_Int32 )
485 0 : bMergeSource = TRUE;
486 5 : else if( psOptions->padfSrcNoDataReal[iBand] > 4294967295.0
487 : && psOptions->eWorkingDataType == GDT_UInt32 )
488 0 : bMergeSource = TRUE;
489 :
490 5 : if( bMergeSource )
491 : psOptions->eWorkingDataType =
492 : GDALDataTypeUnion( psOptions->eWorkingDataType,
493 0 : GDALGetRasterDataType( hSrcBand ) );
494 : }
495 : }
496 : }
497 :
498 : /* -------------------------------------------------------------------- */
499 : /* Default memory available. */
500 : /* */
501 : /* For now we default to 64MB of RAM, but eventually we should */
502 : /* try various schemes to query physical RAM. This can */
503 : /* certainly be done on Win32 and Linux. */
504 : /* -------------------------------------------------------------------- */
505 304 : if( psOptions->dfWarpMemoryLimit == 0.0 )
506 : {
507 267 : psOptions->dfWarpMemoryLimit = 64.0 * 1024*1024;
508 : }
509 :
510 : /* -------------------------------------------------------------------- */
511 : /* Are we doing timings? */
512 : /* -------------------------------------------------------------------- */
513 : bReportTimings = CSLFetchBoolean( psOptions->papszWarpOptions,
514 304 : "REPORT_TIMINGS", FALSE );
515 :
516 : /* -------------------------------------------------------------------- */
517 : /* Support creating cutline from text warpoption. */
518 : /* -------------------------------------------------------------------- */
519 : const char *pszCutlineWKT =
520 304 : CSLFetchNameValue( psOptions->papszWarpOptions, "CUTLINE" );
521 :
522 304 : if( pszCutlineWKT )
523 : {
524 3 : if( OGR_G_CreateFromWkt( (char **) &pszCutlineWKT, NULL,
525 : (OGRGeometryH *) &(psOptions->hCutline) )
526 : != OGRERR_NONE )
527 : {
528 0 : eErr = CE_Failure;
529 : CPLError( CE_Failure, CPLE_AppDefined,
530 0 : "Failed to parse CUTLINE geometry wkt." );
531 : }
532 : else
533 : {
534 : const char *pszBD = CSLFetchNameValue( psOptions->papszWarpOptions,
535 3 : "CUTLINE_BLEND_DIST" );
536 3 : if( pszBD )
537 0 : psOptions->dfCutlineBlendDist = atof(pszBD);
538 : }
539 : }
540 :
541 : /* -------------------------------------------------------------------- */
542 : /* If the options don't validate, then wipe them. */
543 : /* -------------------------------------------------------------------- */
544 304 : if( !ValidateOptions() )
545 0 : eErr = CE_Failure;
546 :
547 304 : if( eErr != CE_None )
548 0 : WipeOptions();
549 :
550 304 : return eErr;
551 : }
552 :
553 : /************************************************************************/
554 : /* GDALCreateWarpOperation() */
555 : /************************************************************************/
556 :
557 : /**
558 : * @see GDALWarpOperation::Initialize()
559 : */
560 :
561 0 : GDALWarpOperationH GDALCreateWarpOperation(
562 : const GDALWarpOptions *psNewOptions )
563 : {
564 : GDALWarpOperation *poOperation;
565 :
566 0 : poOperation = new GDALWarpOperation;
567 0 : if ( poOperation->Initialize( psNewOptions ) != CE_None )
568 : {
569 0 : delete poOperation;
570 0 : return NULL;
571 : }
572 :
573 0 : return (GDALWarpOperationH)poOperation;
574 : }
575 :
576 : /************************************************************************/
577 : /* GDALDestroyWarpOperation() */
578 : /************************************************************************/
579 :
580 : /**
581 : * @see GDALWarpOperation::~GDALWarpOperation()
582 : */
583 :
584 0 : void GDALDestroyWarpOperation( GDALWarpOperationH hOperation )
585 : {
586 0 : if ( hOperation )
587 0 : delete static_cast<GDALWarpOperation *>(hOperation);
588 0 : }
589 :
590 : /************************************************************************/
591 : /* ChunkAndWarpImage() */
592 : /************************************************************************/
593 :
594 : /**
595 : * \fn CPLErr GDALWarpOperation::ChunkAndWarpImage(
596 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize );
597 : *
598 : * This method does a complete warp of the source image to the destination
599 : * image for the indicated region with the current warp options in effect.
600 : * Progress is reported to the installed progress monitor, if any.
601 : *
602 : * This function will subdivide the region and recursively call itself
603 : * until the total memory required to process a region chunk will all fit
604 : * in the memory pool defined by GDALWarpOptions::dfWarpMemoryLimit.
605 : *
606 : * Once an appropriate region is selected GDALWarpOperation::WarpRegion()
607 : * is invoked to do the actual work.
608 : *
609 : * @param nDstXOff X offset to window of destination data to be produced.
610 : * @param nDstYOff Y offset to window of destination data to be produced.
611 : * @param nDstXSize Width of output window on destination file to be produced.
612 : * @param nDstYSize Height of output window on destination file to be produced.
613 : *
614 : * @return CE_None on success or CE_Failure if an error occurs.
615 : */
616 :
617 : struct WarpChunk {
618 : int dx, dy, dsx, dsy;
619 : int sx, sy, ssx, ssy;
620 : };
621 :
622 0 : static int OrderWarpChunk(const void* _a, const void *_b)
623 : {
624 0 : const WarpChunk* a = (const WarpChunk* )_a;
625 0 : const WarpChunk* b = (const WarpChunk* )_b;
626 0 : if (a->dy < b->dy)
627 0 : return -1;
628 0 : else if (a->dy > b->dy)
629 0 : return 1;
630 0 : else if (a->dx < b->dx)
631 0 : return -1;
632 0 : else if (a->dx > b->dx)
633 0 : return 1;
634 : else
635 0 : return 0;
636 : }
637 :
638 264 : CPLErr GDALWarpOperation::ChunkAndWarpImage(
639 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize )
640 :
641 : {
642 : /* -------------------------------------------------------------------- */
643 : /* Collect the list of chunks to operate on. */
644 : /* -------------------------------------------------------------------- */
645 264 : WipeChunkList();
646 264 : CollectChunkList( nDstXOff, nDstYOff, nDstXSize, nDstYSize );
647 :
648 : /* Sort chucks from top to bottom, and for equal y, from left to right */
649 264 : qsort(panChunkList, nChunkListCount, sizeof(WarpChunk), OrderWarpChunk);
650 :
651 : /* -------------------------------------------------------------------- */
652 : /* Total up output pixels to process. */
653 : /* -------------------------------------------------------------------- */
654 : int iChunk;
655 264 : double dfTotalPixels = 0;
656 :
657 528 : for( iChunk = 0; iChunk < nChunkListCount; iChunk++ )
658 : {
659 264 : int *panThisChunk = panChunkList + iChunk*8;
660 264 : double dfChunkPixels = panThisChunk[2] * (double) panThisChunk[3];
661 :
662 264 : dfTotalPixels += dfChunkPixels;
663 : }
664 :
665 : /* -------------------------------------------------------------------- */
666 : /* Process them one at a time, updating the progress */
667 : /* information for each region. */
668 : /* -------------------------------------------------------------------- */
669 264 : double dfPixelsProcessed=0.0;
670 :
671 528 : for( iChunk = 0; iChunk < nChunkListCount; iChunk++ )
672 : {
673 264 : int *panThisChunk = panChunkList + iChunk*8;
674 264 : double dfChunkPixels = panThisChunk[2] * (double) panThisChunk[3];
675 : CPLErr eErr;
676 :
677 264 : dfProgressBase = dfPixelsProcessed / dfTotalPixels;
678 264 : dfProgressScale = dfChunkPixels / dfTotalPixels;
679 :
680 : eErr = WarpRegion( panThisChunk[0], panThisChunk[1],
681 : panThisChunk[2], panThisChunk[3],
682 : panThisChunk[4], panThisChunk[5],
683 264 : panThisChunk[6], panThisChunk[7] );
684 :
685 264 : if( eErr != CE_None )
686 0 : return eErr;
687 :
688 264 : dfPixelsProcessed += dfChunkPixels;
689 : }
690 :
691 264 : WipeChunkList();
692 :
693 264 : psOptions->pfnProgress( 1.00001, "", psOptions->pProgressArg );
694 :
695 264 : return CE_None;
696 : }
697 :
698 : /************************************************************************/
699 : /* GDALChunkAndWarpImage() */
700 : /************************************************************************/
701 :
702 : /**
703 : * @see GDALWarpOperation::ChunkAndWarpImage()
704 : */
705 :
706 0 : CPLErr GDALChunkAndWarpImage( GDALWarpOperationH hOperation,
707 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize )
708 : {
709 0 : VALIDATE_POINTER1( hOperation, "GDALChunkAndWarpImage", CE_Failure );
710 :
711 : return ( (GDALWarpOperation *)hOperation )->
712 0 : ChunkAndWarpImage( nDstXOff, nDstYOff, nDstXSize, nDstYSize );
713 : }
714 :
715 : /************************************************************************/
716 : /* ChunkThreadMain() */
717 : /************************************************************************/
718 :
719 1 : static void ChunkThreadMain( void *pThreadData )
720 :
721 : {
722 1 : void *hThreadMutex = ((void **) pThreadData)[0];
723 : GDALWarpOperation *poOperation =
724 1 : (GDALWarpOperation *) (((void **) pThreadData)[1]);
725 1 : int *panChunkInfo = (int *) (((void **) pThreadData)[2]);
726 :
727 1 : if( !CPLAcquireMutex( hThreadMutex, 2.0 ) )
728 : {
729 : CPLError( CE_Failure, CPLE_AppDefined,
730 0 : "Failed to acquire thread mutex in ChunkThreadMain()." );
731 0 : return;
732 : }
733 :
734 : CPLErr eErr;
735 :
736 : eErr =
737 : poOperation->WarpRegion( panChunkInfo[0], panChunkInfo[1],
738 : panChunkInfo[2], panChunkInfo[3],
739 : panChunkInfo[4], panChunkInfo[5],
740 1 : panChunkInfo[6], panChunkInfo[7] );
741 :
742 : /* Return error. */
743 1 : ((void **) pThreadData)[2] = (void *) (long) eErr;
744 :
745 : /* Marks that we are done. */
746 1 : ((void **) pThreadData)[1] = NULL;
747 :
748 : /* Release mutex so parent knows we are done. */
749 1 : CPLReleaseMutex( hThreadMutex );
750 : }
751 :
752 : /************************************************************************/
753 : /* ChunkAndWarpMulti() */
754 : /************************************************************************/
755 :
756 : /**
757 : * \fn CPLErr GDALWarpOperation::ChunkAndWarpMulti(
758 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize );
759 : *
760 : * This method does a complete warp of the source image to the destination
761 : * image for the indicated region with the current warp options in effect.
762 : * Progress is reported to the installed progress monitor, if any.
763 : *
764 : * Externally this method operates the same as ChunkAndWarpImage(), but
765 : * internally this method uses multiple threads to interleave input/output
766 : * for one region while the processing is being done for another.
767 : *
768 : * @param nDstXOff X offset to window of destination data to be produced.
769 : * @param nDstYOff Y offset to window of destination data to be produced.
770 : * @param nDstXSize Width of output window on destination file to be produced.
771 : * @param nDstYSize Height of output window on destination file to be produced.
772 : *
773 : * @return CE_None on success or CE_Failure if an error occurs.
774 : */
775 :
776 1 : CPLErr GDALWarpOperation::ChunkAndWarpMulti(
777 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize )
778 :
779 : {
780 1 : hThread1Mutex = CPLCreateMutex();
781 1 : hThread2Mutex = CPLCreateMutex();
782 1 : hIOMutex = CPLCreateMutex();
783 1 : hWarpMutex = CPLCreateMutex();
784 :
785 1 : CPLReleaseMutex( hThread1Mutex );
786 1 : CPLReleaseMutex( hThread2Mutex );
787 1 : CPLReleaseMutex( hIOMutex );
788 1 : CPLReleaseMutex( hWarpMutex );
789 :
790 : /* -------------------------------------------------------------------- */
791 : /* Collect the list of chunks to operate on. */
792 : /* -------------------------------------------------------------------- */
793 1 : WipeChunkList();
794 1 : CollectChunkList( nDstXOff, nDstYOff, nDstXSize, nDstYSize );
795 :
796 : /* Sort chucks from top to bottom, and for equal y, from left to right */
797 1 : qsort(panChunkList, nChunkListCount, sizeof(WarpChunk), OrderWarpChunk);
798 :
799 : /* -------------------------------------------------------------------- */
800 : /* Process them one at a time, updating the progress */
801 : /* information for each region. */
802 : /* -------------------------------------------------------------------- */
803 : void * volatile papThreadDataList[6] = { hThread1Mutex, NULL, NULL,
804 1 : hThread2Mutex, NULL, NULL };
805 : int iChunk;
806 1 : double dfPixelsProcessed=0.0, dfTotalPixels = nDstXSize*(double)nDstYSize;
807 1 : CPLErr eErr = CE_None;
808 :
809 3 : for( iChunk = 0; iChunk < nChunkListCount+1; iChunk++ )
810 : {
811 2 : int iThread = iChunk % 2;
812 :
813 : /* -------------------------------------------------------------------- */
814 : /* Launch thread for this chunk. */
815 : /* -------------------------------------------------------------------- */
816 2 : if( iChunk < nChunkListCount )
817 : {
818 1 : int *panThisChunk = panChunkList + iChunk*8;
819 1 : double dfChunkPixels = panThisChunk[2] * (double) panThisChunk[3];
820 :
821 1 : dfProgressBase = dfPixelsProcessed / dfTotalPixels;
822 1 : dfProgressScale = dfChunkPixels / dfTotalPixels;
823 :
824 1 : dfPixelsProcessed += dfChunkPixels;
825 :
826 1 : papThreadDataList[iThread*3+1] = (void *) this;
827 1 : papThreadDataList[iThread*3+2] = (void *) panThisChunk;
828 :
829 1 : CPLDebug( "GDAL", "Start chunk %d.", iChunk );
830 1 : if( CPLCreateThread( ChunkThreadMain,
831 : (void *) (papThreadDataList + iThread*3)) == -1 )
832 : {
833 : CPLError( CE_Failure, CPLE_AppDefined,
834 0 : "CPLCreateThread() failed in ChunkAndWarpMulti()" );
835 0 : return CE_Failure;
836 : }
837 :
838 : /* Eventually we need a mechanism to ensure we wait for this
839 : thread to acquire the IO mutex before proceeding. */
840 1 : if( iChunk == 0 )
841 1 : CPLSleep( 0.25 );
842 : }
843 :
844 :
845 : /* -------------------------------------------------------------------- */
846 : /* Wait for previous chunks thread to complete. */
847 : /* -------------------------------------------------------------------- */
848 2 : if( iChunk > 0 )
849 : {
850 1 : iThread = (iChunk-1) % 2;
851 :
852 : /* Wait for thread to finish. */
853 2 : while( papThreadDataList[iThread*3+1] != NULL )
854 : {
855 0 : if( CPLAcquireMutex( papThreadDataList[iThread*3+0], 1.0 ) )
856 0 : CPLReleaseMutex( papThreadDataList[iThread*3+0] );
857 : }
858 :
859 1 : CPLDebug( "GDAL", "Finished chunk %d.", iChunk-1 );
860 :
861 1 : eErr = (CPLErr) (long) papThreadDataList[iThread*3+2];
862 :
863 1 : if( eErr != CE_None )
864 0 : break;
865 : }
866 : }
867 :
868 : /* -------------------------------------------------------------------- */
869 : /* Wait for all threads to complete. */
870 : /* -------------------------------------------------------------------- */
871 : int iThread;
872 3 : for(iThread = 0; iThread < 2; iThread ++)
873 : {
874 4 : while( papThreadDataList[iThread*3+1] != NULL )
875 : {
876 0 : if( CPLAcquireMutex( papThreadDataList[iThread*3+0], 1.0 ) )
877 0 : CPLReleaseMutex( papThreadDataList[iThread*3+0] );
878 : }
879 : }
880 :
881 1 : WipeChunkList();
882 :
883 1 : return eErr;
884 : }
885 :
886 : /************************************************************************/
887 : /* GDALChunkAndWarpMulti() */
888 : /************************************************************************/
889 :
890 : /**
891 : * @see GDALWarpOperation::ChunkAndWarpMulti()
892 : */
893 :
894 0 : CPLErr GDALChunkAndWarpMulti( GDALWarpOperationH hOperation,
895 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize )
896 : {
897 0 : VALIDATE_POINTER1( hOperation, "GDALChunkAndWarpMulti", CE_Failure );
898 :
899 : return ( (GDALWarpOperation *)hOperation )->
900 0 : ChunkAndWarpMulti( nDstXOff, nDstYOff, nDstXSize, nDstYSize );
901 : }
902 :
903 : /************************************************************************/
904 : /* WipeChunkList() */
905 : /************************************************************************/
906 :
907 834 : void GDALWarpOperation::WipeChunkList()
908 :
909 : {
910 834 : CPLFree( panChunkList );
911 834 : panChunkList = NULL;
912 834 : nChunkListCount = 0;
913 834 : nChunkListMax = 0;
914 834 : }
915 :
916 : /************************************************************************/
917 : /* CollectChunkList() */
918 : /************************************************************************/
919 :
920 265 : CPLErr GDALWarpOperation::CollectChunkList(
921 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize )
922 :
923 : {
924 : /* -------------------------------------------------------------------- */
925 : /* Compute the bounds of the input area corresponding to the */
926 : /* output area. */
927 : /* -------------------------------------------------------------------- */
928 : int nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize;
929 : CPLErr eErr;
930 :
931 : eErr = ComputeSourceWindow( nDstXOff, nDstYOff, nDstXSize, nDstYSize,
932 265 : &nSrcXOff, &nSrcYOff, &nSrcXSize, &nSrcYSize );
933 :
934 265 : if( eErr != CE_None )
935 0 : return eErr;
936 :
937 : /* -------------------------------------------------------------------- */
938 : /* If we are allowed to drop no-source regons, do so now if */
939 : /* appropriate. */
940 : /* -------------------------------------------------------------------- */
941 265 : if( (nSrcXSize == 0 || nSrcYSize == 0)
942 : && CSLFetchBoolean( psOptions->papszWarpOptions, "SKIP_NOSOURCE",0 ))
943 0 : return CE_None;
944 :
945 : /* -------------------------------------------------------------------- */
946 : /* Based on the types of masks in use, how many bits will each */
947 : /* source pixel cost us? */
948 : /* -------------------------------------------------------------------- */
949 : int nSrcPixelCostInBits;
950 :
951 : nSrcPixelCostInBits =
952 : GDALGetDataTypeSize( psOptions->eWorkingDataType )
953 265 : * psOptions->nBandCount;
954 :
955 265 : if( psOptions->pfnSrcDensityMaskFunc != NULL )
956 0 : nSrcPixelCostInBits += 32; /* float mask */
957 :
958 265 : if( psOptions->papfnSrcPerBandValidityMaskFunc != NULL
959 : || psOptions->padfSrcNoDataReal != NULL )
960 5 : nSrcPixelCostInBits += psOptions->nBandCount; /* bit/band mask */
961 :
962 265 : if( psOptions->pfnSrcValidityMaskFunc != NULL )
963 0 : nSrcPixelCostInBits += 1; /* bit mask */
964 :
965 : /* -------------------------------------------------------------------- */
966 : /* What about the cost for the destination. */
967 : /* -------------------------------------------------------------------- */
968 : int nDstPixelCostInBits;
969 :
970 : nDstPixelCostInBits =
971 : GDALGetDataTypeSize( psOptions->eWorkingDataType )
972 265 : * psOptions->nBandCount;
973 :
974 265 : if( psOptions->pfnDstDensityMaskFunc != NULL )
975 0 : nDstPixelCostInBits += 32;
976 :
977 265 : if( psOptions->padfDstNoDataReal != NULL
978 : || psOptions->pfnDstValidityMaskFunc != NULL )
979 1 : nDstPixelCostInBits += psOptions->nBandCount;
980 :
981 : /* -------------------------------------------------------------------- */
982 : /* Does the cost of the current rectangle exceed our memory */
983 : /* limit? If so, split the destination along the longest */
984 : /* dimension and recurse. */
985 : /* -------------------------------------------------------------------- */
986 : double dfTotalMemoryUse;
987 :
988 : dfTotalMemoryUse =
989 : (((double) nSrcPixelCostInBits) * nSrcXSize * nSrcYSize
990 265 : + ((double) nDstPixelCostInBits) * nDstXSize * nDstYSize) / 8.0;
991 :
992 :
993 265 : int nBlockXSize = 1, nBlockYSize = 1;
994 265 : if (psOptions->hDstDS)
995 : {
996 : GDALGetBlockSize(GDALGetRasterBand(psOptions->hDstDS, 1),
997 265 : &nBlockXSize, &nBlockYSize);
998 : }
999 :
1000 265 : if( dfTotalMemoryUse > psOptions->dfWarpMemoryLimit
1001 : && (nDstXSize > 2 || nDstYSize > 2) )
1002 : {
1003 : /* If the region width is greater than the region height, */
1004 : /* cut in half in the width. Do this only if each half part */
1005 : /* is at least as wide as the block width */
1006 0 : if( nDstXSize > nDstYSize &&
1007 : (nDstXSize / 2 >= nBlockXSize || nDstYSize == 1) )
1008 : {
1009 0 : int nChunk1 = nDstXSize / 2;
1010 :
1011 : /* Try to stick on target block boundaries */
1012 0 : if (nChunk1 > nBlockXSize)
1013 0 : nChunk1 = (nChunk1 / nBlockXSize) * nBlockXSize;
1014 :
1015 0 : int nChunk2 = nDstXSize - nChunk1;
1016 :
1017 : eErr = CollectChunkList( nDstXOff, nDstYOff,
1018 0 : nChunk1, nDstYSize );
1019 :
1020 0 : if( eErr == CE_None )
1021 : {
1022 : eErr = CollectChunkList( nDstXOff+nChunk1, nDstYOff,
1023 0 : nChunk2, nDstYSize );
1024 : }
1025 : }
1026 : else
1027 : {
1028 0 : int nChunk1 = nDstYSize / 2;
1029 0 : int nChunk2 = nDstYSize - nChunk1;
1030 :
1031 : /* Try to stick on target block boundaries */
1032 0 : if (nChunk1 > nBlockYSize)
1033 0 : nChunk1 = (nChunk1 / nBlockYSize) * nBlockYSize;
1034 :
1035 : eErr = CollectChunkList( nDstXOff, nDstYOff,
1036 0 : nDstXSize, nChunk1 );
1037 :
1038 0 : if( eErr == CE_None )
1039 : {
1040 : eErr = CollectChunkList( nDstXOff, nDstYOff+nChunk1,
1041 0 : nDstXSize, nChunk2 );
1042 : }
1043 : }
1044 :
1045 0 : return eErr;
1046 : }
1047 :
1048 : /* -------------------------------------------------------------------- */
1049 : /* OK, everything fits, so add to the chunk list. */
1050 : /* -------------------------------------------------------------------- */
1051 265 : if( nChunkListCount == nChunkListMax )
1052 : {
1053 265 : nChunkListMax = nChunkListMax * 2 + 1;
1054 : panChunkList = (int *)
1055 265 : CPLRealloc(panChunkList,sizeof(int)*nChunkListMax*8 );
1056 : }
1057 :
1058 265 : panChunkList[nChunkListCount*8+0] = nDstXOff;
1059 265 : panChunkList[nChunkListCount*8+1] = nDstYOff;
1060 265 : panChunkList[nChunkListCount*8+2] = nDstXSize;
1061 265 : panChunkList[nChunkListCount*8+3] = nDstYSize;
1062 265 : panChunkList[nChunkListCount*8+4] = nSrcXOff;
1063 265 : panChunkList[nChunkListCount*8+5] = nSrcYOff;
1064 265 : panChunkList[nChunkListCount*8+6] = nSrcXSize;
1065 265 : panChunkList[nChunkListCount*8+7] = nSrcYSize;
1066 :
1067 265 : nChunkListCount++;
1068 :
1069 265 : return CE_None;
1070 : }
1071 :
1072 :
1073 : /************************************************************************/
1074 : /* WarpRegion() */
1075 : /************************************************************************/
1076 :
1077 : /**
1078 : * \fn CPLErr GDALWarpOperation::WarpRegion(int nDstXOff, int nDstYOff,
1079 : int nDstXSize, int nDstYSize,
1080 : int nSrcXOff=0, int nSrcYOff=0,
1081 : int nSrcXSize=0, int nSrcYSize=0 );
1082 : *
1083 : * This method requests the indicated region of the output file be generated.
1084 : *
1085 : * Note that WarpRegion() will produce the requested area in one low level warp
1086 : * operation without verifying that this does not exceed the stated memory
1087 : * limits for the warp operation. Applications should take care not to call
1088 : * WarpRegion() on too large a region! This function
1089 : * is normally called by ChunkAndWarpImage(), the normal entry point for
1090 : * applications. Use it instead if staying within memory constraints is
1091 : * desired.
1092 : *
1093 : * Progress is reported from 0.0 to 1.0 for the indicated region.
1094 : *
1095 : * @param nDstXOff X offset to window of destination data to be produced.
1096 : * @param nDstYOff Y offset to window of destination data to be produced.
1097 : * @param nDstXSize Width of output window on destination file to be produced.
1098 : * @param nDstYSize Height of output window on destination file to be produced.
1099 : * @param nSrcXOff source window X offset (computed if window all zero)
1100 : * @param nSrcYOff source window Y offset (computed if window all zero)
1101 : * @param nSrcXSize source window X size (computed if window all zero)
1102 : * @param nSrcYSize source window Y size (computed if window all zero)
1103 : *
1104 : * @return CE_None on success or CE_Failure if an error occurs.
1105 : */
1106 :
1107 265 : CPLErr GDALWarpOperation::WarpRegion( int nDstXOff, int nDstYOff,
1108 : int nDstXSize, int nDstYSize,
1109 : int nSrcXOff, int nSrcYOff,
1110 : int nSrcXSize, int nSrcYSize )
1111 :
1112 : {
1113 : CPLErr eErr;
1114 : int iBand;
1115 :
1116 : /* -------------------------------------------------------------------- */
1117 : /* Acquire IO mutex. */
1118 : /* -------------------------------------------------------------------- */
1119 265 : if( hIOMutex != NULL )
1120 : {
1121 1 : if( !CPLAcquireMutex( hIOMutex, 600.0 ) )
1122 : {
1123 : CPLError( CE_Failure, CPLE_AppDefined,
1124 0 : "Failed to acquire IOMutex in WarpRegion()." );
1125 0 : return CE_Failure;
1126 : }
1127 : }
1128 :
1129 265 : ReportTiming( NULL );
1130 :
1131 : /* -------------------------------------------------------------------- */
1132 : /* Allocate the output buffer. */
1133 : /* -------------------------------------------------------------------- */
1134 : void *pDstBuffer;
1135 265 : int nWordSize = GDALGetDataTypeSize(psOptions->eWorkingDataType)/8;
1136 265 : int nBandSize = nWordSize * nDstXSize * nDstYSize;
1137 :
1138 265 : if (nDstXSize > INT_MAX / nDstYSize ||
1139 : nDstXSize * nDstYSize > INT_MAX / (nWordSize * psOptions->nBandCount))
1140 : {
1141 : CPLError( CE_Failure, CPLE_AppDefined,
1142 : "Integer overflow : nDstXSize=%d, nDstYSize=%d",
1143 0 : nDstXSize, nDstYSize);
1144 0 : return CE_Failure;
1145 : }
1146 :
1147 265 : pDstBuffer = VSIMalloc( nBandSize * psOptions->nBandCount );
1148 265 : if( pDstBuffer == NULL )
1149 : {
1150 : CPLError( CE_Failure, CPLE_OutOfMemory,
1151 : "Out of memory allocating %d byte destination buffer.",
1152 0 : nBandSize * psOptions->nBandCount );
1153 0 : return CE_Failure;
1154 : }
1155 :
1156 : /* -------------------------------------------------------------------- */
1157 : /* If the INIT_DEST option is given the initialize the output */
1158 : /* destination buffer to the indicated value without reading it */
1159 : /* from the hDstDS. This is sometimes used to optimize */
1160 : /* operation to a new output file ... it doesn't have to */
1161 : /* written out and read back for nothing. */
1162 : /* NOTE:The following code is 99% similar in gdalwarpoperation.cpp and */
1163 : /* vrtwarped.cpp. Be careful to keep it in sync ! */
1164 : /* -------------------------------------------------------------------- */
1165 : const char *pszInitDest = CSLFetchNameValue( psOptions->papszWarpOptions,
1166 265 : "INIT_DEST" );
1167 :
1168 265 : if( pszInitDest != NULL && !EQUAL(pszInitDest, "") )
1169 : {
1170 : char **papszInitValues =
1171 247 : CSLTokenizeStringComplex( pszInitDest, ",", FALSE, FALSE );
1172 247 : int nInitCount = CSLCount(papszInitValues);
1173 :
1174 506 : for( iBand = 0; iBand < psOptions->nBandCount; iBand++ )
1175 : {
1176 : double adfInitRealImag[2];
1177 : GByte *pBandData;
1178 259 : const char *pszBandInit = papszInitValues[MIN(iBand,nInitCount-1)];
1179 :
1180 260 : if( EQUAL(pszBandInit,"NO_DATA")
1181 : && psOptions->padfDstNoDataReal != NULL )
1182 : {
1183 1 : adfInitRealImag[0] = psOptions->padfDstNoDataReal[iBand];
1184 1 : adfInitRealImag[1] = psOptions->padfDstNoDataImag[iBand];
1185 : }
1186 : else
1187 : {
1188 : CPLStringToComplex( pszBandInit,
1189 258 : adfInitRealImag + 0, adfInitRealImag + 1);
1190 : }
1191 :
1192 259 : pBandData = ((GByte *) pDstBuffer) + iBand * nBandSize;
1193 :
1194 259 : if( psOptions->eWorkingDataType == GDT_Byte )
1195 : memset( pBandData,
1196 117 : MAX(0,MIN(255,(int)adfInitRealImag[0])),
1197 175 : nBandSize);
1198 402 : else if( adfInitRealImag[0] == 0.0 && adfInitRealImag[1] == 0 )
1199 : {
1200 201 : memset( pBandData, 0, nBandSize );
1201 : }
1202 0 : else if( adfInitRealImag[1] == 0.0 )
1203 : {
1204 : GDALCopyWords( &adfInitRealImag, GDT_Float64, 0,
1205 : pBandData,psOptions->eWorkingDataType,nWordSize,
1206 0 : nDstXSize * nDstYSize );
1207 : }
1208 : else
1209 : {
1210 : GDALCopyWords( &adfInitRealImag, GDT_CFloat64, 0,
1211 : pBandData,psOptions->eWorkingDataType,nWordSize,
1212 0 : nDstXSize * nDstYSize );
1213 : }
1214 : }
1215 :
1216 247 : CSLDestroy( papszInitValues );
1217 : }
1218 :
1219 : /* -------------------------------------------------------------------- */
1220 : /* If we aren't doing fixed initialization of the output buffer */
1221 : /* then read it from disk so we can overlay on existing imagery. */
1222 : /* -------------------------------------------------------------------- */
1223 265 : if( pszInitDest == NULL )
1224 : {
1225 : eErr = GDALDatasetRasterIO( psOptions->hDstDS, GF_Read,
1226 : nDstXOff, nDstYOff, nDstXSize, nDstYSize,
1227 : pDstBuffer, nDstXSize, nDstYSize,
1228 : psOptions->eWorkingDataType,
1229 : psOptions->nBandCount,
1230 : psOptions->panDstBands,
1231 18 : 0, 0, 0 );
1232 :
1233 18 : if( eErr != CE_None )
1234 : {
1235 0 : CPLFree( pDstBuffer );
1236 0 : return eErr;
1237 : }
1238 :
1239 18 : ReportTiming( "Output buffer read" );
1240 : }
1241 :
1242 : /* -------------------------------------------------------------------- */
1243 : /* Perform the warp. */
1244 : /* -------------------------------------------------------------------- */
1245 : eErr = WarpRegionToBuffer( nDstXOff, nDstYOff, nDstXSize, nDstYSize,
1246 : pDstBuffer, psOptions->eWorkingDataType,
1247 265 : nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize );
1248 :
1249 : /* -------------------------------------------------------------------- */
1250 : /* Write the output data back to disk if all went well. */
1251 : /* -------------------------------------------------------------------- */
1252 265 : if( eErr == CE_None )
1253 : {
1254 : eErr = GDALDatasetRasterIO( psOptions->hDstDS, GF_Write,
1255 : nDstXOff, nDstYOff, nDstXSize, nDstYSize,
1256 : pDstBuffer, nDstXSize, nDstYSize,
1257 : psOptions->eWorkingDataType,
1258 : psOptions->nBandCount,
1259 : psOptions->panDstBands,
1260 265 : 0, 0, 0 );
1261 265 : if( CSLFetchBoolean( psOptions->papszWarpOptions, "WRITE_FLUSH",
1262 : FALSE ) )
1263 : {
1264 0 : GDALFlushCache( psOptions->hDstDS );
1265 : }
1266 265 : ReportTiming( "Output buffer write" );
1267 : }
1268 :
1269 : /* -------------------------------------------------------------------- */
1270 : /* Cleanup and return. */
1271 : /* -------------------------------------------------------------------- */
1272 265 : VSIFree( pDstBuffer );
1273 :
1274 265 : if( hIOMutex != NULL )
1275 1 : CPLReleaseMutex( hIOMutex );
1276 :
1277 265 : return eErr;
1278 : }
1279 :
1280 : /************************************************************************/
1281 : /* GDALWarpRegion() */
1282 : /************************************************************************/
1283 :
1284 : /**
1285 : * @see GDALWarpOperation::WarpRegion()
1286 : */
1287 :
1288 0 : CPLErr GDALWarpRegion( GDALWarpOperationH hOperation,
1289 : int nDstXOff, int nDstYOff,
1290 : int nDstXSize, int nDstYSize,
1291 : int nSrcXOff, int nSrcYOff,
1292 : int nSrcXSize, int nSrcYSize )
1293 :
1294 : {
1295 0 : VALIDATE_POINTER1( hOperation, "GDALWarpRegion", CE_Failure );
1296 :
1297 : return ( (GDALWarpOperation *)hOperation )->
1298 : WarpRegion( nDstXOff, nDstYOff, nDstXSize, nDstYSize,
1299 0 : nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize);
1300 : }
1301 :
1302 : /************************************************************************/
1303 : /* WarpRegionToBuffer() */
1304 : /************************************************************************/
1305 :
1306 : /**
1307 : * \fn CPLErr GDALWarpOperation::WarpRegionToBuffer(
1308 : int nDstXOff, int nDstYOff,
1309 : int nDstXSize, int nDstYSize,
1310 : void *pDataBuf,
1311 : GDALDataType eBufDataType,
1312 : int nSrcXOff=0, int nSrcYOff=0,
1313 : int nSrcXSize=0, int nSrcYSize=0 );
1314 : *
1315 : * This method requests that a particular window of the output dataset
1316 : * be warped and the result put into the provided data buffer. The output
1317 : * dataset doesn't even really have to exist to use this method as long as
1318 : * the transformation function in the GDALWarpOptions is setup to map to
1319 : * a virtual pixel/line space.
1320 : *
1321 : * This method will do the whole region in one chunk, so be wary of the
1322 : * amount of memory that might be used.
1323 : *
1324 : * @param nDstXOff X offset to window of destination data to be produced.
1325 : * @param nDstYOff Y offset to window of destination data to be produced.
1326 : * @param nDstXSize Width of output window on destination file to be produced.
1327 : * @param nDstYSize Height of output window on destination file to be produced.
1328 : * @param pDataBuf the data buffer to place result in, of type eBufDataType.
1329 : * @param eBufDataType the type of the output data buffer. For now this
1330 : * must match GDALWarpOptions::eWorkingDataType.
1331 : * @param nSrcXOff source window X offset (computed if window all zero)
1332 : * @param nSrcYOff source window Y offset (computed if window all zero)
1333 : * @param nSrcXSize source window X size (computed if window all zero)
1334 : * @param nSrcYSize source window Y size (computed if window all zero)
1335 : *
1336 : * @return CE_None on success or CE_Failure if an error occurs.
1337 : */
1338 :
1339 330 : CPLErr GDALWarpOperation::WarpRegionToBuffer(
1340 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize,
1341 : void *pDataBuf, GDALDataType eBufDataType,
1342 : int nSrcXOff, int nSrcYOff, int nSrcXSize, int nSrcYSize )
1343 :
1344 : {
1345 330 : CPLErr eErr = CE_None;
1346 : int i;
1347 330 : int nWordSize = GDALGetDataTypeSize(psOptions->eWorkingDataType)/8;
1348 :
1349 : (void) eBufDataType;
1350 : CPLAssert( eBufDataType == psOptions->eWorkingDataType );
1351 :
1352 : /* -------------------------------------------------------------------- */
1353 : /* If not given a corresponding source window compute one now. */
1354 : /* -------------------------------------------------------------------- */
1355 330 : if( nSrcXSize == 0 && nSrcYSize == 0 )
1356 : {
1357 : eErr = ComputeSourceWindow( nDstXOff, nDstYOff, nDstXSize, nDstYSize,
1358 : &nSrcXOff, &nSrcYOff,
1359 65 : &nSrcXSize, &nSrcYSize );
1360 :
1361 65 : if( eErr != CE_None )
1362 0 : return eErr;
1363 : }
1364 :
1365 : /* -------------------------------------------------------------------- */
1366 : /* Prepare a WarpKernel object to match this operation. */
1367 : /* -------------------------------------------------------------------- */
1368 330 : GDALWarpKernel oWK;
1369 :
1370 330 : oWK.eResample = psOptions->eResampleAlg;
1371 330 : oWK.nBands = psOptions->nBandCount;
1372 330 : oWK.eWorkingDataType = psOptions->eWorkingDataType;
1373 :
1374 330 : oWK.pfnTransformer = psOptions->pfnTransformer;
1375 330 : oWK.pTransformerArg = psOptions->pTransformerArg;
1376 :
1377 330 : oWK.pfnProgress = psOptions->pfnProgress;
1378 330 : oWK.pProgress = psOptions->pProgressArg;
1379 330 : oWK.dfProgressBase = dfProgressBase;
1380 330 : oWK.dfProgressScale = dfProgressScale;
1381 :
1382 330 : oWK.papszWarpOptions = psOptions->papszWarpOptions;
1383 :
1384 330 : oWK.padfDstNoDataReal = psOptions->padfDstNoDataReal;
1385 :
1386 : /* -------------------------------------------------------------------- */
1387 : /* Setup the source buffer. */
1388 : /* */
1389 : /* Eventually we may need to take advantage of pixel */
1390 : /* interleaved reading here. */
1391 : /* -------------------------------------------------------------------- */
1392 330 : oWK.nSrcXOff = nSrcXOff;
1393 330 : oWK.nSrcYOff = nSrcYOff;
1394 330 : oWK.nSrcXSize = nSrcXSize;
1395 330 : oWK.nSrcYSize = nSrcYSize;
1396 :
1397 330 : if (nSrcXSize != 0 && nSrcYSize != 0 &&
1398 : (nSrcXSize > INT_MAX / nSrcYSize ||
1399 : nSrcXSize * nSrcYSize > INT_MAX / (nWordSize * psOptions->nBandCount)))
1400 : {
1401 : CPLError( CE_Failure, CPLE_AppDefined,
1402 : "Integer overflow : nSrcXSize=%d, nSrcYSize=%d",
1403 0 : nSrcXSize, nSrcYSize);
1404 0 : return CE_Failure;
1405 : }
1406 :
1407 : oWK.papabySrcImage = (GByte **)
1408 330 : CPLCalloc(sizeof(GByte*),psOptions->nBandCount);
1409 330 : oWK.papabySrcImage[0] = (GByte *)
1410 330 : VSIMalloc( nWordSize * nSrcXSize * nSrcYSize * psOptions->nBandCount );
1411 :
1412 330 : if( nSrcXSize != 0 && nSrcYSize != 0 && oWK.papabySrcImage[0] == NULL )
1413 : {
1414 : CPLError( CE_Failure, CPLE_OutOfMemory,
1415 : "Failed to allocate %d byte source buffer.",
1416 0 : nWordSize * nSrcXSize * nSrcYSize * psOptions->nBandCount );
1417 0 : eErr = CE_Failure;
1418 : }
1419 :
1420 694 : for( i = 0; i < psOptions->nBandCount && eErr == CE_None; i++ )
1421 364 : oWK.papabySrcImage[i] = ((GByte *) oWK.papabySrcImage[0])
1422 364 : + nWordSize * nSrcXSize * nSrcYSize * i;
1423 :
1424 330 : if( eErr == CE_None && nSrcXSize > 0 && nSrcYSize > 0 )
1425 : eErr =
1426 : GDALDatasetRasterIO( psOptions->hSrcDS, GF_Read,
1427 : nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
1428 330 : oWK.papabySrcImage[0], nSrcXSize, nSrcYSize,
1429 : psOptions->eWorkingDataType,
1430 : psOptions->nBandCount, psOptions->panSrcBands,
1431 660 : 0, 0, 0 );
1432 :
1433 330 : ReportTiming( "Input buffer read" );
1434 :
1435 : /* -------------------------------------------------------------------- */
1436 : /* Initialize destination buffer. */
1437 : /* -------------------------------------------------------------------- */
1438 330 : oWK.nDstXOff = nDstXOff;
1439 330 : oWK.nDstYOff = nDstYOff;
1440 330 : oWK.nDstXSize = nDstXSize;
1441 330 : oWK.nDstYSize = nDstYSize;
1442 :
1443 : oWK.papabyDstImage = (GByte **)
1444 330 : CPLCalloc(sizeof(GByte*),psOptions->nBandCount);
1445 :
1446 694 : for( i = 0; i < psOptions->nBandCount && eErr == CE_None; i++ )
1447 : {
1448 364 : oWK.papabyDstImage[i] = ((GByte *) pDataBuf)
1449 364 : + i * nDstXSize * nDstYSize * nWordSize;
1450 : }
1451 :
1452 : /* -------------------------------------------------------------------- */
1453 : /* Eventually we need handling for a whole bunch of the */
1454 : /* validity and density masks here. */
1455 : /* -------------------------------------------------------------------- */
1456 :
1457 : /* TODO */
1458 :
1459 : /* -------------------------------------------------------------------- */
1460 : /* Generate a source density mask if we have a source alpha band */
1461 : /* -------------------------------------------------------------------- */
1462 330 : if( eErr == CE_None && psOptions->nSrcAlphaBand > 0 )
1463 : {
1464 : CPLAssert( oWK.pafDstDensity == NULL );
1465 :
1466 7 : eErr = CreateKernelMask( &oWK, 0, "UnifiedSrcDensity" );
1467 :
1468 7 : if( eErr == CE_None )
1469 : eErr =
1470 : GDALWarpSrcAlphaMasker( psOptions,
1471 : psOptions->nBandCount,
1472 : psOptions->eWorkingDataType,
1473 : oWK.nSrcXOff, oWK.nSrcYOff,
1474 : oWK.nSrcXSize, oWK.nSrcYSize,
1475 : oWK.papabySrcImage,
1476 7 : TRUE, oWK.pafUnifiedSrcDensity );
1477 : }
1478 :
1479 : /* -------------------------------------------------------------------- */
1480 : /* Generate a source density mask if we have a source cutline. */
1481 : /* -------------------------------------------------------------------- */
1482 330 : if( eErr == CE_None && psOptions->hCutline != NULL )
1483 : {
1484 6 : if( oWK.pafUnifiedSrcDensity == NULL )
1485 : {
1486 6 : int j = oWK.nSrcXSize * oWK.nSrcYSize;
1487 :
1488 6 : eErr = CreateKernelMask( &oWK, 0, "UnifiedSrcDensity" );
1489 :
1490 6 : if( eErr == CE_None )
1491 : {
1492 60006 : for( j = oWK.nSrcXSize * oWK.nSrcYSize - 1; j >= 0; j-- )
1493 60000 : oWK.pafUnifiedSrcDensity[j] = 1.0;
1494 : }
1495 : }
1496 :
1497 6 : if( eErr == CE_None )
1498 : eErr =
1499 : GDALWarpCutlineMasker( psOptions,
1500 : psOptions->nBandCount,
1501 : psOptions->eWorkingDataType,
1502 : oWK.nSrcXOff, oWK.nSrcYOff,
1503 : oWK.nSrcXSize, oWK.nSrcYSize,
1504 : oWK.papabySrcImage,
1505 6 : TRUE, oWK.pafUnifiedSrcDensity );
1506 : }
1507 :
1508 : /* -------------------------------------------------------------------- */
1509 : /* Generate a destination density mask if we have a destination */
1510 : /* alpha band. */
1511 : /* -------------------------------------------------------------------- */
1512 330 : if( eErr == CE_None && psOptions->nDstAlphaBand > 0 )
1513 : {
1514 : CPLAssert( oWK.pafDstDensity == NULL );
1515 :
1516 9 : eErr = CreateKernelMask( &oWK, i, "DstDensity" );
1517 :
1518 9 : if( eErr == CE_None )
1519 : eErr =
1520 : GDALWarpDstAlphaMasker( psOptions,
1521 : psOptions->nBandCount,
1522 : psOptions->eWorkingDataType,
1523 : oWK.nDstXOff, oWK.nDstYOff,
1524 : oWK.nDstXSize, oWK.nDstYSize,
1525 : oWK.papabyDstImage,
1526 9 : TRUE, oWK.pafDstDensity );
1527 : }
1528 :
1529 : /* -------------------------------------------------------------------- */
1530 : /* If we have source nodata values create, or update the */
1531 : /* validity mask. */
1532 : /* -------------------------------------------------------------------- */
1533 330 : if( eErr == CE_None && psOptions->padfSrcNoDataReal != NULL )
1534 : {
1535 12 : for( i = 0; i < psOptions->nBandCount && eErr == CE_None; i++ )
1536 : {
1537 6 : eErr = CreateKernelMask( &oWK, i, "BandSrcValid" );
1538 6 : if( eErr == CE_None )
1539 : {
1540 : double adfNoData[2];
1541 :
1542 6 : adfNoData[0] = psOptions->padfSrcNoDataReal[i];
1543 6 : adfNoData[1] = psOptions->padfSrcNoDataImag[i];
1544 :
1545 : eErr =
1546 : GDALWarpNoDataMasker( adfNoData, 1,
1547 : psOptions->eWorkingDataType,
1548 : oWK.nSrcXOff, oWK.nSrcYOff,
1549 : oWK.nSrcXSize, oWK.nSrcYSize,
1550 : &(oWK.papabySrcImage[i]),
1551 6 : FALSE, oWK.papanBandSrcValid[i] );
1552 : }
1553 : }
1554 :
1555 : /* -------------------------------------------------------------------- */
1556 : /* If UNIFIED_SRC_NODATA is set, then compute a unified input */
1557 : /* pixel mask if and only if all bands nodata is true. That */
1558 : /* is, we only treat a pixel as nodata if all bands match their */
1559 : /* respective nodata values. */
1560 : /* -------------------------------------------------------------------- */
1561 6 : if( CSLFetchBoolean( psOptions->papszWarpOptions, "UNIFIED_SRC_NODATA",
1562 : FALSE )
1563 : && eErr == CE_None )
1564 : {
1565 1 : int nBytesInMask = (oWK.nSrcXSize * oWK.nSrcYSize + 31) / 8;
1566 : int iWord;
1567 :
1568 1 : eErr = CreateKernelMask( &oWK, i, "UnifiedSrcValid" );
1569 :
1570 1 : memset( oWK.panUnifiedSrcValid, 0, nBytesInMask );
1571 :
1572 2 : for( i = 0; i < psOptions->nBandCount; i++ )
1573 : {
1574 14 : for( iWord = nBytesInMask/4 - 1; iWord >= 0; iWord-- )
1575 : oWK.panUnifiedSrcValid[iWord] |=
1576 13 : oWK.papanBandSrcValid[i][iWord];
1577 1 : CPLFree( oWK.papanBandSrcValid[i] );
1578 1 : oWK.papanBandSrcValid[i] = NULL;
1579 : }
1580 :
1581 1 : CPLFree( oWK.papanBandSrcValid );
1582 1 : oWK.papanBandSrcValid = NULL;
1583 : }
1584 : }
1585 :
1586 : /* -------------------------------------------------------------------- */
1587 : /* If we have destination nodata values create, or update the */
1588 : /* validity mask. We clear the DstValid for any pixel that we */
1589 : /* do no have valid data in *any* of the source bands. */
1590 : /* */
1591 : /* Note that we don't support any concept of unified nodata on */
1592 : /* the destination image. At some point that should be added */
1593 : /* and then this logic will be significantly different. */
1594 : /* -------------------------------------------------------------------- */
1595 330 : if( eErr == CE_None && psOptions->padfDstNoDataReal != NULL )
1596 : {
1597 2 : GUInt32 *panBandMask = NULL, *panMergedMask = NULL;
1598 2 : int nMaskWords = (oWK.nDstXSize * oWK.nDstYSize + 31)/32;
1599 :
1600 2 : eErr = CreateKernelMask( &oWK, 0, "DstValid" );
1601 2 : if( eErr == CE_None )
1602 : {
1603 2 : panBandMask = (GUInt32 *) CPLMalloc(nMaskWords*4);
1604 2 : panMergedMask = (GUInt32 *) CPLCalloc(nMaskWords,4);
1605 : }
1606 :
1607 2 : if( eErr == CE_None && panBandMask != NULL )
1608 : {
1609 : int iBand, iWord;
1610 :
1611 4 : for( iBand = 0; iBand < psOptions->nBandCount; iBand++ )
1612 : {
1613 : double adfNoData[2];
1614 :
1615 2 : memset( panBandMask, 0xff, nMaskWords * 4 );
1616 :
1617 2 : adfNoData[0] = psOptions->padfDstNoDataReal[iBand];
1618 2 : adfNoData[1] = psOptions->padfDstNoDataImag[iBand];
1619 :
1620 : eErr =
1621 : GDALWarpNoDataMasker( adfNoData, 1,
1622 : psOptions->eWorkingDataType,
1623 : oWK.nDstXOff, oWK.nDstYOff,
1624 : oWK.nDstXSize, oWK.nDstYSize,
1625 : oWK.papabyDstImage + iBand,
1626 2 : FALSE, panBandMask );
1627 :
1628 2064 : for( iWord = nMaskWords - 1; iWord >= 0; iWord-- )
1629 2062 : panMergedMask[iWord] |= panBandMask[iWord];
1630 : }
1631 2 : CPLFree( panBandMask );
1632 :
1633 2064 : for( iWord = nMaskWords - 1; iWord >= 0; iWord-- )
1634 2062 : oWK.panDstValid[iWord] &= panMergedMask[iWord];
1635 :
1636 2 : CPLFree( panMergedMask );
1637 : }
1638 : }
1639 :
1640 : /* -------------------------------------------------------------------- */
1641 : /* Release IO Mutex, and acquire warper mutex. */
1642 : /* -------------------------------------------------------------------- */
1643 330 : if( hIOMutex != NULL )
1644 : {
1645 1 : CPLReleaseMutex( hIOMutex );
1646 1 : if( !CPLAcquireMutex( hWarpMutex, 600.0 ) )
1647 : {
1648 : CPLError( CE_Failure, CPLE_AppDefined,
1649 0 : "Failed to acquire WarpMutex in WarpRegion()." );
1650 0 : return CE_Failure;
1651 : }
1652 : }
1653 :
1654 : /* -------------------------------------------------------------------- */
1655 : /* Optional application provided prewarp chunk processor. */
1656 : /* -------------------------------------------------------------------- */
1657 330 : if( eErr == CE_None && psOptions->pfnPreWarpChunkProcessor != NULL )
1658 : eErr = psOptions->pfnPreWarpChunkProcessor(
1659 0 : (void *) &oWK, psOptions->pPreWarpProcessorArg );
1660 :
1661 : /* -------------------------------------------------------------------- */
1662 : /* Perform the warp. */
1663 : /* -------------------------------------------------------------------- */
1664 330 : if( eErr == CE_None )
1665 : {
1666 330 : eErr = oWK.PerformWarp();
1667 330 : ReportTiming( "In memory warp operation" );
1668 : }
1669 :
1670 : /* -------------------------------------------------------------------- */
1671 : /* Optional application provided postwarp chunk processor. */
1672 : /* -------------------------------------------------------------------- */
1673 330 : if( eErr == CE_None && psOptions->pfnPostWarpChunkProcessor != NULL )
1674 : eErr = psOptions->pfnPostWarpChunkProcessor(
1675 0 : (void *) &oWK, psOptions->pPostWarpProcessorArg );
1676 :
1677 : /* -------------------------------------------------------------------- */
1678 : /* Release Warp Mutex, and acquire io mutex. */
1679 : /* -------------------------------------------------------------------- */
1680 330 : if( hIOMutex != NULL )
1681 : {
1682 1 : CPLReleaseMutex( hWarpMutex );
1683 1 : if( !CPLAcquireMutex( hIOMutex, 600.0 ) )
1684 : {
1685 : CPLError( CE_Failure, CPLE_AppDefined,
1686 0 : "Failed to acquire IOMutex in WarpRegion()." );
1687 0 : return CE_Failure;
1688 : }
1689 : }
1690 :
1691 : /* -------------------------------------------------------------------- */
1692 : /* Write destination alpha if available. */
1693 : /* -------------------------------------------------------------------- */
1694 330 : if( eErr == CE_None && psOptions->nDstAlphaBand > 0 )
1695 : {
1696 : eErr =
1697 : GDALWarpDstAlphaMasker( psOptions,
1698 : -psOptions->nBandCount,
1699 : psOptions->eWorkingDataType,
1700 : oWK.nDstXOff, oWK.nDstYOff,
1701 : oWK.nDstXSize, oWK.nDstYSize,
1702 : oWK.papabyDstImage,
1703 9 : TRUE, oWK.pafDstDensity );
1704 : }
1705 :
1706 : /* -------------------------------------------------------------------- */
1707 : /* Cleanup. */
1708 : /* -------------------------------------------------------------------- */
1709 330 : CPLFree( oWK.papabySrcImage[0] );
1710 330 : CPLFree( oWK.papabySrcImage );
1711 330 : CPLFree( oWK.papabyDstImage );
1712 :
1713 330 : if( oWK.papanBandSrcValid != NULL )
1714 : {
1715 10 : for( i = 0; i < oWK.nBands; i++ )
1716 5 : CPLFree( oWK.papanBandSrcValid[i] );
1717 5 : CPLFree( oWK.papanBandSrcValid );
1718 : }
1719 330 : CPLFree( oWK.panUnifiedSrcValid );
1720 330 : CPLFree( oWK.pafUnifiedSrcDensity );
1721 330 : CPLFree( oWK.panDstValid );
1722 330 : CPLFree( oWK.pafDstDensity );
1723 :
1724 330 : return eErr;
1725 : }
1726 :
1727 : /************************************************************************/
1728 : /* GDALWarpRegionToBuffer() */
1729 : /************************************************************************/
1730 :
1731 : /**
1732 : * @see GDALWarpOperation::WarpRegionToBuffer()
1733 : */
1734 :
1735 0 : CPLErr GDALWarpRegionToBuffer( GDALWarpOperationH hOperation,
1736 : int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize,
1737 : void *pDataBuf, GDALDataType eBufDataType,
1738 : int nSrcXOff, int nSrcYOff, int nSrcXSize, int nSrcYSize )
1739 :
1740 : {
1741 0 : VALIDATE_POINTER1( hOperation, "GDALWarpRegionToBuffer", CE_Failure );
1742 :
1743 : return ( (GDALWarpOperation *)hOperation )->
1744 : WarpRegionToBuffer( nDstXOff, nDstYOff, nDstXSize, nDstYSize,
1745 : pDataBuf, eBufDataType,
1746 0 : nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize );
1747 : }
1748 :
1749 : /************************************************************************/
1750 : /* CreateKernelMask() */
1751 : /* */
1752 : /* If mask does not yet exist, create it. Supported types are */
1753 : /* the name of the variable in question. That is */
1754 : /* "BandSrcValid", "UnifiedSrcValid", "UnifiedSrcDensity", */
1755 : /* "DstValid", and "DstDensity". */
1756 : /************************************************************************/
1757 :
1758 31 : CPLErr GDALWarpOperation::CreateKernelMask( GDALWarpKernel *poKernel,
1759 : int iBand, const char *pszType )
1760 :
1761 : {
1762 : void **ppMask;
1763 : int nXSize, nYSize, nBitsPerPixel, nDefault;
1764 :
1765 : /* -------------------------------------------------------------------- */
1766 : /* Get particulars of mask to be updated. */
1767 : /* -------------------------------------------------------------------- */
1768 31 : if( EQUAL(pszType,"BandSrcValid") )
1769 : {
1770 6 : if( poKernel->papanBandSrcValid == NULL )
1771 : poKernel->papanBandSrcValid = (GUInt32 **)
1772 6 : CPLCalloc( sizeof(void*),poKernel->nBands);
1773 :
1774 6 : ppMask = (void **) &(poKernel->papanBandSrcValid[iBand]);
1775 6 : nXSize = poKernel->nSrcXSize;
1776 6 : nYSize = poKernel->nSrcYSize;
1777 6 : nBitsPerPixel = 1;
1778 6 : nDefault = 0xff;
1779 : }
1780 25 : else if( EQUAL(pszType,"UnifiedSrcValid") )
1781 : {
1782 1 : ppMask = (void **) &(poKernel->panUnifiedSrcValid);
1783 1 : nXSize = poKernel->nSrcXSize;
1784 1 : nYSize = poKernel->nSrcYSize;
1785 1 : nBitsPerPixel = 1;
1786 1 : nDefault = 0xff;
1787 : }
1788 24 : else if( EQUAL(pszType,"UnifiedSrcDensity") )
1789 : {
1790 13 : ppMask = (void **) &(poKernel->pafUnifiedSrcDensity);
1791 13 : nXSize = poKernel->nSrcXSize;
1792 13 : nYSize = poKernel->nSrcYSize;
1793 13 : nBitsPerPixel = 32;
1794 13 : nDefault = 0;
1795 : }
1796 11 : else if( EQUAL(pszType,"DstValid") )
1797 : {
1798 2 : ppMask = (void **) &(poKernel->panDstValid);
1799 2 : nXSize = poKernel->nDstXSize;
1800 2 : nYSize = poKernel->nDstYSize;
1801 2 : nBitsPerPixel = 1;
1802 2 : nDefault = 0xff;
1803 : }
1804 9 : else if( EQUAL(pszType,"DstDensity") )
1805 : {
1806 9 : ppMask = (void **) &(poKernel->pafDstDensity);
1807 9 : nXSize = poKernel->nDstXSize;
1808 9 : nYSize = poKernel->nDstYSize;
1809 9 : nBitsPerPixel = 32;
1810 9 : nDefault = 0;
1811 : }
1812 : else
1813 : {
1814 : CPLError( CE_Failure, CPLE_AppDefined,
1815 : "Internal error in CreateKernelMask(%s).",
1816 0 : pszType );
1817 0 : return CE_Failure;
1818 : }
1819 :
1820 : /* -------------------------------------------------------------------- */
1821 : /* Allocate if needed. */
1822 : /* -------------------------------------------------------------------- */
1823 31 : if( *ppMask == NULL )
1824 : {
1825 : int nBytes;
1826 :
1827 31 : if( nBitsPerPixel == 32 )
1828 22 : nBytes = nXSize * nYSize * 4;
1829 : else
1830 9 : nBytes = (nXSize * nYSize + 31) / 8;
1831 :
1832 31 : *ppMask = VSIMalloc( nBytes );
1833 :
1834 31 : if( *ppMask == NULL )
1835 : {
1836 : CPLError( CE_Failure, CPLE_OutOfMemory,
1837 : "Out of memory allocating %d bytes for %s mask.",
1838 0 : nBytes, pszType );
1839 0 : return CE_Failure;
1840 : }
1841 :
1842 31 : memset( *ppMask, nDefault, nBytes );
1843 : }
1844 :
1845 31 : return CE_None;
1846 : }
1847 :
1848 :
1849 :
1850 : /************************************************************************/
1851 : /* ComputeSourceWindow() */
1852 : /************************************************************************/
1853 :
1854 330 : CPLErr GDALWarpOperation::ComputeSourceWindow(int nDstXOff, int nDstYOff,
1855 : int nDstXSize, int nDstYSize,
1856 : int *pnSrcXOff, int *pnSrcYOff,
1857 : int *pnSrcXSize, int *pnSrcYSize)
1858 :
1859 : {
1860 : /* -------------------------------------------------------------------- */
1861 : /* Figure out whether we just want to do the usual "along the */
1862 : /* edge" sampling, or using a grid. The grid usage is */
1863 : /* important in some weird "inside out" cases like WGS84 to */
1864 : /* polar stereographic around the pole. Also figure out the */
1865 : /* sampling rate. */
1866 : /* -------------------------------------------------------------------- */
1867 : double dfStepSize;
1868 330 : int nSampleMax, nStepCount = 21, bUseGrid;
1869 330 : int *pabSuccess = NULL;
1870 : double *padfX, *padfY, *padfZ;
1871 : int nSamplePoints;
1872 : double dfRatio;
1873 :
1874 330 : if( CSLFetchNameValue( psOptions->papszWarpOptions,
1875 : "SAMPLE_STEPS" ) != NULL )
1876 : {
1877 : nStepCount =
1878 : atoi(CSLFetchNameValue( psOptions->papszWarpOptions,
1879 0 : "SAMPLE_STEPS" ));
1880 0 : nStepCount = MAX(2,nStepCount);
1881 : }
1882 :
1883 330 : dfStepSize = 1.0 / (nStepCount-1);
1884 :
1885 : bUseGrid = CSLFetchBoolean( psOptions->papszWarpOptions,
1886 330 : "SAMPLE_GRID", FALSE );
1887 :
1888 : TryAgainWithGrid:
1889 335 : nSamplePoints = 0;
1890 335 : if( bUseGrid )
1891 : {
1892 5 : if (nStepCount > INT_MAX / nStepCount)
1893 : {
1894 : CPLError( CE_Failure, CPLE_AppDefined,
1895 0 : "Too many steps : %d", nStepCount);
1896 0 : return CE_Failure;
1897 : }
1898 5 : nSampleMax = nStepCount * nStepCount;
1899 : }
1900 : else
1901 : {
1902 330 : if (nStepCount > INT_MAX / 4)
1903 : {
1904 : CPLError( CE_Failure, CPLE_AppDefined,
1905 0 : "Too many steps : %d", nStepCount);
1906 0 : return CE_Failure;
1907 : }
1908 330 : nSampleMax = nStepCount * 4;
1909 : }
1910 :
1911 335 : pabSuccess = (int *) VSIMalloc2(sizeof(int), nSampleMax);
1912 335 : padfX = (double *) VSIMalloc2(sizeof(double) * 3, nSampleMax);
1913 335 : if (pabSuccess == NULL || padfX == NULL)
1914 : {
1915 0 : CPLFree( padfX );
1916 0 : CPLFree( pabSuccess );
1917 0 : return CE_Failure;
1918 : }
1919 335 : padfY = padfX + nSampleMax;
1920 335 : padfZ = padfX + nSampleMax * 2;
1921 :
1922 : /* -------------------------------------------------------------------- */
1923 : /* Setup sample points on a grid pattern throughout the area. */
1924 : /* -------------------------------------------------------------------- */
1925 335 : if( bUseGrid )
1926 : {
1927 : double dfRatioY;
1928 :
1929 110 : for( dfRatioY = 0.0;
1930 : dfRatioY <= 1.0 + dfStepSize*0.5;
1931 : dfRatioY += dfStepSize )
1932 : {
1933 2310 : for( dfRatio = 0.0;
1934 : dfRatio <= 1.0 + dfStepSize*0.5;
1935 : dfRatio += dfStepSize )
1936 : {
1937 2205 : padfX[nSamplePoints] = dfRatio * nDstXSize + nDstXOff;
1938 2205 : padfY[nSamplePoints] = dfRatioY * nDstYSize + nDstYOff;
1939 2205 : padfZ[nSamplePoints++] = 0.0;
1940 : }
1941 : }
1942 : }
1943 : /* -------------------------------------------------------------------- */
1944 : /* Setup sample points all around the edge of the output raster. */
1945 : /* -------------------------------------------------------------------- */
1946 : else
1947 : {
1948 7260 : for( dfRatio = 0.0; dfRatio <= 1.0 + dfStepSize*0.5; dfRatio += dfStepSize )
1949 : {
1950 : // Along top
1951 6930 : padfX[nSamplePoints] = dfRatio * nDstXSize + nDstXOff;
1952 6930 : padfY[nSamplePoints] = nDstYOff;
1953 6930 : padfZ[nSamplePoints++] = 0.0;
1954 :
1955 : // Along bottom
1956 6930 : padfX[nSamplePoints] = dfRatio * nDstXSize + nDstXOff;
1957 6930 : padfY[nSamplePoints] = nDstYOff + nDstYSize;
1958 6930 : padfZ[nSamplePoints++] = 0.0;
1959 :
1960 : // Along left
1961 6930 : padfX[nSamplePoints] = nDstXOff;
1962 6930 : padfY[nSamplePoints] = dfRatio * nDstYSize + nDstYOff;
1963 6930 : padfZ[nSamplePoints++] = 0.0;
1964 :
1965 : // Along right
1966 6930 : padfX[nSamplePoints] = nDstXSize + nDstXOff;
1967 6930 : padfY[nSamplePoints] = dfRatio * nDstYSize + nDstYOff;
1968 6930 : padfZ[nSamplePoints++] = 0.0;
1969 : }
1970 : }
1971 :
1972 : CPLAssert( nSamplePoints == nSampleMax );
1973 :
1974 : /* -------------------------------------------------------------------- */
1975 : /* Transform them to the input pixel coordinate space */
1976 : /* -------------------------------------------------------------------- */
1977 335 : if( !psOptions->pfnTransformer( psOptions->pTransformerArg,
1978 : TRUE, nSamplePoints,
1979 : padfX, padfY, padfZ, pabSuccess ) )
1980 : {
1981 0 : CPLFree( padfX );
1982 0 : CPLFree( pabSuccess );
1983 :
1984 : CPLError( CE_Failure, CPLE_AppDefined,
1985 : "GDALWarperOperation::ComputeSourceWindow() failed because\n"
1986 0 : "the pfnTransformer failed." );
1987 0 : return CE_Failure;
1988 : }
1989 :
1990 : /* -------------------------------------------------------------------- */
1991 : /* Collect the bounds, ignoring any failed points. */
1992 : /* -------------------------------------------------------------------- */
1993 335 : double dfMinXOut=0.0, dfMinYOut=0.0, dfMaxXOut=0.0, dfMaxYOut=0.0;
1994 335 : int bGotInitialPoint = FALSE;
1995 335 : int nFailedCount = 0, i;
1996 :
1997 30260 : for( i = 0; i < nSamplePoints; i++ )
1998 : {
1999 29925 : if( !pabSuccess[i] )
2000 : {
2001 1341 : nFailedCount++;
2002 1341 : continue;
2003 : }
2004 :
2005 28584 : if( !bGotInitialPoint )
2006 : {
2007 334 : bGotInitialPoint = TRUE;
2008 334 : dfMinXOut = dfMaxXOut = padfX[i];
2009 334 : dfMinYOut = dfMaxYOut = padfY[i];
2010 : }
2011 : else
2012 : {
2013 28250 : dfMinXOut = MIN(dfMinXOut,padfX[i]);
2014 28250 : dfMinYOut = MIN(dfMinYOut,padfY[i]);
2015 28250 : dfMaxXOut = MAX(dfMaxXOut,padfX[i]);
2016 28250 : dfMaxYOut = MAX(dfMaxYOut,padfY[i]);
2017 : }
2018 : }
2019 :
2020 335 : CPLFree( padfX );
2021 335 : CPLFree( pabSuccess );
2022 :
2023 : /* -------------------------------------------------------------------- */
2024 : /* If we got any failures when not using a grid, we should */
2025 : /* really go back and try again with the grid. Sorry for the */
2026 : /* goto. */
2027 : /* -------------------------------------------------------------------- */
2028 335 : if( !bUseGrid && nFailedCount > 0 )
2029 : {
2030 5 : bUseGrid = TRUE;
2031 5 : goto TryAgainWithGrid;
2032 : }
2033 :
2034 : /* -------------------------------------------------------------------- */
2035 : /* If we get hardly any points (or none) transforming, we give */
2036 : /* up. */
2037 : /* -------------------------------------------------------------------- */
2038 330 : if( nFailedCount > nSamplePoints - 5 )
2039 : {
2040 : CPLError( CE_Failure, CPLE_AppDefined,
2041 : "Too many points (%d out of %d) failed to transform,\n"
2042 : "unable to compute output bounds.",
2043 0 : nFailedCount, nSamplePoints );
2044 0 : return CE_Failure;
2045 : }
2046 :
2047 330 : if( nFailedCount > 0 )
2048 : CPLDebug( "GDAL",
2049 : "GDALWarpOperation::ComputeSourceWindow() %d out of %d points failed to transform.",
2050 5 : nFailedCount, nSamplePoints );
2051 :
2052 : /* -------------------------------------------------------------------- */
2053 : /* How much of a window around our source pixel might we need */
2054 : /* to collect data from based on the resampling kernel? Even */
2055 : /* if the requested central pixel falls off the source image, */
2056 : /* we may need to collect data if some portion of the */
2057 : /* resampling kernel could be on-image. */
2058 : /* -------------------------------------------------------------------- */
2059 330 : int nResWinSize = GWKGetFilterRadius(psOptions->eResampleAlg);
2060 :
2061 : /* -------------------------------------------------------------------- */
2062 : /* Allow addition of extra sample pixels to source window to */
2063 : /* avoid missing pixels due to sampling error. In fact, */
2064 : /* fallback to adding a bit to the window if any points failed */
2065 : /* to transform. */
2066 : /* -------------------------------------------------------------------- */
2067 330 : if( CSLFetchNameValue( psOptions->papszWarpOptions,
2068 : "SOURCE_EXTRA" ) != NULL )
2069 : {
2070 : nResWinSize += atoi(
2071 0 : CSLFetchNameValue( psOptions->papszWarpOptions, "SOURCE_EXTRA" ));
2072 : }
2073 330 : else if( nFailedCount > 0 )
2074 5 : nResWinSize += 10;
2075 :
2076 : /* -------------------------------------------------------------------- */
2077 : /* return bounds. */
2078 : /* -------------------------------------------------------------------- */
2079 330 : *pnSrcXOff = MAX(0,(int) floor( dfMinXOut ) - nResWinSize );
2080 330 : *pnSrcYOff = MAX(0,(int) floor( dfMinYOut ) - nResWinSize );
2081 330 : *pnSrcXOff = MIN(*pnSrcXOff,GDALGetRasterXSize(psOptions->hSrcDS));
2082 330 : *pnSrcYOff = MIN(*pnSrcYOff,GDALGetRasterYSize(psOptions->hSrcDS));
2083 :
2084 : *pnSrcXSize = MIN( GDALGetRasterXSize(psOptions->hSrcDS) - *pnSrcXOff,
2085 330 : ((int) ceil( dfMaxXOut )) - *pnSrcXOff + nResWinSize );
2086 : *pnSrcYSize = MIN( GDALGetRasterYSize(psOptions->hSrcDS) - *pnSrcYOff,
2087 330 : ((int) ceil( dfMaxYOut )) - *pnSrcYOff + nResWinSize );
2088 330 : *pnSrcXSize = MAX(0,*pnSrcXSize);
2089 330 : *pnSrcYSize = MAX(0,*pnSrcYSize);
2090 :
2091 :
2092 330 : return CE_None;
2093 : }
2094 :
2095 : /************************************************************************/
2096 : /* ReportTiming() */
2097 : /************************************************************************/
2098 :
2099 1208 : void GDALWarpOperation::ReportTiming( const char * pszMessage )
2100 :
2101 : {
2102 1208 : if( !bReportTimings )
2103 1208 : return;
2104 :
2105 0 : unsigned long nNewTime = VSITime(NULL);
2106 :
2107 0 : if( pszMessage != NULL )
2108 : {
2109 : CPLDebug( "WARP_TIMING", "%s: %lds",
2110 0 : pszMessage, (long)(nNewTime - nLastTimeReported) );
2111 : }
2112 :
2113 0 : nLastTimeReported = nNewTime;
2114 : }
2115 :
|
