1 /******************************************************************************
2 * $Id: vrtfilters.cpp 28053 2014-12-04 09:31:07Z rouault $
3 *
4 * Project: Virtual GDAL Datasets
5 * Purpose: Implementation of some filter types.
6 * Author: Frank Warmerdam <warmerdam@pobox.com>
7 *
8 ******************************************************************************
9 * Copyright (c) 2003, Frank Warmerdam <warmerdam@pobox.com>
10 * Copyright (c) 2008-2013, Even Rouault <even dot rouault at mines-paris dot org>
11 *
12 * Permission is hereby granted, free of charge, to any person obtaining a
13 * copy of this software and associated documentation files (the "Software"),
14 * to deal in the Software without restriction, including without limitation
15 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 * and/or sell copies of the Software, and to permit persons to whom the
17 * Software is furnished to do so, subject to the following conditions:
18 *
19 * The above copyright notice and this permission notice shall be included
20 * in all copies or substantial portions of the Software.
21 *
22 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 * DEALINGS IN THE SOFTWARE.
29 ****************************************************************************/
30
31 #include "vrtdataset.h"
32 #include "cpl_minixml.h"
33 #include "cpl_string.h"
34
35 CPL_CVSID("$Id: vrtfilters.cpp 28053 2014-12-04 09:31:07Z rouault $");
36
37 /************************************************************************/
38 /* ==================================================================== */
39 /* VRTFilteredSource */
40 /* ==================================================================== */
41 /************************************************************************/
42
43 /************************************************************************/
44 /* VRTFilteredSource() */
45 /************************************************************************/
46
VRTFilteredSource()47 VRTFilteredSource::VRTFilteredSource()
48
49 {
50 nExtraEdgePixels = 0;
51 nSupportedTypesCount = 1;
52 aeSupportedTypes[0] = GDT_Float32;
53 }
54
55 /************************************************************************/
56 /* ~VRTFilteredSource() */
57 /************************************************************************/
58
~VRTFilteredSource()59 VRTFilteredSource::~VRTFilteredSource()
60
61 {
62 }
63
64 /************************************************************************/
65 /* SetExtraEdgePixels() */
66 /************************************************************************/
67
SetExtraEdgePixels(int nEdgePixels)68 void VRTFilteredSource::SetExtraEdgePixels( int nEdgePixels )
69
70 {
71 nExtraEdgePixels = nEdgePixels;
72 }
73
74 /************************************************************************/
75 /* SetFilteringDataTypesSupported() */
76 /************************************************************************/
77
SetFilteringDataTypesSupported(int nTypeCount,GDALDataType * paeTypes)78 void VRTFilteredSource::SetFilteringDataTypesSupported( int nTypeCount,
79 GDALDataType *paeTypes)
80
81 {
82 if( nTypeCount >
83 (int) (sizeof(aeSupportedTypes)/sizeof(GDALDataType)) )
84 {
85 CPLAssert( FALSE );
86 nTypeCount = (int)
87 (sizeof(aeSupportedTypes)/sizeof(GDALDataType));
88 }
89
90 nSupportedTypesCount = nTypeCount;
91 memcpy( aeSupportedTypes, paeTypes, sizeof(GDALDataType) * nTypeCount );
92 }
93
94 /************************************************************************/
95 /* IsTypeSupported() */
96 /************************************************************************/
97
IsTypeSupported(GDALDataType eTestType)98 int VRTFilteredSource::IsTypeSupported( GDALDataType eTestType )
99
100 {
101 int i;
102
103 for( i = 0; i < nSupportedTypesCount; i++ )
104 {
105 if( eTestType == aeSupportedTypes[i] )
106 return TRUE;
107 }
108
109 return FALSE;
110 }
111
112 /************************************************************************/
113 /* RasterIO() */
114 /************************************************************************/
115
116 CPLErr
RasterIO(int nXOff,int nYOff,int nXSize,int nYSize,void * pData,int nBufXSize,int nBufYSize,GDALDataType eBufType,GSpacing nPixelSpace,GSpacing nLineSpace,GDALRasterIOExtraArg * psExtraArg)117 VRTFilteredSource::RasterIO( int nXOff, int nYOff, int nXSize, int nYSize,
118 void *pData, int nBufXSize, int nBufYSize,
119 GDALDataType eBufType,
120 GSpacing nPixelSpace,
121 GSpacing nLineSpace,
122 GDALRasterIOExtraArg* psExtraArg )
123
124 {
125 /* -------------------------------------------------------------------- */
126 /* For now we don't support filtered access to non-full */
127 /* resolution requests. Just collect the data directly without */
128 /* any operator. */
129 /* -------------------------------------------------------------------- */
130 if( nBufXSize != nXSize || nBufYSize != nYSize )
131 {
132 return VRTComplexSource::RasterIO( nXOff, nYOff, nXSize, nYSize,
133 pData, nBufXSize, nBufYSize,
134 eBufType, nPixelSpace, nLineSpace, psExtraArg );
135 }
136
137 // The window we will actually request from the source raster band.
138 double dfReqXOff, dfReqYOff, dfReqXSize, dfReqYSize;
139 int nReqXOff, nReqYOff, nReqXSize, nReqYSize;
140
141 // The window we will actual set _within_ the pData buffer.
142 int nOutXOff, nOutYOff, nOutXSize, nOutYSize;
143
144 if( !GetSrcDstWindow( nXOff, nYOff, nXSize, nYSize, nBufXSize, nBufYSize,
145 &dfReqXOff, &dfReqYOff, &dfReqXSize, &dfReqYSize,
146 &nReqXOff, &nReqYOff, &nReqXSize, &nReqYSize,
147 &nOutXOff, &nOutYOff, &nOutXSize, &nOutYSize ) )
148 return CE_None;
149
150 pData = ((GByte *)pData)
151 + nPixelSpace * nOutXOff
152 + nLineSpace * nOutYOff;
153
154 /* -------------------------------------------------------------------- */
155 /* Determine the data type we want to request. We try to match */
156 /* the source or destination request, and if both those fail we */
157 /* fallback to the first supported type at least as expressive */
158 /* as the request. */
159 /* -------------------------------------------------------------------- */
160 GDALDataType eOperDataType = GDT_Unknown;
161 int i;
162
163 if( IsTypeSupported( eBufType ) )
164 eOperDataType = eBufType;
165
166 if( eOperDataType == GDT_Unknown
167 && IsTypeSupported( poRasterBand->GetRasterDataType() ) )
168 eOperDataType = poRasterBand->GetRasterDataType();
169
170 if( eOperDataType == GDT_Unknown )
171 {
172 for( i = 0; i < nSupportedTypesCount; i++ )
173 {
174 if( GDALDataTypeUnion( aeSupportedTypes[i], eBufType )
175 == aeSupportedTypes[i] )
176 {
177 eOperDataType = aeSupportedTypes[i];
178 }
179 }
180 }
181
182 if( eOperDataType == GDT_Unknown )
183 {
184 eOperDataType = aeSupportedTypes[0];
185
186 for( i = 1; i < nSupportedTypesCount; i++ )
187 {
188 if( GDALGetDataTypeSize( aeSupportedTypes[i] )
189 > GDALGetDataTypeSize( eOperDataType ) )
190 {
191 eOperDataType = aeSupportedTypes[i];
192 }
193 }
194 }
195
196 /* -------------------------------------------------------------------- */
197 /* Allocate the buffer of data into which our imagery will be */
198 /* read, with the extra edge pixels as well. This will be the */
199 /* source data fed into the filter. */
200 /* -------------------------------------------------------------------- */
201 int nPixelOffset, nLineOffset;
202 int nExtraXSize = nOutXSize + 2 * nExtraEdgePixels;
203 int nExtraYSize = nOutYSize + 2 * nExtraEdgePixels;
204 GByte *pabyWorkData;
205
206 // FIXME? : risk of multiplication overflow
207 pabyWorkData = (GByte *)
208 VSICalloc( nExtraXSize * nExtraYSize,
209 (GDALGetDataTypeSize(eOperDataType) / 8) );
210
211 if( pabyWorkData == NULL )
212 {
213 CPLError( CE_Failure, CPLE_OutOfMemory,
214 "Work buffer allocation failed." );
215 return CE_Failure;
216 }
217
218 nPixelOffset = GDALGetDataTypeSize( eOperDataType ) / 8;
219 nLineOffset = nPixelOffset * nExtraXSize;
220
221 /* -------------------------------------------------------------------- */
222 /* Allocate the output buffer if the passed in output buffer is */
223 /* not of the same type as our working format, or if the passed */
224 /* in buffer has an unusual organization. */
225 /* -------------------------------------------------------------------- */
226 GByte *pabyOutData;
227
228 if( nPixelSpace != nPixelOffset || nLineSpace != nLineOffset
229 || eOperDataType != eBufType )
230 {
231 pabyOutData = (GByte *)
232 VSIMalloc3(nOutXSize, nOutYSize, nPixelOffset );
233
234 if( pabyOutData == NULL )
235 {
236 CPLError( CE_Failure, CPLE_OutOfMemory,
237 "Work buffer allocation failed." );
238 return CE_Failure;
239 }
240 }
241 else
242 pabyOutData = (GByte *) pData;
243
244 /* -------------------------------------------------------------------- */
245 /* Figure out the extended window that we want to load. Note */
246 /* that we keep track of the file window as well as the amount */
247 /* we will need to edge fill past the edge of the source dataset. */
248 /* -------------------------------------------------------------------- */
249 int nTopFill=0, nLeftFill=0, nRightFill=0, nBottomFill=0;
250 int nFileXOff, nFileYOff, nFileXSize, nFileYSize;
251
252 nFileXOff = nReqXOff - nExtraEdgePixels;
253 nFileYOff = nReqYOff - nExtraEdgePixels;
254 nFileXSize = nExtraXSize;
255 nFileYSize = nExtraYSize;
256
257 if( nFileXOff < 0 )
258 {
259 nLeftFill = -nFileXOff;
260 nFileXOff = 0;
261 nFileXSize -= nLeftFill;
262 }
263
264 if( nFileYOff < 0 )
265 {
266 nTopFill = -nFileYOff;
267 nFileYOff = 0;
268 nFileYSize -= nTopFill;
269 }
270
271 if( nFileXOff + nFileXSize > poRasterBand->GetXSize() )
272 {
273 nRightFill = nFileXOff + nFileXSize - poRasterBand->GetXSize();
274 nFileXSize -= nRightFill;
275 }
276
277 if( nFileYOff + nFileYSize > poRasterBand->GetYSize() )
278 {
279 nBottomFill = nFileYOff + nFileYSize - poRasterBand->GetYSize();
280 nFileYSize -= nBottomFill;
281 }
282
283 /* -------------------------------------------------------------------- */
284 /* Load the data. */
285 /* -------------------------------------------------------------------- */
286 CPLErr eErr;
287
288 eErr =
289 VRTComplexSource::RasterIOInternal( nFileXOff, nFileYOff, nFileXSize, nFileYSize,
290 pabyWorkData
291 + nLineOffset * nTopFill
292 + nPixelOffset * nLeftFill,
293 nFileXSize, nFileYSize, eOperDataType,
294 nPixelOffset, nLineOffset, psExtraArg );
295
296 if( eErr != CE_None )
297 {
298 if( pabyWorkData != pData )
299 VSIFree( pabyWorkData );
300
301 return eErr;
302 }
303
304 /* -------------------------------------------------------------------- */
305 /* Fill in missing areas. Note that we replicate the edge */
306 /* valid values out. We don't using "mirroring" which might be */
307 /* more suitable for some times of filters. We also don't mark */
308 /* these pixels as "nodata" though perhaps we should. */
309 /* -------------------------------------------------------------------- */
310 if( nLeftFill != 0 || nRightFill != 0 )
311 {
312 for( i = nTopFill; i < nExtraYSize - nBottomFill; i++ )
313 {
314 if( nLeftFill != 0 )
315 GDALCopyWords( pabyWorkData + nPixelOffset * nLeftFill
316 + i * nLineOffset, eOperDataType, 0,
317 pabyWorkData + i * nLineOffset, eOperDataType,
318 nPixelOffset, nLeftFill );
319
320 if( nRightFill != 0 )
321 GDALCopyWords( pabyWorkData + i * nLineOffset
322 + nPixelOffset * (nExtraXSize - nRightFill - 1),
323 eOperDataType, 0,
324 pabyWorkData + i * nLineOffset
325 + nPixelOffset * (nExtraXSize - nRightFill),
326 eOperDataType, nPixelOffset, nRightFill );
327 }
328 }
329
330 for( i = 0; i < nTopFill; i++ )
331 {
332 memcpy( pabyWorkData + i * nLineOffset,
333 pabyWorkData + nTopFill * nLineOffset,
334 nLineOffset );
335 }
336
337 for( i = nExtraYSize - nBottomFill; i < nExtraYSize; i++ )
338 {
339 memcpy( pabyWorkData + i * nLineOffset,
340 pabyWorkData + (nExtraYSize - nBottomFill - 1) * nLineOffset,
341 nLineOffset );
342 }
343
344 /* -------------------------------------------------------------------- */
345 /* Filter the data. */
346 /* -------------------------------------------------------------------- */
347 eErr = FilterData( nOutXSize, nOutYSize, eOperDataType,
348 pabyWorkData, pabyOutData );
349
350 VSIFree( pabyWorkData );
351 if( eErr != CE_None )
352 {
353 if( pabyOutData != pData )
354 VSIFree( pabyOutData );
355
356 return eErr;
357 }
358
359 /* -------------------------------------------------------------------- */
360 /* Copy from work buffer to target buffer. */
361 /* -------------------------------------------------------------------- */
362 if( pabyOutData != pData )
363 {
364 for( i = 0; i < nOutYSize; i++ )
365 {
366 GDALCopyWords( pabyOutData + i * (nPixelOffset * nOutXSize),
367 eOperDataType, nPixelOffset,
368 ((GByte *) pData) + i * nLineSpace,
369 eBufType, nPixelSpace, nOutXSize );
370 }
371
372 VSIFree( pabyOutData );
373 }
374
375 return CE_None;
376 }
377
378 /************************************************************************/
379 /* ==================================================================== */
380 /* VRTKernelFilteredSource */
381 /* ==================================================================== */
382 /************************************************************************/
383
384 /************************************************************************/
385 /* VRTKernelFilteredSource() */
386 /************************************************************************/
387
VRTKernelFilteredSource()388 VRTKernelFilteredSource::VRTKernelFilteredSource()
389
390 {
391 GDALDataType aeSupTypes[] = { GDT_Float32 };
392 padfKernelCoefs = NULL;
393 nKernelSize = 0;
394 bNormalized = FALSE;
395
396 SetFilteringDataTypesSupported( 1, aeSupTypes );
397 }
398
399 /************************************************************************/
400 /* ~VRTKernelFilteredSource() */
401 /************************************************************************/
402
~VRTKernelFilteredSource()403 VRTKernelFilteredSource::~VRTKernelFilteredSource()
404
405 {
406 CPLFree( padfKernelCoefs );
407 }
408
409 /************************************************************************/
410 /* SetNormalized() */
411 /************************************************************************/
412
SetNormalized(int bNormalizedIn)413 void VRTKernelFilteredSource::SetNormalized( int bNormalizedIn )
414
415 {
416 bNormalized = bNormalizedIn;
417 }
418
419 /************************************************************************/
420 /* SetKernel() */
421 /************************************************************************/
422
SetKernel(int nNewKernelSize,double * padfNewCoefs)423 CPLErr VRTKernelFilteredSource::SetKernel( int nNewKernelSize,
424 double *padfNewCoefs )
425
426 {
427 if( nNewKernelSize < 1 || (nNewKernelSize % 2) != 1 )
428 {
429 CPLError( CE_Failure, CPLE_AppDefined,
430 "Illegal filtering kernel size %d, must be odd positive number.",
431 nNewKernelSize );
432 return CE_Failure;
433 }
434
435 CPLFree( padfKernelCoefs );
436 nKernelSize = nNewKernelSize;
437
438 padfKernelCoefs = (double *)
439 CPLMalloc(sizeof(double) * nKernelSize * nKernelSize );
440 memcpy( padfKernelCoefs, padfNewCoefs,
441 sizeof(double) * nKernelSize * nKernelSize );
442
443 SetExtraEdgePixels( (nNewKernelSize - 1) / 2 );
444
445 return CE_None;
446 }
447
448 /************************************************************************/
449 /* FilterData() */
450 /************************************************************************/
451
452 CPLErr VRTKernelFilteredSource::
FilterData(int nXSize,int nYSize,GDALDataType eType,GByte * pabySrcData,GByte * pabyDstData)453 FilterData( int nXSize, int nYSize, GDALDataType eType,
454 GByte *pabySrcData, GByte *pabyDstData )
455
456 {
457 /* -------------------------------------------------------------------- */
458 /* Validate data type. */
459 /* -------------------------------------------------------------------- */
460 if( eType != GDT_Float32 )
461 {
462 CPLError( CE_Failure, CPLE_AppDefined,
463 "Unsupported data type (%s) in VRTKernelFilteredSource::FilterData()",
464 GDALGetDataTypeName( eType ) );
465 return CE_Failure;
466 }
467
468 CPLAssert( nExtraEdgePixels*2 + 1 == nKernelSize ||
469 (nKernelSize == 0 && nExtraEdgePixels == 0) );
470
471 /* -------------------------------------------------------------------- */
472 /* Float32 case. */
473 /* -------------------------------------------------------------------- */
474 if( eType == GDT_Float32 )
475 {
476 int iX, iY;
477
478 int bHasNoData;
479 float fNoData = (float) poRasterBand->GetNoDataValue(&bHasNoData);
480
481 for( iY = 0; iY < nYSize; iY++ )
482 {
483 for( iX = 0; iX < nXSize; iX++ )
484 {
485 int iYY, iKern = 0;
486 double dfSum = 0.0, dfKernSum = 0.0;
487 float fResult;
488 int iIndex = (iY+nKernelSize/2 ) * (nXSize+2*nExtraEdgePixels) + iX + nKernelSize/2;
489 float fCenter = ((float *)pabySrcData)[iIndex];
490
491 // Check if center srcpixel is NoData
492 if(!bHasNoData || fCenter != fNoData)
493 {
494 for( iYY = 0; iYY < nKernelSize; iYY++ )
495 {
496 int i;
497 float *pafData = ((float *)pabySrcData)
498 + (iY+iYY) * (nXSize+2*nExtraEdgePixels) + iX;
499
500 for( i = 0; i < nKernelSize; i++, pafData++, iKern++ )
501 {
502 if(!bHasNoData || *pafData != fNoData)
503 {
504 dfSum += *pafData * padfKernelCoefs[iKern];
505 dfKernSum += padfKernelCoefs[iKern];
506 }
507 }
508 }
509 if( bNormalized )
510 {
511 if( dfKernSum != 0.0 )
512 fResult = (float) (dfSum / dfKernSum);
513 else
514 fResult = 0.0;
515 }
516 else
517 fResult = (float) dfSum;
518
519 ((float *) pabyDstData)[iX + iY * nXSize] = fResult;
520 }
521 else
522 ((float *) pabyDstData)[iX + iY * nXSize] = fNoData;
523 }
524 }
525 }
526
527 return CE_None;
528 }
529
530 /************************************************************************/
531 /* XMLInit() */
532 /************************************************************************/
533
XMLInit(CPLXMLNode * psTree,const char * pszVRTPath)534 CPLErr VRTKernelFilteredSource::XMLInit( CPLXMLNode *psTree,
535 const char *pszVRTPath )
536
537 {
538 CPLErr eErr = VRTFilteredSource::XMLInit( psTree, pszVRTPath );
539 int nNewKernelSize, i, nCoefs;
540 double *padfNewCoefs;
541
542 if( eErr != CE_None )
543 return eErr;
544
545 nNewKernelSize = atoi(CPLGetXMLValue(psTree,"Kernel.Size","0"));
546
547 if( nNewKernelSize == 0 )
548 return CE_None;
549
550 char **papszCoefItems =
551 CSLTokenizeString( CPLGetXMLValue(psTree,"Kernel.Coefs","") );
552
553 nCoefs = CSLCount(papszCoefItems);
554
555 if( nCoefs != nNewKernelSize * nNewKernelSize )
556 {
557 CSLDestroy( papszCoefItems );
558 CPLError( CE_Failure, CPLE_AppDefined,
559 "Got wrong number of filter kernel coefficients (%s).\n"
560 "Expected %d, got %d.",
561 CPLGetXMLValue(psTree,"Kernel.Coefs",""),
562 nNewKernelSize * nNewKernelSize, nCoefs );
563 return CE_Failure;
564 }
565
566 padfNewCoefs = (double *) CPLMalloc(sizeof(double) * nCoefs);
567
568 for( i = 0; i < nCoefs; i++ )
569 padfNewCoefs[i] = CPLAtof(papszCoefItems[i]);
570
571 eErr = SetKernel( nNewKernelSize, padfNewCoefs );
572
573 CPLFree( padfNewCoefs );
574 CSLDestroy( papszCoefItems );
575
576 SetNormalized( atoi(CPLGetXMLValue(psTree,"Kernel.normalized","0")) );
577
578 return eErr;
579 }
580
581 /************************************************************************/
582 /* SerializeToXML() */
583 /************************************************************************/
584
SerializeToXML(const char * pszVRTPath)585 CPLXMLNode *VRTKernelFilteredSource::SerializeToXML( const char *pszVRTPath )
586
587 {
588 CPLXMLNode *psSrc = VRTFilteredSource::SerializeToXML( pszVRTPath );
589 CPLXMLNode *psKernel;
590 char *pszKernelCoefs;
591 int iCoef, nCoefCount = nKernelSize * nKernelSize;
592
593 if( psSrc == NULL )
594 return NULL;
595
596 CPLFree( psSrc->pszValue );
597 psSrc->pszValue = CPLStrdup("KernelFilteredSource" );
598
599 if( nKernelSize == 0 )
600 return psSrc;
601
602 psKernel = CPLCreateXMLNode( psSrc, CXT_Element, "Kernel" );
603
604 if( bNormalized )
605 CPLCreateXMLNode(
606 CPLCreateXMLNode( psKernel, CXT_Attribute, "normalized" ),
607 CXT_Text, "1" );
608 else
609 CPLCreateXMLNode(
610 CPLCreateXMLNode( psKernel, CXT_Attribute, "normalized" ),
611 CXT_Text, "0" );
612
613 pszKernelCoefs = (char *) CPLMalloc(nCoefCount * 32);
614
615 strcpy( pszKernelCoefs, "" );
616 for( iCoef = 0; iCoef < nCoefCount; iCoef++ )
617 CPLsprintf( pszKernelCoefs + strlen(pszKernelCoefs),
618 "%.8g ", padfKernelCoefs[iCoef] );
619
620 CPLSetXMLValue( psKernel, "Size", CPLSPrintf( "%d", nKernelSize ) );
621 CPLSetXMLValue( psKernel, "Coefs", pszKernelCoefs );
622
623 CPLFree( pszKernelCoefs );
624
625 return psSrc;
626 }
627
628 /************************************************************************/
629 /* VRTParseFilterSources() */
630 /************************************************************************/
631
VRTParseFilterSources(CPLXMLNode * psChild,const char * pszVRTPath)632 VRTSource *VRTParseFilterSources( CPLXMLNode *psChild, const char *pszVRTPath )
633
634 {
635 VRTSource *poSrc;
636
637 if( EQUAL(psChild->pszValue,"KernelFilteredSource") )
638 {
639 poSrc = new VRTKernelFilteredSource();
640 if( poSrc->XMLInit( psChild, pszVRTPath ) == CE_None )
641 return poSrc;
642 else
643 delete poSrc;
644 }
645
646 return NULL;
647 }
648
649