[gdal-dev] Geolocation arrays - location interpretation

[Agram, Piyush S (334D)]

Hi Even,
     Here is my first cut at fixing the half pixel shifts:
https://github.com/piyushrpt/gdal/commit/ead79550e7120767187a0b9c836b5d861223d997

With these changes, the shifts are down to sub 0.2 pixels. Still not sure if I got the math right. Feel like the errors should be smaller for the ideal simulated datasets that I’m using.

Piyush

From: Even Rouault <even.rouault at spatialys.com>
Date: Thursday, September 28, 2017 at 3:36 PM
To: "Agram, Piyush S (334D)" <Piyush.Agram at jpl.nasa.gov>
Cc: "gdal-dev at lists.osgeo.org" <gdal-dev at lists.osgeo.org>
Subject: Re: [gdal-dev] Geolocation arrays - location interpretation


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/5b21caaf/attachment-0001.html>