Opened 16 years ago

Closed 15 years ago

#2445 closed defect (fixed)

GDALAutoCreateWarpedVRT can cause a divide by zero in GDALWarpKernel::PerformWarp

Reported by: mpd Owned by: warmerdam
Priority: high Milestone: 1.5.4
Component: GDAL_Raster Version: unspecified
Severity: normal Keywords: warper
Cc:

Description

GDALWarpKernel::PerformWarp has these two lines:

    dfXScale = (double)nDstXSize / nSrcXSize;
    dfYScale = (double)nDstYSize / nSrcYSize;

Using an auto-warped VRT dataset can lead to nSrcYSize being 0, and hence a divide by zero.

To reproduce, add the following after the GDALOpenShared in gdal_translate.cpp (after also adding #include "gdalwarper.h" at the top of the file):

    {
        if ( GDALGetGCPCount(hDataset)>0 )
        {
            fprintf(stdout,"Creating VRTDataset\n");
            hDataset=GDALAutoCreateWarpedVRT(hDataset,NULL,NULL,GRA_NearestNeighbour,0.0,NULL);
            if( hDataset == NULL )
            {
                fprintf( stderr,
                         "GDALAutoCreateWarpedVRT failed - %d\n%s\n",
                         CPLGetLastErrorNo(), CPLGetLastErrorMsg() );
                GDALDestroyDriverManager();
                exit( 1 );
            }
        }
    }

Then try using the hacked gdal_translate on the attached file (an old GDAL sample).

Attachments (1)

ottawa_patch_small.tif.bz2 (955.8 KB ) - added by mpd 16 years ago.

Download all attachments as: .zip

Change History (14)

by mpd, 16 years ago

Attachment: ottawa_patch_small.tif.bz2 added

comment:1 by warmerdam, 16 years ago

Cc: warmerdam added
Keywords: warper added
Milestone: 1.5.3
Owner: changed from warmerdam to Mateusz Łoskot
Priority: normalhigh

comment:2 by Mateusz Łoskot, 16 years ago

Status: newassigned

comment:3 by Mateusz Łoskot, 16 years ago

Unfortunately, I'm not able to reproduce this problem.

I followed the instructions above to patch gdal_translate.cpp and run it against the attached file, but no error occurred:

$ gdal_translate ottawa_patch_small.tif out.tif
GDAL: GDALOpen(ottawa_patch_small.tif) succeeds as GTiff.
Creating VRTDataset
Input file size is 1126, 805
0GDAL: GDALOpen(out.tif) succeeds as GTiff.
GDAL: GDALDatasetCopyWholeRaster(): 805 line swath, bInterleave=0
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=0,0,476x103 Dst=0,0,512x128
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=451,0,444x168 Dst=512,0,512x128
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=895,5,0x228 Dst=1024,0,512x128
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=0,35,452x230 Dst=0,128,512x128
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=428,102,467x229 Dst=512,128,512x128
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=885,167,10x229 Dst=1024,128,512x128
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=0,197,429x230 Dst=0,256,512x128
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=404,264,482x229 Dst=512,256,512x128
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=862,330,33x228 Dst=1024,256,512x128
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=0,360,405x230 Dst=0,384,512x128
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=380,426,483x230 Dst=512,384,512x128
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=838,492,57x229 Dst=1024,384,512x128
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=0,522,381x230 Dst=0,512,512x128
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=356,589,483x229 Dst=512,512,512x128
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=815,655,80x228 Dst=1024,512,512x128
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=0,684,357x229 Dst=0,640,512x128
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=332,751,484x162 Dst=512,640,512x128
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=791,817,104x96 Dst=1024,640,512x128
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=0,846,333x67 Dst=0,768,512x128
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=308,913,484x0 Dst=512,768,512x128
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=767,913,128x0 Dst=1024,768,512x128
...10...20...30...40...50...60...70...80...90...100 - done.
GDAL: GDALClose(out.tif)
GDAL: GDALClose()
Open GDAL Datasets:
  1 S GTiff  -1227987248 895x913x1 ottawa_patch_small.tif
GDAL: GDALDeregister_GTiff() called.

comment:4 by warmerdam, 16 years ago

Mateusz,

You might want to try putting an "assert( nSrcXSize != 0 );" near the line mentioned to see if we are getting a divide by zero. The actual result of a div-by-zero may vary by platform.

comment:5 by Mateusz Łoskot, 16 years ago

Frank,

You are perfectly right the assert should catch the problem:

$ gdal_translate ottawa_patch_small.tif out.tif
GDAL: GDALOpen(ottawa_patch_small.tif) succeeds as GTiff.
Creating VRTDataset
Input file size is 1126, 805
0GDAL: GDALOpen(out.tif) succeeds as GTiff.
GDAL: GDALDatasetCopyWholeRaster(): 805 line swath, bInterleave=0
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=0,0,476x103 Dst=0,0,512x128
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=451,0,444x168 Dst=512,0,512x128
ERROR 7: Assertion `0 != nSrcXSize' failed
in file `gdalwarpkernel.cpp', line 531

Aborted (core dumped)

in reply to:  5 comment:6 by mpd, 16 years ago

Alternatively you can add something like this to turn the divide-by-zero to an exception (on Windows):

_controlfp(~(_EM_ZERODIVIDE|_EM_OVERFLOW|_EM_INVALID),_MCW_EM);

comment:7 by Mateusz Łoskot, 16 years ago

I was trying to track the problem with debugger and find source of Zero dimension. Actually, I've found it but I'm not enough experienced with warping algorithm to tell what is the right solution on that. Here is what I've found:

The processing fails after 3rd call of VRTWarpedDataset::ProcessBlock function with following parameters:

VRTWarpedDataset::ProcessBlock (this=0x8098d60, iBlockX=2, iBlockY=0) at vrtwarped.cpp:807

More specifically, the 3rd call of ProcessBlock makes a call to GDALWarpOperation::WarpRegionToBuffer function, passing nSrcXSize = 0:

GDALWarpOperation::WarpRegionToBuffer (this=0x8098ec0, nDstXOff=1024, nDstYOff=0, nDstXSize=512, nDstYSize=128, pDataBuf=0x8144710, eBufDataType=GDT_UInt16, nSrcXOff=895, nSrcYOff=5, 
    nSrcXSize=0, nSrcYSize=228) at gdalwarpoperation.cpp:1561 

This causes WarpRegionToBuffer to compute source window and this computation leads to incorrect value of X dimention, equal to Zero:

...
XXX: ----- SOURCE OF ZERO X SIZE (3rd call of ProcessBlock) -----
XXX: GDALGetRasterXSize(psOptions->hSrcDS) is 895 and *pnSrcXOff is 895
XXX: GDALGetRasterXSize(psOptions->hSrcDS) - *pnSrcXOff = 0
XXX: ceil(dfMaxXOut) - *pnSrcXOff + nResWinSize = 497
XXX: pnSrcXSize = 0
XXX: pnSrcYSize = 228
XXX: ----- END OF SOURCE OF ZERO X SIZE -----
GDAL: Computed source window: 895,5,0x228
XXX: GDALWarpKernel ctor
ERROR 7: Assertion `0 != nSrcXSize' failed
in file `gdalwarpkernel.cpp', line 549

As I said above, I'm unfortunately unable to suggest a solution but I sense that this message from the above is important:

GDALGetRasterXSize(psOptions->hSrcDS) is 895 and *pnSrcXOff is 895

and if I got it correctly, the WarpRegionToBuffer function should detect situation when X or Y dimension is equal to X or Y offset, and handle it somehow, probably skip such case(s).


For easier understanding of the explanation above, here is complete output of execution of gdal_translate (patched to use autowarped VR) with my debug messages marked using XXX:

Starting program: /home/mloskot/dev/gdal/_svn/trunk/gdal/apps/gdal_translate ottawa_patch_small.tif out.tif
[Thread debugging using libthread_db enabled]
[New Thread -1227917616 (LWP 18097)]
GDAL: GDALOpen(ottawa_patch_small.tif) succeeds as GTiff.
Creating VRTDataset
XXX: Pixels: 1126
XXX: Lines: 805
Input file size is 1126, 805
0GDAL: GDALOpen(out.tif) succeeds as GTiff.
GDAL: GDALDatasetCopyWholeRaster(): 805 line swath, bInterleave=0
XXX: ----- SOURCE OF ZERO X SIZE -----
XXX: GDALGetRasterXSize(psOptions->hSrcDS) is 895 and *pnSrcXOff is 0
XXX: GDALGetRasterXSize(psOptions->hSrcDS) - *pnSrcXOff = 895
XXX: ceil(dfMaxXOut) - *pnSrcXOff + nResWinSize = 476
XXX: pnSrcXSize = 476
XXX: pnSrcYSize = 103
XXX: ----- END OF SOURCE OF ZERO X SIZE -----
GDAL: Computed source window: 0,0,476x103
XXX: GDALWarpKernel ctor
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=0,0,476x103 Dst=0,0,512x128
XXX: ----- SOURCE OF ZERO X SIZE -----
XXX: GDALGetRasterXSize(psOptions->hSrcDS) is 895 and *pnSrcXOff is 451
XXX: GDALGetRasterXSize(psOptions->hSrcDS) - *pnSrcXOff = 444
XXX: ceil(dfMaxXOut) - *pnSrcXOff + nResWinSize = 483
XXX: pnSrcXSize = 444
XXX: pnSrcYSize = 168
XXX: ----- END OF SOURCE OF ZERO X SIZE -----
GDAL: Computed source window: 451,0,444x168
XXX: GDALWarpKernel ctor
GDAL: GDALWarpKernel()::GWKNearestNoMasksShort()
Src=451,0,444x168 Dst=512,0,512x128
XXX: ----- SOURCE OF ZERO X SIZE (3rd call of ProcessBlock) -----
XXX: GDALGetRasterXSize(psOptions->hSrcDS) is 895 and *pnSrcXOff is 895
XXX: GDALGetRasterXSize(psOptions->hSrcDS) - *pnSrcXOff = 0
XXX: ceil(dfMaxXOut) - *pnSrcXOff + nResWinSize = 497
XXX: pnSrcXSize = 0
XXX: pnSrcYSize = 228
XXX: ----- END OF SOURCE OF ZERO X SIZE -----
GDAL: Computed source window: 895,5,0x228
XXX: GDALWarpKernel ctor
ERROR 7: Assertion `0 != nSrcXSize' failed
in file `gdalwarpkernel.cpp', line 549

And this is full backtrace after assertion is caught (during 3rd call of VRTWarpedDataset::ProcessBlock function):

(gdb) bt
#0  0xffffe410 in __kernel_vsyscall ()
#1  0xb6e91df0 in raise () from /lib/tls/i686/cmov/libc.so.6
#2  0xb6e93641 in abort () from /lib/tls/i686/cmov/libc.so.6
#3  0xb7c004e3 in CPLErrorV (eErrClass=CE_Fatal, err_no=7, fmt=0xb7ea3aa8 "Assertion `%s' failed\nin file `%s', line %d\n", args=0xbfdfc3dc "�q���m��%\002") at cpl_error.cpp:233
#4  0xb7c00524 in CPLError (eErrClass=CE_Fatal, err_no=7, fmt=0xb7ea3aa8 "Assertion `%s' failed\nin file `%s', line %d\n") at cpl_error.cpp:133
#5  0xb7c0056f in _CPLAssert (pszExpression=0xb7ea71c8 "0 != nSrcXSize", pszFile=0xb7ea6df0 "gdalwarpkernel.cpp", iLine=549) at cpl_error.cpp:737
#6  0xb7c2f198 in GDALWarpKernel::PerformWarp (this=0xbfdfc508) at gdalwarpkernel.cpp:549
#7  0xb7c31333 in GDALWarpOperation::WarpRegionToBuffer (this=0x8098ec0, nDstXOff=1024, nDstYOff=0, nDstXSize=512, nDstYSize=128, pDataBuf=0x8144710, eBufDataType=GDT_UInt16, nSrcXOff=895, nSrcYOff=5, 
    nSrcXSize=0, nSrcYSize=228) at gdalwarpoperation.cpp:1561
