i.rectify: add thin plate spline coordinate transformation

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@56566 15284696-431f-4ddb-bdfa-cd5b030d7da7

imagery/i.rectify/cp.c

@@ -7,7 +7,6 @@ int get_control_points(char *group, int order	/* THIS HAS BEEN ADDED WITH THE CR
     )
 {
     char msg[200];
-    struct Control_Points cp;
 
     if (!I_get_control_points(group, &cp))
 	exit(0);
@@ -15,30 +14,58 @@ int get_control_points(char *group, int order	/* THIS HAS BEEN ADDED WITH THE CR
     sprintf(msg, _("Control Point file for group <%s@%s> - "),
 	    group, G_mapset());
 
-    switch (I_compute_georef_equations(&cp, E12, N12, E21, N21, order)) {
-    case 0:
-	sprintf(&msg[strlen(msg)],
-		_("Not enough active control points for current order, %d are required."),
-		(order + 1) * (order + 2) / 2);
-	break;
-    case -1:
-	strcat(msg, _("Poorly placed control points."));
-	strcat(msg, _(" Can not generate the transformation equation."));
-	break;
-    case -2:
-	strcat(msg, _("Not enough memory to solve for transformation equation"));
-	break;
-    case -3:
-	strcat(msg, _("Invalid order"));
-	break;
-    default:
-	/* COMMENTED OUT WHEN SUPPORT FOR 3rd ORDER WAS ADDED BY 'CRS'
-	   E12a = E12[0]; E12b = E12[1]; E12c = E12[2];
-	   N12a = N12[0]; N12b = N12[1]; N12c = N12[2];
-	   E21a = E21[0]; E21b = E21[1]; E21c = E21[2];
-	   N21a = N21[0]; N21b = N21[1]; N21c = N21[2];
-	 */
-	return 1;
+    if (order == 0) {
+	switch (I_compute_georef_equations_tps(&cp, &E12_t, &N12_t, &E21_t, &N21_t)) {
+	case 0:
+	    sprintf(&msg[strlen(msg)],
+		    _("Not enough active control points for thin plate spline."));
+	    break;
+	case -1:
+	    strcat(msg, _("Poorly placed control points."));
+	    strcat(msg, _(" Can not generate the transformation equation."));
+	    break;
+	case -2:
+	    strcat(msg, _("Not enough memory to solve for transformation equation"));
+	    break;
+	case -3:
+	    strcat(msg, _("Invalid order"));
+	    break;
+	default:
+	    /* COMMENTED OUT WHEN SUPPORT FOR 3rd ORDER WAS ADDED BY 'CRS'
+	       E12a = E12[0]; E12b = E12[1]; E12c = E12[2];
+	       N12a = N12[0]; N12b = N12[1]; N12c = N12[2];
+	       E21a = E21[0]; E21b = E21[1]; E21c = E21[2];
+	       N21a = N21[0]; N21b = N21[1]; N21c = N21[2];
+	     */
+	    return 1;
+	}
+    }
+    else {
+	switch (I_compute_georef_equations(&cp, E12, N12, E21, N21, order)) {
+	case 0:
+	    sprintf(&msg[strlen(msg)],
+		    _("Not enough active control points for current order, %d are required."),
+		    (order + 1) * (order + 2) / 2);
+	    break;
+	case -1:
+	    strcat(msg, _("Poorly placed control points."));
+	    strcat(msg, _(" Can not generate the transformation equation."));
+	    break;
+	case -2:
+	    strcat(msg, _("Not enough memory to solve for transformation equation"));
+	    break;
+	case -3:
+	    strcat(msg, _("Invalid order"));
+	    break;
+	default:
+	    /* COMMENTED OUT WHEN SUPPORT FOR 3rd ORDER WAS ADDED BY 'CRS'
+	       E12a = E12[0]; E12b = E12[1]; E12c = E12[2];
+	       N12a = N12[0]; N12b = N12[1]; N12c = N12[2];
+	       E21a = E21[0]; E21b = E21[1]; E21c = E21[2];
+	       N21a = N21[0]; N21b = N21[1]; N21c = N21[2];
+	     */
+	    return 1;
+	}
     }
     G_fatal_error(msg);
 

imagery/i.rectify/get_wind.c

@@ -10,13 +10,19 @@ int georef_window(struct Cell_head *w1, struct Cell_head *w2, int order, double
     } nw, ne, se, sw;
 
     /* extends */
-    I_georef(w1->west, w1->north, &e, &n, E12, N12, order);
+    if (order == 0)
+	I_georef_tps(w1->west, w1->north, &e, &n, E12_t, N12_t, &cp, 1);
+    else
+	I_georef(w1->west, w1->north, &e, &n, E12, N12, order);
     w2->north = w2->south = n;
     w2->west = w2->east = e;
     nw.n = n;
     nw.e = e;
 
-    I_georef(w1->east, w1->north, &e, &n, E12, N12, order);
+    if (order == 0)
+	I_georef_tps(w1->east, w1->north, &e, &n, E12_t, N12_t, &cp, 1);
+    else
+	I_georef(w1->east, w1->north, &e, &n, E12, N12, order);
     ne.n = n;
     ne.e = e;
     if (n > w2->north)
@@ -28,7 +34,10 @@ int georef_window(struct Cell_head *w1, struct Cell_head *w2, int order, double
     if (e < w2->west)
 	w2->west = e;
 
-    I_georef(w1->west, w1->south, &e, &n, E12, N12, order);
+    if (order == 0)
+	I_georef_tps(w1->west, w1->south, &e, &n, E12_t, N12_t, &cp, 1);
+    else
+	I_georef(w1->west, w1->south, &e, &n, E12, N12, order);
     sw.n = n;
     sw.e = e;
     if (n > w2->north)
@@ -40,7 +49,10 @@ int georef_window(struct Cell_head *w1, struct Cell_head *w2, int order, double
     if (e < w2->west)
 	w2->west = e;
 
-    I_georef(w1->east, w1->south, &e, &n, E12, N12, order);
+    if (order == 0)
+	I_georef_tps(w1->east, w1->south, &e, &n, E12_t, N12_t, &cp, 1);
+    else
+	I_georef(w1->east, w1->south, &e, &n, E12, N12, order);
     se.n = n;
     se.e = e;
     if (n > w2->north)

imagery/i.rectify/global.h

@@ -48,6 +48,11 @@ struct menu
 extern double E12[10], N12[10];
 extern double E21[10], N21[10];
 
+extern double *E12_t, *N12_t;
+extern double *E21_t, *N21_t;
+
+extern struct Control_Points cp;
+
 /* DELETED WITH CRS MODIFICATIONS
    extern double E12a, E12b, E12c, N12a, N12b, N12c;
    extern double E21a, E21b, E21c, N21a, N21b, N21c;

imagery/i.rectify/main.c

@@ -44,6 +44,11 @@ func interpolate;
 double E12[10], N12[10];
 double E21[10], N21[10];
 
+double *E12_t, *N12_t;
+double *E21_t, *N21_t;
+
+struct Control_Points cp;
+
 /* DELETED WITH CRS MODIFICATIONS
    double E12a, E12b, E12c, N12a, N12b, N12c;
    double E21a, E21b, E21c, N21a, N21b, N21c;
@@ -85,7 +90,7 @@ int main(int argc, char *argv[])
      *mem,			/* amount of memory for cache */
      *interpol;			/* interpolation method:
 				   nearest neighbor, bilinear, cubic */
-    struct Flag *c, *a;
+    struct Flag *c, *a, *t;
     struct GModule *module;
 
     G_gisinit(argv[0]);
@@ -151,6 +156,10 @@ int main(int argc, char *argv[])
     a->key = 'a';
     a->description = _("Rectify all raster maps in group");
 
+    t = G_define_flag();
+    t->key = 't';
+    t->description = _("Use thin plate spline");
+
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
 
@@ -186,9 +195,11 @@ int main(int argc, char *argv[])
 	k++;
     }
 
-    if (order < 1 || order > 3)  /* MAXORDER in lib/imagery/georef.c */
+    if (!t->answer && (order < 1 || order > 3))  /* MAXORDER in lib/imagery/georef.c */
 	G_fatal_error(_("Invalid order (%d); please enter 1 to %d"), order,
 		      3);
+    if (t->answer)
+	order = 0;
 
     /* determine the number of files in this group */
     if (I_get_group_ref(group, &ref) <= 0) {

imagery/i.rectify/rectify.c

@@ -72,7 +72,10 @@ int rectify(char *name, char *mapset, char *result, int order, char *interp_meth
 	    e1 = target_window.west + (col + 0.5) * target_window.ew_res;
 
 	    /* backwards transformation of target cell center */
-	    I_georef(e1, n1, &ex, &nx, E21, N21, order);
+	    if (order == 0)
+		I_georef_tps(e1, n1, &ex, &nx, E21_t, N21_t, &cp, 0);
+	    else
+		I_georef(e1, n1, &ex, &nx, E21, N21, order);
 
 	    /* convert to row/column indices of source raster */
 	    row_idx = (cellhd.north - nx) / cellhd.ns_res;