← Back to branch summary

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

« back to all changes in this revision

Viewing changes to raster/r.in.lidar/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.in.lidar/main.c
 
1
/****************************************************************************
 
2
 *
 
3
 * MODULE:       r.in.Lidar
 
4
 *               
 
5
 * AUTHOR(S):    Markus Metz
 
6
 *               Based on r.in.xyz by Hamish Bowman, Volker Wichmann
 
7
 *
 
8
 * PURPOSE:      Imports LAS LiDAR point clouds to a raster map using 
 
9
 *               aggregate statistics.
 
10
 *
 
11
 * COPYRIGHT:    (C) 2011 Markus Metz and the The GRASS Development Team
 
12
 *
 
13
 *               This program is free software under the GNU General Public
 
14
 *               License (>=v2). Read the file COPYING that comes with GRASS
 
15
 *               for details.
 
16
 *
 
17
 *****************************************************************************/
 
18
 
 
19
#include <stdio.h>
 
20
#include <stdlib.h>
 
21
#include <string.h>
 
22
#include <unistd.h>
 
23
#include <math.h>
 
24
#include <sys/types.h>
 
25
#include <grass/gis.h>
 
26
#include <grass/raster.h>
 
27
#include <grass/gprojects.h>
 
28
#include <grass/glocale.h>
 
29
#include <liblas/capi/liblas.h>
 
30
#include "local_proto.h"
 
31
 
 
32
struct node
 
33
{
 
34
    int next;
 
35
    double z;
 
36
};
 
37
 
 
38
#define SIZE_INCREMENT 10
 
39
int num_nodes = 0;
 
40
int max_nodes = 0;
 
41
struct node *nodes;
 
42
 
 
43
#define LAS_ALL 0
 
44
#define LAS_FIRST 1
 
45
#define LAS_LAST 2
 
46
#define LAS_MID 3
 
47
 
 
48
int new_node(void)
 
49
{
 
50
    int n = num_nodes++;
 
51
 
 
52
    if (num_nodes >= max_nodes) {
 
53
        max_nodes += SIZE_INCREMENT;
 
54
        nodes = G_realloc(nodes, (size_t)max_nodes * sizeof(struct node));
 
55
    }
 
56
 
 
57
    return n;
 
58
}
 
59
 
 
60
 
 
61
/* add node to sorted, single linked list 
 
62
 * returns id if head has to be saved to index array, otherwise -1 */
 
63
int add_node(int head, double z)
 
64
{
 
65
    int node_id, last_id, newnode_id, head_id;
 
66
 
 
67
    head_id = head;
 
68
    node_id = head_id;
 
69
    last_id = head_id;
 
70
 
 
71
    while (node_id != -1 && nodes[node_id].z < z) {
 
72
        last_id = node_id;
 
73
        node_id = nodes[node_id].next;
 
74
    }
 
75
 
 
76
    /* end of list, simply append */
 
77
    if (node_id == -1) {
 
78
        newnode_id = new_node();
 
79
        nodes[newnode_id].next = -1;
 
80
        nodes[newnode_id].z = z;
 
81
        nodes[last_id].next = newnode_id;
 
82
        return -1;
 
83
    }
 
84
    else if (node_id == head_id) {      /* pole position, insert as head */
 
85
        newnode_id = new_node();
 
86
        nodes[newnode_id].next = head_id;
 
87
        head_id = newnode_id;
 
88
        nodes[newnode_id].z = z;
 
89
        return (head_id);
 
90
    }
 
91
    else {                      /* somewhere in the middle, insert */
 
92
        newnode_id = new_node();
 
93
        nodes[newnode_id].z = z;
 
94
        nodes[newnode_id].next = node_id;
 
95
        nodes[last_id].next = newnode_id;
 
96
        return -1;
 
97
    }
 
98
}
 
99
 
 
100
int main(int argc, char *argv[])
 
