← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to raster/r.terraflow/water.cc

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
raster/r.terraflow/water.cc
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
 
 
20
 
#include <assert.h>
21
 
#include <iostream>
22
 
using namespace std;
23
 
 
24
 
#include <grass/iostream/ami.h>
25
 
 
26
 
 
27
 
#include "3scan.h"
28
 
#include "water.h"
29
 
#include "streamutils.h"
30
 
#include "sortutils.h"
31
 
 
32
 
 
33
 
#define WATER_DEBUG if(0)
34
 
#define XXX if(0)
35
 
 
36
 
 
37
 
char *
38
 
labelElevType::printLabel(const  labelElevType  &p) {
39
 
  static char buf[8];
40
 
  sprintf(buf, CCLABEL_FMT, p.label);
41
 
  return buf;
42
 
}
43
 
 
44
 
 
45
 
 
46
 
 
47
 
 
48
 
 
49
 
/* smaller elevation, depth is smaller priority */
50
 
int 
51
 
fillPriority::compare(const fillPriority &a, const fillPriority &b) {
52
 
  if(a.el < b.el) return -1;
53
 
  if(a.el > b.el) return 1;
54
 
  
55
 
  if(a.depth < b.depth) return -1;
56
 
  if(a.depth > b.depth) return 1;
57
 
  
58
 
  if(a.i < b.i) return -1;
59
 
  if(a.i > b.i) return 1;
60
 
  
61
 
  if(a.j < b.j) return -1;
62
 
  if(a.j > b.j) return 1;
63
 
  
64
 
  return 0;
65
 
}
66
 
 
67
 
int
68
 
fillPriority::qscompare(const void *a, const void *b) {
69
 
  fillPriority *x = (fillPriority*)a;
70
 
  fillPriority *y = (fillPriority*)b;
71
 
  return compare(*x, *y);
72
 
}
73
 
 
74
 
int 
75
 
operator<(const fillPriority &a, const fillPriority &b) {
76
 
  if(a.el < b.el) return 1;
77
 
  if(a.el > b.el) return 0;
78
 
  
79
 
  if(a.depth < b.depth) return 1;
80
 
  if(a.depth > b.depth) return 0;
81
 
  
82
 
  if(a.i < b.i) return 1;
83
 
  if(a.i > b.i) return 0;
84
 
  
85
 
  if(a.j < b.j) return 1;
86
 
  if(a.j > b.j) return 0;
87
 
  
88
 
  return 0;
89
 
}
90
 
 
91
 
int
92
 
operator<=(const fillPriority &a, const fillPriority &b) {
93
 
  if(a.el < b.el) return 1;
94
 
  if(a.el > b.el) return 0;
95
 
  
96
 
  if(a.depth < b.depth) return 1;
97
 
  if(a.depth > b.depth) return 0;
98
 
  
99
 
  if(a.i < b.i) return 1;
100
 
  if(a.i > b.i) return 0;
101
 
  
102
 
  if(a.j < b.j) return 1;
103
 
  if(a.j > b.j) return 0;
104
 
  
105
 
  return 1;
106
 
}
107
 
 
108
 
int 
109
 
operator>(const fillPriority &a, const fillPriority &b) {
110
 
  if(a.el < b.el) return 0;
111
 
  if(a.el > b.el) return 1;
112
 
  
113
 
  if(a.depth < b.depth) return 0;
114
 
  if(a.depth > b.depth) return 1;
115
 
  
116
 
  if(a.i < b.i) return 0;
117
 
  if(a.i > b.i) return 1;
118
 
  
119
 
  if(a.j < b.j) return 0;
120
 
  if(a.j > b.j) return 1;
121
 
  
122
 
  return 0;
123
 
}
124
 
 
125
 
 
126
 
int 
127
 
operator==(const fillPriority &a, const fillPriority &b) {
128
 
  return (a.el == b.el) 
129
 
        && (a.depth == b.depth) 
130
 
        && (a.i == b.i)
131
 
        && (a.j == b.j);
132
 
}
133
 
 
134
 
 
135
 
