[GRASS-SVN] r52153 - in grass-addons/grass7/raster: . r.ht
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Jun 20 01:00:10 PDT 2012
Author: wenzeslaus
Date: 2012-06-20 01:00:10 -0700 (Wed, 20 Jun 2012)
New Revision: 52153
Added:
grass-addons/grass7/raster/r.ht/
grass-addons/grass7/raster/r.ht/Makefile
grass-addons/grass7/raster/r.ht/extract_line.cpp
grass-addons/grass7/raster/r.ht/extract_line.h
grass-addons/grass7/raster/r.ht/hough.cpp
grass-addons/grass7/raster/r.ht/hough.h
grass-addons/grass7/raster/r.ht/houghtransform.cpp
grass-addons/grass7/raster/r.ht/houghtransform.h
grass-addons/grass7/raster/r.ht/main.cpp
grass-addons/grass7/raster/r.ht/matrix.cpp
grass-addons/grass7/raster/r.ht/matrix.h
grass-addons/grass7/raster/r.ht/r.ht.html
Log:
extracting line segments using hough transform: initial version
Added: grass-addons/grass7/raster/r.ht/Makefile
===================================================================
--- grass-addons/grass7/raster/r.ht/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.ht/Makefile 2012-06-20 08:00:10 UTC (rev 52153)
@@ -0,0 +1,20 @@
+# fix this relative to include/
+# or use absolute path to the GRASS source code
+MODULE_TOPDIR = ../..
+
+PGM = r.ht
+
+LIBES = $(VECTORLIB) $(GISLIB) $(RASTERLIB) $(GMATHLIB)
+DEPENDENCIES = $(VECTORDEP) $(GISDEP) $(RASTERDEP)
+
+EXTRA_INC = $(VECT_INC)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+EXTRA_CFLAGS = $(VECT_CFLAGS) -Wno-sign-compare -Wall -Wextra -O0
+
+LINK = $(CXX)
+
+#ifneq ($(strip $(CXX)),)
+default: cmd
+#endif
Added: grass-addons/grass7/raster/r.ht/extract_line.cpp
===================================================================
--- grass-addons/grass7/raster/r.ht/extract_line.cpp (rev 0)
+++ grass-addons/grass7/raster/r.ht/extract_line.cpp 2012-06-20 08:00:10 UTC (rev 52153)
@@ -0,0 +1,257 @@
+#include "extract_line.h"
+
+#include "matrix.h"
+
+#include <stdio.h>
+
+#include <cmath>
+#include <vector>
+#include <list>
+
+using namespace matrix;
+
+struct Point
+{
+ int x;
+ int y;
+};
+
+double segmentLenght(const Segment& segment)
+{
+ double x = segment.first.first - segment.second.first;
+ double y = segment.first.second - segment.second.second;
+ return std::sqrt(x*x + y*y);
+}
+
+bool checkBounds(int y, int x, int cols, int rows)
+{
+ if ( y < 0 || y >= cols || x < 0 || x >= rows )
+ return false;
+ return true;
+}
+
+bool isData(const Matrix &I, const std::vector<int> &y, const std::vector<int> &x, const int cols, const int rows)
+{
+ for (size_t k = 0; k < y.size(); k++)
+ {
+ if (checkBounds(y.at(k), x.at(k), cols, rows) && I(x.at(k), y.at(k)))
+ {
+ return true;
+ }
+ }
+ return false;
+}
+
+void remove_points(LineCoordinates& lineCoordinates, const std::vector<int>& indicesI, const std::vector<int>& indicesJ)
+{
+ for (size_t i = 0; i < indicesI.size(); ++i)
+ {
+ std::pair<int, int> coords;
+ coords.second = indicesI[i];
+ coords.first = indicesJ[i];
+ lineCoordinates.remove(coords);
+ }
+}
+
+void find_indices(std::vector<int>& indicesI, std::vector<int>& indicesJ,
+ const int shift, const bool xflag, const int x, const int y,
+ const int line_width)
+{
+ if (xflag)
+ {
+ int tmp1, tmp2;
+ indicesJ[0] = x;
+ tmp1 = y >> shift;
+ tmp2 = tmp1;
+ indicesI[0] = tmp1;
+ for (int l = 1; l < line_width; l += 2)
+ {
+ indicesJ[l] = x;
+ indicesJ[l + 1] = x;
+ indicesI[l] = --tmp1;
+ indicesI[l + 1] = ++tmp2;
+ }
+ }
+ else
+ {
+ int tmp1, tmp2;
+ indicesI[0] = y;
+ tmp1 = x >> shift;
+ tmp2 = tmp1;
+ indicesJ[0] = tmp1;
+ for (int l = 1; l < line_width; l += 2)
+ {
+ indicesI[l] = y;
+ indicesI[l + 1] = y;
+ indicesJ[l] = --tmp1;
+ indicesJ[l + 1] = ++tmp2;
+ }
+ }
+}
+
+bool segmentContainsPoint(Segment segment, std::pair<int, int> p, bool xFlag, int tol)
+{
+ int min, max;
+ if (xFlag)
+ {
+ min = segment.first.first;
+ max = segment.second.first;
+ if (segment.first.first > segment.second.first)
+ {
+ min = segment.second.first;
+ max = segment.first.first;
+ }
+
+ if (p.first > min + tol && p.first <= max - tol)
+ return true;
+ }
+ else
+ {
+ min = segment.first.second;
+ max = segment.second.second;
+ if (segment.first.second > segment.second.second)
+ {
+ min = segment.second.second;
+ max = segment.first.second;
+ }
+ if (p.second > min + tol && p.second <= max - tol)
+ return true;
+ }
+ return false;
+}
+
+void extract(const Matrix& I,
+ float orient, const int lineGap, const int lineLength,
+ LineCoordinates lineCoordinates, SegmentList& segments)
+{
+ const int rows = I.rows ();
+ const int cols = I.columns ();
+
+ float rho = 1.0;
+ float irho = 1./rho;
+
+ bool xflag;
+
+ std::vector<Point> line_end(2);
+ SegmentList newSegments;
+
+ // from the current point walk in each direction
+ // along the found line and extract the line segment
+ float a = -(float)(sin(orient) * irho); //-trigtab[theta_n*2+1];
+ float b = (float)(cos(orient) * irho); //trigtab[theta_n*2];
+
+ int dx0, dy0;
+ int line_width = 3;
+ const int shift = 16;
+
+ std::vector<int> indicesI(line_width);
+ std::vector<int> indicesJ(line_width);
+
+ int lineNum = 0;
+ while (!lineCoordinates.empty())
+ {
+ lineNum++;
+ std::pair<int, int> startCoordinates = lineCoordinates.front();
+
+ int x0 = startCoordinates.first;
+ int y0 = startCoordinates.second;
+
+ if ( fabs(a) > fabs(b) )
+ {
+ xflag = 1;
+ dx0 = a > 0 ? 1 : -1;
+ dy0 = (int)( b*(1 << shift)/fabs(a) + 0.5);
+ y0 = (y0 << shift) + (1 << (shift-1));
+ }
+ else
+ {
+ xflag = 0;
+ dy0 = b > 0 ? 1 : -1;
+ dx0 = (int)( a*(1 << shift)/fabs(b) + 0.5);
+ x0 = (x0 << shift) + (1 << (shift-1));
+ }
+
+ int numOfGaps = 0;
+ bool useLine = true;
+ for ( int k = 0; k < 2; k++ )
+ {
+ int gap = 0, x = x0, y = y0, dx = dx0, dy = dy0;
+
+ if ( k > 0 )
+ dx = -dx, dy = -dy;
+
+ // walk along the line using fixed-point arithmetics,
+ // stop at the image border or in case of too big gap
+ for ( ;; x += dx, y += dy )
+ {
+ find_indices(indicesI, indicesJ, shift, xflag, x, y, line_width);
+
+ remove_points(lineCoordinates, indicesI, indicesJ);
+
+ // for each non-zero point:
+ // update line end,
+ // clear the mask element
+ // reset the gap
+ if ( isData(I, indicesI, indicesJ, cols, rows) )
+ {
+ if (gap > 5)
+ {
+ numOfGaps++;
+ if (numOfGaps > 5)
+ {
+ useLine = false;
+ break;
+ }
+ }
+ gap = 0;
+ line_end[k].y = indicesI[0];
+ line_end[k].x = indicesJ[0];
+ }
+ else if ( ++gap > lineGap )
+ {
+ break;
+ }
+ }
+ }
+
+ bool good_line = abs(line_end[1].x - line_end[0].x) >= lineLength ||
+ abs(line_end[1].y - line_end[0].y) >= lineLength;
+
+ if (good_line && useLine)
+ {
+ Segment newSegment;
+ newSegment.first.first = line_end[0].x;
+ newSegment.first.second = line_end[0].y;
+ newSegment.second.first = line_end[1].x;
+ newSegment.second.second = line_end[1].y;
+ int limit = 5;
+ bool add = true;
+ if (!newSegments.empty())
+ {
+ add = false;
+ const Segment lastSegment = newSegments.back();
+ if (!segmentContainsPoint(lastSegment, newSegment.first, xflag, limit) &&
+ !segmentContainsPoint(lastSegment, newSegment.second, xflag, limit) &&
+ !segmentContainsPoint(newSegment, lastSegment.first, xflag, limit) &&
+ !segmentContainsPoint(newSegment, lastSegment.second, xflag, limit))
+ {
+ add = true;
+ }
+
+ if (add == false)
+ {
+ if (segmentLenght(lastSegment) < segmentLenght(newSegment))
+ {
+ newSegments.pop_back();
+ add = true;
+ }
+ }
+ }
+ if (add)
+ {
+ newSegments.push_back(newSegment);
+ }
+ }
+ }
+ segments.insert(segments.end(), newSegments.begin(), newSegments.end());
+}
Added: grass-addons/grass7/raster/r.ht/extract_line.h
===================================================================
--- grass-addons/grass7/raster/r.ht/extract_line.h (rev 0)
+++ grass-addons/grass7/raster/r.ht/extract_line.h 2012-06-20 08:00:10 UTC (rev 52153)
@@ -0,0 +1,16 @@
+#ifndef EXTRACT_LINE_H
+#define EXTRACT_LINE_H
+
+#include "matrix.h"
+
+#include <vector>
+#include <list>
+
+typedef std::pair<int, int> Coordinates;
+typedef std::pair< Coordinates, Coordinates > Segment;
+typedef std::vector< Segment > SegmentList;
+typedef std::list< std::pair<int, int> > LineCoordinates;
+
+void extract(const matrix::Matrix& I, const float orient, const int lineGap, const int lineLength, LineCoordinates lineCoordinates, SegmentList& segments);
+
+#endif // EXTRACT_LINE_H
Added: grass-addons/grass7/raster/r.ht/hough.cpp
===================================================================
--- grass-addons/grass7/raster/r.ht/hough.cpp (rev 0)
+++ grass-addons/grass7/raster/r.ht/hough.cpp 2012-06-20 08:00:10 UTC (rev 52153)
@@ -0,0 +1,179 @@
+#include "houghtransform.h"
+
+#include"extract_line.h"
+#include "matrix.h"
+
+extern "C" {
+namespace grass {
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+#include <grass/gmath.h>
+}
+}
+
+using grass::DCELL;
+using grass::CELL;
+using grass::G__calloc;
+using grass::Cell_head;
+using grass::Map_info;
+
+using grass::Vect_new_cats_struct;
+using grass::Vect_new_line_struct;
+
+using grass::Rast_allocate_d_input_buf;
+using grass::Rast_open_old;
+using grass::Rast_get_row;
+using grass::Rast_close;
+
+using grass::G_gettext;
+using grass::G_fatal_error;
+using grass::G_debug;
+using grass::G_free;
+
+/** Loads map into memory.
+
+ \param[out] mat map in a matrix (row order), field have to be allocated
+ */
+void read_raster_map(const char *name, const char *mapset, int nrows,
+ int ncols, Matrix& mat)
+{
+
+ int r, c;
+
+ int map_fd;
+
+ DCELL *row_buffer;
+
+ DCELL cell_value;
+
+ row_buffer = Rast_allocate_d_input_buf();
+
+ /* 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(r, c) = cell_value;
+ else
+ mat(r, c) = 0.0;
+ }
+ }
+ G_free(row_buffer);
+
+ Rast_close(map_fd);
+}
+
+/**
+ \param cellhd raster map header, used for converting rows/cols to easting/northing
+ */
+void create_vector_map(const char * name, const SegmentList& segments,
+ const Cell_head* cellhd)
+{
+ struct Map_info Map;
+ Vect_open_new(&Map, name, 0);
+
+ struct grass::line_cats *Cats;
+ Cats = Vect_new_cats_struct();
+
+ for (size_t i = 0; i < segments.size(); ++i)
+ {
+ const Segment& seg = segments[i];
+
+ struct grass::line_pnts *points = Vect_new_line_struct();
+
+ double y1 = Rast_row_to_northing(seg.first.first, cellhd);
+ double x1 = Rast_col_to_easting(seg.first.second, cellhd);
+ double y2 = Rast_row_to_northing(seg.second.first, cellhd);
+ double x2 = Rast_col_to_easting(seg.second.second, cellhd);
+
+
+ Vect_cat_set(Cats, 1, i+1); // cat is segment number (counting from one)
+
+ Vect_append_point(points, x1, y1, 0);
+ Vect_append_point(points, x2, y2, 0);
+
+ Vect_write_line(&Map, GV_LINE, points, Cats);
+
+ }
+
+ Vect_build(&Map);
+ Vect_close(&Map);
+}
+
+void extract_line_segments(const Matrix &I,
+ const HoughTransform::Peaks& peaks,
+ const HoughTransform::TracebackMap& houghMap,
+ int gap,
+ int minSegmentLength,
+ SegmentList& segments)
+{
+ for (size_t i = 0; i < peaks.size(); ++i)
+ {
+ const HoughTransform::Peak& peak = peaks[i];
+
+ double theta = peak.coordinates.second;
+
+ HoughTransform::TracebackMap::const_iterator coordsIt =
+ houghMap.find(peak.coordinates);
+ if (coordsIt != houghMap.end())
+ {
+ const HoughTransform::CoordinatesList& lineCoordinates = coordsIt->second;
+
+ extract(I, (theta-90)/180*M_PI, gap, minSegmentLength, lineCoordinates, segments);
+ }
+ else
+ {
+ // logic error
+ }
+ }
+}
+
+void hough_peaks(int maxPeaks, int threshold, int sizeOfNeighbourhood,
+ int gap, int minSegmentLength,
+ const char *name, const char *mapset, size_t nrows, size_t ncols,
+ const char *anglesMapName, int angleWidth,
+ const char *result)
+{
+ Matrix I(nrows, ncols);
+ read_raster_map(name, mapset, nrows, ncols, I);
+
+ HoughTransform hough(I);
+
+ if (anglesMapName != NULL)
+ {
+ Matrix angles(nrows, ncols);
+ read_raster_map(anglesMapName, mapset, nrows, ncols, angles);
+ hough.compute(angles, angleWidth);
+ }
+ else
+ {
+ hough.compute();
+ }
+
+ hough.findPeaks(maxPeaks, threshold, sizeOfNeighbourhood);
+
+ const HoughTransform::Peaks& peaks = hough.getPeaks();
+ const HoughTransform::TracebackMap& houghMap = hough.getHoughMap();
+ SegmentList segments;
+
+ extract_line_segments(I, peaks, houghMap, gap, minSegmentLength, segments);
+
+ Cell_head cellhd;
+ Rast_get_window(&cellhd);
+ create_vector_map(result, segments, &cellhd);
+}
Added: grass-addons/grass7/raster/r.ht/hough.h
===================================================================
--- grass-addons/grass7/raster/r.ht/hough.h (rev 0)
+++ grass-addons/grass7/raster/r.ht/hough.h 2012-06-20 08:00:10 UTC (rev 52153)
@@ -0,0 +1,11 @@
+#ifndef HOUGH_H
+#define HOUGH_H
+
+extern "C" {
+#include <grass/gis.h>
+#include <grass/raster.h>
+}
+
+void hough_peaks(int maxPeaks, int threshold, int sizeOfNeighbourhood, int gap, int minSegmentLength, const char * name, const char* mapset, size_t nrows, size_t ncols, const char * angleMapName, int angleWidth, const char * result);
+
+#endif // HOUGH_H
Added: grass-addons/grass7/raster/r.ht/houghtransform.cpp
===================================================================
--- grass-addons/grass7/raster/r.ht/houghtransform.cpp (rev 0)
+++ grass-addons/grass7/raster/r.ht/houghtransform.cpp 2012-06-20 08:00:10 UTC (rev 52153)
@@ -0,0 +1,239 @@
+#include "houghtransform.h"
+
+/* helpers */
+
+struct SortByX
+{
+ bool operator()(const HoughTransform::Coordinates& c1, const HoughTransform::Coordinates& c2) const
+ {
+
+ return c1.first < c2.first;
+ }
+};
+
+struct SortByY
+{
+ bool operator()(const HoughTransform::Coordinates& c1, const HoughTransform::Coordinates& c2) const
+ {
+
+ return c1.second < c2.second;
+ }
+};
+
+
+/* ctors */
+
+HoughTransform::HoughTransform(const Matrix &matrix):
+ mOriginalMatrix(matrix)
+{
+ mNumR = mOriginalMatrix.rows();
+ mNumC = mOriginalMatrix.columns();
+
+ mThetas = ColumnVector(matrix::Range(-M_PI/2.0, M_PI/2.0, M_PI/180.0).matrix_value ());
+
+ const float diag_length = std::sqrt(mNumR*mNumR + mNumC*mNumC);
+ mNumBins = ceil(diag_length) - 1;
+
+ mHoughMatrix = Matrix(mNumBins, mThetas.length(), 0.0);
+
+ c2 = ceil(mNumC/2.);
+ r2 = ceil(mNumR/2.);
+
+ bins0 = 1 - ceil(mNumBins/2.0);
+}
+
+/* functions */
+
+void HoughTransform::compute()
+{
+ for (int x = 0; x < mNumR; x++)
+ {
+ for (int y = 0; y < mNumC; y++)
+ {
+ if (mOriginalMatrix(x, y) == 1)
+ {
+ computeHoughForXY(x, y, 0, mThetas.length());
+ }
+ }
+ }
+
+}
+
+void HoughTransform::computeHoughForXY(int x, int y, size_t minIndex, size_t maxIndex)
+{
+ for (size_t i = minIndex; i < maxIndex; i++)
+ {
+ const double theta = mThetas(i);
+ const double cT = cos(theta);
+ const double sT = sin(theta);
+ const int rho = (int) floor(cT*(x - c2) + sT*(y - r2) + 0.5);
+ const int bin = (int) (rho - bins0);
+ if ((bin > 0) && (bin < mNumBins))
+ {
+ mHoughMatrix(bin, i)++;
+ mHoughMap[Coordinates(bin, i)].push_back(Coordinates(x, y));
+ mOriginalMap[Coordinates(x, y)].push_back(Coordinates(bin, i));
+ }
+ }
+}
+
+/**
+ Expecting angles to be in degrees in range [0,180).
+
+ Size of agles matrix must be the same as size of input edge matrix.
+ mNumR != angles.rows() && mNumC != angles.columns()
+ */
+void HoughTransform::compute(const Matrix& angles, double angleWith)
+{
+ for (int x = 0; x < mNumR; x++)
+ {
+ for (int y = 0; y < mNumC; y++)
+ {
+ if (mOriginalMatrix(x, y) == 1)
+ {
+ double angle = angles(x, y);
+ if (angle < 0)
+ angle = std::abs(angle);
+ else if (angle > 0)
+ angle = 180 - angle;
+
+ // converting angle [0,180) to index
+ // in internal table of angles
+ // assuming that table size is 180
+ // TODO: angleIndex = mThetas.length() * angle/180
+ int angleIndex = angle;
+ int angleShift = angleWith/2 + 0.5;
+
+ // FIXME: magic number
+ int minIndex = angleIndex - angleShift;
+ int maxIndex = angleIndex + angleShift;
+
+ if (minIndex < 0) {
+ computeHoughForXY(x, y, 0, maxIndex);
+ minIndex = 180 + minIndex;
+ computeHoughForXY(x, y, minIndex, 180);
+ }
+ else if (maxIndex > 180)
+ {
+ computeHoughForXY(x, y, minIndex, 180);
+ maxIndex = maxIndex - 180;
+ computeHoughForXY(x, y, 0, maxIndex);
+ }
+ else
+ {
+ computeHoughForXY(x, y, minIndex, maxIndex);
+ }
+ }
+ }
+ }
+}
+
+void HoughTransform::findPeaks(int maxPeakNumber, int threshold, int sizeOfNeighbourhood)
+{
+ int maxIt = mHoughMatrix.rows() * mHoughMatrix.columns();
+ int peaksFound = 0;
+ for (int i = 0; i < maxIt; i++)
+ {
+ Coordinates coordinates;
+ int maxValue = findMax(mHoughMatrix, coordinates);
+
+ if (maxValue < threshold || peaksFound == maxPeakNumber )
+ break;
+
+ CoordinatesList neighbours(neighbourhood(coordinates, sizeOfNeighbourhood));
+
+ Coordinates beginLine, endLine;
+
+ removePeakEffect(neighbours, beginLine, endLine);
+
+ Peak peak(coordinates, maxValue, beginLine, endLine);
+ mPeaks.push_back(peak);
+ peaksFound++;
+ }
+}
+
+HoughTransform::CoordinatesList HoughTransform::neighbourhood(Coordinates &coordinates, const int sizeOfNeighbourhood)
+{
+ CoordinatesList list;
+ const int x = coordinates.first;
+ const int y = coordinates.second;
+ for (int i = -sizeOfNeighbourhood; i <= sizeOfNeighbourhood; i++)
+ {
+ for (int j = -sizeOfNeighbourhood; j <= sizeOfNeighbourhood; j++)
+ {
+ list.push_back(Coordinates(x + i, y + j));
+ }
+ }
+ return list;
+}
+
+void HoughTransform::removePeakEffect(const CoordinatesList &neighbours, Coordinates &beginLine, Coordinates &endLine)
+{
+ CoordinatesList lineList;
+ for (CoordinatesList::const_iterator it = neighbours.begin(), end = neighbours.end(); it != end; ++it)
+ {
+ CoordinatesList &toRemove = mHoughMap[*it];
+
+ for (CoordinatesList::const_iterator it1 = toRemove.begin(), end = toRemove.end(); it1 != end; ++it1)
+ {
+ if ( mOriginalMap.find(*it1) != mOriginalMap.end() )
+ {
+ CoordinatesList &toDecrease = mOriginalMap[*it1];
+
+ for (CoordinatesList::const_iterator it2 = toDecrease.begin(), end = toDecrease.end(); it2 != end; ++it2)
+ {
+ mHoughMatrix(it2->first, it2->second)--;
+ }
+ mOriginalMap.erase(*it1);
+
+ lineList.push_back(*it1);
+ }
+ }
+ }
+ int angle = neighbours.front().second;
+ // what to do if find endpoints finds nothing reasonable?
+ findEndPoints( lineList, beginLine, endLine, angle);
+}
+
+bool HoughTransform::findEndPoints(CoordinatesList list, Coordinates &beginLine, Coordinates &endLine, const int angle )
+{
+ if (angle > 45 && angle <= 135)
+ {
+ list.sort(SortByY());
+ }
+ else
+ {
+ list.sort(SortByX());
+ }
+
+ beginLine = list.front();
+ endLine = list.back();
+
+ return true;
+}
+
+int HoughTransform::findMax(const Matrix& matrix, Coordinates &coordinates)
+{
+ size_t rowIndex = 0;
+ size_t colIndex;
+ std::vector<size_t> colIndexes;
+
+ std::vector<double> rowMax = matrix.row_max(colIndexes);
+
+ int max = 0;
+ for (size_t i = 0; i < rowMax.size(); i++)
+ {
+ int tmp = rowMax[i];
+ if (tmp > max)
+ {
+ max = tmp;
+ rowIndex = i;
+ }
+ }
+ colIndex = colIndexes[rowIndex];
+
+ coordinates.first = rowIndex;
+ coordinates.second = colIndex;
+
+ return max;
+}
Added: grass-addons/grass7/raster/r.ht/houghtransform.h
===================================================================
--- grass-addons/grass7/raster/r.ht/houghtransform.h (rev 0)
+++ grass-addons/grass7/raster/r.ht/houghtransform.h 2012-06-20 08:00:10 UTC (rev 52153)
@@ -0,0 +1,103 @@
+#ifndef HOUGHTRANSFORM_H
+#define HOUGHTRANSFORM_H
+
+#include "matrix.h"
+
+#include <cmath>
+#include <vector>
+#include <list>
+#include <set>
+#include <map>
+#include <limits>
+#include <algorithm>
+
+#include <stdio.h>
+
+using namespace matrix;
+
+class HoughTransform
+{
+public:
+
+ /* typedefs */
+
+ template <typename T>
+ class MyVector: public std::vector<T>
+ {
+ public:
+ template <typename Comp>
+ void sort(Comp comp)
+ {
+ std::sort(this->begin(), this->end(), comp);
+ }
+ operator std::list<T>() const
+ {
+ std::list<T> tmp;
+ tmp.insert(tmp.begin(), this->begin(), this->end());
+ return tmp;
+ }
+ };
+
+ typedef std::pair<int, int> Coordinates;
+ typedef MyVector<Coordinates> CoordinatesList;
+ typedef std::map<Coordinates, CoordinatesList> TracebackMap;
+
+ struct Peak
+ {
+ Peak(Coordinates coordinates, int value, Coordinates begin, Coordinates end)
+ : coordinates(coordinates), value(value), beginLine(begin), endLine(end)
+ {}
+ Coordinates coordinates;
+ int value;
+ Coordinates beginLine;
+ Coordinates endLine;
+ };
+
+ typedef std::vector<Peak> Peaks;
+
+ /* functions */
+
+ HoughTransform(const Matrix &matrix);
+
+ void compute();
+ void compute(const Matrix& angles, double angleWith);
+ void findPeaks(int maxPeakNumber, int threshold,
+ const int sizeOfNeighbourhood);
+
+ /* getters */
+
+ const Matrix & getHoughMatrix() { return mHoughMatrix; }
+ const Matrix & getOrigMatrix() { return mOriginalMatrix; }
+ const Peaks & getPeaks() const { return mPeaks; }
+ const TracebackMap & getHoughMap() { return mHoughMap; }
+
+private:
+
+ /* functions */
+
+ CoordinatesList neighbourhood(Coordinates &coordinates, const int sizeOfNeighbourhood);
+ void removePeakEffect(const CoordinatesList &neighbours, Coordinates &beginLine, Coordinates &endLine);
+ bool findEndPoints(CoordinatesList list, Coordinates &beginLine, Coordinates &endLine, const int angle );
+ int findMax(const Matrix& matrix, Coordinates &coordinates);
+ void addToMaps(const std::vector<std::pair<Coordinates, Coordinates> > &pairs);
+ void computeHoughForXY(int x, int y, size_t minIndex, size_t maxIndex);
+
+ /* data members */
+
+ const Matrix& mOriginalMatrix;
+ Matrix mHoughMatrix;
+ TracebackMap mOriginalMap;
+ TracebackMap mHoughMap;
+ Peaks mPeaks;
+
+ int mNumR;
+ int mNumC;
+ ColumnVector mThetas;
+ int mNumBins;
+ int c2;
+ int r2;
+ int bins0;
+};
+
+
+#endif // HOUGHTRANSFORM_H
Added: grass-addons/grass7/raster/r.ht/main.cpp
===================================================================
--- grass-addons/grass7/raster/r.ht/main.cpp (rev 0)
+++ grass-addons/grass7/raster/r.ht/main.cpp 2012-06-20 08:00:10 UTC (rev 52153)
@@ -0,0 +1,171 @@
+
+/****************************************************************************
+ *
+ * MODULE: r.ht
+ * AUTHOR(S): Anna Kratochvilova - kratochanna gmail.com
+ * Vaclav Petras - wenzeslaus gmail.com
+ *
+ * PURPOSE: Detects line segments using Hough transform.
+ *
+ * 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 "hough.h"
+
+extern "C" {
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+#include <grass/gmath.h>
+}
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+/**
+
+ \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 */
+
+ char *result; /* output raster name */
+
+ int nrows, ncols;
+
+ struct GModule *module; /* GRASS module for parsing arguments */
+
+ /* options */
+ struct Option *input, *output, *anglesOption,
+ *maxLinesOption, *maxGapOption, *minSegmentLengthOption,
+ *angleWidthOption;
+
+ /* 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(_("hought"));
+ G_add_keyword(_(""));
+ 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_V_OUTPUT);
+
+ anglesOption = G_define_standard_option(G_OPT_R_INPUT);
+ anglesOption->key = "angles";
+ anglesOption->required = NO;
+ anglesOption->description = _("Approximate number of line segments.");
+
+ // this option will become max peaks number to find in HT
+ maxLinesOption = G_define_option();
+ maxLinesOption->key = "lines_number";
+ maxLinesOption->type = TYPE_INTEGER;
+ maxLinesOption->required = NO;
+ maxLinesOption->multiple = NO;
+ maxLinesOption->description = _("Approximate number of line segments."
+ " Actually, this option will become"
+ " maximal number of line candidates"
+ " detected in Hough transform image."
+ " Final number of line segments can be"
+ " smaller or greater."
+ );
+ maxLinesOption->answer = const_cast<char *>("20");
+
+ maxGapOption = G_define_option();
+ maxGapOption->key = "max_gap";
+ maxGapOption->type = TYPE_INTEGER;
+ maxGapOption->required = NO;
+ maxGapOption->multiple = NO;
+ maxGapOption->description = _("Maximum gap in pixels");
+ maxGapOption->answer = const_cast<char *>("5");
+
+ minSegmentLengthOption = G_define_option();
+ minSegmentLengthOption->key = "min_segment_length";
+ minSegmentLengthOption->type = TYPE_INTEGER;
+ minSegmentLengthOption->required = NO;
+ minSegmentLengthOption->multiple = NO;
+ minSegmentLengthOption->description = _("Minimal length of line segment");
+ minSegmentLengthOption->answer = const_cast<char *>("50");
+
+ angleWidthOption = G_define_option();
+ angleWidthOption->key = "angle_width";
+ angleWidthOption->type = TYPE_INTEGER;
+ angleWidthOption->required = NO;
+ angleWidthOption->multiple = NO;
+ angleWidthOption->description = _("Width of circle sector (only when you provide angle map).");
+ angleWidthOption->answer = const_cast<char *>("5");
+
+ /* options and flags parser */
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ /* stores options and flags to variables */
+ result = output->answer;
+ name = input->answer;
+
+ int maxPeaks = atoi(maxLinesOption->answer);
+ int threshold = 10;
+ int gap = atoi(maxGapOption->answer);
+ int minSegmentLength = atoi(minSegmentLengthOption->answer);
+ int sizeOfNeighbourhood = 1;
+
+ int angleWidth = atoi(angleWidthOption->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);
+
+ // 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 = Rast_window_rows();
+
+ ncols = Rast_window_cols();
+
+ /* **** */
+
+ hough_peaks(maxPeaks, threshold, sizeOfNeighbourhood, gap, minSegmentLength, name, mapset, nrows, ncols, anglesOption->answer, angleWidth, result);
+
+ /* **** */
+
+ /* memory cleanup */
+ G_free(name);
+
+ /* add command line incantation to history file */
+ // Rast_short_history(templName, "raster", &history);
+ // Rast_command_history(&history);
+ // Rast_write_history(templName, &history);
+
+ exit(EXIT_SUCCESS);
+}
Added: grass-addons/grass7/raster/r.ht/matrix.cpp
===================================================================
Added: grass-addons/grass7/raster/r.ht/matrix.h
===================================================================
--- grass-addons/grass7/raster/r.ht/matrix.h (rev 0)
+++ grass-addons/grass7/raster/r.ht/matrix.h 2012-06-20 08:00:10 UTC (rev 52153)
@@ -0,0 +1,221 @@
+#ifndef MATRIX_H
+#define MATRIX_H
+
+#include <vector>
+#include <algorithm>
+#include <limits>
+#include <cmath>
+
+/* mimics Octave matrices and other classes */
+
+namespace matrix {
+
+class Matrix
+{
+public:
+ Matrix(size_t r, size_t c, double val = 0)
+ {
+ resize(r, c, val);
+ }
+ Matrix() {}
+
+ void resize(size_t r, size_t c, double val = 0)
+ {
+ std::vector<double> row;
+ row.resize(c, val);
+ mat.resize(r, row);
+ }
+
+ double& operator ()(size_t r, size_t c)
+ {
+ if (r > rows() || c > columns())
+ {
+ // FIXME: fatal error
+ }
+ return mat[r][c];
+ }
+ const double& operator ()(size_t r, size_t c) const
+ {
+ return mat[r][c];
+ }
+ size_t rows() const { return mat.size(); }
+ size_t columns() const
+ {
+ if (!mat.empty())
+ return mat[0].size();
+ return 0;
+ }
+
+ std::vector<double> row_max(std::vector<size_t>& colIndexes) const
+ {
+ std::vector<double> ret;
+ for (size_t i = 0; i < rows(); ++i)
+ {
+ std::vector<double>::const_iterator maxe = std::max_element(mat[i].begin(), mat[i].end());
+ double max = *maxe;
+ size_t maxi = maxe - mat[i].begin();
+
+ ret.push_back(max);
+ colIndexes.push_back(maxi);
+ }
+ return ret;
+ }
+
+private:
+ std::vector< std::vector<double> > mat;
+};
+
+class ColumnVector
+{
+public:
+ ColumnVector(Matrix mat)
+ {
+ for (size_t i = 0; i < mat.columns(); ++i)
+ {
+ vec.push_back(mat(0, i));
+ }
+ }
+ ColumnVector() {}
+
+ double& operator ()(size_t i)
+ {
+ return vec[i];
+ }
+ const double& operator ()(size_t i) const
+ {
+ return vec[i];
+ }
+ size_t length() { return vec.size(); }
+
+private:
+ std::vector<double> vec;
+};
+
+class RowVector
+{
+public:
+ RowVector(Matrix mat)
+ {
+ for (size_t i = 0; i < mat.columns(); ++i)
+ {
+ vec.push_back(mat(0, i));
+ }
+ }
+ RowVector() {}
+
+ double& operator ()(size_t i)
+ {
+ return vec[i];
+ }
+ const double& operator ()(size_t i) const
+ {
+ return vec[i];
+ }
+ size_t length() { return vec.size(); }
+
+ RowVector operator -(double d)
+ {
+ RowVector ret(*this);
+ for (size_t i = 0; i < ret.length(); ++i)
+ {
+ ret(i) = ret(i) - d;
+ }
+ return ret;
+ }
+
+private:
+ std::vector<double> vec;
+};
+
+class Range
+{
+public:
+ Range (double b, double l, double i)
+ : rng_base(b), rng_limit(l), rng_inc(i), rng_nelem(nelem_internal())
+ { }
+
+ Range (double b, double l)
+ : rng_base(b), rng_limit(l), rng_inc(1),
+ rng_nelem(nelem_internal()) { }
+
+ Matrix matrix_value() const
+ {
+ Matrix cache;
+ //if (rng_nelem > 0 && cache.nelem () == 0)
+ //{
+ cache.resize(1, rng_nelem);
+ double b = rng_base;
+ double increment = rng_inc;
+ for (size_t i = 0; i < rng_nelem; i++)
+ cache(0, i) = b + i * increment;
+
+ // On some machines (x86 with extended precision floating point
+ // arithmetic, for example) it is possible that we can overshoot
+ // the limit by approximately the machine precision even though
+ // we were very careful in our calculation of the number of
+ // elements.
+
+ if ((rng_inc > 0 && cache(0, rng_nelem-1) > rng_limit)
+ || (rng_inc < 0 && cache(0, rng_nelem-1) < rng_limit))
+ cache(0, rng_nelem-1) = rng_limit;
+ // }
+ return cache;
+ }
+
+private:
+ /* from Octave source codes */
+ size_t nelem_internal() const
+ {
+ size_t retval = -1;
+
+ if (rng_inc == 0
+ || (rng_limit > rng_base && rng_inc < 0)
+ || (rng_limit < rng_base && rng_inc > 0))
+ {
+ retval = 0;
+ }
+ else
+ {
+ // double ct = 3.0 * std::numeric_limits<double>::epsilon(); // FIXME: not used
+ // was DBL_EPSILON;
+
+ // double tmp = tfloor ((rng_limit - rng_base + rng_inc) / rng_inc, ct); // octave code
+ double tmp = floor ((rng_limit - rng_base + rng_inc) / rng_inc);
+
+ size_t n_elt = (tmp > 0.0 ? static_cast<size_t> (tmp) : 0);
+
+ // If the final element that we would compute for the range is
+ // equal to the limit of the range, or is an adjacent floating
+ // point number, accept it. Otherwise, try a range with one
+ // fewer element. If that fails, try again with one more
+ // element.
+ //
+ // I'm not sure this is very good, but it seems to work better than
+ // just using tfloor as above. For example, without it, the
+ // expression 1.8:0.05:1.9 fails to produce the expected result of
+ // [1.8, 1.85, 1.9].
+
+ // octave code
+ // if (! teq (rng_base + (n_elt - 1) * rng_inc, rng_limit))
+ // {
+ // if (teq (rng_base + (n_elt - 2) * rng_inc, rng_limit))
+ // n_elt--;
+ // else if (teq (rng_base + n_elt * rng_inc, rng_limit))
+ // n_elt++;
+ // }
+
+ retval = (n_elt >= std::numeric_limits<size_t>::max () - 1) ? -1 : n_elt;
+ }
+
+ return retval;
+ }
+
+ double rng_base;
+ double rng_limit;
+ double rng_inc;
+ size_t rng_nelem;
+};
+
+}
+
+#endif // MATRIX_H
Added: grass-addons/grass7/raster/r.ht/r.ht.html
===================================================================
--- grass-addons/grass7/raster/r.ht/r.ht.html (rev 0)
+++ grass-addons/grass7/raster/r.ht/r.ht.html 2012-06-20 08:00:10 UTC (rev 52153)
@@ -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