← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to vector/lidar/v.surf.bspline/crosscorr.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
vector/lidar/v.surf.bspline/crosscorr.c
1
 
 
2
 
/***********************************************************************
3
 
 *
4
 
 * MODULE:       v.surf.bspline
5
 
 *
6
 
 * AUTHOR(S):    Roberto Antolin
7
 
 *
8
 
 * PURPOSE:      Spline Interpolation and cross correlation
9
 
 *
10
 
 * COPYRIGHT:    (C) 2006 by Politecnico di Milano -
11
 
 *                           Polo Regionale di Como
12
 
 *
13
 
 *               This program is free software under the
14
 
 *               GNU General Public License (>=v2).
15
 
 *               Read the file COPYING that comes with GRASS
16
 
 *               for details.
17
 
 *
18
 
 **************************************************************************/
19
 
 
20
 
 /*INCLUDES*/
21
 
#include <stdlib.h>
22
 
#include <string.h>
23
 
#include <math.h>
24
 
#include <grass/gis.h>
25
 
#include <grass/Vect.h>
26
 
#include <grass/dbmi.h>
27
 
#include <grass/glocale.h>
28
 
#include <grass/config.h>
29
 
#include <grass/PolimiFunct.h>
30
 
#include "bspline.h"
31
 
#define NDATA_MAX 100
32
 
#define PARAM_LAMBDA 6
33
 
#define PARAM_SPLINE 0
34
 
#define SWAP(a, b) {double t=a; a=b; b=t;}
35
 
 
36
 
/*-------------------------------------------------------------------------------------------*/
37
 
int cross_correlation(struct Map_info *Map, double passWE, double passNS)
38
 
    /*
39
 
       Map: Vector map from which cross-crorrelation will take values
40
 
       passWE: spline step in West-East direction
41
 
       passNS: spline step in North-South direction
42
 
 
43
 
       RETURN:
44
 
       TRUE on success
45
 
       FALSE on failure
46
 
     */
47
 