int 
136
 
operator!=(const fillPriority &a, const fillPriority &b) {
137
 
  return (a.el != b.el) 
138
 
        || (a.depth != b.depth) 
139
 
        || (a.i != b.i)
140
 
        || (a.j != b.j);        
141
 
}
142
 
 
143
 
 
144
 
ostream&
145
 
operator << (ostream& s, const fillPriority &p) {
146
 
  return s << "[fillPriority el=" << p.el 
147
 
                   << ", d=" << p.depth << ", "
148
 
                   << p.i  << ","
149
 
                   << p.j << "]";
150
 
}
151
 
 
152
 
 
153
 
 
154
 
/* ********************************************************************** */
155
 
 
156
 
 
157
 
ostream&
158
 
operator << (ostream& s, const labelElevType &p) {
159
 
  return s << (ijBaseType)p << " "
160
 
                   << "el=" << p.el << ", "
161
 
                   << p.label;
162
 
}
163
 
 
164
 
/* ********************************************************************** */
165
 
 
166
 
char *
167
 
waterType::printLabel(const waterType &p) {
168
 
  static char buf[8];
169
 
  sprintf(buf, CCLABEL_FMT, p.label);
170
 
  return buf;
171
 
}
172
 
 
173
 
/* ********************************************************************** */
174
 
 
175
 
 
176
 
char *
177
 
boundaryType::print(const boundaryType &p) {
178
 
  static char buf[4];
179
 
  if(p.isValid()) {
180
 
        buf[0] = '1';
181
 
  } else {
182
 
        buf[0] = '0';
183
 
  }
184
 
  buf[1] = '\0';
185
 
  
186
 
  return buf;
187
 
}
188
 
 
189
 
 
190
 
ostream&
191
 
operator << (ostream& s, const boundaryType &p) {
192
 
  return s << "[boundaryType "  
193
 
                   << (labelElevType)p << ", "
194
 
                   << p.label2 << "]";
195
 
}
196
 
 
197
 
 
198
 
 
199
 
/* ********************************************************************** */
200
 
/* ********************************************************************** */
201
 
 
202
 
class waterWindower {
203
 
private:
204
 
  AMI_STREAM<waterWindowType> *waterWindows;
205
 
public:
206
 
  waterWindower(AMI_STREAM<waterWindowType> *str) :
207
 
        waterWindows(str) {};
208
 
  void processWindow(dimension_type i, dimension_type j, 
209
 
                     waterGridType &point,
210
 
                     waterWindowBaseType *a,
211
 
                     waterWindowBaseType *b,
212
 
                     waterWindowBaseType *c);
213
 
};
214
 
 
215
 
void
216
 
waterWindower::processWindow(dimension_type i, dimension_type j, 
217
 
                             waterGridType &point,
218
 
                             waterWindowBaseType *a,
219
 
                             waterWindowBaseType *b,
220
 
                             waterWindowBaseType *c) {
221
 
  waterWindowType win = waterWindowType(i, j, point.getLabel(), a, b, c);
222
 
  AMI_err ae = waterWindows->write_item(win);
223
 
  assert(ae == AMI_ERROR_NO_ERROR);
224
 
}
225
 
 
226
 
void
227
 
createWaterWindows(AMI_STREAM<waterGridType> *mergedWaterStr, 
228
 
                   const dimension_type nrows, const dimension_type ncols,
229
 
                   AMI_STREAM<waterWindowType> *waterWindows) {
230
 
  stats->comment("creating windows", opt->verbose);
231
 
  waterWindower winfo(waterWindows);
232
 
  waterWindowBaseType nodata;
233
 
  assert(mergedWaterStr->stream_len() > 0);
234
 
  stats->comment("warning: using slower scan", opt->verbose);
235
 
  scan3(*mergedWaterStr, nrows, ncols, nodata, winfo);
236
 
}
237
 
 
238
 
 
239
 
