← Back to branch summary

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

« back to all changes in this revision

Viewing changes to alg/gdalwarpoperation.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/gdalwarpoperation.cpp
1
1
/******************************************************************************
2
 
 * $Id: gdalwarpoperation.cpp 19066 2010-03-13 20:26:18Z rouault $
 
2
 * $Id: gdalwarpoperation.cpp 22888 2011-08-07 13:06:36Z rouault $
3
3
 *
4
4
 * Project:  High Performance Image Reprojector
5
5
 * Purpose:  Implementation of the GDALWarpOperation class.
32
32
#include "cpl_multiproc.h"
33
33
#include "ogr_api.h"
34
34
 
35
 
CPL_CVSID("$Id: gdalwarpoperation.cpp 19066 2010-03-13 20:26:18Z rouault $");
 
35
CPL_CVSID("$Id: gdalwarpoperation.cpp 22888 2011-08-07 13:06:36Z rouault $");
36
36
 
37
37
/* Defined in gdalwarpkernel.cpp */
38
38
int GWKGetFilterRadius(GDALResampleAlg eResampleAlg);
121
121
{
122
122
    psOptions = NULL;
123
123
 
124
 
    dfProgressBase = 0.0;
125
 
    dfProgressScale = 1.0;
126
 
 
127
124
    hThread1Mutex = NULL;
128
125
    hThread2Mutex = NULL;
129
126
    hIOMutex = NULL;
341
338
        }
342
339
    }
343
340
 
 
341
    if( psOptions->nSrcAlphaBand > 0)
 
342
    {
 
343
        if ( psOptions->hSrcDS == NULL ||
 
344
             psOptions->nSrcAlphaBand > GDALGetRasterCount(psOptions->hSrcDS) )
 
345
        {
 
346
            CPLError( CE_Failure, CPLE_IllegalArg,
 
347
                      "nSrcAlphaBand = %d ... out of range for dataset.",
 
348
                      psOptions->nSrcAlphaBand );
 
349
            return FALSE;
 
350
        }
 
351
    }
 
352
 
 
353
    if( psOptions->nDstAlphaBand > 0)
 
354
    {
 
355
        if ( psOptions->hDstDS == NULL ||
 
356
             psOptions->nDstAlphaBand > GDALGetRasterCount(psOptions->hDstDS) )
 
357
        {
 
358
            CPLError( CE_Failure, CPLE_IllegalArg,
 
359
                      "nDstAlphaBand = %d ... out of range for dataset.",
 
360
                      psOptions->nDstAlphaBand );
 
361
            return FALSE;
 
362
        }
 
363
    }
 
364
 
344
365
    if( psOptions->nSrcAlphaBand > 0 
345
366
        && psOptions->pfnSrcDensityMaskFunc != NULL )
346
367
    {
674
695
        double dfChunkPixels = panThisChunk[2] * (double) panThisChunk[3];
675
696
        CPLErr eErr;
676
697
 
677
 
        dfProgressBase = dfPixelsProcessed / dfTotalPixels;
678
 
        dfProgressScale = dfChunkPixels / dfTotalPixels;
 
698
        double dfProgressBase = dfPixelsProcessed / dfTotalPixels;
 
699
        double dfProgressScale = dfChunkPixels / dfTotalPixels;
679
700
 
680
701
        eErr = WarpRegion( panThisChunk[0], panThisChunk[1], 
681
702
                           panThisChunk[2], panThisChunk[3],
682
703
                           panThisChunk[4], panThisChunk[5],
683
 
                           panThisChunk[6], panThisChunk[7] );
 
704
                           panThisChunk[6], panThisChunk[7],
 
705
                           dfProgressBase, dfProgressScale);
684
706
 
685
707
        if( eErr != CE_None )
686
708
            return eErr;
716
738
/*                          ChunkThreadMain()                           */
717
739
/************************************************************************/
718
740
 
 
741
typedef struct
 
742
{
 
743
    void              *hThreadMutex;
 
744
    GDALWarpOperation *poOperation;
 
745
    int               *panChunkInfo;
 
746
    int                bFinished;
 
747
    CPLErr             eErr;
 
748
    double             dfProgressBase;
 
749
    double             dfProgressScale;
 
750
} ChunkThreadData;
 