#8  0xb7bb4449 in VRTWarpedDataset::ProcessBlock (this=0x8098d08, iBlockX=2, iBlockY=0) at vrtwarped.cpp:841
#9  0xb7bb4626 in VRTWarpedRasterBand::IReadBlock (this=0x8098de8, nBlockXOff=2, nBlockYOff=0, pImage=0x8124708) at vrtwarped.cpp:943
#10 0xb7bf2284 in GDALRasterBand::GetLockedBlockRef (this=0x8098de8, nXBlockOff=2, nYBlockOff=0, bJustInitialize=0) at gdalrasterband.cpp:1109
#11 0xb7bfa6d4 in GDALRasterBand::IRasterIO (this=0x8098de8, eRWFlag=GF_Read, nXOff=0, nYOff=0, nXSize=1126, nYSize=805, pData=0xb6bac008, nBufXSize=1126, nBufYSize=805, eBufType=GDT_UInt16, nPixelSpace=2, 
    nLineSpace=2252) at rasterio.cpp:212
#12 0xb7bd8139 in GDALDataset::IRasterIO (this=0x8098d08, eRWFlag=GF_Read, nXOff=0, nYOff=0, nXSize=1126, nYSize=805, pData=0xb6bac008, nBufXSize=1126, nBufYSize=805, eBufType=GDT_UInt16, nBandCount=1, 
    panBandMap=0xbfdfc948, nPixelSpace=2, nLineSpace=2252, nBandSpace=1812860) at gdaldataset.cpp:1313
