← Back to branch summary

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

« back to all changes in this revision

Viewing changes to raster/r.flow/main.c

  • 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.flow/main.c
 
1
/*
 
2
 **  Original Algorithm:    H. Mitasova, L. Mitas, J. Hofierka, M. Zlocha 
 
3
 **  New GRASS Implementation:  J. Caplan, M. Ruesink  1995
 
4
 **
 
5
 **  US Army Construction Engineering Research Lab, University of Illinois 
 
6
 **
 
7
 **  Copyright  J. Caplan, H. Mitasova, L. Mitas, M.Ruesink, J. Hofierka, 
 
8
 **     M. Zlocha  1995
 
9
 **
 
10
 **This program is free software; you can redistribute it and/or
 
11
 **modify it under the terms of the GNU General Public License
 
12
 **as published by the Free Software Foundation; either version 2
 
13
 **of the License, or (at your option) any later version.
 
14
 **
 
15
 **This program is distributed in the hope that it will be useful,
 
16
 **but WITHOUT ANY WARRANTY; without even the implied warranty of
 
17
 **MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
18
 **GNU General Public License for more details.
 
19
 **
 
20
 **You should have received a copy of the GNU General Public License
 
21
 **along with this program; if not, write to the Free Software
 
22
 **Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
 
23
 **
 
24
 **  version 13 for GRASS5.0
 
25
 **  FP related bugs in slope length output fixed by Helena oct. 1999)
 
26
 **  Update MN: commented line 387
 
27
 */
 
28
 
 
29
#include <grass/gis.h>
 
30
#include <grass/raster.h>
 
31
#include <grass/glocale.h>
 
32
#include "r.flow.h"
 
33
#include "mem.h"
 
34
#include "io.h"
 
35
#include "aspect.h"
 
36
#include "precomp.h"
 
37
 
 
38
#define HORIZ   1               /* \            */
 
39
#define VERT    0               /* |            */
 
40
#define EAST    1               /* |            */
 
41
#define WEST    0               /* |_ magic     */
 
42
#define NORTH   1               /* |  numbers   */
 
43
#define SOUTH   0               /* |            */
 
44
#define ROW     1               /* |            */
 
45
#define COL     0               /* /            */
 
46
 
 
47
CELL v;                         /* address for segment retrieval macros */
 
48
 
 
49
/* heap memory */
 
50
struct Cell_head region;        /* resolution and boundaries            */
 
51
struct Map_info fl;             /* output vector file header            */
 
52
struct BM *bitbar;              /* space-efficient barrier matrix       */
 
53
int lgfd;                       /* output length file descriptor        */
 
54
char string[1024];              /* space for strings                    */
 
55
layer el, as, ds;               /* elevation, aspect, density (accumulation) */
 
56
double *ew_dist;                /* east-west distances for rows         */
 
57
double *epsilon[2];             /* quantization errors for rows         */
 
58
 
 
59
/* command-line parameters */
 
60
params parm;
 
61
 
 
62
typedef struct
 
63
{
 
64
    int row, col;               /* current matrix address       */
 
65
}
 
66
addr;
 
67
 
 
68
typedef int bbox[2][2];         /* current bounding box         */
 
69
 
 
70
typedef struct
 
71
{
 
72
    double x, y, z;             /* exact earth coordinates      */
 
73
    double theta;               /* aspect                       */
 
74
    double r, c;                /* cell matrix coordinates      */
 
75
}
 
76
point;
 
77
 
 
78
typedef struct
 
79
{
 
80
    double *px, *py;            /* x/y point arrays             */
 
81
    int index;                  /* index of next point          */
 
82
}
 
83
flowline;
 
84
 
 
85
/***************************** CALCULATION ******************************/
 
86
 
 
87
/*
 
88
 * height_angle_bounding_box: averages matrix values at sub between
 
89
 *      floor(cut) and ceil(cut), based on proximity; adjusts bbox
 
90
 * globals r: o, z, parm
 
91
 * params  w: p, b
 
92
 */
 
93
static void
 
94
height_angle_bounding_box(int sub, double cut, int horiz, point * p, bbox b)
 