751
 
 
752
 
719
753
static void ChunkThreadMain( void *pThreadData )
720
754
 
721
755
{
722
 
    void *hThreadMutex = ((void **) pThreadData)[0];
723
 
    GDALWarpOperation *poOperation = 
724
 
        (GDALWarpOperation *) (((void **) pThreadData)[1]);
725
 
    int *panChunkInfo = (int *) (((void **) pThreadData)[2]);
 
756
    volatile ChunkThreadData* psData = (volatile ChunkThreadData*) pThreadData;
726
757
 
727
 
    if( !CPLAcquireMutex( hThreadMutex, 2.0 ) )
 
758
    if( !CPLAcquireMutex( psData->hThreadMutex, 2.0 ) )
728
759
    {
729
760
        CPLError( CE_Failure, CPLE_AppDefined, 
730
761
                  "Failed to acquire thread mutex in ChunkThreadMain()." );
731
762
        return;
732
763
    }
733
764
 
734
 
    CPLErr eErr;
735
 
 
736
 
    eErr = 
737
 
        poOperation->WarpRegion( panChunkInfo[0], panChunkInfo[1], 
 
765
    int *panChunkInfo = psData->panChunkInfo;
 
766
    psData->eErr = psData->poOperation->WarpRegion(
 
767
                                 panChunkInfo[0], panChunkInfo[1],
738
768
                                 panChunkInfo[2], panChunkInfo[3], 
739
769
                                 panChunkInfo[4], panChunkInfo[5], 
740
 
                                 panChunkInfo[6], panChunkInfo[7] );
741
 
 
742
 
    /* Return error. */
743
 
    ((void **) pThreadData)[2] = (void *) (long) eErr;
 
770
                                 panChunkInfo[6], panChunkInfo[7],
 
771
                                 psData->dfProgressBase,
 
772
                                 psData->dfProgressScale);
744
773
 
745
774
    /* Marks that we are done. */
746
 
    ((void **) pThreadData)[1] = NULL;
 
775
    psData->bFinished = TRUE;
747
776
 
748
777
    /* Release mutex so parent knows we are done. */
749
 
    CPLReleaseMutex( hThreadMutex );
 
778
    CPLReleaseMutex( psData->hThreadMutex );
750
779
}
751
780
 
752
781
/************************************************************************/
800
829
/*      Process them one at a time, updating the progress               */
801
830
/*      information for each region.                                    */
802
831
/* -------------------------------------------------------------------- */
803
 
    void * volatile papThreadDataList[6] = { hThread1Mutex, NULL, NULL, 
804
 
                                            hThread2Mutex, NULL, NULL };
 
832
    ChunkThreadData volatile asThreadData[2];
 
833
    memset((void*)&asThreadData, 0, sizeof(asThreadData));
 
834
    asThreadData[0].hThreadMutex = hThread1Mutex;
 
835
    asThreadData[0].poOperation = this;
 
836
    asThreadData[0].bFinished = TRUE;
 
837
    asThreadData[1].hThreadMutex = hThread2Mutex;
 
838
    asThreadData[1].poOperation = this;
 
839
    asThreadData[1].bFinished = TRUE;
 
840
 
805
841
    int iChunk;
806
842
    double dfPixelsProcessed=0.0, dfTotalPixels = nDstXSize*(double)nDstYSize;
807
843
    CPLErr eErr = CE_None;
818
854
            int *panThisChunk = panChunkList + iChunk*8;
819
855
            double dfChunkPixels = panThisChunk[2] * (double) panThisChunk[3];
820
856
 
821
 
            dfProgressBase = dfPixelsProcessed / dfTotalPixels;
822
 
            dfProgressScale = dfChunkPixels / dfTotalPixels;
 
857
            asThreadData[iThread].dfProgressBase = dfPixelsProcessed / dfTotalPixels;
 
858
            asThreadData[iThread].dfProgressScale = dfChunkPixels / dfTotalPixels;
823
859
 