101
{
 
102
    int out_fd;
 
103
    char *infile, *outmap;
 
104
    int percent;
 
105
    int method = -1;
 
106
    int bin_n, bin_min, bin_max, bin_sum, bin_sumsq, bin_index;
 
107
    double zrange_min, zrange_max, d_tmp;
 
108
    unsigned long estimated_lines;
 
109
 
 
110
    RASTER_MAP_TYPE rtype;
 
111
    struct History history;
 
112
    char title[64];
 
113
    void *n_array, *min_array, *max_array, *sum_array, *sumsq_array,
 
114
        *index_array;
 
115
    void *raster_row, *ptr;
 
116
    struct Cell_head region;
 
117
    int rows, cols;             /* scan box size */
 
118
    int row, col;               /* counters */
 
119
 
 
120
    int pass, npasses;
 
121
    unsigned long line, line_total;
 
122
    unsigned int counter;
 
123
    char buff[BUFFSIZE];
 
124
    double x, y, z;
 
125
    double pass_north, pass_south;
 
126
    int arr_row, arr_col;
 
127
    unsigned long count, count_total;
 
128
    int skipme, i;
 
129
    int point_class;
 
130
 
 
131
    double min = 0.0 / 0.0;     /* init as nan */
 
132
    double max = 0.0 / 0.0;     /* init as nan */
 
133
    double zscale = 1.0;
 
134
    size_t offset, n_offset;
 
135
    int n = 0;
 
136
    double sum = 0.;
 
137
    double sumsq = 0.;
 
138
    double variance, mean, skew, sumdev;
 
139
    int pth = 0;
 
140
    double trim = 0.0;
 
141
    double res = 0.0;
 
142
 
 
143
    int j, k;
 
144
    int head_id, node_id;
 
145
    int r_low, r_up;
 
146
 
 
147
    struct GModule *module;
 
148
    struct Option *input_opt, *output_opt, *percent_opt, *type_opt, *filter_opt, *class_opt;
 
149
    struct Option *method_opt, *zrange_opt, *zscale_opt;
 
150
    struct Option *trim_opt, *pth_opt, *res_opt;
 
151
    struct Flag *print_flag, *scan_flag, *shell_style, *over_flag, *extents_flag, *intens_flag;
 
152
 
 
153
    /* LAS */
 
154
    LASReaderH LAS_reader;
 
155
    LASHeaderH LAS_header;
 
156
    LASSRSH LAS_srs;
 
157
    LASPointH LAS_point;
 
158
    int return_filter;
 
159
 
 
160
    struct Key_Value *loc_proj_info = NULL, *loc_proj_units = NULL;
 
161
    struct Key_Value *proj_info, *proj_units;
 
162
    const char *projstr;
 
163
    struct Cell_head cellhd, loc_wind;
 
164
 
 
165
    unsigned int n_filtered;
 
166
 
 
167
    G_gisinit(argv[0]);
 
168
 
 
169
    module = G_define_module();
 
170
    G_add_keyword(_("raster"));
 
171
    G_add_keyword(_("import"));
 
172
    G_add_keyword(_("LIDAR"));
 
173
    module->description =
 
174
        _("Creates a raster map from LAS LiDAR points using univariate statistics.");
 
175
 
 
176
    input_opt = G_define_standard_option(G_OPT_F_INPUT);
 
177
    input_opt->label = _("LAS input file");
 
178
    input_opt->description = _("LiDAR input files in LAS format (*.las or *.laz)");
 
179
 
 
180
    output_opt = G_define_standard_option(G_OPT_R_OUTPUT);
 
181
 
 
182
    method_opt = G_define_option();
 
183
    method_opt->key = "method";
 
184
    method_opt->type = TYPE_STRING;
 
185
    method_opt->required = NO;
 
186
    method_opt->description = _("Statistic to use for raster values");
 
187
    method_opt->options =
 
188
        "n,min,max,range,sum,mean,stddev,variance,coeff_var,median,percentile,skewness,trimmean";
 
189
    method_opt->answer = "mean";
 
190
    method_opt->guisection = _("Statistic");
 
191
 
 
192
    type_opt = G_define_option();
 
193
    type_opt->key = "type";
 
194
    type_opt->type = TYPE_STRING;
 
195
    type_opt->required = NO;
 
196
    type_opt->options = "CELL,FCELL,DCELL";
 
197
    type_opt->answer = "FCELL";
 
198
    type_opt->description = _("Storage type for resultant raster map");
 
199
 
 
200
    zrange_opt = G_define_option();
 
201
    zrange_opt->key = "zrange";
 
202
    zrange_opt->type = TYPE_DOUBLE;
 
203
    zrange_opt->required = NO;
 
204
    zrange_opt->key_desc = "min,max";
 
205
    zrange_opt->description = _("Filter range for z data (min,max)");
 
206
 
 
207
    zscale_opt = G_define_option();
 
208
    zscale_opt->key = "zscale";
 
209
    zscale_opt->type = TYPE_DOUBLE;
 
210
    zscale_opt->required = NO;
 
211
    zscale_opt->answer = "1.0";
 
212
    zscale_opt->description = _("Scale to apply to z data");
 
213
 
 
214
    percent_opt = G_define_option();
 
215
    percent_opt->key = "percent";
 
216
    percent_opt->type = TYPE_INTEGER;
 
217
    percent_opt->required = NO;
 
218
    percent_opt->answer = "100";
 
219
    percent_opt->options = "1-100";
 
220
    percent_opt->description = _("Percent of map to keep in memory");
 
221
 
 
222
    /* I would prefer to call the following "percentile", but that has too
 
223
     * much namespace overlap with the "percent" option above */
 
224
    pth_opt = G_define_option();
 
225
    pth_opt->key = "pth";
 
226
    pth_opt->type = TYPE_INTEGER;
 
227
    pth_opt->required = NO;
 
228
    pth_opt->options = "1-100";
 
229
    pth_opt->description = _("pth percentile of the values");
 
230
    pth_opt->guisection = _("Statistic");
 
231
 
 
232
    trim_opt = G_define_option();
 
233
    trim_opt->key = "trim";
 
234
    trim_opt->type = TYPE_DOUBLE;
 
235
    trim_opt->required = NO;
 
236
    trim_opt->options = "0-50";
 
237
    trim_opt->description =
 
238
        _("Discard <trim> percent of the smallest and <trim> percent of the largest observations");
 
239
    trim_opt->guisection = _("Statistic");
 
240
 
 
241
    res_opt = G_define_option();
 
242
    res_opt->key = "resolution";
 
243
    res_opt->type = TYPE_DOUBLE;
 
244
    res_opt->required = NO;
 
245
    res_opt->description =
 
246
        _("Output raster resolution");
 
247
 
 
248
    filter_opt = G_define_option();
 
249
    filter_opt->key = "return_filter";
 
250
    filter_opt->type = TYPE_STRING;
 
251
    filter_opt->required = NO;
 
252
    filter_opt->label = _("Only import points of selected return type");
 
253
    filter_opt->description = _("If not specified, all points are imported");
 
254
    filter_opt->options = "first,last,mid";
 
255
 
 
256
    class_opt = G_define_option();
 
257
    class_opt->key = "class_filter";
 
258
    class_opt->type = TYPE_INTEGER;
 
259
    class_opt->multiple = YES;
 
260
    class_opt->required = NO;
 
261
    class_opt->label = _("Only import points of selected class(es)");
 
262
    class_opt->description = _("Input is comma separated integers. "
 
263
                               "If not specified, all points are imported.");
 
264
 
 
265
    print_flag = G_define_flag();
 
266
    print_flag->key = 'p';
 
267
    print_flag->description =
 
268
        _("Print LAS file info and exit");
 
269
    print_flag->suppress_required = YES;
 
270
 
 
271
    extents_flag = G_define_flag();
 
272
    extents_flag->key = 'e';
 
273
    extents_flag->description =
 
274
        _("Extend region extents based on new dataset");
 
275
 
 
276
    over_flag = G_define_flag();
 
277
    over_flag->key = 'o';
 
278
    over_flag->description =
 
279
        _("Override dataset projection (use location's projection)");
 
280
 
 
281
    scan_flag = G_define_flag();
 
282
    scan_flag->key = 's';
 
283
    scan_flag->description = _("Scan data file for extent then exit");
 
284
    scan_flag->suppress_required = YES;
 
285
 
 
286
    shell_style = G_define_flag();
 
287
    shell_style->key = 'g';
 
288
    shell_style->description =
 
289
        _("In scan mode, print using shell script style");
 
290
 
 
291
    intens_flag = G_define_flag();
 
292
    intens_flag->key = 'i';
 
293
    intens_flag->description =
 
294
        _("Import intensity values rather than z values");
 
295
 
 
296
    if (G_parser(argc, argv))
 
297
        exit(EXIT_FAILURE);
 
298
 
 
299
 
 
300
    /* parse input values */
 
301
    infile = input_opt->answer;
 
302
    outmap = output_opt->answer;
 
303
 
 
304
    if (shell_style->answer && !scan_flag->answer) {
 
305
        scan_flag->answer = 1; /* pointer not int, so set = shell_style->answer ? */
 
306
    }
 
307
 
 
308
    /* Don't crash on cmd line if file not found */
 
309
    if (access(infile, F_OK) != 0) {
 
310
        G_fatal_error(_("Input file <%s> does not exist"), infile);
 
311
    }
 
312
    /* Open LAS file*/
 
313
    LAS_reader = LASReader_Create(infile);
 
314
    if (LAS_reader == NULL) {
 
315
        G_fatal_error(_("Unable to open file <%s>"), infile);
 
316
    }
 
317
    
 
318
    LAS_header = LASReader_GetHeader(LAS_reader);
 
319
    if  (LAS_header == NULL) {
 
320
        G_fatal_error(_("Input file <%s> is not a LAS LiDAR point cloud"),
 
321
                        infile);
 
322
    }
 
323
 
 
324
    LAS_srs = LASHeader_GetSRS(LAS_header);
 
325
 
 
326
    /* Print LAS header */
 
327
    if (print_flag->answer) {
 
328
        /* print... */
 
329
        print_lasinfo(LAS_header, LAS_srs);
 
330
 
 
331
        LASSRS_Destroy(LAS_srs);
 
332
        LASHeader_Destroy(LAS_header);
 
333
        LASReader_Destroy(LAS_reader);
 
334
 
 
335
        exit(EXIT_SUCCESS);
 
336
    }
 
337
 
 
338
    return_filter = LAS_ALL;
 
339
    if (filter_opt->answer) {
 
340
        if (strcmp(filter_opt->answer, "first") == 0)
 
341
            return_filter = LAS_FIRST;
 
342
        else if (strcmp(filter_opt->answer, "last") == 0)
 
343
            return_filter = LAS_LAST;
 
344
        else if (strcmp(filter_opt->answer, "mid") == 0)
 
345
            return_filter = LAS_MID;
 
346
        else
 
347
            G_fatal_error(_("Unknown filter option <%s>"), filter_opt->answer);
 
348
    }
 
349
 
 
350
    /* Fetch input map projection in GRASS form. */
 
351
    proj_info = NULL;
 
352
    proj_units = NULL;
 
353
    projstr = LASSRS_GetWKT_CompoundOK(LAS_srs);
 
354
 
 
355
    if (TRUE) {
 
356
        int err = 0;
 
357
        char error_msg[8192];
 
358
 
 
359
        /* Projection only required for checking so convert non-interactively */
 
360
        if (GPJ_wkt_to_grass(&cellhd, &proj_info,
 
361
                             &proj_units, projstr, 0) < 0)
 
362
            G_warning(_("Unable to convert input map projection information to "
 
363
                       "GRASS format for checking"));
 
364
        
 
365
        /* Does the projection of the current location match the dataset? */
 
366
        /* G_get_window seems to be unreliable if the location has been changed */
 
367
        G_get_set_window(&loc_wind); /* TODO: v.in.lidar uses G_get_default_window() */
 
368
        /* fetch LOCATION PROJ info */
 
369
        if (loc_wind.proj != PROJECTION_XY) {
 
370
            loc_proj_info = G_get_projinfo();
 
371
            loc_proj_units = G_get_projunits();
 
372
        }
 
373
 
 
374
        if (over_flag->answer) {
 
375
            cellhd.proj = loc_wind.proj;
 
376
            cellhd.zone = loc_wind.zone;
 
377
            G_message(_("Over-riding projection check"));
 
378
        }
 
379
        else if (loc_wind.proj != cellhd.proj
 
380
                 || (err =
 
381
                     G_compare_projections(loc_proj_info, loc_proj_units,
 
382
                                           proj_info, proj_units)) != TRUE) {
 
383
            int i_value;
 
384
 
 
385
            strcpy(error_msg,
 
386
                   _("Projection of dataset does not"
 
387
                     " appear to match current location.\n\n"));
 
388
 
 
389
            /* TODO: output this info sorted by key: */
 
390
            if (loc_wind.proj != cellhd.proj || err != -2) {
 
391
                if (loc_proj_info != NULL) {
 
392
                    strcat(error_msg, _("GRASS LOCATION PROJ_INFO is:\n"));
 
393
                    for (i_value = 0; i_value < loc_proj_info->nitems;
 
394
                         i_value++)
 
395
                        sprintf(error_msg + strlen(error_msg), "%s: %s\n",
 
396
                                loc_proj_info->key[i_value],
 
397
                                loc_proj_info->value[i_value]);
 
398
                    strcat(error_msg, "\n");
 
399
                }
 
400
 
 
401
                if (proj_info != NULL) {
 
402
                    strcat(error_msg, _("Import dataset PROJ_INFO is:\n"));
 
403
                    for (i_value = 0; i_value < proj_info->nitems; i_value++)
 
404
                        sprintf(error_msg + strlen(error_msg), "%s: %s\n",
 
405
                                proj_info->key[i_value],
 
406
                                proj_info->value[i_value]);
 
407
                }
 
408
                else {
 
409
                    strcat(error_msg, _("Import dataset PROJ_INFO is:\n"));
 
410
                    if (cellhd.proj == PROJECTION_XY)
 
411
                        sprintf(error_msg + strlen(error_msg),
 
412
                                "Dataset proj = %d (unreferenced/unknown)\n",
 
413
                                cellhd.proj);
 
414
                    else if (cellhd.proj == PROJECTION_LL)
 
415
                        sprintf(error_msg + strlen(error_msg),
 
416
                                "Dataset proj = %d (lat/long)\n",
 
417
                                cellhd.proj);
 
418
                    else if (cellhd.proj == PROJECTION_UTM)
 
419
                        sprintf(error_msg + strlen(error_msg),
 
420
                                "Dataset proj = %d (UTM), zone = %d\n",
 
421
                                cellhd.proj, cellhd.zone);
 
422
                    else
 
423
                        sprintf(error_msg + strlen(error_msg),
 
424
                                "Dataset proj = %d (unknown), zone = %d\n",
 
425
                                cellhd.proj, cellhd.zone);
 
426
                }
 
427
            }
 
428
            else {
 
429
                if (loc_proj_units != NULL) {
 
430
                    strcat(error_msg, "GRASS LOCATION PROJ_UNITS is:\n");
 
431
                    for (i_value = 0; i_value < loc_proj_units->nitems;
 
432
                         i_value++)
 
433
                        sprintf(error_msg + strlen(error_msg), "%s: %s\n",
 
434
                                loc_proj_units->key[i_value],
 
435
                                loc_proj_units->value[i_value]);
 
436
                    strcat(error_msg, "\n");
 
437
                }
 
438
 
 
439
                if (proj_units != NULL) {
 
440
                    strcat(error_msg, "Import dataset PROJ_UNITS is:\n");
 
441
                    for (i_value = 0; i_value < proj_units->nitems; i_value++)
 
442
                        sprintf(error_msg + strlen(error_msg), "%s: %s\n",
 
443
                                proj_units->key[i_value],
 
444
                                proj_units->value[i_value]);
 
445
                }
 
446
            }
 
447
            sprintf(error_msg + strlen(error_msg),
 
448
                    _("\nYou can use the -o flag to %s to override this projection check.\n"),
 
449
                    G_program_name());
 
450
            strcat(error_msg,
 
451
                   _("Consider generating a new location with 'location' parameter"
 
452
                    " from input data set.\n"));
 
453
            G_fatal_error(error_msg);
 
454
        }
 
455
        else if (!shell_style->answer) {
 
456
            G_message(_("Projection of input dataset and current location "
 
457
                        "appear to match"));
 
458
        }
 
459
    }
 
460
 
 
461
    percent = atoi(percent_opt->answer);
 
462
    zscale = atof(zscale_opt->answer);
 
463
 
 
464
    /* parse zrange */
 
465
    if (zrange_opt->answer != NULL) {
 
466
        if (zrange_opt->answers[0] == NULL)
 
467
            G_fatal_error(_("Invalid zrange"));
 
468
 
 
469
        sscanf(zrange_opt->answers[0], "%lf", &zrange_min);
 
470
        sscanf(zrange_opt->answers[1], "%lf", &zrange_max);
 
471
 
 
472
        if (zrange_min > zrange_max) {
 
473
            d_tmp = zrange_max;
 
474
            zrange_max = zrange_min;
 
475
            zrange_min = d_tmp;
 
476
        }
 
477
    }
 
478
 
 
479
    /* figure out what maps we need in memory */
 
480
    /*  n               n
 
481
       min              min
 
482
       max              max
 
483
       range            min max         max - min
 
484
       sum              sum
 
485
       mean             sum n           sum/n
 
486
       stddev           sum sumsq n     sqrt((sumsq - sum*sum/n)/n)
 
487
       variance         sum sumsq n     (sumsq - sum*sum/n)/n
 
488
       coeff_var        sum sumsq n     sqrt((sumsq - sum*sum/n)/n) / (sum/n)
 
489
       median           n               array index to linked list
 
490
       percentile       n               array index to linked list
 
491
       skewness         n               array index to linked list
 
492
       trimmean         n               array index to linked list
 
493
     */
 
494
    bin_n = FALSE;
 
495
    bin_min = FALSE;
 
496
    bin_max = FALSE;
 
497
    bin_sum = FALSE;
 
498
    bin_sumsq = FALSE;
 
499
    bin_index = FALSE;
 
500
 
 
501
    n_array = NULL;
 
502
    min_array = NULL;
 
503
    max_array = NULL;
 
504
    sum_array = NULL;
 
505
    sumsq_array = NULL;
 
506
    index_array = NULL;
 
507
    
 
508
    if (strcmp(method_opt->answer, "n") == 0) {
 
509
        method = METHOD_N;
 
510
        bin_n = TRUE;
 
511
    }
 
512
    if (strcmp(method_opt->answer, "min") == 0) {
 
513
        method = METHOD_MIN;
 
514
        bin_min = TRUE;
 
515
    }
 
516
    if (strcmp(method_opt->answer, "max") == 0) {
 
517
        method = METHOD_MAX;
 
518
        bin_max = TRUE;
 
519
    }
 
520
    if (strcmp(method_opt->answer, "range") == 0) {
 
521
        method = METHOD_RANGE;
 
522
        bin_min = TRUE;
 
523
        bin_max = TRUE;
 
524
    }
 
525
    if (strcmp(method_opt->answer, "sum") == 0) {
 
526
        method = METHOD_SUM;
 
527
        bin_sum = TRUE;
 
528
    }
 
529
    if (strcmp(method_opt->answer, "mean") == 0) {
 
530
        method = METHOD_MEAN;
 
531
        bin_sum = TRUE;
 
532
        bin_n = TRUE;
 
533
    }
 
534
    if (strcmp(method_opt->answer, "stddev") == 0) {
 
535
        method = METHOD_STDDEV;
 
536
        bin_sum = TRUE;
 
537
        bin_sumsq = TRUE;
 
538
        bin_n = TRUE;
 
539
    }
 
540
    if (strcmp(method_opt->answer, "variance") == 0) {
 
541
        method = METHOD_VARIANCE;
 
542
        bin_sum = TRUE;
 
543
        bin_sumsq = TRUE;
 
544
        bin_n = TRUE;
 
545
    }
 
546
    if (strcmp(method_opt->answer, "coeff_var") == 0) {
 
547
        method = METHOD_COEFF_VAR;
 
548
        bin_sum = TRUE;
 
549
        bin_sumsq = TRUE;
 
550
        bin_n = TRUE;
 
551
    }
 
552
    if (strcmp(method_opt->answer, "median") == 0) {
 
553
        method = METHOD_MEDIAN;
 
554
        bin_index = TRUE;
 
555
    }
 
556
    if (strcmp(method_opt->answer, "percentile") == 0) {
 
557
        if (pth_opt->answer != NULL)
 
558
            pth = atoi(pth_opt->answer);
 
559
        else
 
560
            G_fatal_error(_("Unable to calculate percentile without the pth option specified!"));
 
561
        method = METHOD_PERCENTILE;
 
562
        bin_index = TRUE;
 
563
    }
 
564
    if (strcmp(method_opt->answer, "skewness") == 0) {
 
565
        method = METHOD_SKEWNESS;
 
566
        bin_index = TRUE;
 
567
    }
 
568
    if (strcmp(method_opt->answer, "trimmean") == 0) {
 
569
        if (trim_opt->answer != NULL)
 
570
            trim = atof(trim_opt->answer) / 100.0;
 
571
        else
 
572
            G_fatal_error(_("Unable to calculate trimmed mean without the trim option specified!"));
 
573
        method = METHOD_TRIMMEAN;
 
574
        bin_index = TRUE;
 
575
    }
 
576
 
 
577
    if (strcmp("CELL", type_opt->answer) == 0)
 
578
        rtype = CELL_TYPE;
 
579
    else if (strcmp("DCELL", type_opt->answer) == 0)
 
580
        rtype = DCELL_TYPE;
 
581
    else
 
582
        rtype = FCELL_TYPE;
 
583
 
 
584
    if (method == METHOD_N)
 
585
        rtype = CELL_TYPE;
 
586
 
 
587
 
 
588
    Rast_get_window(&region);
 
589
 
 
590
 
 
591
    if (scan_flag->answer || extents_flag->answer) {
 
592
        if (zrange_opt->answer)
 
593
            G_warning(_("zrange will not be taken into account during scan"));
 
594
 
 
595
        scan_bounds(LAS_reader, shell_style->answer, extents_flag->answer,
 
596
                    zscale, &region);
 
597
 
 
598
        if (!extents_flag->answer) {
 
599
            LASHeader_Destroy(LAS_header);
 
600
            LASReader_Destroy(LAS_reader);
 
601
 
 
602
            exit(EXIT_SUCCESS);
 
603
        }
 
604
    }
 
605
    
 
606
    if (res_opt->answer) {
 
607
        /* align to resolution */
 
608
        res = atof(res_opt->answer);
 
609
 
 
610
        if (!G_scan_resolution(res_opt->answer, &res, region.proj))
 
611
            G_fatal_error(_("Invalid input <%s=%s>"), res_opt->key, res_opt->answer);
 
612
 
 
613
        if (res <= 0)
 
614
            G_fatal_error(_("Option '%s' must be > 0.0"), res_opt->key);
 
615
        
 
616
        region.ns_res = region.ew_res = res;
 
617
 
 
618
        region.north = ceil(region.north / res) * res;
 
619
        region.south = floor(region.south / res) * res;
 
620
        region.east = ceil(region.east / res) * res;
 
621
        region.west = floor(region.west / res) * res;
 
622
 
 
623
        G_adjust_Cell_head(&region, 0, 0);
 
624
    }
 
625
    else if (extents_flag->answer) {
 
626
        /* align to current region */
 
627
        Rast_align_window(&region, &loc_wind);
 
628
    }
 
629
    Rast_set_output_window(&region);
 
630
 
 
631
    rows = (int)(region.rows * (percent / 100.0));
 
632
    cols = region.cols;
 
633
 
 
634
    G_debug(2, "region.n=%f  region.s=%f  region.ns_res=%f", region.north,
 
635
            region.south, region.ns_res);
 
636
    G_debug(2, "region.rows=%d  [box_rows=%d]  region.cols=%d", region.rows,
 
637
            rows, region.cols);
 
638
 
 
639
    npasses = (int)ceil(1.0 * region.rows / rows);
 
640
 
 
641
    if (!scan_flag->answer) {
 
642
        /* check if rows * (cols + 1) go into a size_t */
 
643
        if (sizeof(size_t) < 8) {
 
644
            double dsize = rows * (cols + 1);
 
645
            
 
646
            if (dsize != (size_t)rows * (cols + 1))
 
647
                G_fatal_error(_("Unable to process the hole map at once. "
 
648
                                "Please set the %s option to some value lower than 100."),
 
649
                                percent_opt->key);
 
650
        }
 
651
        /* allocate memory (test for enough before we start) */
 
652
        if (bin_n)
 
653
            n_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(CELL_TYPE));
 
654
        if (bin_min)
 
655
            min_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
 
656
        if (bin_max)
 
657
            max_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
 
658
        if (bin_sum)
 
659
            sum_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
 
660
        if (bin_sumsq)
 
661
            sumsq_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
 
662
        if (bin_index)
 
663
            index_array =
 
664
                G_calloc((size_t)rows * (cols + 1), Rast_cell_size(CELL_TYPE));
 
665
 
 
666
        /* and then free it again */
 
667
        if (bin_n)
 
668
            G_free(n_array);
 
669
        if (bin_min)
 
670
            G_free(min_array);
 
671
        if (bin_max)
 
672
            G_free(max_array);
 
673
        if (bin_sum)
 
674
            G_free(sum_array);
 
675
        if (bin_sumsq)
 
676
            G_free(sumsq_array);
 
677
        if (bin_index)
 
678
            G_free(index_array);
 
679
 
 
680
        /** end memory test **/
 
681
    }
 