95
{
 
96
    int c, f = (int)cut;
 
97
    double r = cut - (double)f;
 
98
    double a1, a2, a, d;
 
99
 
 
100
    b[horiz][horiz] = sub - 1;
 
101
    b[horiz][!horiz] = sub + 1;
 
102
    b[!horiz][horiz] = f + 1;
 
103
    c = (b[!horiz][!horiz] = f - (r == 0.)) + 1;
 
104
 
 
105
    if (horiz) {
 
106
        a1 = (double)aspect(sub, f);
 
107
        a2 = (double)aspect(sub, c);
 
108
        p->z = (double)get(el, sub, f) * (1. - r) +
 
109
            (double)get(el, sub, c) * r;
 
110
    }
 
111
    else {
 
112
        a1 = (double)aspect(f, sub);
 
113
        a2 = (double)aspect(c, sub);
 
114
        p->z = (double)get(el, f, sub) * (1. - r) +
 
115
            (double)get(el, c, sub) * r;
 
116
    }
 
117
 
 
118
    if (!(a1 == UNDEF || a2 == UNDEF) &&
 
119
        !(Rast_is_d_null_value(&a1) || Rast_is_d_null_value(&a2)))
 
120
        /*    if (!(Rast_is_d_null_value(&a1) || Rast_is_d_null_value(&a2))) */
 
121
    {
 
122
        if ((d = a1 - a2) >= D_PI || d <= -D_PI) {
 
123
            if (a2 > D_PI)
 
124
                a2 -= D2_PI;
 
125
            else
 
126
                a1 -= D2_PI;
 
127
        }
 
128
        a = r * a2 + (1. - r) * a1;
 
129
        p->theta = a + (a < 0.) * D2_PI;
 
130
    }
 
131
    else
 
132
        p->theta = UNDEF;
 
133
 
 
134
    return;
 
135
}
 
136
 
 
137
/*
 
138
 * sloping: returns elevation condition for continuing current line
 
139
 */
 
140
#define sloping(z1, z2) (z1 > z2)
 
141
 
 
142
/*
 
143
 * on_map: returns map boundary condition for continuing current line
 
144
 * globals r:  region
 
145
 */
 
146
static int on_map(int sub, double cut, int horiz)
 
147
{
 
148
    return
 
149
        (sub >= 0 && cut >= 0.0 &&
 
150
         ((horiz && sub < region.rows && cut <= (double)(region.cols - 1)) ||
 
151
          (!horiz && sub < region.cols && cut <= (double)(region.rows - 1))));
 
152
}
 
153
 
 
154
/*
 
155
 * add_to_line: puts a new point on the end of the current flowline
 
156
 * globals r:  parm
 
157
 * params  w:  f
 
158
 */
 
159
static void add_to_line(point * p, flowline * f)
 
160
{
 
161
    if (parm.flout) {
 
162
        f->px[f->index] = (double)p->x;
 
163
        f->py[f->index] = (double)p->y;
 
164
    }
 
165
    ++f->index;
 
166
}
 
167
 
 
168
/*
 
169
 * rectify: correct quantization problems (designed for speed, not elegance)
 
170
 */
 
171
static double rectify(double delta, double bd[2], double e)
 
172
{
 
173
    if (delta > 0.) {
 
174
        if (delta > bd[1] + e)
 
175
            return delta;
 
176
    }
 
177
    else {
 
178
        if (delta < bd[0] - e)
 
179
            return delta;
 
180
    }
 
181
    if (delta < bd[1] - e)
 
182
        if (delta > bd[0] + e)
 
183
            return delta;
 
184
        else
 
185
            return bd[0];
 
186
    else
 
187
        return bd[1];
 
188
}
 
189
 
 
190
/*
 
191
 * next_point: computes next point based on current point, z, and o
 
192
 *      returns continuation condition
 
193
 * globals r:  region, bitbar, parm
 
194
 * params  w:  p, a, l
 
195
 * globals w:  density
 
196
 */
 
197
static int next_point(point * p,        /* current/next point               */
 
198
                      addr * a, /* current/next matrix address      */
 
199
                      bbox b,   /* current/next bounding box        */
 
200
                      double *l /* current/eventual length          */
 
201
    )
 
