[PATCH] RPC transform: pixel center vs pixel corner

[pablo]

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:

1. 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.

2. 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.

[warmerdam]

In regard to ​http://lists.osgeo.org/pipermail/gdal-dev/2015-June/042087.html, I just wanted to take an alternate position to Evens, and go with what I think Pablo is suggesting.

I feel that the RPC coefficients have well establish meanings from the NITF spec, and file formats like _rpc.txt. I assume they are center-of-pixel oriented. I would *not* want the RPC metadata we keep (ie ​https://trac.osgeo.org/gdal/wiki/rfc22_rpc) to have a different meaning for the pixel/line locations. So I would suggest we should not need to transform the RPC coefficients when they are imported - instead it is the evaluator which needs to adapt between the RPC pixel/line model and the usual GDAL interpretation.

[Even Rouault]

I've just a bit rearranged fix_rpc_dem_interpolation.patch to avoid doing -0.5 and then +0.5 in the nearest neighbour case. gdal_datamodel.doc and RFC 22 edited to clarify the meaning of LINE_OFF and SAMP_OFF.

trunk r29416, branches/2.0 r29417 "RPC: fix off-by-half pixel computation of (pixel, line), and in bilinear and bicubic RPC DEM interpolation; fix off-by-one pixel registration for Pleaiades RPC (#5993, derived from patches by Pablo d'Angelo)"

trunk r29418, branches/2.0 r29419 "Adjust vrtwarp_8 expected results due to #5993"

[Even Rouault]

trunk r29420 (so only GDAL 2.1) "gdalwarp: if RPC_DEM warping option is specified, use exact transformer by default (#5993)"