682
 
 
683
 
 
684
    /* open output map */
 
685
    out_fd = Rast_open_new(outmap, rtype);
 
686
 
 
687
    estimated_lines = LASHeader_GetPointRecordsCount(LAS_header);
 
688
 
 
689
    /* allocate memory for a single row of output data */
 
690
    raster_row = Rast_allocate_output_buf(rtype);
 
691
 
 
692
    G_message(_("Reading data ..."));
 
693
 
 
694
    count_total = line_total = 0;
 
695
 
 
696
    /* init northern border */
 
697
    pass_south = region.north;
 
698
 
 
699
    /* main binning loop(s) */
 
700
    for (pass = 1; pass <= npasses; pass++) {
 
701
        LASError LAS_error;
 
702
 
 
703
        if (npasses > 1)
 
704
            G_message(_("Pass #%d (of %d) ..."), pass, npasses);
 
705
 
 
706
        LAS_error = LASReader_Seek(LAS_reader, 0);
 
707
        
 
708
        if (LAS_error != LE_None)
 
709
            G_fatal_error(_("Could not rewind input file"));
 
710
 
 
711
        /* figure out segmentation */
 
712
        pass_north = pass_south;  /* exact copy to avoid fp errors */
 
713
        pass_south = region.north - pass * rows * region.ns_res;
 
714
        if (pass == npasses) {
 
715
            rows = region.rows - (pass - 1) * rows;
 
716
            pass_south = region.south; /* exact copy to avoid fp errors */
 
717
        }
 
718
 
 
719
        G_debug(2, "pass=%d/%d  pass_n=%f  pass_s=%f  rows=%d",
 
720
                pass, npasses, pass_north, pass_south, rows);
 
721
 
 
722
 
 
723
        if (bin_n) {
 
724
            G_debug(2, "allocating n_array");
 
725
            n_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(CELL_TYPE));
 
726
            blank_array(n_array, rows, cols, CELL_TYPE, 0);
 
727
        }
 
