GRASS Programmer's Manual  6.5.svn(2014)-r66266
transform.c
Go to the documentation of this file.
1
20 /****************************************************************
21 note: uses sqrt() from math library
22 *****************************************************************
23 Points from one system may be converted into the second by
24 use of one of the two equation routines.
25
26 transform_a_into_b (ax,ay,bx,by)
27
28  double ax,ay; input point from system a
29  double *bx,*by; resultant point in system b
30
31 transform_b_into_a (bx,by,ax,ay)
32
33  double bx,by; input point from system b
34  double *ax,*ay; resultant point in system a
35 *****************************************************************
36 Residual analysis on the equation can be run to test how well
37 the equations work. Either test how well b is predicted by a
38 or vice versa.
39
40 residuals_a_predicts_b (ax,ay,bx,by,use,n,residuals,rms)
41 residuals_b_predicts_a (ax,ay,bx,by,use,n,residuals,rms)
42
43  double ax[], ay[]; coordinate from system a
44  double bx[], by[]; coordinate from system b
45  char use[]; use point flags
46  int n; number of points in ax,ay,bx,by
47  double residual[] residual error for each point
48  double *rms; overall root mean square error
49 ****************************************************************/
50
51 #include <stdio.h>
52 #include <math.h>
53 #include <grass/libtrans.h>
54
55 /* the coefficients */
56 static double A0, A1, A2, A3, A4, A5;
57 static double B0, B1, B2, B3, B4, B5;
58
59 /* function prototypes */
60 static int resid(double *, double *, double *, double *, int *, int, double *,
61  double *, int);
62
63
90 int compute_transformation_coef(double ax[], double ay[], double bx[],
91  double by[], int *use, int n)
92 {
93  int i;
94  int j;
95  int count;
96  double aa[3];
97  double aar[3];
98  double bb[3];
99  double bbr[3];
100
101  double cc[3][3];
102  double x;
103
104  count = 0;
105  for (i = 0; i < n; i++)
106  if (use[i])
107  count++;
108  if (count < 4)
109  return -2; /* must have at least 4 points */
110
111  for (i = 0; i < 3; i++) {
112  aa[i] = bb[i] = 0.0;
113
114  for (j = 0; j < 3; j++)
115  cc[i][j] = 0.0;
116  }
117
118  for (i = 0; i < n; i++) {
119  if (!use[i])
120  continue; /* skip this point */
121  cc[0][0] += 1;
122  cc[0][1] += bx[i];
123  cc[0][2] += by[i];
124
125  cc[1][1] += bx[i] * bx[i];
126  cc[1][2] += bx[i] * by[i];
127  cc[2][2] += by[i] * by[i];
128
129  aa[0] += ay[i];
130  aa[1] += ay[i] * bx[i];
131  aa[2] += ay[i] * by[i];
132
133  bb[0] += ax[i];
134  bb[1] += ax[i] * bx[i];
135  bb[2] += ax[i] * by[i];
136  }
137
138  cc[1][0] = cc[0][1];
139  cc[2][0] = cc[0][2];
140  cc[2][1] = cc[1][2];
141
142  /* aa and bb are solved */
143  if (inverse(cc) < 0)
144  return (-1);
145  if (m_mult(cc, aa, aar) < 0 || m_mult(cc, bb, bbr) < 0)
146  return (-1);
147
148  /* the equation coefficients */
149  B0 = aar[0];
150  B1 = aar[1];
151  B2 = aar[2];
152
153  B3 = bbr[0];
154  B4 = bbr[1];
155  B5 = bbr[2];
156
157  /* the inverse equation */
158  x = B2 * B4 - B1 * B5;
159
160  if (!x)
161  return (-1);
162
163  A0 = (B1 * B3 - B0 * B4) / x;
164  A1 = -B1 / x;
165  A2 = B4 / x;
166  A3 = (B0 * B5 - B2 * B3) / x;
167  A4 = B2 / x;
168  A5 = -B5 / x;
169
170  return 1;
171 }
172
173
174 int transform_a_into_b(double ax, double ay, double *bx, double *by)
175 {
176  *by = A0 + A1 * ax + A2 * ay;
177  *bx = A3 + A4 * ax + A5 * ay;
178
179  return 0;
180 }
181
182
183 int transform_b_into_a(double bx, double by, double *ax, double *ay)
184 {
185  *ay = B0 + B1 * bx + B2 * by;
186  *ax = B3 + B4 * bx + B5 * by;
187
188  return 0;
189 }
190
191 /**************************************************************
192 These routines are internal to this source code
193
194 solve (a, b)
195  double a[3][3]
196  double b[3]
197
198  equation solver used by compute_transformation_coef()
199 **************************************************************/
200
201 /* #define abs(xx) (xx >= 0 ? xx : -xx) */
202 /* #define N 3 */
203
204
205 int residuals_a_predicts_b(double ax[], double ay[], double bx[], double by[],
206  int use[], int n, double residuals[], double *rms)
207 {
208  resid(ax, ay, bx, by, use, n, residuals, rms, 1);
209
210  return 0;
211 }
212
213
214 int residuals_b_predicts_a(double ax[], double ay[], double bx[], double by[],
215  int use[], int n, double residuals[], double *rms)
216 {
217  resid(ax, ay, bx, by, use, n, residuals, rms, 0);
218
219  return 0;
220 }
221
222
232 {
233  fprintf(stdout, "\nTransformation Matrix\n");
234  fprintf(stdout, "| xoff a b |\n");
235  fprintf(stdout, "| yoff d e |\n");
236  fprintf(stdout, "-------------------------------------------\n");
237  fprintf(stdout, "%f %f %f \n", -B3, B2, -B5);
238  fprintf(stdout, "%f %f %f \n", -B0, -B1, B4);
239  fprintf(stdout, "-------------------------------------------\n");
240
241  return 1;
242 }
243
244
245 static int resid(double ax[], double ay[], double bx[], double by[],
246  int use[], int n, double residuals[], double *rms, int atob)
247 {
248  double x, y;
249  int i;
250  int count;
251  double sum;
252  double delta;
253  double dx, dy;
254
255  count = 0;
256  sum = 0.0;
257  for (i = 0; i < n; i++) {
258  if (!use[i])
259  continue;
260
261  count++;
262  if (atob) {
263  transform_a_into_b(ax[i], ay[i], &x, &y);
264  dx = x - bx[i];
265  dy = y - by[i];
266  }
267  else {
268  transform_b_into_a(bx[i], by[i], &x, &y);
269  dx = x - ax[i];
270  dy = y - ay[i];
271  }
272
273  delta = dx * dx + dy * dy;
274  residuals[i] = sqrt(delta);
275  sum += delta;
276  }
277  *rms = sqrt(sum / count);
278
279  return 0;
280 }
int transform_a_into_b(double ax, double ay, double *bx, double *by)
Definition: transform.c:174
int inverse(double m[N][N])
Definition: inverse.c:17
int transform_b_into_a(double bx, double by, double *ax, double *ay)
Definition: transform.c:183
int count
int y
Definition: plot.c:34
int compute_transformation_coef(double ax[], double ay[], double bx[], double by[], int *use, int n)
Definition: transform.c:90
int residuals_b_predicts_a(double ax[], double ay[], double bx[], double by[], int use[], int n, double residuals[], double *rms)
Definition: transform.c:214
int m_mult(double a[N][N], double b[N], double c[N])
Definition: m_mult.c:12
int residuals_a_predicts_b(double ax[], double ay[], double bx[], double by[], int use[], int n, double residuals[], double *rms)
Definition: transform.c:205
int n