← Back to branch summary

~ubuntu-branches/debian/sid/gdal/sid

« back to all changes in this revision

Viewing changes to alg/gdalwarper.cpp

  • Committer: Package Import Robot
  • Author(s): Francesco Paolo Lovergine
  • Date: 2012-05-07 15:04:42 UTC
  • mfrom: (5.5.16 experimental)
  • Revision ID: package-import@ubuntu.com-20120507150442-2eks97loeh6rq005
Tags: 1.9.0-1
* Ready for sid, starting transition.
* All symfiles updated to latest builds.
* Added dh_numpy call in debian/rules to depend on numpy ABI.
* Policy bumped to 3.9.3, no changes required.

Show diffs side-by-side

added added

removed removed

Lines of Context:
alg/gdalwarper.cpp
1
1
/******************************************************************************
2
 
 * $Id: gdalwarper.cpp 20188 2010-08-06 15:04:44Z rouault $
 
2
 * $Id: gdalwarper.cpp 22888 2011-08-07 13:06:36Z rouault $
3
3
 *
4
4
 * Project:  High Performance Image Reprojector
5
5
 * Purpose:  Implementation of high level convenience APIs for warper.
31
31
#include "cpl_string.h"
32
32
#include "cpl_minixml.h"
33
33
#include "ogr_api.h"
 
34
#include "gdal_priv.h"
34
35
 
35
 
CPL_CVSID("$Id: gdalwarper.cpp 20188 2010-08-06 15:04:44Z rouault $");
 
36
CPL_CVSID("$Id: gdalwarper.cpp 22888 2011-08-07 13:06:36Z rouault $");
36
37
 
37
38
/************************************************************************/
38
39
/*                         GDALReprojectImage()                         */
47
48
 * implement the reprojection, and will default a variety of other 
48
49
 * warp options. 
49
50
 *
50
 
 * By default all bands are transferred, with no masking or nodata values
51
 
 * in effect.  No metadata, projection info, or color tables are transferred 
 
51
 * No metadata, projection info, or color tables are transferred
52
52
 * to the output file. 
53
53
 *
54
54
 * @param hSrcDS the source image file. 
158
158
        GDALRasterBandH hBand = GDALGetRasterBand( hSrcDS, iBand+1 );
159
159
        int             bGotNoData = FALSE;
160
160
        double          dfNoDataValue;
161
 
        
 
161
 
 
162
        if (GDALGetRasterColorInterpretation(hBand) == GCI_AlphaBand)
 
163
        {
 
164
            psWOptions->nSrcAlphaBand = iBand + 1;
 
165
        }
 
166
 
162
167
        dfNoDataValue = GDALGetRasterNoDataValue( hBand, &bGotNoData );
163
168
        if( bGotNoData )
164
169
        {
180
185
 
181
186
            psWOptions->padfSrcNoDataReal[iBand] = dfNoDataValue;
182
187
        }
 
188
 
 
189
        hBand = GDALGetRasterBand( hDstDS, iBand+1 );
 
190
        if (hBand && GDALGetRasterColorInterpretation(hBand) == GCI_AlphaBand)
 
191
        {
 
192
            psWOptions->nDstAlphaBand = iBand + 1;
 
193
        }
183
194
    }
184
195
 
185
196
/* -------------------------------------------------------------------- */
409
420
          float fNoData = (float) padfNoData[0];
410
421
          float *pafData = (float *) *ppImageData;
411
422
          int iOffset;
412
 
 
413
 
          // nothing to do if value is out of range.
414
 
          if( padfNoData[1] != 0.0 )
415
 
              return CE_None;
416
 
 
417
 
          for( iOffset = nXSize*nYSize-1; iOffset >= 0; iOffset-- )