728
        if (bin_min) {
 
729
            G_debug(2, "allocating min_array");
 
730
            min_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
 
731
            blank_array(min_array, rows, cols, rtype, -1);      /* fill with NULLs */
 
732
        }
 
733
        if (bin_max) {
 
734
            G_debug(2, "allocating max_array");
 
735
            max_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
 
736
            blank_array(max_array, rows, cols, rtype, -1);      /* fill with NULLs */
 
737
        }
 
738
        if (bin_sum) {
 
739
            G_debug(2, "allocating sum_array");
 
740
            sum_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
 
741
            blank_array(sum_array, rows, cols, rtype, 0);
 
742
        }
 
743
        if (bin_sumsq) {
 
744
            G_debug(2, "allocating sumsq_array");
 
745
            sumsq_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
 
746
            blank_array(sumsq_array, rows, cols, rtype, 0);
 
747
        }
 
748
        if (bin_index) {
 
749
            G_debug(2, "allocating index_array");
 
750
            index_array =
 
751
                G_calloc((size_t)rows * (cols + 1), Rast_cell_size(CELL_TYPE));
 
752
            blank_array(index_array, rows, cols, CELL_TYPE, -1);        /* fill with NULLs */
 
753
        }
 
754
 
 
755
        line = 0;
 
756
        count = 0;
 
