From fa1e07d6637d844972902059785adf06a725bfbf Mon Sep 17 00:00:00 2001 From: Markus Neteler Date: Thu, 19 Nov 2020 15:07:07 +0100 Subject: [PATCH] r.geomorphon: Run grass_indent.sh on the source (#1102) * r.geomorphon: Run grass_indent.sh on the source --- raster/r.geomorphon/geom.c | 144 ++--- raster/r.geomorphon/local_proto.h | 48 +- raster/r.geomorphon/main.c | 980 +++++++++++++++--------------- raster/r.geomorphon/memory.c | 95 +-- raster/r.geomorphon/multires.c | 12 +- raster/r.geomorphon/pattern.c | 208 +++---- 6 files changed, 744 insertions(+), 743 deletions(-) diff --git a/raster/r.geomorphon/geom.c b/raster/r.geomorphon/geom.c index 8d3fe3a63e6..6c7401e0fed 100644 --- a/raster/r.geomorphon/geom.c +++ b/raster/r.geomorphon/geom.c @@ -1,7 +1,7 @@ #include "local_proto.h" -/* static double dirs[8] = { 0.7854, 0., 5.4978, 4.7124, 3.9270, 3.1416, 2.3562, 1.5708 };*/ /* radians */ -static double sins[8] = { 0.7071067812, 0, -0.7071067812, -1, -0.7071067812, 0, 0.7071067812, 1 }; /* sinus */ -static double coss[8] = { 0.7071067812, 1, 0.7071067812, 0, -0.7071067812, -1, -0.7071067812, 0 }; /* cosinus */ + /* static double dirs[8] = { 0.7854, 0., 5.4978, 4.7124, 3.9270, 3.1416, 2.3562, 1.5708 }; *//* radians */ +static double sins[8] = { 0.7071067812, 0, -0.7071067812, -1, -0.7071067812, 0, 0.7071067812, 1 }; /* sinus */ +static double coss[8] = { 0.7071067812, 1, 0.7071067812, 0, -0.7071067812, -1, -0.7071067812, 0 }; /* cosinus */ /* DIRS in DEGREES from NORTH: 45,0,315,270,225,180,135,90 */ @@ -21,25 +21,25 @@ unsigned int ternary_rotate(unsigned int value) int i, j, k; for (i = 0; i < 8; i++) { - pattern[i] = value % 3; - rev_pattern[7 - i] = value % 3; - value /= 3; + pattern[i] = value % 3; + rev_pattern[7 - i] = value % 3; + value /= 3; } for (j = 0; j < 8; j++) { - power = 1; - tmp_code = 0; - tmp_rev_code = 0; - for (i = 0; i < 8; i++) { - k = (i - j) < 0 ? j - 8 : j; - tmp_pattern[i] = pattern[i - k]; - tmp_rev_pattern[i] = rev_pattern[i - k]; - tmp_code += tmp_pattern[i] * power; - tmp_rev_code += tmp_rev_pattern[i] * power; - power *= 3; - } - code = tmp_code < code ? tmp_code : code; - rev_code = tmp_rev_code < rev_code ? tmp_rev_code : rev_code; + power = 1; + tmp_code = 0; + tmp_rev_code = 0; + for (i = 0; i < 8; i++) { + k = (i - j) < 0 ? j - 8 : j; + tmp_pattern[i] = pattern[i - k]; + tmp_rev_pattern[i] = rev_pattern[i - k]; + tmp_code += tmp_pattern[i] * power; + tmp_rev_code += tmp_rev_pattern[i] * power; + power *= 3; + } + code = tmp_code < code ? tmp_code : code; + rev_code = tmp_rev_code < rev_code ? tmp_rev_code : rev_code; } return code < rev_code ? code : rev_code; } @@ -50,17 +50,17 @@ int determine_form(int num_minus, int num_plus) * simple approach to determine form pattern */ const FORMS forms[9][9] = { - /* minus ------------- plus ---------------- */ - /* 0 1 2 3 4 5 6 7 8 */ - /* 0 */ {FL, FL, FL, FS, FS, VL, VL, VL, PT}, - /* 1 */ {FL, FL, FS, FS, FS, VL, VL, VL, __}, - /* 2 */ {FL, SH, SL, SL, CN, CN, VL, __, __}, - /* 3 */ {SH, SH, SL, SL, SL, CN, __, __, __}, - /* 4 */ {SH, SH, CV, SL, SL, __, __, __, __}, - /* 5 */ {RI, RI, CV, CV, __, __, __, __, __}, - /* 6 */ {RI, RI, RI, __, __, __, __, __, __}, - /* 7 */ {RI, RI, __, __, __, __, __, __, __}, - /* 8 */ {PK, __, __, __, __, __, __, __, __}, + /* minus ------------- plus ---------------- */ + /* 0 1 2 3 4 5 6 7 8 */ + /* 0 */ {FL, FL, FL, FS, FS, VL, VL, VL, PT}, + /* 1 */ {FL, FL, FS, FS, FS, VL, VL, VL, __}, + /* 2 */ {FL, SH, SL, SL, CN, CN, VL, __, __}, + /* 3 */ {SH, SH, SL, SL, SL, CN, __, __, __}, + /* 4 */ {SH, SH, CV, SL, SL, __, __, __, __}, + /* 5 */ {RI, RI, CV, CV, __, __, __, __, __}, + /* 6 */ {RI, RI, RI, __, __, __, __, __, __}, + /* 7 */ {RI, RI, __, __, __, __, __, __, __}, + /* 8 */ {PK, __, __, __, __, __, __, __, __}, }; /* legend: @@ -86,14 +86,14 @@ int determine_binary(int *pattern, int sign) unsigned char binary = 0, result = 255, test = 0; for (i = 0, n = 1; i < 8; i++, n *= 2) - binary += (pattern[i] == sign) ? n : 0; + binary += (pattern[i] == sign) ? n : 0; /* rotate */ for (i = 0; i < 8; ++i) { - if ((i &= 7) == 0) - test = binary; - else - test = (binary << i) | (binary >> (8 - i)); - result = (result < test) ? result : test; + if ((i &= 7) == 0) + test = binary; + else + test = (binary << i) | (binary >> (8 - i)); + result = (result < test) ? result : test; } return (int)result; } @@ -106,11 +106,11 @@ int rotate(unsigned char binary) /* rotate */ for (i = 0; i < 8; ++i) { - if ((i &= 7) == 0) - test = binary; - else - test = (binary << i) | (binary >> (8 - i)); - result = (result < test) ? result : test; + if ((i &= 7) == 0) + test = binary; + else + test = (binary << i) | (binary >> (8 - i)); + result = (result < test) ? result : test; } return (int)result; } @@ -122,7 +122,7 @@ int determine_ternary(int *pattern) int power, i; for (i = 0, power = 1; i < 8; ++i, power *= 3) - ternary_code += (pattern[i] + 1) * power; + ternary_code += (pattern[i] + 1) * power; return global_ternary_codes[ternary_code]; } @@ -133,7 +133,7 @@ float intensity(float *elevation, int pattern_size) int i; for (i = 0; i < 8; i++) - sum_elevation -= elevation[i]; + sum_elevation -= elevation[i]; return sum_elevation / (float)pattern_size; } @@ -145,19 +145,19 @@ float exposition(float *elevation) int i; for (i = 0; i < 8; i++) - max = fabs(elevation[i]) > fabs(max) ? elevation[i] : max; + max = fabs(elevation[i]) > fabs(max) ? elevation[i] : max; return -max; } float range(float *elevation) { /* calculate relative difference in visible range of central cell */ - float max = -90000000000., min = 9000000000000.; /* should be enough */ + float max = -90000000000., min = 9000000000000.; /* should be enough */ int i; for (i = 0; i < 8; i++) { - max = elevation[i] > max ? elevation[i] : max; - min = elevation[i] < min ? elevation[i] : min; + max = elevation[i] > max ? elevation[i] : max; + min = elevation[i] < min ? elevation[i] : min; } return max - min; @@ -171,11 +171,11 @@ float variance(float *elevation, int pattern_size) int i; for (i = 0; i < 8; i++) - sum_elevation += elevation[i]; + sum_elevation += elevation[i]; sum_elevation /= (float)pattern_size; for (i = 0; i < 8; i++) - variance += - ((sum_elevation - elevation[i]) * (sum_elevation - elevation[i])); + variance += + ((sum_elevation - elevation[i]) * (sum_elevation - elevation[i])); return variance / (float)pattern_size; } @@ -188,14 +188,14 @@ int radial2cartesian(PATTERN * pattern) int i; for (i = 0; i < 8; ++i) - if (pattern->distance > 0) { - pattern->x[i] = pattern->distance[i] * sins[i]; - pattern->y[i] = pattern->distance[i] * coss[i]; - } - else { - pattern->x[i] = 0; - pattern->y[i] = 0; - } + if (pattern->distance > 0) { + pattern->x[i] = pattern->distance[i] * sins[i]; + pattern->y[i] = pattern->distance[i] * coss[i]; + } + else { + pattern->x[i] = 0; + pattern->y[i] = 0; + } return 0; } @@ -205,15 +205,15 @@ float extends(PATTERN * pattern, int pattern_size) float area = 0; for (i = 0, j = 1; i < 8; ++i, ++j) { - j = j < 8 ? j : 0; - area += - (pattern->x[i] * pattern->y[j] - pattern->x[j] * pattern->y[i]); + j = j < 8 ? j : 0; + area += + (pattern->x[i] * pattern->y[j] - pattern->x[j] * pattern->y[i]); } return fabs(area) / 2.; } int shape(PATTERN * pattern, int pattern_size, float *azimuth, - float *elongation, float *width) + float *elongation, float *width) { /* calculates azimuth, elongation and width of geomorphon's polygon */ int i; @@ -226,10 +226,10 @@ int shape(PATTERN * pattern, int pattern_size, float *azimuth, double rxmin, rxmax, rymin, rymax; for (i = 0; i < 8; ++i) { - avg_y += pattern->y[i]; - avg_x += pattern->x[i]; - avg_x_square += pattern->x[i] * pattern->x[i]; - avg_x_y += pattern->x[i] * pattern->y[i]; + avg_y += pattern->y[i]; + avg_x += pattern->x[i]; + avg_x_square += pattern->x[i] * pattern->x[i]; + avg_x_y += pattern->x[i] * pattern->y[i]; } avg_y /= (float)pattern_size; avg_x /= (float)pattern_size; @@ -245,12 +245,12 @@ int shape(PATTERN * pattern, int pattern_size, float *azimuth, rxmin = rxmax = pattern->x[0] * cosine - pattern->y[0] * sine; rymin = rymax = pattern->x[0] * sine + pattern->y[0] * cosine; for (i = 1; i < 8; ++i) { - rx = pattern->x[i] * cosine - pattern->y[i] * sine; - ry = pattern->x[i] * sine + pattern->y[i] * cosine; - rxmin = rx < rxmin ? rx : rxmin; - rxmax = rx > rxmax ? rx : rxmax; - rymin = ry < rymin ? ry : rymin; - rymax = ry > rymax ? ry : rymax; + rx = pattern->x[i] * cosine - pattern->y[i] * sine; + ry = pattern->x[i] * sine + pattern->y[i] * cosine; + rxmin = rx < rxmin ? rx : rxmin; + rxmax = rx > rxmax ? rx : rxmax; + rymin = ry < rymin ? ry : rymin; + rymax = ry > rymax ? ry : rymax; } rx = (rxmax - rxmin); ry = (rymax - rymin); diff --git a/raster/r.geomorphon/local_proto.h b/raster/r.geomorphon/local_proto.h index e88a93259f1..75150be3031 100644 --- a/raster/r.geomorphon/local_proto.h +++ b/raster/r.geomorphon/local_proto.h @@ -15,11 +15,11 @@ #endif -#ifndef PI2 /* PI/2 */ +#ifndef PI2 /* PI/2 */ #define PI2 (2*atan(1)) #endif -#ifndef PI4 /* PI/4 */ +#ifndef PI4 /* PI/4 */ #define PI4 (atan(1)) #endif @@ -27,7 +27,7 @@ #define PI (4*atan(1)) #endif -#ifndef M2PI /* 2*PI */ +#ifndef M2PI /* 2*PI */ #define M2PI (8*atan(1)) #endif @@ -60,11 +60,11 @@ typedef struct char elevname[150]; RASTER_MAP_TYPE raster_type; FCELL **elev; - int fd; /* file descriptor */ + int fd; /* file descriptor */ } MAPS; typedef struct -{ /* struct is used both for interface and output */ +{ /* struct is used both for interface and output */ char *name; int required; char *description; @@ -90,24 +90,24 @@ typedef struct int pattern[8]; float elevation[8]; double distance[8]; - double x[8], y[8]; /* cartesian coordinates of geomorphon */ + double x[8], y[8]; /* cartesian coordinates of geomorphon */ } PATTERN; typedef enum { - ZERO, /* zero cats do not accept zero category */ - FL, /* flat */ - PK, /* peak, summit */ - RI, /* ridge */ - SH, /* shoulder */ - CV, /* convex slope */ - SL, /* slope */ - CN, /* concave slope */ - FS, /* footslope */ - VL, /* valley */ - PT, /* pit, depression */ - __, /* error */ - CNT /* counter */ + ZERO, /* zero cats do not accept zero category */ + FL, /* flat */ + PK, /* peak, summit */ + RI, /* ridge */ + SH, /* shoulder */ + CV, /* convex slope */ + SL, /* slope */ + CN, /* concave slope */ + FS, /* footslope */ + VL, /* valley */ + PT, /* pit, depression */ + __, /* error */ + CNT /* counter */ } FORMS; typedef struct @@ -134,7 +134,7 @@ GLOBAL int nrows, ncols, row_radius_size, row_buffer_size; GLOBAL int search_cells, skip_cells, step_cells, start_cells; GLOBAL int num_of_steps; GLOBAL double search_distance, skip_distance, flat_distance; -GLOBAL double start_distance, step_distance; /* multiresolution mode */ +GLOBAL double start_distance, step_distance; /* multiresolution mode */ GLOBAL double flat_threshold, flat_threshold_height; GLOBAL double H, V; GLOBAL struct Cell_head window; @@ -164,7 +164,7 @@ float exposition(float *elevation); float range(float *elevation); float variance(float *elevation, int n); int shape(PATTERN * pattern, int pattern_size, float *azimuth, - float *elongation, float *width); + float *elongation, float *width); float extends(PATTERN * pattern, int pattern_size); int radial2cartesian(PATTERN *); @@ -172,8 +172,8 @@ int radial2cartesian(PATTERN *); int reset_multi_patterns(void); int calc_multi_patterns(int row, int cur_row, int col); int update_pattern(int k, int i, - double zenith_height, double zenith_distance, - double zenith_angle, double nadir_height, - double nadir_distance, double nadir_angle); + double zenith_height, double zenith_distance, + double zenith_angle, double nadir_height, + double nadir_distance, double nadir_angle); #endif diff --git a/raster/r.geomorphon/main.c b/raster/r.geomorphon/main.c index 7219493ebab..fd6217b009e 100644 --- a/raster/r.geomorphon/main.c +++ b/raster/r.geomorphon/main.c @@ -23,73 +23,73 @@ #include "local_proto.h" typedef enum { i_dem, o_forms, o_ternary, o_positive, o_negative, o_intensity, - o_exposition, + o_exposition, o_range, o_variance, o_elongation, o_azimuth, o_extend, o_width, io_size } outputs; int main(int argc, char **argv) { - IO rasters[] = { /* rasters stores output buffers */ - {"elevation", YES, "Name of input elevation raster map", "input", UNKNOWN, -1, NULL}, /* WARNING: this one map is input */ - {"forms", NO, "Most common geomorphic forms", "Patterns", CELL_TYPE, - -1, NULL}, - {"ternary", NO, "Code of ternary patterns", "Patterns", CELL_TYPE, -1, - NULL}, - {"positive", NO, "Code of binary positive patterns", "Patterns", - CELL_TYPE, -1, NULL}, - {"negative", NO, "Code of binary negative patterns", "Patterns", - CELL_TYPE, -1, NULL}, - {"intensity", NO, - "Rasters containing mean relative elevation of the form", "Geometry", - FCELL_TYPE, -1, NULL}, - {"exposition", NO, - "Rasters containing maximum difference between extend and central cell", - "Geometry", FCELL_TYPE, -1, NULL}, - {"range", NO, - "Rasters containing difference between max and min elevation of the form extend", - "Geometry", FCELL_TYPE, -1, NULL}, - {"variance", NO, "Rasters containing variance of form boundary", - "Geometry", FCELL_TYPE, -1, NULL}, - {"elongation", NO, "Rasters containing local elongation", "Geometry", - FCELL_TYPE, -1, NULL}, - {"azimuth", NO, "Rasters containing local azimuth of the elongation", - "Geometry", FCELL_TYPE, -1, NULL}, - {"extend", NO, "Rasters containing local extend (area) of the form", - "Geometry", FCELL_TYPE, -1, NULL}, - {"width", NO, "Rasters containing local width of the form", - "Geometry", FCELL_TYPE, -1, NULL} - }; /* adding more maps change IOSIZE macro */ - - CATCOLORS ccolors[CNT] = { /* colors and cats for forms */ - {ZERO, 0, 0, 0, "forms"}, - {FL, 220, 220, 220, "flat"}, - {PK, 56, 0, 0, "summit"}, - {RI, 200, 0, 0, "ridge"}, - {SH, 255, 80, 20, "shoulder"}, - {CV, 250, 210, 60, "spur"}, - {SL, 255, 255, 60, "slope"}, - {CN, 180, 230, 20, "hollow"}, - {FS, 60, 250, 150, "footslope"}, - {VL, 0, 0, 255, "valley"}, - {PT, 0, 0, 56, "depression"}, - {__, 255, 0, 255, "ERROR"} + IO rasters[] = { /* rasters stores output buffers */ + {"elevation", YES, "Name of input elevation raster map", "input", UNKNOWN, -1, NULL}, /* WARNING: this one map is input */ + {"forms", NO, "Most common geomorphic forms", "Patterns", CELL_TYPE, + -1, NULL}, + {"ternary", NO, "Code of ternary patterns", "Patterns", CELL_TYPE, -1, + NULL}, + {"positive", NO, "Code of binary positive patterns", "Patterns", + CELL_TYPE, -1, NULL}, + {"negative", NO, "Code of binary negative patterns", "Patterns", + CELL_TYPE, -1, NULL}, + {"intensity", NO, + "Rasters containing mean relative elevation of the form", "Geometry", + FCELL_TYPE, -1, NULL}, + {"exposition", NO, + "Rasters containing maximum difference between extend and central cell", + "Geometry", FCELL_TYPE, -1, NULL}, + {"range", NO, + "Rasters containing difference between max and min elevation of the form extend", + "Geometry", FCELL_TYPE, -1, NULL}, + {"variance", NO, "Rasters containing variance of form boundary", + "Geometry", FCELL_TYPE, -1, NULL}, + {"elongation", NO, "Rasters containing local elongation", "Geometry", + FCELL_TYPE, -1, NULL}, + {"azimuth", NO, "Rasters containing local azimuth of the elongation", + "Geometry", FCELL_TYPE, -1, NULL}, + {"extend", NO, "Rasters containing local extend (area) of the form", + "Geometry", FCELL_TYPE, -1, NULL}, + {"width", NO, "Rasters containing local width of the form", + "Geometry", FCELL_TYPE, -1, NULL} + }; /* adding more maps change IOSIZE macro */ + + CATCOLORS ccolors[CNT] = { /* colors and cats for forms */ + {ZERO, 0, 0, 0, "forms"}, + {FL, 220, 220, 220, "flat"}, + {PK, 56, 0, 0, "summit"}, + {RI, 200, 0, 0, "ridge"}, + {SH, 255, 80, 20, "shoulder"}, + {CV, 250, 210, 60, "spur"}, + {SL, 255, 255, 60, "slope"}, + {CN, 180, 230, 20, "hollow"}, + {FS, 60, 250, 150, "footslope"}, + {VL, 0, 0, 255, "valley"}, + {PT, 0, 0, 56, "depression"}, + {__, 255, 0, 255, "ERROR"} }; struct GModule *module; struct Option - *opt_input, - *opt_output[io_size], - *par_search_radius, - *par_skip_radius, - *par_flat_threshold, - *par_flat_distance, - *par_multi_prefix, *par_multi_step, *par_multi_start; + *opt_input, + *opt_output[io_size], + *par_search_radius, + *par_skip_radius, + *par_flat_threshold, + *par_flat_distance, + *par_multi_prefix, *par_multi_step, *par_multi_start; struct Flag *flag_units, *flag_extended; struct History history; int i; - int meters = 0, multires = 0, extended = 0; /* flags */ + int meters = 0, multires = 0, extended = 0; /* flags */ int row, cur_row, col; int pattern_size; double max_resolution; @@ -97,458 +97,458 @@ int main(int argc, char **argv) G_gisinit(argv[0]); - { /* interface parameters */ - module = G_define_module(); - module->description = - _("Calculates geomorphons (terrain forms) and associated geometry using machine vision approach."); + { /* interface parameters */ + module = G_define_module(); + module->description = + _("Calculates geomorphons (terrain forms) and associated geometry using machine vision approach."); G_add_keyword(_("raster")); - G_add_keyword(_("geomorphons")); - G_add_keyword(_("terrain patterns")); - G_add_keyword(_("machine vision geomorphometry")); - - opt_input = G_define_standard_option(G_OPT_R_ELEV); - opt_input->key = rasters[0].name; - opt_input->required = rasters[0].required; - opt_input->description = _(rasters[0].description); - - for (i = 1; i < io_size; ++i) { /* WARNING: loop starts from one, zero is for input */ - opt_output[i] = G_define_standard_option(G_OPT_R_OUTPUT); - opt_output[i]->key = rasters[i].name; - opt_output[i]->required = NO; - opt_output[i]->description = _(rasters[i].description); - opt_output[i]->guisection = _(rasters[i].gui); - } - - par_search_radius = G_define_option(); - par_search_radius->key = "search"; - par_search_radius->type = TYPE_INTEGER; - par_search_radius->answer = "3"; - par_search_radius->required = YES; - par_search_radius->description = _("Outer search radius"); - - par_skip_radius = G_define_option(); - par_skip_radius->key = "skip"; - par_skip_radius->type = TYPE_INTEGER; - par_skip_radius->answer = "0"; - par_skip_radius->required = YES; - par_skip_radius->description = _("Inner search radius"); - - par_flat_threshold = G_define_option(); - par_flat_threshold->key = "flat"; - par_flat_threshold->type = TYPE_DOUBLE; - par_flat_threshold->answer = "1"; - par_flat_threshold->required = YES; - par_flat_threshold->description = _("Flatness threshold (degrees)"); - - par_flat_distance = G_define_option(); - par_flat_distance->key = "dist"; - par_flat_distance->type = TYPE_DOUBLE; - par_flat_distance->answer = "0"; - par_flat_distance->required = YES; - par_flat_distance->description = - _("Flatness distance, zero for none"); - - par_multi_prefix = G_define_option(); - par_multi_prefix->key = "prefix"; - par_multi_prefix->type = TYPE_STRING; - par_multi_prefix->description = - _("Prefix for maps resulting from multiresolution approach"); - par_multi_prefix->guisection = _("Multires"); - - par_multi_step = G_define_option(); - par_multi_step->key = "step"; - par_multi_step->type = TYPE_DOUBLE; - par_multi_step->answer = "0"; - par_multi_step->description = - _("Distance step for every iteration (zero to omit)"); - par_multi_step->guisection = _("Multires"); - - par_multi_start = G_define_option(); - par_multi_start->key = "start"; - par_multi_start->type = TYPE_DOUBLE; - par_multi_start->answer = "0"; - par_multi_start->description = - _("Distance where search will start in multiple mode (zero to omit)"); - par_multi_start->guisection = _("Multires"); - - flag_units = G_define_flag(); - flag_units->key = 'm'; - flag_units->description = - _("Use meters to define search units (default is cells)"); - - flag_extended = G_define_flag(); - flag_extended->key = 'e'; - flag_extended->description = _("Use extended form correction"); - - if (G_parser(argc, argv)) - exit(EXIT_FAILURE); + G_add_keyword(_("geomorphons")); + G_add_keyword(_("terrain patterns")); + G_add_keyword(_("machine vision geomorphometry")); + + opt_input = G_define_standard_option(G_OPT_R_ELEV); + opt_input->key = rasters[0].name; + opt_input->required = rasters[0].required; + opt_input->description = _(rasters[0].description); + + for (i = 1; i < io_size; ++i) { /* WARNING: loop starts from one, zero is for input */ + opt_output[i] = G_define_standard_option(G_OPT_R_OUTPUT); + opt_output[i]->key = rasters[i].name; + opt_output[i]->required = NO; + opt_output[i]->description = _(rasters[i].description); + opt_output[i]->guisection = _(rasters[i].gui); + } + + par_search_radius = G_define_option(); + par_search_radius->key = "search"; + par_search_radius->type = TYPE_INTEGER; + par_search_radius->answer = "3"; + par_search_radius->required = YES; + par_search_radius->description = _("Outer search radius"); + + par_skip_radius = G_define_option(); + par_skip_radius->key = "skip"; + par_skip_radius->type = TYPE_INTEGER; + par_skip_radius->answer = "0"; + par_skip_radius->required = YES; + par_skip_radius->description = _("Inner search radius"); + + par_flat_threshold = G_define_option(); + par_flat_threshold->key = "flat"; + par_flat_threshold->type = TYPE_DOUBLE; + par_flat_threshold->answer = "1"; + par_flat_threshold->required = YES; + par_flat_threshold->description = _("Flatness threshold (degrees)"); + + par_flat_distance = G_define_option(); + par_flat_distance->key = "dist"; + par_flat_distance->type = TYPE_DOUBLE; + par_flat_distance->answer = "0"; + par_flat_distance->required = YES; + par_flat_distance->description = + _("Flatness distance, zero for none"); + + par_multi_prefix = G_define_option(); + par_multi_prefix->key = "prefix"; + par_multi_prefix->type = TYPE_STRING; + par_multi_prefix->description = + _("Prefix for maps resulting from multiresolution approach"); + par_multi_prefix->guisection = _("Multires"); + + par_multi_step = G_define_option(); + par_multi_step->key = "step"; + par_multi_step->type = TYPE_DOUBLE; + par_multi_step->answer = "0"; + par_multi_step->description = + _("Distance step for every iteration (zero to omit)"); + par_multi_step->guisection = _("Multires"); + + par_multi_start = G_define_option(); + par_multi_start->key = "start"; + par_multi_start->type = TYPE_DOUBLE; + par_multi_start->answer = "0"; + par_multi_start->description = + _("Distance where search will start in multiple mode (zero to omit)"); + par_multi_start->guisection = _("Multires"); + + flag_units = G_define_flag(); + flag_units->key = 'm'; + flag_units->description = + _("Use meters to define search units (default is cells)"); + + flag_extended = G_define_flag(); + flag_extended->key = 'e'; + flag_extended->description = _("Use extended form correction"); + + if (G_parser(argc, argv)) + exit(EXIT_FAILURE); } - { /* calculate parameters */ - int num_outputs = 0; - double search_radius, skip_radius, start_radius, step_radius; - double ns_resolution; - - multires = (par_multi_prefix->answer) ? 1 : 0; - for (i = 1; i < io_size; ++i) /* check for outputs */ - if (opt_output[i]->answer) { - if (G_legal_filename(opt_output[i]->answer) < 0) - G_fatal_error(_("<%s> is an illegal file name"), - opt_output[i]->answer); - num_outputs++; - } - if (!num_outputs && !multires) - G_fatal_error(_("At least one output is required, e.g. %s"), - opt_output[o_forms]->key); - - meters = (flag_units->answer != 0); - extended = (flag_extended->answer != 0); - nrows = Rast_window_rows(); - ncols = Rast_window_cols(); - Rast_get_window(&window); - G_begin_distance_calculations(); - - if (G_projection() == PROJECTION_LL) { /* for LL max_res should be NS */ - ns_resolution = - G_distance(0, Rast_row_to_northing(0, &window), 0, - Rast_row_to_northing(1, &window)); - max_resolution = ns_resolution; - } - else { - max_resolution = MAX(window.ns_res, window.ew_res); /* max_resolution MORE meters per cell */ - ns_resolution = window.ns_res; - } - - /* search distance */ - search_radius = atof(par_search_radius->answer); - search_cells = - meters ? (int)(search_radius / max_resolution) : search_radius; - if (search_cells < 1) - G_fatal_error(_("Search radius size must cover at least 1 cell")); - row_radius_size = - meters ? ceil(search_radius / ns_resolution) : search_radius; - row_buffer_size = row_radius_size * 2 + 1; - search_distance = - (meters) ? search_radius : ns_resolution * search_cells; - - /* skip distance */ - skip_radius = atof(par_skip_radius->answer); - skip_cells = - meters ? (int)(skip_radius / max_resolution) : skip_radius; - if (skip_cells >= search_cells) - G_fatal_error(_("Skip radius size must be at least 1 cell lower than radius")); - skip_distance = (meters) ? skip_radius : ns_resolution * skip_cells; - - /* flatness parameters */ - flat_threshold = atof(par_flat_threshold->answer); - if (flat_threshold <= 0.) - G_fatal_error(_("Flatness threshold must be grater than 0")); - flat_threshold = DEGREE2RAD(flat_threshold); - - flat_distance = atof(par_flat_distance->answer); - flat_distance = - (meters) ? flat_distance : ns_resolution * flat_distance; - flat_threshold_height = tan(flat_threshold) * flat_distance; - if ((flat_distance > 0 && flat_distance <= skip_distance) || - flat_distance >= search_distance) { - G_warning(_("Flatness distance should be between skip and search radius. Otherwise ignored")); - flat_distance = 0; - } - if (multires) { - start_radius = atof(par_multi_start->answer); - start_cells = - meters ? (int)(start_radius / max_resolution) : start_radius; - if (start_cells <= skip_cells) - start_cells = skip_cells + 1; - start_distance = - (meters) ? start_radius : ns_resolution * start_cells; - - step_radius = atof(par_multi_step->answer); - step_cells = - meters ? (int)(step_radius / max_resolution) : step_radius; - step_distance = - (meters) ? step_radius : ns_resolution * step_cells; - if (step_distance < ns_resolution) - G_fatal_error(_("For multiresolution mode step must be greater than or equal to resolution of one cell")); - - if (G_legal_filename(par_multi_prefix->answer) < 0 || - strlen(par_multi_prefix->answer) > 19) - G_fatal_error(_("<%s> is an incorrect prefix"), - par_multi_prefix->answer); - strcpy(prefix, par_multi_prefix->answer); - strcat(prefix, "_"); - num_of_steps = (int)ceil(search_distance / step_distance); - } /* end multires preparation */ - - /* print information about distances */ - G_verbose_message("Search distance m: %f, cells: %d", search_distance, - search_cells); - G_verbose_message("Skip distance m: %f, cells: %d", skip_distance, - skip_cells); - G_verbose_message("Flat threshold distance m: %f, height: %f", flat_distance, - flat_threshold_height); - G_verbose_message("%s version", (extended) ? "Extended" : "Basic"); - if (multires) { - G_verbose_message - ("Multiresolution mode: search start at: m: %f, cells: %d", - start_distance, start_cells); - G_verbose_message - ("Multiresolution mode: search step is: m: %f, number of steps %d", - step_distance, num_of_steps); - G_verbose_message("Prefix for output: %s", prefix); - } + { /* calculate parameters */ + int num_outputs = 0; + double search_radius, skip_radius, start_radius, step_radius; + double ns_resolution; + + multires = (par_multi_prefix->answer) ? 1 : 0; + for (i = 1; i < io_size; ++i) /* check for outputs */ + if (opt_output[i]->answer) { + if (G_legal_filename(opt_output[i]->answer) < 0) + G_fatal_error(_("<%s> is an illegal file name"), + opt_output[i]->answer); + num_outputs++; + } + if (!num_outputs && !multires) + G_fatal_error(_("At least one output is required, e.g. %s"), + opt_output[o_forms]->key); + + meters = (flag_units->answer != 0); + extended = (flag_extended->answer != 0); + nrows = Rast_window_rows(); + ncols = Rast_window_cols(); + Rast_get_window(&window); + G_begin_distance_calculations(); + + if (G_projection() == PROJECTION_LL) { /* for LL max_res should be NS */ + ns_resolution = + G_distance(0, Rast_row_to_northing(0, &window), 0, + Rast_row_to_northing(1, &window)); + max_resolution = ns_resolution; + } + else { + max_resolution = MAX(window.ns_res, window.ew_res); /* max_resolution MORE meters per cell */ + ns_resolution = window.ns_res; + } + + /* search distance */ + search_radius = atof(par_search_radius->answer); + search_cells = + meters ? (int)(search_radius / max_resolution) : search_radius; + if (search_cells < 1) + G_fatal_error(_("Search radius size must cover at least 1 cell")); + row_radius_size = + meters ? ceil(search_radius / ns_resolution) : search_radius; + row_buffer_size = row_radius_size * 2 + 1; + search_distance = + (meters) ? search_radius : ns_resolution * search_cells; + + /* skip distance */ + skip_radius = atof(par_skip_radius->answer); + skip_cells = + meters ? (int)(skip_radius / max_resolution) : skip_radius; + if (skip_cells >= search_cells) + G_fatal_error(_("Skip radius size must be at least 1 cell lower than radius")); + skip_distance = (meters) ? skip_radius : ns_resolution * skip_cells; + + /* flatness parameters */ + flat_threshold = atof(par_flat_threshold->answer); + if (flat_threshold <= 0.) + G_fatal_error(_("Flatness threshold must be grater than 0")); + flat_threshold = DEGREE2RAD(flat_threshold); + + flat_distance = atof(par_flat_distance->answer); + flat_distance = + (meters) ? flat_distance : ns_resolution * flat_distance; + flat_threshold_height = tan(flat_threshold) * flat_distance; + if ((flat_distance > 0 && flat_distance <= skip_distance) || + flat_distance >= search_distance) { + G_warning(_("Flatness distance should be between skip and search radius. Otherwise ignored")); + flat_distance = 0; + } + if (multires) { + start_radius = atof(par_multi_start->answer); + start_cells = + meters ? (int)(start_radius / max_resolution) : start_radius; + if (start_cells <= skip_cells) + start_cells = skip_cells + 1; + start_distance = + (meters) ? start_radius : ns_resolution * start_cells; + + step_radius = atof(par_multi_step->answer); + step_cells = + meters ? (int)(step_radius / max_resolution) : step_radius; + step_distance = + (meters) ? step_radius : ns_resolution * step_cells; + if (step_distance < ns_resolution) + G_fatal_error(_("For multiresolution mode step must be greater than or equal to resolution of one cell")); + + if (G_legal_filename(par_multi_prefix->answer) < 0 || + strlen(par_multi_prefix->answer) > 19) + G_fatal_error(_("<%s> is an incorrect prefix"), + par_multi_prefix->answer); + strcpy(prefix, par_multi_prefix->answer); + strcat(prefix, "_"); + num_of_steps = (int)ceil(search_distance / step_distance); + } /* end multires preparation */ + + /* print information about distances */ + G_verbose_message("Search distance m: %f, cells: %d", search_distance, + search_cells); + G_verbose_message("Skip distance m: %f, cells: %d", skip_distance, + skip_cells); + G_verbose_message("Flat threshold distance m: %f, height: %f", + flat_distance, flat_threshold_height); + G_verbose_message("%s version", (extended) ? "Extended" : "Basic"); + if (multires) { + G_verbose_message + ("Multiresolution mode: search start at: m: %f, cells: %d", + start_distance, start_cells); + G_verbose_message + ("Multiresolution mode: search step is: m: %f, number of steps %d", + step_distance, num_of_steps); + G_verbose_message("Prefix for output: %s", prefix); + } } /* generate global ternary codes */ for (i = 0; i < 6561; ++i) - global_ternary_codes[i] = ternary_rotate(i); + global_ternary_codes[i] = ternary_rotate(i); /* open DEM */ strcpy(elevation.elevname, opt_input->answer); open_map(&elevation); if (!multires) { - PATTERN *pattern; - PATTERN patterns[4]; - void *pointer_buf; - double search_dist = search_distance; - double skip_dist = skip_distance; - double flat_dist = flat_distance; - double area_of_octagon = - 4 * (search_distance * search_distance) * sin(DEGREE2RAD(45.)); - - cell_step = 1; - /* prepare outputs */ - for (i = 1; i < io_size; ++i) - if (opt_output[i]->answer) { - rasters[i].fd = - Rast_open_new(opt_output[i]->answer, - rasters[i].out_data_type); - rasters[i].buffer = - Rast_allocate_buf(rasters[i].out_data_type); - } - - /* main loop */ - for (row = 0; row < nrows; ++row) { - G_percent(row, nrows, 2); - cur_row = (row < row_radius_size) ? row : - ((row >= - nrows - row_radius_size - 1) ? row_buffer_size - (nrows - - row - - 1) : - row_radius_size); - - if (row > (row_radius_size) && - row < nrows - (row_radius_size + 1)) - shift_buffers(row); - for (col = 0; col < ncols; ++col) { - /* on borders forms ussualy are innatural. */ - if (row < (skip_cells + 1) || row > nrows - (skip_cells + 2) - || col < (skip_cells + 1) || - col > ncols - (skip_cells + 2) || - Rast_is_f_null_value(&elevation.elev[cur_row][col])) { - /* set outputs to NULL and do nothing if source value is null or border */ - for (i = 1; i < io_size; ++i) - if (opt_output[i]->answer) { - pointer_buf = rasters[i].buffer; - switch (rasters[i].out_data_type) { - case CELL_TYPE: - Rast_set_c_null_value(&((CELL *) pointer_buf) - [col], 1); - break; - case FCELL_TYPE: - Rast_set_f_null_value(&((FCELL *) pointer_buf) - [col], 1); - break; - case DCELL_TYPE: - Rast_set_d_null_value(&((DCELL *) pointer_buf) - [col], 1); - break; - default: - G_fatal_error(_("Unknown output data type")); - } - } - continue; - } /* end null value */ - { - int cur_form, small_form; - - search_distance = search_dist; - skip_distance = skip_dist; - flat_distance = flat_dist; - pattern_size = - calc_pattern(&patterns[0], row, cur_row, col); - pattern = &patterns[0]; - cur_form = - determine_form(pattern->num_negatives, - pattern->num_positives); - - /* correction of forms */ - if (extended && search_distance > 10 * max_resolution) { - /* 1) remove extensive innatural forms: ridges, peaks, shoulders and footslopes */ - if ((cur_form == 4 || cur_form == 8 || cur_form == 2 - || cur_form == 3)) { - search_distance = - (search_dist / 2. < - 4 * max_resolution) ? 4 * - max_resolution : search_dist / 2.; - skip_distance = 0; - flat_distance = 0; - pattern_size = - calc_pattern(&patterns[1], row, cur_row, col); - pattern = &patterns[1]; - small_form = - determine_form(pattern->num_negatives, - pattern->num_positives); - if (cur_form == 4 || cur_form == 8) - cur_form = (small_form == 1) ? 1 : cur_form; - if (cur_form == 2 || cur_form == 3) - cur_form = small_form; - } - /* 3) Depressions */ - - } /* end of correction */ - pattern = &patterns[0]; - if (opt_output[o_forms]->answer) - ((CELL *) rasters[o_forms].buffer)[col] = cur_form; - } - - if (opt_output[o_ternary]->answer) - ((CELL *) rasters[o_ternary].buffer)[col] = - determine_ternary(pattern->pattern); - if (opt_output[o_positive]->answer) - ((CELL *) rasters[o_positive].buffer)[col] = - rotate(pattern->positives); - if (opt_output[o_negative]->answer) - ((CELL *) rasters[o_negative].buffer)[col] = - rotate(pattern->negatives); - if (opt_output[o_intensity]->answer) - ((FCELL *) rasters[o_intensity].buffer)[col] = - intensity(pattern->elevation, pattern_size); - if (opt_output[o_exposition]->answer) - ((FCELL *) rasters[o_exposition].buffer)[col] = - exposition(pattern->elevation); - if (opt_output[o_range]->answer) - ((FCELL *) rasters[o_range].buffer)[col] = - range(pattern->elevation); - if (opt_output[o_variance]->answer) - ((FCELL *) rasters[o_variance].buffer)[col] = - variance(pattern->elevation, pattern_size); - - /* used only for next four shape functions */ - if (opt_output[o_elongation]->answer || - opt_output[o_azimuth]->answer || - opt_output[o_extend]->answer || - opt_output[o_width]->answer) { - float azimuth, elongation, width; - - radial2cartesian(pattern); - shape(pattern, pattern_size, &azimuth, &elongation, - &width); - if (opt_output[o_azimuth]->answer) - ((FCELL *) rasters[o_azimuth].buffer)[col] = azimuth; - if (opt_output[o_elongation]->answer) - ((FCELL *) rasters[o_elongation].buffer)[col] = - elongation; - if (opt_output[o_width]->answer) - ((FCELL *) rasters[o_width].buffer)[col] = width; - } - if (opt_output[o_extend]->answer) - ((FCELL *) rasters[o_extend].buffer)[col] = - extends(pattern, pattern_size) / area_of_octagon; - - } /* end for col */ - - /* write existing outputs */ - for (i = 1; i < io_size; ++i) - if (opt_output[i]->answer) - Rast_put_row(rasters[i].fd, rasters[i].buffer, - rasters[i].out_data_type); - } - G_percent(row, nrows, 2); /* end main loop */ - - /* finish and close */ - free_map(elevation.elev, row_buffer_size + 1); - for (i = 1; i < io_size; ++i) - if (opt_output[i]->answer) { - G_free(rasters[i].buffer); - Rast_close(rasters[i].fd); - Rast_short_history(opt_output[i]->answer, "raster", &history); - Rast_command_history(&history); - Rast_write_history(opt_output[i]->answer, &history); - } - - if (opt_output[o_forms]->answer) - write_form_cat_colors(opt_output[o_forms]->answer, ccolors); - if (opt_output[o_intensity]->answer) - write_contrast_colors(opt_output[o_intensity]->answer); - if (opt_output[o_exposition]->answer) - write_contrast_colors(opt_output[o_exposition]->answer); - if (opt_output[o_range]->answer) - write_contrast_colors(opt_output[o_range]->answer); - - G_done_msg(" "); - exit(EXIT_SUCCESS); - } /* end of multiresolution */ + PATTERN *pattern; + PATTERN patterns[4]; + void *pointer_buf; + double search_dist = search_distance; + double skip_dist = skip_distance; + double flat_dist = flat_distance; + double area_of_octagon = + 4 * (search_distance * search_distance) * sin(DEGREE2RAD(45.)); + + cell_step = 1; + /* prepare outputs */ + for (i = 1; i < io_size; ++i) + if (opt_output[i]->answer) { + rasters[i].fd = + Rast_open_new(opt_output[i]->answer, + rasters[i].out_data_type); + rasters[i].buffer = + Rast_allocate_buf(rasters[i].out_data_type); + } + + /* main loop */ + for (row = 0; row < nrows; ++row) { + G_percent(row, nrows, 2); + cur_row = (row < row_radius_size) ? row : + ((row >= + nrows - row_radius_size - 1) ? row_buffer_size - (nrows - + row - + 1) : + row_radius_size); + + if (row > (row_radius_size) && + row < nrows - (row_radius_size + 1)) + shift_buffers(row); + for (col = 0; col < ncols; ++col) { + /* on borders forms ussualy are innatural. */ + if (row < (skip_cells + 1) || row > nrows - (skip_cells + 2) + || col < (skip_cells + 1) || + col > ncols - (skip_cells + 2) || + Rast_is_f_null_value(&elevation.elev[cur_row][col])) { + /* set outputs to NULL and do nothing if source value is null or border */ + for (i = 1; i < io_size; ++i) + if (opt_output[i]->answer) { + pointer_buf = rasters[i].buffer; + switch (rasters[i].out_data_type) { + case CELL_TYPE: + Rast_set_c_null_value(&((CELL *) pointer_buf) + [col], 1); + break; + case FCELL_TYPE: + Rast_set_f_null_value(&((FCELL *) pointer_buf) + [col], 1); + break; + case DCELL_TYPE: + Rast_set_d_null_value(&((DCELL *) pointer_buf) + [col], 1); + break; + default: + G_fatal_error(_("Unknown output data type")); + } + } + continue; + } /* end null value */ + { + int cur_form, small_form; + + search_distance = search_dist; + skip_distance = skip_dist; + flat_distance = flat_dist; + pattern_size = + calc_pattern(&patterns[0], row, cur_row, col); + pattern = &patterns[0]; + cur_form = + determine_form(pattern->num_negatives, + pattern->num_positives); + + /* correction of forms */ + if (extended && search_distance > 10 * max_resolution) { + /* 1) remove extensive innatural forms: ridges, peaks, shoulders and footslopes */ + if ((cur_form == 4 || cur_form == 8 || cur_form == 2 + || cur_form == 3)) { + search_distance = + (search_dist / 2. < + 4 * max_resolution) ? 4 * + max_resolution : search_dist / 2.; + skip_distance = 0; + flat_distance = 0; + pattern_size = + calc_pattern(&patterns[1], row, cur_row, col); + pattern = &patterns[1]; + small_form = + determine_form(pattern->num_negatives, + pattern->num_positives); + if (cur_form == 4 || cur_form == 8) + cur_form = (small_form == 1) ? 1 : cur_form; + if (cur_form == 2 || cur_form == 3) + cur_form = small_form; + } + /* 3) Depressions */ + + } /* end of correction */ + pattern = &patterns[0]; + if (opt_output[o_forms]->answer) + ((CELL *) rasters[o_forms].buffer)[col] = cur_form; + } + + if (opt_output[o_ternary]->answer) + ((CELL *) rasters[o_ternary].buffer)[col] = + determine_ternary(pattern->pattern); + if (opt_output[o_positive]->answer) + ((CELL *) rasters[o_positive].buffer)[col] = + rotate(pattern->positives); + if (opt_output[o_negative]->answer) + ((CELL *) rasters[o_negative].buffer)[col] = + rotate(pattern->negatives); + if (opt_output[o_intensity]->answer) + ((FCELL *) rasters[o_intensity].buffer)[col] = + intensity(pattern->elevation, pattern_size); + if (opt_output[o_exposition]->answer) + ((FCELL *) rasters[o_exposition].buffer)[col] = + exposition(pattern->elevation); + if (opt_output[o_range]->answer) + ((FCELL *) rasters[o_range].buffer)[col] = + range(pattern->elevation); + if (opt_output[o_variance]->answer) + ((FCELL *) rasters[o_variance].buffer)[col] = + variance(pattern->elevation, pattern_size); + + /* used only for next four shape functions */ + if (opt_output[o_elongation]->answer || + opt_output[o_azimuth]->answer || + opt_output[o_extend]->answer || + opt_output[o_width]->answer) { + float azimuth, elongation, width; + + radial2cartesian(pattern); + shape(pattern, pattern_size, &azimuth, &elongation, + &width); + if (opt_output[o_azimuth]->answer) + ((FCELL *) rasters[o_azimuth].buffer)[col] = azimuth; + if (opt_output[o_elongation]->answer) + ((FCELL *) rasters[o_elongation].buffer)[col] = + elongation; + if (opt_output[o_width]->answer) + ((FCELL *) rasters[o_width].buffer)[col] = width; + } + if (opt_output[o_extend]->answer) + ((FCELL *) rasters[o_extend].buffer)[col] = + extends(pattern, pattern_size) / area_of_octagon; + + } /* end for col */ + + /* write existing outputs */ + for (i = 1; i < io_size; ++i) + if (opt_output[i]->answer) + Rast_put_row(rasters[i].fd, rasters[i].buffer, + rasters[i].out_data_type); + } + G_percent(row, nrows, 2); /* end main loop */ + + /* finish and close */ + free_map(elevation.elev, row_buffer_size + 1); + for (i = 1; i < io_size; ++i) + if (opt_output[i]->answer) { + G_free(rasters[i].buffer); + Rast_close(rasters[i].fd); + Rast_short_history(opt_output[i]->answer, "raster", &history); + Rast_command_history(&history); + Rast_write_history(opt_output[i]->answer, &history); + } + + if (opt_output[o_forms]->answer) + write_form_cat_colors(opt_output[o_forms]->answer, ccolors); + if (opt_output[o_intensity]->answer) + write_contrast_colors(opt_output[o_intensity]->answer); + if (opt_output[o_exposition]->answer) + write_contrast_colors(opt_output[o_exposition]->answer); + if (opt_output[o_range]->answer) + write_contrast_colors(opt_output[o_range]->answer); + + G_done_msg(" "); + exit(EXIT_SUCCESS); + } /* end of multiresolution */ if (multires) { - PATTERN *multi_patterns; - MULTI multiple_output[5]; /* ten form maps + all forms */ - char *postfixes[] = { "scale_300", "scale_100", "scale_50", "scale_20", "scale_10" }; /* in pixels */ - num_of_steps = 5; - multi_patterns = G_malloc(num_of_steps * sizeof(PATTERN)); - /* prepare outputs */ - for (i = 0; i < 5; ++i) { - multiple_output[i].forms_buffer = Rast_allocate_buf(CELL_TYPE); - strcpy(multiple_output[i].name, prefix); - strcat(multiple_output[i].name, postfixes[i]); - multiple_output[i].fd = - Rast_open_new(multiple_output[i].name, CELL_TYPE); - } - - /* main loop */ - for (row = 0; row < nrows; ++row) { - G_percent(row, nrows, 2); - cur_row = (row < row_radius_size) ? row : - ((row >= - nrows - row_radius_size - 1) ? row_buffer_size - (nrows - - row - - 1) : - row_radius_size); - - if (row > (row_radius_size) && - row < nrows - (row_radius_size + 1)) - shift_buffers(row); - for (col = 0; col < ncols; ++col) { - if (row < (skip_cells + 1) || row > nrows - (skip_cells + 2) - || col < (skip_cells + 1) || - col > ncols - (skip_cells + 2) || - Rast_is_f_null_value(&elevation.elev[cur_row][col])) { - for (i = 0; i < num_of_steps; ++i) - Rast_set_c_null_value(&multiple_output[i]. - forms_buffer[col], 1); - continue; - } - cell_step = 10; - calc_pattern(&multi_patterns[0], row, cur_row, col); - } - - for (i = 0; i < num_of_steps; ++i) - Rast_put_row(multiple_output[i].fd, - multiple_output[i].forms_buffer, CELL_TYPE); - - } - G_percent(row, nrows, 2); /* end main loop */ - - for (i = 0; i < num_of_steps; ++i) { - G_free(multiple_output[i].forms_buffer); - Rast_close(multiple_output[i].fd); - Rast_short_history(multiple_output[i].name, "raster", &history); - Rast_command_history(&history); - Rast_write_history(multiple_output[i].name, &history); - } - G_message("Multiresolution Done!"); - exit(EXIT_SUCCESS); + PATTERN *multi_patterns; + MULTI multiple_output[5]; /* ten form maps + all forms */ + char *postfixes[] = { "scale_300", "scale_100", "scale_50", "scale_20", "scale_10" }; /* in pixels */ + num_of_steps = 5; + multi_patterns = G_malloc(num_of_steps * sizeof(PATTERN)); + /* prepare outputs */ + for (i = 0; i < 5; ++i) { + multiple_output[i].forms_buffer = Rast_allocate_buf(CELL_TYPE); + strcpy(multiple_output[i].name, prefix); + strcat(multiple_output[i].name, postfixes[i]); + multiple_output[i].fd = + Rast_open_new(multiple_output[i].name, CELL_TYPE); + } + + /* main loop */ + for (row = 0; row < nrows; ++row) { + G_percent(row, nrows, 2); + cur_row = (row < row_radius_size) ? row : + ((row >= + nrows - row_radius_size - 1) ? row_buffer_size - (nrows - + row - + 1) : + row_radius_size); + + if (row > (row_radius_size) && + row < nrows - (row_radius_size + 1)) + shift_buffers(row); + for (col = 0; col < ncols; ++col) { + if (row < (skip_cells + 1) || row > nrows - (skip_cells + 2) + || col < (skip_cells + 1) || + col > ncols - (skip_cells + 2) || + Rast_is_f_null_value(&elevation.elev[cur_row][col])) { + for (i = 0; i < num_of_steps; ++i) + Rast_set_c_null_value(&multiple_output[i].forms_buffer + [col], 1); + continue; + } + cell_step = 10; + calc_pattern(&multi_patterns[0], row, cur_row, col); + } + + for (i = 0; i < num_of_steps; ++i) + Rast_put_row(multiple_output[i].fd, + multiple_output[i].forms_buffer, CELL_TYPE); + + } + G_percent(row, nrows, 2); /* end main loop */ + + for (i = 0; i < num_of_steps; ++i) { + G_free(multiple_output[i].forms_buffer); + Rast_close(multiple_output[i].fd); + Rast_short_history(multiple_output[i].name, "raster", &history); + Rast_command_history(&history); + Rast_write_history(multiple_output[i].name, &history); + } + G_message("Multiresolution Done!"); + exit(EXIT_SUCCESS); } } diff --git a/raster/r.geomorphon/memory.c b/raster/r.geomorphon/memory.c index 5e764cf26e6..4e7b1ea3c3f 100644 --- a/raster/r.geomorphon/memory.c +++ b/raster/r.geomorphon/memory.c @@ -11,25 +11,25 @@ int open_map(MAPS * rast) mapset = (char *)G_find_raster2(rast->elevname, ""); if (mapset == NULL) - G_fatal_error(_("Raster map <%s> not found"), rast->elevname); + G_fatal_error(_("Raster map <%s> not found"), rast->elevname); rast->fd = Rast_open_old(rast->elevname, mapset); Rast_get_cellhd(rast->elevname, mapset, &cellhd); rast->raster_type = Rast_map_type(rast->elevname, mapset); if (window.ew_res < cellhd.ew_res || window.ns_res < cellhd.ns_res) - G_warning(_("Region resolution shoudn't be lesser than map %s resolution. Run g.region raster=%s to set proper resolution"), - rast->elevname, rast->elevname); + G_warning(_("Region resolution shoudn't be lesser than map %s resolution. Run g.region raster=%s to set proper resolution"), + rast->elevname, rast->elevname); tmp_buf = Rast_allocate_buf(rast->raster_type); rast->elev = (FCELL **) G_malloc((row_buffer_size + 1) * sizeof(FCELL *)); for (row = 0; row < row_buffer_size + 1; ++row) { - rast->elev[row] = Rast_allocate_buf(FCELL_TYPE); - Rast_get_row(rast->fd, tmp_buf, row, rast->raster_type); - for (col = 0; col < ncols; ++col) - get_cell(col, rast->elev[row], tmp_buf, rast->raster_type); - } /* end elev */ + rast->elev[row] = Rast_allocate_buf(FCELL_TYPE); + Rast_get_row(rast->fd, tmp_buf, row, rast->raster_type); + for (col = 0; col < ncols; ++col) + get_cell(col, rast->elev[row], tmp_buf, rast->raster_type); + } /* end elev */ G_free(tmp_buf); return 0; @@ -41,25 +41,25 @@ int get_cell(int col, float *buf_row, void *buf, RASTER_MAP_TYPE raster_type) switch (raster_type) { case CELL_TYPE: - if (Rast_is_null_value(&((CELL *) buf)[col], CELL_TYPE)) - Rast_set_f_null_value(&buf_row[col], 1); - else - buf_row[col] = (FCELL) ((CELL *) buf)[col]; - break; + if (Rast_is_null_value(&((CELL *) buf)[col], CELL_TYPE)) + Rast_set_f_null_value(&buf_row[col], 1); + else + buf_row[col] = (FCELL) ((CELL *) buf)[col]; + break; case FCELL_TYPE: - if (Rast_is_null_value(&((FCELL *) buf)[col], FCELL_TYPE)) - Rast_set_f_null_value(&buf_row[col], 1); - else - buf_row[col] = (FCELL) ((FCELL *) buf)[col]; - break; + if (Rast_is_null_value(&((FCELL *) buf)[col], FCELL_TYPE)) + Rast_set_f_null_value(&buf_row[col], 1); + else + buf_row[col] = (FCELL) ((FCELL *) buf)[col]; + break; case DCELL_TYPE: - if (Rast_is_null_value(&((DCELL *) buf)[col], DCELL_TYPE)) - Rast_set_f_null_value(&buf_row[col], 1); - else - buf_row[col] = (FCELL) ((DCELL *) buf)[col]; - break; + if (Rast_is_null_value(&((DCELL *) buf)[col], DCELL_TYPE)) + Rast_set_f_null_value(&buf_row[col], 1); + else + buf_row[col] = (FCELL) ((DCELL *) buf)[col]; + break; } return 0; @@ -76,15 +76,15 @@ int shift_buffers(int row) tmp_elev_buf = elevation.elev[0]; for (i = 1; i < row_buffer_size + 1; ++i) - elevation.elev[i - 1] = elevation.elev[i]; + elevation.elev[i - 1] = elevation.elev[i]; elevation.elev[row_buffer_size] = tmp_elev_buf; Rast_get_row(elevation.fd, tmp_buf, row + row_radius_size + 1, - elevation.raster_type); + elevation.raster_type); for (col = 0; col < ncols; ++col) - get_cell(col, elevation.elev[row_buffer_size], tmp_buf, - elevation.raster_type); + get_cell(col, elevation.elev[row_buffer_size], tmp_buf, + elevation.raster_type); G_free(tmp_buf); return 0; @@ -95,7 +95,7 @@ int free_map(FCELL ** map, int n) int i; for (i = 0; i < n; ++i) - G_free(map[i]); + G_free(map[i]); G_free(map); return 0; } @@ -109,15 +109,15 @@ int write_form_cat_colors(char *raster, CATCOLORS * ccolors) Rast_init_colors(&colors); for (i = 1; i < CNT; ++i) - Rast_add_color_rule(&ccolors[i].cat, ccolors[i].r, ccolors[i].g, - ccolors[i].b, &ccolors[i].cat, ccolors[i].r, - ccolors[i].g, ccolors[i].b, &colors, CELL_TYPE); + Rast_add_color_rule(&ccolors[i].cat, ccolors[i].r, ccolors[i].g, + ccolors[i].b, &ccolors[i].cat, ccolors[i].r, + ccolors[i].g, ccolors[i].b, &colors, CELL_TYPE); Rast_write_colors(raster, G_mapset(), &colors); Rast_free_colors(&colors); Rast_init_cats("Forms", &cats); for (i = 1; i < CNT; ++i) - Rast_set_cat(&ccolors[i].cat, &ccolors[i].cat, ccolors[i].label, - &cats, CELL_TYPE); + Rast_set_cat(&ccolors[i].cat, &ccolors[i].cat, ccolors[i].label, + &cats, CELL_TYPE); Rast_write_cats(raster, &cats); Rast_free_cats(&cats); return 0; @@ -126,28 +126,29 @@ int write_form_cat_colors(char *raster, CATCOLORS * ccolors) int write_contrast_colors(char *raster) { struct Colors colors; + /* struct Categories cats; */ - FCOLORS fcolors[9] = { /* colors for positive openness */ - {-2500, 0, 0, 50, NULL}, - {-100, 0, 0, 56, NULL}, - {-15, 0, 56, 128, NULL}, - {-3, 0, 128, 255, NULL}, - {0, 255, 255, 255, NULL}, - {3, 255, 128, 0, NULL}, - {15, 128, 56, 0, NULL}, - {100, 56, 0, 0, NULL}, - {2500, 50, 0, 0, NULL} + FCOLORS fcolors[9] = { /* colors for positive openness */ + {-2500, 0, 0, 50, NULL}, + {-100, 0, 0, 56, NULL}, + {-15, 0, 56, 128, NULL}, + {-3, 0, 128, 255, NULL}, + {0, 255, 255, 255, NULL}, + {3, 255, 128, 0, NULL}, + {15, 128, 56, 0, NULL}, + {100, 56, 0, 0, NULL}, + {2500, 50, 0, 0, NULL} }; int i; Rast_init_colors(&colors); for (i = 0; i < 8; ++i) - Rast_add_d_color_rule(&fcolors[i].cat, fcolors[i].r, fcolors[i].g, - fcolors[i].b, &fcolors[i + 1].cat, - fcolors[i + 1].r, fcolors[i + 1].g, - fcolors[i + 1].b, &colors); + Rast_add_d_color_rule(&fcolors[i].cat, fcolors[i].r, fcolors[i].g, + fcolors[i].b, &fcolors[i + 1].cat, + fcolors[i + 1].r, fcolors[i + 1].g, + fcolors[i + 1].b, &colors); Rast_write_colors(raster, G_mapset(), &colors); Rast_free_colors(&colors); /* diff --git a/raster/r.geomorphon/multires.c b/raster/r.geomorphon/multires.c index cc3a21c4aba..da2f9c5a46c 100644 --- a/raster/r.geomorphon/multires.c +++ b/raster/r.geomorphon/multires.c @@ -8,14 +8,14 @@ int pattern_matching(int *pattern) int sign = -1; for (i = 0, n = 1; i < 8; i++, n *= 2) - binary += (pattern[i] == sign) ? n : 0; + binary += (pattern[i] == sign) ? n : 0; /* rotate */ for (i = 0; i < 8; ++i) { - if ((i &= 7) == 0) - test = binary; - else - test = (binary << i) | (binary >> (8 - i)); - result = (result < test) ? result : test; + if ((i &= 7) == 0) + test = binary; + else + test = (binary << i) | (binary >> (8 - i)); + result = (result < test) ? result : test; } return ((result & source) == source) ? 1 : 0; diff --git a/raster/r.geomorphon/pattern.c b/raster/r.geomorphon/pattern.c index 969158c5c7c..1693ff42ae6 100644 --- a/raster/r.geomorphon/pattern.c +++ b/raster/r.geomorphon/pattern.c @@ -28,118 +28,118 @@ int calc_pattern(PATTERN * pattern, int row, int cur_row, int col) pattern->negatives = 0; for (i = 0; i < 8; ++i) { - /* reset patterns */ - pattern->pattern[i] = 0; - pattern->elevation[i] = 0.; - pattern->distance[i] = 0.; - j = skip_cells + 1; - zenith_angle = -(PI2); - nadir_angle = PI2; + /* reset patterns */ + pattern->pattern[i] = 0; + pattern->elevation[i] = 0.; + pattern->distance[i] = 0.; + j = skip_cells + 1; + zenith_angle = -(PI2); + nadir_angle = PI2; - if (cur_row + j * nextr[i] < 0 || - cur_row + j * nextr[i] > row_buffer_size - 1 || - col + j * nextc[i] < 0 || col + j * nextc[i] > ncols - 1) - continue; /* border: current cell is on the end of DEM */ - if (Rast_is_f_null_value - (&elevation.elev[cur_row + nextr[i]][col + nextc[i]])) - continue; /* border: next value is null, line-of-sight does not exists */ - pattern_size++; /* line-of-sight exists, continue calculate visibility */ + if (cur_row + j * nextr[i] < 0 || + cur_row + j * nextr[i] > row_buffer_size - 1 || + col + j * nextc[i] < 0 || col + j * nextc[i] > ncols - 1) + continue; /* border: current cell is on the end of DEM */ + if (Rast_is_f_null_value + (&elevation.elev[cur_row + nextr[i]][col + nextc[i]])) + continue; /* border: next value is null, line-of-sight does not exists */ + pattern_size++; /* line-of-sight exists, continue calculate visibility */ - target_northing = - Rast_row_to_northing(row + j * nextr[i] + 0.5, &window); - target_easting = - Rast_col_to_easting(col + j * nextc[i] + 0.5, &window); - cur_distance = - G_distance(cur_easting, cur_northing, target_easting, - target_northing); + target_northing = + Rast_row_to_northing(row + j * nextr[i] + 0.5, &window); + target_easting = + Rast_col_to_easting(col + j * nextc[i] + 0.5, &window); + cur_distance = + G_distance(cur_easting, cur_northing, target_easting, + target_northing); - while (cur_distance < search_distance) { - if (cur_row + j * nextr[i] < 0 || - cur_row + j * nextr[i] > row_buffer_size - 1 || - col + j * nextc[i] < 0 || col + j * nextc[i] > ncols - 1) - break; /* reached end of DEM (cols) or buffer (rows) */ + while (cur_distance < search_distance) { + if (cur_row + j * nextr[i] < 0 || + cur_row + j * nextr[i] > row_buffer_size - 1 || + col + j * nextc[i] < 0 || col + j * nextc[i] > ncols - 1) + break; /* reached end of DEM (cols) or buffer (rows) */ - height = - elevation.elev[cur_row + j * nextr[i]][col + j * nextc[i]] - - center_height; - angle = atan2(height, cur_distance); + height = + elevation.elev[cur_row + j * nextr[i]][col + j * nextc[i]] - + center_height; + angle = atan2(height, cur_distance); - if (angle > zenith_angle) { - zenith_angle = angle; - zenith_height = height; - zenith_distance = cur_distance; - } - if (angle < nadir_angle) { - nadir_angle = angle; - nadir_height = height; - nadir_distance = cur_distance; - } - j += cell_step; - /* j++; */ /* go to next cell */ - target_northing = - Rast_row_to_northing(row + j * nextr[i] + 0.5, &window); - target_easting = - Rast_col_to_easting(col + j * nextc[i] + 0.5, &window); - cur_distance = - G_distance(cur_easting, cur_northing, target_easting, - target_northing); - } /* end line of sight */ + if (angle > zenith_angle) { + zenith_angle = angle; + zenith_height = height; + zenith_distance = cur_distance; + } + if (angle < nadir_angle) { + nadir_angle = angle; + nadir_height = height; + nadir_distance = cur_distance; + } + j += cell_step; + /* j++; *//* go to next cell */ + target_northing = + Rast_row_to_northing(row + j * nextr[i] + 0.5, &window); + target_easting = + Rast_col_to_easting(col + j * nextc[i] + 0.5, &window); + cur_distance = + G_distance(cur_easting, cur_northing, target_easting, + target_northing); + } /* end line of sight */ - /* original paper version */ - /* zenith_angle=PI2-zenith_angle; - nadir_angle=PI2+nadir_angle; - if(fabs(zenith_angle-nadir_angle) > flat_threshold) { - if((nadir_angle-zenith_angle) > 0) { - patterns->pattern[i]=1; - patterns->elevation[i]=nadir_height; - patterns->distance[i]=nadir_distance; - patterns->num_positives++; - } else { - patterns->pattern[i]=-1; - patterns->elevation[i]=zenith_height; - patterns->distance[i]=zenith_distance; - patterns->num_negatives++; - } - } else { - patterns->distance[i]=search_distance; - } - */ - /* this is used to lower flat threshold if distance exceed flat_distance parameter */ - zenith_threshold = (flat_distance > 0 && - flat_distance < - zenith_distance) ? atan2(flat_threshold_height, - zenith_distance) : - flat_threshold; - nadir_threshold = (flat_distance > 0 && - flat_distance < - nadir_distance) ? atan2(flat_threshold_height, - nadir_distance) : - flat_threshold; + /* original paper version */ + /* zenith_angle=PI2-zenith_angle; + nadir_angle=PI2+nadir_angle; + if(fabs(zenith_angle-nadir_angle) > flat_threshold) { + if((nadir_angle-zenith_angle) > 0) { + patterns->pattern[i]=1; + patterns->elevation[i]=nadir_height; + patterns->distance[i]=nadir_distance; + patterns->num_positives++; + } else { + patterns->pattern[i]=-1; + patterns->elevation[i]=zenith_height; + patterns->distance[i]=zenith_distance; + patterns->num_negatives++; + } + } else { + patterns->distance[i]=search_distance; + } + */ + /* this is used to lower flat threshold if distance exceed flat_distance parameter */ + zenith_threshold = (flat_distance > 0 && + flat_distance < + zenith_distance) ? atan2(flat_threshold_height, + zenith_distance) : + flat_threshold; + nadir_threshold = (flat_distance > 0 && + flat_distance < + nadir_distance) ? atan2(flat_threshold_height, + nadir_distance) : + flat_threshold; - if (zenith_angle > zenith_threshold) - pattern->positives += i; - if (nadir_angle < -nadir_threshold) - pattern->negatives += i; + if (zenith_angle > zenith_threshold) + pattern->positives += i; + if (nadir_angle < -nadir_threshold) + pattern->negatives += i; - if (fabs(zenith_angle) > zenith_threshold || - fabs(nadir_angle) > nadir_threshold) { - if (fabs(nadir_angle) < fabs(zenith_angle)) { - pattern->pattern[i] = 1; - pattern->elevation[i] = zenith_height; /* A CHANGE! */ - pattern->distance[i] = zenith_distance; - pattern->num_positives++; - } - if (fabs(nadir_angle) > fabs(zenith_angle)) { - pattern->pattern[i] = -1; - pattern->elevation[i] = nadir_height; /* A CHANGE! */ - pattern->distance[i] = nadir_distance; - pattern->num_negatives++; - } - } - else { - pattern->distance[i] = search_distance; - } + if (fabs(zenith_angle) > zenith_threshold || + fabs(nadir_angle) > nadir_threshold) { + if (fabs(nadir_angle) < fabs(zenith_angle)) { + pattern->pattern[i] = 1; + pattern->elevation[i] = zenith_height; /* A CHANGE! */ + pattern->distance[i] = zenith_distance; + pattern->num_positives++; + } + if (fabs(nadir_angle) > fabs(zenith_angle)) { + pattern->pattern[i] = -1; + pattern->elevation[i] = nadir_height; /* A CHANGE! */ + pattern->distance[i] = nadir_distance; + pattern->num_negatives++; + } + } + else { + pattern->distance[i] = search_distance; + } - } /* end for */ + } /* end for */ return pattern_size; }