824
860
            dfPixelsProcessed += dfChunkPixels;
825
 
            
826
 
            papThreadDataList[iThread*3+1] = (void *) this;
827
 
            papThreadDataList[iThread*3+2] = (void *) panThisChunk;
 
861
 
 
862
            asThreadData[iThread].panChunkInfo = panThisChunk;
 
863
            asThreadData[iThread].bFinished = FALSE;
828
864
 
829
865
            CPLDebug( "GDAL", "Start chunk %d.", iChunk );
830
 
            if( CPLCreateThread( ChunkThreadMain, 
831
 
                          (void *) (papThreadDataList + iThread*3)) == -1 )
 
866
            if( CPLCreateThread( ChunkThreadMain, (void*) &asThreadData[iThread]) == -1 )
832
867
            {
833
868
                CPLError( CE_Failure, CPLE_AppDefined, 
834
869
                          "CPLCreateThread() failed in ChunkAndWarpMulti()" );
850
885
            iThread = (iChunk-1) % 2;
851
886
            
852
887
            /* Wait for thread to finish. */
853
 
            while( papThreadDataList[iThread*3+1] != NULL )
 
888
            while( !(asThreadData[iThread].bFinished) )
854
889
            {
855
 
                if( CPLAcquireMutex( papThreadDataList[iThread*3+0], 1.0 ) )
856
 
                    CPLReleaseMutex( papThreadDataList[iThread*3+0] );
 
890
                if( CPLAcquireMutex( asThreadData[iThread].hThreadMutex, 1.0 ) )
 
891
                    CPLReleaseMutex( asThreadData[iThread].hThreadMutex );
857
892
            }
858
893
 
859
894
            CPLDebug( "GDAL", "Finished chunk %d.", iChunk-1 );
860
895
 
861
 
            eErr = (CPLErr) (long) (GIntBig) papThreadDataList[iThread*3+2];
 
896
            eErr = asThreadData[iThread].eErr;
862
897
 
863
898
            if( eErr != CE_None )
864
899
                break;
871
906
    int iThread;
872
907
    for(iThread = 0; iThread < 2; iThread ++)
873
908
    {
874
 
        while( papThreadDataList[iThread*3+1] != NULL )
 
909
        while( !(asThreadData[iThread].bFinished) )
875
910
        {
876
 
            if( CPLAcquireMutex( papThreadDataList[iThread*3+0], 1.0 ) )
877
 
                CPLReleaseMutex( papThreadDataList[iThread*3+0] );
 
911
            if( CPLAcquireMutex( asThreadData[iThread].hThreadMutex, 1.0 ) )
 
912
                CPLReleaseMutex( asThreadData[iThread].hThreadMutex );
878
913
        }
879
914
    }
880
915
 
953
988
        * psOptions->nBandCount;
954
989
 
955
990
    if( psOptions->pfnSrcDensityMaskFunc != NULL )
956
 
        nSrcPixelCostInBits += 32; /* float mask */
 
991
        nSrcPixelCostInBits += 32; /* ?? float mask */
 
992
 
 
993
    GDALRasterBandH hSrcBand = NULL;
 
994
    if( psOptions->nBandCount > 0 )
 
995
        hSrcBand = GDALGetRasterBand(psOptions->hSrcDS,
 
996
                                     psOptions->panSrcBands[0]);
 
997
 
 
998
    if( psOptions->nSrcAlphaBand > 0 || psOptions->hCutline != NULL )
 
999
        nSrcPixelCostInBits += 32; /* UnifiedSrcDensity float mask */
 
1000
    else if (hSrcBand != NULL && (GDALGetMaskFlags(hSrcBand) & GMF_PER_DATASET))
 
1001
        nSrcPixelCostInBits += 1; /* UnifiedSrcValid bit mask */
957
1002
 
958
1003
    if( psOptions->papfnSrcPerBandValidityMaskFunc != NULL 
959
1004
        || psOptions->padfSrcNoDataReal != NULL )
978
1023
        || psOptions->pfnDstValidityMaskFunc != NULL )
979
1024
        nDstPixelCostInBits += psOptions->nBandCount;