757
        counter = 0;
 
758
        G_percent_reset();
 
759
 
 
760
        while ((LAS_point = LASReader_GetNextPoint(LAS_reader)) != NULL) {
 
761
            line++;
 
762
            counter++;
 
763
 
 
764
            if (counter == 100000) {    /* speed */
 
765
                if (line < estimated_lines)
 
766
                    G_percent(line, estimated_lines, 3);
 
767
                counter = 0;
 
768
            }
 
769
 
 
770
            if (!LASPoint_IsValid(LAS_point)) {
 
771
                continue;
 
772
            }
 
773
 
 
774
            x = LASPoint_GetX(LAS_point);
 
775
            y = LASPoint_GetY(LAS_point);
 
776
            if (intens_flag->answer)
 
777
                /* use z variable here to allow for scaling of intensity below */
 
778
                z = LASPoint_GetIntensity(LAS_point);
 
779
            else
 
780
                z = LASPoint_GetZ(LAS_point);
 
781
 
 
782
        if (return_filter != LAS_ALL) {
 
783
            int return_no = LASPoint_GetReturnNumber(LAS_point);
 
784
            int n_returns = LASPoint_GetNumberOfReturns(LAS_point);
 
785
            skipme = 1;
 
786
 
 
787
            switch (return_filter) {
 
788
            case LAS_FIRST:
 
789
                if (return_no == 1)
 
790
                    skipme = 0;
 
791
                break;
 
792
            case LAS_MID:
 
793
                if (return_no > 1 && return_no < n_returns)
 
794
                    skipme = 0;
 
795
                break;
 
796
            case LAS_LAST:
 
797
                if (n_returns > 1 && return_no == n_returns)
 
798
                    skipme = 0;
 
799
                break;
 
800
            }
 
801
 
 
802
            if (skipme) {
 
803
                n_filtered++;
 
804
                continue;
 
805
            }
 
806
        }
 
807
        if (class_opt->answer) {
 
808
            point_class = (int) LASPoint_GetClassification(LAS_point);
 
809
            i = 0;
 
810
            skipme = TRUE;
 
811
            while (class_opt->answers[i]) {
 
812
                if (point_class == atoi(class_opt->answers[i])) {
 
813
                    skipme = FALSE;
 
814
                    break;
 
815
                }
 
816
                i++;
 
817
            }
 
818
            if (skipme) {
 
819
                continue;
 
820
            }
 
821
        }
 
822
 
 
823
            if (y <= pass_south || y > pass_north) {
 
824
                continue;
 
825
            }
 
826
            if (x < region.west || x >= region.east) {
 
827
                continue;
 
828
            }
 
829
 
 
830
            z = z * zscale;
 
831
 
 
832
            if (zrange_opt->answer) {
 
833
                if (z < zrange_min || z > zrange_max) {
 
834
                    continue;
 
835
                }
 
836
            }
 
837
 
 
838
            count++;
 
839
            /*          G_debug(5, "x: %f, y: %f, z: %f", x, y, z); */
 
840
 
 
841
            /* find the bin in the current array box */
 
842
            arr_row = (int)((pass_north - y) / region.ns_res);
 
843
            arr_col = (int)((x - region.west) / region.ew_res);
 
844
 
 
845
            if (bin_n)
 
846
                update_n(n_array, cols, arr_row, arr_col);
 
847
            if (bin_min)
 
848
                update_min(min_array, cols, arr_row, arr_col, rtype, z);
 
849
            if (bin_max)
 
850
                update_max(max_array, cols, arr_row, arr_col, rtype, z);
 
851
            if (bin_sum)
 
852
                update_sum(sum_array, cols, arr_row, arr_col, rtype, z);
 
853
            if (bin_sumsq)
 
854
                update_sumsq(sumsq_array, cols, arr_row, arr_col, rtype, z);
 
855
            if (bin_index) {
 
856
                ptr = index_array;
 
857
                ptr =
 
858
                    G_incr_void_ptr(ptr,
 
859
                                    ((arr_row * cols) +
 
860
                                     arr_col) * Rast_cell_size(CELL_TYPE));
 
861
 
 
862
                if (Rast_is_null_value(ptr, CELL_TYPE)) {       /* first node */
 
863
                    head_id = new_node();
 
864
                    nodes[head_id].next = -1;
 
865
                    nodes[head_id].z = z;
 
866
                    Rast_set_c_value(ptr, head_id, CELL_TYPE);  /* store index to head */
 
867
                }
 
868
                else {          /* head is already there */
 
869
 
 
870
                    head_id = Rast_get_c_value(ptr, CELL_TYPE); /* get index to head */
 
871
                    head_id = add_node(head_id, z);
 
872
                    if (head_id != -1)
 
873
                        Rast_set_c_value(ptr, head_id, CELL_TYPE);      /* store index to head */
 
874
                }
 
875
            }
 
876
        }                       /* while !EOF */
 
877
 
 
878
        G_percent(1, 1, 1);     /* flush */
 
879
        G_debug(2, "pass %d finished, %lu coordinates in box", pass, count);
 
880
        count_total += count;
 
881
        line_total += line;
 
882
 
 
883
        /* calc stats and output */
 
884
        G_message(_("Writing to map ..."));
 
