1 #include <stdlib.h>
2 #include <stdio.h>
3 #include <string.h>
4 #include <math.h>
5 #include <limits.h>
6 #include <float.h>
7 #include <grass/gis.h>
8 #include <grass/raster.h>
9 #include "tinf.h"
10
11 /* get the slope between two cells and return a slope direction */
check(CELL newdir,CELL * dir,void * center,void * edge,double cnst,double * oldslope)12 void check(CELL newdir, CELL * dir, void *center, void *edge, double cnst,
13 double *oldslope)
14 {
15 double newslope;
16
17 /* always discharge to a null boundary */
18 if (is_null(edge)) {
19 *oldslope = DBL_MAX;
20 *dir = newdir;
21 }
22 else {
23 newslope = slope(center, edge, cnst);
24 if (newslope == *oldslope) {
25 *dir += newdir;
26 }
27 else if (newslope > *oldslope) {
28 *oldslope = newslope;
29 *dir = newdir;
30 }
31 }
32
33 return;
34
35 }
36
37 /* process one row, filling single-cell pits */
fill_row(int nl,int ns,struct band3 * bnd)38 int fill_row(int nl, int ns, struct band3 *bnd)
39 {
40 int j, offset, inc, rc;
41 void *min;
42 char *center;
43 char *edge;
44
45 inc = bpe();
46
47 min = G_malloc(bpe());
48
49 rc = 0;
50 for (j = 1; j < ns - 1; j += 1) {
51 offset = j * bpe();
52 center = bnd->b[1] + offset;
53 if (is_null(center))
54 return rc;
55
56 edge = bnd->b[0] + offset;
57 min = edge - inc;
58 min = get_min(min, edge);
59 min = get_min(min, edge + inc);
60
61 min = get_min(min, center - inc);
62 min = get_min(min, center + inc);
63
64 edge = bnd->b[2] + offset;
65 min = get_min(min, edge - inc);
66 min = get_min(min, edge);
67 min = get_min(min, edge + inc);
68
69 if (get_min(center, min) == center) {
70 rc = 1;
71 memcpy(center, min, bpe());
72 }
73
74 }
75 return rc;
76 }
77
78 /* determine the flow direction at each cell on one row */
build_one_row(int i,int nl,int ns,struct band3 * bnd,CELL * dir)79 void build_one_row(int i, int nl, int ns, struct band3 *bnd, CELL * dir)
80 {
81 int j, inc;
82 size_t offset;
83 CELL sdir;
84 double curslope;
85 char *center;
86 char *edge;
87
88 inc = bpe();
89
90 for (j = 0; j < ns; j += 1) {
91 offset = j * bpe();
92 center = bnd->b[1] + offset;
93 if (is_null(center)) {
94 Rast_set_c_null_value(dir + j, 1);
95 continue;
96 }
97
98 sdir = 0;
99 curslope = HUGE_VAL;
100 if (i == 0) {
101 sdir = 128;
102 }
103 else if (i == nl - 1) {
104 sdir = 8;
105 }
106 else if (j == 0) {
107 sdir = 32;
108 }
109 else if (j == ns - 1) {
110 sdir = 2;
111 }
112 else {
113 curslope = -HUGE_VAL;
114
115 /* check one row back */
116 edge = bnd->b[0] + offset;
117 check(64, &sdir, center, edge - inc, 1.4142136, &curslope);
118 check(128, &sdir, center, edge, 1., &curslope);
119 check(1, &sdir, center, edge + inc, 1.4142136, &curslope);
120
121 /* check this row */
122 check(32, &sdir, center, center - inc, 1., &curslope);
123 check(2, &sdir, center, center + inc, 1., &curslope);
124
125 /* check one row forward */
126 edge = bnd->b[2] + offset;
127 check(16, &sdir, center, edge - inc, 1.4142136, &curslope);
128 check(8, &sdir, center, edge, 1., &curslope);
129 check(4, &sdir, center, edge + inc, 1.4142136, &curslope);
130 }
131
132 if (curslope == 0.)
133 sdir = -sdir;
134 else if (curslope < 0.)
135 sdir = -256;
136 dir[j] = sdir;
137 }
138 return;
139 }
140
filldir(int fe,int fd,int nl,struct band3 * bnd)141 void filldir(int fe, int fd, int nl, struct band3 *bnd)
142 {
143 int i;
144 size_t bufsz;
145 CELL *dir;
146
147 /* fill single-cell depressions, except on outer rows and columns */
148 lseek(fe, 0, SEEK_SET);
149 advance_band3(fe, bnd);
150 advance_band3(fe, bnd);
151 for (i = 1; i < nl - 1; i += 1) {
152 lseek(fe, (off_t) (i + 1) * bnd->sz, SEEK_SET);
153 advance_band3(fe, bnd);
154 if (fill_row(nl, bnd->ns, bnd)) {
155 lseek(fe, (off_t) i * bnd->sz, SEEK_SET);
156 write(fe, bnd->b[1], bnd->sz);
157 }
158 }
159 /* why on the last row? it's an outer row */
160 #if 0
161 advance_band3(0, bnd);
162 if (fill_row(nl, bnd->ns, bnd)) {
163 lseek(fe, (off_t) i * bnd->sz, SEEK_SET);
164 write(fe, bnd->b[1], bnd->sz);
165 }
166 #endif
167 /* determine the flow direction in each cell. On outer rows and columns
168 * the flow direction is always directly out of the map */
169
170 dir = G_calloc(bnd->ns, sizeof(CELL));
171 bufsz = bnd->ns * sizeof(CELL);
172
173 lseek(fe, 0, SEEK_SET);
174 lseek(fd, 0, SEEK_SET);
175 advance_band3(fe, bnd);
176 for (i = 0; i < nl; i += 1) {
177 advance_band3(fe, bnd);
178 build_one_row(i, nl, bnd->ns, bnd, dir);
179 write(fd, dir, bufsz);
180 }
181 /* why this extra row ? */
182 #if 0
183 advance_band3(fe, bnd);
184 build_one_row(i, nl, bnd->ns, bnd, dir);
185 write(fd, dir, bufsz);
186 #endif
187
188 G_free(dir);
189
190 return;
191 }
192