[GRASS-SVN] r52150 - in grass-addons/grass7/raster: . r.canny
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Jun 20 00:58:15 PDT 2012
Author: wenzeslaus
Date: 2012-06-20 00:58:15 -0700 (Wed, 20 Jun 2012)
New Revision: 52150
Added:
grass-addons/grass7/raster/r.canny/
grass-addons/grass7/raster/r.canny/Makefile
grass-addons/grass7/raster/r.canny/canny.c
grass-addons/grass7/raster/r.canny/canny.h
grass-addons/grass7/raster/r.canny/gauss.c
grass-addons/grass7/raster/r.canny/gauss.h
grass-addons/grass7/raster/r.canny/main.c
grass-addons/grass7/raster/r.canny/r.canny.html
Log:
canny edge detector: initial version
Added: grass-addons/grass7/raster/r.canny/Makefile
===================================================================
--- grass-addons/grass7/raster/r.canny/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.canny/Makefile 2012-06-20 07:58:15 UTC (rev 52150)
@@ -0,0 +1,12 @@
+# fix this relative to include/
+# or use absolute path to the GRASS source code
+MODULE_TOPDIR = ../..
+
+PGM = r.canny
+
+LIBES = $(GISLIB) $(RASTERLIB) $(GMATHLIB)
+DEPENDENCIES = $(GISDEP) $(RASTERDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass-addons/grass7/raster/r.canny/canny.c
===================================================================
--- grass-addons/grass7/raster/r.canny/canny.c (rev 0)
+++ grass-addons/grass7/raster/r.canny/canny.c 2012-06-20 07:58:15 UTC (rev 52150)
@@ -0,0 +1,344 @@
+/*
+ Canny Edge Detector Implementation
+ http://www.tomgibara.com/computer-vision/canny-edge-detector
+ Public Domain.
+
+ NOTE from Java source code: The elements of the method below (specifically the technique for
+ non-maximal suppression and the technique for gradient computation)
+ are derived from an implementation posted in the following forum (with the
+ clear intent of others using the code):
+ http://forum.java.sun.com/thread.jspa?threadID=546211&start=45&tstart=0
+ My code effectively mimics the algorithm exhibited above.
+ Since I don't know the providence of the code that was posted it is a
+ possibility (though I think a very remote one) that this code violates
+ someone's intellectual property rights. If this concerns you feel free to
+ contact me for an alternative, though less efficient, implementation.
+
+ NOTE - Citation of forum Terms of Use accoding to http://archive.org/:
+ Java Technology Forums
+ Sun.com Terms of Use
+ http://www.sun.com/termsofuse.html
+
+ 4. CONTENT SUBMITTED TO SUN
+
+ 4.1 Sun does not claim ownership of the Content You place on the Website and
+ shall have no obligation of any kind with respect to such Content.
+ Unless otherwise stated herein, or in Sun's Privacy Policy,
+ any Content You provide in connection with this Website shall be deemed
+ to be provided on a nonconfidential basis. Sun shall be free to use
+ or disseminate such Content on an unrestricted basis for any purpose,
+ and You grant Sun and all other users of the Website an irrevocable,
+ worldwide, royalty-free, nonexclusive license to use, reproduce, modify,
+ distribute, transmit, display, perform, adapt,
+ resell and publish such Content (including in digital form).
+ You represent and warrant that you have proper authorization
+ for the worldwide transfer and processing among Sun, its affiliates,
+ and third-party providers of any information
+ that You may provide on the Website.
+ */
+
+#include "canny.h"
+
+#include <math.h>
+
+void computeXGradients(DCELL * diffKernel, DCELL * yConv, DCELL * xGradient,
+ int rows, int cols, int kernelWidth)
+{
+ int initX = kernelWidth - 1;
+
+ int maxX = cols - (kernelWidth - 1);
+
+ int initY = cols * (kernelWidth - 1);
+
+ int maxY = cols * (rows - (kernelWidth - 1));
+
+ int x;
+
+ int y;
+
+ int i;
+
+ for (x = initX; x < maxX; x++) {
+ for (y = initY; y < maxY; y += cols) {
+ float sum = 0.;
+
+ int index = x + y;
+
+ for (i = 1; i < kernelWidth; i++) {
+ sum += diffKernel[i] * (yConv[index - i] - yConv[index + i]);
+ }
+ xGradient[index] = sum;
+ }
+ }
+}
+
+
+void computeYGradients(DCELL * diffKernel, DCELL * xConv, DCELL * yGradient,
+ int rows, int cols, int kernelWidth)
+{
+ int initY = cols * (kernelWidth - 1);
+
+ int maxY = cols * (rows - (kernelWidth - 1));
+
+ int x;
+
+ int y;
+
+ int i;
+
+ for (x = kernelWidth; x < cols - kernelWidth; x++) {
+ for (y = initY; y < maxY; y += cols) {
+ float sum = 0.0;
+
+ int index = x + y;
+
+ int yOffset = cols;
+
+ for (i = 1; i < kernelWidth; i++) {
+ sum +=
+ diffKernel[i] * (xConv[index - yOffset] -
+ xConv[index + yOffset]);
+ yOffset += cols;
+ }
+ yGradient[index] = sum;
+ }
+ }
+}
+
+static float custom_hypot(float x, float y)
+{
+ float t;
+
+ x = fabs(x);
+ y = fabs(y);
+
+ if (x < y) {
+ t = x;
+ x = y;
+ y = t;
+ }
+
+ if (x == 0.0)
+ return 0.0;
+ if (y == 0.0)
+ return x;
+ return x * sqrt(1 + (y / x) * (y / x));
+}
+
+/*
+ * An explanation of what's happening here, for those who want
+ * to understand the source: This performs the "non-maximal
+ * supression" phase of the Canny edge detection in which we
+ * need to compare the gradient magnitude to that in the
+ * direction of the gradient; only if the value is a local
+ * maximum do we consider the point as an edge candidate.
+ *
+ * We need to break the comparison into a number of different
+ * cases depending on the gradient direction so that the
+ * appropriate values can be used. To avoid computing the
+ * gradient direction, we use two simple comparisons: first we
+ * check that the partial derivatives have the same sign (1)
+ * and then we check which is larger (2). As a consequence, we
+ * have reduced the problem to one of four identical cases that
+ * each test the central gradient magnitude against the values at
+ * two points with 'identical support'; what this means is that
+ * the geometry required to accurately interpolate the magnitude
+ * of gradient function at those points has an identical
+ * geometry (upto right-angled-rotation/reflection).
+ *
+ * When comparing the central gradient to the two interpolated
+ * values, we avoid performing any divisions by multiplying both
+ * sides of each inequality by the greater of the two partial
+ * derivatives. The common comparand is stored in a temporary
+ * variable (3) and reused in the mirror case (4).
+ *
+ */
+static int isLocalMax(float xGrad, float yGrad, float gradMag,
+ float neMag, float seMag, float swMag, float nwMag,
+ float nMag, float eMag, float sMag, float wMag)
+{
+ float tmp, tmp1, tmp2;
+
+ if (xGrad * yGrad <= 0.0f) {
+ if (fabs(xGrad) >= fabs(yGrad)) {
+ tmp = fabs(xGrad * gradMag);
+ tmp1 = fabs(yGrad * neMag - (xGrad + yGrad) * eMag) /*(3) */ ;
+ tmp2 = fabs(yGrad * swMag - (xGrad + yGrad) * wMag) /*(4) */ ;
+ }
+ else {
+ tmp = fabs(yGrad * gradMag);
+ tmp1 = fabs(xGrad * neMag - (yGrad + xGrad) * nMag) /*(3) */ ;
+ tmp2 = fabs(xGrad * swMag - (yGrad + xGrad) * sMag) /*(4) */ ;
+ }
+ }
+ else {
+ if (fabs(xGrad) >= fabs(yGrad) /*(2) */ ) {
+ tmp = fabs(xGrad * gradMag);
+ tmp1 = fabs(yGrad * seMag + (xGrad - yGrad) * eMag) /*(3) */ ;
+ tmp2 = fabs(yGrad * nwMag + (xGrad - yGrad) * wMag) /*(4) */ ;
+ }
+ else {
+ tmp = fabs(yGrad * gradMag);
+ tmp1 = fabs(xGrad * seMag + (yGrad - xGrad) * sMag) /*(3) */ ;
+ tmp2 = fabs(xGrad * nwMag + (yGrad - xGrad) * nMag) /*(4) */ ;
+ }
+ }
+ if (tmp >= tmp1 && tmp > tmp2) {
+ return 1;
+ }
+ return 0;
+}
+
+void nonmaxSuppresion(DCELL * xGradient, DCELL * yGradient, DCELL * magnitude,
+ DCELL * angle,
+ int rows, int cols, int kernelWidth,
+ float magnitudeScale, float magnitudeLimit)
+{
+ int initX = kernelWidth;
+
+ int maxX = cols - kernelWidth;
+
+ int initY = cols * kernelWidth;
+
+ int maxY = cols * (rows - kernelWidth);
+
+ int x;
+
+ int y;
+
+ int MAGNITUDE_MAX = magnitudeScale * magnitudeLimit;
+
+ for (x = initX; x < maxX; x++) {
+ for (y = initY; y < maxY; y += cols) {
+ int index = x + y;
+
+ int indexN = index - cols;
+
+ int indexS = index + cols;
+
+ int indexW = index - 1;
+
+ int indexE = index + 1;
+
+ int indexNW = indexN - 1;
+
+ int indexNE = indexN + 1;
+
+ int indexSW = indexS - 1;
+
+ int indexSE = indexS + 1;
+
+ float xGrad = xGradient[index];
+
+ float yGrad = yGradient[index];
+
+ float gradMag = custom_hypot(xGrad, yGrad);
+
+ /* perform non-maximal supression */
+ float nMag = custom_hypot(xGradient[indexN], yGradient[indexN]);
+
+ float sMag = custom_hypot(xGradient[indexS], yGradient[indexS]);
+
+ float wMag = custom_hypot(xGradient[indexW], yGradient[indexW]);
+
+ float eMag = custom_hypot(xGradient[indexE], yGradient[indexE]);
+
+ float neMag =
+ custom_hypot(xGradient[indexNE], yGradient[indexNE]);
+
+ float seMag =
+ custom_hypot(xGradient[indexSE], yGradient[indexSE]);
+
+ float swMag =
+ custom_hypot(xGradient[indexSW], yGradient[indexSW]);
+
+ float nwMag =
+ custom_hypot(xGradient[indexNW], yGradient[indexNW]);
+
+ if (isLocalMax(xGrad, yGrad, gradMag, neMag, seMag, swMag, nwMag,
+ nMag, eMag, sMag, wMag)) {
+ magnitude[index] =
+ gradMag >=
+ magnitudeLimit ? MAGNITUDE_MAX : (int)(magnitudeScale *
+ gradMag);
+ /*
+ NOTE: The orientation of the edge is not employed by this
+ implementation. It is a simple matter to compute it at
+ this point as: Math.atan2(yGrad, xGrad);
+ */
+ if (angle != NULL)
+ {
+ angle[index] = atan2(yGrad, xGrad) * 180 / M_PI;;
+ }
+ }
+ else {
+ magnitude[index] = 0;
+ }
+ }
+ }
+}
+
+static void follow(DCELL * edges, DCELL * magnitude, int x1, int y1, int i1,
+ int threshold, int rows, int cols)
+{
+ int x0 = x1 == 0 ? x1 : x1 - 1;
+
+ int x2 = x1 == cols - 1 ? x1 : x1 + 1;
+
+ int y0 = y1 == 0 ? y1 : y1 - 1;
+
+ int y2 = y1 == rows - 1 ? y1 : y1 + 1;
+
+ int x;
+
+ int y;
+
+ edges[i1] = magnitude[i1];
+ for (x = x0; x <= x2; x++) {
+ for (y = y0; y <= y2; y++) {
+ int i2 = x + y * cols;
+
+ if ((y != y1 || x != x1) && edges[i2] == 0
+ && magnitude[i2] >= threshold) {
+ follow(edges, magnitude, x, y, i2, threshold, rows, cols);
+ return;
+ }
+ }
+ }
+}
+
+/* edges.fill(0) */
+void performHysteresis(DCELL * edges, DCELL * magnitude, int low, int high,
+ int rows, int cols)
+{
+ /*
+ NOTE: this implementation reuses the data array to store both
+ luminance data from the image, and edge intensity from the processing.
+ This is done for memory efficiency, other implementations may wish
+ to separate these functions.
+ */
+
+ int x;
+
+ int y;
+
+ int offset = 0;
+
+ for (y = 0; y < rows; y++) {
+ for (x = 0; x < cols; x++) {
+ if (edges[offset] == 0 && magnitude[offset] >= high) {
+ follow(edges, magnitude, x, y, offset, low, rows, cols);
+ }
+ offset++;
+ }
+ }
+}
+
+void thresholdEdges(DCELL * edges, int rows, int cols)
+{
+ int i;
+
+ for (i = 0; i < rows * cols; i++) {
+ edges[i] = edges[i] > 0 ? 1 : 0;
+ }
+}
Added: grass-addons/grass7/raster/r.canny/canny.h
===================================================================
--- grass-addons/grass7/raster/r.canny/canny.h (rev 0)
+++ grass-addons/grass7/raster/r.canny/canny.h 2012-06-20 07:58:15 UTC (rev 52150)
@@ -0,0 +1,21 @@
+#ifndef CANNY_H
+#define CANNY_H
+
+#include <grass/gis.h>
+
+void computeXGradients(DCELL * diffKernel, DCELL * yConv, DCELL * xGradient,
+ int rows, int cols, int kernelWidth);
+
+void computeYGradients(DCELL * diffKernel, DCELL * xConv, DCELL * yGradient,
+ int rows, int cols, int kernelWidth);
+
+void nonmaxSuppresion(DCELL * xGradient, DCELL * yGradient, DCELL * magnitude,
+ DCELL *angle, int rows, int cols,
+ int kernelWidth, float magnitudeScale, float magnitudeLimit);
+
+void performHysteresis(DCELL * edges, DCELL * magnitude,
+ int low, int high, int rows, int cols);
+
+void thresholdEdges(DCELL * edges, int rows, int cols);
+
+#endif /* CANNY_H */
Added: grass-addons/grass7/raster/r.canny/gauss.c
===================================================================
--- grass-addons/grass7/raster/r.canny/gauss.c (rev 0)
+++ grass-addons/grass7/raster/r.canny/gauss.c 2012-06-20 07:58:15 UTC (rev 52150)
@@ -0,0 +1,76 @@
+#include "gauss.h"
+
+#include <math.h>
+
+#include <grass/gis.h>
+
+int getKernelWidth(const float sigma, float gaussianCutOff)
+{
+ return ceil(sqrt(-2 * sigma * sigma * log(gaussianCutOff)));
+}
+
+static float gaussian(float x, float sigma)
+{
+ return exp(-(x * x) / (2.0 * sigma * sigma));
+}
+
+void gaussKernel(DCELL * gaussKernel, DCELL * diffKernel,
+ int kernelWidth, float kernelRadius)
+{
+ int kwidth;
+
+ for (kwidth = 0; kwidth < kernelWidth; kwidth++) {
+ float g1 = gaussian(kwidth, kernelRadius);
+
+ float g2 = gaussian(kwidth - 0.5, kernelRadius);
+
+ float g3 = gaussian(kwidth + 0.5, kernelRadius);
+
+ gaussKernel[kwidth] =
+ (g1 + g2 +
+ g3) / 3. / (2.0 * (float)M_PI * kernelRadius * kernelRadius);
+ diffKernel[kwidth] = g3 - g2;
+ }
+}
+
+void gaussConvolution(DCELL * image, DCELL * kernel, DCELL * xConv,
+ DCELL * yConv, int rows, int cols, int kernelWidth)
+{
+ int x, y;
+
+ int initX = kernelWidth - 1;
+
+ int maxX = cols - (kernelWidth - 1);
+
+ int initY = cols * (kernelWidth - 1);
+
+ int maxY = cols * (rows - (kernelWidth - 1));
+
+ //perform convolution in x and y directions
+ for (x = initX; x < maxX; x++) {
+ for (y = initY; y < maxY; y += cols) {
+ int index = x + y;
+
+ float sumX = image[index] * kernel[0];
+
+ float sumY = sumX;
+
+ int xOffset = 1;
+
+ int yOffset = cols;
+
+ for (; xOffset < kernelWidth;) {
+ sumY +=
+ kernel[xOffset] * (image[index - yOffset] +
+ image[index + yOffset]);
+ sumX +=
+ kernel[xOffset] * (image[index - xOffset] +
+ image[index + xOffset]);
+ yOffset += cols;
+ xOffset++;
+ }
+ yConv[index] = sumY;
+ xConv[index] = sumX;
+ }
+ }
+}
Added: grass-addons/grass7/raster/r.canny/gauss.h
===================================================================
--- grass-addons/grass7/raster/r.canny/gauss.h (rev 0)
+++ grass-addons/grass7/raster/r.canny/gauss.h 2012-06-20 07:58:15 UTC (rev 52150)
@@ -0,0 +1,13 @@
+#ifndef GAUSS_H
+#define GAUSS_H
+
+#include <grass/gis.h>
+
+int getKernelWidth(const float sigma, float gaussianCutOff);
+
+void gaussKernel(DCELL * gaussKernel, DCELL * diffKernel,
+ int kernelWidth, float kernelRadius);
+void gaussConvolution(DCELL * image, DCELL * kernel, DCELL * xConv,
+ DCELL * yConv, int rows, int cols, int kernelWidth);
+
+#endif /* GAUSS_H */
Added: grass-addons/grass7/raster/r.canny/main.c
===================================================================
--- grass-addons/grass7/raster/r.canny/main.c (rev 0)
+++ grass-addons/grass7/raster/r.canny/main.c 2012-06-20 07:58:15 UTC (rev 52150)
@@ -0,0 +1,363 @@
+
+/****************************************************************************
+ *
+ * MODULE: r.canny
+ * AUTHOR(S): Anna Kratochvilova - kratochanna gmail.com
+ * Vaclav Petras - wenzeslaus gmail.com
+ *
+ * PURPOSE: Edge detection usig Canny algorithm.
+ *
+ * COPYRIGHT: (C) 2012 by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+
+
+#include <grass/imagery.h>
+#include <grass/gmath.h>
+
+#include <math.h>
+
+#include "canny.h"
+#include "gauss.h"
+
+/** Loads map into memory.
+
+ \param[out] mat map in a matrix (row order), field have to be allocated
+ */
+static void readMap(const char *name, const char *mapset, int nrows,
+ int ncols, DCELL * mat)
+{
+
+ int r, c;
+
+ int map_fd;
+
+ int check_reading;
+
+ DCELL *row_buffer;
+
+ DCELL cell_value;
+
+ row_buffer = Rast_allocate_d_input_buf();
+ check_reading = 0;
+
+
+ /* load map */
+ map_fd = Rast_open_old(name, mapset);
+ if (map_fd < 0) {
+ G_fatal_error(_("Error opening first raster map <%s>"), name);
+ }
+
+ G_debug(1, "fd %d %s %s", map_fd, name, mapset);
+
+ // if ((first_map_R_type =
+ // Rast_map_type(templName, mapset)) < 0)
+ // G_fatal_error(_("Error getting first raster map type"));
+
+ for (r = 0; r < nrows; r++) {
+ Rast_get_row(map_fd, row_buffer, r, DCELL_TYPE);
+
+ for (c = 0; c < ncols; c++) {
+ cell_value = row_buffer[c];
+ if (!Rast_is_d_null_value(&cell_value))
+ mat[(ncols * (r)) + c] = cell_value;
+ else
+ mat[(ncols * (r)) + c] = 0.0;
+
+ if (mat[(ncols * (r)) + c])
+ check_reading = 1;
+ }
+ }
+ G_free(row_buffer);
+
+ if (!check_reading)
+ G_fatal_error(_("Nothing read from map %d"), check_reading);
+
+ Rast_close(map_fd);
+}
+
+/** Writes map from memory into the file.
+
+ \param[in] map map in a matrix (row order)
+ */
+static void writeMap(const char *name, int nrows, int ncols, DCELL * map)
+{
+ unsigned char *outrast; /* output buffer */
+
+ int outfd;
+
+ outfd = Rast_open_new(name, DCELL_TYPE); // FIXME: using both open old and open new
+ int r, c;
+
+ outrast = Rast_allocate_buf(DCELL_TYPE);
+ for (r = 0; r < nrows; r++) {
+ for (c = 0; c < ncols; c++) {
+ int index = r * ncols + c;
+
+ DCELL value = map[index];
+
+ ((DCELL *) outrast)[c] = value;
+ }
+ Rast_put_row(outfd, outrast, DCELL_TYPE);
+ }
+ G_free(outrast);
+
+ Rast_close(outfd);
+}
+
+
+/**
+
+ \todo Floats are used instead of doubles.
+ \todo Be able to work with FCELL (and CELL?).
+ */
+int main(int argc, char *argv[])
+{
+ struct Cell_head cell_head; /* it stores region information,
+ and header information of rasters */
+ char *name; /* input raster name */
+
+ char *mapset; /* mapset name */
+
+ int kernelWidth;
+
+ float kernelRadius;
+
+ char *result; /* output raster name */
+ char *anglesMapName;
+
+ static const float GAUSSIAN_CUT_OFF = 0.005;
+
+ static const float MAGNITUDE_SCALE = 100.;
+
+ static const float MAGNITUDE_LIMIT = 1000.;
+
+ float lowThreshold, highThreshold, low, high;
+
+ int nrows, ncols, dim_2;
+
+// struct History history; /* holds meta-data (title, comments,..) */
+ struct GModule *module; /* GRASS module for parsing arguments */
+
+ /* options */
+ struct Option *input, *output, *angleOutput,
+ *lowThresholdOption, *highThresholdOption, *sigmaOption;
+
+ int r;
+
+ /* initialize GIS environment */
+ G_gisinit(argv[0]); /* reads grass env, stores program name to G_program_name() */
+
+ /* initialize module */
+ module = G_define_module();
+ G_add_keyword(_("raster"));
+ G_add_keyword(_("canny"));
+ G_add_keyword(_("edge detection"));
+ module->description =
+ _("Canny edge detector. Region shall be set to input map."
+ "Can work only on small images since map is loaded into memory.");
+
+ /* Define the different options as defined in gis.h */
+ input = G_define_standard_option(G_OPT_R_INPUT);
+
+ output = G_define_standard_option(G_OPT_R_OUTPUT);
+
+ angleOutput = G_define_standard_option(G_OPT_R_OUTPUT);
+ angleOutput->key = "angles_map";
+ angleOutput->required = NO;
+ angleOutput->description = _("Map with angles");
+ angleOutput->guisection = _("Outputs");
+
+ lowThresholdOption = G_define_option();
+ lowThresholdOption->key = "low_threshold";
+ lowThresholdOption->type = TYPE_DOUBLE;
+ lowThresholdOption->required = NO;
+ lowThresholdOption->multiple = NO;
+ lowThresholdOption->description = _("Low treshold for edges in Canny");
+ lowThresholdOption->answer = "3";
+ // lowThresholdOption->options = "1-10";
+
+ highThresholdOption = G_define_option();
+ highThresholdOption->key = "high_threshold";
+ highThresholdOption->type = TYPE_DOUBLE;
+ highThresholdOption->required = NO;
+ highThresholdOption->multiple = NO;
+ highThresholdOption->description = _("High treshold for edges in Canny");
+ highThresholdOption->answer = "10";
+ // lowThresholdOption->options = "1-10";
+
+ sigmaOption = G_define_option();
+ sigmaOption->key = "sigma";
+ sigmaOption->type = TYPE_DOUBLE;
+ sigmaOption->required = NO;
+ sigmaOption->multiple = NO;
+ sigmaOption->description = _("Kernel radius");
+ sigmaOption->answer = "2";
+
+ /* options and flags parser */
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ lowThreshold = atof(lowThresholdOption->answer);
+ highThreshold = atof(highThresholdOption->answer);
+
+ low = (int)((lowThreshold * MAGNITUDE_SCALE) + 0.5);
+ high = (int)((highThreshold * MAGNITUDE_SCALE) + 0.5);
+
+
+ kernelRadius = atoi(sigmaOption->answer);
+
+ result = output->answer;
+
+ /* stores options and flags to variables */
+ name = input->answer;
+
+ anglesMapName = angleOutput->answer;
+
+ /* returns NULL if the map was not found in any mapset,
+ * mapset name otherwise */
+ mapset = (char *)G_find_raster2(name, "");
+ if (mapset == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), name);
+
+ /* determine the inputmap type (CELL/FCELL/DCELL) */
+ //data_type = Rast_map_type(name, mapset);
+
+ /* Rast_open_old - returns file destriptor (>0) */
+ // infd = Rast_open_old(name, mapset);
+
+
+
+ // struct Cell_head templCellhd;
+
+ // Rast_get_cellhd(name, mapset, &cellhd);
+ // Rast_get_cellhd(first_map_R_name, first_map_R_mapset, &cellhd_zoom1);
+
+ /* controlling, if we can open input raster */
+ Rast_get_cellhd(name, mapset, &cell_head);
+
+ G_debug(3, "number of rows %d", cell_head.rows);
+
+ nrows = cell_head.rows;
+
+ ncols = cell_head.cols;
+
+ dim_2 = nrows * ncols;
+
+ DCELL *mat1;
+
+ /* Memory allocation for map_1: */
+ mat1 = (DCELL *) G_calloc((dim_2), sizeof(DCELL));
+
+
+ // FIXME: it is necessary?
+ for (r = 0; r < dim_2; r++) {
+ mat1[r] = 0.0;
+ }
+
+ // if ((first_map_R_type =
+ // Rast_map_type(templName, mapset)) < 0)
+ // G_fatal_error(_("Error getting first raster map type"));
+
+
+ readMap(name, mapset, nrows, ncols, mat1);
+
+ /* **** */
+
+ kernelWidth = getKernelWidth(kernelRadius, GAUSSIAN_CUT_OFF);
+
+ DCELL *kernel;
+
+ DCELL *diffKernel;
+
+ kernel = (DCELL *) G_calloc((kernelWidth), sizeof(DCELL));
+ diffKernel = (DCELL *) G_calloc((kernelWidth), sizeof(DCELL));
+ gaussKernel(kernel, diffKernel, kernelWidth, kernelRadius);
+
+
+ DCELL *yConv = (DCELL *) G_calloc((dim_2), sizeof(DCELL));
+
+ DCELL *xConv = (DCELL *) G_calloc((dim_2), sizeof(DCELL));
+
+ for (r = 0; r < dim_2; r++) {
+ yConv[r] = xConv[r] = 0;
+ }
+ gaussConvolution(mat1, kernel, xConv, yConv, nrows, ncols, kernelWidth);
+
+
+ DCELL *yGradient = (DCELL *) G_calloc((dim_2), sizeof(DCELL));
+
+ DCELL *xGradient = (DCELL *) G_calloc((dim_2), sizeof(DCELL));
+
+ for (r = 0; r < dim_2; r++) {
+ yGradient[r] = xGradient[r] = 0;
+ }
+
+
+ computeXGradients(diffKernel, yConv, xGradient, nrows, ncols,
+ kernelWidth);
+ computeYGradients(diffKernel, xConv, yGradient, nrows, ncols,
+ kernelWidth);
+
+
+ DCELL *magnitude = (DCELL *) G_calloc((dim_2), sizeof(DCELL));
+
+ DCELL *angle = NULL;
+ if (anglesMapName != NULL)
+ {
+ angle = (DCELL *) G_calloc((dim_2), sizeof(DCELL));
+
+ for (r = 0; r < dim_2; r++) {
+ angle[r] = 0;
+ }
+ }
+
+ nonmaxSuppresion(xGradient, yGradient, magnitude, angle,
+ nrows, ncols, kernelWidth,
+ MAGNITUDE_SCALE, MAGNITUDE_LIMIT);
+
+
+ DCELL *edges = (DCELL *) G_calloc((dim_2), sizeof(DCELL));
+
+ for (r = 0; r < dim_2; r++) {
+ edges[r] = 0;
+ }
+
+ performHysteresis(edges, magnitude, low, high, nrows, ncols);
+
+ thresholdEdges(edges, nrows, ncols);
+
+ writeMap(result, nrows, ncols, edges);
+
+ if (angle != NULL)
+ {
+ writeMap(anglesMapName, nrows, ncols, angle);
+ }
+
+ /* **** */
+
+ /* memory cleanup */
+ G_free(kernel);
+ G_free(diffKernel);
+ G_free(name);
+
+
+// /* add command line incantation to history file */
+// Rast_short_history(result, "raster", &history);
+// Rast_command_history(&history);
+// Rast_write_history(result, &history);
+
+
+ exit(EXIT_SUCCESS);
+}
Added: grass-addons/grass7/raster/r.canny/r.canny.html
===================================================================
--- grass-addons/grass7/raster/r.canny/r.canny.html (rev 0)
+++ grass-addons/grass7/raster/r.canny/r.canny.html 2012-06-20 07:58:15 UTC (rev 52150)
@@ -0,0 +1,18 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.example</em> does practically do nothing, except
+for illustrating GRASS raster programming. It copies
+over an existing raster map to a new raster map.
+See the source code for details.
+
+<h2>SEE ALSO</h2>
+
+<em><a href="r.stats.html">r.stats</A></em><br>
+<em><a href="http://grass.itc.it/devel/index.php#prog">GRASS Programmer's Manual</A></em>
+
+
+<h2>AUTHOR</h2>
+
+GRASS Development Team
+
+<p><i>Last changed: $Date: 2011-09-29 21:18:47 +0200 (Thu, 29 Sep 2011) $</i>
More information about the grass-commit
mailing list