/* ********************************************************************** */
240
 
/* ********************************************************************** */
241
 
 
242
 
 
243
 
/*
244
 
 * push labels to upslope neighbors 
245
 
 */
246
 
void
247
 
generateWatersheds(AMI_STREAM<waterWindowType> **waterWindows,
248
 
                                   const dimension_type nrows, const dimension_type ncols,
249
 
                                   AMI_STREAM<labelElevType> *labeledWater, 
250
 
                                   AMI_STREAM<boundaryType> *boundaryStr) {
251
 
  AMI_err ae;
252
 
  waterWindowType *winp, prevWin;
253
 
  assert(prevWin.getDepth() == DEPTH_INITIAL);
254
 
  EMPQueueAdaptive<fillPLabel, fillPriority> *pq;
255
 
 
256
 
  stats->comment("generateWatersheds", opt->verbose);
257
 
 
258
 
  assert((*waterWindows)->stream_len() == (nrows * ncols));
259
 
 
260
 
  WATER_DEBUG cout << "sort waterWindowsStream (by priority): ";
261
 
  sort(waterWindows, priorityCmpWaterWindowType());
262
 
  
263
 
  pq = new EMPQueueAdaptive<fillPLabel, fillPriority>();
264
 
  
265
 
/*   if(GETOPT("alwaysUseExternalPQ")) { */
266
 
/*      pq->makeExternal(); */
267
 
/*    } */
268
 
/*    if(GETOPT("useDebugPQ")) { */
269
 
/*      pq->makeExternalDebug(); */
270
 
/*    } */
271
 
  
272
 
  stats->comment("starting generate watersheds main loop", opt->verbose);
273
 
  
274
 
  assert((*waterWindows)->stream_len() == (nrows * ncols));
275
 
  /* not really in a grid, so row, col are not valid (but count correct) */
276
 
  for(dimension_type row=0; row<nrows; row++) {
277
 
    for(dimension_type col=0; col<ncols; col++) {
278
 
      ae = (*waterWindows)->read_item(&winp);
279
 
      assert(ae == AMI_ERROR_NO_ERROR);
280
 
      
281
 
      /* make sure it's sorted; prevWin default constructor should be ok */
282
 
      assert(winp->getPriority() > prevWin.getPriority());
283
 
      prevWin = *winp;
284
 
      
285
 
      XXX winp->sanityCheck();
286
 
      /* cout << "--- PROC: " << *winp << endl; */
287
 
      /* get my label(s) */
288
 
      fillPLabel plabel;                /* from the PQ */
289
 
      fillPriority prio;
290
 
      cclabel_type label = winp->getLabel();
291
 
#ifndef NDEBUG
292
 
      {
293
 
        /* check to see if things are coming out of the pq in
294
 
           order. just peek at the next one */
295
 
        fillPLabel tmp;
296
 
        XXX winp->sanityCheck();
297
 
        pq->min(tmp);
298
 
        /* XXX pq->verify(); */
299
 
        XXX winp->sanityCheck();
300
 
        assert(pq->is_empty() || winp->getPriority() <= tmp.getPriority());
301
 
      }
302
 
#endif
303
 
      while(pq->min(plabel) &&
304
 
                        ((prio=plabel.getPriority()) == winp->getPriority())) {
305
 
                /* XXX pq->verify(); */
306
 
                XXX winp->sanityCheck();
307
 
                pq->extract_min(plabel);
308
 
                /* XXX pq->verify(); */
309
 
                XXX winp->sanityCheck();
310
 
                if(label == LABEL_UNDEF) label = plabel.getLabel();
311
 
      } 
312
 
      /* no label! assign a new one */
313
 
      if((label==LABEL_UNDEF) && (!is_nodata(winp->getElevation()))) {
314
 
#ifndef NDEBUG
315
 
        {
316
 
          /* check to see if things are coming out of the pq in
317
 
             order. just peek at the next one */
318
 
          fillPLabel tmp;               
319
 
          XXX winp->sanityCheck();
320
 
          pq->min(tmp);
321
 
          /* XXX pq->verify(); */
322
 
          XXX winp->sanityCheck();
323
 
          assert(pq->is_empty() || winp->getPriority() <= tmp.getPriority());
324
 
        }
325
 
#endif
326
 
        if(IS_BOUNDARY(winp->i,winp->j,nrows, ncols)) {  /* edge of grid */
327
 
          assert(!is_nodata(winp->getElevation()));
328
 
          label = LABEL_BOUNDARY; /* reserved for watersheds draining
329
 
                                     out of grid */
330
 
        } else {
331
 
          label = labelFactory::getNewLabel();
332
 
        }
333
 
      }
334
 
      winp->setLabel(label);
335
 
      
336
 
      /* push label to 'upslope' neighbors.  let's assume that the
337
 
       * edges cause no problems, since they have no directions...  */
338
 
      if(label != LABEL_UNDEF) {
339
 
        int k=0;
340
 
        for(int i=-1; i<2; i++) {
341
 
          for(int j=-1; j<2; j++) {
342
 
            assert(k==linear(i,j));
343
 
            if(!is_nodata(winp->getElevation(k))
344
 
               && winp->drainsFrom(i,j)) { /* points to me */
345
 
              assert(i || j);
346
 
              prio = fillPriority(winp->getElevation(k),
347
 
                                  winp->getDepth(k),
348
 
                                  winp->i + i, winp->j + j);
349
 
#ifndef NDEBUG
350
 
              /* dont insert if preceeds us */
351
 
              if(winp->getPriority() < prio) {
352
 
                fillPLabel plabel(prio, label);
353
 
                pq->insert(plabel);
354
 
              } else {                  /* trying to send a label backward */
355
 
                cerr << "WARNING: time travel attempted" << endl; 
356
 
                cerr << "inst priority is " << prio << endl;
357
 
                cerr << "source is " << *winp << "; prio=" 
358
 
                     << winp->getPriority() << endl;
359
 
                assert(0);
360
 
              }
361
 
#else
362
 
              fillPLabel plabel(prio, label);
363
 
              pq->insert(plabel);
364
 
#endif
365
 
            }
366
 
            k++;
367
 
          }
368
 
        }
369
 
      }
370
 
      
371
 
      /* write myself to output */
372
 
      ae = labeledWater->write_item(winp->getCenter());
373
 
      assert(ae == AMI_ERROR_NO_ERROR);
374
 
    }
375
 
  }
376
 
  
377
 
  assert(pq->is_empty());
378
 
  delete pq;
379
 
  
380
 
  stats->comment("done with generate watersheds", opt->verbose);
381
 
}
382
 
 
383
 
 
384
 
 
385
 