418
 
          {
419
 
              if( pafData[iOffset] == fNoData )
 
423
          int bIsNoDataNan = CPLIsNan(fNoData);
 
424
 
 
425
          // nothing to do if value is out of range.
 
426
          if( padfNoData[1] != 0.0 )
 
427
              return CE_None;
 
428
 
 
429
          for( iOffset = nXSize*nYSize-1; iOffset >= 0; iOffset-- )
 
430
          {
 
431
              float fVal = pafData[iOffset];
 
432
              if( (bIsNoDataNan && CPLIsNan(fVal)) || (!bIsNoDataNan && ARE_REAL_EQUAL(fVal, fNoData)) )
 
433
              {
 
434
                  panValidityMask[iOffset>>5] &= ~(0x01 << (iOffset & 0x1f));
 
435
              }
 
436
          }
 
437
      }
 
438
      break;
 
439
      
 
440
      case GDT_Float64:
 
441
      {
 
442
          double dfNoData = padfNoData[0];
 
443
          double *padfData = (double *) *ppImageData;
 
444
          int iOffset;
 
445
          int bIsNoDataNan = CPLIsNan(dfNoData);
 
446
 
 
447
          // nothing to do if value is out of range.
 
448
          if( padfNoData[1] != 0.0 )
 
449
              return CE_None;
 
450
 
 
451
          for( iOffset = nXSize*nYSize-1; iOffset >= 0; iOffset-- )
 
452
          {
 
453
              double dfVal = padfData[iOffset];
 
454
              if( (bIsNoDataNan && CPLIsNan(dfVal)) || (!bIsNoDataNan && ARE_REAL_EQUAL(dfVal, dfNoData)) )
420
455
              {
421
456
                  panValidityMask[iOffset>>5] &= ~(0x01 << (iOffset & 0x1f));
422
457
              }
429
464
          double  *padfWrk;
430
465
          int     iLine, iPixel;
431
466
          int     nWordSize = GDALGetDataTypeSize(eType)/8;
 
467
          
 
468
          int bIsNoDataRealNan = CPLIsNan(padfNoData[0]);
 
469
          int bIsNoDataImagNan = CPLIsNan(padfNoData[1]);
432
470
 
433
471
          padfWrk = (double *) CPLMalloc(nXSize * sizeof(double) * 2);
434
472
          for( iLine = 0; iLine < nYSize; iLine++ )
439
477
              
440
478
              for( iPixel = 0; iPixel < nXSize; iPixel++ )
441
479
              {
442
 
                  if( padfWrk[iPixel*2] == padfNoData[0]
443
 
                      && padfWrk[iPixel*2+1] == padfNoData[1] )
 
480
                  if( ((bIsNoDataRealNan && CPLIsNan(padfWrk[iPixel*2])) ||
 
481
                       (!bIsNoDataRealNan && ARE_REAL_EQUAL(padfWrk[iPixel*2], padfNoData[0])))
 
482
                      && ((bIsNoDataImagNan && CPLIsNan(padfWrk[iPixel*2+1])) ||
 
483
                          (!bIsNoDataImagNan && ARE_REAL_EQUAL(padfWrk[iPixel*2+1], padfNoData[1]))) )
444
484
                  {
445
485
                      int iOffset = iPixel + iLine * nXSize;
446
486
                      
520
560
}
521
561
 
522
562
/************************************************************************/
 
563
/*                       GDALWarpSrcMaskMasker()                        */
 
564
/*                                                                      */
 
565
/*      GDALMaskFunc for reading source simple 8bit validity mask       */
 
566
/*      information and building a one bit validity mask.               */
 
567
/************************************************************************/
 
568
 
 
569
CPLErr 
 
570
GDALWarpSrcMaskMasker( void *pMaskFuncArg, int nBandCount, GDALDataType eType, 
 
571
                       int nXOff, int nYOff, int nXSize, int nYSize,
 
572
                       GByte ** /*ppImageData */,
 
573
                       int bMaskIsFloat, void *pValidityMask )
 
574
 
 
575
{
 
576
    GDALWarpOptions *psWO = (GDALWarpOptions *) pMaskFuncArg;
 
577
    GUInt32  *panMask = (GUInt32 *) pValidityMask;
 
578
 
 
579
/* -------------------------------------------------------------------- */
 
580
/*      Do some minimal checking.                                       */
 
581
/* -------------------------------------------------------------------- */
 
582
    if( bMaskIsFloat )
 
583
    {
 
584
        CPLAssert( FALSE );
 
585
        return CE_Failure;
 
586
    }
 
587
 
 
588
    if( psWO == NULL )
 
589
    {
 
590
        CPLAssert( FALSE );
 
591
        return CE_Failure;
 
592
    }
 
593
 
 
594
/* -------------------------------------------------------------------- */
 
595
/*      Allocate a temporary buffer to read mask byte data into.        */
 
596
/* -------------------------------------------------------------------- */
 
597
    GByte *pabySrcMask;
 
598
 
 
599
    pabySrcMask = (GByte *) VSIMalloc2(nXSize,nYSize);
 
600
    if( pabySrcMask == NULL )
 
601
    {
 
602
        CPLError( CE_Failure, CPLE_OutOfMemory,
 
603
                  "Failed to allocate pabySrcMask (%dx%d) in GDALWarpSrcMaskMasker()", 
 
604
                  nXSize, nYSize );
 
605
        return CE_Failure;
 
606
    }
 
607
 
 
608
/* -------------------------------------------------------------------- */
 
609
/*      Fetch our mask band.                                            */
 
610
/* -------------------------------------------------------------------- */
 
611
    GDALRasterBandH hSrcBand, hMaskBand = NULL;
 
612
 
 
613
    hSrcBand = GDALGetRasterBand( psWO->hSrcDS, psWO->panSrcBands[0] );
 
614
    if( hSrcBand != NULL )
 
615
        hMaskBand = GDALGetMaskBand( hSrcBand );
 
616
 
 
617
    if( hMaskBand == NULL )
 
618
    {
 
619
        CPLAssert( FALSE );
 
620
        return CE_Failure;
 
621
    }
 
622
 
 
623
/* -------------------------------------------------------------------- */
 
624
/*      Read the mask band.                                             */
 
625
/* -------------------------------------------------------------------- */
 
626
    CPLErr eErr;
 
627
 
 
628
    eErr = GDALRasterIO( hMaskBand, GF_Read, nXOff, nYOff, nXSize, nYSize, 
 
629
                         pabySrcMask, nXSize, nYSize, GDT_Byte, 0, 0 );
 
630
 
 
631
    if( eErr != CE_None )
 
632
    {
 
633
        CPLFree( pabySrcMask );
 
634
        return eErr;
 
635
    }
 
636
 
 
637
/* -------------------------------------------------------------------- */
 
638
/*      Pack into 1 bit per pixel for validity.                         */
 
639
/* -------------------------------------------------------------------- */
 
640
    for( int iPixel = nXSize * nYSize - 1; iPixel >= 0; iPixel-- )
 
641
    {                                    
 
642
        if( pabySrcMask[iPixel] == 0 )
 
643
            panMask[iPixel>>5] &= ~(0x01 << (iPixel & 0x1f));
 
644
    }
 
645
 
 
646
    CPLFree( pabySrcMask );
 
647
 
 
648
    return CE_None;
 
649
}
 
650
 
 
651
/************************************************************************/
523
652
/*                       GDALWarpDstAlphaMasker()                       */
524
653
/*                                                                      */
525
654
/*      GDALMaskFunc for reading or writing the destination simple      */
599
728
    else
600
729
    {
601
730
        for( iPixel = nXSize * nYSize - 1; iPixel >= 0; iPixel-- )
602
 
            pafMask[iPixel] = (int) ( pafMask[iPixel] * 255.1 );
 
731
            pafMask[iPixel] = (float)(int) ( pafMask[iPixel] * 255.1 );
603
732
        
604
733
        // Write data.
605
734
 
688
817
 *
689
818
 * - CUTLINE: This may contain the WKT geometry for a cutline.  It will
690
819
 * be converted into a geometry by GDALWarpOperation::Initialize() and assigned
691
 
 * to the GDALWarpOptions hCutline field.
 
820
 * to the GDALWarpOptions hCutline field. The coordinates must be expressed
 
821
 * in source pixel/line coordinates. Note: this is different from the assumptions
 
822
 * made for the -cutline option of the gdalwarp utility !
692
823
 *
693
824
 * - CUTLINE_BLEND_DIST: This may be set with a distance in pixels which
694
825
 * will be assigned to the dfCutlineBlendDist field in the GDALWarpOptions.
753
884
 
754
885
 
755
886
#define COPY_MEM(target,type,count)                                     \
756
 
   if( (psSrcOptions->target) != NULL && (count) != 0 )                 \
 
887
   do { if( (psSrcOptions->target) != NULL && (count) != 0 )            \
757
888
   {                                                                    \
758
 
       (psDstOptions->target) = (type *) CPLMalloc(sizeof(type)*count); \
 
889
       (psDstOptions->target) = (type *) CPLMalloc(sizeof(type)*(count)); \
759
890
       memcpy( (psDstOptions->target), (psSrcOptions->target),          \
760
 
               sizeof(type) * count );                                  \
761
 
   }
 
891
               sizeof(type) * (count) );                                        \
 
892
   } \
 
893
   else \
 
894
       (psDstOptions->target) = NULL; } while(0)
762
895
 
763
896
/************************************************************************/
764
897
/*                        GDALCloneWarpOptions()                        */
784
917
    COPY_MEM( padfDstNoDataImag, double, psSrcOptions->nBandCount );
785
918
    COPY_MEM( papfnSrcPerBandValidityMaskFunc, GDALMaskFunc, 
786
919
              psSrcOptions->nBandCount );
 
920
    psDstOptions->papSrcPerBandValidityMaskFuncArg = NULL;
787
921
 
788
922
    if( psSrcOptions->hCutline != NULL )
789
923
        psDstOptions->hCutline = 
927
1061
                CXT_Text, CPLString().Printf( "%d", psWO->panDstBands[i] ) );
928
1062
        
929
1063
        if( psWO->padfSrcNoDataReal != NULL )
930
 
            CPLCreateXMLElementAndValue( 
931
 
                psBand, "SrcNoDataReal", 
932
 
                CPLString().Printf( "%.16g", psWO->padfSrcNoDataReal[i] ) );
 
1064
        {
 
1065
            if (CPLIsNan(psWO->padfSrcNoDataReal[i]))
 
1066
                CPLCreateXMLElementAndValue(psBand, "SrcNoDataReal", "nan");
 
1067
            else
 
1068
                CPLCreateXMLElementAndValue( 
 
1069
                    psBand, "SrcNoDataReal", 
 
1070
                    CPLString().Printf( "%.16g", psWO->padfSrcNoDataReal[i] ) );
 
1071
        }
933
1072
 
934
1073
        if( psWO->padfSrcNoDataImag != NULL )
935
 
            CPLCreateXMLElementAndValue( 
936
 
                psBand, "SrcNoDataImag", 
937
 
                CPLString().Printf( "%.16g", psWO->padfSrcNoDataImag[i] ) );
 
1074
        {
 
1075
            if (CPLIsNan(psWO->padfSrcNoDataImag[i]))
 
1076
                CPLCreateXMLElementAndValue(psBand, "SrcNoDataImag", "nan");
 
1077
            else
 
1078
                CPLCreateXMLElementAndValue( 
 
1079
                    psBand, "SrcNoDataImag", 
 
1080
                    CPLString().Printf( "%.16g", psWO->padfSrcNoDataImag[i] ) );
 
1081
        }
938
1082
 
939
1083
        if( psWO->padfDstNoDataReal != NULL )
940
 
            CPLCreateXMLElementAndValue( 
941
 
                psBand, "DstNoDataReal", 
942
 
                CPLString().Printf( "%.16g", psWO->padfDstNoDataReal[i] ) );
 
1084
        {
 
1085
            if (CPLIsNan(psWO->padfDstNoDataReal[i]))
 
1086
                CPLCreateXMLElementAndValue(psBand, "DstNoDataReal", "nan");
 
1087
            else
 
1088
                CPLCreateXMLElementAndValue( 
 
1089
                    psBand, "DstNoDataReal", 
 
1090
                    CPLString().Printf( "%.16g", psWO->padfDstNoDataReal[i] ) );
 
1091
        }
943
1092
 
944
1093
        if( psWO->padfDstNoDataImag != NULL )
945
 
            CPLCreateXMLElementAndValue( 
946
 
                psBand, "DstNoDataImag", 
947
 
                CPLString().Printf( "%.16g", psWO->padfDstNoDataImag[i] ) );
 
1094
        {
 
1095
            if (CPLIsNan(psWO->padfDstNoDataImag[i]))
 
1096
                CPLCreateXMLElementAndValue(psBand, "DstNoDataImag", "nan");
 
1097
            else
 
1098
                CPLCreateXMLElementAndValue( 
 
1099
                    psBand, "DstNoDataImag", 
 
1100
                    CPLString().Printf( "%.16g", psWO->padfDstNoDataImag[i] ) );
 
1101
        }
948
1102
    }
