1
2 /***********************************************************************
3 *
4 * spread.c in ../r.spread
5 *
6 * This is a raster version of Dijkstra's shortest path algorithm
7 * that is suited for simulating elliptical spread phenomena. It
8 * 1) starts from each spread origin (stored in
9 * a linked list of type costHa) - spread();
10 * 2) selects appropriate cells as links for the current
11 * spread cell and stored in a linked list of type
12 * cell_ptrHa - select() ;
13 * A) caculates the cumulative cost (time) of the
14 * end cell of each link - calculate();
15 * B) compares this new cumulative cost (time) with
16 * the previously computed cumulative time/cost,
17 * if there is any, of the same cell - update();
18 * C) puts this cell into a min-heap and puts the
19 * new cumulative cost (time) together with UTM
20 * coordinates in the cumulative cost (time)
21 * map, x (East) map and y (North) map if
22 * there is no previous cumulative cost (time);
23 * otherwise, if the new cumulative cost (time)
24 * is less, replaces with it both in the heap
25 * and the output maps - update().
26 * 3) gets the first cell in the min-heap, which is the
27 * cell with the least cumulative cost (time), and
28 * repeats Step 2 until the heap is empty or desired
29 * simulated cumulative cost (time) is reached - spread().
30 *
31 ***********************************************************************/
32 #include <stdio.h>
33 #include <stdlib.h>
34 #include <math.h>
35 #include <grass/gis.h>
36 #include <grass/raster.h>
37 #include "cmd_line.h"
38 #include "costHa.h"
39 #include "cell_ptrHa.h"
40 #include "local_proto.h"
41
42 #ifndef PI
43 #define PI M_PI
44 #endif
45 #define DATA(map, r, c) (map)[(r) * ncols + (c)]
46
47 /*#define DEBUG */
48
49 extern CELL *map_max, *map_base, *map_dir, *map_visit;
50 extern CELL *map_x_out, *map_y_out;
51 extern float *map_out;
52 extern int BARRIER;
53 extern int nrows, ncols;
54
55 /* extern float PI; */
56 extern long heap_len;
57 extern struct costHa *heap;
58 extern struct Cell_head window;
59
60 struct cell_ptrHa *front_cell = NULL, *rear_cell = NULL;
61
spread(void)62 void spread(void)
63 {
64 float min_cost;
65 int ros_max, ros_base, dir;
66 int row, col;
67 int cell_count = 0, ncells = 0;
68 struct cell_ptrHa *to_cell, *old_to_cell;
69 struct costHa *pres_cell;
70
71 /* initialize using arbitrary value, this value is never used except for debug */
72 min_cost = 0;
73
74 ncells = nrows * ncols;
75 G_message
76 ("Finding spread time - number of cells visited in percentage ... %3d%%",
77 0);
78 pres_cell = (struct costHa *)G_malloc(sizeof(struct costHa));
79 get_minHa(heap, pres_cell, heap_len);
80 G_debug(2, "begin spread: cost(%d,%d)=%f", pres_cell->row, pres_cell->col,
81 pres_cell->min_cost);
82 G_debug(2,
83 " heap_len=%ld pres_cell->min_cost=%f time_lag=%d",
84 heap_len, pres_cell->min_cost, time_lag);
85 while (heap_len-- > 0 && pres_cell->min_cost < init_time + time_lag + 1.0) {
86 ros_max = DATA(map_max, pres_cell->row, pres_cell->col);
87 ros_base = DATA(map_base, pres_cell->row, pres_cell->col);
88 dir = DATA(map_dir, pres_cell->row, pres_cell->col);
89
90 /*Select end cells of links of the present cell */
91 select_linksB(pres_cell, least / 2, comp_dens);
92
93 #ifdef DEBUG
94 to_cell = front_cell;
95 while (to_cell != NULL) {
96 printf("(%d,%d) ", to_cell->row, to_cell->col);
97 to_cell = to_cell->next;
98 }
99 #endif
100 /*Get a cell in the list each time, and compute culmulative costs
101 *via the current spread cell*/
102 to_cell = front_cell;
103 while (to_cell != NULL) {
104 /*calculate cumulative costs,
105 *function returns -1 if detected a barrier */
106 if (cumulative
107 (pres_cell, to_cell, ros_max, ros_base, dir,
108 &min_cost) == -1) {
109 old_to_cell = to_cell;
110 to_cell = to_cell->next;
111 front_cell = to_cell;
112 G_free(old_to_cell);
113 continue;
114 }
115
116 G_debug(2, " finish a link: cost(%d,%d)->(%d,%d)=%f",
117 pres_cell->row, pres_cell->col, to_cell->row,
118 to_cell->col, min_cost);
119 /*update the cumulative time/cost */
120 update(pres_cell, to_cell->row, to_cell->col, to_cell->angle,
121 min_cost);
122 old_to_cell = to_cell;
123 to_cell = to_cell->next;
124 front_cell = to_cell;
125 G_free(old_to_cell);
126 }
127
128 /*compute spotting fires */
129 if (spotting)
130 spot(pres_cell, dir);
131
132 /*mark a visited cell */
133 DATA(map_visit, pres_cell->row, pres_cell->col) = YES;
134 #if 0
135 if (display)
136 draw_a_cell(pres_cell->row, pres_cell->col,
137 (int)pres_cell->min_cost);
138 #endif
139
140 cell_count++;
141 if ((100 * cell_count / ncells) % 2 == 0 &&
142 (100 * (cell_count + (int)(0.009 * ncells)) / ncells) % 2 == 0) {
143 G_percent(cell_count, ncells, 2);
144 }
145
146 get_minHa(heap, pres_cell, heap_len);
147 G_debug(2,
148 "in while: heap_len=%ld pres_cell->min_cost=%f time_lag=%d",
149 heap_len, pres_cell->min_cost, time_lag);
150 } /*end 'while (heap_len-- >0)' */
151 G_free(pres_cell);
152
153 /*Assign min_cost values to un-reached area */
154 for (row = 0; row < nrows; row++) {
155 for (col = 0; col < ncols; col++) {
156 if (!DATA(map_visit, row, col)) {
157 DATA(map_out, row, col) = (float)BARRIER;
158 if (x_out)
159 DATA(map_x_out, row, col) = 0;
160 if (y_out)
161 DATA(map_y_out, row, col) = 0;
162 }
163 }
164 }
165 G_debug(2, "end spread");
166 } /*end spread () */
167
168
169 /******* function computing cumulative spread time/cost, ***************
170 ******* good for both adjacent cell links and non-adjacent cell links */
171
172 int
cumulative(struct costHa * pres_cell,struct cell_ptrHa * to_cell,int ros_max,int ros_base,int dir,float * min_cost)173 cumulative(struct costHa *pres_cell, struct cell_ptrHa *to_cell,
174 int ros_max, int ros_base, int dir, float *min_cost)
175 {
176 float ros, xros, cost;
177 float xstep_len;
178 float cos_angle, sin_angle;
179 int xrow, xcol, xsteps, count;
180
181 /*most of the actions below calculate the cumulative time/cost,
182 *from the current spread cell, of the end cell of one link*/
183
184 sin_angle = sin(to_cell->angle);
185 cos_angle = cos(to_cell->angle);
186
187 if (abs(pres_cell->row - to_cell->row) >
188 abs(pres_cell->col - to_cell->col)) {
189 xsteps = abs(pres_cell->row - to_cell->row);
190 xstep_len = 1 / cos_angle;
191 if (xstep_len < 0.0)
192 xstep_len = -xstep_len;
193 }
194 else {
195 xsteps = abs(pres_cell->col - to_cell->col);
196 xstep_len = 1 / sin_angle;
197 if (xstep_len < 0.0)
198 xstep_len = -xstep_len;
199 }
200
201 /*ROS value based on a 'from_cell', (elliptical cases) */
202 ros =
203 ros_base / (1 -
204 (1 - ros_base / (float)ros_max) * cos(to_cell->angle -
205 dir % 360 * PI /
206 180));
207
208 /*the next cell */
209 xrow = pres_cell->row - xstep_len * cos_angle + 0.5;
210 xcol = pres_cell->col + xstep_len * sin_angle + 0.5;
211
212 cost = 0.0;
213 count = 1;
214 while (count <= xsteps) {
215 /*Can't go through a barrer in a path */
216 if (DATA(map_base, xrow, xcol) <= 0)
217 return -1;
218
219 /*ROS value based on current 'to_cell', (elliptical cases) */
220 xros =
221 DATA(map_base, xrow,
222 xcol) / (1 - (1 -
223 DATA(map_base, xrow,
224 xcol) / (float)DATA(map_max, xrow,
225 xcol)) *
226 cos(to_cell->angle -
227 DATA(map_dir, xrow, xcol) % 360 * PI / 180));
228 /*Calculate cost to this cell */
229 cost =
230 cost + 0.5 * (xstep_len * window.ns_res / ros +
231 xstep_len * window.ns_res / xros);
232
233 /*Update temp cell along the path, and counter */
234 ros = xros;
235 xrow = pres_cell->row - count * xstep_len * cos_angle + 0.5;
236 xcol = pres_cell->col + count * xstep_len * sin_angle + 0.5;
237 count++;
238 } /*end'while (count<= ..)' */
239 G_debug(2, " in cumulatvie() cost=%.2f pre min_cost=%.2f",
240 cost, *min_cost);
241 /*from the origin, cumulative time/cost of the end cell of one link */
242 *min_cost = pres_cell->min_cost + cost;
243 G_debug(2, " in cumulatvie() post min_cost=%.2f",
244 *min_cost);
245
246 return 0;
247 }
248
249
250 /****** function for updating the cumulative cost/time, possibaly ********
251 ****** back path x,y coordinates, both in the output(s) and the heap ********/
252
253 void
update(struct costHa * pres_cell,int row,int col,double angle,float min_cost)254 update(struct costHa *pres_cell, int row, int col, double angle,
255 float min_cost)
256 {
257 if (DATA(map_out, row, col) < -1.0) {
258 G_debug(2, " insert: out(%d,%d)=%f min_cost=%f", row, col,
259 DATA(map_out, row, col), min_cost);
260 DATA(map_out, row, col) = min_cost;
261 if (x_out)
262 DATA(map_x_out, row, col) = pres_cell->col;
263 if (y_out)
264 DATA(map_y_out, row, col) = pres_cell->row;
265
266 insertHa(min_cost, angle, row, col, heap, &heap_len);
267 #if 0
268 if (display && min_cost < init_time + time_lag + 1.0)
269 draw_a_burning_cell(row, col);
270 #endif
271 }
272 else {
273 if (DATA(map_out, row, col) > min_cost + 0.001) {
274 G_debug(2, " replace: out(%d,%d)=%f min_cost=%f", row, col,
275 DATA(map_out, row, col), min_cost);
276 DATA(map_out, row, col) = min_cost;
277 if (x_out)
278 DATA(map_x_out, row, col) = pres_cell->col;
279 if (y_out)
280 DATA(map_y_out, row, col) = pres_cell->row;
281
282 replaceHa(min_cost, angle, row, col, heap, &heap_len);
283 #if 0
284 if (display && min_cost < init_time + time_lag + 1.0)
285 draw_a_burning_cell(row, col);
286 #endif
287 }
288 }
289 }
290