/* ********************************************************************** */
386
 
 
387
 
 
388
 
class boundaryDetector {
389
 
private:
390
 
  const dimension_type nrows, ncols;
391
 
  AMI_STREAM<boundaryType> *boundaryStr;
392
 
  void processPair(labelElevType &pt, 
393
 
                   dimension_type i, dimension_type j, labelElevType &n);
394
 
public:
395
 
  boundaryDetector(AMI_STREAM<boundaryType> *str, 
396
 
                                   const dimension_type gnrows, const dimension_type gncols) 
397
 
    : nrows(gnrows), ncols(gncols), boundaryStr(str) {};
398
 
  
399
 
  void processWindow(dimension_type i, dimension_type j, 
400
 
                     labelElevType &point,
401
 
                     labelElevType *a,
402
 
                     labelElevType *b,
403
 
                     labelElevType *c);
404
 
};
405
 
 
406
 
template<class T>
407
 
T mymax(T a, T b) {
408
 
  return (a>b?a:b);
409
 
}
410
 
 
411
 
void
412
 
boundaryDetector::processPair(labelElevType &pt, 
413
 
                              dimension_type i, dimension_type j,
414
 
                              labelElevType &n) {
415
 
  if(n.getLabel() != LABEL_UNDEF && pt.getLabel() != n.getLabel()) {
416
 
    boundaryType bt(pt, mymax(pt.getElevation(), n.getElevation()), 
417
 
                    n.getLabel());
418
 
    AMI_err ae = boundaryStr->write_item(bt);
419
 
    assert(ae == AMI_ERROR_NO_ERROR);
420
 
  } else if(IS_BOUNDARY(i,j,nrows, ncols) && pt.getLabel() != LABEL_BOUNDARY) {
421
 
    /* this part makes sure that regions on the grid edge
422
 
       are considered 'boundary' */
423
 
    boundaryType bt(pt, LABEL_BOUNDARY);
424
 
    AMI_err ae = boundaryStr->write_item(bt);
425
 
    assert(ae == AMI_ERROR_NO_ERROR);
426
 
  }
427
 
}
428
 
 
429
 