885
        for (row = 0; row < rows; row++) {
 
886
 
 
887
            switch (method) {
 
888
            case METHOD_N:      /* n is a straight copy */
 
889
                Rast_raster_cpy(raster_row,
 
890
                             n_array +
 
891
                             (row * cols * Rast_cell_size(CELL_TYPE)), cols,
 
892
                             CELL_TYPE);
 
893
                break;
 
894
 
 
895
            case METHOD_MIN:
 
896
                Rast_raster_cpy(raster_row,
 
897
                             min_array + (row * cols * Rast_cell_size(rtype)),
 
898
                             cols, rtype);
 
899
                break;
 
900
 
 
901
            case METHOD_MAX:
 
902
                Rast_raster_cpy(raster_row,
 
903
                             max_array + (row * cols * Rast_cell_size(rtype)),
 
904
                             cols, rtype);
 
905
                break;
 
906
 
 
907
            case METHOD_SUM:
 
908
                Rast_raster_cpy(raster_row,
 
909
                             sum_array + (row * cols * Rast_cell_size(rtype)),
 
910
                             cols, rtype);
 
911
                break;
 
912
 
 
913
            case METHOD_RANGE:  /* (max-min) */
 
914
                ptr = raster_row;
 
915
                for (col = 0; col < cols; col++) {
 
916
                    offset = (row * cols + col) * Rast_cell_size(rtype);
 
917
                    min = Rast_get_d_value(min_array + offset, rtype);
 
918
                    max = Rast_get_d_value(max_array + offset, rtype);
 
919
                    Rast_set_d_value(ptr, max - min, rtype);
 
920
                    ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
 
921
                }
 
922
                break;
 
923
 
 
924
            case METHOD_MEAN:   /* (sum / n) */
 
925
                ptr = raster_row;
 
926
                for (col = 0; col < cols; col++) {
 
927
                    offset = (row * cols + col) * Rast_cell_size(rtype);
 
928
                    n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
 
929
                    n = Rast_get_c_value(n_array + n_offset, CELL_TYPE);
 
930
                    sum = Rast_get_d_value(sum_array + offset, rtype);
 
931
 
 
932
                    if (n == 0)
 
933
                        Rast_set_null_value(ptr, 1, rtype);
 
934
                    else
 
935
                        Rast_set_d_value(ptr, (sum / n), rtype);
 
936
 
 
937
                    ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
 
938
                }
 
939
                break;
 
940
 
 
941
            case METHOD_STDDEV: /*  sqrt(variance)        */
 
942
            case METHOD_VARIANCE:       /*  (sumsq - sum*sum/n)/n */
 
943
            case METHOD_COEFF_VAR:      /*  100 * stdev / mean    */
 
944
                ptr = raster_row;
 
945
                for (col = 0; col < cols; col++) {
 
946
                    offset = (row * cols + col) * Rast_cell_size(rtype);
 
947
                    n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
 
948
                    n = Rast_get_c_value(n_array + n_offset, CELL_TYPE);
 
949
                    sum = Rast_get_d_value(sum_array + offset, rtype);
 
950
                    sumsq = Rast_get_d_value(sumsq_array + offset, rtype);
 
951
 
 
952
                    if (n == 0)
 
953
                        Rast_set_null_value(ptr, 1, rtype);
 
954
                    else if (n == 1)
 
955
                        Rast_set_d_value(ptr, 0.0, rtype);
 
956
                    else {
 
957
                        variance = (sumsq - sum * sum / n) / n;
 
958
                        if (variance < GRASS_EPSILON)
 
959
                            variance = 0.0;
 
960
 
 
961
                        /* nan test */
 
962
                        if (variance != variance)
 
963
                            Rast_set_null_value(ptr, 1, rtype);
 
964
                        else {
 
965
 
 
966
                            if (method == METHOD_STDDEV)
 
967
                                variance = sqrt(variance);
 
968
 
 
969
                            else if (method == METHOD_COEFF_VAR)
 
970
                                variance = 100 * sqrt(variance) / (sum / n);
 
971
 
 
972
                            /* nan test */
 
973
                            if (variance != variance)
 
974
                                variance = 0.0; /* OK for n > 0 ?*/
 
975
 
 
976
                            Rast_set_d_value(ptr, variance, rtype);
 
977
                        }
 
978
 
 
979
                    }
 
980
                    ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
 
981
                }
 
982
                break;
 
983
            case METHOD_MEDIAN: /* median, if only one point in cell we will use that */
 
984
                ptr = raster_row;
 
985
                for (col = 0; col < cols; col++) {
 
986
                    n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
 
987
                    if (Rast_is_null_value(index_array + n_offset, CELL_TYPE))  /* no points in cell */
 
988
                        Rast_set_null_value(ptr, 1, rtype);
 
989
                    else {      /* one or more points in cell */
 
990
 
 
991
                        head_id =
 
992
                            Rast_get_c_value(index_array + n_offset,
 
993
                                                 CELL_TYPE);
 
994
                        node_id = head_id;
 
995
 
 
996
                        n = 0;
 
997
 
 
998
                        while (node_id != -1) { /* count number of points in cell */
 
999
                            n++;
 
1000
                            node_id = nodes[node_id].next;
 
1001
                        }
 
1002
 
 
1003
                        if (n == 1)     /* only one point, use that */
 
1004
                            Rast_set_d_value(ptr, nodes[head_id].z,
 
1005
                                                 rtype);
 
1006
                        else if (n % 2 != 0) {  /* odd number of points: median_i = (n + 1) / 2 */
 
1007
                            n = (n + 1) / 2;
 
1008
                            node_id = head_id;
 
1009
                            for (j = 1; j < n; j++)     /* get "median element" */
 
1010
                                node_id = nodes[node_id].next;
 
1011
 
 
1012
                            Rast_set_d_value(ptr, nodes[node_id].z,
 
1013
                                                 rtype);
 
1014
                        }
 
1015
                        else {  /* even number of points: median = (val_below + val_above) / 2 */
 
1016
 
 
1017
                            z = (n + 1) / 2.0;
 
1018
                            n = floor(z);
 
1019
                            node_id = head_id;
 
1020
                            for (j = 1; j < n; j++)     /* get element "below" */
 
1021
                                node_id = nodes[node_id].next;
 
1022
 
 
1023
                            z = (nodes[node_id].z +
 
1024
                                 nodes[nodes[node_id].next].z) / 2;
 
1025
                            Rast_set_d_value(ptr, z, rtype);
 
1026
                        }
 
1027
                    }
 
1028
                    ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
 
1029
                }
 
1030
                break;
 
1031
            case METHOD_PERCENTILE:     /* rank = (pth*(n+1))/100; interpolate linearly */
 
1032
                ptr = raster_row;
 
1033
                for (col = 0; col < cols; col++) {
 
1034
                    n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
 
1035
                    if (Rast_is_null_value(index_array + n_offset, CELL_TYPE))  /* no points in cell */
 
1036
                        Rast_set_null_value(ptr, 1, rtype);
 
1037
                    else {
 
1038
                        head_id =
 
1039
                            Rast_get_c_value(index_array + n_offset,
 
1040
                                                 CELL_TYPE);
 
1041
                        node_id = head_id;
 
1042
                        n = 0;
 
1043
 
 
1044
                        while (node_id != -1) { /* count number of points in cell */
 
1045
                            n++;
 
1046
                            node_id = nodes[node_id].next;
 
1047
                        }
 
1048
 
 
1049
                        z = (pth * (n + 1)) / 100.0;
 
1050
                        r_low = floor(z);       /* lower rank */
 
1051
                        if (r_low < 1)
 
1052
                            r_low = 1;
 
1053
                        else if (r_low > n)
 
1054
                            r_low = n;
 
1055
 
 
1056
                        r_up = ceil(z); /* upper rank */
 
1057
                        if (r_up > n)
 
1058
                            r_up = n;
 
1059
 
 
1060
                        node_id = head_id;
 
1061
                        for (j = 1; j < r_low; j++)     /* search lower value */
 
1062
                            node_id = nodes[node_id].next;
 
1063
 
 
1064
                        z = nodes[node_id].z;   /* save lower value */
 
1065
                        node_id = head_id;
 
1066
                        for (j = 1; j < r_up; j++)      /* search upper value */
 
1067
                            node_id = nodes[node_id].next;
 
1068
 
 
1069
                        z = (z + nodes[node_id].z) / 2;
 
1070
                        Rast_set_d_value(ptr, z, rtype);
 
1071
                    }
 
1072
                    ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
 
1073
                }
 
1074
                break;
 
1075
            case METHOD_SKEWNESS:       /* skewness = sum(xi-mean)^3/(N-1)*s^3 */
 
1076
                ptr = raster_row;
 
1077
                for (col = 0; col < cols; col++) {
 
1078
                    n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
 
1079
                    if (Rast_is_null_value(index_array + n_offset, CELL_TYPE))  /* no points in cell */
 
1080
                        Rast_set_null_value(ptr, 1, rtype);
 
1081
                    else {
 
1082
                        head_id =
 
1083
                            Rast_get_c_value(index_array + n_offset,
 
1084
                                                 CELL_TYPE);
 
1085
                        node_id = head_id;
 
1086
 
 
1087
                        n = 0;  /* count */
 
1088
                        sum = 0.0;      /* sum */
 
1089
                        sumsq = 0.0;    /* sum of squares */
 
1090
                        sumdev = 0.0;   /* sum of (xi - mean)^3 */
 
1091
                        skew = 0.0;     /* skewness */
 
1092
 
 
1093
                        while (node_id != -1) {
 
1094
                            z = nodes[node_id].z;
 
1095
                            n++;
 
1096
                            sum += z;
 
1097
                            sumsq += (z * z);
 
1098
                            node_id = nodes[node_id].next;
 
1099
                        }
 
1100
 
 
1101
                        if (n > 1) {    /* if n == 1, skew is "0.0" */
 
1102
                            mean = sum / n;
 
1103
                            node_id = head_id;
 
1104
                            while (node_id != -1) {
 
1105
                                z = nodes[node_id].z;
 
1106
                                sumdev += pow((z - mean), 3);
 
1107
                                node_id = nodes[node_id].next;
 
1108
                            }
 
1109
 
 
1110
                            variance = (sumsq - sum * sum / n) / n;
 
1111
                            if (variance < GRASS_EPSILON)
 
1112
                                skew = 0.0;
 
1113
                            else
 
1114
                                skew =
 
1115
                                    sumdev / ((n - 1) *
 
1116
                                              pow(sqrt(variance), 3));
 
1117
                        }
 
1118
                        Rast_set_d_value(ptr, skew, rtype);
 
1119
                    }
 
1120
                    ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
 
1121
                }
 
1122
                break;
 
1123
            case METHOD_TRIMMEAN:
 
1124
                ptr = raster_row;
 
1125
                for (col = 0; col < cols; col++) {
 
1126
                    n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
 
1127
                    if (Rast_is_null_value(index_array + n_offset, CELL_TYPE))  /* no points in cell */
 
1128
                        Rast_set_null_value(ptr, 1, rtype);
 
1129
                    else {
 
1130
                        head_id =
 
1131
                            Rast_get_c_value(index_array + n_offset,
 
1132
                                                 CELL_TYPE);
 
1133
 
 
1134
                        node_id = head_id;
 
1135
                        n = 0;
 
1136
                        while (node_id != -1) { /* count number of points in cell */
 
1137
                            n++;
 
1138
                            node_id = nodes[node_id].next;
 
1139
                        }
 
1140
 
 
1141
                        if (1 == n)
 
1142
                            mean = nodes[head_id].z;
 
1143
                        else {
 
1144
                            k = floor(trim * n + 0.5);  /* number of ranks to discard on each tail */
 
1145
 
 
1146
                            if (k > 0 && (n - 2 * k) > 0) {     /* enough elements to discard */
 
1147
                                node_id = head_id;
 
1148
                                for (j = 0; j < k; j++) /* move to first rank to consider */
 
1149
                                    node_id = nodes[node_id].next;
 
1150
 
 
1151
                                j = k + 1;
 
1152
                                k = n - k;
 
1153
                                n = 0;
 
1154
                                sum = 0.0;
 
1155
 
 
1156
                                while (j <= k) {        /* get values in interval */
 
1157
                                    n++;
 
1158
                                    sum += nodes[node_id].z;
 
1159
                                    node_id = nodes[node_id].next;
 
1160
                                    j++;
 
1161
                                }
 
1162
                            }
 
1163
                            else {
 
1164
                                node_id = head_id;
 
1165
                                n = 0;
 
1166
                                sum = 0.0;
 
1167
                                while (node_id != -1) {
 
1168
                                    n++;
 
1169
                                    sum += nodes[node_id].z;
 
1170
                                    node_id = nodes[node_id].next;
 
1171
                                }
 
1172
                            }
 
1173
                            mean = sum / n;
 
1174
                        }
 
1175
                        Rast_set_d_value(ptr, mean, rtype);
 
1176
                    }
 
1177
                    ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
 
1178
                }
 
1179
                break;
 
1180
 
 
1181
            default:
 
1182
                G_fatal_error("?");
 
1183
            }
 
1184
 
 
1185
            /* write out line of raster data */
 
1186
            Rast_put_row(out_fd, raster_row, rtype);
 
1187
        }
 
