[gdal-dev] Geolocation arrays - location interpretation

[Even Rouault]

On jeudi 28 septembre 2017 22:18:31 CEST Agram, Piyush S (334D) wrote:
> Hi Even,
>     I implemented the bilinear weighting like we discussed and that took out
> the systematic pixel shift from 1 to 0.5 (sign depending on orientation).
> You can see the changes here:
> https://github.com/piyushrpt/gdal/commit/b199cbae9ea15e31c787c1480c342bb84de
> bf774
 
> I also tried playing with the upsampling factor (1.3 to 4.0) and the
> discrepancy numbers were consistent after 1.5. Probably, the minimum needs
> to be sqrt(2) for square pixels.
 Using the weighting fixed the truncation
> / rounding error and the results seem consistent when I move the output
> extent around. 
> The systematic 0.5 pixel shift suggests that there is still an issue. Is my
> interpretation that the GeoLocTransformer is supposed to operate on pixel
> center coordinates correct?
 i.e – if my geoloc arrays were regular grid
> and I feed in the coordinates to the transformer, they return pixel centers
> – 0.5,0.5 etc. I’m going to try implementing some round trip checks to see
> if that shows something? 

Yes I agree there must be a 0.5 pixel shift issue when looking at the forward path in 
GDALGeoLocTransform(), ie the one taken by if( !bDstToSrc )

So we have currently

           const double dfGeoLocPixel =
                (padfX[i] - psTransform->dfPIXEL_OFFSET)
                / psTransform->dfPIXEL_STEP;
            const double dfGeoLocLine =
                (padfY[i] - psTransform->dfLINE_OFFSET)
                / psTransform->dfLINE_STEP;


For example take the case of offset = 0 and step = 1 to simplify

Imagine than padfX[0] = 0.5 and padfY[0] = 0.5 (so center of the top left pixel)

which will lead to iX = 0 and iY = 10

Now we are going to take a bilinear resampling at lines 788-800 (current trunk)

out_padfX[0] = 0.25 * (geoloc[0] + geoloc[1] + geoloc[width] + geoloc[width+1])

whereas it should be just out_padfX[0] = geoloc[0] . No resampling needed.

So perhaps something like the following ( substracting 0.5 ) would work ?

           const double dfGeoLocPixel =
                (padfX[i] - 0.5 - psTransform->dfPIXEL_OFFSET)
                / psTransform->dfPIXEL_STEP;
            const double dfGeoLocLine =
                (padfY[i] - 0.5 - psTransform->dfLINE_OFFSET)
                / psTransform->dfLINE_STEP;

and probably something equivalent in the other direction

Perhaps we should have some metadata item to indicate the convention and apply or not this 
0.5 shift
Like GEOLOC_COORD_CONVENTION=PIXEL_CENTER (default if not specified) or 
GEOLOC_COORD_CONVENTION=PIXEL_TOP_LEFT
But I'd expect most (all ?) datasources to use a pixel_center convention.

Which format/driver are you using ?

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20170929/698b0f61/attachment.html>