202
{
 
203
    int sub;
 
204
    double cut;
 
205
    int horiz;
 
206
    int semi;
 
207
    double length, delta;
 
208
    double deltaz;
 
209
    double tangent;
 
210
 
 
211
    double oldz = p->z;
 
212
    double oldtheta = p->theta;
 
213
    double oldr = p->r;
 
214
    double oldc = p->c;
 
215
 
 
216
    addr ads;
 
217
    double bdy[2], bdx[2];
 
218
 
 
219
    ads = *a;
 
220
    bdy[SOUTH] = (double)(oldr - b[ROW][SOUTH]) * region.ns_res;
 
221
    bdy[NORTH] = (double)(oldr - b[ROW][NORTH]) * region.ns_res;
 
222
    bdx[WEST] = (double)(b[COL][WEST] - oldc) * ew_dist[ads.row];
 
223
    bdx[EAST] = (double)(b[COL][EAST] - oldc) * ew_dist[ads.row];
 
224
 
 
225
    semi = oldtheta < 90 || oldtheta >= 270;
 
226
    tangent = tan(oldtheta * DEG2RAD);
 
227
 
 
228
    if (oldtheta != 90 && oldtheta != 270 &&    /* north/south */
 
229
        (delta = (bdy[semi] * tangent)) < bdx[EAST] && delta > bdx[WEST]) {
 
230
        delta = rectify(delta, bdx, epsilon[HORIZ][ads.row]);
 
231
        p->x += delta;
 
232
        p->y += bdy[semi];
 
233
        p->r = (double)b[ROW][semi];
 
234
        p->c += delta / ew_dist[ads.row];
 
235
        a->row = b[ROW][semi];
 
236
        a->col = ROUND(p->c);
 
237
        sub = b[ROW][semi];
 
238
        cut = p->c;
 
239
        horiz = HORIZ;
 
240
        if (parm.lgout)
 
241
            length = hypot(delta, bdy[semi]);
 
242
    }
 
243
    else {                      /*  east/west  */
 
244
 
 
245
        semi = oldtheta < 180;
 
246
        if (oldtheta == 90 || oldtheta == 270)
 
247
            delta = 0;
 
248
        else {
 
249
            /* I don't know if this is right case.
 
250
             * Anyway, should be avoid from dividing by zero.
 
251
             * Any hydrologic idea?
 
252
             */
 
253
            if (tangent == 0.0)
 
254
                tangent = 0.000001;
 
255
 
 
256
            delta = bdx[semi] / tangent;
 
257
        }
 
258
 
 
259
        delta = rectify(delta, bdy, epsilon[VERT][ads.row]);
 
260
        p->y += delta;
 
261
        p->x += bdx[semi];
 
262
        p->r -= delta / region.ns_res;
 
263
        p->c = (double)b[COL][semi];
 
264
        a->row = ROUND(p->r);
 
265
        a->col = b[COL][semi];
 
266
        sub = b[COL][semi];
 
267
        cut = p->r;
 
268
        horiz = VERT;
 
269
        if (parm.lgout)
 
270
            length = hypot(bdx[semi], delta);
 
271
    }
 
272
 
 
273
    if (on_map(sub, cut, horiz) &&
 
274
        (height_angle_bounding_box(sub, cut, horiz, p, b),
 
275
         sloping(oldz, p->z)) &&
 
276
        !(parm.barin && BM_get(bitbar, a->col, a->row))) {
 
277
        if (parm.dsout && (ads.row != a->row || ads.col != a->col))
 
278
            put(ds, a->row, a->col, get(ds, a->row, a->col) + 1);
 
279
        /*      if (parm.lgout) 
 
280
         *      *l += parm.l3d ? hypot(length, oldz - p->z) : length; this did not work, helena*/
 
281
        if (parm.lgout) {
 
282
            if (parm.l3d) {
 
283
                deltaz = oldz - p->z;   /*fix by helena Dec. 06 */
 
284
                *l += hypot(length, deltaz);
 
285
            }
 
286
            else
 
287
                *l += length;
 
288
        }
 
289
        return 1;
 
290
    }
 
291
 
 
292
    return 0;
 
293
}
 
294
 
 
295
/*
 
296
 * calculate: create a flowline for each cell
 
297
 * globals r: region, bitbar, parm, lgfd
 
298
 */
 
299
static void calculate(void)
 