1188
 
 
1189
        /* free memory */
 
1190
        if (bin_n)
 
1191
            G_free(n_array);
 
1192
        if (bin_min)
 
1193
            G_free(min_array);
 
1194
        if (bin_max)
 
1195
            G_free(max_array);
 
1196
        if (bin_sum)
 
1197
            G_free(sum_array);
 
1198
        if (bin_sumsq)
 
1199
            G_free(sumsq_array);
 
1200
        if (bin_index) {
 
1201
            G_free(index_array);
 
1202
            G_free(nodes);
 
1203
            num_nodes = 0;
 
1204
            max_nodes = 0;
 
1205
            nodes = NULL;
 
1206
        }
 
1207
 
 
1208
    }                           /* passes loop */
 
1209
 
 
1210
    G_percent(1, 1, 1);         /* flush */
 
1211
    G_free(raster_row);
 
1212
 
 
1213
    /* close input LAS file */
 
1214
    LASHeader_Destroy(LAS_header);
 
1215
    LASReader_Destroy(LAS_reader);
 
1216
 
 
1217
    /* close raster file & write history */
 
1218
    Rast_close(out_fd);
 
1219
 
 
1220
    sprintf(title, "Raw x,y,z data binned into a raster grid by cell %s",
 
1221
            method_opt->answer);
 
1222
    Rast_put_cell_title(outmap, title);
 
1223
 
 
1224
    Rast_short_history(outmap, "raster", &history);
 
1225
    Rast_command_history(&history);
 
1226
    Rast_set_history(&history, HIST_DATSRC_1, infile);
 
1227
    Rast_write_history(outmap, &history);
 
1228
 
 
1229
 
 
1230
    sprintf(buff, _("%lu points found in region."), count_total);
 
1231
    G_done_msg(buff);
 
1232
    G_debug(1, "Processed %lu points.", line_total);
 
1233
 
 
1234
    exit(EXIT_SUCCESS);
 
1235
 
 
1236
}
 
1237
 
 
1238
void print_lasinfo(LASHeaderH LAS_header, LASSRSH LAS_srs)
 
