1
2 /*****************************************************************************/
3
4 /*** ***/
5
6 /*** process() ***/
7
8 /*** Reads in a raster map row by row for processing. ***/
9
10 /*** Jo Wood, Project ASSIST, V1.0 7th February 1993 ***/
11
12 /*** ***/
13
14 /*****************************************************************************/
15
16 #include <stdlib.h>
17 #include <grass/gis.h>
18 #include <grass/raster.h>
19 #include <grass/glocale.h>
20 #include <grass/gmath.h>
21 #include "param.h"
22 #include "nrutil.h"
23
24
process(void)25 void process(void)
26 {
27
28 /*--------------------------------------------------------------------------*/
29 /* INITIALISE */
30
31 /*--------------------------------------------------------------------------*/
32
33
34 DCELL *row_in, /* Buffer large enough to hold `wsize' */
35 *row_out = NULL, /* raster rows. When GRASS reads in a */
36 /* raster row, each element is of type */
37 /* DCELL */
38 *window_ptr, /* Stores local terrain window. */
39 centre, /* Elevation of central cell in window. */
40 *window_cell; /* Stores last read cell in window. */
41
42 CELL *featrow_out = NULL; /* store features in CELL */
43
44 struct Cell_head region; /* Structure to hold region information */
45
46 int nrows, /* Will store the current number of */
47 ncols, /* rows and columns in the raster. */
48 row, col, /* Counts through each row and column */
49 /* of the input raster. */
50 wind_row, /* Counts through each row and column */
51 wind_col, /* of the local neighbourhood window. */
52 *index_ptr; /* Row permutation vector for LU decomp. */
53
54 double **normal_ptr, /* Cross-products matrix. */
55 *obs_ptr, /* Observed vector. */
56 temp; /* Unused */
57
58 double *weight_ptr; /* Weighting matrix for observed values. */
59 int found_null; /* If null was found among window cells */
60
61 /*--------------------------------------------------------------------------*/
62 /* GET RASTER AND WINDOW DETAILS */
63
64 /*--------------------------------------------------------------------------*/
65
66 G_get_window(®ion); /* Fill out the region structure (the */
67 /* geographical limits etc.) */
68
69 nrows = Rast_window_rows(); /* Find out the number of rows and */
70 ncols = Rast_window_cols(); /* columns of the raster. */
71
72
73 if ((region.ew_res / region.ns_res >= 1.01) || /* If EW and NS resolns are */
74 (region.ns_res / region.ew_res >= 1.01)) { /* >1% different, warn user. */
75 G_warning(_("E-W and N-S grid resolutions are different. Taking average."));
76 resoln = (region.ns_res + region.ew_res) / 2;
77 }
78 else
79 resoln = region.ns_res;
80
81
82 /*--------------------------------------------------------------------------*/
83 /* RESERVE MEMORY TO HOLD Z VALUES AND MATRICES */
84
85 /*--------------------------------------------------------------------------*/
86
87 row_in = (DCELL *) G_malloc(ncols * sizeof(DCELL) * wsize);
88 /* Reserve `wsize' rows of memory. */
89
90 if (mparam != FEATURE) {
91 row_out = Rast_allocate_buf(DCELL_TYPE); /* Initialise output row buffer. */
92 Rast_set_d_null_value(row_out, ncols);
93 }
94 else {
95 featrow_out = Rast_allocate_buf(CELL_TYPE); /* Initialise output row buffer. */
96 Rast_set_c_null_value(featrow_out, ncols);
97 }
98
99 window_ptr = (DCELL *) G_malloc(SQR(wsize) * sizeof(DCELL));
100 /* Reserve enough memory for local wind. */
101
102 weight_ptr = (double *)G_malloc(SQR(wsize) * sizeof(double));
103 /* Reserve enough memory weights matrix. */
104
105 normal_ptr = dmatrix(0, 5, 0, 5); /* Allocate memory for 6*6 matrix */
106 index_ptr = ivector(0, 5); /* and for 1D vector holding indices */
107 obs_ptr = dvector(0, 5); /* and for 1D vector holding observed z */
108
109
110 /* ---------------------------------------------------------------- */
111 /* - CALCULATE LEAST SQUARES COEFFICIENTS - */
112 /* ---------------------------------------------------------------- */
113
114 /*--- Calculate weighting matrix. ---*/
115
116 find_weight(weight_ptr);
117
118 /* Initial coefficients need only be found once since they are
119 constant for any given window size. The only element that
120 changes is the observed vector (RHS of normal equations). */
121
122 /*--- Find normal equations in matrix form. ---*/
123
124 find_normal(normal_ptr, weight_ptr);
125
126
127 /*--- Apply LU decomposition to normal equations. ---*/
128
129 if (constrained) {
130 G_ludcmp(normal_ptr, 5, index_ptr, &temp);
131 /* To constrain the quadtratic
132 through the central cell, ignore
133 the calculations involving the
134 coefficient f. Since these are
135 all in the last row and column of
136 the matrix, simply redimension. */
137 /* disp_matrix(normal_ptr,obs_ptr,obs_ptr,5);
138 */
139 }
140
141 else {
142 G_ludcmp(normal_ptr, 6, index_ptr, &temp);
143 /* disp_matrix(normal_ptr,obs_ptr,obs_ptr,6);
144 */
145 }
146
147
148 /*--------------------------------------------------------------------------*/
149 /* PROCESS INPUT RASTER AND WRITE OUT RASTER LINE BY LINE */
150
151 /*--------------------------------------------------------------------------*/
152
153 if (mparam != FEATURE)
154 for (wind_row = 0; wind_row < EDGE; wind_row++)
155 Rast_put_row(fd_out, row_out, DCELL_TYPE); /* Write out the edge cells as NULL. */
156 else
157 for (wind_row = 0; wind_row < EDGE; wind_row++)
158 Rast_put_row(fd_out, featrow_out, CELL_TYPE); /* Write out the edge cells as NULL. */
159
160 for (wind_row = 0; wind_row < wsize - 1; wind_row++)
161 Rast_get_row(fd_in, row_in + (wind_row * ncols), wind_row,
162 DCELL_TYPE);
163 /* Read in enough of the first rows to */
164 /* allow window to be examined. */
165
166 for (row = EDGE; row < (nrows - EDGE); row++) {
167 G_percent(row + 1, nrows - EDGE, 2);
168
169 Rast_get_row(fd_in, row_in + ((wsize - 1) * ncols), row + EDGE,
170 DCELL_TYPE);
171
172 for (col = EDGE; col < (ncols - EDGE); col++) {
173 /* Find central z value */
174 centre = *(row_in + EDGE * ncols + col);
175 /* Test for no data and propagate */
176 if (Rast_is_d_null_value(¢re)) {
177 if (mparam == FEATURE)
178 Rast_set_c_null_value(featrow_out + col, 1);
179 else {
180 Rast_set_d_null_value(row_out + col, 1);
181 }
182 continue;
183 }
184 found_null = FALSE;
185 for (wind_row = 0; wind_row < wsize; wind_row++) {
186 if (found_null)
187 break;
188 for (wind_col = 0; wind_col < wsize; wind_col++) {
189
190 /* Express all window values relative */
191 /* to the central elevation. */
192 window_cell = row_in + (wind_row * ncols) + col + wind_col - EDGE;
193 /* Test for no data and propagate */
194 if (Rast_is_d_null_value(window_cell)) {
195 if (mparam == FEATURE)
196 Rast_set_c_null_value(featrow_out + col, 1);
197 else {
198 Rast_set_d_null_value(row_out + col, 1);
199 }
200 found_null = TRUE;
201 break;
202 }
203 *(window_ptr + (wind_row * wsize) + wind_col) =
204 *(window_cell) - centre;
205 }
206 }
207
208 if (found_null)
209 continue;
210
211 /*--- Use LU back substitution to solve normal equations. ---*/
212 find_obs(window_ptr, obs_ptr, weight_ptr);
213
214 /* disp_wind(window_ptr);
215 disp_matrix(normal_ptr,obs_ptr,obs_ptr,6);
216 */
217
218 if (constrained) {
219 G_lubksb(normal_ptr, 5, index_ptr, obs_ptr);
220 /*
221 disp_matrix(normal_ptr,obs_ptr,obs_ptr,5);
222 */
223 }
224
225 else {
226 G_lubksb(normal_ptr, 6, index_ptr, obs_ptr);
227 /*
228 disp_matrix(normal_ptr,obs_ptr,obs_ptr,6);
229 */
230
231 }
232
233 /*--- Calculate terrain parameter based on quad. coefficients. ---*/
234 if (mparam == FEATURE)
235 *(featrow_out + col) = (CELL) feature(obs_ptr);
236 else
237 *(row_out + col) = param(mparam, obs_ptr);
238
239 if (mparam == ELEV)
240 *(row_out + col) += centre; /* Add central elevation back */
241 }
242
243 if (mparam != FEATURE)
244 Rast_put_row(fd_out, row_out, DCELL_TYPE); /* Write the row buffer to the output */
245 /* raster. */
246 else /* write FEATURE to CELL */
247 Rast_put_row(fd_out, featrow_out, CELL_TYPE); /* Write the row buffer to the output */
248 /* raster. */
249
250 /* 'Shuffle' rows down one, and read in */
251 /* one new row. */
252 for (wind_row = 0; wind_row < wsize - 1; wind_row++)
253 for (col = 0; col < ncols; col++)
254 *(row_in + (wind_row * ncols) + col) =
255 *(row_in + ((wind_row + 1) * ncols) + col);
256 }
257
258 if (mparam != FEATURE)
259 Rast_set_d_null_value(row_out, ncols);
260 else
261 Rast_set_c_null_value(featrow_out, ncols);
262 for (wind_row = 0; wind_row < EDGE; wind_row++) {
263 if (mparam != FEATURE)
264 Rast_put_row(fd_out, row_out, DCELL_TYPE); /* Write out the edge cells as NULL. */
265 else
266 Rast_put_row(fd_out, featrow_out, CELL_TYPE); /* Write out the edge cells as NULL. */
267 }
268
269 /*--------------------------------------------------------------------------*/
270 /* FREE MEMORY USED TO STORE RASTER ROWS, LOCAL WINDOW AND MATRICES */
271
272 /*--------------------------------------------------------------------------*/
273
274 G_free(row_in);
275 if (mparam != FEATURE)
276 G_free(row_out);
277 else
278 G_free(featrow_out);
279
280 G_free(window_ptr);
281 free_dmatrix(normal_ptr, 0, 5, 0, 5);
282 free_dvector(obs_ptr, 0, 5);
283 free_ivector(index_ptr, 0, 5);
284 }
285