void
430
 
boundaryDetector::processWindow(dimension_type i, dimension_type j, 
431
 
                                labelElevType &point,
432
 
                                labelElevType *a,
433
 
                                labelElevType *b,
434
 
                                labelElevType *c) {
435
 
  if(point.getLabel() == LABEL_UNDEF) return;
436
 
  /* NODATA_FIX */
437
 
  /* don't use the nodata as the boundary. */
438
 
  /* if(point.getLabel() == LABEL_NODATA) return; */
439
 
  assert(point.getLabel() != LABEL_NODATA);
440
 
  
441
 
  for(int k=0; k<3; k++) {
442
 
        processPair(point, i, j, a[k]);
443
 
        processPair(point, i, j, b[k]);
444
 
        processPair(point, i, j, c[k]);
445
 
  }
446
 
  /* processPair(point, i, j, b[0]); */
447
 
}
448
 
 
449
 
 
450
 
/* ********************************************************************** */
451
 
 
452
 
void
453
 
findBoundaries(AMI_STREAM<labelElevType> *labeledWater,
454
 
               const dimension_type nrows, const dimension_type ncols,   
455
 
               AMI_STREAM<boundaryType> *boundaryStr) {
456
 
  stats->comment("creating windows", opt->verbose);
457
 
  boundaryDetector det(boundaryStr, nrows, ncols);
458
 
  /* cerr << "WARNING: using scan3 instead of scan2" << endl; */
459
 
  scan3(*labeledWater, nrows, ncols, labelElevType(), det);
460
 
  
461
 
  /*  NODATA_FIX 
462
 
          assert(LABEL_BOUNDARY < LABEL_NODATA);
463
 
          boundaryType bt(-1, -1, ELEVATION_MIN, LABEL_BOUNDARY, LABEL_NODATA);
464
 
          AMI_err ae = boundaryStr->write_item(bt);
465
 
          assert(ae == AMI_ERROR_NO_ERROR);
466
 
  */
467
 
}
468
 
  
469
 
 
470
 
/* ********************************************************************** */
471
 
 
472
 
int
473
 
compressedWaterWindowBaseType::computeDelta(waterWindowBaseType *center,
474
 
                                            int index,
475
 
                                            waterWindowBaseType *p) const{
476
 
  if(center->el != p->el) {
477
 
    assert(p->depth == 1 || center->el > p->el);
478
 
    return 0;
479
 
  }
480
 
  if(index > 7) return 0;               /* we store our depth elsewhere */
481
 
  
482
 
  int d = p->depth - center->depth + 1;
483
 
  assert(d >= 0);
484
 
#ifndef NDEBUG
485
 
  if(d>2) {
486
 
        cerr << "whoops - assertion failure" << endl;
487
 
        cerr << "center = " << *center << endl;
488
 
        cerr << "p = " << *p << endl;
489
 
        cerr << "this = " << *this << endl;
490
 
  }
491
 
#endif
492
 
  assert(d <= 2);
493
 
  return d<<(2*index);
494
 
}
495
 
 
496
 
