Geoloc transformer: fix inverse transform to be exact #5520
Conversation
|
CC @piyushrpt as you've worked in the past in that area of the code |
|
The GDAL project highly values your contribution and would love to see this work merged! Unfortunately this PR has not had any activity in the last 21 days and is being automatically marked as "stale". If you think this pull request should be merged, please check - that all unit tests are passing - that all comments by reviewers have been addressed - that there is enough information for reviewers, in particular
|
|
Thanks @rouault - you are right about |
|
Thanks again for all the updates on this. Only comment is that |
Despite a few past attempts to improve it, the inverse transformer relying on a backmap was far from being exact, with errors frequently above 1 pixel. This commit fixes it by using an (exact, non-interative) inverse bilinear interpolation method (https://stackoverflow.com/a/812077) during backmap generation itself (when the geolocation array expresses an affine transformation, the backmap itself is exact), but also when inverse transforming a coordinate, to fine tune the initial solution obtained by bilinear interpolation of the backmap. The FSHIFT and ISHIFT offsets, that were probably more empirical than are no longer needed. The oversampling factor [0.1,2] (default: 1.3) of the backmap can now be controlled with the GEOLOC_BACKMAP_OVERSAMPLE_FACTOR transformer option (or the GDAL_GEOLOC_BACKMAP_OVERSAMPLE_FACTOR configuration option). I also noticed that the referencing convention of the geolocation array was not clearly documented. Looking at the forward method, the existing convention is a top-left corner one. This means that the right-most column and bottom-most column of a raster (assuming a 1:1 match between the raster and geolocation array) are not within in a cell of the geolocation array where we can naturally apply bilinear interpolation, thus interpolation from the nearest cell was applied. For a number of drivers, the georeferencing will be more using a pixel-center convention, hence we add support in the GEOLOCATION metadata domain for a GEOREFERENCING_CONVENTION item that defaults to TOP_LEFT_CORNER but can be set to PIXEL_CENTER (and we set it to PIXEL_CENTER in the netCDF driver where this is the most likely convention). The forward and inverse transformers take into account that to apply a half-pixel shift to go back/from the top-left corner convention used internally. There are also quite a lot of precaution to deal with cells of the geolocation array that span the antimeridian, to avoid erroneous interpolation to be done (we don't want averaging longitudes close to -180 deg with ones close to +180 deg), or missed pixels. Those cases are particular triggered by the below mentionned dataset that includes the north Pole and thus has discontinuity in longitudes in the geolocation array around +/- 180 deg. There's a performance cost for those improvements. gdalwarp applied on a Sentinel 5P product (S5P_TEST_L2__NO2____20190509T220707_20190509T234837_08137_01_010400_20200220T091343) of size 450x3245 goes from 1.26 s to 3.52 s. This commit also includes an alternate method for the inverse transform that no longer uses a backmap, but a quadtree indexing the quadrilaterals (in georeferenced coordinates) of each cell of the geolocation array. It gives even more exact results than the above improvements (no issue with holes and hole filling at high latitudes), and make the code simpler, but it uses more RAM and is slower: 8.14 s on the above mentioned test case. It can be selected by setting the GDAL_GEOLOC_INVERSE_METHOD configuration option to QUADTREE.
I've ammended the commit to document the option in GDALCreateGenImgProjTransformer2() close to other transformer options |
I believe that with the new method even if the backmap is stored on float, the refinement step (exact inverse bilinear interpolation) is done with double precision, so the inverse geoloc tranformer should have double precision accuracy. |
…loc array is too big Add a GDAL_GEOLOC_USE_TEMP_DATASETS=YES/NO config option (GEOLOC_USE_TEMP_DATASETS transformer option) to control whether to use the temporary GTiff or the in-memory C array. Temporary GTiff is used by default if the geolocation array is larger than 16 megapixels Adds also a few speed optimizations to recover some of the performance drop of the exact inverse bilinear transformer. Runtime of the S5P_TEST_L2__NO2____20190509T220707_20190509T234837_08137_01_010400_20200220T091343 product is 2.8 s with the in-memory C array and 3.5 s for the temporary GTiff dataset.
|
New commit added: |
|
I'm slowly picking through this one @rouault 😆 |
|
@nyalldawson will you have some time to complete your review ? I'd like this to go in 3.5.0 whose release candidate is due this week |
Despite a few past attempts to improve it, the inverse transformer
relying on a backmap was far from being exact, with errors frequently
above 1 pixel.
This commit fixes it by using an (exact, non-interative) inverse bilinear
interpolation method (https://stackoverflow.com/a/812077) during backmap
generation itself (when the geolocation array expresses an affine
transformation, the backmap itself is exact), but also when inverse
transforming a coordinate, to fine tune the initial solution obtained
by bilinear interpolation of the backmap.
The FSHIFT and ISHIFT offsets, that were probably more empirical than
are no longer needed.
The oversampling factor [0.1,2] (default: 1.3) of the backmap can now be
controlled with the GEOLOC_BACKMAP_OVERSAMPLE_FACTOR transformer option
(or the GDAL_GEOLOC_BACKMAP_OVERSAMPLE_FACTOR configuration option).
I also noticed that the referencing convention of the geolocation array
was not clearly documented. Looking at the forward method, the existing
convention is a top-left corner one. This means that the
right-most column and bottom-most column of a raster (assuming a 1:1
match between the raster and geolocation array) are not within in a cell
of the geolocation array where we can naturally apply bilinear interpolation,
thus interpolation from the nearest cell was applied.
For a number of drivers, the georeferencing will be more using a pixel-center
convention, hence we add support in the GEOLOCATION metadata domain for a
GEOREFERENCING_CONVENTION item that defaults to TOP_LEFT_CORNER but can be
set to PIXEL_CENTER (and we set it to PIXEL_CENTER in the netCDF driver
where this is the most likely convention). The forward and inverse
transformers take into account that to apply a half-pixel shift to go
back/from the top-left corner convention used internally.
There are also quite a lot of precaution to deal with cells of the geolocation
array that span the antimeridian, to avoid erroneous interpolation to be
done (we don't want averaging longitudes close to -180 deg with ones
close to +180 deg), or missed pixels. Those cases are particular triggered
by the below mentionned dataset that includes the north Pole and thus has
discontinuity in longitudes in the geolocation array around +/- 180 deg.
There's a performance cost for those improvements.
gdalwarp applied on a Sentinel 5P product (S5P_TEST_L2__NO2____20190509T220707_20190509T234837_08137_01_010400_20200220T091343)
of size 450x3245 goes from 1.26 s to 3.52 s.
This commit also includes an alternate method for the inverse transform
that no longer uses a backmap, but a quadtree indexing the
quadrilaterals (in georeferenced coordinates) of each cell of the
geolocation array. It gives even more exact results than the above
improvements (no issue with holes and hole filling at high latitudes),
and make the code simpler, but it uses more RAM and is slower: 8.14 s on
the above mentioned test case. It can be selected by setting the GDAL_GEOLOC_INVERSE_METHOD
configuration option to QUADTREE.