#13 0xb7bd93d0 in GDALDataset::RasterIO (this=0x8098d08, eRWFlag=GF_Read, nXOff=0, nYOff=0, nXSize=1126, nYSize=805, pData=0xb6bac008, nBufXSize=1126, nBufYSize=805, eBufType=GDT_UInt16, nBandCount=1, 
    panBandMap=0xbfdfc948, nPixelSpace=2, nLineSpace=2252, nBandSpace=1812860) at gdaldataset.cpp:1506
#14 0xb7bf81e8 in GDALDatasetCopyWholeRaster (hSrcDS=0x8098d08, hDstDS=0x805d1b0, papszOptions=0x0, pfnProgress=0xb7bc698e <GDALTermProgress>, pProgressData=0x0) at rasterio.cpp:1882
#15 0xb79e418a in GTiffDataset::CreateCopy (pszFilename=0x8056860 "out.tif", poSrcDS=0x8098d08, bStrict=0, papszOptions=0x0, pfnProgress=0xb7bc698e <GDALTermProgress>, pProgressData=0x0) at geotiff.cpp:5449
#16 0xb7bde6bd in GDALDriver::CreateCopy (this=0x804f780, pszFilename=0x8056860 "out.tif", poSrcDS=0x8098d08, bStrict=0, papszOptions=0x0, pfnProgress=0xb7bc698e <GDALTermProgress>, pProgressData=0x0)
    at gdaldriver.cpp:558
