r.geomorphon: moved to trunk from Addons

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@71632 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/Makefile

@@ -23,6 +23,7 @@ SUBDIRS = \
 	r.external.out \
 	r.fill.dir \
 	r.flow \
+	r.geomorphon \
 	r.grow.distance \
 	r.gwflow \
 	r.his \

raster/r.geomorphon/Makefile

@@ -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

raster/r.geomorphon/geom.c

@@ -0,0 +1,259 @@
+#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;
+
+    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;
+}

raster/r.geomorphon/local_proto.h

@@ -0,0 +1,179 @@
+#ifndef __LOCAL_PROTO_H__
+#define __LOCAL_PROTO_H__
+
+#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))
+
+#undef MIN
+#undef MAX
+#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 */
+GLOBAL 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);
+int write_contrast_colors(char *);
+
+/* 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);
+int radial2cartesian(PATTERN *);
+
+/* 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);
+
+#endif

raster/r.geomorphon/main.c

@@ -0,0 +1,553 @@
+
+/****************************************************************************
+*
+* 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 */
+	{"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_treshold,
+	*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 row, cur_row, col;
+    int pattern_size;
+    double max_resolution;
+    char prefix[20];
+
+    G_gisinit(argv[0]);
+
+    {				/* 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_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_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);
+
+    /* 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 */
+
+    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);
+    }
+
+}

raster/r.geomorphon/memory.c

@@ -0,0 +1,161 @@
+#include "local_proto.h"
+
+int open_map(MAPS * rast)
+{
+
+    int row, col;
+    char *mapset;
+    struct Cell_head cellhd;
+    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 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 */
+
+    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;
+
+    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;
+}

raster/r.geomorphon/multires.c

@@ -0,0 +1,22 @@
+#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;
+
+}

raster/r.geomorphon/pattern.c

@@ -0,0 +1,145 @@
+#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;
+}

raster/r.geomorphon/testsuite/test_r_geom.py

@@ -0,0 +1,72 @@
+"""
+Name:       r.geomorphon tests
+Purpose:    Tests r.geomorphon input parsing.
+            Uses NC Basic data set.
+
+Author:     Luca Delucchi, Markus Neteler
+Copyright:  (C) 2017 by Luca Delucchi, Markus Neteler and the GRASS Development Team
+Licence:    This program is free software under the GNU General Public
+            License (>=v2). Read the file COPYING that comes with GRASS
+            for details.
+"""
+
+from grass.gunittest.case import TestCase
+from grass.gunittest.main import test
+from grass.script.core import read_command
+
+synth_out = \
+"""1	flat
+3	ridge
+4	shoulder
+6	slope
+8	footslope
+9	valley
+"""
+
+ele_out = \
+"""1	flat
+2	summit
+3	ridge
+4	shoulder
+5	spur
+6	slope
+7	hollow
+8	footslope
+9	valley
+10	depression
+"""
+
+class TestClipling(TestCase):
+    inele = 'elevation'
+    insint = 'synthetic_dem'
+    outele = 'ele_geomorph'
+    outsint = 'synth_geomorph'
+
+    @classmethod
+    def setUpClass(cls):
+        """Ensures expected computational region and generated data"""
+        cls.use_temp_region()
+        cls.runModule('g.region', raster=cls.inele)
+        cls.runModule('r.mapcalc', expression="{ou} = sin(x() / 5.0) + (sin(x() / 5.0) * 100.0 + 200)".format(ou=cls.insint))
+
+    @classmethod
+    def tearDownClass(cls):
+        """Remove the temporary region and generated data"""
+        cls.runModule('g.remove', flags='f', type='raster',
+                      name=(cls.insint, cls.outele, cls.outsint))
+        cls.del_temp_region()
+    
+    def test_ele(self):
+        self.runModule('r.geomorphon', elevation=self.inele, forms=self.outele,
+                      search=10)
+        category = read_command('r.category', map=self.outele)
+        self.assertEqual(first=ele_out, second=category)
+
+    def test_sint(self):
+        self.runModule('r.geomorphon', elevation=self.insint,
+                       forms=self.outsint, search=10)
+        category = read_command('r.category', map=self.outsint)
+        self.assertEqual(first=synth_out, second=category)
+    
+if __name__ == '__main__':
+    test()

raster/r.param.scale/r.param.scale.html

@@ -104,11 +104,9 @@ tables for all output files (presently only on features).
 </ul>
 
 <h2>SEE ALSO</h2>
-<!-- not ported to GRASS 6 due to non-GPLness of numerical recipes.
-<i><a href="d.param.scale.html">d.param.scale</a></i>
--->
 
 <em>
+  <a href="r.geomorphon.html">r.geomorphon</a>,
   <a href="r.slope.aspect.html">r.slope.aspect</a>
 </em>