[Gdal-dev] GDAL GCPsToGeoTransform Problem

[Eric Dönges]

Am 17.02.2005 um 16:52 schrieb Gillian Walter:

> Eric Dönges wrote:
>> While I am reasonably sure the GCPsToGeoTransform code is 
>> mathematically correct (the
>> equations were determined with a symbolic algebra package), it is 
>> probably not optimal
>> from a numerical standpoint. The problem is that if you have lots of 
>> GCPs with large
>> coordinates, the calculations may overflow 'mere' 64 bit double 
>> arithmetic (not being
>> a mathematician, I wrote the code for clarity, not numerical 
>> stability). Looking
>> at the equations, I do not see an obvious solution for avoiding the 
>> rather large
>> temporary values during the calculation.
>
> If large intermediate values are the problem, maybe the gcps could be 
> shifted by an offset, and then shifted back?  ie. subtract mean(pixel 
> values), mean(line values), mean(x values), mean(y values) from each 
> gcp, do a fit to the shifted gcps so you'd end up with (where gt is 
> the geotransform):
>
> x - xmean = gt[0] + gt[1](p - pmean) + gt[2](l - lmean)
> y - ymean = gt[3] + gt[4](p - pmean) + gt[5](l - lmean)
>
> then the real gt[0] and gt[3] would be:
>
> gt[0] actual = gt[0] + xmean - gt[1]*pmean - gt[2]*lmean
> gt[3] actual = gt[3] + ymean - gt[4]*pmean - gt[5]*lmean
>
> or something like that.   The fit code seems to have a lot of sums of 
> squares of pixel/line/x/y values, and this should make them smaller.  
> You might also be able to do some kind of normalization by spread to 
> make them even smaller if necessary.
> Gillian

I'm not sure if this will help with floating point arithmetic, since 
the way I
understand floating point arithmetic (and I may be completely wrong 
here !), the
problem is not absolute magnitude (IEEE doubles go to something like 
10^108 or so,
I believe), but relative magnitude between numbers. In our case, we 
would need to
loose the products of pixel/line/x/y values, i.e. try to rework the 
equations so that
the intermediate results are not several orders of magnitude larger or 
smaller than
the input values. A possibly simpler solution might be to simply switch 
to using
'long double'; however I am not sure how well this data type is 
supported on
different compilers.

Another thought I had is that the original poster 
(<dlegault at vantpoint.com>) should
check the return value of GCPsToGeoTransform - it may be that the 
complete set
of GCPs does not result in a valid transform (which is possible if some 
of the GCPs
are poorly placed) and thus the function returns gibberish.

With kind regards,
Eric