← Back to branch summary

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

« back to all changes in this revision

Viewing changes to imagery/i.topo.corr/correction.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.topo.corr/correction.c
11
11
#include <stdio.h>
12
12
#include <stdlib.h>
13
13
#include <math.h>
 
14
#include <float.h>
14
15
#include <unistd.h>
15
16
#include <grass/gis.h>
 
17
#include <grass/raster.h>
16
18
#include <grass/glocale.h>
17
19
 
18
20
#include "local_proto.h"
19
21
 
20
22
void eval_tcor(int method, Gfile * out, Gfile * cosi, Gfile * band,
21
 
               double zenith)
 
23
               double zenith, int do_scale)
22
24
{
23
25
    int row, col, nrows, ncols;
24
26
    void *pref, *pcos;
25
27
 
26
 
    double cos_z, cos_i, ref_i;
 
28
    double cos_z, cos_i, ref_i, result;
27
29
    double n, sx, sxx, sy, sxy, tx, ty;
28
30
    double a, m, cka, ckb, kk;
29
31
 
 
32
    double imin, imax, omin, omax, factor;    /* for scaling to input */
30
33
 
31
 
    nrows = G_window_rows();
32
 
    ncols = G_window_cols();
 
34
    nrows = Rast_window_rows();
 
35
    ncols = Rast_window_cols();
33
36
 
34
37
    cos_z = cos(D2R * zenith);
 
38
    
 
39
    imin = omin = DBL_MAX;
 
40
    imax = omax = -DBL_MAX;
 
41
    factor = 1;
35
42
 
36
43
    /* Calculating regression */
37
44
    if (method > NON_LAMBERTIAN) {
39
46
        for (row = 0; row < nrows; row++) {
40
47
            G_percent(row, nrows, 2);
41
48
 
42
 
            G_get_raster_row(band->fd, band->rast, row, DCELL_TYPE);
43
 
            G_get_raster_row(cosi->fd, cosi->rast, row, cosi->type);
 
49
            Rast_get_row(band->fd, band->rast, row, band->type);
 
50
            Rast_get_row(cosi->fd, cosi->rast, row, cosi->type);
 
51
            
 
52
            pref = band->rast;
 
53
            pcos = cosi->rast;
44
54
 
45
55
            for (col = 0; col < ncols; col++) {
46
 
                switch (cosi->type) {
47
 
                case FCELL_TYPE:
48
 
                    pcos = (void *)((FCELL *) cosi->rast + col);
49
 
                    cos_i = (double)((FCELL *) cosi->rast)[col];
50
 
                    break;
51
 
                case DCELL_TYPE:
52
 
                    pcos = (void *)((DCELL *) cosi->rast + col);
53
 
                    cos_i = (double)((DCELL *) cosi->rast)[col];
54
 
                    break;
55
 
                default:
56
 
                    pcos = NULL;
57
 
                    cos_i = 0.;
58
 
                    break;
59
 
                }
60
 
                pref = (void *)((DCELL *) band->rast + col);
 
56
                
 
57
                cos_i = Rast_get_d_value(pcos, cosi->type);
61
58
 
62
 
                if (!G_is_null_value(pref, DCELL_TYPE) &&
63
 
                    !G_is_null_value(pcos, cosi->type)) {
64
 
                    ref_i = (double)*((DCELL *) pref);
 
59
                if (!Rast_is_null_value(pref, band->type) &&
 
60
                    !Rast_is_null_value(pcos, cosi->type)) {
 
61
                    ref_i = Rast_get_d_value(pref, band->type);
 
62
                    
 
63
                    if (imin > ref_i)
 
64
                        imin = ref_i;
 
65
                    if (imax < ref_i)
 
66
                        imax = ref_i;
 
67
                    
65
68
                    switch (method) {
66
69
                    case MINNAERT:
67
70
                        if (cos_i > 0. && cos_z > 0. && ref_i > 0.) {
87
90
                        break;
88
91
                    }
89
92
                }
 
93
                pref = G_incr_void_ptr(pref, Rast_cell_size(band->type));
 
94
                pcos = G_incr_void_ptr(pcos, Rast_cell_size(cosi->type));
90
95
            }
91
96
        }
92
97
        m = (n == 0.) ? 1. : (n * sxy - sx * sy) / (n * sxx - sx * sx);
97
102
    case MINNAERT:
98
103
        cka = ckb = 0.;
99
104
        kk = m;
100
 
        G_message(_("Minnaert constant = %lf"), kk);
 
105
        G_message("Minnaert constant = %lf", kk);
101
106
        break;
102
107
    case C_CORRECT:
103
108
        cka = ckb = a / m; /* Richter changes to m/a */
104
109
        kk = 1.;
105
 
        G_message(_("C-factor constant = %lf (a=%.4f; m=%.4f)"), cka, a, m);
 
110
        G_message("C-factor constant = %lf (a=%.4f; m=%.4f)", cka, a, m);
106
111
        break;
107
112
    case PERCENT:
108
113
        cka = 2. - cos_z;
113
118
        cka = ckb = 0.;
114
119
        kk = 1.;
115
120
    }
 
121
 
 
122
    if (do_scale) {
 
123
        /* Topographic correction, pass 1 */
 
124
        for (row = 0; row < nrows; row++) {
 
125
            G_percent(row, nrows, 2);
 
126
 
 
127
            Rast_get_row(band->fd, band->rast, row, band->type);
 
128
            Rast_get_row(cosi->fd, cosi->rast, row, cosi->type);
 
129
 
 
130
            pref = band->rast;
 
131
            pcos = cosi->rast;
 
132
            
 
133
            for (col = 0; col < ncols; col++) {
 
134
 
 
135
                cos_i = Rast_get_d_value(pcos, cosi->type);
 
136
 
 
137
                if (!Rast_is_null_value(pref, band->type) &&
 
138
                    !Rast_is_null_value(pcos, cosi->type)) {
 
139
 
 
140
                    ref_i = Rast_get_d_value(pref, band->type);
 
141
                    result = (DCELL) (ref_i * pow((cos_z + cka) /
 
142
                                      (cos_i + ckb), kk));
 
143
                    G_debug(3,
 
144
                            "Old val: %f, cka: %f, cos_i: %f, ckb: %f, kk: %f, New val: %f",
 
145
                            ref_i, cka, cos_i, ckb, kk, result);
 
146
 
 
147
                    if (omin > result)
 
148
                        omin = result;
 
149
                    if (omax < result)
 
150
                        omax = result;
 
151
 
 
152
                }
 
153
                pref = G_incr_void_ptr(pref, Rast_cell_size(band->type));
 
154
                pcos = G_incr_void_ptr(pcos, Rast_cell_size(cosi->type));
 
155
            }
 
156
        }
 
157
        G_percent(1, 1, 1);
 
158
        factor = (imax - imin) / (omax - omin);
 
159
    }
116
160
    /* Topographic correction */
117
161
    for (row = 0; row < nrows; row++) {
118
162
        G_percent(row, nrows, 2);
119
163
 
120
 
        G_get_raster_row(band->fd, band->rast, row, band->type);
121
 
        G_get_raster_row(cosi->fd, cosi->rast, row, cosi->type);
 
164
        Rast_get_row(band->fd, band->rast, row, band->type);
 
165
        Rast_get_row(cosi->fd, cosi->rast, row, cosi->type);
 
166
 
 
167
        pref = band->rast;
 
168
        pcos = cosi->rast;
 
169
        
 
170
        Rast_set_null_value(out->rast, ncols, DCELL_TYPE);
122
171
 
123
172
        for (col = 0; col < ncols; col++) {
124
 
            switch (cosi->type) {
125
 
            case FCELL_TYPE:
126
 
                pcos = (void *)((FCELL *) cosi->rast + col);
127
 
                cos_i = (double)*((FCELL *) pcos);
128
 
                break;
129
 
            case DCELL_TYPE:
130
 
                pcos = (void *)((DCELL *) cosi->rast + col);
131
 
                cos_i = (double)*((DCELL *) pcos);
132
 
                break;
133
 
            default:
134
 
                pcos = NULL;
135
 
                cos_i = 0.;
136
 
                break;
137
 
            }
138
 
            pref = (void *)((DCELL *) band->rast + col);
139
 
 
140
 
            if (pcos == NULL ||
141
 
                G_is_null_value(pref, DCELL_TYPE) ||
142
 
                G_is_null_value(pcos, cosi->type)) {
143
 
                G_set_null_value((DCELL *) out->rast + col, 1, DCELL_TYPE);
144
 
            }
145
 
            else {
146
 
                ref_i = (double)*((DCELL *) pref);
147
 
                ((DCELL *) out->rast)[col] =
148
 
                    (DCELL) (ref_i * pow((cos_z + cka) / (cos_i + ckb), kk));
 
173
 
 
174
            cos_i = Rast_get_d_value(pcos, cosi->type);
 
175
 
 
176
            if (!Rast_is_null_value(pref, band->type) &&
 
177
                !Rast_is_null_value(pcos, cosi->type)) {
 
178
 
 
179
                ref_i = Rast_get_d_value(pref, band->type);
 
180
                result = (DCELL) (ref_i * pow((cos_z + cka) /
 
181
                                  (cos_i + ckb), kk));
 
182
 
 
183
                if (do_scale)
 
184
                    result = (result - omin) * factor + imin;
 
185
 
 
186
                ((DCELL *) out->rast)[col] = result;
149
187
                G_debug(3,
150
188
                        "Old val: %f, cka: %f, cos_i: %f, ckb: %f, kk: %f, New val: %f",
151
189
                        ref_i, cka, cos_i, ckb, kk,
152
190
                        ((DCELL *) out->rast)[col]);
153
191
            }
 
192
            pref = G_incr_void_ptr(pref, Rast_cell_size(band->type));
 
193
            pcos = G_incr_void_ptr(pcos, Rast_cell_size(cosi->type));
154
194
        }
155
 
        G_put_raster_row(out->fd, out->rast, out->type);
 
195
        Rast_put_row(out->fd, out->rast, out->type);
156
196
    }
 
197
    G_percent(1, 1, 1);
157
198
}