300
{
 
301
    point pts;
 
302
    addr ads;
 
303
    bbox bbs;
 
304
    flowline fls;
 
305
    int row, col;
 
306
 
 
307
    /*    double     x, y, length, xstep, ystep, roffset, coffset; */
 
308
    double x, y, length, xstep, ystep;
 
309
    FCELL *lg = Rast_allocate_f_buf();
 
310
    struct line_pnts *points = Vect_new_line_struct();
 
311
    struct line_cats *cats = Vect_new_cats_struct();
 
312
    int loopstep = (!parm.dsout && !parm.lgout && parm.flout) ? parm.skip : 1;
 
313
 
 
314
    G_important_message(_("Calculating..."));
 
315
 
 
316
    fls.px = (double *)G_calloc(parm.bound, sizeof(double));
 
317
    fls.py = (double *)G_calloc(parm.bound, sizeof(double));
 
318
 
 
319
    ystep = region.ns_res * (double)loopstep;
 
320
 
 
321
    /* FIXME - allow seed to be specified for repeatability */
 
322
    G_srand48_auto();
 
323
 
 
324
    for (row = 0, y = (double)region.north - (region.ns_res * .5);
 
325
         row < region.rows; row += loopstep, y -= ystep) {
 
326
        xstep = ew_dist[row] * (double)loopstep;
 
327
        G_percent(row, region.rows, 2);
 
328
 
 
329
        for (col = 0, x = (double)region.west + (ew_dist[row] * .5);
 
330
             col < region.cols; col += loopstep, x += xstep) {
 
331
 
 
332
            length = 0.0;
 
333
            fls.index = 0;
 
334
 
 
335
            if (!(parm.barin && BM_get(bitbar, col, row))) {
 
336
#ifdef OFFSET
 
337
                /* disabled by helena June 2005 */
 
338
                roffset = parm.offset * (double)region.ew_res
 
339
                    * ((2. * G_drand48()) - 1.);
 
340
                coffset = parm.offset * (double)region.ns_res
 
341
                    * ((2. * G_drand48()) - 1.);
 
342
#endif
 
343
                pts.x = x;
 
344
                pts.y = y;
 
345
                pts.z = (double)get(el, row, col);
 
346
                pts.theta = (double)aspect(row, col);
 
347
                pts.r = (double)row;    /* + roffset; */
 
348
                pts.c = (double)col;    /*+ coffset; */
 
349
 
 
350
                ads.row = row;
 
351
                ads.col = col;
 
352
 
 
353
#ifdef OFFSET
 
354
                G_debug(3, "dx: %f  x: %f %f  row: %f %f\n",
 
355
                        roffset, x, pts.x, (double)row, pts.r);
 
356
                G_debug(3, "dy: %f  y: %f %f  col: %f %f\n",
 
357
                        roffset, y, pts.y, (double)col, pts.c);
 
358
#endif
 
359
 
 
360
                bbs[ROW][SOUTH] = row + 1;
 
361
                bbs[ROW][NORTH] = row - 1;
 
362
                bbs[COL][WEST] = col - 1;
 
363
                bbs[COL][EAST] = col + 1;
 
364
 
 
365
                do
 
366
                    add_to_line(&pts, &fls);
 
367
                while (fls.index <= parm.bound &&
 
368
                       (pts.z != UNDEFZ && pts.theta >= 0 && pts.theta <= 360)
 
369
                       &&
 
370
                       /*  (!Rast_is_d_null_value(&pts.z) && pts.theta != UNDEF) && */
 
371
                       next_point(&pts, &ads, bbs, &length));
 
372
            }
 
373
 
 
374
            if (fls.index > 1 && parm.flout &&
 
375
                (loopstep == parm.skip ||
 
376
                 !(row % parm.skip || col % parm.skip))) {
 
377
                Vect_copy_xyz_to_pnts(points, fls.px, fls.py, NULL,
 
378
                                      fls.index);
 
379
                Vect_write_line(&fl, GV_LINE, points, cats);
 
380
            }
 
381
 
 
382
            if (parm.lgout)
 
383
                lg[col] = (float)length;
 
384
        }
 
385
 
 
386
        if (parm.lgout)
 
387
            Rast_put_f_row(lgfd, lg);
 
388
    }
 
389
    G_percent (1, 1, 1);
 
390
    
 
391
    G_free(fls.px);
 
392
    G_free(fls.py);
 
393
    /*    G_free (fls); *//* commented 19/10/99 MN */
 
394
    G_free(lg);
 
395
    Vect_destroy_line_struct(points);
 
396
    Vect_destroy_cats_struct(cats);
 
397
 
 
398
    if (parm.lgout)
 
399
        Rast_close(lgfd);
 
400
}
 
