[gdal-dev] Overview scales

[Even Rouault]

Le mercredi 06 juillet 2016 12:52:15, Wilfried Karel a écrit :
> Dear Even,
> 
> I'm not referring to a specific interpolation method. However, for
> checking GDAL's behaviour, nearest-neighbor interpolation is surely the
> easiest way.

I guess the effect is the most distracting with nearest neighbor, as other 
methods do interpolations and thus you cannot map a target pixel to a single 
source pixel, and the effect of taking a non exact sampling factor has less 
visible consequences as the contribution of fractional source pixel is well 
taken into account (at least for those using the generic convolution 
infrastructure: cubic, cubicspline, bilinear and lanczos. Wouldn't be the case 
for Gauss, which apparently due to recent tickets may have its own problems, 
or Mode / Average)

> And yes, it may seem arbitrary whether to select only even
> rows/columns or only odd rows/columns of the base image when
> downsampling by a factor of 2 (although I'd prefer to select even rows,
> as that means that the image coordinate systems of overviews only differ
> by scales, and not also by shifts - if indices are 0-based and if the
> image coordinate systems' origins are at the center of the upper/left
> pixel).
> 
> However, what seems not arbitrary to me, but is quite an obstacle for me
> (and possibly others who would like to use GDAL for doing
> photogrammetry), is GDAL's behavior to not use strictly constant
> sampling distances, as can be seen in the following example of 1 line of
> a base image with 11 columns:
> 
> 246   7   2 212  93  80 231 174 134  69  24
> 
> downsampling by a decimation factor of 2, using nearest-neighbor
> interpolation (via GDALDataset::BuildOverviews) results in:
> 
> 246       2      93     231 174      69
> 
> 2 neighboring pixels in the base image are both mapped to the overview,
> while all other source pixels are separated by a non-mapped source
> pixel. That means that there is no linear transformation between image
> coordinates in those 2 images. What I would strongly welcome, at least
> for the very common decimation factors of 2, would be the following line
> in the overview:
> 
> 246       2      93     231     134      24
> 
> In that case, when transforming image coordinates from the overview to
> the base image or vice versa, simply the scale factor of 2 would need to
> be considered.
> 
> Looking into GDAL, it seems to me that (relative) image scales are
> always derived from the base image / overview resolutions. While this
> may be a good way to handle overviews with arbitrary scales (probably
> created with other tools), I find this very distracting when the goal is
> to simply have scale factors of exactly 2.

OK I see your point.

> While decimation factors for GDALDataset::BuildOverviews can be
> specified as integers only, the actual scale ratios may result as e.g.
> 1.8999 or so. This ratio is applied and then casted to
> integer-coordinates of the source image grid.

A criterion could be that the rounded scale ratio shouldn't account for a 
difference of more than one pixel. Something like :

dfFactor = (float)nFullResSize / nOvrSize
int nClosestIntFactor = (int)(dfFactor + 0.5)
if( ABS(nOvrSize *  nClosestIntFactor - nFullResSize) <= 1 )
	dfFactor = nClosestIntFactor

For non nearest neighbour methods, some care should be taken to ensure that 
for the pixels at the right and bottom of the image, all the source pixels are 
well taken into account. I think that this must be already the case, but 
should be checked.

> Hence the 2 neighboring
> rows/columns close to the horizontal and vertical center lines of the
> source image that are both mapped to the overview image.
> 
> I think that scale factors of 2 are the by far most common case. Thus,
> I'd vote for rounding scale factors that are "close" to 2 to exactly 2.
> However, as I can see it, the derivation of scales by dividing the
> resolutions of 2 images is done in multiple places.
> Hence, this change
> would need to be done in all those places, I believe.

In gcore/rasterio.cpp, in GDALRasterBand::IRasterIO() there's similar logic as 
well for the nearest resampling logic. Other methods end up using overview 
building internally (but there would be likely changes here to do in factor 
computation)

And we would perhaps want the warper alg in alg/gdalwarpkernel.cpp to be 
updated as well. Not sure how much feasible it is however. Would only make 
sense in the cases the warper is used for simple resampling, ie that there's 
no reprojection or GCP/RPC fitting, but the warper kernel doesn't know the 
properties of the src <--> target pixel transformer. It is just a blackbox for 
it. Hum, we could probably skip the warper update, as using it for just 
resampling is overkill since GDAL 2.0 gdal_translate supports the -r 
parameter.

So yes there's a bunch of code change, and also test suite "golden" results to 
be updated. That'd need a champion.

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com