980
1025
 
 
1026
    if( psOptions->nDstAlphaBand > 0 )
 
1027
        nDstPixelCostInBits += 32; /* DstDensity float mask */
 
1028
 
981
1029
/* -------------------------------------------------------------------- */
982
1030
/*      Does the cost of the current rectangle exceed our memory        */
983
1031
/*      limit? If so, split the destination along the longest           */
1084
1132
 * \fn CPLErr GDALWarpOperation::WarpRegion(int nDstXOff, int nDstYOff, 
1085
1133
                                            int nDstXSize, int nDstYSize,
1086
1134
                                            int nSrcXOff=0, int nSrcYOff=0,
1087
 
                                            int nSrcXSize=0, int nSrcYSize=0 );
 
1135
                                            int nSrcXSize=0, int nSrcYSize=0,
 
1136
                                            double dfProgressBase = 0,
 
1137
                                            double dfProgressScale = 1);
1088
1138
 *
1089
1139
 * This method requests the indicated region of the output file be generated.
1090
1140
 * 
1096
1146
 * applications.  Use it instead if staying within memory constraints is
1097
1147
 * desired. 
1098
1148
 *
1099
 
 * Progress is reported from 0.0 to 1.0 for the indicated region. 
 
1149
 * Progress is reported from dfProgressBase to dfProgressBase + dfProgressScale
 
1150
 * for the indicated region.
1100
1151
 *
1101
1152
 * @param nDstXOff X offset to window of destination data to be produced.
1102
1153
 * @param nDstYOff Y offset to window of destination data to be produced.
1106
1157
 * @param nSrcYOff source window Y offset (computed if window all zero)
1107
1158
 * @param nSrcXSize source window X size (computed if window all zero)
1108
1159
 * @param nSrcYSize source window Y size (computed if window all zero)
 
1160
 * @param dfProgressBase minimum progress value reported
 
1161
 * @param dfProgressScale value such as dfProgressBase + dfProgressScale is the
 
1162
 *                        maximum progress value reported
1109
1163
 *
1110
1164
 * @return CE_None on success or CE_Failure if an error occurs.
1111
1165
 */
1113
1167
CPLErr GDALWarpOperation::WarpRegion( int nDstXOff, int nDstYOff, 
1114
1168
                                      int nDstXSize, int nDstYSize,
1115
1169
                                      int nSrcXOff, int nSrcYOff,
1116
 
                                      int nSrcXSize, int nSrcYSize )
 
1170
                                      int nSrcXSize, int nSrcYSize,
 
1171
                                      double dfProgressBase,
 
1172
                                      double dfProgressScale)
1117
1173
 
1118
1174
{
1119
1175
    CPLErr eErr;
1147
1203
        CPLError( CE_Failure, CPLE_AppDefined,
1148
1204
                  "Integer overflow : nDstXSize=%d, nDstYSize=%d",
1149
1205
                  nDstXSize, nDstYSize);
 
1206
        if( hIOMutex != NULL )
 
1207
            CPLReleaseMutex( hIOMutex );
1150
1208
        return CE_Failure;
1151
1209
    }
1152
1210
 
1156
1214
        CPLError( CE_Failure, CPLE_OutOfMemory,
1157
1215
                  "Out of memory allocating %d byte destination buffer.",
1158
1216
                  nBandSize * psOptions->nBandCount );
 
1217
        if( hIOMutex != NULL )
 
1218
            CPLReleaseMutex( hIOMutex );
1159
1219
        return CE_Failure;
1160
1220
    }
1161
1221
 
1201
1261
                memset( pBandData, 
1202
1262
                        MAX(0,MIN(255,(int)adfInitRealImag[0])), 
1203
1263
                        nBandSize);
1204
 
            else if( adfInitRealImag[0] == 0.0 && adfInitRealImag[1] == 0 )
 
1264
            else if( !CPLIsNan(adfInitRealImag[0]) && adfInitRealImag[0] == 0.0 &&
 
1265
                     !CPLIsNan(adfInitRealImag[1]) && adfInitRealImag[1] == 0.0 )
