1 /* File: correction.c
2 *
3 * AUTHOR: E. Jorge Tizado, Spain 2010
4 *
5 * COPYRIGHT: (c) 2007-10 E. Jorge Tizado
6 * This program is free software under the GNU General Public
7 * License (>=v2). Read the file COPYING that comes with GRASS
8 * for details.
9 */
10
11 #include <stdio.h>
12 #include <stdlib.h>
13 #include <math.h>
14 #include <float.h>
15 #include <unistd.h>
16 #include <grass/gis.h>
17 #include <grass/raster.h>
18 #include <grass/glocale.h>
19
20 #include "local_proto.h"
21
eval_tcor(int method,Gfile * out,Gfile * cosi,Gfile * band,double zenith,int do_scale)22 void eval_tcor(int method, Gfile * out, Gfile * cosi, Gfile * band,
23 double zenith, int do_scale)
24 {
25 int row, col, nrows, ncols;
26 void *pref, *pcos;
27
28 double cos_z, cos_i, ref_i, result;
29 double n, sx, sxx, sy, sxy, tx, ty;
30 double a, m, cka, ckb, kk;
31
32 double imin, imax, omin, omax, factor; /* for scaling to input */
33
34 nrows = Rast_window_rows();
35 ncols = Rast_window_cols();
36
37 cos_z = cos(D2R * zenith);
38
39 imin = omin = DBL_MAX;
40 imax = omax = -DBL_MAX;
41 factor = 1;
42
43 /* Calculating regression */
44 if (method > NON_LAMBERTIAN) {
45 n = sx = sxx = sy = sxy = 0.;
46 for (row = 0; row < nrows; row++) {
47 G_percent(row, nrows, 2);
48
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;
54
55 for (col = 0; col < ncols; col++) {
56
57 cos_i = Rast_get_d_value(pcos, cosi->type);
58
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
68 switch (method) {
69 case MINNAERT:
70 if (cos_i > 0. && cos_z > 0. && ref_i > 0.) {
71 n++;
72 /* tx = log(cos_i / cos_z) */
73 /* cos_z is constant then m not changes */
74 tx = log(cos_i);
75 ty = log(ref_i);
76 sx += tx;
77 sxx += tx * tx;
78 sy += ty;
79 sxy += tx * ty;
80 }
81 break;
82 case C_CORRECT:
83 {
84 n++;
85 sx += cos_i;
86 sxx += cos_i * cos_i;
87 sy += ref_i;
88 sxy += cos_i * ref_i;
89 }
90 break;
91 }
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));
95 }
96 }
97 m = (n == 0.) ? 1. : (n * sxy - sx * sy) / (n * sxx - sx * sx);
98 a = (n == 0.) ? 0. : (sy - m * sx) / n;
99 }
100 /* Calculating Constants */
101 switch (method) {
102 case MINNAERT:
103 cka = ckb = 0.;
104 kk = m;
105 G_message("Minnaert constant = %lf", kk);
106 break;
107 case C_CORRECT:
108 cka = ckb = a / m; /* Richter changes to m/a */
109 kk = 1.;
110 G_message("C-factor constant = %lf (a=%.4f; m=%.4f)", cka, a, m);
111 break;
112 case PERCENT:
113 cka = 2. - cos_z;
114 ckb = 1.;
115 kk = 1.;
116 break;
117 default: /* COSINE */
118 cka = ckb = 0.;
119 kk = 1.;
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 }
160 /* Topographic correction */
161 for (row = 0; row < nrows; row++) {
162 G_percent(row, nrows, 2);
163
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);
171
172 for (col = 0; col < ncols; col++) {
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;
187 G_debug(3,
188 "Old val: %f, cka: %f, cos_i: %f, ckb: %f, kk: %f, New val: %f",
189 ref_i, cka, cos_i, ckb, kk,
190 ((DCELL *) out->rast)[col]);
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));
194 }
195 Rast_put_row(out->fd, out->rast, out->type);
196 }
197 G_percent(1, 1, 1);
198 }
199