Index: ogr2ogr.cpp
===================================================================
RCS file: /cvs/maptools/cvsroot/gdal/ogr/ogr2ogr.cpp,v
retrieving revision 1.28
diff -u -r1.28 ogr2ogr.cpp
--- ogr2ogr.cpp	14 Apr 2005 14:20:24 -0000	1.28
+++ ogr2ogr.cpp	6 Jul 2005 04:13:58 -0000
@@ -117,6 +117,10 @@
 #include "ogrsf_frmts.h"
 #include "cpl_conv.h"
 #include "cpl_string.h"
+#include <gdal_alg.h>
+#include <gdal.h>
+#include <stdio.h>
+
 
 CPL_CVSID("$Id: ogr2ogr.cpp,v 1.28 2005/04/14 14:20:24 fwarmerdam Exp $");
 
@@ -133,10 +137,16 @@
                            char **papszSelFields,
                            int bAppend, int eGType );
 
+static void *CreateGCPTransform (const char *file);
+
+static OGRGeometry *warpGeometry (
+    OGRGeometry *poGeometry, void *poTransform );
+
 static int bSkipFailures = FALSE;
 static int nGroupTransactions = 200;
 static int bPreserveFID = FALSE;
 static int nFIDToFetch = OGRNullFID;
+static void *pTPSxform = NULL;
 
 /************************************************************************/
 /*                                main()                                */
@@ -303,6 +313,10 @@
             papszSelFields = CSLTokenizeStringComplex(pszSelect, " ,", 
                                                       FALSE, FALSE );
         }
+        else if( EQUAL(papszArgv[iArg],"-gcp") && papszArgv[iArg+1] != NULL )
+        {
+            pTPSxform = CreateGCPTransform(papszArgv[++iArg]);
+        }
         else if( papszArgv[iArg][0] == '-' )
         {
             Usage();
@@ -806,6 +820,15 @@
             }
         }
 
+	if (pTPSxform && poDstFeature->GetGeometryRef() != NULL ) {
+	    OGRGeometry *poDstGeom = poDstFeature->StealGeometry();
+	    OGRGeometry *poDstGeom2;
+
+	    poDstGeom2 = warpGeometry( poDstGeom, pTPSxform );
+	    poDstFeature->SetGeometryDirectly(poDstGeom2);
+	}
+
+
         if( poDstFeature->GetGeometryRef() != NULL && bForceToPolygon )
         {
             poDstFeature->SetGeometryDirectly( 
@@ -842,3 +865,113 @@
     return TRUE;
 }
 
+static void warpPoint (OGRPoint *poPoint, void *pXformer)
+{
+	double x = poPoint->getX(),
+	       y = poPoint->getY(),	
+	       z = poPoint->getZ();
+	int r = 0;
+	// printf("%lf %lf %lf -> ", x, y, z);
+	GDALTPSTransform(pXformer, FALSE, 1, &x, &y, &z, &r );
+	// printf("%lf %lf %lf\n", x, y, z);
+	if (!r) {
+	    printf("transform failed!\n");
+	    return;
+	}
+	poPoint->setX(x);	
+	poPoint->setY(y);	
+	poPoint->setZ(z);	
+}
+
+static void warpCurve (OGRCurve *poOut, OGRCurve *poIn, void *pXformer)
+{
+    OGRPoint *poPoint = new OGRPoint();
+    OGRLineString *poIn2 = (OGRLineString *) poIn;
+    OGRLineString *poOut2 = (OGRLineString *) poOut;
+    for (int n = 0; n < poIn2->getNumPoints(); n++) {
+	poIn2->getPoint(n, poPoint);
+	warpPoint(poPoint, pXformer);
+	poOut2->addPoint(poPoint);
+    }
+    delete poPoint;
+}
+
+static OGRGeometry *warpGeometry (OGRGeometry *poGeometry, void *poTransform)
+{
+    OGRGeometry *poGeomOut = NULL;
+    OGRwkbGeometryType eGeomType = poGeometry->getGeometryType();
+
+    if (eGeomType == wkbPoint) {
+	OGRPoint *poPoint = (OGRPoint *) poGeometry->clone();
+	warpPoint(poPoint, poTransform);
+	poGeomOut = (OGRGeometry *)poPoint;
+    }
+    if (eGeomType == wkbLineString) {
+	OGRLineString *poLineIn = (OGRLineString *) poGeometry;
+	OGRLineString *poLineOut = new OGRLineString();
+	warpCurve(poLineOut, poLineIn, poTransform);
+	poGeomOut = (OGRGeometry *)poLineOut;
+    }
+    else if (eGeomType == wkbPolygon) {
+	OGRPolygon *poPolyIn = (OGRPolygon *) poGeometry;
+	OGRPolygon *poPolyOut = new OGRPolygon();
+	OGRLinearRing *poRingIn = poPolyIn->getExteriorRing();
+	OGRLinearRing *poRingOut = new OGRLinearRing();
+
+	//printf("-.");
+	warpCurve( poRingOut, poRingIn, poTransform );
+	poPolyOut->addRingDirectly( poRingOut );
+	//printf(".-");
+
+	for (int n = 0; n < poPolyIn->getNumInteriorRings(); n++) {
+	    poRingIn = poPolyIn->getInteriorRing(n);
+	    poRingOut = new OGRLinearRing();
+	    warpCurve( poRingOut, poRingIn, poTransform );
+	    poPolyOut->addRingDirectly( poRingOut );
+	}
+
+	poGeomOut = (OGRGeometry *)poPolyOut;
+
+    }
+    return poGeomOut;
+}
+
+static void *CreateGCPTransform (const char *file) {
+    int n = 0, r = 0;
+    GDAL_GCP gcps[32];
+    FILE *gcpfile;
+    void *poTransform;
+
+    gcpfile = fopen(file, "r");
+    if (!gcpfile) {
+	printf("Can't open GCP file %s.\n", file);
+	exit(0);
+    }
+
+    while (!feof(gcpfile)) {
+	r = fscanf(gcpfile, "%lf %lf %lf %lf", 
+		&(gcps[n].dfGCPX),	&(gcps[n].dfGCPY),
+		&(gcps[n].dfGCPPixel),	&(gcps[n].dfGCPLine) );
+	if (r < 4)
+	    continue;
+/*	printf( "GCP %lf %lf %lf %lf\n", 
+		gcps[n].dfGCPX,	gcps[n].dfGCPX,
+		gcps[n].dfGCPPixel,	gcps[n].dfGCPLine ); */
+	gcps[n].dfGCPZ = 0;
+	gcps[n].pszId = gcps[n].pszInfo = "";
+	if (n++ == 32) {
+	    printf("only 32 GCPs supported, sorry :-/\n");
+	    break;
+	}
+    }
+    fclose(gcpfile);
+
+    poTransform = GDALCreateTPSTransformer( n, gcps, 0 );
+    if (poTransform == NULL) {
+	printf("Couldn't generate thin plate spline transform from GCPs.\n");
+	exit(1);
+    }
+
+    return poTransform;
+}
+