[gdal-dev] GSoC Image Correlator

[Even Rouault]

Selon Frank Warmerdam <warmerdam at pobox.com>:

> Folks,
>
> Андрей Мигаль (Andrew) has completed his Google Summer of Code
> image
> correlator project and I have moderately adapted the code and checked
> into GDAL trunk svn.  Some information is available at:
>
>    http://trac.osgeo.org/gdal/wiki/Correlator
>    http://correlatorgsoc2012.blogspot.com/
>    https://github.com/migal-drew/GDAL-correlator
>    http://trac.osgeo.org/gdal/changeset/24868
>
> While I have exposed a function with the following signature, I have
> not yet checked in any changes that would use this from an existing
> commandline program.  My plan is to expose it via SWIG and write some
> demonstration Python scripts using it:
>
> GDAL_GCP CPL_DLL *
> GDALComputeMatchingPoints( GDALDatasetH hFirstImage,
>                             GDALDatasetH hSecondImage,
>                             char **papszOptions,
>                             int *pnGCPCount );
>
>
> I have successfully use it to perform an image match between a UTM image
> and the same image warped into WGS84.  I then wrote out the GCPS and use
> them to warp the the one to the other and got a good result.
>
> I will also document the GDALComputeMatchingPoints() function based on the
> extensive documentation already prepared by Andrew for lower level functions.
>
> Andrew has done a great job this summer, and I hope we will continue to
> evolve and exploit the new functionality developed.

I've done a quick review of the code and indeed it is impressively
documented ! I guess I should find some paper that explains the maths behind for
a better understanding. Perhaps Andrew has someuseful links ?

Looking at gdalmatching.cpp, I'm wondering if GDALComputeMatchingPoints()
shouldn't operate on bands directly rather than on datasets. Because I see that
the first processing step is to compute a luminosity buffer from the R,G,B
bands. This could perhaps be done in a different API, like :

GDALRasterBandH GDALCreateLuminosityBandFromRGB(GDALDatasetH hDS);

void GDALFreeLuminosityBand(GDALRasterBandH );

(could be implemented internally by a VRT dataset that uses a custom pixel
function)

That you would call before GDALComputeMatchingPoints(), and that you will not
need to call if you operate directly on single band datasets.

I've noted that GDALSimpleSURF::ConvertRGBToLuminosity() does a normalization of
the values to [0,1] by dividing by 255, so that assumes that the input range is
[0,255]. Is that a strong requirement that the values are in the [0,1] range
after ConvertRGBToLuminosity() ? I guess there are use cases where satellite
imagery could be in other datatypes, and that dividing by 255 isn't the
appropriate normalization value. But perhaps the normalization must still be
done in order to be able to compare properly the 2 images.

so the API could be :

GDAL_GCP CPL_DLL *
GDALComputeMatchingPoints( GDALRasterBandH hFirstBand,
                           double dfMinFirstBand, double dfMaxFirstBand,
                           GDALDatasetH hSecondImage,
                           double dfMinSecondBand, double dfMaxSecondBand,
                           char **papszOptions,
                           int *pnGCPCount )

dfMinFirstBand and dfMaxFirstBand would be used to do the normalization to
[0,1]. (We could deduce them from statistics, but sometimes you perhaps want to
clamp input values because the extremes are not significant : error measures,
etc...).

(We could argue that the normalization could be done before calling
GDALComputeMatchingPoints(), but that's perhaps not necesserary).

The other important point is that the algorithms assume that each image can
entirely fits in RAM. I've the feeling by looking at gdal_octave.cpp that the
maths done effectively imply to have the full image (or that modifying it to be
able to operate on regions might require a significant rewrite). So we should
probably at least catch std::bad_alloc exceptions in GDALComputeMatchingPoints()

Best regards,

Even