1239
{
 
1240
    char *las_srs_proj4 = LASSRS_GetProj4(LAS_srs);
 
1241
    int las_point_format = LASHeader_GetDataFormatId(LAS_header);
 
1242
 
 
1243
    fprintf(stdout, "\nUsing LAS Library Version '%s'\n\n",
 
1244
                    LAS_GetFullVersion());
 
1245
    fprintf(stdout, "LAS File Version:                  %d.%d\n",
 
1246
                    LASHeader_GetVersionMajor(LAS_header),
 
1247
                    LASHeader_GetVersionMinor(LAS_header));
 
1248
    fprintf(stdout, "System ID:                         '%s'\n",
 
1249
                    LASHeader_GetSystemId(LAS_header));
 
1250
    fprintf(stdout, "Generating Software:               '%s'\n",
 
1251
                    LASHeader_GetSoftwareId(LAS_header));
 
1252
    fprintf(stdout, "File Creation Day/Year:            %d/%d\n",
 
1253
                    LASHeader_GetCreationDOY(LAS_header),
 
1254
                    LASHeader_GetCreationYear(LAS_header));
 
1255
    fprintf(stdout, "Point Data Format:                 %d\n",
 
1256
                    las_point_format);
 
1257
    fprintf(stdout, "Number of Point Records:           %d\n",
 
1258
                    LASHeader_GetPointRecordsCount(LAS_header));
 
1259
    fprintf(stdout, "Number of Points by Return:        %d %d %d %d %d\n",
 
1260
                    LASHeader_GetPointRecordsByReturnCount(LAS_header, 0),
 
1261
                    LASHeader_GetPointRecordsByReturnCount(LAS_header, 1),
 
1262
                    LASHeader_GetPointRecordsByReturnCount(LAS_header, 2),
 
1263
                    LASHeader_GetPointRecordsByReturnCount(LAS_header, 3),
 
1264
                    LASHeader_GetPointRecordsByReturnCount(LAS_header, 4));
 
1265
    fprintf(stdout, "Scale Factor X Y Z:                %g %g %g\n",
 
1266
                    LASHeader_GetScaleX(LAS_header),
 
1267
                    LASHeader_GetScaleY(LAS_header),
 
1268
                    LASHeader_GetScaleZ(LAS_header));
 
1269
    fprintf(stdout, "Offset X Y Z:                      %g %g %g\n",
 
1270
                    LASHeader_GetOffsetX(LAS_header),
 
1271
                    LASHeader_GetOffsetY(LAS_header),
 
1272
                    LASHeader_GetOffsetZ(LAS_header));
 
1273
    fprintf(stdout, "Min X Y Z:                         %g %g %g\n",
 
1274
                    LASHeader_GetMinX(LAS_header),
 
1275
                    LASHeader_GetMinY(LAS_header),
 
1276
                    LASHeader_GetMinZ(LAS_header));
 
1277
    fprintf(stdout, "Max X Y Z:                         %g %g %g\n",
 
1278
                    LASHeader_GetMaxX(LAS_header),
 
1279
                    LASHeader_GetMaxY(LAS_header),
 
1280
                    LASHeader_GetMaxZ(LAS_header));
 
1281
    if (las_srs_proj4 && strlen(las_srs_proj4) > 0) {
 
1282
        fprintf(stdout, "Spatial Reference:\n");
 
1283
        fprintf(stdout, "%s\n", las_srs_proj4);
 
1284
    }
 
1285
    else {
 
1286
        fprintf(stdout, "Spatial Reference:                 None\n");
 
1287
    }
 
1288
    
 
1289
    fprintf(stdout, "\nData Fields:\n");
 
1290
    fprintf(stdout, "  'X'\n  'Y'\n  'Z'\n  'Intensity'\n  'Return Number'\n");
 
1291
    fprintf(stdout, "  'Number of Returns'\n  'Scan Direction'\n");
 
1292
    fprintf(stdout, "  'Flighline Edge'\n  'Classification'\n  'Scan Angle Rank'\n");
 
1293
    fprintf(stdout, "  'User Data'\n  'Point Source ID'\n");
 
1294
    if (las_point_format == 1 || las_point_format == 3 || las_point_format == 4 || las_point_format == 5) {
 
1295
        fprintf(stdout, "  'GPS Time'\n");
 
1296
    }
 
1297
    if (las_point_format == 2 || las_point_format == 3 || las_point_format == 5) {
 
1298
        fprintf(stdout, "  'Red'\n  'Green'\n  'Blue'\n");
 
1299
    }
 
1300
    fprintf(stdout, "\n");
 
1301
    fflush(stdout);
 
1302
 
 
1303
    return;
 
1304
}
 
1305
 
 
1306
 
 
1307
int scan_bounds(LASReaderH LAS_reader, int shell_style, int extents,
 
1308
                double zscale, struct Cell_head *region)
 
1309
{
 
1310
    unsigned long line;
 
1311
    int first;
 
1312
    double min_x, max_x, min_y, max_y, min_z, max_z;
 
1313
    double x, y, z;
 
1314
    LASPointH LAS_point;
 
1315
 
 
1316
    line = 0;
 
1317
    first = TRUE;
 
1318
    
 
1319
    /* init to nan in case no points are found */
 
1320
    min_x = max_x = min_y = max_y = min_z = max_z = 0.0 / 0.0;
 
1321
 
 
1322
    G_verbose_message(_("Scanning data ..."));
 
1323
    
 
1324
    LASReader_Seek(LAS_reader, 0);
 
1325
 
 
1326
    while ((LAS_point = LASReader_GetNextPoint(LAS_reader)) != NULL) {
 
1327
        line++;
 
1328
 
 
1329
        if (!LASPoint_IsValid(LAS_point)) {
 
1330
            continue;
 
1331
        }
 
1332
 
 
1333
        x = LASPoint_GetX(LAS_point);
 
1334
        y = LASPoint_GetY(LAS_point);
 
1335
        z = LASPoint_GetZ(LAS_point);
 
1336
 
 
1337
        if (first) {
 
1338
            min_x = x;
 
1339
            max_x = x;
 
1340
            min_y = y;
 
1341
            max_y = y;
 
1342
            min_z = z;
 
1343
            max_z = z;
 
1344
            first = FALSE;
 
1345
        }
 
1346
        else {
 
1347
            if (x < min_x)
 
1348
                min_x = x;
 
1349
            if (x > max_x)
 
1350
                max_x = x;
 
1351
            if (y < min_y)
 
1352
                min_y = y;
 
1353
            if (y > max_y)
 
1354
                max_y = y;
 
1355
            if (z < min_z)
 
1356
                min_z = z;
 
1357
            if (z > max_z)
 
1358
                max_z = z;
 
1359
        }
 
1360
    }
 
1361
 
 
1362
    if (!extents) {
 
1363
        if (!shell_style) {
 
1364
            fprintf(stderr, _("Range:     min         max\n"));
 
1365
            fprintf(stdout, "x: %11f %11f\n", min_x, max_x);
 
1366
            fprintf(stdout, "y: %11f %11f\n", min_y, max_y);
 
1367
            fprintf(stdout, "z: %11f %11f\n", min_z * zscale, max_z * zscale);
 
1368
        }
 
1369
        else
 
1370
            fprintf(stdout, "n=%f s=%f e=%f w=%f b=%f t=%f\n",
 
1371
                    max_y, min_y, max_x, min_x, min_z * zscale, max_z * zscale);
 
1372
 
 
1373
        G_debug(1, "Processed %lu points.", line);
 
1374
        G_debug(1, "region template: g.region n=%f s=%f e=%f w=%f",
 
1375
                max_y, min_y, max_x, min_x);
 
1376
    }
 
1377
    else {
 
1378
        region->east = max_x;
 
1379
        region->west = min_x;
 
1380
        region->north = max_y;
 
1381
        region->south = min_y;
 
1382
    }
 
1383
 
 
1384
    return 0;
 
1385
}