1 /****************************************************************************
2 *
3 * MODULE: r.terraflow
4 *
5 * COPYRIGHT (C) 2007 Laura Toma
6 *
7 * This program is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 *****************************************************************************/
18
19 #include <stdlib.h>
20 #include <ctype.h>
21 #include <time.h>
22
23 #include <string>
24 #include "fill.h"
25 #include "common.h"
26 #include "water.h"
27 #include "sortutils.h"
28 #include "streamutils.h"
29 #include "filldepr.h"
30 #include "grid.h"
31
32 /* globals in common.H
33
34 extern statsRecorder *stats; stats file
35 extern userOptions *opt; command-line options
36 extern struct Cell_head *region; header of the region
37 extern dimension_type nrows, ncols;
38 */
39
40
41 #define FILL_DEBUG if(0)
42 #define FILL_SAVEALL if(0)
43
44
45 /* defined in this module */
46 void recordWatersheds(AMI_STREAM<labelElevType> *labeledWater);
47 long assignDirections(AMI_STREAM<plateauStats> *statstr,
48 AMI_STREAM<plateauType> *platstr,
49 AMI_STREAM<waterType> *waterstr);
50 void assignFinalDirections(AMI_STREAM<plateauStats> *statstr,
51 AMI_STREAM<plateauType> *platstr,
52 AMI_STREAM<waterType> *waterstr);
53 template<class T1, class T2, class T3, class T4, class FUN>
54 void mergeStreamGridGrid(AMI_STREAM<T1> *grid1,
55 AMI_STREAM<T2> *grid2,
56 dimension_type rows, dimension_type cols,
57 AMI_STREAM<T3> *str,
58 FUN fo,
59 AMI_STREAM<T4> *outStream);
60 void merge2waterBase(AMI_STREAM<waterType> *unsortedWaterStr,
61 AMI_STREAM<direction_type> *dirStr,
62 AMI_STREAM<elevation_type> *elStr,
63 AMI_STREAM<waterWindowBaseType> *merge);
64 AMI_STREAM<waterGridType>*
65 merge2waterGrid(AMI_STREAM<waterType> *unsortedWaterStr,
66 AMI_STREAM<direction_type> *dirStr,
67 AMI_STREAM<elevation_type> *elStr);
68 AMI_STREAM<boundaryType> *
69 findBoundariesMain(AMI_STREAM<labelElevType> *labeledWater);
70
71
72
73
74 /* ********************************************************************** */
75 /* some helper classes */
76 /* ********************************************************************** */
77
78 class printElevation {
79 public:
operator ()(const elevation_type & p)80 char *operator()(const elevation_type &p) {
81 static char buf[20];
82 sprintf(buf, "%.1f", (float)p);
83 return buf;
84 }
85 };
86
87 class printDirection {
88 public:
operator ()(const direction_type & p)89 char *operator()(const direction_type &p) {
90 static char buf[20];
91 sprintf(buf, "%3d", p);
92 return buf;
93 }
operator ()(const waterWindowBaseType & p)94 char *operator()(const waterWindowBaseType &p) {
95 static char buf[20];
96 sprintf(buf, "%3d", p.dir);
97 return buf;
98 }
99 #if(0)
operator ()(const waterWindowBaseType & p)100 char *operator()(const waterWindowBaseType &p) {
101 static char buf[3];
102 buf[0] = directionSymbol(p.dir);
103 buf[1] = '\0';
104 return buf;
105 }
106 #endif
107 };
108
109 class printLabel {
110 public:
operator ()(const labelElevType & p)111 char *operator()(const labelElevType&p) {
112 static char buf[8];
113 sprintf(buf, CCLABEL_FMT, p.getLabel());
114 return buf;
115 }
operator ()(const waterGridType & p)116 char *operator()(const waterGridType &p) {
117 static char buf[8];
118 sprintf(buf, CCLABEL_FMT, p.getLabel());
119 return buf;
120 }
operator ()(const waterType & p)121 char *operator()(const waterType &p) {
122 static char buf[8];
123 sprintf(buf, CCLABEL_FMT, p.getLabel());
124 return buf;
125 }
126 };
127
128 class printDepth {
129 public:
operator ()(const waterGridType & p)130 char *operator()(const waterGridType &p) {
131 static char buf[3];
132 sprintf(buf, "%1u", p.depth);
133 return buf;
134 }
135 };
136
137
138
139
140 char *
verbosedir(std::string s)141 verbosedir(std::string s) {
142 static char buf[BUFSIZ];
143 sprintf(buf, "dump/%s", s.c_str());
144 return buf;
145 }
146
147
148
149 /* ---------------------------------------------------------------------- */
150 /* fill the terrain (if necessary) and compute flow direction stream;
151 elstr must exist and contain grid data before call, filledstr and
152 dirstr are created; elstr is deleted and replaced with the
153 classified elstr, which has boundary nodata distinguished from
154 inner nodata */
155 AMI_STREAM<waterWindowBaseType>*
computeFlowDirections(AMI_STREAM<elevation_type> * & elstr,AMI_STREAM<elevation_type> * & filledstr,AMI_STREAM<direction_type> * & dirstr,AMI_STREAM<labelElevType> * & labeledWater)156 computeFlowDirections(AMI_STREAM<elevation_type>*& elstr,
157 AMI_STREAM<elevation_type>*& filledstr,
158 AMI_STREAM<direction_type>*& dirstr,
159 AMI_STREAM<labelElevType> *& labeledWater) {
160
161 Rtimer rt, rtTotal;
162 AMI_STREAM<elevation_type> *elstr_reclass=NULL;
163 AMI_STREAM<ElevationWindow > *winstr=NULL;
164 AMI_STREAM<plateauStats> *statstr=NULL;
165 AMI_STREAM<plateauType> *platstr=NULL;
166 AMI_STREAM<waterType> *waterstr=NULL;
167 AMI_STREAM<waterGridType> *mergedWaterStr=NULL;
168 AMI_STREAM<boundaryType> *boundaryStr=NULL;
169 AMI_STREAM<waterWindowType> *waterWindows=NULL;
170
171 rt_start(rtTotal);
172 assert(elstr && filledstr == NULL && dirstr == NULL && labeledWater == NULL);
173 if (stats) {
174 stats->comment("------------------------------");
175 stats->comment("COMPUTING FLOW DIRECTIONS");
176
177 /* adjust nodata -- boundary nodata distinguished from inner
178 nodata */
179 stats->comment("classifying nodata (inner & boundary)");
180 }
181
182 elstr_reclass = classifyNodata(elstr);
183 delete elstr;
184 elstr = elstr_reclass;
185
186
187 /* ---------------------------------------------------------------------- */
188 /* find the plateaus. */
189 /* ---------------------------------------------------------------------- */
190 if (stats) {
191 stats->comment("----------", opt->verbose);
192 stats->comment("assigning preliminary directions");
193 }
194
195 rt_start(rt);
196 dirstr = new AMI_STREAM<direction_type>;
197 winstr = new AMI_STREAM<ElevationWindow>();
198 statstr = new AMI_STREAM<plateauStats>;
199
200 platstr = findPlateaus(elstr, nrows, ncols, nodataType::ELEVATION_NODATA,
201 winstr, dirstr, statstr);
202
203 delete winstr; /* not used; not made */
204 rt_stop(rt);
205
206 if (stats) {
207 stats->recordTime("findingPlateaus", rt);
208 stats->recordLength("plateaus", platstr);
209 stats->recordLength("plateau stats", statstr);
210 }
211 FILL_SAVEALL {
212 /* printStream(*stats, statstr); */
213 FILL_DEBUG cout << "sort plateauStr (by ij): ";
214 AMI_STREAM<plateauType> *tmp = sort(platstr, ijCmpPlateauType());
215 printStream2Grid(tmp, nrows, ncols,
216 verbosedir("label1.asc"), plateauType::printLabel);
217 delete tmp;
218 }
219
220
221 /* ---------------------------------------------------------------------- */
222 /* assign labels and directions & BFS ordering. depressions have
223 labels, but no direction information.
224 */
225 /* ---------------------------------------------------------------------- */
226 rt_start(rt);
227 waterstr = new AMI_STREAM<waterType>();
228 assignDirections(statstr, platstr, waterstr);
229 delete platstr;
230 delete statstr;
231 rt_stop(rt);
232 if (stats) {
233 stats->recordTime("assigning directions", rt);
234 *stats << "maxWatershedCount=" << labelFactory::getLabelCount() << endl;
235 }
236
237 rt_start(rt);
238 mergedWaterStr = merge2waterGrid(waterstr, dirstr, elstr);
239 delete dirstr;
240 delete waterstr;
241 rt_stop(rt);
242 if (stats) {
243 stats->recordTime("merging", rt);
244 stats->recordLength("mergedWaterStr", mergedWaterStr);
245 }
246
247 /* ---------------------------------------------------------------------- */
248 /* watershed analysis */
249 /* IN: mergedWaterStr, labelFactory::... */
250 /* ---------------------------------------------------------------------- */
251 if (stats) {
252 stats->comment("--------------", opt->verbose);
253 stats->comment("generating watersheds and watershed graph");
254 }
255
256 rt_start(rt);
257 waterWindows = new AMI_STREAM<waterWindowType>();
258 createWaterWindows(mergedWaterStr, nrows, ncols, waterWindows);
259 delete mergedWaterStr;
260 rt_stop(rt);
261 if (stats) {
262 stats->recordTime("creating windows", rt);
263 stats->recordLength("waterWindows", waterWindows);
264 }
265
266 /* ---------------------------------------------------------------------- */
267 rt_start(rt);
268 labeledWater = new AMI_STREAM<labelElevType>();
269 boundaryStr = new AMI_STREAM<boundaryType>();
270 generateWatersheds(&waterWindows, nrows, ncols, labeledWater, boundaryStr);
271
272 /* do we need to make boundaries here?? */
273 delete waterWindows;
274 /* cout << "bogus boundary length = " << boundaryStr->stream_len() << endl;*/
275 assert(boundaryStr->stream_len() == 0);
276 delete boundaryStr;
277
278 assert(labeledWater->stream_len() == nrows * ncols);
279 rt_stop(rt);
280 if (stats)
281 stats->recordTime("generating watersheds", rt);
282
283 /* ---------------------------------------------------------------------- */
284 /* find boundaries */
285 /* ---------------------------------------------------------------------- */
286 FILL_DEBUG cerr << "sort labeledWater (by ij):";
287 sort(&labeledWater, ijCmpLabelElevType());
288
289 #ifdef SAVE_ASCII
290 cerr << "saving WATERSHED grid as watershed_grid\n";
291 printStream2Grid(labeledWater, nrows, ncols,
292 "watershed.asc", labelElevType::printLabel);
293 #endif
294 boundaryStr = findBoundariesMain(labeledWater);
295
296
297 /* ---------------------------------------------------------------------- */
298 /* filling */
299 /* valid streams are: boundaryStr, labeledWater */
300 /* ---------------------------------------------------------------------- */
301 rt_start(rt);
302 elevation_type *raise;
303 /*find the raise elevations */
304
305 FILL_DEBUG cerr << "sort boundaryStr (by elev): ";
306 sort(&boundaryStr, elevCmpBoundaryType());
307
308 raise = fill_depression(boundaryStr, labelFactory::getLabelCount());
309 delete boundaryStr;
310 rt_stop(rt);
311 if (stats)
312 stats->recordTime("filling depressions", rt);
313
314 /*fill the terrain*/
315 rt_start(rt);
316 filledstr = new AMI_STREAM<elevation_type>();
317 commit_fill(labeledWater, raise, labelFactory::getLabelCount(), filledstr);
318 assert(filledstr->stream_len() == nrows * ncols);
319 delete [] raise;
320 rt_stop(rt);
321 if (stats) {
322 stats->recordTime("updating filled grid", rt);
323 stats->recordLength("filledstr", filledstr);
324
325 /* ---------------------------------------------------------------------- */
326 /* find plateaus again and reassign directions */
327 /* ---------------------------------------------------------------------- */
328
329 stats->comment("------------------------------");
330 stats->comment("REASSIGNING DIRECTIONS");
331 }
332
333 rt_start(rt);
334 winstr = NULL;
335 dirstr = new AMI_STREAM<direction_type>();
336 statstr = new AMI_STREAM<plateauStats>();
337 platstr = findPlateaus(filledstr, nrows, ncols, nodataType::ELEVATION_NODATA,
338 winstr, dirstr, statstr);
339 rt_stop(rt);
340 if (stats) {
341 stats->recordTime("findingPlateaus2", rt);
342 stats->recordLength("final plateaus", platstr);
343 stats->recordLength("final plateau stats", statstr);
344 }
345 FILL_SAVEALL {
346 FILL_DEBUG cout << "sort plateauStr (by ij): ";
347 AMI_STREAM<plateauType> *tmp = sort(platstr, ijCmpPlateauType());
348 printStream2Grid(tmp, nrows, ncols,
349 verbosedir("plateaus.asc"), plateauType::printLabel);
350 delete tmp;
351 }
352
353 /* assign final directions */
354 rt_start(rt);
355 waterstr = new AMI_STREAM<waterType>();
356 long dc = assignDirections(statstr, platstr, waterstr);
357 if(dc && stats) {
358 *stats << "WARNING: " << dc << " depressions (islands) detected\n";
359 }
360 delete platstr;
361 delete statstr;
362 rt_stop(rt);
363 if (stats)
364 stats->recordTime("final directions", rt);
365
366 /* merge */
367 rt_start(rt);
368 AMI_STREAM<waterWindowBaseType> *flowStream;
369 /*STREAM_TO_OPTION(flowStream, "flowStream");*/
370 char path[BUFSIZ];
371 char* base_dir = getenv(STREAM_TMPDIR);
372 assert(base_dir);
373 sprintf(path, "%s/flowStream", base_dir);
374 flowStream = new AMI_STREAM<waterWindowBaseType>(path);
375 /*flowStream->persist(PERSIST_PERSISTENT); */
376 if (stats)
377 stats->comment("creating flowStream: ", flowStream->sprint());
378
379 merge2waterBase(waterstr, dirstr, filledstr, flowStream);
380 delete waterstr;
381 rt_stop(rt);
382 if (stats)
383 stats->recordTime("merge water,dir,elev to flow", rt);
384 rt_stop(rtTotal);
385
386 #ifdef SAVE_ASCII
387 /*write grids as ascii file */
388 printGridStream(filledstr, nrows, ncols,
389 "filled_elev.asc", printElevation());
390 printGridStream(flowStream, nrows, ncols,
391 "direction.asc", printDirection());
392 #endif
393
394 if (stats) {
395 stats->recordTime("Total compute flow direction running time", rtTotal);
396 stats->comment("compute flow directions done.");
397 }
398
399 return flowStream;
400 }
401
402
403 /* ---------------------------------------------------------------------- */
404 void
recordWatersheds(AMI_STREAM<labelElevType> * labeledWater)405 recordWatersheds(AMI_STREAM<labelElevType> *labeledWater) {
406 AMI_err ae;
407 long labelCount = 0;
408 AMI_STREAM<labelElevType> *tmp;
409
410 if (stats)
411 *stats << "--- watershed stats" << endl;
412 FILL_DEBUG cout << "sort labeledWater (by wat label): ";
413 tmp = sort(labeledWater, labelCmpLabelElevType());
414
415 labelElevType *p;
416 cclabel_type prev(LABEL_UNDEF);
417 tmp->seek(0);
418 while((ae = tmp->read_item(&p)) == AMI_ERROR_NO_ERROR) {
419 if(p->getLabel() != prev) {
420 labelCount++;
421 prev = p->getLabel();
422 }
423 }
424 assert(ae == AMI_ERROR_END_OF_STREAM);
425
426 if (stats) {
427 *stats << "watershed count = " << labelCount << endl;
428 *stats << "---" << endl;
429 stats->flush();
430 }
431
432 delete tmp;
433 }
434
435
436
437
438
439 /* ********************************************************************** */
440 /* assign directions to plateaus that have sinks;
441 * reassign labels to depressions (don't drain out).
442 * all plateaus are written out to the waterstr. */
443 long
assignDirections(AMI_STREAM<plateauStats> * statstr,AMI_STREAM<plateauType> * platstr,AMI_STREAM<waterType> * waterstr)444 assignDirections(AMI_STREAM<plateauStats> *statstr,
445 AMI_STREAM<plateauType> *platstr,
446 AMI_STREAM<waterType> *waterstr) {
447 size_t fmem;
448 AMI_err ae;
449 plateauStats *ps;
450
451 if (stats) {
452 stats->comment("----------", opt->verbose);
453 stats->comment("assigning directions on plateaus");
454 }
455
456 labelFactory::reset(); /* we are relabeling now */
457
458 statstr->seek(0);
459 platstr->seek(0);
460 fmem = getAvailableMemory();
461 long depressionCount=0;
462 long spillCount=0;
463 while((ae = statstr->read_item(&ps)) == AMI_ERROR_NO_ERROR) {
464 if(ps->size*sizeof(gridElement) > fmem) {
465 cerr << "WARNING: grid larger than memory (ignored)" << endl;
466 }
467 assert(ps->label != LABEL_NODATA);
468 if(ps->hasSpill) {
469 spillCount++;
470 grid *platGrid = new grid(ps->iMin, ps->jMin, ps->iMax, ps->jMax,
471 ps->size, ps->label);
472 platGrid->load(*platstr);
473 platGrid->assignDirections(opt->d8 ? 1 : 0);
474 platGrid->save(*waterstr); /* this doesn't save labels */
475 delete platGrid;
476 } else {
477 /* depression - just give contiguous labels only */
478 depressionCount++;
479 cclabel_type label = labelFactory::getNewLabel();
480 for(int i=0; i<ps->size; i++) {
481 plateauType *pt;
482 platstr->read_item(&pt);
483 pt->cclabel = label; /* assign new label */
484 waterType wt(*pt); /* write it out */
485 ae = waterstr->write_item(wt);
486 assert(ae == AMI_ERROR_NO_ERROR);
487 }
488 }
489 }
490 if (stats) {
491 *stats << "depression count = " << depressionCount << endl;
492 *stats << "spill count = " << spillCount << endl;
493 }
494
495 return depressionCount;
496 }
497
498
499
500
501 /* ********************************************************************** */
502 /* assign directions to plateaus that have sinks;
503 * check that there are no depressions.
504 */
505 void
assignFinalDirections(AMI_STREAM<plateauStats> * statstr,AMI_STREAM<plateauType> * platstr,AMI_STREAM<waterType> * waterstr)506 assignFinalDirections(AMI_STREAM<plateauStats> *statstr,
507 AMI_STREAM<plateauType> *platstr,
508 AMI_STREAM<waterType> *waterstr) {
509 AMI_err ae;
510 plateauStats *ps;
511
512 if (stats)
513 stats->comment("assigning final directions");
514
515 statstr->seek(0);
516 platstr->seek(0);
517 while((ae = statstr->read_item(&ps)) == AMI_ERROR_NO_ERROR) {
518
519 if(ps->hasSpill) {
520 grid *platGrid = new grid(ps->iMin, ps->jMin, ps->iMax, ps->jMax,
521 ps->size, ps->label);
522 platGrid->load(*platstr);
523 platGrid->assignDirections(opt->d8 ? 1 : 0);
524 platGrid->save(*waterstr); /* this doesn't save labels */
525 delete platGrid;
526 } else {
527 /* could be legitimate */
528 cerr << "WARNING: depression detected: " << *ps << endl;
529 continue;
530 }
531 }
532 };
533
534
535
536 /* ********************************************************************** */
537 class directionElevationMerger {
538 public:
operator ()(elevation_type el,direction_type dir,waterType p)539 waterGridType operator()(elevation_type el, direction_type dir,
540 waterType p) {
541 /* check that no (boundary) nodata values got in here */
542 assert(el != nodataType::ELEVATION_BOUNDARY);
543 assert(!is_nodata(el)); /* p should be a valid grid cell */
544 return waterGridType(el, p.getDirection(), p.getLabel(), p.getDepth());
545 }
operator ()(elevation_type el,direction_type dir)546 waterGridType operator()(elevation_type el, direction_type dir) {
547 waterGridType wg(el, dir);
548 if(el == nodataType::ELEVATION_BOUNDARY) { /* hack XXX (approved) */
549 wg.setLabel(LABEL_BOUNDARY);
550 }
551 /* nodata => boundary or undef */
552 assert(!is_nodata(el) ||
553 (wg.getLabel()==LABEL_BOUNDARY || wg.getLabel()==LABEL_UNDEF));
554 return wg;
555 }
556 };
557
558
559
560 /* ---------------------------------------------------------------------- */
561 AMI_STREAM<waterGridType> *
merge2waterGrid(AMI_STREAM<waterType> * unsortedWaterStr,AMI_STREAM<direction_type> * dirStr,AMI_STREAM<elevation_type> * elStr)562 merge2waterGrid(AMI_STREAM<waterType> *unsortedWaterStr,
563 AMI_STREAM<direction_type> *dirStr,
564 AMI_STREAM<elevation_type> *elStr) {
565 AMI_STREAM<waterType> *waterStr;
566
567
568 FILL_DEBUG cout << "sort waterStr (by ij): ";
569 waterStr = sort(unsortedWaterStr, ijCmpWaterType());
570
571 FILL_SAVEALL printStream2Grid(waterStr, nrows, ncols,
572 verbosedir("platlabels.asc"), printLabel());
573
574 AMI_STREAM<waterGridType> *mergedWaterStr = new AMI_STREAM<waterGridType>();
575 mergeStreamGridGrid(elStr, dirStr,
576 nrows, ncols,
577 waterStr,
578 directionElevationMerger(),
579 mergedWaterStr);
580 delete waterStr;
581 FILL_SAVEALL printGridStream(mergedWaterStr, nrows, ncols,
582 verbosedir("mergedlabels.asc"), printLabel());
583
584 assert(mergedWaterStr->stream_len());
585 return mergedWaterStr;
586 }
587
588
589
590 /* ---------------------------------------------------------------------- */
591 void
merge2waterBase(AMI_STREAM<waterType> * unsortedWaterStr,AMI_STREAM<direction_type> * dirStr,AMI_STREAM<elevation_type> * elStr,AMI_STREAM<waterWindowBaseType> * merge)592 merge2waterBase(AMI_STREAM<waterType> *unsortedWaterStr,
593 AMI_STREAM<direction_type> *dirStr,
594 AMI_STREAM<elevation_type> *elStr,
595 AMI_STREAM<waterWindowBaseType> *merge) {
596 AMI_STREAM<waterType> *waterStr;
597 FILL_DEBUG cout << "sort waterStr (by ij): ";
598 waterStr = sort(unsortedWaterStr, ijCmpWaterType());
599 mergeStreamGridGrid(elStr, dirStr,
600 nrows, ncols,
601 waterStr,
602 directionElevationMerger(),
603 merge);
604 delete waterStr;
605 }
606
607
608
609
610 /* ---------------------------------------------------------------------- */
611 /*
612 * merge 2 grids and a stream together to form a new grid.
613 * str should be sorted in ij order
614 */
615 template<class T1, class T2, class T3, class T4, class FUN>
616 void
mergeStreamGridGrid(AMI_STREAM<T1> * grid1,AMI_STREAM<T2> * grid2,dimension_type rows,dimension_type cols,AMI_STREAM<T3> * str,FUN fo,AMI_STREAM<T4> * outStream)617 mergeStreamGridGrid(AMI_STREAM<T1> *grid1,
618 AMI_STREAM<T2> *grid2,
619 dimension_type rows, dimension_type cols,
620 AMI_STREAM<T3> *str,
621 FUN fo,
622 AMI_STREAM<T4> *outStream) {
623 T1 *t1p;
624 T2 *t2p;
625 T3 *t3p;
626 AMI_err aeUpd, ae;
627
628
629 grid1->seek(0);
630 grid2->seek(0);
631 str->seek(0);
632 aeUpd = str->read_item(&t3p);
633 assert(aeUpd == AMI_ERROR_NO_ERROR || aeUpd == AMI_ERROR_END_OF_STREAM);
634
635 for(dimension_type row = 0; row < rows; row++) {
636 for(dimension_type col = 0; col < cols; col++) {
637 ae = grid1->read_item(&t1p);
638 assert(ae == AMI_ERROR_NO_ERROR);
639 ae = grid2->read_item(&t2p);
640 assert(ae == AMI_ERROR_NO_ERROR);
641
642 T4 t4;
643 if(aeUpd == AMI_ERROR_NO_ERROR && t3p->i == row && t3p->j == col) {
644 /* cell present in stream */
645 t4 = fo(*t1p, *t2p, *t3p);
646 aeUpd = str->read_item(&t3p);
647 assert(aeUpd == AMI_ERROR_NO_ERROR ||
648 aeUpd == AMI_ERROR_END_OF_STREAM);
649 } else {
650 /* not in stream */
651 t4 = fo(*t1p, *t2p);
652 }
653 ae = outStream->write_item(t4);
654 assert(ae == AMI_ERROR_NO_ERROR);
655 }
656 /*assert(outStream->stream_len() == (row+1) * cols); */
657 }
658 assert(outStream->stream_len() == rows * cols);
659 return;
660 }
661
662
663
664 /* ---------------------------------------------------------------------- */
665 /* make boundaryStr from labeledWater */
666 AMI_STREAM<boundaryType> *
findBoundariesMain(AMI_STREAM<labelElevType> * labeledWater)667 findBoundariesMain(AMI_STREAM<labelElevType> *labeledWater) {
668 AMI_STREAM<boundaryType> *boundaryStr;
669 Rtimer rt;
670
671 rt_start(rt);
672 boundaryStr = new AMI_STREAM<boundaryType>();
673 FILL_SAVEALL printGridStream(labeledWater, nrows, ncols,
674 verbosedir("labels.asc"), printLabel());
675
676 findBoundaries(labeledWater, nrows, ncols, boundaryStr);
677 if (stats)
678 stats->recordLength("all boundaries", boundaryStr);
679
680 FILL_SAVEALL {
681 FILL_DEBUG cout << "sort boundaryStr (by ij): ";
682 sort(&boundaryStr, ijCmpBoundaryType());
683 removeDuplicatesEx(&boundaryStr, ijCmpBoundaryType());
684 printStream2Grid(boundaryStr, nrows, ncols,
685 verbosedir("boundary.asc"), boundaryType::print);
686 }
687 FILL_DEBUG cout << "sort boundaryStr (by wat label): ";
688 sort(&boundaryStr, waterCmpBoundaryType());
689 removeDuplicatesEx(&boundaryStr, boundaryCmpBoundaryType());
690
691 rt_stop(rt);
692 if (stats) {
693 stats->recordTime("generating boundaries", rt);
694 stats->recordLength("boundary stream", boundaryStr);
695 }
696
697 /*if(GETOPT("veryfillverbose")) printStream(cout, boundaryStr);
698 SAVE_ON_OPTION(boundaryStr, "saveBoundaryStream");
699 */
700 return boundaryStr;
701 }
702
703
704