[GRASS-SVN] r51549 - grass/trunk/vector/v.rectify

[svn_grass at osgeo.org]

Author: mmetz
Date: 2012-04-26 10:44:32 -0700 (Thu, 26 Apr 2012)
New Revision: 51549

Modified:
   grass/trunk/vector/v.rectify/orthorot.c
Log:
code cleanup

Modified: grass/trunk/vector/v.rectify/orthorot.c
===================================================================
--- grass/trunk/vector/v.rectify/orthorot.c	2012-04-26 16:21:52 UTC (rev 51548)
+++ grass/trunk/vector/v.rectify/orthorot.c	2012-04-26 17:44:32 UTC (rev 51549)
@@ -41,7 +41,7 @@
 static int calcscale(struct Control_Points_3D *, double *);
 
 void transpose_matrix(int, int, double **, double **);
-void matmult(int, int, double **, double **, double **);
+void matmult(int, int, int, double **, double **, double **);
 void copy_matrix(int, int, double **, double **);
 void init_matrix(int, int, double **);
 void scale_matrix(int, int, double, double **, double **);
@@ -293,6 +293,7 @@
     double **src_mat_T = NULL;
     double **dest_mat = NULL;
     double **dest_mat_T = NULL;
+    double **src_dest_mat = NULL;
     double *S_vec = NULL;
     double **R_mat = NULL;
     double **R_mat_T = NULL;
@@ -301,10 +302,8 @@
     double **mat_nm1 = NULL;
     double **mat_nm2 = NULL;
     double **mat_nn1 = NULL;
-    double **mat_nn2 = NULL;
     double **E_mat = NULL;
     double **P_mat = NULL;
-    double **P_mat_T = NULL;
     double **Q_mat = NULL;
     double *D_vec = NULL;
     double *one_vec = NULL;
@@ -344,6 +343,7 @@
 
     src_mat_T = G_alloc_matrix(n, m);
     dest_mat_T = G_alloc_matrix(n, m);
+    src_dest_mat = G_alloc_matrix(n, n);
     R_mat = G_alloc_matrix(n, n);
     R_mat_T = G_alloc_matrix(n, n);
 
@@ -352,11 +352,9 @@
     mat_nm1 = G_alloc_matrix(n, m);
     mat_nm2 = G_alloc_matrix(n, m);
     mat_nn1 = G_alloc_matrix(n, n);
-    mat_nn2 = G_alloc_matrix(n, n);
 
     E_mat = G_alloc_matrix(m, m);
     P_mat = G_alloc_matrix(ndims, ndims);
-    P_mat_T = G_alloc_matrix(ndims, ndims);
     Q_mat = G_alloc_matrix(ndims, ndims);
 
     transpose_matrix(m, n, dest_mat, dest_mat_T);
@@ -372,35 +370,35 @@
 	}
     }
 
-    matmult(n, m, dest_mat_T, E_mat, mat_nm1);
-    matmult(n, m, mat_nm1, src_mat, mat_nm2);
-    copy_matrix(ndims, ndims, mat_nm2, P_mat);
-    copy_matrix(ndims, ndims, mat_nm2, mat_nn1);
+    matmult(n, m, m, dest_mat_T, E_mat, mat_nm1);
+    matmult(n, m, n, mat_nm1, src_mat, src_dest_mat);
+    copy_matrix(n, n, src_dest_mat, P_mat);
+    copy_matrix(n, n, src_dest_mat, mat_nn1);
 
-    status = G_math_svduv(D_vec, mat_nn1, P_mat, ndims, Q_mat, ndims);
+    status = G_math_svduv(D_vec, mat_nn1, P_mat, n, Q_mat, n);
 
     if (status == 0)
 	status = MSUCCESS;
 
-    transpose_matrix(n, n, P_mat, P_mat_T);
+    transpose_matrix(n, n, P_mat, mat_nn1);
 
     /* rotation matrix */
-    matmult(n, n, Q_mat, P_mat_T, R_mat_T);
+    matmult(n, n, n, Q_mat, mat_nn1, R_mat_T);
     transpose_matrix(n, n, R_mat_T, R_mat);
 
     /* scale */
-    matmult(n, n, mat_nm2, R_mat_T, mat_nn2);
-    trace1 = trace(n, n, mat_nn2);
+    matmult(n, n, n, src_dest_mat, R_mat_T, mat_nn1);
+    trace1 = trace(n, n, mat_nn1);
 
     transpose_matrix(m, n, src_mat, src_mat_T);
-    matmult(n, m, src_mat_T, E_mat, mat_nm1);
-    matmult(n, m, mat_nm1, src_mat, mat_nm2);
-    trace2 = trace(n, n, mat_nm2);
+    matmult(n, m, m, src_mat_T, E_mat, mat_nm1);
+    matmult(n, m, n, mat_nm1, src_mat, mat_nn1);
+    trace2 = trace(n, n, mat_nn1);
 
     OR[14] = trace1 / trace2;
 
     /* shifts */
-    matmult(m, n, src_mat, R_mat_T, mat_mn1);
+    matmult(m, n, n, src_mat, R_mat_T, mat_mn1);
     scale_matrix(m, n, OR[14], mat_mn1, mat_mn2);
     subtract_matrix(m, n, dest_mat, mat_mn2, mat_mn1);
     scale_matrix(m, n, 1.0 / m, mat_mn1, mat_mn2);
@@ -426,10 +424,10 @@
     G_free_matrix(src_mat_T);
     G_free_matrix(dest_mat);
     G_free_matrix(dest_mat_T);
+    G_free_matrix(src_dest_mat);
     G_free_vector(D_vec);
     G_free_matrix(E_mat);
     G_free_matrix(P_mat);
-    G_free_matrix(P_mat_T);
     G_free_matrix(Q_mat);
     G_free_matrix(R_mat);
     G_free_matrix(R_mat_T);
@@ -438,7 +436,6 @@
     G_free_matrix(mat_nm1);
     G_free_matrix(mat_nm2);
     G_free_matrix(mat_nn1);
-    G_free_matrix(mat_nn2);
     G_free_vector(S_vec);
     G_free_vector(one_vec);
 
@@ -518,13 +515,16 @@
     }
 }
 
-void matmult(int m, int n, double **mat1, double **mat2, double **mat3)
+void matmult(int m, int n, int p, double **mat1, double **mat2, double **mat3)
 {
     int i, j, k;
     double sum;
 
+    /* input   mat1: m x n */
+    /* input   mat2: n x p */
+    /* output  mat3: m x p */
     for (i = 0; i < m; i++) {
-	for (j = 0; j < n; j++) {
+	for (j = 0; j < p; j++) {
 	    sum = 0.0;
 	    for (k = 0; k < n; k++) {
 		sum += mat1[i][k] * mat2[k][j];