← Back to branch summary

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

« back to all changes in this revision

Viewing changes to imagery/i.smap/model.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:
imagery/i.smap/model.c
 
1
#include <grass/gis.h>
 
2
#include <grass/raster.h>
 
3
#include <grass/glocale.h>
 
4
#include <grass/imagery.h>
 
5
#include <grass/gmath.h>
 
6
#include "bouman.h"
 
7
#include "region.h"
 
8
 
 
9
#define PI M_PI
 
10
 
 
11
 
 
12
void extract_init(struct SigSet *S)
 
13
{
 
14
    int m;
 
15
    int i;
 
16
    int b1, b2;
 
17
    int nbands;
 
18
    double *lambda;
 
19
    double **tmp_mat;
 
20
    struct ClassSig *C;
 
21
    struct SubSig *SubS;
 
22
 
 
23
    nbands = S->nbands;
 
24
    /* allocate scratch memory */
 
25
    lambda = G_alloc_vector(nbands);
 
26
    tmp_mat = G_alloc_matrix(nbands, nbands);
 
27
 
 
28
 
 
29
    /* invert matrix and compute constant for each subclass */
 
30
 
 
31
    /* for each class */
 
32
    for (m = 0; m < S->nclasses; m++) {
 
33
        C = &(S->ClassSig[m]);
 
34
 
 
35
        /* for each subclass */
 
36
        for (i = 0; i < C->nsubclasses; i++) {
 
37
            SubS = &(C->SubSig[i]);
 
38
 
 
39
            /* Test for symetric  matrix */
 
40
            for (b1 = 0; b1 < nbands; b1++)
 
41
                for (b2 = 0; b2 < nbands; b2++) {
 
42
                    if (SubS->R[b1][b2] != SubS->R[b2][b1])
 
43
                        G_warning(_("Nonsymetric covariance for class %d subclass %d"),
 
44
                                  m + 1, i + 1);
 
45
 
 
46
                    SubS->Rinv[b1][b2] = SubS->R[b1][b2];
 
47
                    tmp_mat[b1][b2] = SubS->R[b1][b2];
 
48
 
 
49
                }
 
50
 
 
51
            /* Test for positive definite matrix */
 
52
            G_math_eigen(tmp_mat, lambda, nbands);
 
53
            for (b1 = 0; b1 < nbands; b1++) {
 
54
                if (lambda[b1] <= 0.0)
 
55
                    G_warning(_("Nonpositive eigenvalues for class %d subclass %d"),
 
56
                              m + 1, i + 1);
 
57
            }
 
58
 
 
59
            /* Precomputes the cnst */
 
60
            SubS->cnst = (-nbands / 2.0) * log(2 * PI);
 
61
            for (b1 = 0; b1 < nbands; b1++) {
 
62
                SubS->cnst += -0.5 * log(lambda[b1]);
 
63
            }
 
64
 
 
65
            /* Precomputes the inverse of tex->R */
 
66
            invert(SubS->Rinv, nbands);
 
67
        }
 
68
    }
 
69
    G_free_vector(lambda);
 
70
    G_free_matrix(tmp_mat);
 
71
}
 
72
 
 
73
 
 
74
void extract(DCELL *** img,     /* multispectral image, img[band][i][j] */
 
75
             struct Region *region,     /* region to extract */
 
76
             LIKELIHOOD *** ll, /* log likelihood, ll[i][j][class] */
 
77
             struct SigSet *S   /* class signatures */
 
78
    )
 
79
{
 
80
    int i, j;                   /* column and row indexes */
 
81
    int m;                      /* class index */
 
82
    int k;                      /* subclass index */
 
83
    int b1, b2;                 /* spectral index */
 
84
    int no_data;                /* no data flag */
 
85
    int max_nsubclasses;        /* maximum number of subclasses */
 
86
    int nbands;                 /* number of spectral bands */
 
87
    double *subll;              /* log likelihood of subclasses */
 
88
    double *diff;
 
89
    double maxlike = 0.0L;
 
90
    double subsum;
 
91
    struct ClassSig *C;
 
92
    struct SubSig *SubS;
 
93
 
 
94
    nbands = S->nbands;
 
95
 
 
96
    /* determine the maximum number of subclasses */
 
97
    max_nsubclasses = 0;
 
98
    for (m = 0; m < S->nclasses; m++)
 
99
        if (S->ClassSig[m].nsubclasses > max_nsubclasses)
 
100
            max_nsubclasses = S->ClassSig[m].nsubclasses;
 
101
 
 
102
    /* allocate memory */
 
103
    diff = (double *)G_malloc(nbands * sizeof(double));
 
104
    subll = (double *)G_malloc(max_nsubclasses * sizeof(double));
 
105
 
 
106
    /* Compute log likelihood at each pixel and for every class. */
 
107
 
 
108
    /* for each pixel */
 
109
    for (i = region->ymin; i < region->ymax; i++)
 
110
        for (j = region->xmin; j < region->xmax; j++) {
 
111
 
 
112
            /* Check for no data condition */
 
113
            no_data = 1;
 
114
            for (b1 = 0; (b1 < nbands) && no_data; b1++)
 
115
                no_data = no_data && (Rast_is_d_null_value(&img[b1][i][j]));
 
116
 
 
117
            if (no_data) {
 
118
                for (m = 0; m < S->nclasses; m++)
 
119
                    ll[i][j][m] = 0.0;
 
120
            }
 
121
            else {
 
122
                /* for each class */
 
123
                for (m = 0; m < S->nclasses; m++) {
 
124
                    C = &(S->ClassSig[m]);
 
125
 
 
126
                    /* compute log likelihood for each subclass */
 
127
                    for (k = 0; k < C->nsubclasses; k++) {
 
128
                        SubS = &(C->SubSig[k]);
 
129
                        subll[k] = SubS->cnst;
 
130
                        for (b1 = 0; b1 < nbands; b1++) {
 
131
                            diff[b1] = img[b1][i][j] - SubS->means[b1];
 
132
                            subll[k] -=
 
133
                                0.5 * diff[b1] * diff[b1] *
 
134
                                SubS->Rinv[b1][b1];
 
135
                        }
 
136
                        for (b1 = 0; b1 < nbands; b1++)
 
137
                            for (b2 = b1 + 1; b2 < nbands; b2++)
 
138
                                subll[k] -=
 
139
                                    diff[b1] * diff[b2] * SubS->Rinv[b1][b2];
 
140
                    }
 
141
 
 
142
                    /* shortcut for one subclass */
 
143
                    if (C->nsubclasses == 1) {
 
144
                        ll[i][j][m] = subll[0];
 
145
                    }
 
146
                    /* compute mixture likelihood */
 
147
                    else {
 
148
                        /* find the most likely subclass */
 
149
                        for (k = 0; k < C->nsubclasses; k++) {
 
150
                            if (k == 0)
 
151
                                maxlike = subll[k];
 
152
                            if (subll[k] > maxlike)
 
153
                                maxlike = subll[k];
 
154
                        }
 
155
 
 
156
                        /* Sum weighted subclass likelihoods */
 
157
                        subsum = 0;
 
158
                        for (k = 0; k < C->nsubclasses; k++)
 
159
                            subsum +=
 
160
                                exp(subll[k] - maxlike) * C->SubSig[k].pi;
 
161
 
 
162
                        ll[i][j][m] = log(subsum) + maxlike;
 
163
                    }
 
164
                }
 
165
            }
 
166
        }
 
167
    G_free((char *)diff);
 
168
    G_free((char *)subll);
 
169
}