1205
1266
            {
1206
1267
                memset( pBandData, 0, nBandSize );
1207
1268
            }
1208
 
            else if( adfInitRealImag[1] == 0.0 )
 
1269
            else if( !CPLIsNan(adfInitRealImag[1]) && adfInitRealImag[1] == 0.0 )
1209
1270
            {
1210
1271
                GDALCopyWords( &adfInitRealImag, GDT_Float64, 0, 
1211
1272
                               pBandData,psOptions->eWorkingDataType,nWordSize,
1239
1300
        if( eErr != CE_None )
1240
1301
        {
1241
1302
            CPLFree( pDstBuffer );
 
1303
            if( hIOMutex != NULL )
 
1304
                CPLReleaseMutex( hIOMutex );
1242
1305
            return eErr;
1243
1306
        }
1244
1307
 
1250
1313
/* -------------------------------------------------------------------- */
1251
1314
    eErr = WarpRegionToBuffer( nDstXOff, nDstYOff, nDstXSize, nDstYSize, 
1252
1315
                               pDstBuffer, psOptions->eWorkingDataType, 
1253
 
                               nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize );
 
1316
                               nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
 
1317
                               dfProgressBase, dfProgressScale);
1254
1318
 
1255
1319
/* -------------------------------------------------------------------- */
1256
1320
/*      Write the output data back to disk if all went well.            */
1264
1328
                                    psOptions->nBandCount, 
1265
1329
                                    psOptions->panDstBands,
1266
1330
                                    0, 0, 0 );
1267
 
        if( CSLFetchBoolean( psOptions->papszWarpOptions, "WRITE_FLUSH", 
 
1331
        if( eErr == CE_None &&
 
1332
            CSLFetchBoolean( psOptions->papszWarpOptions, "WRITE_FLUSH",
1268
1333
                             FALSE ) )                                  
1269
1334
        {
 
1335
            CPLErr eOldErr = CPLGetLastErrorType();
 
1336
            CPLString osLastErrMsg = CPLGetLastErrorMsg();
1270
1337
            GDALFlushCache( psOptions->hDstDS );
 
1338
            CPLErr eNewErr = CPLGetLastErrorType();
 
1339
            if (eNewErr != eOldErr || osLastErrMsg.compare(CPLGetLastErrorMsg()) != 0)
 
1340
                eErr = CE_Failure;
1271
1341
        }
1272
1342
        ReportTiming( "Output buffer write" );
1273
1343
    }
1316
1386
                                  void *pDataBuf, 
1317
1387
                                  GDALDataType eBufDataType,
1318
1388
                                  int nSrcXOff=0, int nSrcYOff=0,
1319
 
                                  int nSrcXSize=0, int nSrcYSize=0 );
 
1389
                                  int nSrcXSize=0, int nSrcYSize=0,
 
1390
                                  double dfProgressBase = 0,
 
1391
                                  double dfProgressScale = 1 );
1320
1392
 * 
1321
1393
 * This method requests that a particular window of the output dataset
1322
1394
 * be warped and the result put into the provided data buffer.  The output
1338
1410
 * @param nSrcYOff source window Y offset (computed if window all zero)
1339
1411
 * @param nSrcXSize source window X size (computed if window all zero)
1340
1412
 * @param nSrcYSize source window Y size (computed if window all zero)
 
1413
 * @param dfProgressBase minimum progress value reported
 
1414
 * @param dfProgressScale value such as dfProgressBase + dfProgressScale is the
 
1415
 *                        maximum progress value reported
1341
1416
 *
1342
1417
 * @return CE_None on success or CE_Failure if an error occurs.
1343
1418
 */
1345
1420
CPLErr GDALWarpOperation::WarpRegionToBuffer( 
1346
1421
    int nDstXOff, int nDstYOff, int nDstXSize, int nDstYSize, 
1347
1422
    void *pDataBuf, GDALDataType eBufDataType,
1348
 
    int nSrcXOff, int nSrcYOff, int nSrcXSize, int nSrcYSize )
 
1423
    int nSrcXOff, int nSrcYOff, int nSrcXSize, int nSrcYSize,
 
