[PATCH] RPC transform: pixel center vs pixel corner
Hello GDAL developers,
I have noticed that the output of the RPC based gdalwarp does not match the output of our internal rpc orthorectification tools.
I thus digged into the source code of gdal 2.0.0 rc1 and found several places where I think the "upper left corner of upper left pixel is 0,0" convention used by the rest of GDAL is not properly handled:
- When using the RPC_DEM option, DEM row, column values are computed
using the geotransform, but the interpolation weights for bicubic interpolation are computed as (dX - int(dX)), which is valid for "origin is at center of pixel" convention. My understanding of the warp kernel code is that it computes the bilinear interpolation weights with "origin is at center of pixel" convention. I have changed the gdal_rpc.cpp code to work in a similar way.
Interestingly the "nearest neighbour" was implemented correctly, so I had to adjust that to fit these changes, too. This was because a simple floor(x_in_ul_corner_coords) is the correct case for nearest neighbour.
I am not so certain about the cubic convolution approach, as I our rpc ortho tool does not implement cubic convolution elevation interpolation, so I cannot compare against another independent implementation. I tried comparing the gdalwarpkernel.cpp implementation in the warp kernel with that in gdal_rpc.cpp, but they use different formulation, and it is not obvious if they compute the interpolation in the same way.
When comparing the images orthorectified with bilinear and cubic convolution DEM interpolation, one can see some smaller, local changes, but the result is quite similar, so I believe the coordinate origin changes are appropriate for cubic convolution, too.
The patch fix_rpc_dem_interpolation.patch contains these changes. This is especially important when using coarse resolution DEMs like SRTM when orthorectifying images with larger off nadir angles in hilly regions.
- The Digital globe RPC coefficients use the "center of upper left
pixel is 0,0" convention, and Pleiades (and SPOT-6,7) uses "center of upper left pixel is 1,1" convention, see https://www.orfeo-toolbox.org/otb-4-2-and-monteverdi2-0-8-are-out/ The gdal RPCTransformPoint() function does not consider this, and returns the image coordinates in "center of pixel" coordinates, instead of the "upper left corner coordinates" expected by the rest of GDAL.
I assume that the DG convention is most commonly used, and thus add 0.5 to the pixel coordinates in RPCTransformPoint, see fix_RPCTransformPoint.patch
The fix_PleiadesRPC.patch does subtract 1 pixel from the RPC offsets, when reading them using the GDALMDReaderPleiades::LoadRPCXmlFile() reader, so that the RPCs used by GDAL follow the Digital Globe and NITF conventions.
I had to adjust the tests in gdalautotest to deal with these changed convention, see the gdalautotest_rpc_corner_changes.patch.
After these patches, the output of gdalwarp fits our internal RPC orthorectification tool very well. A difference images shows a few pixels with a difference of -1 and 1, but that is likely due to slight rounding and floating point computation issues.
Additionally, I noticed that -et 0 is required for good orthorectification, as the default of -et=0.1 leads to strange artifacts with errors of more than a few pixels in hilly areas. Probably the error detection cannot handle irregular transformations well. I wonder if gdalwarp should default to -et=0 if RPC_DEM is set.
Change History (8)
comment:1 by , 9 years ago
|RPC transform: pixel center vs pixel corner → [PATCH] RPC transform: pixel center vs pixel corner