949
1103
 
950
1104
/* -------------------------------------------------------------------- */
1091
1245
 
1092
1246
    psWO->nBandCount = 0;
1093
1247
    
1094
 
    for( psBand=psBandTree->psChild; psBand != NULL; psBand = psBand->psNext )
 
1248
    if (psBandTree)
 
1249
        psBand = psBandTree->psChild;
 
1250
    else
 
1251
        psBand = NULL;
 
1252
 
 
1253
    for( ; psBand != NULL; psBand = psBand->psNext )
1095
1254
    {
1096
1255
        if( psBand->eType != CXT_Element 
1097
1256
            || !EQUAL(psBand->pszValue,"BandMapping") )
1105
1264
/* ==================================================================== */
1106
1265
    int iBand = 0;
1107
1266
 
1108
 
    for( psBand=psBandTree->psChild; psBand != NULL; psBand = psBand->psNext )
 
1267
    if (psBandTree)
 
1268
        psBand = psBandTree->psChild;
 
1269
    else
 
1270
        psBand = NULL;
 
1271
 
 
1272
    for( ; psBand != NULL; psBand = psBand->psNext )
1109
1273
    {
1110
1274
        if( psBand->eType != CXT_Element 
1111
1275
            || !EQUAL(psBand->pszValue,"BandMapping") )
1146
1310
                psWO->padfSrcNoDataReal = 
1147
1311
                    (double *) CPLCalloc(sizeof(double),psWO->nBandCount);
1148
1312
 
1149
 
            psWO->padfSrcNoDataReal[iBand] = atof(pszValue);
 
1313
            psWO->padfSrcNoDataReal[iBand] = CPLAtofM(pszValue);
1150
1314
        }
1151
1315
        
1152
1316
        pszValue = CPLGetXMLValue(psBand,"SrcNoDataImag",NULL);
1156
1320
                psWO->padfSrcNoDataImag = 
1157
1321
                    (double *) CPLCalloc(sizeof(double),psWO->nBandCount);
1158
1322
 
1159
 
            psWO->padfSrcNoDataImag[iBand] = atof(pszValue);
 
1323
            psWO->padfSrcNoDataImag[iBand] = CPLAtofM(pszValue);
1160
1324
        }
1161
1325
        
1162
1326
/* -------------------------------------------------------------------- */
1169
1333
                psWO->padfDstNoDataReal = 
1170
1334
                    (double *) CPLCalloc(sizeof(double),psWO->nBandCount);
1171
1335
 
1172
 
            psWO->padfDstNoDataReal[iBand] = atof(pszValue);
 
1336
            psWO->padfDstNoDataReal[iBand] = CPLAtofM(pszValue);
1173
1337
        }
1174
1338
        
1175
1339
        pszValue = CPLGetXMLValue(psBand,"DstNoDataImag",NULL);
1179
1343
                psWO->padfDstNoDataImag = 
1180
1344
                    (double *) CPLCalloc(sizeof(double),psWO->nBandCount);
1181
1345
 
1182
 
            psWO->padfDstNoDataImag[iBand] = atof(pszValue);
 
1346
            psWO->padfDstNoDataImag[iBand] = CPLAtofM(pszValue);
1183
1347
        }
1184
1348
        
1185
1349
        iBand++;