1424
    double dfProgressBase, double dfProgressScale)
1349
1425
 
1350
1426
{
1351
1427
    CPLErr eErr = CE_None;
1465
1541
/* -------------------------------------------------------------------- */
1466
1542
/*      Generate a source density mask if we have a source alpha band   */
1467
1543
/* -------------------------------------------------------------------- */
1468
 
    if( eErr == CE_None && psOptions->nSrcAlphaBand > 0 )
 
1544
    if( eErr == CE_None && psOptions->nSrcAlphaBand > 0 &&
 
1545
        nSrcXSize > 0 && nSrcYSize > 0 )
1469
1546
    {
1470
1547
        CPLAssert( oWK.pafDstDensity == NULL );
1471
1548
 
1481
1558
                                        oWK.papabySrcImage,
1482
1559
                                        TRUE, oWK.pafUnifiedSrcDensity );
1483
1560
    }
1484
 
    
 
1561
 
1485
1562
/* -------------------------------------------------------------------- */
1486
1563
/*      Generate a source density mask if we have a source cutline.     */
1487
1564
/* -------------------------------------------------------------------- */
1488
 
    if( eErr == CE_None && psOptions->hCutline != NULL )
 
1565
    if( eErr == CE_None && psOptions->hCutline != NULL  &&
 
1566
        nSrcXSize > 0 && nSrcYSize > 0 )
1489
1567
    {
1490
1568
        if( oWK.pafUnifiedSrcDensity == NULL )
1491
1569
        {
1536
1614
/*      If we have source nodata values create, or update the           */
1537
1615
/*      validity mask.                                                  */
1538
1616
/* -------------------------------------------------------------------- */
1539
 
    if( eErr == CE_None && psOptions->padfSrcNoDataReal != NULL )
 
1617
    if( eErr == CE_None && psOptions->padfSrcNoDataReal != NULL &&
 
1618
        nSrcXSize > 0 && nSrcYSize > 0 )
1540
1619
    {
1541
1620
        for( i = 0; i < psOptions->nBandCount && eErr == CE_None; i++ )
1542
1621
        {
1590
1669
    }
1591
1670
 
1592
1671
/* -------------------------------------------------------------------- */
 
1672
/*      Generate a source validity mask if we have a source mask for     */
 
1673
/*      the whole input dataset (and didn't already treat it as         */
 
1674
/*      alpha band).                                                    */
 
1675
/* -------------------------------------------------------------------- */
 
1676
    GDALRasterBandH hSrcBand = NULL;
 
1677
    if( psOptions->nBandCount > 0 )
 
1678
        hSrcBand = GDALGetRasterBand(psOptions->hSrcDS,
 
1679
                                     psOptions->panSrcBands[0]);
 
1680
    
 
1681
    if( eErr == CE_None 
 
1682
        && oWK.pafUnifiedSrcDensity == NULL 
 
1683
        && (GDALGetMaskFlags(hSrcBand) & GMF_PER_DATASET) &&
 
1684
        nSrcXSize > 0 && nSrcYSize > 0 )
 
1685
 
 
1686
    {
 
1687
        eErr = CreateKernelMask( &oWK, 0, "UnifiedSrcValid" );
 
1688
        
 
1689
        if( eErr == CE_None )
 
1690
            eErr = 
 
1691
                GDALWarpSrcMaskMasker( psOptions, 
 
1692
                                       psOptions->nBandCount, 
 
1693
                                       psOptions->eWorkingDataType,
 
1694
                                       oWK.nSrcXOff, oWK.nSrcYOff, 
 
1695
                                       oWK.nSrcXSize, oWK.nSrcYSize,
 
1696
                                       oWK.papabySrcImage,
 
1697
                                       FALSE, oWK.panUnifiedSrcValid );
 
1698
    }
 
1699
    
 
1700
/* -------------------------------------------------------------------- */
1593
1701
/*      If we have destination nodata values create, or update the      */
1594
1702
/*      validity mask.  We clear the DstValid for any pixel that we     */
1595
1703
/*      do no have valid data in *any* of the source bands.             */