← Back to branch summary

~ubuntu-branches/ubuntu/wily/grass/wily

« back to all changes in this revision

Viewing changes to imagery/i.smap/decimate.c

Tags: 7.0.0~rc1+ds1-1~exp1
http://bugs.debian.org/776349
* New upstream release candidate.
* Repack upstream tarball, remove precompiled Python objects.
* Add upstream metadata.
* Update gbp.conf and Vcs-Git URL to use the experimental branch.
* Update watch file for GRASS 7.0.
* Drop build dependencies for Tcl/Tk, add build dependencies:
  python-numpy, libnetcdf-dev, netcdf-bin, libblas-dev, liblapack-dev
* Update Vcs-Browser URL to use cgit instead of gitweb.
* Update paths to use grass70.
* Add configure options: --with-netcdf, --with-blas, --with-lapack,
  remove --with-tcltk-includes.
* Update patches for GRASS 7.
* Update copyright file, changes:
  - Update copyright years
  - Group files by license
  - Remove unused license sections
* Add patches for various typos.
* Fix desktop file with patch instead of d/rules.
* Use minimal dh rules.
* Bump Standards-Version to 3.9.6, no changes.
* Use dpkg-maintscript-helper to replace directories with symlinks.
  (closes: #776349)
* Update my email to use @debian.org address.

Show diffs side-by-side

added added

removed removed

Lines of Context:
imagery/i.smap/decimate.c
 
1
#include <grass/gis.h>
 
2
#include "bouman.h"
 
3
#include "region.h"
 
4
 
 
5
 
 
6
static void decimate(LIKELIHOOD ***, struct Region *, int, LIKELIHOOD ***,
 
7
                     double);
 
8
static void up_ll(LIKELIHOOD *, int, double, LIKELIHOOD *);
 
9
 
 
10
 
 
11
void make_pyramid(LIKELIHOOD **** ll_pym,       /* log likelihood pyramid, ll_pym[scale][i][j][class] */
 
12
                  struct Region *region,        /* specifies image subregion */
 
13
                  int M,        /* number of classes */
 
14
                  double *alpha /* decimation parameters */
 
15
    )
 
16
{
 
17
    int D;
 
18
    int wd, ht;
 
19
    struct Region region_buff;
 
20
 
 
21
    /* save region stucture */
 
22
    copy_reg(region, &region_buff);
 
23
 
 
24
    D = 0;
 
25
    reg_to_wdht(region, &wd, &ht);
 
26
    while ((wd > 2) && (ht > 2)) {
 
27
        G_debug(1, "D = %d  alpha = %f; 1-alpha = %f", D, alpha[D],
 
28
                1 - alpha[D]);
 
29
        decimate(ll_pym[D], region, M, ll_pym[D + 1], alpha[D]);
 
30
        dec_reg(region);
 
31
        reg_to_wdht(region, &wd, &ht);
 
32
        D++;
 
33
    }
 
34
 
 
35
    /* restore region structure */
 
36
    copy_reg(&region_buff, region);
 
37
}
 
38
 
 
39
static void decimate(
 
40
                        /* decimate statistics ll1 to form ll2 */
 
41
                        LIKELIHOOD *** ll1,     /* log likelihood ll1[i][j][class], at fine resolution */
 
42
                        struct Region *region1, /* specifies image subregion */
 
43
                        int M,  /* number of classes */
 
44
                        LIKELIHOOD *** ll2,     /* loglikelihood ll1[i][j][class], at coarse resolution */
 
45
                        double alpha    /* transition parameter */
 
46
    )
 
47
{
 
48
    struct Region *region2, reg_spc;    /* coarse resolution region */
 
49
    int wflag, hflag;           /* flags indicate odd number of pixels */
 
50
    int i, j, m;
 
51
    LIKELIHOOD *node;           /* coarse resolution point */
 
52
    LIKELIHOOD *pt1, *pt2, *pt3, *pt4;  /* fine resolution neighbors */
 
53
 
 
54
    region2 = &reg_spc;
 
55
 
 
56
    copy_reg(region1, region2);
 
57
    dec_reg(region2);
 
58
 
 
59
    wflag = region1->xmax & 1;
 
60
    hflag = region1->ymax & 1;
 
61
 
 
62
    for (i = region2->ymin; i < region2->ymax; i++)
 
63
        for (j = region2->xmin; j < region2->xmax; j++) {
 
64
            pt1 = ll1[2 * i][2 * j];
 
65
            pt2 = ll1[2 * i][2 * j + 1];
 
66
            pt3 = ll1[2 * i + 1][2 * j];
 
67
            pt4 = ll1[2 * i + 1][2 * j + 1];
 
68
 
 
69
            node = ll2[i][j];
 
70
            for (m = 0; m < M; m++)
 
71
                node[m] = 0.0;
 
72
 
 
73
            up_ll(pt1, M, alpha, node);
 
74
            up_ll(pt2, M, alpha, node);
 
75
            up_ll(pt3, M, alpha, node);
 
76
            up_ll(pt4, M, alpha, node);
 
77
        }
 
78
 
 
79
    if (wflag) {
 
80
        for (i = region2->ymin; i < region2->ymax; i++) {
 
81
            node = ll2[i][region2->xmax - 1];
 
82
            for (m = 0; m < M; m++)
 
83
                node[m] = 0.0;
 
84
 
 
85
            pt1 = ll1[2 * i][region1->xmax - 1];
 
86
            pt2 = ll1[2 * i + 1][region1->xmax - 1];
 
87
            up_ll(pt1, M, alpha, node);
 
88
            up_ll(pt2, M, alpha, node);
 
89
        }
 
90
    }
 
91
 
 
92
    if (hflag) {
 
93
        for (j = region2->xmin; j < region2->xmax; j++) {
 
94
            node = ll2[region2->ymax - 1][j];
 
95
            for (m = 0; m < M; m++)
 
96
                node[m] = 0.0;
 
97
 
 
98
            pt1 = ll1[region1->ymax - 1][2 * j];
 
99
            pt2 = ll1[region1->ymax - 1][2 * j + 1];
 
100
            up_ll(pt1, M, alpha, node);
 
101
            up_ll(pt2, M, alpha, node);
 
102
        }
 
103
    }
 
104
 
 
105
    if (hflag && wflag) {
 
106
        node = ll2[region2->ymax - 1][region2->xmax - 1];
 
107
        for (m = 0; m < M; m++)
 
108
            node[m] = 0.0;
 
109
 
 
110
        pt1 = ll1[region1->ymax - 1][region1->xmax - 1];
 
111
        up_ll(pt1, M, alpha, node);
 
112
    }
 
113
}
 
114
 
 
115
 
 
116
static void up_ll(LIKELIHOOD * pt1,     /* array of log likelihood values, pt1[class] */
 
117
                  int M, double alpha, LIKELIHOOD * pt2 /* array of log likelihood values, pt2[class] */
 
118
    )
 
119
{
 
120
    static int m;
 
121
    static double sum, max, cprob[256];
 
122
 
 
123
    if (alpha != 1.0) {
 
124
        max = pt1[0];
 
125
        for (m = 1; m < M; m++)
 
126
            if (pt1[m] > max)
 
127
                max = pt1[m];
 
128
 
 
129
        sum = 0;
 
130
        for (m = 0; m < M; m++) {
 
131
            cprob[m] = exp(pt1[m] - max);
 
132
            sum += cprob[m];
 
133
        }
 
134
 
 
135
        sum = ((1 - alpha) * sum) / M;
 
136
 
 
137
        for (m = 0; m < M; m++)
 
138
            pt2[m] += log(sum + alpha * cprob[m]) + max;
 
139
    }
 
140
    else {
 
141
        for (m = 0; m < M; m++)
 
142
            pt2[m] += pt1[m];
 
143
    }
 
144
}
 
145
 
 
146
char ***get_pyramid(int w0, int h0, size_t size)
 
147
{
 
148
    char ***pym;
 
149
    int D, wn, hn;
 
150
 
 
151
    D = levels(w0, h0);
 
152
    pym = (char ***)G_malloc((D + 1) * sizeof(char *));
 
153
 
 
154
    D = 0;
 
155
    wn = w0;
 
156
    hn = h0;
 
157
    pym[D] = (char **)get_img(wn, hn, size);
 
158
    while ((wn > 2) && (hn > 2)) {
 
159
        D++;
 
160
        wn = wn / 2;
 
161
        hn = hn / 2;
 
162
        pym[D] = (char **)get_img(wn, hn, size);
 
163
    }
 
164
 
 
165
    return (pym);
 
166
}
 
167
 
 
168
 
 
169
void free_pyramid(char *pym, int wd, int ht)
 
170
{
 
171
    unsigned char ***pt;
 
172
    int i, D;
 
173
 
 
174
    pt = (unsigned char ***)pym;
 
175
    D = levels(wd, ht);
 
176
    for (i = 0; i <= D; i++)
 
177
        free_img(pt[i]);
 
178
    G_free(pym);
 
179
}
 
180
 
 
181
char ****get_cubic_pyramid(int w0, int h0, int M, size_t size)
 
182
{
 
183
    char ****pym;
 
184
    int D, wn, hn;
 
185
 
 
186
    D = levels(w0, h0);
 
187
    pym = (char ****)G_malloc((D + 1) * sizeof(char *));
 
188
 
 
189
    D = 0;
 
190
    wn = w0;
 
191
    hn = h0;
 
192
    pym[D] = (char ***)multialloc(size, 3, hn, wn, M);
 
193
    while ((wn > 2) && (hn > 2)) {
 
194
        D++;
 
195
        wn = wn / 2;
 
196
        hn = hn / 2;
 
197
        pym[D] = (char ***)multialloc(size, 3, hn, wn, M);
 
198
    }
 
199
 
 
200
    return (pym);
 
201
}
 
202
 
 
203
 
 
204
void free_cubic_pyramid(char *pym, int wd, int ht, int M)
 
205
{
 
206
    int i, D;
 
207
    char **pt;
 
208
 
 
209
    pt = (char **)pym;
 
210
    D = levels(wd, ht);
 
211
    for (i = 0; i <= D; i++)
 
212
        multifree(pt[i], M);
 
213
    G_free(pym);
 
214
}
 
215
 
 
216
int levels(int wd, int ht)
 
217
{
 
218
    int D;
 
219
 
 
220
    D = 0;
 
221
    while ((wd > 2) && (ht > 2)) {
 
222
        D++;
 
223
        wd = wd / 2;
 
224
        ht = ht / 2;
 
225
    }
 
226
    return (D);
 
227
}