Changeset 66595


Ignore:
Timestamp:
Oct 24, 2015, 8:42:35 PM (9 years ago)
Author:
wenzeslaus
Message:

r.in.lidar: separate binning code, split to functions

Location:
grass/trunk/raster/r.in.lidar
Files:
2 added
1 edited

Legend:

Unmodified
Added
Removed
  • grass/trunk/raster/r.in.lidar/main.c

    r66594 r66595  
    3333#include "local_proto.h"
    3434#include "rast_segment.h"
    35 
    36 struct node
    37 {
    38     int next;
    39     double z;
    40 };
    41 
    42 #define SIZE_INCREMENT 10
    43 int num_nodes = 0;
    44 int max_nodes = 0;
    45 struct node *nodes;
     35#include "point_binning.h"
    4636
    4737#define LAS_ALL 0
     
    4939#define LAS_LAST 2
    5040#define LAS_MID 3
    51 
    52 int new_node(void)
    53 {
    54     int n = num_nodes++;
    55 
    56     if (num_nodes >= max_nodes) {
    57         max_nodes += SIZE_INCREMENT;
    58         nodes = G_realloc(nodes, (size_t)max_nodes * sizeof(struct node));
    59     }
    60 
    61     return n;
    62 }
    63 
    64 
    65 /* add node to sorted, single linked list
    66  * returns id if head has to be saved to index array, otherwise -1 */
    67 int add_node(int head, double z)
    68 {
    69     int node_id, last_id, newnode_id, head_id;
    70 
    71     head_id = head;
    72     node_id = head_id;
    73     last_id = head_id;
    74 
    75     while (node_id != -1 && nodes[node_id].z < z) {
    76         last_id = node_id;
    77         node_id = nodes[node_id].next;
    78     }
    79 
    80     /* end of list, simply append */
    81     if (node_id == -1) {
    82         newnode_id = new_node();
    83         nodes[newnode_id].next = -1;
    84         nodes[newnode_id].z = z;
    85         nodes[last_id].next = newnode_id;
    86         return -1;
    87     }
    88     else if (node_id == head_id) {      /* pole position, insert as head */
    89         newnode_id = new_node();
    90         nodes[newnode_id].next = head_id;
    91         head_id = newnode_id;
    92         nodes[newnode_id].z = z;
    93         return (head_id);
    94     }
    95     else {                      /* somewhere in the middle, insert */
    96         newnode_id = new_node();
    97         nodes[newnode_id].z = z;
    98         nodes[newnode_id].next = node_id;
    99         nodes[last_id].next = newnode_id;
    100         return -1;
    101     }
    102 }
    10341
    10442int main(int argc, char *argv[])
     
    10745    char *infile, *outmap;
    10846    int percent;
    109     int method = -1;
    110     int bin_n, bin_min, bin_max, bin_sum, bin_sumsq, bin_index;
    11147    double zrange_min, zrange_max, d_tmp;
    11248    unsigned long estimated_lines;
     
    11652    char title[64];
    11753    SEGMENT base_segment;
    118     void *n_array, *min_array, *max_array, *sum_array, *sumsq_array,
    119         *index_array, *base_array;
    120     void *raster_row, *ptr;
     54    struct PointBinning point_binning;
     55    void *base_array;
     56    void *raster_row;
    12157    struct Cell_head region;
    12258    struct Cell_head input_region;
    12359    int rows, cols;             /* scan box size */
    124     int row, col;               /* counters */
     60    int row;            /* counters */
    12561
    12662    int pass, npasses;
     
    13571    int point_class;
    13672
    137     double min = 0.0 / 0.0;     /* init as nan */
    138     double max = 0.0 / 0.0;     /* init as nan */
    13973    double zscale = 1.0;
    140     size_t offset, n_offset;
    141     int n = 0;
    142     double sum = 0.;
    143     double sumsq = 0.;
    144     double variance, mean, skew, sumdev;
    145     int pth = 0;
    146     double trim = 0.0;
    14774    double res = 0.0;
    14875
    149     int j, k;
    150     int head_id, node_id;
    151     int r_low, r_up;
     76    struct BinIndex bin_index_nodes;
     77    bin_index_nodes.num_nodes = 0;
     78    bin_index_nodes.max_nodes = 0;
     79    bin_index_nodes.nodes = 0;
    15280
    15381    struct GModule *module;
     
    397325    }
    398326
    399     /* figure out what maps we need in memory */
    400     /*  n               n
    401        min              min
    402        max              max
    403        range            min max         max - min
    404        sum              sum
    405        mean             sum n           sum/n
    406        stddev           sum sumsq n     sqrt((sumsq - sum*sum/n)/n)
    407        variance         sum sumsq n     (sumsq - sum*sum/n)/n
    408        coeff_var        sum sumsq n     sqrt((sumsq - sum*sum/n)/n) / (sum/n)
    409        median           n               array index to linked list
    410        percentile       n               array index to linked list
    411        skewness         n               array index to linked list
    412        trimmean         n               array index to linked list
    413      */
    414     bin_n = FALSE;
    415     bin_min = FALSE;
    416     bin_max = FALSE;
    417     bin_sum = FALSE;
    418     bin_sumsq = FALSE;
    419     bin_index = FALSE;
    420 
    421     n_array = NULL;
    422     min_array = NULL;
    423     max_array = NULL;
    424     sum_array = NULL;
    425     sumsq_array = NULL;
    426     index_array = NULL;
     327    point_binning_set(&point_binning, method_opt->answer, pth_opt->answer, trim_opt->answer);
     328
    427329    base_array = NULL;
    428    
    429     if (strcmp(method_opt->answer, "n") == 0) {
    430         method = METHOD_N;
    431         bin_n = TRUE;
    432     }
    433     if (strcmp(method_opt->answer, "min") == 0) {
    434         method = METHOD_MIN;
    435         bin_min = TRUE;
    436     }
    437     if (strcmp(method_opt->answer, "max") == 0) {
    438         method = METHOD_MAX;
    439         bin_max = TRUE;
    440     }
    441     if (strcmp(method_opt->answer, "range") == 0) {
    442         method = METHOD_RANGE;
    443         bin_min = TRUE;
    444         bin_max = TRUE;
    445     }
    446     if (strcmp(method_opt->answer, "sum") == 0) {
    447         method = METHOD_SUM;
    448         bin_sum = TRUE;
    449     }
    450     if (strcmp(method_opt->answer, "mean") == 0) {
    451         method = METHOD_MEAN;
    452         bin_sum = TRUE;
    453         bin_n = TRUE;
    454     }
    455     if (strcmp(method_opt->answer, "stddev") == 0) {
    456         method = METHOD_STDDEV;
    457         bin_sum = TRUE;
    458         bin_sumsq = TRUE;
    459         bin_n = TRUE;
    460     }
    461     if (strcmp(method_opt->answer, "variance") == 0) {
    462         method = METHOD_VARIANCE;
    463         bin_sum = TRUE;
    464         bin_sumsq = TRUE;
    465         bin_n = TRUE;
    466     }
    467     if (strcmp(method_opt->answer, "coeff_var") == 0) {
    468         method = METHOD_COEFF_VAR;
    469         bin_sum = TRUE;
    470         bin_sumsq = TRUE;
    471         bin_n = TRUE;
    472     }
    473     if (strcmp(method_opt->answer, "median") == 0) {
    474         method = METHOD_MEDIAN;
    475         bin_index = TRUE;
    476     }
    477     if (strcmp(method_opt->answer, "percentile") == 0) {
    478         if (pth_opt->answer != NULL)
    479             pth = atoi(pth_opt->answer);
    480         else
    481             G_fatal_error(_("Unable to calculate percentile without the pth option specified!"));
    482         method = METHOD_PERCENTILE;
    483         bin_index = TRUE;
    484     }
    485     if (strcmp(method_opt->answer, "skewness") == 0) {
    486         method = METHOD_SKEWNESS;
    487         bin_index = TRUE;
    488     }
    489     if (strcmp(method_opt->answer, "trimmean") == 0) {
    490         if (trim_opt->answer != NULL)
    491             trim = atof(trim_opt->answer) / 100.0;
    492         else
    493             G_fatal_error(_("Unable to calculate trimmed mean without the trim option specified!"));
    494         method = METHOD_TRIMMEAN;
    495         bin_index = TRUE;
    496     }
    497330
    498331    if (strcmp("CELL", type_opt->answer) == 0)
     
    503336        rtype = FCELL_TYPE;
    504337
    505     if (method == METHOD_N)
     338    if (point_binning.method == METHOD_N)
    506339        rtype = CELL_TYPE;
    507340
     
    585418
    586419    if (!scan_flag->answer) {
    587         /* check if rows * (cols + 1) go into a size_t */
    588         if (sizeof(size_t) < 8) {
    589             double dsize = rows * (cols + 1);
    590            
    591             if (dsize != (size_t)rows * (cols + 1))
     420        if (!check_rows_cols_fit_to_size_t(rows, cols))
    592421                G_fatal_error(_("Unable to process the hole map at once. "
    593                                 "Please set the %s option to some value lower than 100."),
     422                        "Please set the '%s' option to some value lower than 100."),
    594423                                percent_opt->key);
     424        point_binning_memory_test(&point_binning, rows, cols, rtype);
    595425        }
    596         /* allocate memory (test for enough before we start) */
    597         if (bin_n)
    598             n_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(CELL_TYPE));
    599         if (bin_min)
    600             min_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
    601         if (bin_max)
    602             max_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
    603         if (bin_sum)
    604             sum_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
    605         if (bin_sumsq)
    606             sumsq_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
    607         if (bin_index)
    608             index_array =
    609                 G_calloc((size_t)rows * (cols + 1), Rast_cell_size(CELL_TYPE));
    610         /* TODO: perhaps none of them needs to be freed */
    611 
    612         /* and then free it again */
    613         if (bin_n)
    614             G_free(n_array);
    615         if (bin_min)
    616             G_free(min_array);
    617         if (bin_max)
    618             G_free(max_array);
    619         if (bin_sum)
    620             G_free(sum_array);
    621         if (bin_sumsq)
    622             G_free(sumsq_array);
    623         if (bin_index)
    624             G_free(index_array);
    625 
    626         /** end memory test **/
    627     }
    628 
    629426
    630427    /* open output map */
     
    673470                pass, npasses, pass_north, pass_south, rows);
    674471
    675 
    676         if (bin_n) {
    677             G_debug(2, "allocating n_array");
    678             n_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(CELL_TYPE));
    679             blank_array(n_array, rows, cols, CELL_TYPE, 0);
    680         }
    681         if (bin_min) {
    682             G_debug(2, "allocating min_array");
    683             min_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
    684             blank_array(min_array, rows, cols, rtype, -1);      /* fill with NULLs */
    685         }
    686         if (bin_max) {
    687             G_debug(2, "allocating max_array");
    688             max_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
    689             blank_array(max_array, rows, cols, rtype, -1);      /* fill with NULLs */
    690         }
    691         if (bin_sum) {
    692             G_debug(2, "allocating sum_array");
    693             sum_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
    694             blank_array(sum_array, rows, cols, rtype, 0);
    695         }
    696         if (bin_sumsq) {
    697             G_debug(2, "allocating sumsq_array");
    698             sumsq_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
    699             blank_array(sumsq_array, rows, cols, rtype, 0);
    700         }
    701         if (bin_index) {
    702             G_debug(2, "allocating index_array");
    703             index_array =
    704                 G_calloc((size_t)rows * (cols + 1), Rast_cell_size(CELL_TYPE));
    705             blank_array(index_array, rows, cols, CELL_TYPE, -1);        /* fill with NULLs */
    706         }
     472    point_binning_allocate(&point_binning, rows, cols, rtype);
    707473
    708474        line = 0;
     
    815581            /*          G_debug(5, "x: %f, y: %f, z: %f", x, y, z); */
    816582
    817             if (bin_n)
    818                 update_n(n_array, cols, arr_row, arr_col);
    819             if (bin_min)
    820                 update_min(min_array, cols, arr_row, arr_col, rtype, z);
    821             if (bin_max)
    822                 update_max(max_array, cols, arr_row, arr_col, rtype, z);
    823             if (bin_sum)
    824                 update_sum(sum_array, cols, arr_row, arr_col, rtype, z);
    825             if (bin_sumsq)
    826                 update_sumsq(sumsq_array, cols, arr_row, arr_col, rtype, z);
    827             if (bin_index) {
    828                 ptr = index_array;
    829                 ptr =
    830                     G_incr_void_ptr(ptr,
    831                                     ((arr_row * cols) +
    832                                      arr_col) * Rast_cell_size(CELL_TYPE));
    833 
    834                 if (Rast_is_null_value(ptr, CELL_TYPE)) {       /* first node */
    835                     head_id = new_node();
    836                     nodes[head_id].next = -1;
    837                     nodes[head_id].z = z;
    838                     Rast_set_c_value(ptr, head_id, CELL_TYPE);  /* store index to head */
    839                 }
    840                 else {          /* head is already there */
    841 
    842                     head_id = Rast_get_c_value(ptr, CELL_TYPE); /* get index to head */
    843                     head_id = add_node(head_id, z);
    844                     if (head_id != -1)
    845                         Rast_set_c_value(ptr, head_id, CELL_TYPE);      /* store index to head */
    846                 }
    847             }
     583            update_value(&point_binning, &bin_index_nodes, cols, arr_row, arr_col, rtype, z);
    848584        }                       /* while !EOF */
    849585
     
    856592        G_message(_("Writing to map ..."));
    857593        for (row = 0; row < rows; row++) {
    858 
    859             switch (method) {
    860             case METHOD_N:      /* n is a straight copy */
    861                 Rast_raster_cpy(raster_row,
    862                              n_array +
    863                              (row * cols * Rast_cell_size(CELL_TYPE)), cols,
    864                              CELL_TYPE);
    865                 break;
    866 
    867             case METHOD_MIN:
    868                 Rast_raster_cpy(raster_row,
    869                              min_array + (row * cols * Rast_cell_size(rtype)),
    870                              cols, rtype);
    871                 break;
    872 
    873             case METHOD_MAX:
    874                 Rast_raster_cpy(raster_row,
    875                              max_array + (row * cols * Rast_cell_size(rtype)),
    876                              cols, rtype);
    877                 break;
    878 
    879             case METHOD_SUM:
    880                 Rast_raster_cpy(raster_row,
    881                              sum_array + (row * cols * Rast_cell_size(rtype)),
    882                              cols, rtype);
    883                 break;
    884 
    885             case METHOD_RANGE:  /* (max-min) */
    886                 ptr = raster_row;
    887                 for (col = 0; col < cols; col++) {
    888                     offset = (row * cols + col) * Rast_cell_size(rtype);
    889                     min = Rast_get_d_value(min_array + offset, rtype);
    890                     max = Rast_get_d_value(max_array + offset, rtype);
    891                     Rast_set_d_value(ptr, max - min, rtype);
    892                     ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
    893                 }
    894                 break;
    895 
    896             case METHOD_MEAN:   /* (sum / n) */
    897                 ptr = raster_row;
    898                 for (col = 0; col < cols; col++) {
    899                     offset = (row * cols + col) * Rast_cell_size(rtype);
    900                     n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
    901                     n = Rast_get_c_value(n_array + n_offset, CELL_TYPE);
    902                     sum = Rast_get_d_value(sum_array + offset, rtype);
    903 
    904                     if (n == 0)
    905                         Rast_set_null_value(ptr, 1, rtype);
    906                     else
    907                         Rast_set_d_value(ptr, (sum / n), rtype);
    908 
    909                     ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
    910                 }
    911                 break;
    912 
    913             case METHOD_STDDEV: /*  sqrt(variance)        */
    914             case METHOD_VARIANCE:       /*  (sumsq - sum*sum/n)/n */
    915             case METHOD_COEFF_VAR:      /*  100 * stdev / mean    */
    916                 ptr = raster_row;
    917                 for (col = 0; col < cols; col++) {
    918                     offset = (row * cols + col) * Rast_cell_size(rtype);
    919                     n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
    920                     n = Rast_get_c_value(n_array + n_offset, CELL_TYPE);
    921                     sum = Rast_get_d_value(sum_array + offset, rtype);
    922                     sumsq = Rast_get_d_value(sumsq_array + offset, rtype);
    923 
    924                     if (n == 0)
    925                         Rast_set_null_value(ptr, 1, rtype);
    926                     else if (n == 1)
    927                         Rast_set_d_value(ptr, 0.0, rtype);
    928                     else {
    929                         variance = (sumsq - sum * sum / n) / n;
    930                         if (variance < GRASS_EPSILON)
    931                             variance = 0.0;
    932 
    933                         /* nan test */
    934                         if (variance != variance)
    935                             Rast_set_null_value(ptr, 1, rtype);
    936                         else {
    937 
    938                             if (method == METHOD_STDDEV)
    939                                 variance = sqrt(variance);
    940 
    941                             else if (method == METHOD_COEFF_VAR)
    942                                 variance = 100 * sqrt(variance) / (sum / n);
    943 
    944                             /* nan test */
    945                             if (variance != variance)
    946                                 variance = 0.0; /* OK for n > 0 ?*/
    947 
    948                             Rast_set_d_value(ptr, variance, rtype);
    949                         }
    950 
    951                     }
    952                     ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
    953                 }
    954                 break;
    955             case METHOD_MEDIAN: /* median, if only one point in cell we will use that */
    956                 ptr = raster_row;
    957                 for (col = 0; col < cols; col++) {
    958                     n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
    959                     if (Rast_is_null_value(index_array + n_offset, CELL_TYPE))  /* no points in cell */
    960                         Rast_set_null_value(ptr, 1, rtype);
    961                     else {      /* one or more points in cell */
    962 
    963                         head_id =
    964                             Rast_get_c_value(index_array + n_offset,
    965                                                  CELL_TYPE);
    966                         node_id = head_id;
    967 
    968                         n = 0;
    969 
    970                         while (node_id != -1) { /* count number of points in cell */
    971                             n++;
    972                             node_id = nodes[node_id].next;
    973                         }
    974 
    975                         if (n == 1)     /* only one point, use that */
    976                             Rast_set_d_value(ptr, nodes[head_id].z,
    977                                                  rtype);
    978                         else if (n % 2 != 0) {  /* odd number of points: median_i = (n + 1) / 2 */
    979                             n = (n + 1) / 2;
    980                             node_id = head_id;
    981                             for (j = 1; j < n; j++)     /* get "median element" */
    982                                 node_id = nodes[node_id].next;
    983 
    984                             Rast_set_d_value(ptr, nodes[node_id].z,
    985                                                  rtype);
    986                         }
    987                         else {  /* even number of points: median = (val_below + val_above) / 2 */
    988 
    989                             z = (n + 1) / 2.0;
    990                             n = floor(z);
    991                             node_id = head_id;
    992                             for (j = 1; j < n; j++)     /* get element "below" */
    993                                 node_id = nodes[node_id].next;
    994 
    995                             z = (nodes[node_id].z +
    996                                  nodes[nodes[node_id].next].z) / 2;
    997                             Rast_set_d_value(ptr, z, rtype);
    998                         }
    999                     }
    1000                     ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
    1001                 }
    1002                 break;
    1003             case METHOD_PERCENTILE:     /* rank = (pth*(n+1))/100; interpolate linearly */
    1004                 ptr = raster_row;
    1005                 for (col = 0; col < cols; col++) {
    1006                     n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
    1007                     if (Rast_is_null_value(index_array + n_offset, CELL_TYPE))  /* no points in cell */
    1008                         Rast_set_null_value(ptr, 1, rtype);
    1009                     else {
    1010                         head_id =
    1011                             Rast_get_c_value(index_array + n_offset,
    1012                                                  CELL_TYPE);
    1013                         node_id = head_id;
    1014                         n = 0;
    1015 
    1016                         while (node_id != -1) { /* count number of points in cell */
    1017                             n++;
    1018                             node_id = nodes[node_id].next;
    1019                         }
    1020 
    1021                         z = (pth * (n + 1)) / 100.0;
    1022                         r_low = floor(z);       /* lower rank */
    1023                         if (r_low < 1)
    1024                             r_low = 1;
    1025                         else if (r_low > n)
    1026                             r_low = n;
    1027 
    1028                         r_up = ceil(z); /* upper rank */
    1029                         if (r_up > n)
    1030                             r_up = n;
    1031 
    1032                         node_id = head_id;
    1033                         for (j = 1; j < r_low; j++)     /* search lower value */
    1034                             node_id = nodes[node_id].next;
    1035 
    1036                         z = nodes[node_id].z;   /* save lower value */
    1037                         node_id = head_id;
    1038                         for (j = 1; j < r_up; j++)      /* search upper value */
    1039                             node_id = nodes[node_id].next;
    1040 
    1041                         z = (z + nodes[node_id].z) / 2;
    1042                         Rast_set_d_value(ptr, z, rtype);
    1043                     }
    1044                     ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
    1045                 }
    1046                 break;
    1047             case METHOD_SKEWNESS:       /* skewness = sum(xi-mean)^3/(N-1)*s^3 */
    1048                 ptr = raster_row;
    1049                 for (col = 0; col < cols; col++) {
    1050                     n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
    1051                     if (Rast_is_null_value(index_array + n_offset, CELL_TYPE))  /* no points in cell */
    1052                         Rast_set_null_value(ptr, 1, rtype);
    1053                     else {
    1054                         head_id =
    1055                             Rast_get_c_value(index_array + n_offset,
    1056                                                  CELL_TYPE);
    1057                         node_id = head_id;
    1058 
    1059                         n = 0;  /* count */
    1060                         sum = 0.0;      /* sum */
    1061                         sumsq = 0.0;    /* sum of squares */
    1062                         sumdev = 0.0;   /* sum of (xi - mean)^3 */
    1063                         skew = 0.0;     /* skewness */
    1064 
    1065                         while (node_id != -1) {
    1066                             z = nodes[node_id].z;
    1067                             n++;
    1068                             sum += z;
    1069                             sumsq += (z * z);
    1070                             node_id = nodes[node_id].next;
    1071                         }
    1072 
    1073                         if (n > 1) {    /* if n == 1, skew is "0.0" */
    1074                             mean = sum / n;
    1075                             node_id = head_id;
    1076                             while (node_id != -1) {
    1077                                 z = nodes[node_id].z;
    1078                                 sumdev += pow((z - mean), 3);
    1079                                 node_id = nodes[node_id].next;
    1080                             }
    1081 
    1082                             variance = (sumsq - sum * sum / n) / n;
    1083                             if (variance < GRASS_EPSILON)
    1084                                 skew = 0.0;
    1085                             else
    1086                                 skew =
    1087                                     sumdev / ((n - 1) *
    1088                                               pow(sqrt(variance), 3));
    1089                         }
    1090                         Rast_set_d_value(ptr, skew, rtype);
    1091                     }
    1092                     ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
    1093                 }
    1094                 break;
    1095             case METHOD_TRIMMEAN:
    1096                 ptr = raster_row;
    1097                 for (col = 0; col < cols; col++) {
    1098                     n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
    1099                     if (Rast_is_null_value(index_array + n_offset, CELL_TYPE))  /* no points in cell */
    1100                         Rast_set_null_value(ptr, 1, rtype);
    1101                     else {
    1102                         head_id =
    1103                             Rast_get_c_value(index_array + n_offset,
    1104                                                  CELL_TYPE);
    1105 
    1106                         node_id = head_id;
    1107                         n = 0;
    1108                         while (node_id != -1) { /* count number of points in cell */
    1109                             n++;
    1110                             node_id = nodes[node_id].next;
    1111                         }
    1112 
    1113                         if (1 == n)
    1114                             mean = nodes[head_id].z;
    1115                         else {
    1116                             k = floor(trim * n + 0.5);  /* number of ranks to discard on each tail */
    1117 
    1118                             if (k > 0 && (n - 2 * k) > 0) {     /* enough elements to discard */
    1119                                 node_id = head_id;
    1120                                 for (j = 0; j < k; j++) /* move to first rank to consider */
    1121                                     node_id = nodes[node_id].next;
    1122 
    1123                                 j = k + 1;
    1124                                 k = n - k;
    1125                                 n = 0;
    1126                                 sum = 0.0;
    1127 
    1128                                 while (j <= k) {        /* get values in interval */
    1129                                     n++;
    1130                                     sum += nodes[node_id].z;
    1131                                     node_id = nodes[node_id].next;
    1132                                     j++;
    1133                                 }
    1134                             }
    1135                             else {
    1136                                 node_id = head_id;
    1137                                 n = 0;
    1138                                 sum = 0.0;
    1139                                 while (node_id != -1) {
    1140                                     n++;
    1141                                     sum += nodes[node_id].z;
    1142                                     node_id = nodes[node_id].next;
    1143                                 }
    1144                             }
    1145                             mean = sum / n;
    1146                         }
    1147                         Rast_set_d_value(ptr, mean, rtype);
    1148                     }
    1149                     ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
    1150                 }
    1151                 break;
    1152 
    1153             default:
    1154                 G_fatal_error("?");
    1155             }
    1156 
     594        write_values(&point_binning, &bin_index_nodes, raster_row, row, cols, rtype);
    1157595            /* write out line of raster data */
    1158596            Rast_put_row(out_fd, raster_row, rtype);
     
    1160598
    1161599        /* free memory */
    1162         if (bin_n)
    1163             G_free(n_array);
    1164         if (bin_min)
    1165             G_free(min_array);
    1166         if (bin_max)
    1167             G_free(max_array);
    1168         if (bin_sum)
    1169             G_free(sum_array);
    1170         if (bin_sumsq)
    1171             G_free(sumsq_array);
    1172         if (bin_index) {
    1173             G_free(index_array);
    1174             G_free(nodes);
    1175             num_nodes = 0;
    1176             max_nodes = 0;
    1177             nodes = NULL;
    1178         }
     600        point_binning_free(&point_binning, &bin_index_nodes);
    1179601    }                           /* passes loop */
    1180602    if (base_array)
Note: See TracChangeset for help on using the changeset viewer.