#17 0xb7bde884 in GDALCreateCopy (hDriver=0x804f780, pszFilename=0x8056860 "out.tif", hSrcDS=0x8098d08, bStrict=0, papszOptions=0x0, pfnProgress=0xb7bc698e <GDALTermProgress>, pProgressData=0x0)
    at gdaldriver.cpp:597
#18 0x0804b6fb in ProxyMain (argc=3, argv=0x80568c0) at gdal_translate.cpp:607
#19 0x0804c2e3 in main (argc=Cannot access memory at address 0x48b5
) at gdal_translate.cpp:921

Frank, sorry I'm of little help in solving this. I've still not caught the GDAL warper internals enough, so I'd be thankful if you could help me here.

comment:8 by Even Rouault, 16 years ago

The problem is easier to reproduce in fact. No need to patch anything :

gdalwarp -of VRT ottawa_patch_small.tif ottawa_patch_small.vrt
gdal_translate ottawa_patch_small.vrt ottawa_patch_small_warped.tif
As for the fix, skipping GDALWarpKernel::PerformWarp when nSrcXSize == 0
nSrcYSize == 0 might be enough as Mateusz suggested. I've tried it and it seemed OK.

comment:9 by Mateusz Łoskot, 16 years ago

Even, thanks for complementing my analysis.

Frank, would you like to share your opinion?

comment:10 by warmerdam, 16 years ago

My concern is that I think PerformWarp() is still supposed to be executed even if there is no source data to ensure that output density maps and stuff get set. I'm not really sure about this though.

Actually, reviewing GWKGeneralCase() my concern is unfounded. Clearly the output density maps and stuff are already setup properly for a "no action" case before PerformWarp() is called. So bailing out of PerformWarp() when there is no source data should be just fine.

comment:11 by Mateusz Łoskot, 16 years ago

Owner: Mateusz Łoskot removed
Status: assignednew

comment:12 by warmerdam, 15 years ago

Cc: warmerdam removed
Owner: set to warmerdam

I have reviewed and determined that there is little magic here. I have just put an if around the computation of the scaling (used in filter sizing heuristics). The change is in trunk (r15820) and 1.5 branch (r15821).

comment:13 by warmerdam, 15 years ago

Resolution: fixed
Status: newclosed
Note: See TracTickets for help on using tickets.