{
48
 
    int bilin = TRUE;           /*booleans */
49
 
    int nsplx, nsply, nparam_spl, ndata;
50
 
    double *mean, *rms, *stdev, rms_min, stdev_min;
51
 
 
52
 
    /* double lambda[PARAM_LAMBDA] = { 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0 }; */        /* Fixed values (at the moment) */
53
 
    double lambda[PARAM_LAMBDA] = { 0.01, 0.05, 0.1, 0.2, 0.3, 0.4 };   /* Fixed values (at the moment) */
54
 
    /* a more exhaustive search:
55
 
       #define PARAM_LAMBDA 11
56
 
       double lambda[PARAM_LAMBDA] = { 0.0001, 0.0005, 0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1.0, 5.0, 10.0 }; */
57
 
 
58
 
    double *TN, *Q, *parVect;   /* Interpolation and least-square vectors */
59
 
    double **N, **obsVect;      /* Interpolation and least-square matrix */
60
 
 
61
 
    struct Point *observ;
62
 
    struct Stats stat_vect;
63
 
 
64
 
    /*struct line_pnts *points; */
65
 
    /*struct line_cats *Cats; */
66
 
    struct Cell_head region;
67
 
 
68
 
    G_get_window(&region);
69
 
 
70
 
    extern int bspline_field;
71
 
    extern char *bspline_column;
72
 
    dbCatValArray cvarr;
73
 
 
74
 
    G_debug(5,
75
 
            "CrossCorrelation: Some tests using different lambda_i values will be done");
76
 
 
77
 
    ndata = Vect_get_num_lines(Map);
78
 
 
79
 
    if (ndata > NDATA_MAX)
80
 
        G_warning(_("%d are too many points. "
81
 
                    "The cross validation would take too much time."), ndata);
82
 
 
83
 
    /*points = Vect_new_line_struct (); */
84
 
    /*Cats = Vect_new_cats_struct (); */
85
 
 
86
 
    /* Current region is read and points recorded into observ */
87
 
    observ = P_Read_Vector_Region_Map(Map, &region, &ndata, 1024, 1);
88
 
    G_debug(5, "CrossCorrelation: %d points read in region. ", ndata);
89
 
    G_verbose_message(_("%d points read in region"),
90
 
                      ndata);
91
 
 
92
 
    if (ndata > 50)
93
 
        G_warning(_("Maybe, it takes too long. "
94
 
                    "It will depend on how many points you are considering."));
95
 
    else
96
 
        G_debug(5, "CrossCorrelation: It shouldn't take too long.");
97
 
 
98
 
    if (ndata > 0) {            /* If at least one point is in the region */
99
 
        int i, j, lbd;          /* lbd: lambda index */
100
 
        int BW, lbd_min;        /* lbd_min: index where minimun is found */
101
 
        double mean_reg, *obs_mean;
102
 
 
103
 
        int nrec, ctype = 0, verbosity;
104
 
        struct field_info *Fi;
105
 
        dbDriver *driver_cats;
106
 
 
107
 
        mean = G_alloc_vector(PARAM_LAMBDA);    /* Alloc as much mean, rms and stdev values as the total */
108
 
        rms = G_alloc_vector(PARAM_LAMBDA);     /* number of parameter used used for cross validation */
109
 
        stdev = G_alloc_vector(PARAM_LAMBDA);
110
 
 
111
 
        verbosity = G_verbose(); /* store for later reset */
112
 
 
113
 
        /* Working with attributes */
114
 
        if (bspline_field > 0) {
115
 
            db_CatValArray_init(&cvarr);
116
 
 
117
 
            Fi = Vect_get_field(Map, bspline_field);
118
 
            if (Fi == NULL)
119
 
              G_fatal_error(_("Database connection not defined for layer %d"),
120
 
                            bspline_field);
121
 
 
122
 
            driver_cats =
123
 
                db_start_driver_open_database(Fi->driver, Fi->database);
124
 
            G_debug(1, _("CrossCorrelation: driver=%s db=%s"), Fi->driver,
125
 
                    Fi->database);
126
 
 
127
 
            if (driver_cats == NULL)
128
 
                G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
129
 
                              Fi->database, Fi->driver);
130
 
 
131
 
            nrec =
132
 
                db_select_CatValArray(driver_cats, Fi->table, Fi->key,
133
 
                                      bspline_column, NULL, &cvarr);
134
 
            G_debug(3, "nrec = %d", nrec);
135
 
 
136
 
            ctype = cvarr.ctype;
137
 
            if (ctype != DB_C_TYPE_INT && ctype != DB_C_TYPE_DOUBLE)
138
 
                G_fatal_error(_("Column type not supported"));
139
 
 
140
 
            if (nrec < 0)
141
 
                G_fatal_error(_("No records selected from table <%s> "),
142
 
                              Fi->table);
143
 
 
144
 
            G_debug(1, "%d records selected from table",
145
 
                    nrec);
146
 
 
147
 
            db_close_database_shutdown_driver(driver_cats);
148
 
        }
149
 
 
150
 
        /* Setting number of splines as a function of WE and SN spline steps */
151
 
        nsplx = ceil((region.east - region.west) / passWE);
152
 
        nsply = ceil((region.north - region.south) / passNS);
153
 
        nparam_spl = nsplx * nsply;     /* Total number of splines */
154
 
 
155
 
        if (nparam_spl > 22900)
156
 
            G_fatal_error(_("Too many splines (%d x %d). "
157
 
                            "Consider changing spline steps \"sie=\" \"sin=\"."),
158
 
                          nsplx, nsply);
159
 
 
160
 
        BW = P_get_BandWidth(bilin, nsply);
161
 
        /**/
162
 
            /*Least Squares system */
163
 
            N = G_alloc_matrix(nparam_spl, BW); /* Normal matrix */
164
 
        TN = G_alloc_vector(nparam_spl);        /* vector */
165
 
        parVect = G_alloc_vector(nparam_spl);   /* Parameters vector */
166
 
        obsVect = G_alloc_matrix(ndata, 3);     /* Observation vector */
167
 
        Q = G_alloc_vector(ndata);              /* "a priori" var-cov matrix */
168
 
 
169
 
        obs_mean = G_alloc_vector(ndata);
170
 
        stat_vect = alloc_Stats(ndata);
171
 
 
172
 
        for (lbd = 0; lbd < PARAM_LAMBDA; lbd++) {      /* For each lambda value */
173
 
 
174
 
            G_message(_("Beginning cross validation with "
175
 
                        "lambda_i=%.4f ... (%d of %d)"), lambda[lbd],
176
 
                        lbd+1, PARAM_LAMBDA);
177
 
 
178
 
            /*
179
 
               How the cross correlation algorithm is done:
180
 
               For each cycle, only the first ndata-1 "observ" elements are considered for the 
181
 
               interpolation. Within every interpolation mean is calculated to lowering border 
182
 
               errors. The point left out will be used for an estimation. The error between the 
183
 
               estimation and the observation is recorded for further statistics.
184
 
               At the end of the cycle, the last point, that is, the ndata-1 index, and the point 
185
 
               with j index are swapped.
186
 
             */
187
 
            for (j = 0; j < ndata; j++) {       /* Cross Correlation will use all ndata points */
188
 
                double out_x, out_y, out_z;     /* This point is left out */
189
 
 
190
 
                for (i = 0; i < ndata; i++) {   /* Each time, only the first ndata-1 points */
191
 
                    double dval;                /* are considered in the interpolation */
192
 
 
193
 
                    /* Setting obsVect vector & Q matrix */
194
 
                    Q[i] = 1;   /* Q=I */
195
 
                    obsVect[i][0] = observ[i].coordX;
196
 
                    obsVect[i][1] = observ[i].coordY;
197
 
 
198
 
                    if (bspline_field > 0) {
199
 
                        int cat, ival, ret;
200
 
 
201
 
                        /*type = Vect_read_line (Map, points, Cats, observ[i].lineID); */
202
 
                        /*if ( !(type & GV_POINTS ) ) continue; */
203
 
 
204
 
                        /*Vect_cat_get ( Cats, bspline_field, &cat ); */
205
 
                        cat = observ[i].cat;
206
 
 
207
 
                        if (cat < 0)
208
 
                            continue;
209
 
 
210
 
                        if (ctype == DB_C_TYPE_INT) {
211
 
                            ret =
212
 
                                db_CatValArray_get_value_int(&cvarr, cat,
213
 
                                                             &ival);
214
 
                            obsVect[i][2] = ival;
215
 
                            obs_mean[i] = ival;
216
 
                        }
217
 
                        else {  /* DB_C_TYPE_DOUBLE */
218
 
                            ret =
219
 
                                db_CatValArray_get_value_double(&cvarr, cat,
220
 
                                                                &dval);
221
 
                            obsVect[i][2] = dval;
222
 
                            obs_mean[i] = dval;
223
 
                        }
224
 
                        if (ret != DB_OK) {
225
 
                            G_warning(_("No record for point (cat = %d)"),
226
 
                                      cat);
227
 
                            continue;
228
 
                        }
229
 
                    }
230
 
                    else {
231
 
                        obsVect[i][2] = observ[i].coordZ;
232
 
                        obs_mean[i] = observ[i].coordZ;
233
 
                    }
234
 
                }               /* i index */
235
 
 
236
 
                /* Mean calculation for every point less the last one */
237
 
                mean_reg = calc_mean(obs_mean, ndata - 1);
238
 
 
239
 
                for (i = 0; i < ndata; i++)
240
 
                    obsVect[i][2] -= mean_reg;
241
 
 
242
 
                /* This is left out */
243
 
                out_x = observ[ndata - 1].coordX;
244
 
                out_y = observ[ndata - 1].coordY;
245
 
                out_z = obsVect[ndata - 1][2];
246
 
 
247
 
                if (bilin) {    /* Bilinear interpolation */
248
 
                    normalDefBilin(N, TN, Q, obsVect, passWE, passNS, nsplx,
249
 
                                   nsply, region.west, region.south,
250
 
                                   ndata - 1, nparam_spl, BW);
251
 
                    nCorrectGrad(N, lambda[lbd], nsplx, nsply, passWE,
252
 
                                 passNS);
253
 
                }
254
 
                else {          /* Bicubic interpolation */
255
 
                    normalDefBicubic(N, TN, Q, obsVect, passWE, passNS, nsplx,
256
 
                                     nsply, region.west, region.south,
257
 
                                     ndata - 1, nparam_spl, BW);
258
 
                    nCorrectGrad(N, lambda[lbd], nsplx, nsply, passWE,
259
 
                                 passNS);
260
 
                }
261
 
 
262
 
                /* 
263
 
                   if (bilin) interpolation (&interp, P_BILINEAR);
264
 
                   else interpolation (&interp, P_BICUBIC);
265
 
                 */
266
 
                G_set_verbose(G_verbose_min());
267
 
                tcholSolve(N, TN, parVect, nparam_spl, BW);
268
 
                G_set_verbose(verbosity);
269
 
 
270
 
                /* Estimation of j-point */
271
 
                if (bilin)
272
 
                    stat_vect.estima[j] =
273
 
                        dataInterpolateBilin(out_x, out_y, passWE, passNS,
274
 
                                             nsplx, nsply, region.west,
275
 
                                             region.south, parVect);
276
 
 
277
 
                else
278
 
                    stat_vect.estima[j] =
279
 
                        dataInterpolateBilin(out_x, out_y, passWE, passNS,
280
 
                                             nsplx, nsply, region.west,
281
 
                                             region.south, parVect);
282
 
 
283
 
                /* Difference between estimated and observated i-point */
284
 
                stat_vect.error[j] = out_z - stat_vect.estima[j];
285
 
                G_debug(1, "CrossCorrelation: stat_vect.error[%d]  =  %lf", j,
286
 
                        stat_vect.error[j]);
287
 
 
288
 
                /* Once the last value is left out, it is swaped with j-value */
289
 
                observ = swap(observ, j, ndata - 1);
290
 
 
291
 
                G_percent(j, ndata, 2);
292
 
            }
293
 
 
294
 
            mean[lbd] = calc_mean(stat_vect.error, stat_vect.n_points);
295
 
            rms[lbd] =
296
 
                calc_root_mean_square(stat_vect.error, stat_vect.n_points);
297
 
            stdev[lbd] =
298
 
                calc_standard_deviation(stat_vect.error, stat_vect.n_points);
299
 
 
300
 
            G_message(_("Mean = %.5lf"), mean[lbd]);
301
 
            G_message(_("Root Mean Square (RMS) = %.5lf"),
302
 
                      rms[lbd]);
303
 
            G_message("---");
304
 
        }                       /* ENDFOR each lambda value */
305
 
 
306
 
        G_free_matrix(N);
307
 
        G_free_vector(TN);
308
 
        G_free_vector(Q);
309
 
        G_free_matrix(obsVect);
310
 
        G_free_vector(parVect);
311
 
#ifdef nodef
312
 
        /*TODO: if the minimum lambda is wanted, the function declaration must be changed */
313
 
        /* At this moment, consider rms only */
314
 
        rms_min = find_minimum(rms, &lbd_min);
315
 
        stdev_min = find_minimum(stdev, &lbd_min);
316
 
 
317
 
        /* Writing some output */
318
 
        G_message(_("Different number of splines and lambda_i values have "
319
 
                    "been taken for the cross correlation"));
320
 
        G_message(_("The minimum value for the test (rms=%lf) was "
321
 
                    "obtained with: lambda_i = %.3f"),
322
 
                  rms_min,
323
 
                  lambda[lbd_min]);
324
 
 
325
 
        *lambda_min = lambda[lbd_min];
326
 
#endif
327
 
 
328
 
        G_message(_("Table of results:"));
329
 
        fprintf(stdout, _("    lambda |       mean |        rms |\n"));
330
 
        for (lbd = 0; lbd < PARAM_LAMBDA; lbd++) {
331
 
            fprintf(stdout, " %9.5f | %10.4f | %10.4f |\n", lambda[lbd],
332
 
                    mean[lbd], rms[lbd]);
333
 
        }
334
 
        
335
 
        G_free_vector(mean);
336
 
        G_free_vector(rms);
337
 
    }                           /* ENDIF (ndata > 0) */
