[gdal-dev] Overview scales

[Wilfried Karel]

Dear Even,

On 06.07.2016 13:45, Even Rouault wrote:
> 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)
Yes, with interpolation methods that take into account more than a
single pixel of the source image, this effect will not be that obvious.
But that doesn't resolve the problem if one is interested in more than
just visualization, but e.g. in deriving positions from the image
content, as in least-squares image matching at decreasing image pyramid
levels and/or between images with largely different scales.

> 
>> 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.
Unfortunately, I'm surely no GDAL champion. So I hope that one of those
picks up this problem!

Cheers,
Willi.
> 
> Even
> 

-- 
                        _____ _____ _____  Technische Universität Wien
   DI Wilfried Karel   /____// ___//    / Dept. of Geodesy & Geoinform.
  wilfried.karel@     // __ / /__ / // / Research Groups Photogrammetry
 geo.tuwien.ac.at    //__/// /__ / // /   & Remote Sensing
+43 1 58801 12244   /____//____//____/ Gusshausstr. 27-29, A-1040 Wien