1 : /******************************************************************************
2 : * $Id: vrtsources.cpp 23465 2011-12-04 15:54:52Z rouault $
3 : *
4 : * Project: Virtual GDAL Datasets
5 : * Purpose: Implementation of VRTSimpleSource, VRTFuncSource and
6 : * VRTAveragedSource.
7 : * Author: Frank Warmerdam <warmerdam@pobox.com>
8 : *
9 : ******************************************************************************
10 : * Copyright (c) 2001, Frank Warmerdam <warmerdam@pobox.com>
11 : *
12 : * Permission is hereby granted, free of charge, to any person obtaining a
13 : * copy of this software and associated documentation files (the "Software"),
14 : * to deal in the Software without restriction, including without limitation
15 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 : * and/or sell copies of the Software, and to permit persons to whom the
17 : * Software is furnished to do so, subject to the following conditions:
18 : *
19 : * The above copyright notice and this permission notice shall be included
20 : * in all copies or substantial portions of the Software.
21 : *
22 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 : * DEALINGS IN THE SOFTWARE.
29 : ****************************************************************************/
30 :
31 : #include "vrtdataset.h"
32 : #include "gdal_proxy.h"
33 : #include "cpl_minixml.h"
34 : #include "cpl_string.h"
35 :
36 : #include <algorithm>
37 :
38 : CPL_CVSID("$Id: vrtsources.cpp 23465 2011-12-04 15:54:52Z rouault $");
39 :
40 : /************************************************************************/
41 : /* ==================================================================== */
42 : /* VRTSource */
43 : /* ==================================================================== */
44 : /************************************************************************/
45 :
46 631 : VRTSource::~VRTSource()
47 : {
48 631 : }
49 :
50 : /************************************************************************/
51 : /* GetFileList() */
52 : /************************************************************************/
53 :
54 0 : void VRTSource::GetFileList(char*** ppapszFileList, int *pnSize,
55 : int *pnMaxSize, CPLHashSet* hSetFiles)
56 : {
57 0 : }
58 :
59 : /************************************************************************/
60 : /* ==================================================================== */
61 : /* VRTSimpleSource */
62 : /* ==================================================================== */
63 : /************************************************************************/
64 :
65 : /************************************************************************/
66 : /* VRTSimpleSource() */
67 : /************************************************************************/
68 :
69 631 : VRTSimpleSource::VRTSimpleSource()
70 :
71 : {
72 631 : poRasterBand = NULL;
73 631 : poMaskBandMainBand = NULL;
74 631 : bNoDataSet = FALSE;
75 631 : dfNoDataValue = VRT_NODATA_UNSET;
76 631 : }
77 :
78 : /************************************************************************/
79 : /* ~VRTSimpleSource() */
80 : /************************************************************************/
81 :
82 631 : VRTSimpleSource::~VRTSimpleSource()
83 :
84 : {
85 631 : if( poMaskBandMainBand != NULL )
86 : {
87 13 : if (poMaskBandMainBand->GetDataset() != NULL )
88 : {
89 13 : if( poMaskBandMainBand->GetDataset()->GetShared() )
90 12 : GDALClose( (GDALDatasetH) poMaskBandMainBand->GetDataset() );
91 : else
92 1 : poMaskBandMainBand->GetDataset()->Dereference();
93 : }
94 : }
95 618 : else if( poRasterBand != NULL && poRasterBand->GetDataset() != NULL )
96 : {
97 617 : if( poRasterBand->GetDataset()->GetShared() )
98 519 : GDALClose( (GDALDatasetH) poRasterBand->GetDataset() );
99 : else
100 98 : poRasterBand->GetDataset()->Dereference();
101 : }
102 631 : }
103 :
104 : /************************************************************************/
105 : /* SetSrcBand() */
106 : /************************************************************************/
107 :
108 289 : void VRTSimpleSource::SetSrcBand( GDALRasterBand *poNewSrcBand )
109 :
110 : {
111 289 : poRasterBand = poNewSrcBand;
112 289 : }
113 :
114 :
115 : /************************************************************************/
116 : /* SetSrcMaskBand() */
117 : /************************************************************************/
118 :
119 : /* poSrcBand is not the mask band, but the band from which the mask band is taken */
120 5 : void VRTSimpleSource::SetSrcMaskBand( GDALRasterBand *poNewSrcBand )
121 :
122 : {
123 5 : poRasterBand = poNewSrcBand->GetMaskBand();
124 5 : poMaskBandMainBand = poNewSrcBand;
125 5 : }
126 :
127 : /************************************************************************/
128 : /* SetSrcWindow() */
129 : /************************************************************************/
130 :
131 294 : void VRTSimpleSource::SetSrcWindow( int nNewXOff, int nNewYOff,
132 : int nNewXSize, int nNewYSize )
133 :
134 : {
135 294 : nSrcXOff = nNewXOff;
136 294 : nSrcYOff = nNewYOff;
137 294 : nSrcXSize = nNewXSize;
138 294 : nSrcYSize = nNewYSize;
139 294 : }
140 :
141 : /************************************************************************/
142 : /* SetDstWindow() */
143 : /************************************************************************/
144 :
145 294 : void VRTSimpleSource::SetDstWindow( int nNewXOff, int nNewYOff,
146 : int nNewXSize, int nNewYSize )
147 :
148 : {
149 294 : nDstXOff = nNewXOff;
150 294 : nDstYOff = nNewYOff;
151 294 : nDstXSize = nNewXSize;
152 294 : nDstYSize = nNewYSize;
153 294 : }
154 :
155 : /************************************************************************/
156 : /* SetNoDataValue() */
157 : /************************************************************************/
158 :
159 8 : void VRTSimpleSource::SetNoDataValue( double dfNewNoDataValue )
160 :
161 : {
162 8 : if( dfNewNoDataValue == VRT_NODATA_UNSET )
163 : {
164 0 : bNoDataSet = FALSE;
165 0 : dfNoDataValue = VRT_NODATA_UNSET;
166 : }
167 : else
168 : {
169 8 : bNoDataSet = TRUE;
170 8 : dfNoDataValue = dfNewNoDataValue;
171 : }
172 8 : }
173 :
174 : /************************************************************************/
175 : /* SerializeToXML() */
176 : /************************************************************************/
177 :
178 181 : CPLXMLNode *VRTSimpleSource::SerializeToXML( const char *pszVRTPath )
179 :
180 : {
181 : CPLXMLNode *psSrc;
182 : int bRelativeToVRT;
183 : const char *pszRelativePath;
184 : int nBlockXSize, nBlockYSize;
185 :
186 181 : if( poRasterBand == NULL )
187 0 : return NULL;
188 :
189 : GDALDataset *poDS;
190 :
191 181 : if (poMaskBandMainBand)
192 : {
193 5 : poDS = poMaskBandMainBand->GetDataset();
194 5 : if( poDS == NULL || poMaskBandMainBand->GetBand() < 1 )
195 0 : return NULL;
196 : }
197 : else
198 : {
199 176 : poDS = poRasterBand->GetDataset();
200 176 : if( poDS == NULL || poRasterBand->GetBand() < 1 )
201 0 : return NULL;
202 : }
203 :
204 181 : psSrc = CPLCreateXMLNode( NULL, CXT_Element, "SimpleSource" );
205 :
206 : VSIStatBufL sStat;
207 362 : if ( strstr(poDS->GetDescription(), "/vsicurl/http") != NULL ||
208 181 : strstr(poDS->GetDescription(), "/vsicurl/ftp") != NULL )
209 : {
210 : /* Testing the existence of remote ressources can be excruciating */
211 : /* slow, so let's just suppose they exist */
212 0 : pszRelativePath = poDS->GetDescription();
213 0 : bRelativeToVRT = FALSE;
214 : }
215 : /* If this isn't actually a file, don't even try to know if it is */
216 : /* a relative path. It can't be !, and unfortunately */
217 : /* CPLIsFilenameRelative() can only work with strings that are filenames */
218 : /* To be clear NITF_TOC_ENTRY:CADRG_JOG-A_250K_1_0:some_path isn't a relative */
219 : /* file path */
220 181 : else if( VSIStatExL( poDS->GetDescription(), &sStat, VSI_STAT_EXISTS_FLAG ) != 0 )
221 : {
222 1 : pszRelativePath = poDS->GetDescription();
223 1 : bRelativeToVRT = FALSE;
224 : }
225 : else
226 : {
227 : pszRelativePath =
228 180 : CPLExtractRelativePath( pszVRTPath, poDS->GetDescription(),
229 360 : &bRelativeToVRT );
230 : }
231 :
232 181 : CPLSetXMLValue( psSrc, "SourceFilename", pszRelativePath );
233 :
234 : CPLCreateXMLNode(
235 : CPLCreateXMLNode( CPLGetXMLNode( psSrc, "SourceFilename" ),
236 : CXT_Attribute, "relativeToVRT" ),
237 181 : CXT_Text, bRelativeToVRT ? "1" : "0" );
238 :
239 181 : if (poMaskBandMainBand)
240 : CPLSetXMLValue( psSrc, "SourceBand",
241 5 : CPLSPrintf("mask,%d",poMaskBandMainBand->GetBand()) );
242 : else
243 : CPLSetXMLValue( psSrc, "SourceBand",
244 176 : CPLSPrintf("%d",poRasterBand->GetBand()) );
245 :
246 : /* Write a few additional useful properties of the dataset */
247 : /* so that we can use a proxy dataset when re-opening. See XMLInit() */
248 : /* below */
249 : CPLSetXMLValue( psSrc, "SourceProperties.#RasterXSize",
250 181 : CPLSPrintf("%d",poRasterBand->GetXSize()) );
251 : CPLSetXMLValue( psSrc, "SourceProperties.#RasterYSize",
252 181 : CPLSPrintf("%d",poRasterBand->GetYSize()) );
253 : CPLSetXMLValue( psSrc, "SourceProperties.#DataType",
254 181 : GDALGetDataTypeName( poRasterBand->GetRasterDataType() ) );
255 181 : poRasterBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
256 : CPLSetXMLValue( psSrc, "SourceProperties.#BlockXSize",
257 181 : CPLSPrintf("%d",nBlockXSize) );
258 : CPLSetXMLValue( psSrc, "SourceProperties.#BlockYSize",
259 181 : CPLSPrintf("%d",nBlockYSize) );
260 :
261 181 : if( nSrcXOff != -1
262 : || nSrcYOff != -1
263 : || nSrcXSize != -1
264 : || nSrcYSize != -1 )
265 : {
266 : CPLSetXMLValue( psSrc, "SrcRect.#xOff",
267 178 : CPLSPrintf( "%d", nSrcXOff ) );
268 : CPLSetXMLValue( psSrc, "SrcRect.#yOff",
269 178 : CPLSPrintf( "%d", nSrcYOff ) );
270 : CPLSetXMLValue( psSrc, "SrcRect.#xSize",
271 178 : CPLSPrintf( "%d", nSrcXSize ) );
272 : CPLSetXMLValue( psSrc, "SrcRect.#ySize",
273 178 : CPLSPrintf( "%d", nSrcYSize ) );
274 : }
275 :
276 181 : if( nDstXOff != -1
277 : || nDstYOff != -1
278 : || nDstXSize != -1
279 : || nDstYSize != -1 )
280 : {
281 178 : CPLSetXMLValue( psSrc, "DstRect.#xOff", CPLSPrintf( "%d", nDstXOff ) );
282 178 : CPLSetXMLValue( psSrc, "DstRect.#yOff", CPLSPrintf( "%d", nDstYOff ) );
283 178 : CPLSetXMLValue( psSrc, "DstRect.#xSize",CPLSPrintf( "%d", nDstXSize ));
284 178 : CPLSetXMLValue( psSrc, "DstRect.#ySize",CPLSPrintf( "%d", nDstYSize ));
285 : }
286 :
287 181 : return psSrc;
288 : }
289 :
290 : /************************************************************************/
291 : /* XMLInit() */
292 : /************************************************************************/
293 :
294 337 : CPLErr VRTSimpleSource::XMLInit( CPLXMLNode *psSrc, const char *pszVRTPath )
295 :
296 : {
297 : /* -------------------------------------------------------------------- */
298 : /* Prepare filename. */
299 : /* -------------------------------------------------------------------- */
300 337 : char *pszSrcDSName = NULL;
301 337 : CPLXMLNode* psSourceFileNameNode = CPLGetXMLNode(psSrc,"SourceFilename");
302 : const char *pszFilename =
303 337 : psSourceFileNameNode ? CPLGetXMLValue(psSourceFileNameNode,NULL, NULL) : NULL;
304 :
305 337 : if( pszFilename == NULL )
306 : {
307 : CPLError( CE_Warning, CPLE_AppDefined,
308 0 : "Missing <SourceFilename> element in VRTRasterBand." );
309 0 : return CE_Failure;
310 : }
311 :
312 337 : if( pszVRTPath != NULL
313 : && atoi(CPLGetXMLValue( psSourceFileNameNode, "relativetoVRT", "0")) )
314 : {
315 : pszSrcDSName = CPLStrdup(
316 193 : CPLProjectRelativeFilename( pszVRTPath, pszFilename ) );
317 : }
318 : else
319 144 : pszSrcDSName = CPLStrdup( pszFilename );
320 :
321 337 : const char* pszSourceBand = CPLGetXMLValue(psSrc,"SourceBand","1");
322 337 : int nSrcBand = 0;
323 337 : int bGetMaskBand = FALSE;
324 337 : if (EQUALN(pszSourceBand, "mask",4))
325 : {
326 8 : bGetMaskBand = TRUE;
327 8 : if (pszSourceBand[4] == ',')
328 8 : nSrcBand = atoi(pszSourceBand + 5);
329 : else
330 0 : nSrcBand = 1;
331 : }
332 : else
333 329 : nSrcBand = atoi(pszSourceBand);
334 337 : if (!GDALCheckBandCount(nSrcBand, 0))
335 : {
336 : CPLError( CE_Warning, CPLE_AppDefined,
337 0 : "Invalid <SourceBand> element in VRTRasterBand." );
338 0 : CPLFree( pszSrcDSName );
339 0 : return CE_Failure;
340 : }
341 :
342 : /* Newly generated VRT will have RasterXSize, RasterYSize, DataType, */
343 : /* BlockXSize, BlockYSize tags, so that we don't have actually to */
344 : /* open the real dataset immediately, but we can use a proxy dataset */
345 : /* instead. This is particularly usefull when dealing with huge VRT */
346 : /* For example, a VRT with the world coverage of DTED0 (25594 files) */
347 337 : CPLXMLNode* psSrcProperties = CPLGetXMLNode(psSrc,"SourceProperties");
348 337 : int nRasterXSize = 0, nRasterYSize =0;
349 337 : GDALDataType eDataType = (GDALDataType)-1;
350 337 : int nBlockXSize = 0, nBlockYSize = 0;
351 337 : if (psSrcProperties)
352 : {
353 219 : nRasterXSize = atoi(CPLGetXMLValue(psSrcProperties,"RasterXSize","0"));
354 219 : nRasterYSize = atoi(CPLGetXMLValue(psSrcProperties,"RasterYSize","0"));
355 219 : const char *pszDataType = CPLGetXMLValue(psSrcProperties, "DataType", NULL);
356 219 : if( pszDataType != NULL )
357 : {
358 489 : for( int iType = 0; iType < GDT_TypeCount; iType++ )
359 : {
360 489 : const char *pszThisName = GDALGetDataTypeName((GDALDataType)iType);
361 :
362 489 : if( pszThisName != NULL && EQUAL(pszDataType,pszThisName) )
363 : {
364 219 : eDataType = (GDALDataType) iType;
365 219 : break;
366 : }
367 : }
368 : }
369 219 : nBlockXSize = atoi(CPLGetXMLValue(psSrcProperties,"BlockXSize","0"));
370 219 : nBlockYSize = atoi(CPLGetXMLValue(psSrcProperties,"BlockYSize","0"));
371 : }
372 :
373 : GDALDataset *poSrcDS;
374 455 : if (nRasterXSize == 0 || nRasterYSize == 0 || eDataType == (GDALDataType)-1 ||
375 : nBlockXSize == 0 || nBlockYSize == 0)
376 : {
377 : /* -------------------------------------------------------------------- */
378 : /* Open the file (shared). */
379 : /* -------------------------------------------------------------------- */
380 118 : poSrcDS = (GDALDataset *) GDALOpenShared( pszSrcDSName, GA_ReadOnly );
381 : }
382 : else
383 : {
384 : /* -------------------------------------------------------------------- */
385 : /* Create a proxy dataset */
386 : /* -------------------------------------------------------------------- */
387 : int i;
388 219 : GDALProxyPoolDataset* proxyDS = new GDALProxyPoolDataset(pszSrcDSName, nRasterXSize, nRasterYSize, GA_ReadOnly, TRUE);
389 219 : poSrcDS = proxyDS;
390 :
391 : /* Only the information of rasterBand nSrcBand will be accurate */
392 : /* but that's OK since we only use that band afterwards */
393 492 : for(i=1;i<=nSrcBand;i++)
394 273 : proxyDS->AddSrcBandDescription(eDataType, nBlockXSize, nBlockYSize);
395 219 : if (bGetMaskBand)
396 8 : ((GDALProxyPoolRasterBand*)proxyDS->GetRasterBand(nSrcBand))->AddSrcMaskBandDescription(eDataType, nBlockXSize, nBlockYSize);
397 : }
398 :
399 337 : CPLFree( pszSrcDSName );
400 :
401 337 : if( poSrcDS == NULL )
402 1 : return CE_Failure;
403 :
404 : /* -------------------------------------------------------------------- */
405 : /* Get the raster band. */
406 : /* -------------------------------------------------------------------- */
407 :
408 336 : poRasterBand = poSrcDS->GetRasterBand(nSrcBand);
409 336 : if( poRasterBand == NULL )
410 : {
411 0 : if( poSrcDS->GetShared() )
412 0 : GDALClose( (GDALDatasetH) poSrcDS );
413 0 : return CE_Failure;
414 : }
415 336 : if (bGetMaskBand)
416 : {
417 8 : poMaskBandMainBand = poRasterBand;
418 8 : poRasterBand = poRasterBand->GetMaskBand();
419 8 : if( poRasterBand == NULL )
420 0 : return CE_Failure;
421 : }
422 :
423 : /* -------------------------------------------------------------------- */
424 : /* Set characteristics. */
425 : /* -------------------------------------------------------------------- */
426 336 : CPLXMLNode* psSrcRect = CPLGetXMLNode(psSrc,"SrcRect");
427 336 : if (psSrcRect)
428 : {
429 323 : nSrcXOff = atoi(CPLGetXMLValue(psSrcRect,"xOff","-1"));
430 323 : nSrcYOff = atoi(CPLGetXMLValue(psSrcRect,"yOff","-1"));
431 323 : nSrcXSize = atoi(CPLGetXMLValue(psSrcRect,"xSize","-1"));
432 323 : nSrcYSize = atoi(CPLGetXMLValue(psSrcRect,"ySize","-1"));
433 : }
434 : else
435 : {
436 13 : nSrcXOff = nSrcYOff = nSrcXSize = nSrcYSize = -1;
437 : }
438 :
439 336 : CPLXMLNode* psDstRect = CPLGetXMLNode(psSrc,"DstRect");
440 336 : if (psDstRect)
441 : {
442 323 : nDstXOff = atoi(CPLGetXMLValue(psDstRect,"xOff","-1"));
443 323 : nDstYOff = atoi(CPLGetXMLValue(psDstRect,"yOff","-1"));
444 323 : nDstXSize = atoi(CPLGetXMLValue(psDstRect,"xSize","-1"));
445 323 : nDstYSize = atoi(CPLGetXMLValue(psDstRect,"ySize","-1"));
446 : }
447 : else
448 : {
449 13 : nDstXOff = nDstYOff = nDstXSize = nDstYSize = -1;
450 : }
451 :
452 336 : return CE_None;
453 : }
454 :
455 : /************************************************************************/
456 : /* GetFileList() */
457 : /************************************************************************/
458 :
459 5 : void VRTSimpleSource::GetFileList(char*** ppapszFileList, int *pnSize,
460 : int *pnMaxSize, CPLHashSet* hSetFiles)
461 : {
462 : const char* pszFilename;
463 10 : if (poRasterBand != NULL && poRasterBand->GetDataset() != NULL &&
464 5 : (pszFilename = poRasterBand->GetDataset()->GetDescription()) != NULL)
465 : {
466 : /* -------------------------------------------------------------------- */
467 : /* Is the filename even a real filesystem object? */
468 : /* -------------------------------------------------------------------- */
469 5 : if ( strstr(pszFilename, "/vsicurl/http") != NULL ||
470 : strstr(pszFilename, "/vsicurl/ftp") != NULL )
471 : {
472 : /* Testing the existence of remote ressources can be excruciating */
473 : /* slow, so let's just suppose they exist */
474 : }
475 : else
476 : {
477 : VSIStatBufL sStat;
478 5 : if( VSIStatExL( pszFilename, &sStat, VSI_STAT_EXISTS_FLAG ) != 0 )
479 0 : return;
480 : }
481 :
482 : /* -------------------------------------------------------------------- */
483 : /* Is it already in the list ? */
484 : /* -------------------------------------------------------------------- */
485 5 : if( CPLHashSetLookup(hSetFiles, pszFilename) != NULL )
486 0 : return;
487 :
488 : /* -------------------------------------------------------------------- */
489 : /* Grow array if necessary */
490 : /* -------------------------------------------------------------------- */
491 5 : if (*pnSize + 1 >= *pnMaxSize)
492 : {
493 5 : *pnMaxSize = 2 + 2 * (*pnMaxSize);
494 : *ppapszFileList = (char **) CPLRealloc(
495 5 : *ppapszFileList, sizeof(char*) * (*pnMaxSize) );
496 : }
497 :
498 : /* -------------------------------------------------------------------- */
499 : /* Add the string to the list */
500 : /* -------------------------------------------------------------------- */
501 5 : (*ppapszFileList)[*pnSize] = CPLStrdup(pszFilename);
502 5 : (*ppapszFileList)[(*pnSize + 1)] = NULL;
503 5 : CPLHashSetInsert(hSetFiles, (*ppapszFileList)[*pnSize]);
504 :
505 5 : (*pnSize) ++;
506 : }
507 : }
508 :
509 : /************************************************************************/
510 : /* GetBand() */
511 : /************************************************************************/
512 :
513 393 : GDALRasterBand* VRTSimpleSource::GetBand()
514 : {
515 393 : return poMaskBandMainBand ? NULL : poRasterBand;
516 : }
517 :
518 : /************************************************************************/
519 : /* IsSameExceptBandNumber() */
520 : /************************************************************************/
521 :
522 31 : int VRTSimpleSource::IsSameExceptBandNumber(VRTSimpleSource* poOtherSource)
523 : {
524 : return nSrcXOff == poOtherSource->nSrcXOff &&
525 : nSrcYOff == poOtherSource->nSrcYOff &&
526 : nSrcXSize == poOtherSource->nSrcXSize &&
527 : nSrcYSize == poOtherSource->nSrcYSize &&
528 : nDstXOff == poOtherSource->nDstXOff &&
529 : nDstYOff == poOtherSource->nDstYOff &&
530 : nDstXSize == poOtherSource->nDstXSize &&
531 : nDstYSize == poOtherSource->nDstYSize &&
532 : bNoDataSet == poOtherSource->bNoDataSet &&
533 : dfNoDataValue == poOtherSource->dfNoDataValue &&
534 : GetBand() != NULL && poOtherSource->GetBand() != NULL &&
535 : GetBand()->GetDataset() != NULL &&
536 : poOtherSource->GetBand()->GetDataset() != NULL &&
537 31 : EQUAL(GetBand()->GetDataset()->GetDescription(),
538 : poOtherSource->GetBand()->GetDataset()->GetDescription());
539 : }
540 :
541 : /************************************************************************/
542 : /* SrcToDst() */
543 : /* */
544 : /* Note: this is a no-op if the dst window is -1,-1,-1,-1. */
545 : /************************************************************************/
546 :
547 10206 : void VRTSimpleSource::SrcToDst( double dfX, double dfY,
548 : double &dfXOut, double &dfYOut )
549 :
550 : {
551 10206 : dfXOut = ((dfX - nSrcXOff) / nSrcXSize) * nDstXSize + nDstXOff;
552 10206 : dfYOut = ((dfY - nSrcYOff) / nSrcYSize) * nDstYSize + nDstYOff;
553 10206 : }
554 :
555 : /************************************************************************/
556 : /* DstToSrc() */
557 : /* */
558 : /* Note: this is a no-op if the dst window is -1,-1,-1,-1. */
559 : /************************************************************************/
560 :
561 3400 : void VRTSimpleSource::DstToSrc( double dfX, double dfY,
562 : double &dfXOut, double &dfYOut )
563 :
564 : {
565 3400 : dfXOut = ((dfX - nDstXOff) / nDstXSize) * nSrcXSize + nSrcXOff;
566 3400 : dfYOut = ((dfY - nDstYOff) / nDstYSize) * nSrcYSize + nSrcYOff;
567 3400 : }
568 :
569 : /************************************************************************/
570 : /* GetSrcDstWindow() */
571 : /************************************************************************/
572 :
573 : int
574 13525 : VRTSimpleSource::GetSrcDstWindow( int nXOff, int nYOff, int nXSize, int nYSize,
575 : int nBufXSize, int nBufYSize,
576 : int *pnReqXOff, int *pnReqYOff,
577 : int *pnReqXSize, int *pnReqYSize,
578 : int *pnOutXOff, int *pnOutYOff,
579 : int *pnOutXSize, int *pnOutYSize )
580 :
581 : {
582 : int bDstWinSet = nDstXOff != -1 || nDstXSize != -1
583 13525 : || nDstYOff != -1 || nDstYSize != -1;
584 :
585 : #ifdef DEBUG
586 : int bSrcWinSet = nSrcXOff != -1 || nSrcXSize != -1
587 13525 : || nSrcYOff != -1 || nSrcYSize != -1;
588 :
589 13525 : if( bSrcWinSet != bDstWinSet )
590 : {
591 0 : return FALSE;
592 : }
593 : #endif
594 :
595 : /* -------------------------------------------------------------------- */
596 : /* If the input window completely misses the portion of the */
597 : /* virtual dataset provided by this source we have nothing to do. */
598 : /* -------------------------------------------------------------------- */
599 13525 : if( bDstWinSet )
600 : {
601 13500 : if( nXOff >= nDstXOff + nDstXSize
602 : || nYOff >= nDstYOff + nDstYSize
603 : || nXOff + nXSize < nDstXOff
604 : || nYOff + nYSize < nDstYOff )
605 400 : return FALSE;
606 : }
607 :
608 : /* -------------------------------------------------------------------- */
609 : /* This request window corresponds to the whole output buffer. */
610 : /* -------------------------------------------------------------------- */
611 13125 : *pnOutXOff = 0;
612 13125 : *pnOutYOff = 0;
613 13125 : *pnOutXSize = nBufXSize;
614 13125 : *pnOutYSize = nBufYSize;
615 :
616 : /* -------------------------------------------------------------------- */
617 : /* If the input window extents outside the portion of the on */
618 : /* the virtual file that this source can set, then clip down */
619 : /* the requested window. */
620 : /* -------------------------------------------------------------------- */
621 13125 : int bModifiedX = FALSE, bModifiedY = FALSE;
622 13125 : int nRXOff = nXOff;
623 13125 : int nRYOff = nYOff;
624 13125 : int nRXSize = nXSize;
625 13125 : int nRYSize = nYSize;
626 :
627 :
628 13125 : if( bDstWinSet )
629 : {
630 13100 : if( nRXOff < nDstXOff )
631 : {
632 2532 : nRXSize = nRXSize + nRXOff - nDstXOff;
633 2532 : nRXOff = nDstXOff;
634 2532 : bModifiedX = TRUE;
635 : }
636 :
637 13100 : if( nRYOff < nDstYOff )
638 : {
639 16 : nRYSize = nRYSize + nRYOff - nDstYOff;
640 16 : nRYOff = nDstYOff;
641 16 : bModifiedY = TRUE;
642 : }
643 :
644 13100 : if( nRXOff + nRXSize > nDstXOff + nDstXSize )
645 : {
646 2593 : nRXSize = nDstXOff + nDstXSize - nRXOff;
647 2593 : bModifiedX = TRUE;
648 : }
649 :
650 13100 : if( nRYOff + nRYSize > nDstYOff + nDstYSize )
651 : {
652 0 : nRYSize = nDstYOff + nDstYSize - nRYOff;
653 0 : bModifiedY = TRUE;
654 : }
655 : }
656 :
657 : /* -------------------------------------------------------------------- */
658 : /* Translate requested region in virtual file into the source */
659 : /* band coordinates. */
660 : /* -------------------------------------------------------------------- */
661 13125 : double dfScaleX = nSrcXSize / (double) nDstXSize;
662 13125 : double dfScaleY = nSrcYSize / (double) nDstYSize;
663 :
664 13125 : *pnReqXOff = (int) floor((nRXOff - nDstXOff) * dfScaleX + nSrcXOff);
665 13125 : *pnReqYOff = (int) floor((nRYOff - nDstYOff) * dfScaleY + nSrcYOff);
666 :
667 13125 : *pnReqXSize = (int) floor(nRXSize * dfScaleX + 0.5);
668 13125 : *pnReqYSize = (int) floor(nRYSize * dfScaleY + 0.5);
669 :
670 : /* -------------------------------------------------------------------- */
671 : /* Clamp within the bounds of the available source data. */
672 : /* -------------------------------------------------------------------- */
673 13125 : if( *pnReqXOff < 0 )
674 : {
675 0 : *pnReqXSize += *pnReqXOff;
676 0 : *pnReqXOff = 0;
677 :
678 0 : bModifiedX = TRUE;
679 : }
680 :
681 13125 : if( *pnReqYOff < 0 )
682 : {
683 0 : *pnReqYSize += *pnReqYOff;
684 0 : *pnReqYOff = 0;
685 0 : bModifiedY = TRUE;
686 : }
687 :
688 13125 : if( *pnReqXSize == 0 )
689 0 : *pnReqXSize = 1;
690 13125 : if( *pnReqYSize == 0 )
691 16 : *pnReqYSize = 1;
692 :
693 13125 : if( *pnReqXOff + *pnReqXSize > poRasterBand->GetXSize() )
694 : {
695 4 : *pnReqXSize = poRasterBand->GetXSize() - *pnReqXOff;
696 4 : bModifiedX = TRUE;
697 : }
698 :
699 13125 : if( *pnReqYOff + *pnReqYSize > poRasterBand->GetYSize() )
700 : {
701 4 : *pnReqYSize = poRasterBand->GetYSize() - *pnReqYOff;
702 4 : bModifiedY = TRUE;
703 : }
704 :
705 : /* -------------------------------------------------------------------- */
706 : /* Don't do anything if the requesting region is completely off */
707 : /* the source image. */
708 : /* -------------------------------------------------------------------- */
709 13125 : if( *pnReqXOff >= poRasterBand->GetXSize()
710 : || *pnReqYOff >= poRasterBand->GetYSize()
711 : || *pnReqXSize <= 0 || *pnReqYSize <= 0 )
712 : {
713 0 : return FALSE;
714 : }
715 :
716 : /* -------------------------------------------------------------------- */
717 : /* If we haven't had to modify the source rectangle, then the */
718 : /* destination rectangle must be the whole region. */
719 : /* -------------------------------------------------------------------- */
720 13125 : if( !bModifiedX && !bModifiedY )
721 8022 : return TRUE;
722 :
723 : /* -------------------------------------------------------------------- */
724 : /* Now transform this possibly reduced request back into the */
725 : /* destination buffer coordinates in case the output region is */
726 : /* less than the whole buffer. */
727 : /* -------------------------------------------------------------------- */
728 : double dfDstULX, dfDstULY, dfDstLRX, dfDstLRY;
729 : double dfScaleWinToBufX, dfScaleWinToBufY;
730 :
731 5103 : SrcToDst( (double) *pnReqXOff, (double) *pnReqYOff, dfDstULX, dfDstULY );
732 : SrcToDst( *pnReqXOff + *pnReqXSize, *pnReqYOff + *pnReqYSize,
733 5103 : dfDstLRX, dfDstLRY );
734 :
735 5103 : if( bModifiedX )
736 : {
737 5103 : dfScaleWinToBufX = nBufXSize / (double) nXSize;
738 :
739 5103 : *pnOutXOff = (int) ((dfDstULX - nXOff) * dfScaleWinToBufX+0.001);
740 : *pnOutXSize = (int) ((dfDstLRX - nXOff) * dfScaleWinToBufX+0.001)
741 5103 : - *pnOutXOff;
742 :
743 5103 : *pnOutXOff = MAX(0,*pnOutXOff);
744 5103 : if( *pnOutXOff + *pnOutXSize > nBufXSize )
745 0 : *pnOutXSize = nBufXSize - *pnOutXOff;
746 : }
747 :
748 5103 : if( bModifiedY )
749 : {
750 20 : dfScaleWinToBufY = nBufYSize / (double) nYSize;
751 :
752 20 : *pnOutYOff = (int) ((dfDstULY - nYOff) * dfScaleWinToBufY+0.001);
753 : *pnOutYSize = (int) ((dfDstLRY - nYOff) * dfScaleWinToBufY+0.001)
754 20 : - *pnOutYOff;
755 :
756 20 : *pnOutYOff = MAX(0,*pnOutYOff);
757 20 : if( *pnOutYOff + *pnOutYSize > nBufYSize )
758 16 : *pnOutYSize = nBufYSize - *pnOutYOff;
759 : }
760 :
761 5103 : if( *pnOutXSize < 1 || *pnOutYSize < 1 )
762 16 : return FALSE;
763 : else
764 5087 : return TRUE;
765 : }
766 :
767 : /************************************************************************/
768 : /* RasterIO() */
769 : /************************************************************************/
770 :
771 : CPLErr
772 11083 : VRTSimpleSource::RasterIO( int nXOff, int nYOff, int nXSize, int nYSize,
773 : void *pData, int nBufXSize, int nBufYSize,
774 : GDALDataType eBufType,
775 : int nPixelSpace, int nLineSpace )
776 :
777 : {
778 : // The window we will actually request from the source raster band.
779 : int nReqXOff, nReqYOff, nReqXSize, nReqYSize;
780 :
781 : // The window we will actual set _within_ the pData buffer.
782 : int nOutXOff, nOutYOff, nOutXSize, nOutYSize;
783 :
784 11083 : if( !GetSrcDstWindow( nXOff, nYOff, nXSize, nYSize,
785 : nBufXSize, nBufYSize,
786 : &nReqXOff, &nReqYOff, &nReqXSize, &nReqYSize,
787 : &nOutXOff, &nOutYOff, &nOutXSize, &nOutYSize ) )
788 : {
789 416 : return CE_None;
790 : }
791 :
792 : /* -------------------------------------------------------------------- */
793 : /* Actually perform the IO request. */
794 : /* -------------------------------------------------------------------- */
795 : CPLErr eErr;
796 :
797 : eErr =
798 : poRasterBand->RasterIO( GF_Read,
799 : nReqXOff, nReqYOff, nReqXSize, nReqYSize,
800 : ((unsigned char *) pData)
801 : + nOutXOff * nPixelSpace
802 : + nOutYOff * nLineSpace,
803 : nOutXSize, nOutYSize,
804 10667 : eBufType, nPixelSpace, nLineSpace );
805 :
806 10667 : return eErr;
807 : }
808 :
809 : /************************************************************************/
810 : /* GetMinimum() */
811 : /************************************************************************/
812 :
813 9 : double VRTSimpleSource::GetMinimum( int nXSize, int nYSize, int *pbSuccess )
814 : {
815 : // The window we will actually request from the source raster band.
816 : int nReqXOff, nReqYOff, nReqXSize, nReqYSize;
817 :
818 : // The window we will actual set _within_ the pData buffer.
819 : int nOutXOff, nOutYOff, nOutXSize, nOutYSize;
820 :
821 9 : if( !GetSrcDstWindow( 0, 0, nXSize, nYSize,
822 : nXSize, nYSize,
823 : &nReqXOff, &nReqYOff, &nReqXSize, &nReqYSize,
824 : &nOutXOff, &nOutYOff, &nOutXSize, &nOutYSize ) ||
825 : nReqXOff != 0 || nReqYOff != 0 ||
826 : nReqXSize != poRasterBand->GetXSize() ||
827 : nReqYSize != poRasterBand->GetYSize())
828 : {
829 0 : *pbSuccess = FALSE;
830 0 : return 0;
831 : }
832 :
833 9 : return poRasterBand->GetMinimum(pbSuccess);
834 : }
835 :
836 : /************************************************************************/
837 : /* GetMaximum() */
838 : /************************************************************************/
839 :
840 9 : double VRTSimpleSource::GetMaximum( int nXSize, int nYSize, int *pbSuccess )
841 : {
842 : // The window we will actually request from the source raster band.
843 : int nReqXOff, nReqYOff, nReqXSize, nReqYSize;
844 :
845 : // The window we will actual set _within_ the pData buffer.
846 : int nOutXOff, nOutYOff, nOutXSize, nOutYSize;
847 :
848 9 : if( !GetSrcDstWindow( 0, 0, nXSize, nYSize,
849 : nXSize, nYSize,
850 : &nReqXOff, &nReqYOff, &nReqXSize, &nReqYSize,
851 : &nOutXOff, &nOutYOff, &nOutXSize, &nOutYSize ) ||
852 : nReqXOff != 0 || nReqYOff != 0 ||
853 : nReqXSize != poRasterBand->GetXSize() ||
854 : nReqYSize != poRasterBand->GetYSize())
855 : {
856 0 : *pbSuccess = FALSE;
857 0 : return 0;
858 : }
859 :
860 9 : return poRasterBand->GetMaximum(pbSuccess);
861 : }
862 :
863 :
864 : /************************************************************************/
865 : /* DatasetRasterIO() */
866 : /************************************************************************/
867 :
868 56 : CPLErr VRTSimpleSource::DatasetRasterIO(
869 : int nXOff, int nYOff, int nXSize, int nYSize,
870 : void * pData, int nBufXSize, int nBufYSize,
871 : GDALDataType eBufType,
872 : int nBandCount, int *panBandMap,
873 : int nPixelSpace, int nLineSpace, int nBandSpace)
874 : {
875 56 : if (!EQUAL(GetType(), "SimpleSource"))
876 : {
877 : CPLError(CE_Failure, CPLE_NotSupported,
878 0 : "DatasetRasterIO() not implemented for %s", GetType());
879 0 : return CE_Failure;
880 : }
881 :
882 : // The window we will actually request from the source raster band.
883 : int nReqXOff, nReqYOff, nReqXSize, nReqYSize;
884 :
885 : // The window we will actual set _within_ the pData buffer.
886 : int nOutXOff, nOutYOff, nOutXSize, nOutYSize;
887 :
888 56 : if( !GetSrcDstWindow( nXOff, nYOff, nXSize, nYSize,
889 : nBufXSize, nBufYSize,
890 : &nReqXOff, &nReqYOff, &nReqXSize, &nReqYSize,
891 : &nOutXOff, &nOutYOff, &nOutXSize, &nOutYSize ) )
892 : {
893 0 : return CE_None;
894 : }
895 :
896 56 : GDALDataset* poDS = poRasterBand->GetDataset();
897 56 : if (poDS == NULL)
898 0 : return CE_Failure;
899 :
900 : return poDS->RasterIO( GF_Read,
901 : nReqXOff, nReqYOff, nReqXSize, nReqYSize,
902 : ((unsigned char *) pData)
903 : + nOutXOff * nPixelSpace
904 : + nOutYOff * nLineSpace,
905 : nOutXSize, nOutYSize,
906 : eBufType, nBandCount, panBandMap,
907 56 : nPixelSpace, nLineSpace, nBandSpace );
908 : }
909 :
910 : /************************************************************************/
911 : /* ==================================================================== */
912 : /* VRTAveragedSource */
913 : /* ==================================================================== */
914 : /************************************************************************/
915 :
916 : /************************************************************************/
917 : /* VRTAveragedSource() */
918 : /************************************************************************/
919 :
920 2 : VRTAveragedSource::VRTAveragedSource()
921 : {
922 2 : }
923 :
924 : /************************************************************************/
925 : /* SerializeToXML() */
926 : /************************************************************************/
927 :
928 0 : CPLXMLNode *VRTAveragedSource::SerializeToXML( const char *pszVRTPath )
929 :
930 : {
931 0 : CPLXMLNode *psSrc = VRTSimpleSource::SerializeToXML( pszVRTPath );
932 :
933 0 : if( psSrc == NULL )
934 0 : return NULL;
935 :
936 0 : CPLFree( psSrc->pszValue );
937 0 : psSrc->pszValue = CPLStrdup( "AveragedSource" );
938 :
939 0 : return psSrc;
940 : }
941 :
942 : /************************************************************************/
943 : /* RasterIO() */
944 : /************************************************************************/
945 :
946 : CPLErr
947 50 : VRTAveragedSource::RasterIO( int nXOff, int nYOff, int nXSize, int nYSize,
948 : void *pData, int nBufXSize, int nBufYSize,
949 : GDALDataType eBufType,
950 : int nPixelSpace, int nLineSpace )
951 :
952 : {
953 : // The window we will actually request from the source raster band.
954 : int nReqXOff, nReqYOff, nReqXSize, nReqYSize;
955 :
956 : // The window we will actual set _within_ the pData buffer.
957 : int nOutXOff, nOutYOff, nOutXSize, nOutYSize;
958 :
959 50 : if( !GetSrcDstWindow( nXOff, nYOff, nXSize, nYSize, nBufXSize, nBufYSize,
960 : &nReqXOff, &nReqYOff, &nReqXSize, &nReqYSize,
961 : &nOutXOff, &nOutYOff, &nOutXSize, &nOutYSize ) )
962 0 : return CE_None;
963 :
964 : /* -------------------------------------------------------------------- */
965 : /* Allocate a temporary buffer to whole the full resolution */
966 : /* data from the area of interest. */
967 : /* -------------------------------------------------------------------- */
968 : float *pafSrc;
969 :
970 50 : pafSrc = (float *) VSIMalloc3(sizeof(float), nReqXSize, nReqYSize);
971 50 : if( pafSrc == NULL )
972 : {
973 : CPLError( CE_Failure, CPLE_OutOfMemory,
974 0 : "Out of memory allocating working buffer in VRTAveragedSource::RasterIO()." );
975 0 : return CE_Failure;
976 : }
977 :
978 : /* -------------------------------------------------------------------- */
979 : /* Load it. */
980 : /* -------------------------------------------------------------------- */
981 : CPLErr eErr;
982 :
983 : eErr = poRasterBand->RasterIO( GF_Read,
984 : nReqXOff, nReqYOff, nReqXSize, nReqYSize,
985 : pafSrc, nReqXSize, nReqYSize, GDT_Float32,
986 50 : 0, 0 );
987 :
988 50 : if( eErr != CE_None )
989 : {
990 0 : VSIFree( pafSrc );
991 0 : return eErr;
992 : }
993 :
994 : /* -------------------------------------------------------------------- */
995 : /* Do the averaging. */
996 : /* -------------------------------------------------------------------- */
997 100 : for( int iBufLine = nOutYOff; iBufLine < nOutYOff + nOutYSize; iBufLine++ )
998 : {
999 : double dfYDst;
1000 :
1001 50 : dfYDst = (iBufLine / (double) nBufYSize) * nYSize + nYOff;
1002 :
1003 1750 : for( int iBufPixel = nOutXOff;
1004 : iBufPixel < nOutXOff + nOutXSize;
1005 : iBufPixel++ )
1006 : {
1007 : double dfXDst;
1008 : double dfXSrcStart, dfXSrcEnd, dfYSrcStart, dfYSrcEnd;
1009 : int iXSrcStart, iYSrcStart, iXSrcEnd, iYSrcEnd;
1010 :
1011 1700 : dfXDst = (iBufPixel / (double) nBufXSize) * nXSize + nXOff;
1012 :
1013 : // Compute the source image rectangle needed for this pixel.
1014 1700 : DstToSrc( dfXDst, dfYDst, dfXSrcStart, dfYSrcStart );
1015 1700 : DstToSrc( dfXDst+1.0, dfYDst+1.0, dfXSrcEnd, dfYSrcEnd );
1016 :
1017 : // Convert to integers, assuming that the center of the source
1018 : // pixel must be in our rect to get included.
1019 1700 : if (dfXSrcEnd >= dfXSrcStart + 1)
1020 : {
1021 100 : iXSrcStart = (int) floor(dfXSrcStart+0.5);
1022 100 : iXSrcEnd = (int) floor(dfXSrcEnd+0.5);
1023 : }
1024 : else
1025 : {
1026 : /* If the resampling factor is less than 100%, the distance */
1027 : /* between the source pixel is < 1, so we stick to nearest */
1028 : /* neighbour */
1029 1600 : iXSrcStart = (int) floor(dfXSrcStart);
1030 1600 : iXSrcEnd = iXSrcStart + 1;
1031 : }
1032 1700 : if (dfYSrcEnd >= dfYSrcStart + 1)
1033 : {
1034 100 : iYSrcStart = (int) floor(dfYSrcStart+0.5);
1035 100 : iYSrcEnd = (int) floor(dfYSrcEnd+0.5);
1036 : }
1037 : else
1038 : {
1039 1600 : iYSrcStart = (int) floor(dfYSrcStart);
1040 1600 : iYSrcEnd = iYSrcStart + 1;
1041 : }
1042 :
1043 : // Transform into the coordinate system of the source *buffer*
1044 1700 : iXSrcStart -= nReqXOff;
1045 1700 : iYSrcStart -= nReqYOff;
1046 1700 : iXSrcEnd -= nReqXOff;
1047 1700 : iYSrcEnd -= nReqYOff;
1048 :
1049 1700 : double dfSum = 0.0;
1050 1700 : int nPixelCount = 0;
1051 :
1052 3500 : for( int iY = iYSrcStart; iY < iYSrcEnd; iY++ )
1053 : {
1054 1800 : if( iY < 0 || iY >= nReqYSize )
1055 0 : continue;
1056 :
1057 3800 : for( int iX = iXSrcStart; iX < iXSrcEnd; iX++ )
1058 : {
1059 2000 : if( iX < 0 || iX >= nReqXSize )
1060 0 : continue;
1061 :
1062 2000 : float fSampledValue = pafSrc[iX + iY * nReqXSize];
1063 2000 : if (CPLIsNan(fSampledValue))
1064 0 : continue;
1065 :
1066 2000 : if( bNoDataSet && ARE_REAL_EQUAL(fSampledValue, dfNoDataValue))
1067 0 : continue;
1068 :
1069 2000 : nPixelCount++;
1070 2000 : dfSum += pafSrc[iX + iY * nReqXSize];
1071 : }
1072 : }
1073 :
1074 1700 : if( nPixelCount == 0 )
1075 0 : continue;
1076 :
1077 : // Compute output value.
1078 1700 : float dfOutputValue = (float) (dfSum / nPixelCount);
1079 :
1080 : // Put it in the output buffer.
1081 : GByte *pDstLocation;
1082 :
1083 : pDstLocation = ((GByte *)pData)
1084 : + nPixelSpace * iBufPixel
1085 1700 : + nLineSpace * iBufLine;
1086 :
1087 1700 : if( eBufType == GDT_Byte )
1088 0 : *pDstLocation = (GByte) MIN(255,MAX(0,dfOutputValue + 0.5));
1089 : else
1090 : GDALCopyWords( &dfOutputValue, GDT_Float32, 4,
1091 1700 : pDstLocation, eBufType, 8, 1 );
1092 : }
1093 : }
1094 :
1095 50 : VSIFree( pafSrc );
1096 :
1097 50 : return CE_None;
1098 : }
1099 :
1100 : /************************************************************************/
1101 : /* GetMinimum() */
1102 : /************************************************************************/
1103 :
1104 0 : double VRTAveragedSource::GetMinimum( int nXSize, int nYSize, int *pbSuccess )
1105 : {
1106 0 : *pbSuccess = FALSE;
1107 0 : return 0;
1108 : }
1109 :
1110 : /************************************************************************/
1111 : /* GetMaximum() */
1112 : /************************************************************************/
1113 :
1114 0 : double VRTAveragedSource::GetMaximum( int nXSize, int nYSize, int *pbSuccess )
1115 : {
1116 0 : *pbSuccess = FALSE;
1117 0 : return 0;
1118 : }
1119 :
1120 : /************************************************************************/
1121 : /* ==================================================================== */
1122 : /* VRTComplexSource */
1123 : /* ==================================================================== */
1124 : /************************************************************************/
1125 :
1126 : /************************************************************************/
1127 : /* VRTComplexSource() */
1128 : /************************************************************************/
1129 :
1130 61 : VRTComplexSource::VRTComplexSource()
1131 :
1132 : {
1133 61 : bDoScaling = FALSE;
1134 61 : dfScaleOff = 0.0;
1135 61 : dfScaleRatio = 1.0;
1136 :
1137 61 : bNoDataSet = FALSE;
1138 61 : dfNoDataValue = 0.0;
1139 :
1140 61 : padfLUTInputs = NULL;
1141 61 : padfLUTOutputs = NULL;
1142 61 : nLUTItemCount = 0;
1143 61 : nColorTableComponent = 0;
1144 61 : }
1145 :
1146 61 : VRTComplexSource::~VRTComplexSource()
1147 : {
1148 61 : if (padfLUTInputs)
1149 7 : VSIFree( padfLUTInputs );
1150 61 : if (padfLUTOutputs)
1151 7 : VSIFree( padfLUTOutputs );
1152 61 : }
1153 :
1154 : /************************************************************************/
1155 : /* SerializeToXML() */
1156 : /************************************************************************/
1157 :
1158 15 : CPLXMLNode *VRTComplexSource::SerializeToXML( const char *pszVRTPath )
1159 :
1160 : {
1161 15 : CPLXMLNode *psSrc = VRTSimpleSource::SerializeToXML( pszVRTPath );
1162 :
1163 15 : if( psSrc == NULL )
1164 0 : return NULL;
1165 :
1166 15 : CPLFree( psSrc->pszValue );
1167 15 : psSrc->pszValue = CPLStrdup( "ComplexSource" );
1168 :
1169 15 : if( bNoDataSet )
1170 : {
1171 8 : if (CPLIsNan(dfNoDataValue))
1172 0 : CPLSetXMLValue( psSrc, "NODATA", "nan");
1173 : else
1174 : CPLSetXMLValue( psSrc, "NODATA",
1175 8 : CPLSPrintf("%g", dfNoDataValue) );
1176 : }
1177 :
1178 15 : if( bDoScaling )
1179 : {
1180 : CPLSetXMLValue( psSrc, "ScaleOffset",
1181 0 : CPLSPrintf("%g", dfScaleOff) );
1182 : CPLSetXMLValue( psSrc, "ScaleRatio",
1183 0 : CPLSPrintf("%g", dfScaleRatio) );
1184 : }
1185 :
1186 15 : if ( nLUTItemCount )
1187 : {
1188 0 : CPLString osLUT = CPLString().Printf("%g:%g", padfLUTInputs[0], padfLUTOutputs[0]);
1189 : int i;
1190 0 : for ( i = 1; i < nLUTItemCount; i++ )
1191 0 : osLUT += CPLString().Printf(",%g:%g", padfLUTInputs[i], padfLUTOutputs[i]);
1192 0 : CPLSetXMLValue( psSrc, "LUT", osLUT );
1193 : }
1194 :
1195 15 : if ( nColorTableComponent )
1196 : {
1197 : CPLSetXMLValue( psSrc, "ColorTableComponent",
1198 7 : CPLSPrintf("%d", nColorTableComponent) );
1199 : }
1200 :
1201 15 : return psSrc;
1202 : }
1203 :
1204 : /************************************************************************/
1205 : /* XMLInit() */
1206 : /************************************************************************/
1207 :
1208 37 : CPLErr VRTComplexSource::XMLInit( CPLXMLNode *psSrc, const char *pszVRTPath )
1209 :
1210 : {
1211 : CPLErr eErr;
1212 :
1213 : /* -------------------------------------------------------------------- */
1214 : /* Do base initialization. */
1215 : /* -------------------------------------------------------------------- */
1216 37 : eErr = VRTSimpleSource::XMLInit( psSrc, pszVRTPath );
1217 37 : if( eErr != CE_None )
1218 0 : return eErr;
1219 :
1220 : /* -------------------------------------------------------------------- */
1221 : /* Complex parameters. */
1222 : /* -------------------------------------------------------------------- */
1223 37 : if( CPLGetXMLValue(psSrc, "ScaleOffset", NULL) != NULL
1224 : || CPLGetXMLValue(psSrc, "ScaleRatio", NULL) != NULL )
1225 : {
1226 1 : bDoScaling = TRUE;
1227 1 : dfScaleOff = atof(CPLGetXMLValue(psSrc, "ScaleOffset", "0") );
1228 1 : dfScaleRatio = atof(CPLGetXMLValue(psSrc, "ScaleRatio", "1") );
1229 : }
1230 :
1231 37 : if( CPLGetXMLValue(psSrc, "NODATA", NULL) != NULL )
1232 : {
1233 10 : bNoDataSet = TRUE;
1234 10 : dfNoDataValue = CPLAtofM(CPLGetXMLValue(psSrc, "NODATA", "0"));
1235 : }
1236 :
1237 37 : if( CPLGetXMLValue(psSrc, "LUT", NULL) != NULL )
1238 : {
1239 : int nIndex;
1240 7 : char **papszValues = CSLTokenizeString2(CPLGetXMLValue(psSrc, "LUT", ""), ",:", CSLT_ALLOWEMPTYTOKENS);
1241 :
1242 7 : if (nLUTItemCount)
1243 : {
1244 0 : if (padfLUTInputs)
1245 : {
1246 0 : VSIFree( padfLUTInputs );
1247 0 : padfLUTInputs = NULL;
1248 : }
1249 0 : if (padfLUTOutputs)
1250 : {
1251 0 : VSIFree( padfLUTOutputs );
1252 0 : padfLUTOutputs = NULL;
1253 : }
1254 0 : nLUTItemCount = 0;
1255 : }
1256 :
1257 7 : nLUTItemCount = CSLCount(papszValues) / 2;
1258 :
1259 7 : padfLUTInputs = (double *) VSIMalloc2(nLUTItemCount, sizeof(double));
1260 7 : if ( !padfLUTInputs )
1261 : {
1262 0 : CSLDestroy(papszValues);
1263 0 : nLUTItemCount = 0;
1264 0 : return CE_Failure;
1265 : }
1266 :
1267 7 : padfLUTOutputs = (double *) VSIMalloc2(nLUTItemCount, sizeof(double));
1268 7 : if ( !padfLUTOutputs )
1269 : {
1270 0 : CSLDestroy(papszValues);
1271 0 : VSIFree( padfLUTInputs );
1272 0 : padfLUTInputs = NULL;
1273 0 : nLUTItemCount = 0;
1274 0 : return CE_Failure;
1275 : }
1276 :
1277 71 : for ( nIndex = 0; nIndex < nLUTItemCount; nIndex++ )
1278 : {
1279 64 : padfLUTInputs[nIndex] = atof( papszValues[nIndex * 2] );
1280 64 : padfLUTOutputs[nIndex] = atof( papszValues[nIndex * 2 + 1] );
1281 :
1282 : // Enforce the requirement that the LUT input array is monotonically non-decreasing.
1283 64 : if ( nIndex > 0 && padfLUTInputs[nIndex] < padfLUTInputs[nIndex - 1] )
1284 : {
1285 0 : CSLDestroy(papszValues);
1286 0 : VSIFree( padfLUTInputs );
1287 0 : VSIFree( padfLUTOutputs );
1288 0 : padfLUTInputs = NULL;
1289 0 : padfLUTOutputs = NULL;
1290 0 : nLUTItemCount = 0;
1291 0 : return CE_Failure;
1292 : }
1293 : }
1294 :
1295 7 : CSLDestroy(papszValues);
1296 : }
1297 :
1298 37 : if( CPLGetXMLValue(psSrc, "ColorTableComponent", NULL) != NULL )
1299 : {
1300 14 : nColorTableComponent = atoi(CPLGetXMLValue(psSrc, "ColorTableComponent", "0"));
1301 : }
1302 :
1303 37 : return CE_None;
1304 : }
1305 :
1306 : /************************************************************************/
1307 : /* LookupValue() */
1308 : /************************************************************************/
1309 :
1310 : double
1311 58964 : VRTComplexSource::LookupValue( double dfInput )
1312 : {
1313 : // Find the index of the first element in the LUT input array that
1314 : // is not smaller than the input value.
1315 58964 : int i = std::lower_bound(padfLUTInputs, padfLUTInputs + nLUTItemCount, dfInput) - padfLUTInputs;
1316 :
1317 58964 : if (i == 0)
1318 0 : return padfLUTOutputs[0];
1319 :
1320 : // If the index is beyond the end of the LUT input array, the input
1321 : // value is larger than all the values in the array.
1322 58964 : if (i == nLUTItemCount)
1323 4 : return padfLUTOutputs[nLUTItemCount - 1];
1324 :
1325 58960 : if (padfLUTInputs[i] == dfInput)
1326 4859 : return padfLUTOutputs[i];
1327 :
1328 : // Otherwise, interpolate.
1329 108202 : return padfLUTOutputs[i - 1] + (dfInput - padfLUTInputs[i - 1]) *
1330 108202 : ((padfLUTOutputs[i] - padfLUTOutputs[i - 1]) / (padfLUTInputs[i] - padfLUTInputs[i - 1]));
1331 : }
1332 :
1333 : /************************************************************************/
1334 : /* RasterIO() */
1335 : /************************************************************************/
1336 :
1337 : CPLErr
1338 2317 : VRTComplexSource::RasterIO( int nXOff, int nYOff, int nXSize, int nYSize,
1339 : void *pData, int nBufXSize, int nBufYSize,
1340 : GDALDataType eBufType,
1341 : int nPixelSpace, int nLineSpace )
1342 :
1343 : {
1344 : // The window we will actually request from the source raster band.
1345 : int nReqXOff, nReqYOff, nReqXSize, nReqYSize;
1346 :
1347 : // The window we will actual set _within_ the pData buffer.
1348 : int nOutXOff, nOutYOff, nOutXSize, nOutYSize;
1349 :
1350 2317 : if( !GetSrcDstWindow( nXOff, nYOff, nXSize, nYSize, nBufXSize, nBufYSize,
1351 : &nReqXOff, &nReqYOff, &nReqXSize, &nReqYSize,
1352 : &nOutXOff, &nOutYOff, &nOutXSize, &nOutYSize ) )
1353 0 : return CE_None;
1354 :
1355 :
1356 : /* -------------------------------------------------------------------- */
1357 : /* Read into a temporary buffer. */
1358 : /* -------------------------------------------------------------------- */
1359 : float *pafData;
1360 : CPLErr eErr;
1361 2317 : GDALColorTable* poColorTable = NULL;
1362 2317 : int bIsComplex = GDALDataTypeIsComplex(eBufType);
1363 2317 : GDALDataType eWrkDataType = (bIsComplex) ? GDT_CFloat32 : GDT_Float32;
1364 2317 : int nWordSize = GDALGetDataTypeSize(eWrkDataType) / 8;
1365 2317 : int bNoDataSetAndNotNan = bNoDataSet && !CPLIsNan(dfNoDataValue);
1366 :
1367 2317 : if( bDoScaling && bNoDataSet == FALSE && dfScaleRatio == 0)
1368 : {
1369 : /* -------------------------------------------------------------------- */
1370 : /* Optimization when outputing a constant value */
1371 : /* (used by the -addalpha option of gdalbuildvrt) */
1372 : /* -------------------------------------------------------------------- */
1373 0 : pafData = NULL;
1374 : }
1375 : else
1376 : {
1377 2317 : pafData = (float *) VSIMalloc3(nOutXSize,nOutYSize,nWordSize);
1378 2317 : if (pafData == NULL)
1379 : {
1380 0 : CPLError( CE_Failure, CPLE_OutOfMemory, "Out of memory");
1381 0 : return CE_Failure;
1382 : }
1383 : eErr = poRasterBand->RasterIO( GF_Read,
1384 : nReqXOff, nReqYOff, nReqXSize, nReqYSize,
1385 : pafData, nOutXSize, nOutYSize, eWrkDataType,
1386 2317 : nWordSize, nWordSize * nOutXSize );
1387 2317 : if( eErr != CE_None )
1388 : {
1389 0 : CPLFree( pafData );
1390 0 : return eErr;
1391 : }
1392 :
1393 2317 : if (nColorTableComponent != 0)
1394 : {
1395 1610 : poColorTable = poRasterBand->GetColorTable();
1396 1610 : if (poColorTable == NULL)
1397 : {
1398 : CPLError(CE_Failure, CPLE_AppDefined,
1399 0 : "Source band has no color table.");
1400 0 : CPLFree( pafData );
1401 0 : return CE_Failure;
1402 : }
1403 : }
1404 : }
1405 :
1406 : /* -------------------------------------------------------------------- */
1407 : /* Selectively copy into output buffer with nodata masking, */
1408 : /* and/or scaling. */
1409 : /* -------------------------------------------------------------------- */
1410 : int iX, iY;
1411 :
1412 6867 : for( iY = 0; iY < nOutYSize; iY++ )
1413 : {
1414 1381656 : for( iX = 0; iX < nOutXSize; iX++ )
1415 : {
1416 : GByte *pDstLocation;
1417 :
1418 : pDstLocation = ((GByte *)pData)
1419 : + nPixelSpace * (iX + nOutXOff)
1420 1377106 : + nLineSpace * (iY + nOutYOff);
1421 :
1422 2753409 : if (pafData && !bIsComplex)
1423 : {
1424 1377105 : float fResult = pafData[iX + iY * nOutXSize];
1425 1377105 : if( CPLIsNan(dfNoDataValue) && CPLIsNan(fResult) )
1426 20 : continue;
1427 1377085 : if( bNoDataSetAndNotNan && ARE_REAL_EQUAL(fResult, dfNoDataValue) )
1428 782 : continue;
1429 :
1430 1376303 : if (nColorTableComponent)
1431 : {
1432 1275780 : const GDALColorEntry* poEntry = poColorTable->GetColorEntry((int)fResult);
1433 1275780 : if (poEntry)
1434 : {
1435 1275780 : if (nColorTableComponent == 1)
1436 474978 : fResult = poEntry->c1;
1437 800802 : else if (nColorTableComponent == 2)
1438 320401 : fResult = poEntry->c2;
1439 480401 : else if (nColorTableComponent == 3)
1440 320401 : fResult = poEntry->c3;
1441 160000 : else if (nColorTableComponent == 4)
1442 160000 : fResult = poEntry->c4;
1443 : }
1444 : else
1445 : {
1446 : static int bHasWarned = FALSE;
1447 0 : if (!bHasWarned)
1448 : {
1449 0 : bHasWarned = TRUE;
1450 : CPLError(CE_Failure, CPLE_AppDefined,
1451 0 : "No entry %d.", (int)fResult);
1452 : }
1453 0 : continue;
1454 : }
1455 : }
1456 :
1457 1376303 : if( bDoScaling )
1458 10000 : fResult = (float) (fResult * dfScaleRatio + dfScaleOff);
1459 :
1460 1376303 : if (nLUTItemCount)
1461 58964 : fResult = (float) LookupValue( fResult );
1462 :
1463 1376303 : if( eBufType == GDT_Byte )
1464 660422 : *pDstLocation = (GByte) MIN(255,MAX(0,fResult + 0.5));
1465 : else
1466 : GDALCopyWords( &fResult, GDT_Float32, 0,
1467 715881 : pDstLocation, eBufType, 0, 1 );
1468 : }
1469 2 : else if (pafData && bIsComplex)
1470 : {
1471 : float afResult[2];
1472 1 : afResult[0] = pafData[2 * (iX + iY * nOutXSize)];
1473 1 : afResult[1] = pafData[2 * (iX + iY * nOutXSize) + 1];
1474 :
1475 : /* Do not use color table */
1476 :
1477 1 : if( bDoScaling )
1478 : {
1479 1 : afResult[0] = (float) (afResult[0] * dfScaleRatio + dfScaleOff);
1480 1 : afResult[1] = (float) (afResult[1] * dfScaleRatio + dfScaleOff);
1481 : }
1482 :
1483 : /* Do not use LUT */
1484 :
1485 1 : if( eBufType == GDT_Byte )
1486 0 : *pDstLocation = (GByte) MIN(255,MAX(0,afResult[0] + 0.5));
1487 : else
1488 : GDALCopyWords( afResult, GDT_CFloat32, 0,
1489 1 : pDstLocation, eBufType, 0, 1 );
1490 : }
1491 : else
1492 : {
1493 0 : float fResult = (float) dfScaleOff;
1494 :
1495 0 : if (nLUTItemCount)
1496 0 : fResult = (float) LookupValue( fResult );
1497 :
1498 0 : if( eBufType == GDT_Byte )
1499 0 : *pDstLocation = (GByte) MIN(255,MAX(0,fResult + 0.5));
1500 : else
1501 : GDALCopyWords( &fResult, GDT_Float32, 0,
1502 0 : pDstLocation, eBufType, 0, 1 );
1503 : }
1504 :
1505 : }
1506 : }
1507 :
1508 2317 : CPLFree( pafData );
1509 :
1510 2317 : return CE_None;
1511 : }
1512 :
1513 : /************************************************************************/
1514 : /* GetMinimum() */
1515 : /************************************************************************/
1516 :
1517 0 : double VRTComplexSource::GetMinimum( int nXSize, int nYSize, int *pbSuccess )
1518 : {
1519 0 : if (dfScaleOff == 0.0 && dfScaleRatio == 1.0 &&
1520 : nLUTItemCount == 0 && nColorTableComponent == 0)
1521 : {
1522 0 : return VRTSimpleSource::GetMinimum(nXSize, nYSize, pbSuccess);
1523 : }
1524 :
1525 0 : *pbSuccess = FALSE;
1526 0 : return 0;
1527 : }
1528 :
1529 : /************************************************************************/
1530 : /* GetMaximum() */
1531 : /************************************************************************/
1532 :
1533 0 : double VRTComplexSource::GetMaximum( int nXSize, int nYSize, int *pbSuccess )
1534 : {
1535 0 : if (dfScaleOff == 0.0 && dfScaleRatio == 1.0 &&
1536 : nLUTItemCount == 0 && nColorTableComponent == 0)
1537 : {
1538 0 : return VRTSimpleSource::GetMaximum(nXSize, nYSize, pbSuccess);
1539 : }
1540 :
1541 0 : *pbSuccess = FALSE;
1542 0 : return 0;
1543 : }
1544 :
1545 : /************************************************************************/
1546 : /* ==================================================================== */
1547 : /* VRTFuncSource */
1548 : /* ==================================================================== */
1549 : /************************************************************************/
1550 :
1551 : /************************************************************************/
1552 : /* VRTFuncSource() */
1553 : /************************************************************************/
1554 :
1555 0 : VRTFuncSource::VRTFuncSource()
1556 :
1557 : {
1558 0 : pfnReadFunc = NULL;
1559 0 : pCBData = NULL;
1560 0 : fNoDataValue = (float) VRT_NODATA_UNSET;
1561 0 : eType = GDT_Byte;
1562 0 : }
1563 :
1564 : /************************************************************************/
1565 : /* ~VRTFuncSource() */
1566 : /************************************************************************/
1567 :
1568 0 : VRTFuncSource::~VRTFuncSource()
1569 :
1570 : {
1571 0 : }
1572 :
1573 : /************************************************************************/
1574 : /* SerializeToXML() */
1575 : /************************************************************************/
1576 :
1577 0 : CPLXMLNode *VRTFuncSource::SerializeToXML( const char * pszVRTPath )
1578 :
1579 : {
1580 0 : return NULL;
1581 : }
1582 :
1583 : /************************************************************************/
1584 : /* RasterIO() */
1585 : /************************************************************************/
1586 :
1587 : CPLErr
1588 0 : VRTFuncSource::RasterIO( int nXOff, int nYOff, int nXSize, int nYSize,
1589 : void *pData, int nBufXSize, int nBufYSize,
1590 : GDALDataType eBufType,
1591 : int nPixelSpace, int nLineSpace )
1592 :
1593 : {
1594 0 : if( nPixelSpace*8 == GDALGetDataTypeSize( eBufType )
1595 : && nLineSpace == nPixelSpace * nXSize
1596 : && nBufXSize == nXSize && nBufYSize == nYSize
1597 : && eBufType == eType )
1598 : {
1599 : return pfnReadFunc( pCBData,
1600 : nXOff, nYOff, nXSize, nYSize,
1601 0 : pData );
1602 : }
1603 : else
1604 : {
1605 : printf( "%d,%d %d,%d, %d,%d %d,%d %d,%d\n",
1606 : nPixelSpace*8, GDALGetDataTypeSize(eBufType),
1607 : nLineSpace, nPixelSpace * nXSize,
1608 : nBufXSize, nXSize,
1609 : nBufYSize, nYSize,
1610 0 : (int) eBufType, (int) eType );
1611 : CPLError( CE_Failure, CPLE_AppDefined,
1612 0 : "VRTFuncSource::RasterIO() - Irregular request." );
1613 0 : return CE_Failure;
1614 : }
1615 : }
1616 :
1617 : /************************************************************************/
1618 : /* GetMinimum() */
1619 : /************************************************************************/
1620 :
1621 0 : double VRTFuncSource::GetMinimum( int nXSize, int nYSize, int *pbSuccess )
1622 : {
1623 0 : *pbSuccess = FALSE;
1624 0 : return 0;
1625 : }
1626 :
1627 : /************************************************************************/
1628 : /* GetMaximum() */
1629 : /************************************************************************/
1630 :
1631 0 : double VRTFuncSource::GetMaximum( int nXSize, int nYSize, int *pbSuccess )
1632 : {
1633 0 : *pbSuccess = FALSE;
1634 0 : return 0;
1635 : }
1636 :
1637 : /************************************************************************/
1638 : /* VRTParseCoreSources() */
1639 : /************************************************************************/
1640 :
1641 332 : VRTSource *VRTParseCoreSources( CPLXMLNode *psChild, const char *pszVRTPath )
1642 :
1643 : {
1644 : VRTSource * poSource;
1645 :
1646 332 : if( EQUAL(psChild->pszValue,"AveragedSource")
1647 : || (EQUAL(psChild->pszValue,"SimpleSource")
1648 : && EQUALN(CPLGetXMLValue(psChild, "Resampling", "Nearest"),
1649 : "Aver",4)) )
1650 : {
1651 2 : poSource = new VRTAveragedSource();
1652 : }
1653 330 : else if( EQUAL(psChild->pszValue,"SimpleSource") )
1654 : {
1655 298 : poSource = new VRTSimpleSource();
1656 : }
1657 32 : else if( EQUAL(psChild->pszValue,"ComplexSource") )
1658 : {
1659 32 : poSource = new VRTComplexSource();
1660 : }
1661 : else
1662 : {
1663 : CPLError( CE_Failure, CPLE_AppDefined,
1664 0 : "VRTParseCoreSources() - Unknown source : %s", psChild->pszValue );
1665 0 : return NULL;
1666 : }
1667 :
1668 332 : if ( poSource->XMLInit( psChild, pszVRTPath ) == CE_None )
1669 331 : return poSource;
1670 :
1671 1 : delete poSource;
1672 1 : return NULL;
1673 : }
1674 :
|