338
 
    else
339
 
        G_warning(_("No point lies into the current region"));
340
 
 
341
 
    G_free(observ);
342
 
    return TRUE;
343
 
}
344
 
 
345
 
#ifdef nodef
346
 
void interpolation(struct ParamInterp *interp, boolean bilin)
347
 
{
348
 
    if (bilin == P_BILINEAR) {  /* Bilinear interpolation */
349
 
        normalDefBilin(interp->N, interp->TN, interp->Q, interp->obsVect,
350
 
                       interp->passoE, interp->passoN, interp->nsplx,
351
 
                       interp->nsply, interp->region.west,
352
 
                       interp->region.south, interp->ndata,
353
 
                       interp->nparam_spl, interp->BW);
354
 
 
355
 
        nCorrectGrad(interp->N, interp->lambda[lbd], interp->nsplx,
356
 
                     interp->nsply, interp->passoE, interp->passoN);
357
 
    }
358
 
    else {                      /* Bicubic interpolation */
359
 
        normalDefBicubic(interp->N, interp->TN, interp->Q, interp->obsVect,
360
 
                         interp->passoE, interp->passoN, interp->nsplx,
361
 
                         interp->nsply, interp->region.west,
362
 
                         interp->region.south, interp->ndata,
363
 
                         interp->nparam_spl, interp->BW);
364
 
 
365
 
        nCorrectGrad(interp->N, interp->lambda[lbd], interp->nsplx,
366
 
                     interp->nsply, interp->passoE, interp->passoN);
367
 
    }
368
 
    return TRUE;
369
 
}
370
 
