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