[Gdal-dev] gdalwarp: difference between -tps and -order 1

[Frank Warmerdam]

Jan Hartmann wrote:
> Hi,
> 
> I'm georeferencing historical maps that are not geometrically correct. 
> When I warp them with control points to a modern map, using a first 
> order rigid transformation (gdalwarp -order 1), the resulting, 
> georeferenced map will keep its original form. This means that the 
> position of the control points on the georeferenced historical map is 
> different from those on the modern map. Warping with  a thin plate 
> spline (gdalwarp -tps) results in a georeferenced map with the control 
> points exactly on the same place on both maps, which is what I want of 
> course.
> 
> I can see the difference between both methods visually, but I would like 
> also to obtain the numerical coodinates of the control points after 
> warping them with both methods. That way I can compute the differences 
> between both methods and construct a distortion map. So essentially I 
> would like to have a warping program that does not output a raster, but 
> produces the coordinates of the transformed input control points. I 
> already did this for the first order transformation (pretty elementary), 
> but how could this be done for the thin plate spline transformation? 
> Where should I look in the gdal codebase to extract the algorithm? 
> Alternatively, if an expert would program this, how much would that cost?

Jan,

In C++ you can construct your "TPSTransformer" and use it to transform
individual points.  But currently this is not available from command
line utilities or the swig bindings.

This is an example of instantating a TPS transformer:

             psInfo->pSrcTPSTransformArg =
                 GDALCreateTPSTransformer( GDALGetGCPCount( hSrcDS ),
                                           GDALGetGCPs( hSrcDS ), FALSE );

And this is an example of actually transforming a point with it.

         if( !GDALTPSTransform( pTPSTransformArg, FALSE,
                                nPointCount, padfX, padfY, padfZ,
                                panSuccess ) )
             return FALSE;

More complete code is available in:

   gdal/alg/gdaltransformer.cpp
   gdal/alg/gdal_tps.cpp

For a compatriot from the old country I would be willing to prepare
a small command line program for free if we can work out the specifications
offline.

At some point I'd like to expose more of the warper and transformer APIs
via swig, but they involve a lot of callbacks and idiosyncratic data
structures.

Best regards,
-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | President OSGeo, http://osgeo.org