#endif
371
 
 
372
 
double calc_mean(double *values, int nvalues)
373
 
{
374
 
    int i;
375
 
    double sum = .0;
376
 
 
377
 
    if (nvalues == 0)
378
 
        return .0;
379
 
    for (i = 0; i < nvalues; i++)
380
 
        sum += values[i];
381
 
    return sum / nvalues;
382
 
}
383
 
 
384
 
 
385
 
double calc_root_mean_square(double *values, int nvalues)
386
 
{
387
 
    int i;
388
 
    double rms, sum = .0;
389
 
 
390
 
    if (nvalues == 0)
391
 
        return .0;
392
 
 
393
 
    for (i = 0; i < nvalues; i++)
394
 
        sum += pow(values[i], 2) / nvalues;
395
 
 
396
 
    rms = sqrt(sum);
397
 
    return rms;
398
 
 
399
 
}
400
 
 
401
 
double calc_standard_deviation(double *values, int nvalues)
402
 
{
403
 
    double mean, rms, stdev;
404
 
 
405
 
    if (nvalues == 0)
406
 
        return .0;
407
 
 
408
 
    rms = calc_root_mean_square(values, nvalues);
409
 
    mean = calc_mean(values, nvalues);
410
 
 
411
 
    stdev = sqrt(pow(rms, 2) - pow(mean, 2));
412
 
    return stdev;
413
 
}
414
 
 
415
 
