[Gdal-dev] GDAL GCPsToGeoTransform Problem

[Eric Dönges]

Am 17.02.2005 um 15:12 schrieb Frank Warmerdam:

> On Thu, 17 Feb 2005 08:34:51 -0500, David Legault
> <dlegault at vantpoint.com> wrote:
>> then use the gdal.GCPsToGeoTransform( gcps, 1 ) function to retreive 
>> the
>> calculated georeferencing data and I always get this
>>
>> (-1.#IND, -1.#IND, 1.#INF, -1.#IND, -1.#IND, 1.#INF)
>>
>> which tends to make me think there is something wrong with the 
>> calculation
>> routine or the input parameters are not good (or any other reason).
>>
>> But if I choose 3 points in all those (any 3 not aligned) the 
>> function does
>> return some data (which I didn't check to see if it was correct or 
>> not)

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.

With kind regards,
Eric