[Gdal-dev] GDAL GCPsToGeoTransform Problem

[Gillian Walter]

Eric Dönges wrote:

>
> 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.

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

>
> With kind regards,
> Eric
>
> _______________________________________________
> Gdal-dev mailing list
> Gdal-dev at xserve.flids.com
> http://xserve.flids.com/mailman/listinfo/gdal-dev