struct Stats alloc_Stats(int n)
416
 
{
417
 
    double *err, *stm;
418
 
    struct Stats stat;
419
 
 
420
 
    stat.n_points = n;
421
 
    err = (double *)G_calloc(n, sizeof(double));
422
 
    stm = (double *)G_calloc(n, sizeof(double));
423
 
 
424
 
    stat.error = err;
425
 
    stat.estima = stm;
426
 
 
427
 
    return stat;
428
 
}
429
 
 
430
 
double find_minimum(double *values, int *l_min)
431
 
{
432
 
    int l;
433
 
    double min;
434
 
 
435
 
    min = values[0];
436
 
 
437
 
    for (l = 0; l < PARAM_LAMBDA; l++) {
438
 
        if (min > values[l]) {
439
 
            min = values[l];
440
 
            *l_min = l;
441
 
        }
442
 
    }
443
 
    return min;
444
 
}
445
 
 
446
 
struct Point *swap(struct Point *point, int a, int b)
447
 
{
448
 
    /* Once the last value is left out, it is swaped with j-value */
449
 
    SWAP(point[a].coordX, point[b].coordX);
450
 
    SWAP(point[a].coordY, point[b].coordY);
451
 
    SWAP(point[a].coordZ, point[b].coordZ);
452
 
    SWAP(point[a].cat, point[b].cat);
453
 
 
454
 
    return point;
455
 
}