[GRASS-SVN] r61373 - in grass/branches/releasebranch_7_0: . raster3d/r3.stats
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Jul 23 06:22:59 PDT 2014
Author: neteler
Date: 2014-07-23 06:22:59 -0700 (Wed, 23 Jul 2014)
New Revision: 61373
Added:
grass/branches/releasebranch_7_0/raster3d/r3.stats/local_proto.h
grass/branches/releasebranch_7_0/raster3d/r3.stats/support.c
Modified:
grass/branches/releasebranch_7_0/
grass/branches/releasebranch_7_0/raster3d/r3.stats/
grass/branches/releasebranch_7_0/raster3d/r3.stats/main.c
Log:
r3.stats: split main.c up into bite sized chunks (trunk, r60100); percent cosmetics
Property changes on: grass/branches/releasebranch_7_0
___________________________________________________________________
Modified: svn:mergeinfo
- /grass/trunk:59505,59646,60215,60219,60278,60376,60610,60807,60835,60936-60937,61159-61160,61165,61275,61288,61290,61292,61294,61301
+ /grass/trunk:59505,59646,60100,60215,60219,60278,60376,60610,60807,60835,60936-60937,61159-61160,61165,61275,61288,61290,61292,61294,61301
Property changes on: grass/branches/releasebranch_7_0/raster3d/r3.stats
___________________________________________________________________
Modified: svn:ignore
- OBJ.*
+ OBJ.*
*.tmp.html
Copied: grass/branches/releasebranch_7_0/raster3d/r3.stats/local_proto.h (from rev 60100, grass/trunk/raster3d/r3.stats/local_proto.h)
===================================================================
--- grass/branches/releasebranch_7_0/raster3d/r3.stats/local_proto.h (rev 0)
+++ grass/branches/releasebranch_7_0/raster3d/r3.stats/local_proto.h 2014-07-23 13:22:59 UTC (rev 61373)
@@ -0,0 +1,49 @@
+#include <grass/raster3d.h>
+
+#define COMPARE_PRECISION 0.000000001
+
+/* statistic table row struct */
+typedef struct
+{
+ double min, max, vol, perc;
+ int num, count;
+} stat_row;
+
+/* the statistic table struct */
+typedef struct
+{
+ stat_row **table;
+ stat_row *null;
+ int sum_count, nsteps, equal;
+ double sum_vol, sum_perc;
+} stat_table;
+
+/* the equal value struct */
+typedef struct
+{
+ double val; /* the equal value */
+ int count; /* the appearance count */
+} equal_val;
+
+/* an array of groups with equal values */
+typedef struct
+{
+ equal_val **values;
+ int count;
+} equal_val_array;
+
+/* prototypes */
+equal_val_array *alloc_equal_val_array(int count);
+void free_equal_val_array(equal_val_array *vals);
+equal_val_array *add_equal_val_to_array(equal_val_array *array, double val);
+int check_equal_value(equal_val_array * array, double val);
+stat_table *create_stat_table(int nsteps, equal_val_array *values,
+ double min, double max);
+void free_stat_table(stat_table *stats);
+void print_stat_table(stat_table *stats);
+void update_stat_table(stat_table *stats, RASTER3D_Region *region);
+void heapsort_eqvals(equal_val_array *data, int n);
+void downheap_eqvals(equal_val_array *data, int n, int k);
+void check_range_value(stat_table *stats, double value);
+void tree_search_range(stat_table *stats, int left, int right,
+ double value);
Modified: grass/branches/releasebranch_7_0/raster3d/r3.stats/main.c
===================================================================
--- grass/branches/releasebranch_7_0/raster3d/r3.stats/main.c 2014-07-23 13:01:33 UTC (rev 61372)
+++ grass/branches/releasebranch_7_0/raster3d/r3.stats/main.c 2014-07-23 13:22:59 UTC (rev 61373)
@@ -18,529 +18,9 @@
#include <grass/gis.h>
#include <grass/raster3d.h>
#include <grass/glocale.h>
+#include "local_proto.h"
-#define COMPARE_PRECISION 0.000000001
-/*statistic table row struct */
-typedef struct
-{
- double min, max, vol, perc;
- int num, count;
-} stat_row;
-
-/*the statistic table struct */
-typedef struct
-{
- stat_row **table;
- stat_row *null;
- int sum_count, nsteps, equal;
- double sum_vol, sum_perc;
-} stat_table;
-
-/*the equal value struct */
-typedef struct
-{
- double val; /* the equal value */
- int count; /* the appearance count */
-} equal_val;
-
-/*an array of groups with equal values */
-typedef struct
-{
- equal_val **values;
- int count;
-} equal_val_array;
-
-/*prototypes */
-static equal_val_array *alloc_equal_val_array(int count);
-static void free_equal_val_array(equal_val_array * vals);
-static equal_val_array *add_equal_val_to_array(equal_val_array * array, double val);
-static int check_equal_value(equal_val_array * array, double val);
-static stat_table *create_stat_table(int nsteps, equal_val_array * values,
- double min, double max);
-static void free_stat_table(stat_table * stats);
-static void print_stat_table(stat_table * stats);
-static void update_stat_table(stat_table * stats, RASTER3D_Region * region);
-static void heapsort_eqvals(equal_val_array * data, int n);
-static void downheap_eqvals(equal_val_array * data, int n, int k);
-static void check_range_value(stat_table * stats, double value);
-static void tree_search_range(stat_table * stats, int left, int right,
- double value);
-
-/* *************************************************************** */
-/* ***** allocate equal_val_array structure ********************* */
-/* *************************************************************** */
-equal_val_array *alloc_equal_val_array(int count)
-{
- equal_val_array *p;
- int i;
-
-
- p = (equal_val_array *) G_calloc(1, sizeof(equal_val_array));
- p->count = count;
-
- p->values = (equal_val **) G_calloc(p->count, sizeof(equal_val *));
-
- for (i = 0; i < p->count; i++)
- p->values[i] = (equal_val *) G_calloc(1, sizeof(equal_val));
-
-
- return p;
-}
-
-/* *************************************************************** */
-/* ***** add an equal_val to the equal_val_array structure ******* */
-/* *************************************************************** */
-equal_val_array *add_equal_val_to_array(equal_val_array * array, double val)
-{
- equal_val_array *p = array;
- int count;
-
- if (p == NULL) {
- p = alloc_equal_val_array(1);
- p->values[0]->val = val;
- p->values[0]->count = 1;
- G_debug(5, "Create new equal_array with value %g\n", val);
- }
- else {
- /*increase the count */
- count = array->count;
- count++;
-
- /*new memory */
- p->values =
- (equal_val **) G_realloc(p->values, count * sizeof(equal_val *));
- p->values[count - 1] = (equal_val *) G_calloc(1, sizeof(equal_val));
- /*set the new value */
- p->values[count - 1]->val = val;
- p->values[count - 1]->count = 1;
- /*set the new counter */
- p->count = count;
- G_debug(5, "Add new value %g at position %i\n", val, p->count);
- }
-
- return p;
-}
-
-/* *************************************************************** */
-/* ***** check if a value exists in the equal_val_array ********* */
-/* *************************************************************** */
-int check_equal_value(equal_val_array * array, double val)
-{
- int i;
-
- /*search if the new value exists and increase the counter */
- if (array != NULL)
- for (i = 0; i < array->count; i++) {
- if (array->values[i]->val == val) {
- array->values[i]->count++;
- G_debug(5, "found value %g with count %i at pos %i\n", val,
- array->values[i]->count, i);
- return 1;
- }
- }
-
- /*if it does not exists, add it to the array */
- add_equal_val_to_array(array, val);
-
- /*return 1 if found, 0 otherwise */
- return 0;
-}
-
-
-/* *************************************************************** */
-/* ***** Release the memory of a equal_val_array structure ****** */
-/* *************************************************************** */
-void free_equal_val_array(equal_val_array * uvals)
-{
- int i;
-
- for (i = 0; i < uvals->count; i++) {
- G_free(uvals->values[i]);
- }
-
- G_free(uvals->values);
- G_free(uvals);
- uvals = NULL;
-
- return;
-}
-
-/* *************************************************************** */
-/* **** Create the structure which manages the statistical ******* */
-/* **** values for a value range or equal values **************** */
-/* *************************************************************** */
-stat_table *create_stat_table(int nsteps, equal_val_array * eqvals,
- double min, double max)
-{
- stat_table *p;
- int i;
- double step;
-
- /* Memory */
- p = (stat_table *) G_calloc(1, sizeof(stat_table));
-
- p->null = (stat_row *) G_calloc(nsteps, sizeof(stat_row));
- p->table = (stat_row **) G_calloc(nsteps, sizeof(stat_row *));
-
- for (i = 0; i < nsteps; i++)
- p->table[i] = (stat_row *) G_calloc(1, sizeof(stat_row));
-
- /* Some value initializing */
- p->null->min = 0;
- p->null->max = 0;
- p->null->vol = 0;
- p->null->perc = 0;
- p->null->count = 0;
- p->null->num = nsteps + 1;
-
- p->nsteps = nsteps;
- p->sum_count = 0;
- p->sum_vol = 0;
- p->sum_perc = 0;
- p->equal = 0;
-
- /*if no equal values are provided, calculate the range and the length of each step */
- if (!eqvals) {
-
- /*calculate the steplength */
- step = (max - min) / (nsteps);
- p->equal = 0;
-
- p->table[0]->min = min;
- p->table[0]->max = min + step;
- p->table[0]->num = 1;
- p->table[0]->count = 0;
- p->table[0]->vol = 0;
- p->table[0]->perc = 0;
-
- G_debug(3, "Step %i range min %.11lf max %.11lf\n", p->table[0]->num,
- p->table[0]->min, p->table[0]->max);
-
- for (i = 1; i < nsteps; i++) {
- p->table[i]->min = p->table[i - 1]->max;
- p->table[i]->max = p->table[i - 1]->max + step;
- p->table[i]->num = i + 1;
- p->table[i]->count = 0;
- p->table[i]->vol = 0;
- p->table[i]->perc = 0;
- G_debug(5, "Step %i range min %.11lf max %.11lf\n",
- p->table[i]->num, p->table[i]->min, p->table[i]->max);
- }
- /* the last value must be a bit larger */
- p->table[nsteps - 1]->max += COMPARE_PRECISION;
- }
- else { /* Create equal value statistic */
- p->equal = 1;
- for (i = 0; i < eqvals->count; i++) {
- p->table[i]->min = eqvals->values[i]->val; /* equal values have no range, set min and max to the same value */
- p->table[i]->max = eqvals->values[i]->val;
- p->table[i]->num = i + 1;
- p->table[i]->count = eqvals->values[i]->count; /* the appearance count for each equal value */
- p->table[i]->vol = 0;
- p->table[i]->perc = 0;
- G_debug(5, "Unique value %i = %g count %i\n", p->table[i]->num,
- p->table[i]->min, p->table[i]->count);
- }
-
- }
-
- return p;
-}
-
-/* *************************************************************** */
-/* ***** Release the memory of a stat_table structure ************ */
-/* *************************************************************** */
-void free_stat_table(stat_table * stats)
-{
- int i;
-
- for (i = 0; i < stats->nsteps; i++) {
- G_free(stats->table[i]);
- }
-
- G_free(stats->table);
- G_free(stats->null);
-
- G_free(stats);
- stats = NULL;
-
- return;
-}
-
-/* *************************************************************** */
-/* **** Compute the volume, percentage and sums ****************** */
-/* *************************************************************** */
-void update_stat_table(stat_table * stats, RASTER3D_Region * region)
-{
- int i;
- double vol = region->ns_res * region->ew_res * region->tb_res;
- int cellnum = region->rows * region->cols * region->depths;
-
- /*Calculate volume and percentage for each range or equal value and the sum stats */
- for (i = 0; i < stats->nsteps; i++) {
- stats->table[i]->vol = stats->table[i]->count * vol;
- stats->table[i]->perc =
- (double)(100.0 * (double)stats->table[i]->count /
- (double)cellnum);
- stats->sum_count += stats->table[i]->count;
- stats->sum_vol += stats->table[i]->vol;
- stats->sum_perc += stats->table[i]->perc;
- }
-
- /*calculate the null value stats */
- stats->null->vol = stats->null->count * vol;
- stats->null->perc =
- (double)(100.0 * (double)stats->null->count / (double)cellnum);
-
- return;
-}
-
-/* *************************************************************** */
-/* *************************************************************** */
-/* *************************************************************** */
-void print_stat_table(stat_table * stats)
-{
- int i;
-
- if (stats->equal) {
- /* 1234567 012345678901234567 0123456789012 0123456 0123456789 */
- fprintf(stdout,
- " num | value | volume | perc | cell count\n");
-
- for (i = 0; i < stats->nsteps; i++) {
- fprintf(stdout, "%7i %18.6lf %13.3lf %7.5lf %10i\n",
- stats->table[i]->num, stats->table[i]->min,
- stats->table[i]->vol, stats->table[i]->perc,
- stats->table[i]->count);
- }
- fprintf(stdout,
- "%7i * %13.3lf %7.5lf %10i\n",
- stats->null->num, stats->null->vol, stats->null->perc,
- stats->null->count);
-
- fprintf(stdout, "\nNumber of groups with equal values: %i",
- stats->nsteps);
- }
- else {
- /* 1234567 012345678901234567 012345678901234567 0123456789012 0123456 0123456789 */
- fprintf(stdout,
- " num | minimum <= value | value < maximum | volume | perc | cell count\n");
-
- for (i = 0; i < stats->nsteps; i++) {
- fprintf(stdout,
- "%7i %18.9lf %18.9lf %13.3lf %7.5lf %10i\n",
- stats->table[i]->num, stats->table[i]->min,
- stats->table[i]->max, stats->table[i]->vol,
- stats->table[i]->perc, stats->table[i]->count);
-
- }
- fprintf(stdout,
- "%7i * * %13.3lf %7.5lf %10i\n",
- stats->null->num, stats->null->vol, stats->null->perc,
- stats->null->count);
- }
-
- fprintf(stdout,
- "\nSum of non Null cells: \n\tVolume = %13.3lf \n\tPercentage = %7.3lf \n\tCell count = %i\n",
- stats->sum_vol, stats->sum_perc, stats->sum_count);
-
- fprintf(stdout,
- "\nSum of all cells: \n\tVolume = %13.3lf \n\tPercentage = %7.3lf \n\tCell count = %i\n",
- stats->sum_vol + stats->null->vol,
- stats->sum_perc + stats->null->perc,
- stats->sum_count + stats->null->count);
-
-
- return;
-}
-
-/* *************************************************************** */
-/* Make an entry in the statistic table based on range value check */
-/* *************************************************************** */
-void check_range_value(stat_table * stats, double value)
-{
-
- /* this is very expensive for large range arrays! */
- /*
- * int i;
- * for (i = 0; i < stats->nsteps; i++) {
- * if (value >= stats->table[i]->min && value < stats->table[i]->max) {
- * stats->table[i]->count++;
- * }
- * }
- */
-
- /* the much faster tree search is used */
- tree_search_range(stats, 0, stats->nsteps - 1, value);
-
- return;
-}
-
-/* *************************************************************** */
-/* tree search method developed by Soeren Gebbert *************** */
-/* *************************************************************** */
-/*
- * This is a divide and conquer tree search algorithm
- *
- * The algorithm counts how many values are located in each range.
- * It divides therefor the range array in smaller pices.
- *
- * e.g:
- * 0 1 2 3 4 5 6 7 8
- * [[0,1],[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9]]
- * / \
- * [[0,1],[1,2],[2,3],[3,4]] [[4,5],[5,6],[6,7],[7,8],[8,9]]
- * / \ / \
- * [[0,1],[1,2]] [[2,3],[3,4]] [[4,5],[5,6]] [[6,7],[7,8],[8,9]]
- * / \ / \ / \ / \
- * [[0,1]] [[1,2]] [[2,3]] [[3,4]] [[4,5]] [[5,6]] [[6,7]] [[7,8],[8,9]]
- * / \
- * [[7,8]] [[8,9]]
- *
- * Searching the value 5.5 will result in the follwoing steps:
- *
- * value 5.5 in range of
- * 1.) left array index = 0 - right array index = 8 -> range [0 - 9]
- * 2.) left array index = 4 - right array index = 8 -> range [4 - 9]
- * 3.) left array index = 4 - right array index = 5 -> range [4 - 6]
- * 4.) the value is located in range array 5 with a range of [5 - 6]
- * 5.) now increase the value range counter of range array with index 5
- *
- * */
-void tree_search_range(stat_table * stats, int left, int right, double value)
-{
- int size = right - left;
-
- G_debug(5,
- "Search value %g in array size %i left border index %i right border index %i\n",
- value, size, left, right);
-
- /*if left and right are equal */
- if (size == 0) {
- stats->table[left]->count++;
- return;
- }
- else if (size == 1) { /* if the size is one, check directly */
-
- if (value >= stats->table[left]->min &&
- value < stats->table[left]->max) {
- stats->table[left]->count++;
- }
- else if (value >= stats->table[right]->min &&
- value < stats->table[right]->max) {
- stats->table[right]->count++;
- }
- return;
- }
- else {
- if ((size % 2 == 0)) { /*even */
- /*left side */
- right = left + (size) / 2;
-
- if (value >= stats->table[left]->min &&
- value < stats->table[right]->max) {
- tree_search_range(stats, left, right, value);
- return;
- }
-
- /*right side */
- left += (size) / 2;
- right = left + (size) / 2;
-
- if (value >= stats->table[left]->min &&
- value < stats->table[right]->max) {
- tree_search_range(stats, left, right, value);
- return;
- }
- }
- else { /*odd */
- /*left side */
- right = left + (size - 1) / 2;
-
- if (value >= stats->table[left]->min &&
- value < stats->table[right]->max) {
- tree_search_range(stats, left, right, value);
- return;
- }
-
- /*right side */
- left += (size - 1) / 2;
- right = left + (size - 1) / 2 + 1; /*the right array is one col larger */
-
- if (value >= stats->table[left]->min &&
- value < stats->table[right]->max) {
- tree_search_range(stats, left, right, value);
- return;
- }
- }
- }
-
- return;
-}
-
-/* *************************************************************** */
-/* ****** heapsort for equal value arrays of size n ************** */
-/* ** Code based on heapsort from Sebastian Cyris **************** */
-/* *************************************************************** */
-void heapsort_eqvals(equal_val_array * e, int n)
-{
- int k;
- double t;
- int count = 0;
-
- --n;
-
- for (k = n / 2; k >= 0; k--)
- downheap_eqvals(e, n, k);
-
- while (n > 0) {
- t = e->values[0]->val;
- count = e->values[0]->count;
- e->values[0]->val = e->values[n]->val;
- e->values[0]->count = e->values[n]->count;
- e->values[n]->val = t;
- e->values[n]->count = count;
- downheap_eqvals(e, --n, 0);
- }
- return;
-}
-
-
-/* *************************************************************** */
-/* ** Code based on heapsort from Sebastian Cyris **************** */
-/* ** ************************************************************ */
-void downheap_eqvals(equal_val_array * e, int n, int k)
-{
- int j;
- double v;
- int count;
-
- v = e->values[k]->val;
- count = e->values[k]->count;
- while (k <= n / 2) {
- j = k + k;
-
- if (j < n && e->values[j]->val < e->values[j + 1]->val)
- j++;
- if (v >= e->values[j]->val)
- break;
-
- e->values[k]->val = e->values[j]->val;
- e->values[k]->count = e->values[j]->count;
- k = j;
- }
-
- e->values[k]->val = v;
- e->values[k]->count = count;
-
- return;
-}
-
-
-/* *************************************************************** */
-/* ** Main function ********************************************** */
-/* *************************************************************** */
int main(int argc, char *argv[])
{
@@ -627,7 +107,7 @@
eqvals = NULL;
n = 0;
for (z = 0; z < depths; z++) {
- G_percent(z, depths - 1, 10);
+ G_percent(z, depths - 1, 2);
for (y = 0; y < rows; y++) {
for (x = 0; x < cols; x++) {
if (map_type == FCELL_TYPE) {
@@ -686,7 +166,7 @@
n = 0;
for (z = 0; z < depths; z++) {
- G_percent(z, depths - 1, 10);
+ G_percent(z, depths - 1, 2);
for (y = 0; y < rows; y++) {
for (x = 0; x < cols; x++) {
if (map_type == FCELL_TYPE) {
Copied: grass/branches/releasebranch_7_0/raster3d/r3.stats/support.c (from rev 60100, grass/trunk/raster3d/r3.stats/support.c)
===================================================================
--- grass/branches/releasebranch_7_0/raster3d/r3.stats/support.c (rev 0)
+++ grass/branches/releasebranch_7_0/raster3d/r3.stats/support.c 2014-07-23 13:22:59 UTC (rev 61373)
@@ -0,0 +1,485 @@
+/*
+ * r3.stats support functions
+ *
+ * Copyright (C) 2004-2014 by the GRASS Development Team
+ * Author(s): Soeren Gebbert
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ */
+
+#include <grass/gis.h>
+#include "local_proto.h"
+
+
+/* *************************************************************** */
+/* ***** allocate equal_val_array structure ********************* */
+/* *************************************************************** */
+equal_val_array *alloc_equal_val_array(int count)
+{
+ equal_val_array *p;
+ int i;
+
+
+ p = (equal_val_array *) G_calloc(1, sizeof(equal_val_array));
+ p->count = count;
+
+ p->values = (equal_val **) G_calloc(p->count, sizeof(equal_val *));
+
+ for (i = 0; i < p->count; i++)
+ p->values[i] = (equal_val *) G_calloc(1, sizeof(equal_val));
+
+
+ return p;
+}
+
+/* *************************************************************** */
+/* ***** add an equal_val to the equal_val_array structure ******* */
+/* *************************************************************** */
+equal_val_array *add_equal_val_to_array(equal_val_array *array, double val)
+{
+ equal_val_array *p = array;
+ int count;
+
+ if (p == NULL) {
+ p = alloc_equal_val_array(1);
+ p->values[0]->val = val;
+ p->values[0]->count = 1;
+ G_debug(5, "Create new equal_array with value %g\n", val);
+ }
+ else {
+ /*increase the count */
+ count = array->count;
+ count++;
+
+ /*new memory */
+ p->values =
+ (equal_val **) G_realloc(p->values, count * sizeof(equal_val *));
+ p->values[count - 1] = (equal_val *) G_calloc(1, sizeof(equal_val));
+ /*set the new value */
+ p->values[count - 1]->val = val;
+ p->values[count - 1]->count = 1;
+ /*set the new counter */
+ p->count = count;
+ G_debug(5, "Add new value %g at position %i\n", val, p->count);
+ }
+
+ return p;
+}
+
+/* *************************************************************** */
+/* ***** check if a value exists in the equal_val_array ********* */
+/* *************************************************************** */
+int check_equal_value(equal_val_array *array, double val)
+{
+ int i;
+
+ /*search if the new value exists and increase the counter */
+ if (array != NULL)
+ for (i = 0; i < array->count; i++) {
+ if (array->values[i]->val == val) {
+ array->values[i]->count++;
+ G_debug(5, "found value %g with count %i at pos %i\n", val,
+ array->values[i]->count, i);
+ return 1;
+ }
+ }
+
+ /*if it does not exists, add it to the array */
+ add_equal_val_to_array(array, val);
+
+ /*return 1 if found, 0 otherwise */
+ return 0;
+}
+
+
+/* *************************************************************** */
+/* ***** Release the memory of a equal_val_array structure ****** */
+/* *************************************************************** */
+void free_equal_val_array(equal_val_array *uvals)
+{
+ int i;
+
+ for (i = 0; i < uvals->count; i++) {
+ G_free(uvals->values[i]);
+ }
+
+ G_free(uvals->values);
+ G_free(uvals);
+ uvals = NULL;
+
+ return;
+}
+
+/* *************************************************************** */
+/* **** Create the structure which manages the statistical ******* */
+/* **** values for a value range or equal values **************** */
+/* *************************************************************** */
+stat_table *create_stat_table(int nsteps, equal_val_array *eqvals,
+ double min, double max)
+{
+ stat_table *p;
+ int i;
+ double step;
+
+ /* Memory */
+ p = (stat_table *) G_calloc(1, sizeof(stat_table));
+
+ p->null = (stat_row *) G_calloc(nsteps, sizeof(stat_row));
+ p->table = (stat_row **) G_calloc(nsteps, sizeof(stat_row *));
+
+ for (i = 0; i < nsteps; i++)
+ p->table[i] = (stat_row *) G_calloc(1, sizeof(stat_row));
+
+ /* Some value initializing */
+ p->null->min = 0;
+ p->null->max = 0;
+ p->null->vol = 0;
+ p->null->perc = 0;
+ p->null->count = 0;
+ p->null->num = nsteps + 1;
+
+ p->nsteps = nsteps;
+ p->sum_count = 0;
+ p->sum_vol = 0;
+ p->sum_perc = 0;
+ p->equal = 0;
+
+ /*if no equal values are provided, calculate the range and the length of each step */
+ if (!eqvals) {
+
+ /*calculate the steplength */
+ step = (max - min) / (nsteps);
+ p->equal = 0;
+
+ p->table[0]->min = min;
+ p->table[0]->max = min + step;
+ p->table[0]->num = 1;
+ p->table[0]->count = 0;
+ p->table[0]->vol = 0;
+ p->table[0]->perc = 0;
+
+ G_debug(3, "Step %i range min %.11lf max %.11lf\n", p->table[0]->num,
+ p->table[0]->min, p->table[0]->max);
+
+ for (i = 1; i < nsteps; i++) {
+ p->table[i]->min = p->table[i - 1]->max;
+ p->table[i]->max = p->table[i - 1]->max + step;
+ p->table[i]->num = i + 1;
+ p->table[i]->count = 0;
+ p->table[i]->vol = 0;
+ p->table[i]->perc = 0;
+ G_debug(5, "Step %i range min %.11lf max %.11lf\n",
+ p->table[i]->num, p->table[i]->min, p->table[i]->max);
+ }
+ /* the last value must be a bit larger */
+ p->table[nsteps - 1]->max += COMPARE_PRECISION;
+ }
+ else { /* Create equal value statistic */
+ p->equal = 1;
+ for (i = 0; i < eqvals->count; i++) {
+ p->table[i]->min = eqvals->values[i]->val; /* equal values have no range, set min and max to the same value */
+ p->table[i]->max = eqvals->values[i]->val;
+ p->table[i]->num = i + 1;
+ p->table[i]->count = eqvals->values[i]->count; /* the appearance count for each equal value */
+ p->table[i]->vol = 0;
+ p->table[i]->perc = 0;
+ G_debug(5, "Unique value %i = %g count %i\n", p->table[i]->num,
+ p->table[i]->min, p->table[i]->count);
+ }
+
+ }
+
+ return p;
+}
+
+/* *************************************************************** */
+/* ***** Release the memory of a stat_table structure ************ */
+/* *************************************************************** */
+void free_stat_table(stat_table *stats)
+{
+ int i;
+
+ for (i = 0; i < stats->nsteps; i++) {
+ G_free(stats->table[i]);
+ }
+
+ G_free(stats->table);
+ G_free(stats->null);
+
+ G_free(stats);
+ stats = NULL;
+
+ return;
+}
+
+/* *************************************************************** */
+/* **** Compute the volume, percentage and sums ****************** */
+/* *************************************************************** */
+void update_stat_table(stat_table *stats, RASTER3D_Region *region)
+{
+ int i;
+ double vol = region->ns_res * region->ew_res * region->tb_res;
+ int cellnum = region->rows * region->cols * region->depths;
+
+ /*Calculate volume and percentage for each range or equal value and the sum stats */
+ for (i = 0; i < stats->nsteps; i++) {
+ stats->table[i]->vol = stats->table[i]->count * vol;
+ stats->table[i]->perc =
+ (double)(100.0 * (double)stats->table[i]->count /
+ (double)cellnum);
+ stats->sum_count += stats->table[i]->count;
+ stats->sum_vol += stats->table[i]->vol;
+ stats->sum_perc += stats->table[i]->perc;
+ }
+
+ /*calculate the null value stats */
+ stats->null->vol = stats->null->count * vol;
+ stats->null->perc =
+ (double)(100.0 * (double)stats->null->count / (double)cellnum);
+
+ return;
+}
+
+/* *************************************************************** */
+/* *************************************************************** */
+/* *************************************************************** */
+void print_stat_table(stat_table *stats)
+{
+ int i;
+
+ if (stats->equal) {
+ /* 1234567 012345678901234567 0123456789012 0123456 0123456789 */
+ fprintf(stdout,
+ " num | value | volume | perc | cell count\n");
+
+ for (i = 0; i < stats->nsteps; i++) {
+ fprintf(stdout, "%7i %18.6lf %13.3lf %7.5lf %10i\n",
+ stats->table[i]->num, stats->table[i]->min,
+ stats->table[i]->vol, stats->table[i]->perc,
+ stats->table[i]->count);
+ }
+ fprintf(stdout,
+ "%7i * %13.3lf %7.5lf %10i\n",
+ stats->null->num, stats->null->vol, stats->null->perc,
+ stats->null->count);
+
+ fprintf(stdout, "\nNumber of groups with equal values: %i",
+ stats->nsteps);
+ }
+ else {
+ /* 1234567 012345678901234567 012345678901234567 0123456789012 0123456 0123456789 */
+ fprintf(stdout,
+ " num | minimum <= value | value < maximum | volume | perc | cell count\n");
+
+ for (i = 0; i < stats->nsteps; i++) {
+ fprintf(stdout,
+ "%7i %18.9lf %18.9lf %13.3lf %7.5lf %10i\n",
+ stats->table[i]->num, stats->table[i]->min,
+ stats->table[i]->max, stats->table[i]->vol,
+ stats->table[i]->perc, stats->table[i]->count);
+
+ }
+ fprintf(stdout,
+ "%7i * * %13.3lf %7.5lf %10i\n",
+ stats->null->num, stats->null->vol, stats->null->perc,
+ stats->null->count);
+ }
+
+ fprintf(stdout,
+ "\nSum of non Null cells: \n\tVolume = %13.3lf \n\tPercentage = %7.3lf \n\tCell count = %i\n",
+ stats->sum_vol, stats->sum_perc, stats->sum_count);
+
+ fprintf(stdout,
+ "\nSum of all cells: \n\tVolume = %13.3lf \n\tPercentage = %7.3lf \n\tCell count = %i\n",
+ stats->sum_vol + stats->null->vol,
+ stats->sum_perc + stats->null->perc,
+ stats->sum_count + stats->null->count);
+
+
+ return;
+}
+
+/* *************************************************************** */
+/* Make an entry in the statistic table based on range value check */
+/* *************************************************************** */
+void check_range_value(stat_table *stats, double value)
+{
+
+ /* this is very expensive for large range arrays! */
+ /*
+ * int i;
+ * for (i = 0; i < stats->nsteps; i++) {
+ * if (value >= stats->table[i]->min && value < stats->table[i]->max) {
+ * stats->table[i]->count++;
+ * }
+ * }
+ */
+
+ /* the much faster tree search is used */
+ tree_search_range(stats, 0, stats->nsteps - 1, value);
+
+ return;
+}
+
+/* *************************************************************** */
+/* tree search method developed by Soeren Gebbert *************** */
+/* *************************************************************** */
+/*
+ * This is a divide and conquer tree search algorithm
+ *
+ * The algorithm counts how many values are located in each range.
+ * It divides therefor the range array in smaller pices.
+ *
+ * e.g:
+ * 0 1 2 3 4 5 6 7 8
+ * [[0,1],[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9]]
+ * / \
+ * [[0,1],[1,2],[2,3],[3,4]] [[4,5],[5,6],[6,7],[7,8],[8,9]]
+ * / \ / \
+ * [[0,1],[1,2]] [[2,3],[3,4]] [[4,5],[5,6]] [[6,7],[7,8],[8,9]]
+ * / \ / \ / \ / \
+ * [[0,1]] [[1,2]] [[2,3]] [[3,4]] [[4,5]] [[5,6]] [[6,7]] [[7,8],[8,9]]
+ * / \
+ * [[7,8]] [[8,9]]
+ *
+ * Searching the value 5.5 will result in the follwoing steps:
+ *
+ * value 5.5 in range of
+ * 1.) left array index = 0 - right array index = 8 -> range [0 - 9]
+ * 2.) left array index = 4 - right array index = 8 -> range [4 - 9]
+ * 3.) left array index = 4 - right array index = 5 -> range [4 - 6]
+ * 4.) the value is located in range array 5 with a range of [5 - 6]
+ * 5.) now increase the value range counter of range array with index 5
+ *
+ * */
+void tree_search_range(stat_table *stats, int left, int right, double value)
+{
+ int size = right - left;
+
+ G_debug(5,
+ "Search value %g in array size %i left border index %i right border index %i\n",
+ value, size, left, right);
+
+ /*if left and right are equal */
+ if (size == 0) {
+ stats->table[left]->count++;
+ return;
+ }
+ else if (size == 1) { /* if the size is one, check directly */
+
+ if (value >= stats->table[left]->min &&
+ value < stats->table[left]->max) {
+ stats->table[left]->count++;
+ }
+ else if (value >= stats->table[right]->min &&
+ value < stats->table[right]->max) {
+ stats->table[right]->count++;
+ }
+ return;
+ }
+ else {
+ if ((size % 2 == 0)) { /*even */
+ /*left side */
+ right = left + (size) / 2;
+
+ if (value >= stats->table[left]->min &&
+ value < stats->table[right]->max) {
+ tree_search_range(stats, left, right, value);
+ return;
+ }
+
+ /*right side */
+ left += (size) / 2;
+ right = left + (size) / 2;
+
+ if (value >= stats->table[left]->min &&
+ value < stats->table[right]->max) {
+ tree_search_range(stats, left, right, value);
+ return;
+ }
+ }
+ else { /*odd */
+ /*left side */
+ right = left + (size - 1) / 2;
+
+ if (value >= stats->table[left]->min &&
+ value < stats->table[right]->max) {
+ tree_search_range(stats, left, right, value);
+ return;
+ }
+
+ /*right side */
+ left += (size - 1) / 2;
+ right = left + (size - 1) / 2 + 1; /*the right array is one col larger */
+
+ if (value >= stats->table[left]->min &&
+ value < stats->table[right]->max) {
+ tree_search_range(stats, left, right, value);
+ return;
+ }
+ }
+ }
+
+ return;
+}
+
+/* *************************************************************** */
+/* ****** heapsort for equal value arrays of size n ************** */
+/* ** Code based on heapsort from Sebastian Cyris **************** */
+/* *************************************************************** */
+void heapsort_eqvals(equal_val_array *e, int n)
+{
+ int k;
+ double t;
+ int count = 0;
+
+ --n;
+
+ for (k = n / 2; k >= 0; k--)
+ downheap_eqvals(e, n, k);
+
+ while (n > 0) {
+ t = e->values[0]->val;
+ count = e->values[0]->count;
+ e->values[0]->val = e->values[n]->val;
+ e->values[0]->count = e->values[n]->count;
+ e->values[n]->val = t;
+ e->values[n]->count = count;
+ downheap_eqvals(e, --n, 0);
+ }
+ return;
+}
+
+
+/* *************************************************************** */
+/* ** Code based on heapsort from Sebastian Cyris **************** */
+/* ** ************************************************************ */
+void downheap_eqvals(equal_val_array *e, int n, int k)
+{
+ int j;
+ double v;
+ int count;
+
+ v = e->values[k]->val;
+ count = e->values[k]->count;
+ while (k <= n / 2) {
+ j = k + k;
+
+ if (j < n && e->values[j]->val < e->values[j + 1]->val)
+ j++;
+ if (v >= e->values[j]->val)
+ break;
+
+ e->values[k]->val = e->values[j]->val;
+ e->values[k]->count = e->values[j]->count;
+ k = j;
+ }
+
+ e->values[k]->val = v;
+ e->values[k]->count = count;
+
+ return;
+}
More information about the grass-commit
mailing list