401
 
 
402
int main(int argc, char *argv[])
 
403
{
 
404
    struct GModule *module;
 
405
    struct Option *pelevin, *paspin, *pbarin, *pskip, *pbound,
 
406
        *pflout, *plgout, *pdsout;
 
407
    struct Flag *fup, *flg, *fmem;
 
408
    int default_skip, larger, default_bound;
 
409
 
 
410
#ifdef OFFSET
 
411
    char *default_offset_ans, *offset_opt;
 
412
#endif
 
413
    char *default_skip_ans, *default_bound_ans, *skip_opt;
 
414
    struct History history;
 
415
 
 
416
 
 
417
    /* Initialize GIS engine */
 
418
    G_gisinit(argv[0]);
 
419
 
 
420
    module = G_define_module();
 
421
    G_add_keyword(_("raster"));
 
422
    G_add_keyword(_("hydrology"));
 
423
    module->label = _("Constructs flowlines.");
 
424
    module->description =
 
425
        _("Computes flowlines, flowpath lengths, "
 
426
          "and flowaccumulation (contributing areas) from a elevation raster "
 
427
          "map.");
 
428
 
 
429
    pelevin = G_define_standard_option(G_OPT_R_ELEV);
 
430
    
 
431
    paspin = G_define_standard_option(G_OPT_R_INPUT);
 
432
    paspin->key = "aspect";
 
433
    paspin->required = NO;
 
434
    paspin->description = _("Name of input aspect raster map");
 
435
    paspin->guisection = _("Input maps");
 
436
 
 
437
    pbarin = G_define_standard_option(G_OPT_R_INPUT);
 
438
    pbarin->key = "barrier";
 
439
    pbarin->required = NO;
 
440
    pbarin->description = _("Name of input barrier raster map");
 
441
    pbarin->guisection = _("Input maps");
 
442
 
 
443
    pskip = G_define_option();
 
444
    pskip->key = "skip";
 
445
    pskip->type = TYPE_INTEGER;
 
446
    pskip->required = NO;
 
447
    pskip->description = _("Number of cells between flowlines");
 
448
 
 
449
    pbound = G_define_option();
 
450
    pbound->key = "bound";
 
451
    pbound->type = TYPE_INTEGER;
 
452
    pbound->required = NO;
 
453
    pbound->description = _("Maximum number of segments per flowline");
 
454
 
 
455
    pflout = G_define_standard_option(G_OPT_V_OUTPUT);
 
456
    pflout->key = "flowline";
 
457
    pflout->required = NO;
 
458
    pflout->description = _("Name for output flow line vector map");
 
459
    pflout->guisection = _("Output maps");
 
460
 
 
461
    plgout = G_define_standard_option(G_OPT_R_OUTPUT);
 
462
    plgout->key = "flowlength";
 
463
    plgout->required = NO;
 
464
    plgout->description = _("Name for output flow path length raster map");
 
465
    plgout->guisection = _("Output maps");
 
466
 
 
467
    pdsout = G_define_standard_option(G_OPT_R_OUTPUT);
 
468
    pdsout->key = "flowaccumulation";
 
469
    pdsout->required = NO;
 
470
    pdsout->description = _("Name for output flow accumulation raster map");
 
471
    pdsout->guisection = _("Output maps");
 
472
 
 
473
    fup = G_define_flag();
 
474
    fup->key = 'u';
 
475
    fup->description =
 
476
        _("Compute upslope flowlines instead of default downhill flowlines");
 
477
 
 
478
    flg = G_define_flag();
 
479
    flg->key = '3';
 
480
    flg->description = _("3D lengths instead of 2D");
 
481
 
 
482
    fmem = G_define_flag();
 
483
    fmem->key = 'm';
 
484
    fmem->description = _("Use less memory, at a performance penalty");
 
485
 
 
486
    if (G_parser(argc, argv))
 
487
        exit(EXIT_FAILURE);
 
488
 
 
489
    G_get_set_window(&region);
 
490
 
 
491
    larger = ((region.cols < region.rows) ? region.rows : region.cols);
 
492
    default_skip = (larger < 50) ? 1 : (int)(larger / 50);
 
493
 
 
494
    default_skip_ans =
 
495
        G_calloc((int)log10((double)default_skip) + 2, sizeof(char));
 
496
    skip_opt = G_calloc((int)log10((double)larger) + 4, sizeof(char));
 
497
 
 
498
    sprintf(default_skip_ans, "%d", default_skip);
 
499
    sprintf(skip_opt, "1-%d", larger);
 
500
 
 
501
    default_bound = (int)(4. * hypot((double)region.rows,
 
502
                                     (double)region.cols));
 
503
    default_bound_ans =
 
504
        G_calloc((int)log10((double)default_bound) + 4, sizeof(char));
 
505
    sprintf(default_bound_ans, "0-%d", default_bound);
 
506
 
 
507
#ifdef OFFSET
 
508
    /* below fix changed from 0.0 to 1.0 and its effect disabled in 
 
509
     * calc.c, Helena June 2005 */
 
510
 
 
511
    default_offset = 1.0;       /* fixed 20. May 2001 Helena */
 
512
    default_offset_ans =
 
513
        G_calloc((int)log10(default_offset) + 2, sizeof(char));
 
514
    sprintf(default_offset_ans, "%f", default_offset);
 
515
 
 
516
    offset_opt = G_calloc((int)log10(default_offset) + 4, sizeof(char));
 
517
    sprintf(offset_opt, "0.0-500.0");
 
518
#endif
 
519
 
 
520
    if (!pskip->answer)
 
521
        pskip->answer = default_skip_ans;
 
522
 
 
523
    if (!pbound->answer)
 
524
        pbound->answer = default_bound_ans + 2;
 
525
 
 
526
    parm.elevin = pelevin->answer;
 
527
    parm.aspin = paspin->answer;
 
528
    parm.barin = pbarin->answer;
 
529
    parm.skip = atoi(pskip->answer);
 
530
    parm.bound = atoi(pbound->answer);
 
531
    parm.flout = pflout->answer;
 
532
    parm.lgout = plgout->answer;
 
533
    parm.dsout = pdsout->answer;
 
534
    parm.up = fup->answer;
 
535
    parm.l3d = flg->answer;
 
536
    parm.mem = fmem->answer;
 
537
 
 
538
    if (!pflout->answer && !plgout->answer && !pdsout->answer)
 
539
        G_fatal_error(_("You must select one or more output maps (flout, lgout, dsout)"));
 
540
 
 
541
    if (parm.seg)
 
542
        parm.mem = '\0';
 
543
    else if (parm.mem)
 
544
        parm.aspin = NULL;
 
545
 
 
546
    el.name = parm.elevin;
 
547
    as.name = (parm.aspin) ? parm.aspin : "internal aspects";
 
548
 
 
549
    ds.name = parm.dsout;
 
550
    el.row_offset = el.col_offset = 1;
 
551
    as.row_offset = as.col_offset = 0;
 
552
    ds.row_offset = ds.col_offset = 0;
 
553
 
 
554
    if ((G_projection() == PROJECTION_LL))      /* added MN 2005 */
 
555
        G_fatal_error(_("lat/long projection not supported by "
 
556
                        "r.flow. Please use 'r.watershed' for calculating "
 
557
                        "flow accumulation."));
 
558
 
 
559
    if (parm.flout || parm.dsout || parm.lgout) {
 
560
        open_output_files();
 
561
        allocate_heap();
 
562
        read_input_files();
 
563
 
 
564
        precompute();
 
565
        calculate();
 
566
        if (parm.dsout)
 
567
            write_density_file();
 
568
 
 
569
        close_files();
 
570
        deallocate_heap();
 
571
    }
 
572
 
 
573
    if (parm.dsout) {
 
574
        Rast_short_history(parm.dsout, "raster", &history);
 
575
        Rast_command_history(&history);
 
576
        Rast_write_history(parm.dsout, &history);
 
577
    }
 
578
    if (parm.lgout) {
 
579
        Rast_short_history(parm.lgout, "raster", &history);
 
580
        Rast_command_history(&history);
 
581
        Rast_write_history(parm.lgout, &history);
 
582
    }
 
583
 
 
584
    exit(EXIT_SUCCESS);
 
585
}