[GRASS-SVN] r54918 - in grass-addons/grass7/raster: . r.geomorphon
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Feb 5 02:55:21 PST 2013
Author: madi
Date: 2013-02-05 02:55:21 -0800 (Tue, 05 Feb 2013)
New Revision: 54918
Added:
grass-addons/grass7/raster/r.geomorphon/
grass-addons/grass7/raster/r.geomorphon/Makefile
grass-addons/grass7/raster/r.geomorphon/geom.c
grass-addons/grass7/raster/r.geomorphon/geomorphon.png
grass-addons/grass7/raster/r.geomorphon/legend.png
grass-addons/grass7/raster/r.geomorphon/local_proto.h
grass-addons/grass7/raster/r.geomorphon/main.c
grass-addons/grass7/raster/r.geomorphon/memory.c
grass-addons/grass7/raster/r.geomorphon/multires.c
grass-addons/grass7/raster/r.geomorphon/pattern.c
grass-addons/grass7/raster/r.geomorphon/r.geomorphon.html
Log:
committing on behalf of Jarek
Added: grass-addons/grass7/raster/r.geomorphon/Makefile
===================================================================
--- grass-addons/grass7/raster/r.geomorphon/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.geomorphon/Makefile 2013-02-05 10:55:21 UTC (rev 54918)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.geomorphon
+
+LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB)
+DEPENDENCIES = $(RASTERDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Property changes on: grass-addons/grass7/raster/r.geomorphon/Makefile
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/grass7/raster/r.geomorphon/geom.c
===================================================================
--- grass-addons/grass7/raster/r.geomorphon/geom.c (rev 0)
+++ grass-addons/grass7/raster/r.geomorphon/geom.c 2013-02-05 10:55:21 UTC (rev 54918)
@@ -0,0 +1,252 @@
+#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 */
+/* DIRS in DEGREES from NORTH: 45,0,315,270,225,180,135,90 */
+
+unsigned int ternary_rotate (unsigned int value)
+{
+/* this function returns rotated and mirrored
+ * ternary code from any number
+ * function is used to create lookup table with original
+ * terrain patterns (6561) and its rotated and mirrored counterparts (498)*/
+
+unsigned char pattern[8];
+unsigned char rev_pattern[8];
+unsigned char tmp_pattern[8];
+unsigned char tmp_rev_pattern[8];
+unsigned int code=10000, tmp_code, rev_code=10000, tmp_rev_code;
+int power=1;
+int i,j,k;
+int res;
+
+for (i=0;i<8;i++) {
+ 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;
+}
+return code<rev_code ? code : rev_code;
+}
+
+int determine_form(int num_minus, int num_plus)
+{
+/* determine form according number of positives and negatives
+ * 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, __, __, __, __, __, __, __, __},
+ };
+
+/* legend:
+ FL, flat
+ PK, peak, summit
+ RI, ridge
+ SH, shoulder
+ CV, convex
+ SL, slope
+ CN, concave
+ FS, footslope
+ VL, valley
+ PT, pit, depression
+ __ error, impossible
+*/
+ return forms[num_minus][num_plus];
+}
+
+int determine_binary(int* pattern, int sign)
+{
+/* extract binary pattern for zenith (+) or nadir (-) from unrotated ternary pattern */
+int n, i;
+unsigned char binary=0, result=255, test=0;
+
+ for (i=0,n=1;i<8;i++,n*=2)
+ 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;
+ }
+return (int)result;
+}
+
+
+int rotate(unsigned char binary)
+{
+ int i;
+ unsigned char result=255, test=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;
+ }
+return (int)result;
+}
+
+int determine_ternary(int* pattern)
+{
+/* extract rotated and mirrored ternary pattern form unrotated ternary pattern */
+unsigned ternary_code=0;
+int power, i;
+
+ for(i=0,power=1;i<8;++i,power*=3)
+ternary_code+=(pattern[i]+1)*power;
+return global_ternary_codes[ternary_code];
+}
+
+float intensity(float* elevation, int pattern_size)
+{
+/* calculate relative elevation of the central cell against its visibility surround */
+float sum_elevation=0.;
+int i;
+
+ for (i=0;i<8;i++)
+sum_elevation-=elevation[i];
+
+return sum_elevation/(float)pattern_size;
+}
+
+float exposition(float* elevation)
+{
+/* calculate relative elevation of the central cell against its visibility */
+float max=0.;
+int i;
+ for (i=0;i<8;i++)
+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 */
+int i;
+
+for (i=0;i<8;i++) {
+ max=elevation[i]>max?elevation[i]:max;
+ min=elevation[i]<min?elevation[i]:min;
+}
+
+return max-min;
+}
+
+float variance(float* elevation, int pattern_size)
+{
+/* calculate relative variation of the visible neighbourhood of the cell */
+float sum_elevation=0;
+float variance=0;
+int i;
+
+ for (i=0;i<8;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]));
+
+return variance/(float)pattern_size;
+}
+
+int radial2cartesian(PATTERN* pattern)
+{
+/* this function converts radial coordinates of geomorphon
+ * (assuming center as 0,0) to cartezian coordinates
+ * with the begining in the central cell of geomorphon */
+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;
+}
+return 0;
+}
+
+float extends(PATTERN* pattern, int pattern_size)
+{
+int i,j;
+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]);
+}
+return fabs(area)/2.;
+}
+
+int shape(PATTERN* pattern, int pattern_size, float* azimuth, float* elongation, float* width)
+{
+/* calculates azimuth, elongation and width of geomorphon's polygon */
+int i;
+double avg_x=0, avg_y=0;
+double avg_x_y=0;
+double avg_x_square;
+double rx, ry;
+double sine, cosine;
+double result;
+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 /= (float)pattern_size;
+avg_x /= (float)pattern_size;
+avg_x_y /= (float)pattern_size;
+avg_x_square /= (float)pattern_size;
+result=(avg_x_y - avg_x * avg_y) / (avg_x_square - avg_x * avg_x);
+result=atan(result);
+*azimuth=(float)RAD2DEGREE(PI2-result);
+
+/* rotation */
+sine=sin(result);
+cosine=cos(result);
+for (i=0;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=(rxmax-rxmin);
+ry=(rymax-rymin);
+*elongation=rx>ry?(float)rx/ry:(float)ry/rx;
+*width=rx>ry?ry:rx;
+
+return 0;
+}
\ No newline at end of file
Property changes on: grass-addons/grass7/raster/r.geomorphon/geom.c
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/grass7/raster/r.geomorphon/geomorphon.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/raster/r.geomorphon/geomorphon.png
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/grass7/raster/r.geomorphon/legend.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/raster/r.geomorphon/legend.png
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/grass7/raster/r.geomorphon/local_proto.h
===================================================================
--- grass-addons/grass7/raster/r.geomorphon/local_proto.h (rev 0)
+++ grass-addons/grass7/raster/r.geomorphon/local_proto.h 2013-02-05 10:55:21 UTC (rev 54918)
@@ -0,0 +1,162 @@
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <grass/glocale.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+
+#ifdef MAIN
+# define GLOBAL
+#else
+# define GLOBAL extern
+#endif
+
+
+#ifndef PI2 /* PI/2 */
+ #define PI2 (2*atan(1))
+#endif
+
+#ifndef PI4 /* PI/4 */
+ #define PI4 (atan(1))
+#endif
+
+#ifndef PI
+ #define PI (4*atan(1))
+#endif
+
+#ifndef M2PI /* 2*PI */
+ #define M2PI (8*atan(1))
+#endif
+
+#ifndef GRADE
+ #define GRADE (atan(1)/45)
+#endif
+
+
+#ifndef PI2PERCENT
+ #define PI2PERCENT (50/atan(1))
+#endif
+
+#ifndef UNKNOWN
+ #define UNKNOWN -1
+#endif
+
+
+#define DEGREE2RAD(a) ((a)/(180/PI))
+#define RAD2DEGREE(a) ((a)*(180/PI))
+
+
+#define MAX(a,b) ((a) > (b) ? (a) : (b))
+#define MIN(a,b) ((a) < (b) ? (a) : (b))
+
+typedef char* STRING;
+
+typedef struct {
+ char elevname[150];
+ RASTER_MAP_TYPE raster_type;
+ FCELL **elev;
+ int fd; /* file descriptor */
+} MAPS;
+
+typedef struct { /* struct is used both for interface and output */
+ char* name;
+ int required;
+ char* description;
+ char* gui;
+ RASTER_MAP_TYPE out_data_type;
+ int fd;
+ void* buffer;
+} IO;
+
+typedef struct {
+ char name[100];
+ int fd;
+ CELL* forms_buffer;
+} MULTI;
+
+typedef struct {
+ int num_positives;
+ int num_negatives;
+ unsigned char positives;
+ unsigned char negatives;
+ int pattern[8];
+ float elevation[8];
+ double distance[8];
+ 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 */
+} FORMS;
+
+typedef struct {
+ int cat;
+ int r;
+ int g;
+ int b;
+ char* label;
+} CATCOLORS;
+
+typedef struct {
+ double cat;
+ int r;
+ int g;
+ int b;
+ char* label;
+} FCOLORS;
+
+/* variables */
+GLOBAL MAPS elevation;
+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 flat_threshold, flat_threshold_height;
+GLOBAL double H,V;
+GLOBAL struct Cell_head window;
+GLOBAL int cell_step;
+
+/* main */
+unsigned int global_ternary_codes[6562];
+
+/* memory */
+int open_map(MAPS* rast);
+int create_maps(void);
+int shift_buffers(int row);
+int get_cell (int col, float* buf_row, void* buf, RASTER_MAP_TYPE raster_type);
+int free_map(FCELL **map, int n);
+int write_form_cat_colors (char* raster, CATCOLORS* ccolors);
+
+/*geom */
+int calc_pattern(PATTERN* pattern, int row, int cur_row, int col);
+unsigned int ternary_rotate (unsigned int value);
+int determine_form(int num_plus, int num_minus);
+int determine_binary(int* pattern, int sign);
+int determine_ternary(int* pattern);
+int rotate(unsigned char binary);
+float intensity(float* elevation, int pattern_size);
+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 extends(PATTERN* pattern, int pattern_size);
+
+/* multires */
+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);
\ No newline at end of file
Property changes on: grass-addons/grass7/raster/r.geomorphon/local_proto.h
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/grass7/raster/r.geomorphon/main.c
===================================================================
--- grass-addons/grass7/raster/r.geomorphon/main.c (rev 0)
+++ grass-addons/grass7/raster/r.geomorphon/main.c 2013-02-05 10:55:21 UTC (rev 54918)
@@ -0,0 +1,465 @@
+/****************************************************************************
+*
+* MODULE: r.geomorphon
+* AUTHOR(S): Jarek Jasiewicz jarekj amu.edu.pl with collaboration of Tom Stepinski stepintz uc.edu based on idea of Tomasz Stepinski and Jaroslaw Jasiewicz
+*
+* PURPOSE: Calculate terrain forms using machine-vison technique called geomorphons
+* This module allow to calculate standard set of terrain forms
+* using look-up table proposed by authors, calculate patterns (binary and ternary)
+* for every pixel as well as several geomorphometrical parameters
+* This technology is currently capable of "experimental" stage.
+*
+* COPYRIGHT: (C) 2002,2012 by the GRASS Development Team
+* (C) Scientific idea of geomotrphon copyrighted to authors.
+*
+* This program is free software under the GNU General Public
+* License (>=v2). Read the file COPYING that comes with GRASS
+* for details.
+*
+ *****************************************************************************/
+
+#define MAIN
+#include "local_proto.h"
+typedef enum {i_dem, o_forms, o_ternary, o_positive, o_negative, o_intensity, 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 */
+ {"dem",YES,"Input dem","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_treshold,
+ *par_flat_distance,
+ *par_multi_prefix,
+ *par_multi_step,
+ *par_multi_start;
+ struct Flag *flag_units,
+ *flag_extended;
+
+ struct History history;
+
+ int i,j, n;
+ int meters=0, multires=0, extended=0; /* flags */
+ int row,cur_row,col,radius;
+ int pattern_size;
+ double max_resolution;
+ char prefix[20];
+
+ G_gisinit(argv[0]);
+
+{ /* interface parameters */
+ module = G_define_module();
+ module->description =
+ _("Calculate geomorphons (terrain forms)and associated geometry using machine vision approach");
+ G_add_keyword("Geomorphons");
+ G_add_keyword("Terrain patterns");
+ G_add_keyword("Machine vision geomorphometry");
+
+ opt_input = G_define_standard_option(G_OPT_R_INPUT);
+ 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_treshold = G_define_option();
+ par_flat_treshold->key = "flat";
+ par_flat_treshold->type = TYPE_DOUBLE;
+ par_flat_treshold->answer = "1";
+ par_flat_treshold->required = YES;
+ par_flat_treshold->description = _("Flatenss treshold (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 = _("Flatenss 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 serch 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"));
+
+ 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_treshold->answer);
+ if(flat_threshold<=0.)
+ G_fatal_error(_("Flatenss treshold 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(_("Flatenss 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_message("Search distance m: %f, cells: %d", search_distance, search_cells);
+ G_message("Skip distance m: %f, cells: %d", skip_distance, skip_cells);
+ G_message("Flat threshold distance m: %f, height: %f",flat_distance, flat_threshold_height);
+ G_message("%s version",(extended)?"extended":"basic");
+ if(multires) {
+ G_message("Multiresolution mode: search start at: m: %f, cells: %d", start_distance, start_cells);
+ G_message("Multiresolution mode: search step is: m: %f, number of steps %d", step_distance,num_of_steps);
+ G_message("Prefix for output: %s",prefix);
+ }
+}
+
+ /* generate global ternary codes */
+ for(i=0;i<6561;++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;
+ int formA, formB, formC;
+ 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_message("Done!");
+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[11].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);
+}
+
+}
+
+
Property changes on: grass-addons/grass7/raster/r.geomorphon/main.c
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/grass7/raster/r.geomorphon/memory.c
===================================================================
--- grass-addons/grass7/raster/r.geomorphon/memory.c (rev 0)
+++ grass-addons/grass7/raster/r.geomorphon/memory.c 2013-02-05 10:55:21 UTC (rev 54918)
@@ -0,0 +1,152 @@
+#include "local_proto.h"
+
+int open_map(MAPS* rast) {
+
+ int row, col;
+ int fd;
+ char* mapset;
+ struct Cell_head cellhd;
+ int bufsize;
+ void* tmp_buf;
+
+ mapset = (char*)G_find_raster2(rast->elevname, "");
+
+ if (mapset == NULL)
+ 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 rast=%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 */
+
+G_free(tmp_buf);
+return 0;
+}
+
+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;
+
+ 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;
+
+ 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;
+ }
+
+return 0;
+}
+
+int shift_buffers(int row)
+{
+int i;
+int col;
+void* tmp_buf;
+FCELL* tmp_elev_buf, *slope_tmp, *aspect_tmp;
+
+tmp_buf=Rast_allocate_buf(elevation.raster_type);
+tmp_elev_buf=elevation.elev[0];
+
+ for (i = 1; i < row_buffer_size+1; ++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);
+
+ for (col=0;col<ncols;++col)
+get_cell(col, elevation.elev[row_buffer_size], tmp_buf, elevation.raster_type);
+
+G_free(tmp_buf);
+return 0;
+}
+
+int free_map (FCELL **map, int n) {
+ int i;
+ for (i=0;i<n;++i)
+ G_free(map[i]);
+ G_free(map);
+ return 0;
+}
+int write_form_cat_colors (char* raster, CATCOLORS* ccolors) {
+struct Colors colors;
+struct Categories cats;
+int i;
+
+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_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_write_cats(raster, &cats);
+Rast_free_cats(&cats);
+return 0;
+}
+
+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}};
+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_write_colors(raster, G_mapset(), &colors);
+Rast_free_colors(&colors);
+/*
+Rast_init_cats("Forms", &cats);
+ for(i=0;i<8;++i)
+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;
+}
\ No newline at end of file
Property changes on: grass-addons/grass7/raster/r.geomorphon/memory.c
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/grass7/raster/r.geomorphon/multires.c
===================================================================
--- grass-addons/grass7/raster/r.geomorphon/multires.c (rev 0)
+++ grass-addons/grass7/raster/r.geomorphon/multires.c 2013-02-05 10:55:21 UTC (rev 54918)
@@ -0,0 +1,25 @@
+#include "local_proto.h"
+
+
+
+
+int pattern_matching(int* pattern)
+{
+int n, i;
+unsigned char binary=0, result=255, test=0;
+unsigned char source=32;
+int sign=-1;
+
+ for (i=0,n=1;i<8;i++,n*=2)
+ 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;
+ }
+return (result & source == source)? 1:0;
+
+}
\ No newline at end of file
Property changes on: grass-addons/grass7/raster/r.geomorphon/multires.c
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/grass7/raster/r.geomorphon/pattern.c
===================================================================
--- grass-addons/grass7/raster/r.geomorphon/pattern.c (rev 0)
+++ grass-addons/grass7/raster/r.geomorphon/pattern.c 2013-02-05 10:55:21 UTC (rev 54918)
@@ -0,0 +1,123 @@
+#include "local_proto.h"
+/*directions
+ * 3|2|1
+ * 4|0|8
+ * 5|6|7 */
+static int nextr[8] = {-1, -1, -1, 0, 1, 1, 1, 0 };
+static int nextc[8] = {1, 0, -1, -1, -1, 0, 1, 1 };
+
+int calc_pattern(PATTERN* pattern, int row, int cur_row, int col)
+{
+/* calculate parameters of geomorphons and store it in the struct pattern */
+int i,j,pattern_size=0;
+double zenith_angle, nadir_angle, angle;
+double nadir_threshold, zenith_threshold;
+double zenith_height, nadir_height, zenith_distance, nadir_distance;
+double cur_northing, cur_easting, target_northing, target_easting;
+double cur_distance;
+double center_height, height;
+
+/* use distance calculation */
+cur_northing=Rast_row_to_northing(row+0.5, &window);
+cur_easting=Rast_col_to_easting(col+0.5, &window);
+center_height=elevation.elev[cur_row][col];
+pattern->num_positives=0;
+pattern->num_negatives=0;
+pattern->positives=0;
+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;
+
+ 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);
+
+ 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);
+
+ 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_treshold) {
+ 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(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; //ZMIANA!
+ pattern->distance[i]=zenith_distance;
+ pattern->num_positives++;
+ }
+ if(fabs(nadir_angle)>fabs(zenith_angle)) {
+ pattern->pattern[i]=-1;
+ pattern->elevation[i]=nadir_height; //ZMIANA!
+ pattern->distance[i]=nadir_distance;
+ pattern->num_negatives++;
+ }
+ } else {
+ pattern->distance[i]=search_distance;
+ }
+
+} /* end for */
+return pattern_size;
+}
Property changes on: grass-addons/grass7/raster/r.geomorphon/pattern.c
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/grass7/raster/r.geomorphon/r.geomorphon.html
===================================================================
--- grass-addons/grass7/raster/r.geomorphon/r.geomorphon.html (rev 0)
+++ grass-addons/grass7/raster/r.geomorphon/r.geomorphon.html 2013-02-05 10:55:21 UTC (rev 54918)
@@ -0,0 +1,60 @@
+<h2>OPTIONS</h2>
+<DL>
+<DT><b>-m</b></DT>
+<DD>All distance parameters (search, skip, flat distances) are supplied as meters instead of cells (default). To avoid situation when supplied distances is smaller than one cell program first check if supplied distance is longer than one cell in both NS and WE directions. For LatLong projection only NS distance checked, because in latitude angular unit comprise always bigger or equal distance than longitude one. If distance is supplied in cells, For all projections is recalculated into meters according formula: <code>number_of_cells*resolution_along_NS_direction</code>. It is important if geomorphons are calculate for large areas in LatLong projecton.</DD>
+<DT><b>dem</b></DT>
+<DD>Digital elevation model. Data can be of any type and any projection. During calculation DEM is stored as floating point raster.</DD>
+<DT><b>search</b></DT>
+<DD>Determines length on the geodesic distances in all eight directions where line-of-sight is calculated. To speed up calculation is determines only these cells which centers falls into the distance</DD>
+<DT><b>skip</b></DT>
+<DD>Determines length on the geodesic distances at the beginning of calculation all eight directions where line-of-sight is yet calculated. To speed up calculation this distance is always recalculated into number of cell which are skipped at the begining of every line-of-sight and is equal in all direction. This parameter eliminates forms of very small extend, smaller than skip parameter.</DD>
+<DT><b>flat</b></DT>
+<DD>The difference (in degrees) between zenith and nadir line-of-sight which indicate flat direction. If higher threshold produce more flat maps. If resolution of the map is low (more than 1 km per cell) threshold should be very small (much smaller than 1 degree) because on such distance 1 degree of difference means several meters of high difference.</DD>
+<DT><b>dist</b></DT>
+<DD>>Flat distance. This is additional parameter defining the distance above which the threshold starts to decrease to avoid problems with pseudo-flat line-of-sights if real elevation difference appears on the distance where its value is higher DO POPRAWKI </DD>
+<DT><b>form</b></DT>
+<DD>Returns geomorphic map with 10 most popular terrestial forms. Legend for forms, its definition by the number of <em>+</em> and <em>-</em> and its idealized visualisation are presented at the image.
+<center>
+<h3>Forms represented by geomorphons:<h3>
+<img src=legend.png border=1><br>
+</center></DD>
+<DT><b>pattern</b></DT>
+<DD>returns code of one of 498 unique ternary patterns for every cell. The code is a decimal representation o 8-tuple minimalised patterns written in ternary system. Full list of patterns is available in source code directory as patterns.txt. This map can be used to create alternative form classification using supervised approach</DD>
+<DT><b>positive and negative</b></DT>
+<DD>returns codes binary patterns for zenith (positive) and nadir (negative) line of sights. The code is a decimal representation o 8-tuple minimalised patterns written in binary system. Full list of patterns is available in source code directory as patterns.txt</DD>
+</DL>
+<P><em>NOTE: parameters below are very experimental. The usefulness of these parameters are currently under investigation</em></P>
+<DL>
+<DT><b>intensity</b></DT>
+<DD>returns avarage difference between central cell of geomorphon and eight cells in visibility neighbourhood. This parameter shows local (as is visible) exposition/abasment of the form in the terrain</DD>
+<DT><b>range</b></DT>
+<DD>returns difference between minimum and maximum values of visibility neigbourhood.</DD>
+<DT><b>variance</b></DT>
+<DD>returns variance (difference between particular values and mean value) ofvisibility neigbourhood.</DD>
+<DT><b>extend</b></DT>
+<DD>returns area of the polygon created by the 8 points where line-of-sight cuts the terrain (see image in description section).</DD>
+<DT><b>azimuth</b></DT>
+<DD>returns orientation of the poligon constituting geomorphon. This orientation is currentlyb calculated as a orientation of least square fit line to the eight verticles of this polygon.</DD>
+<DT><b>elongation</b></DT>
+<DD>returns proportion between sides of the bounding box rectangle calculated for geomorphon rotated to fit lest square line.</DD>
+<DT><b>width</b></DT>
+<DD>returns lenght of the shorter side of the bounding box rectangle calculated for geomorphon rotated to fit lest square line.</DD>
+
+
+<h2>DESCRIPTION</h2>
+<center>
+<h3>Whar is geomorphon:<h3>
+<img src=geomorphon.png border=1><br>
+</center>
+<P>Geomorphon is a new concept of presentation and analysis of terrain forms. This concept utilises 8-tuple pattern of the visibility neighbourhood and breaks well known limitation of standard calculus approach where all terrain forms cannot be detected in a single window size. The pattern arises from a comparison of a focus pixel with its eight neighbors starting from the one located to the east and continuing counterclockwise producing ternary operator. For example, a tuple {+,-,-,-,0,+,+,+} describes one possible pattern of relative measures {higher, lower, lower, lower, equal, higher, higher, higher} for pixels surrounding the focus pixel. It is important to stress that the visibility neighbors are <b>not necessarily an immediate neighbors</b> of the focus pixel in the grid, but the pixels determined from <b>the line-of-sight</b> principle along the eight principal directions. This principle relates surface relief and horizontal distance by means of so-called zenith an
d nadir angles along the eight principal compass directions. The ternary operator converts the information contained in all the pairs of zenith and nadir angles into the ternary pattern (8-tuple). The result depends on the values of two parameters: search radius (L) and relief threshold (d). The search radius is the maximum allowable distance for calculation of zenith and nadir angles. The relief threshold is a minimum value of difference between LOSs angle (zenith and nadir) that is considered significantly different from the horizon. Two lines-of-sight are necessary due to zenith LOS only, does not detect positive forms correctly.
+<P>There are 38 = 6561 possible ternary patterns (8-tuplets). However by eliminating all patterns that are results of either rotation or reflection of other patterns wa set of 498 patterns remain referred as geomorphons. This is a comprehensive and exhaustive set of idealized landforms that are independent of the size, relief, and orientation of the actual landform.
+<P>Form recognition depends on two free parameters: <b>Search radius</b> and <b>flatness threshold</b>. Using larger values of L and is tantamount to terrain classification from a higher and wider perspective, whereas using smaller values of L and is tantamount to terrain classification from a local point of view. A character of the map depends on the value of L. Using small value of L results in the map that correctly identifies landforms if their size is smaller than L; landforms having larger sizes are broken down into components. Using larger values of L allows simultaneous identification of landforms on variety of sizes in expense of recognition smaller, second-order forms. There are two addational parameters: <b>skip radius</b> used to eliminate impact of small irregularities. On the contrary <b>flatness distance</b> eliminates the impact of very high distance (in meters) of search radius which may not detect elevation difference if this is at very far distance. Import
ant especially with low resolution DEMS.
+<H2>NOTES</H2>
+<P>From computational point of view there are no limitations of input DEM and free parameters used in calculation. However in practice there are some issues on DEM resolution and search radius. Low resolution DEM especially above 1 km per cell requires smaller than default flatness threshold. On the other hand, only forms with high local elevation difference will be detected correctly. It results form fact that on very high distance (of order of kilometers or higher) even relatively high elevation difference will be recognized as flat. For example at the distance of 8 km (8 cells with 1 km resolution DEM) an relative elevation difference of at least 136 m is required to be noticed as non-flat. Flatness distance threshold may be helpful to avoid this problem.
+<h2>SEE ALSO</h2>
+<h2>REFERENCES</h2>
+<P>
+Stepinski, T., Jasiewicz, J., 2011, Geomorphons - a new approach to classification of landform, in : Eds: Hengl, T., Evans, I.S., Wilson, J.P., and Gould, M.,, Proceedings to Geomorphometry 2011, Redlands, 109–112
+Jasiewicz, J., Stepinski, T., 2013, Geomorphons - a pattern recognition approach to classification and mapping of landforms, Geomorphology, vol. 182, 147–156
+<h2>AUTHOR</h2>
+Jarek Jasiewicz, Tomek Stepinski (merit contribution)
Property changes on: grass-addons/grass7/raster/r.geomorphon/r.geomorphon.html
___________________________________________________________________
Added: svn:executable
+ *
More information about the grass-commit
mailing list