1
2 /***********************************************************************
3 *
4 * MODULE: r.resamp.bspline
5 *
6 * AUTHOR(S): Roberto Antolin (v.surf.bspline)
7 * Markus Metz (adapted for r.resamp.bspline)
8 *
9 * PURPOSE: Spline Interpolation and cross correlation
10 *
11 * COPYRIGHT: (C) 2006 by Politecnico di Milano -
12 * Polo Regionale di Como
13 *
14 * This program is free software under the
15 * GNU General Public License (>=v2).
16 * Read the file COPYING that comes with GRASS
17 * for details.
18 *
19 **************************************************************************/
20
21 /*INCLUDES*/
22 #include <grass/config.h>
23 #include <stdlib.h>
24 #include <string.h>
25 #include <math.h>
26 #include "bspline.h"
27
28 #define NDATA_MAX 100
29 #define PARAM_LAMBDA 6
30 #define PARAM_SPLINE 0
31 #define SWAP(a, b) {double t=a; a=b; b=t;}
32
33 /* internal functions */
calc_mean(double * values,int nvalues)34 double calc_mean(double *values, int nvalues)
35 {
36 int i;
37 double sum = .0;
38
39 if (nvalues == 0)
40 return .0;
41 for (i = 0; i < nvalues; i++)
42 sum += values[i];
43 return sum / nvalues;
44 }
45
46
calc_root_mean_square(double * values,int nvalues)47 double calc_root_mean_square(double *values, int nvalues)
48 {
49 int i;
50 double rms, sum = .0;
51
52 if (nvalues == 0)
53 return .0;
54
55 for (i = 0; i < nvalues; i++)
56 sum += pow(values[i], 2) / nvalues;
57
58 rms = sqrt(sum);
59 return rms;
60
61 }
62
calc_standard_deviation(double * values,int nvalues)63 double calc_standard_deviation(double *values, int nvalues)
64 {
65 double mean, rms, stdev;
66
67 if (nvalues == 0)
68 return .0;
69
70 rms = calc_root_mean_square(values, nvalues);
71 mean = calc_mean(values, nvalues);
72
73 stdev = sqrt(pow(rms, 2) - pow(mean, 2));
74 return stdev;
75 }
76
alloc_Stats(int n)77 struct Stats alloc_Stats(int n)
78 {
79 double *err, *stm;
80 struct Stats stat;
81
82 stat.n_points = n;
83 err = (double *)G_calloc(n, sizeof(double));
84 stm = (double *)G_calloc(n, sizeof(double));
85
86 stat.error = err;
87 stat.estima = stm;
88
89 return stat;
90 }
91
find_minimum(double * values,int * l_min)92 double find_minimum(double *values, int *l_min)
93 {
94 int l;
95 double min;
96
97 min = values[0];
98
99 for (l = 0; l < PARAM_LAMBDA; l++) {
100 if (min > values[l]) {
101 min = values[l];
102 *l_min = l;
103 }
104 }
105 return min;
106 }
107
swap(struct Point * point,int a,int b)108 struct Point *swap(struct Point *point, int a, int b)
109 {
110
111 SWAP(point[a].coordX, point[b].coordX); /* Once the last value is let out, it is swap with j-value */
112 SWAP(point[a].coordY, point[b].coordY);
113 SWAP(point[a].coordZ, point[b].coordZ);
114 SWAP(point[a].cat, point[b].cat);
115
116 return point;
117 }
118
119 /*-------------------------------------------------------------------------------------------*/
cross_correlation(SEGMENT * in_seg,struct Cell_head * src_reg,double passWE,double passNS)120 int cross_correlation(SEGMENT *in_seg, struct Cell_head *src_reg, double passWE, double passNS)
121 /*
122 matrix: matrix with raster values
123 passWE: spline step in West-East direction
124 passNS: spline step in North-South direction
125
126 RETURN:
127 TRUE on success
128 FALSE on failure
129 */
130 {
131 int bilin = TRUE; /*booleans */
132 int nsplx, nsply, nparam_spl, ndata;
133 double *mean, *rms, *stdev;
134
135 #ifdef nodef
136 double rms_min, stdev_min;
137 int lbd_min; /* lbd_min: index where minimun lambda index is found */
138 #endif
139
140 /* double lambda[PARAM_LAMBDA] = { 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0 }; */ /* Fixed values (by the moment) */
141 double lambda[PARAM_LAMBDA] = { 0.0001, 0.001, 0.005, 0.01, 0.02, 0.05 }; /* Fixed values (by the moment) */
142
143 double *TN, *Q, *parVect; /* Interpolation and least-square vectors */
144 double **N, **obsVect; /* Interpolation and least-square matrix */
145
146 struct Point *observ;
147 struct Stats stat_vect;
148
149 struct Cell_head region;
150
151 G_get_window(®ion);
152
153
154 G_debug(5,
155 "CrossCorrelation: Some tests using different lambda_i values will be done");
156
157 ndata = region.rows * region.cols;
158
159 if (ndata > NDATA_MAX)
160 G_warning(_("%d are too many cells, recommended are < 100 cells. "
161 "The cross validation would take too much time."), ndata);
162
163 /*points = Vect_new_line_struct (); */
164 /*Cats = Vect_new_cats_struct (); */
165
166 /* Current region is read and points recorded into observ */
167 observ = P_Read_Raster_Region_Map(in_seg, ®ion, src_reg, &ndata, 1024);
168 G_debug(5, "CrossCorrelation: %d points read in region. ", ndata);
169 G_verbose_message(n_("%d point read in region",
170 "%d points read in region", ndata),
171 ndata);
172
173 if (ndata > 50)
174 G_warning(_("Maybe it takes too long. "
175 "Consider reducing the region extents."));
176 else
177 G_debug(5, "CrossCorrelation: It shouldn't take too long.");
178
179 if (ndata > 0) { /* If at least one point is in the region */
180 int i, j, lbd; /* lbd: lambda index */
181 int BW;
182 double mean_reg, *obs_mean;
183 int verbosity;
184
185 mean = G_alloc_vector(PARAM_LAMBDA); /* Alloc as much mean, rms and stdev values as the total */
186 rms = G_alloc_vector(PARAM_LAMBDA); /* number of parameter used used for cross validation */
187 stdev = G_alloc_vector(PARAM_LAMBDA);
188
189 verbosity = G_verbose(); /* store for later reset */
190
191 /* Setting number of splines as a function of WE and SN spline steps */
192 nsplx = ceil((region.east - region.west) / passWE);
193 nsply = ceil((region.north - region.south) / passNS);
194 nparam_spl = nsplx * nsply; /* Total number of splines */
195
196 if (nparam_spl > 22900)
197 G_fatal_error(_("Too many splines (%d x %d). "
198 "Consider changing spline steps \"ew_step=\" \"ns_step=\"."),
199 nsplx, nsply);
200
201 BW = P_get_BandWidth(bilin, nsply);
202 /**/
203 /*Least Squares system */
204 N = G_alloc_matrix(nparam_spl, BW); /* Normal matrix */
205 TN = G_alloc_vector(nparam_spl); /* vector */
206 parVect = G_alloc_vector(nparam_spl); /* Parameters vector */
207 obsVect = G_alloc_matrix(ndata, 3); /* Observation vector */
208 Q = G_alloc_vector(ndata); /* "a priori" var-cov matrix */
209
210 obs_mean = G_alloc_vector(ndata);
211 stat_vect = alloc_Stats(ndata);
212
213 for (lbd = 0; lbd < PARAM_LAMBDA; lbd++) { /* For each lambda value */
214
215 G_message(_("Beginning cross validation with "
216 "lambda_i=%.4f ... (%d of %d)"), lambda[lbd],
217 lbd+1, PARAM_LAMBDA);
218
219 /*
220 How the cross correlation algorithm is done:
221 For each cycle, only the first ndata-1 "observ" elements are considered for the
222 interpolation. Within every interpolation mean is calculated to lowering border
223 errors. The point left out will be used for an estimation. The error between the
224 estimation and the observation is recorded for further statistics.
225 At the end of the cycle, the last point, that is, the ndata-1 index, and the point
226 with j index are swapped.
227 */
228 for (j = 0; j < ndata; j++) { /* Cross Correlation will use all ndata points */
229 double out_x, out_y, out_z; /* This point is left out */
230
231 for (i = 0; i < ndata; i++) { /* Each time, only the first ndata-1 points */
232
233 /* Setting obsVect vector & Q matrix */
234 Q[i] = 1; /* Q=I */
235 obsVect[i][0] = observ[i].coordX;
236 obsVect[i][1] = observ[i].coordY;
237
238 obsVect[i][2] = observ[i].coordZ;
239 obs_mean[i] = observ[i].coordZ;
240 } /* i index */
241
242 /* Mean calculation for every point less the last one */
243 mean_reg = calc_mean(obs_mean, ndata - 1);
244
245 for (i = 0; i < ndata; i++)
246 obsVect[i][2] -= mean_reg;
247
248 /* This is left out */
249 out_x = observ[ndata - 1].coordX;
250 out_y = observ[ndata - 1].coordY;
251 out_z = obsVect[ndata - 1][2];
252
253 if (bilin) { /* Bilinear interpolation */
254 normalDefBilin(N, TN, Q, obsVect, passWE, passNS, nsplx,
255 nsply, region.west, region.south,
256 ndata - 1, nparam_spl, BW);
257 nCorrectGrad(N, lambda[lbd], nsplx, nsply, passWE,
258 passNS);
259 }
260 else { /* Bicubic interpolation */
261 normalDefBicubic(N, TN, Q, obsVect, passWE, passNS, nsplx,
262 nsply, region.west, region.south,
263 ndata - 1, nparam_spl, BW);
264 nCorrectGrad(N, lambda[lbd], nsplx, nsply, passWE,
265 passNS);
266 }
267
268 /*
269 if (bilin) interpolation (&interp, P_BILINEAR);
270 else interpolation (&interp, P_BICUBIC);
271 */
272 G_set_verbose(G_verbose_min());
273 G_math_solver_cholesky_sband(N, parVect, TN, nparam_spl, BW);
274 G_set_verbose(verbosity);
275
276 /* Estimation of j-point */
277 if (bilin)
278 stat_vect.estima[j] =
279 dataInterpolateBilin(out_x, out_y, passWE, passNS,
280 nsplx, nsply, region.west,
281 region.south, parVect);
282
283 else
284 stat_vect.estima[j] =
285 dataInterpolateBilin(out_x, out_y, passWE, passNS,
286 nsplx, nsply, region.west,
287 region.south, parVect);
288
289 /* Difference between estimated and observated i-point */
290 stat_vect.error[j] = out_z - stat_vect.estima[j];
291 G_debug(1, "CrossCorrelation: stat_vect.error[%d] = %lf", j,
292 stat_vect.error[j]);
293
294 /* Once the last value is left out, it is swapped with j-value */
295 observ = swap(observ, j, ndata - 1);
296
297 G_percent(j, ndata, 2);
298 }
299
300 mean[lbd] = calc_mean(stat_vect.error, stat_vect.n_points);
301 rms[lbd] =
302 calc_root_mean_square(stat_vect.error, stat_vect.n_points);
303 stdev[lbd] =
304 calc_standard_deviation(stat_vect.error, stat_vect.n_points);
305
306 G_message(_("Mean = %.5lf"), mean[lbd]);
307 G_message(_("Root Mean Square (RMS) = %.5lf"),
308 rms[lbd]);
309 G_message("---");
310 } /* ENDFOR each lambda value */
311
312 G_free_matrix(N);
313 G_free_vector(TN);
314 G_free_vector(Q);
315 G_free_matrix(obsVect);
316 G_free_vector(parVect);
317 #ifdef nodef
318 /*TODO: if the minimum lambda is wanted, the function declaration must be changed */
319 /* At this moment, consider rms only */
320 rms_min = find_minimum(rms, &lbd_min);
321 stdev_min = find_minimum(stdev, &lbd_min);
322
323 /* Writing some output */
324 G_message(_("Different number of splines and lambda_i values have "
325 "been taken for the cross correlation"));
326 G_message(_("The minimum value for the test (rms=%lf) was "
327 "obtained with: lambda_i = %.3f"),
328 rms_min,
329 lambda[lbd_min]);
330
331 *lambda_min = lambda[lbd_min];
332 #endif
333
334 G_message(_("Table of results:"));
335 fprintf(stdout, _(" lambda | mean | rms |\n"));
336 for (lbd = 0; lbd < PARAM_LAMBDA; lbd++) {
337 fprintf(stdout, " %9.5f | %10.4f | %10.4f |\n", lambda[lbd],
338 mean[lbd], rms[lbd]);
339 }
340
341 G_free_vector(mean);
342 G_free_vector(rms);
343 } /* ENDIF (ndata > 0) */
344 else
345 G_warning(_("No point lies into the current region"));
346
347 G_free(observ);
348 return TRUE;
349 }
350
351 #ifdef nodef
interpolation(struct ParamInterp * interp,boolean bilin)352 void interpolation(struct ParamInterp *interp, boolean bilin)
353 {
354 if (bilin == P_BILINEAR) { /* Bilinear interpolation */
355 normalDefBilin(interp->N, interp->TN, interp->Q, interp->obsVect,
356 interp->stepE, interp->stepN, interp->nsplx,
357 interp->nsply, interp->region.west,
358 interp->region.south, interp->ndata,
359 interp->nparam_spl, interp->BW);
360
361 nCorrectGrad(interp->N, interp->lambda[lbd], interp->nsplx,
362 interp->nsply, interp->stepE, interp->stepN);
363 }
364 else { /* Bicubic interpolation */
365 normalDefBicubic(interp->N, interp->TN, interp->Q, interp->obsVect,
366 interp->stepE, interp->stepN, interp->nsplx,
367 interp->nsply, interp->region.west,
368 interp->region.south, interp->ndata,
369 interp->nparam_spl, interp->BW);
370
371 nCorrectGrad(interp->N, interp->lambda[lbd], interp->nsplx,
372 interp->nsply, interp->stepE, interp->stepN);
373 }
374 return TRUE;
375 }
376 #endif
377