compressedWaterWindowBaseType::compressedWaterWindowBaseType(dimension_type gi,
497
 
                                                             dimension_type gj,
498
 
                                                             waterWindowBaseType *a, 
499
 
                                                             waterWindowBaseType *b, 
500
 
                                                             waterWindowBaseType *c) 
501
 
  : ijBaseType(gi, gj) {
502
 
  
503
 
  for(int i=0; i<3; i++) {
504
 
    el[i] = a[i].el;
505
 
    el[i+3] = b[i].el;
506
 
    el[i+6] = c[i].el;
507
 
  }
508
 
  
509
 
  for(int i=0; i<3; i++) {
510
 
    const direction_type mask_a[] = {2, 4, 8};
511
 
    const direction_type mask_b[] = {1, 0, 16};
512
 
    const direction_type mask_c[] = {128, 64, 32};
513
 
    points.setBit(i, a[i].dir & mask_a[i]);
514
 
    points.setBit(norm(i+3), b[i].dir & mask_b[i]);
515
 
    points.setBit(i+5, c[i].dir & mask_c[i]);
516
 
  }
517
 
  dir = b[1].dir;
518
 
  depth = b[1].depth;
519
 
  depth_delta = 0;
520
 
  
521
 
  /* nodata is not processed. */
522
 
  if(is_nodata(b[1].el)) {
523
 
    return;
524
 
  }
525
 
  
526
 
  for(int i=0; i<3; i++) {
527
 
    depth_delta |= computeDelta(b+1, norm(-1,i-1), a+i);
528
 
    depth_delta |= computeDelta(b+1, norm(0,i-1), b+i);
529
 
    depth_delta |= computeDelta(b+1, norm(1,i-1), c+i);
530
 
  }
531
 
}
532
 
 
533
 
fillPriority
534
 
compressedWaterWindowBaseType::getPriority() const {
535
 
  return fillPriority(getElevation(), getDepth(), i, j);
536
 
}
537
 
 
538
 
 
539
 
bfs_depth_type
540
 
compressedWaterWindowBaseType::getDepth(int k) const {
541
 
  if(getElevation() != getElevation(k)) return DEPTH_INITIAL;
542
 
  return depth + ((depth_delta >> (norm(k)*2)) & 0x3) - 1;
543
 
}
544
 
 
545
 
 
546
 
void 
547
 
compressedWaterWindowBaseType::sanityCheck() {
548
 
  assert(i >= -1);
549
 
  assert(j >= -1);
550
 
  assert(depth >= DEPTH_INITIAL);
551
 
}
552
 
  
553
 
 
554
 
ostream&
555
 
operator<<(ostream& s, const compressedWaterWindowBaseType &p) {
556
 
  return s << "[compressedWaterWindowBaseType " 
557
 
           << p.i << "," << p.j 
558
 
           << " " << directionSymbol(p.getDirection())
559
 
           << " e=" << p.getElevation() 
560
 
           << " d =" << p.getDepth() << "]";
561
 
}
562
 
 
563
 
/* ********************************************************************** */
564
 
 
565
 
labelElevType 
566
 
compressedWaterWindowType::getCenter() const {
567
 
  return labelElevType(i, j, getElevation(), label);
568
 
}
569
 
 
570
 
void 
571
 
compressedWaterWindowType::sanityCheck() {
572
 
  assert(label >= LABEL_UNDEF);
573
 
  compressedWaterWindowBaseType::sanityCheck();
574
 
}
575
 
 
576
 
 
577
 
ostream&
578
 
operator<<(ostream& s, const compressedWaterWindowType &p) {
579
 
  return s << "[compressedWaterWindowType " 
580
 
           << p.i << "," <<  p.j 
581
 
           << "  " << directionSymbol(p.getDirection())
582
 
           << " e=" << p.getElevation()
583
 
           << " d=" << p.getDepth() 
584
 
           << " l=" << p.label;
585
 
}
586
 
 
587