[GRASS-SVN] r52180 - grass-addons/grass7/raster/r.ht

svn_grass at osgeo.org svn_grass at osgeo.org
Thu Jun 21 03:00:25 PDT 2012


Author: wenzeslaus
Date: 2012-06-21 03:00:25 -0700 (Thu, 21 Jun 2012)
New Revision: 52180

Added:
   grass-addons/grass7/raster/r.ht/houghparameters.h
   grass-addons/grass7/raster/r.ht/linesegmentsextractor.cpp
   grass-addons/grass7/raster/r.ht/linesegmentsextractor.h
Removed:
   grass-addons/grass7/raster/r.ht/extract_line.cpp
   grass-addons/grass7/raster/r.ht/extract_line.h
Modified:
   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
Log:
hough: new option line_width, colors for ht image, extract segment class

Deleted: grass-addons/grass7/raster/r.ht/extract_line.cpp
===================================================================
--- grass-addons/grass7/raster/r.ht/extract_line.cpp	2012-06-21 09:23:07 UTC (rev 52179)
+++ grass-addons/grass7/raster/r.ht/extract_line.cpp	2012-06-21 10:00:25 UTC (rev 52180)
@@ -1,257 +0,0 @@
-#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, int gapSize, int maxNumOfGaps, 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 > gapSize)
-                    {
-                        numOfGaps++;
-                        if (numOfGaps > maxNumOfGaps)
-                        {
-                            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());
-}

Deleted: grass-addons/grass7/raster/r.ht/extract_line.h
===================================================================
--- grass-addons/grass7/raster/r.ht/extract_line.h	2012-06-21 09:23:07 UTC (rev 52179)
+++ grass-addons/grass7/raster/r.ht/extract_line.h	2012-06-21 10:00:25 UTC (rev 52180)
@@ -1,16 +0,0 @@
-#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, int gapSize, int maxNumOfGaps, const int lineGap, const int lineLength, LineCoordinates lineCoordinates, SegmentList& segments);
-
-#endif // EXTRACT_LINE_H

Modified: grass-addons/grass7/raster/r.ht/hough.cpp
===================================================================
--- grass-addons/grass7/raster/r.ht/hough.cpp	2012-06-21 09:23:07 UTC (rev 52179)
+++ grass-addons/grass7/raster/r.ht/hough.cpp	2012-06-21 10:00:25 UTC (rev 52180)
@@ -1,6 +1,6 @@
 #include "houghtransform.h"
 
-#include"extract_line.h"
+#include"linesegmentsextractor.h"
 #include "matrix.h"
 
 extern "C" {
@@ -39,6 +39,10 @@
 using grass::G_debug;
 using grass::G_free;
 
+using grass::Colors;
+using grass::FPRange; // FIXME: DCELL/CELL
+using grass::G_mapset;
+
 /** Loads map into memory.
 
   \param[out] mat map in a matrix (row order), field have to be allocated
@@ -85,6 +89,18 @@
     Rast_close(map_fd);
 }
 
+void apply_hough_colors_to_map(const char *name)
+{
+    struct Colors colors;
+    struct FPRange range;
+    DCELL min, max;
+
+    Rast_read_fp_range(name, G_mapset(), &range);
+    Rast_get_fp_range_min_max(&range, &min, &max);
+    Rast_make_grey_scale_colors(&colors, min, max);
+    Rast_write_colors(name, G_mapset(), &colors);
+}
+
 void create_raster_map(const char *name, struct Cell_head *window, const Matrix& mat)
 {
     struct Cell_head original_window;
@@ -145,25 +161,24 @@
     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();
+        struct grass::line_pnts *points = Vect_new_line_struct(); // FIXME: some destroy
 
         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);
@@ -173,8 +188,7 @@
 void extract_line_segments(const Matrix &I,
                            const HoughTransform::Peaks& peaks,
                            const HoughTransform::TracebackMap& houghMap,
-                           int gapSize, int maxNumOfGaps, int gap,
-                           int minSegmentLength,
+                           ExtractParametres extractParametres,
                            SegmentList& segments)
 {
     for (size_t i = 0; i < peaks.size(); ++i)
@@ -189,7 +203,9 @@
         {
             const HoughTransform::CoordinatesList& lineCoordinates = coordsIt->second;
 
-            extract(I, (theta-90)/180*M_PI, gapSize, maxNumOfGaps, gap, minSegmentLength, lineCoordinates, segments);
+            LineSegmentsExtractor extractor(I, extractParametres);
+
+            extractor.extract(lineCoordinates, (theta-90)/180*M_PI, segments);
         }
         else
         {
@@ -198,8 +214,8 @@
     }
 }
 
-void hough_peaks(int maxPeaks, int threshold, double angleWith, int sizeOfNeighbourhood,
-                 int gapSize, int maxNumOfGaps, int gap, int minSegmentLength,
+void hough_peaks(HoughParametres houghParametres,
+                 ExtractParametres extractParametres,
                  const char *name, const char *mapset, size_t nrows, size_t ncols,
                  const char *anglesMapName,
                  const char *houghImageName, const char *result)
@@ -207,33 +223,34 @@
     Matrix I(nrows, ncols);
     read_raster_map(name, mapset, nrows, ncols, I);
 
-    HoughTransform hough(I);
+    HoughTransform hough(I, houghParametres);
 
     if (anglesMapName != NULL)
     {
         Matrix angles(nrows, ncols);
         read_raster_map(anglesMapName, mapset, nrows, ncols, angles);
-        hough.compute(angles, angleWith);
+        hough.compute(angles);
     }
     else
     {
         hough.compute();
     }
 
-    hough.findPeaks(maxPeaks, threshold, sizeOfNeighbourhood);
+    hough.findPeaks();
 
     if (houghImageName != NULL)
     {
         struct Cell_head window;
         Rast_get_cellhd(name, "", &window);
         create_raster_map(houghImageName, &window, hough.getHoughMatrix());
+        apply_hough_colors_to_map(houghImageName);
     }
 
     const HoughTransform::Peaks& peaks = hough.getPeaks();
     const HoughTransform::TracebackMap& houghMap = hough.getHoughMap();
     SegmentList segments;
 
-    extract_line_segments(I, peaks, houghMap, gapSize, maxNumOfGaps, gap, minSegmentLength, segments);
+    extract_line_segments(I, peaks, houghMap, extractParametres, segments);
 
     Cell_head cellhd;
     Rast_get_window(&cellhd);

Modified: grass-addons/grass7/raster/r.ht/hough.h
===================================================================
--- grass-addons/grass7/raster/r.ht/hough.h	2012-06-21 09:23:07 UTC (rev 52179)
+++ grass-addons/grass7/raster/r.ht/hough.h	2012-06-21 10:00:25 UTC (rev 52180)
@@ -1,11 +1,13 @@
 #ifndef HOUGH_H
 #define HOUGH_H
 
+#include "houghparameters.h"
+
 extern "C" {
 #include <grass/gis.h>
 #include <grass/raster.h>
 }
 
-void hough_peaks(int maxPeaks, int threshold, double angleWith, int sizeOfNeighbourhood, int gapSize, int maxNumOfGaps, int gap, int minSegmentLength, const char *name, const char *mapset, size_t nrows, size_t ncols, const char *angleMapName, const char *houghImageName, const char *result);
+void hough_peaks(HoughParametres houghParametres, ExtractParametres extractParametres, const char *name, const char *mapset, size_t nrows, size_t ncols, const char *angleMapName, const char *houghImageName, const char *result);
 
 #endif // HOUGH_H

Added: grass-addons/grass7/raster/r.ht/houghparameters.h
===================================================================
--- grass-addons/grass7/raster/r.ht/houghparameters.h	                        (rev 0)
+++ grass-addons/grass7/raster/r.ht/houghparameters.h	2012-06-21 10:00:25 UTC (rev 52180)
@@ -0,0 +1,21 @@
+#ifndef HOUGHPARAMETERS_H
+#define HOUGHPARAMETERS_H
+
+struct HoughParametres
+{
+    double angleWidth;
+    int maxPeaksNum;
+    int threshold;
+    int sizeOfNeighbourhood;
+};
+
+struct ExtractParametres
+{
+    int gapSize;
+    int maxNumOfGaps;
+    int maxGap;
+    int lineLength;
+    int lineWidth;
+};
+
+#endif // HOUGHPARAMETERS_H

Modified: grass-addons/grass7/raster/r.ht/houghtransform.cpp
===================================================================
--- grass-addons/grass7/raster/r.ht/houghtransform.cpp	2012-06-21 09:23:07 UTC (rev 52179)
+++ grass-addons/grass7/raster/r.ht/houghtransform.cpp	2012-06-21 10:00:25 UTC (rev 52180)
@@ -23,8 +23,8 @@
 
 /* ctors */
 
-HoughTransform::HoughTransform(const Matrix &matrix):
-    mOriginalMatrix(matrix)
+HoughTransform::HoughTransform(const Matrix &matrix, const HoughParametres& parametrs)
+    : mOriginalMatrix(matrix), mParams(parametrs)
 {
     mNumR = mOriginalMatrix.rows();
     mNumC = mOriginalMatrix.columns();
@@ -192,18 +192,21 @@
     }
     int angle = neighbours.front().second;
     // what to do if find endpoints finds nothing reasonable?
-    findEndPoints( lineList, beginLine, endLine, angle);
+    findEndPoints(lineList, beginLine, endLine, angle);
 }
 
-bool HoughTransform::findEndPoints(CoordinatesList list, Coordinates &beginLine, Coordinates &endLine, const int angle )
+/** \param list[in, out] will be sorted by y cooridinate
+  if angle is in range (45, 135] or by x otherwise
+  */
+bool HoughTransform::findEndPoints(CoordinatesList& list, Coordinates &beginLine, Coordinates &endLine, const int angle)
 {
     if (angle > 45 && angle <= 135)
     {
-        list.sort(SortByY());
+        std::sort(list.begin(), list.end(), SortByY());
     }
     else
     {
-        list.sort(SortByX());
+        std::sort(list.begin(), list.end(), SortByX());
     }
 
     beginLine = list.front();

Modified: grass-addons/grass7/raster/r.ht/houghtransform.h
===================================================================
--- grass-addons/grass7/raster/r.ht/houghtransform.h	2012-06-21 09:23:07 UTC (rev 52179)
+++ grass-addons/grass7/raster/r.ht/houghtransform.h	2012-06-21 10:00:25 UTC (rev 52180)
@@ -1,6 +1,7 @@
 #ifndef HOUGHTRANSFORM_H
 #define HOUGHTRANSFORM_H
 
+#include "houghparameters.h"
 #include "matrix.h"
 
 #include <cmath>
@@ -19,17 +20,24 @@
 {
 public:
 
-    /* typedefs */
+    /* types */
 
+    /** This is class to provide compatible syntax with \c std::list. */
     template <typename T>
-    class MyVector: public std::vector<T>
+    class Vector: public std::vector<T>
     {
     public:
+
+        /**
+          Maybe it would be better to replace by an overload of sort function for list.
+          And use this function for list.
+          */
         template <typename Comp>
         void sort(Comp comp)
         {
             std::sort(this->begin(), this->end(), comp);
         }
+        /** only for convenience, with some compiler optimalisation overhead shouldn't be so high  */
         operator std::list<T>() const
         {
             std::list<T> tmp;
@@ -39,7 +47,7 @@
     };
 
     typedef std::pair<int, int> Coordinates;
-    typedef MyVector<Coordinates> CoordinatesList;
+    typedef Vector<Coordinates> CoordinatesList;
     typedef std::map<Coordinates, CoordinatesList> TracebackMap;
 
     struct Peak
@@ -57,10 +65,18 @@
 
     /* functions */
 
-    HoughTransform(const Matrix &matrix);
+    HoughTransform(const Matrix &matrix, const HoughParametres &parametrs);
 
     void compute();
+    void compute(const Matrix& angles)
+    {
+        compute(angles, mParams.angleWidth);
+    }
     void compute(const Matrix& angles, double angleWith);
+    void findPeaks()
+    {
+        findPeaks(mParams.maxPeaksNum, mParams.threshold, mParams.sizeOfNeighbourhood);
+    }
     void findPeaks(int maxPeakNumber, int threshold,
                    const int sizeOfNeighbourhood);
 
@@ -75,11 +91,11 @@
 
     /* functions */
 
+    // some of these functions can be changed to non-member
     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 );
+    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 */
@@ -97,6 +113,8 @@
     int c2;
     int r2;
     int bins0;
+
+    HoughParametres mParams;
 };
 
 

Copied: grass-addons/grass7/raster/r.ht/linesegmentsextractor.cpp (from rev 52167, grass-addons/grass7/raster/r.ht/extract_line.cpp)
===================================================================
--- grass-addons/grass7/raster/r.ht/linesegmentsextractor.cpp	                        (rev 0)
+++ grass-addons/grass7/raster/r.ht/linesegmentsextractor.cpp	2012-06-21 10:00:25 UTC (rev 52180)
@@ -0,0 +1,260 @@
+#include "linesegmentsextractor.h"
+
+#include "matrix.h"
+
+#include <stdio.h>
+
+#include <cmath>
+#include <vector>
+#include <list>
+
+using namespace matrix;
+
+struct Point
+{
+    int x;
+    int y;
+};
+
+/* non-member functions */
+
+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(const Segment& segment, const 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;
+}
+
+/* member functions */
+
+void LineSegmentsExtractor::extract(LineCoordinates lineCoordinates,
+                                    float orient,
+                                    SegmentList& segments)
+{
+    const int rows = mImage.rows();
+    const int cols = mImage.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;
+    const int shift = 16;
+
+    std::vector<int> indicesI(lineWidth);
+    std::vector<int> indicesJ(lineWidth);
+
+    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, lineWidth);
+
+                remove_points(lineCoordinates, indicesI, indicesJ);
+
+                // for each non-zero point:
+                //    update line end,
+                //    clear the mask element
+                //    reset the gap
+                if ( isData(mImage, indicesI, indicesJ, cols, rows) )
+                {
+                    if (gap > gapSize)
+                    {
+                        numOfGaps++;
+                        if (numOfGaps > maxNumOfGaps)
+                        {
+                            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());
+}

Copied: grass-addons/grass7/raster/r.ht/linesegmentsextractor.h (from rev 52167, grass-addons/grass7/raster/r.ht/extract_line.h)
===================================================================
--- grass-addons/grass7/raster/r.ht/linesegmentsextractor.h	                        (rev 0)
+++ grass-addons/grass7/raster/r.ht/linesegmentsextractor.h	2012-06-21 10:00:25 UTC (rev 52180)
@@ -0,0 +1,47 @@
+#ifndef LILESEGMENTSEXTRACTOR_H
+#define LILESEGMENTSEXTRACTOR_H
+
+
+#include "houghparameters.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;
+
+class LineSegmentsExtractor
+{
+public:
+    /**
+      \param image should exist during existence of LineSegmentsExtractor
+      \param lineCoordinates will be copied to internal variable
+      */
+    LineSegmentsExtractor(const matrix::Matrix& image,
+                          const ExtractParametres& parametres)
+        : mImage(image),
+          gapSize(parametres.gapSize), maxNumOfGaps(parametres.maxNumOfGaps),
+          lineGap(parametres.maxGap), lineLength(parametres.lineLength),
+          lineWidth(parametres.lineWidth)
+    {}
+
+    /**
+      \param lineCoordinates list of lines
+      \param orient orientation of lines
+      \param[out] segments extracted segments will be added to this list
+      */
+    void extract(LineCoordinates lineCoordinates, const float orient, SegmentList& segments);
+
+private:
+    const matrix::Matrix &mImage;
+    int gapSize;
+    int maxNumOfGaps;
+    int lineGap;
+    int lineLength;
+    int lineWidth;
+};
+
+#endif // LILESEGMENTSEXTRACTOR_H

Modified: grass-addons/grass7/raster/r.ht/main.cpp
===================================================================
--- grass-addons/grass7/raster/r.ht/main.cpp	2012-06-21 09:23:07 UTC (rev 52179)
+++ grass-addons/grass7/raster/r.ht/main.cpp	2012-06-21 10:00:25 UTC (rev 52180)
@@ -52,7 +52,8 @@
     struct Option *input, *output, *anglesOption, *houghImageNameOption,
             *angleWidthOption,
             *minGapOption, *maxNumberOfGapsOption,
-        *maxLinesOption, *maxGapOption, *minSegmentLengthOption;
+        *maxLinesOption, *maxGapOption, *minSegmentLengthOption,
+            *lineWidthOption;
 
     /* initialize GIS environment */
     G_gisinit(argv[0]);		/* reads grass env, stores program name to G_program_name() */
@@ -136,6 +137,14 @@
     minSegmentLengthOption->description = _("Minimal length of line segment");
     minSegmentLengthOption->answer = const_cast<char *>("50");
 
+    lineWidthOption = G_define_option();
+    lineWidthOption->key = "line_width";
+    lineWidthOption->type = TYPE_INTEGER;
+    lineWidthOption->required = NO;
+    lineWidthOption->multiple = NO;
+    lineWidthOption->description = _("Expected width of line (used for searching segments)");
+    lineWidthOption->answer = const_cast<char *>("3");
+
     /* options and flags parser */
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
@@ -144,15 +153,19 @@
     result = output->answer;
     name = input->answer;
 
-    int maxPeaks = atoi(maxLinesOption->answer);
-    int threshold = 10;
-    int angleWidth = atoi(angleWidthOption->answer);
-    int gapSize = atoi(minGapOption->answer);
-    int gap = atoi(maxGapOption->answer);
-    int maxNumOfGaps = atoi(maxNumberOfGapsOption->answer);
-    int minSegmentLength = atoi(minSegmentLengthOption->answer);
-    int sizeOfNeighbourhood = 1;
+    ExtractParametres extractParametres;
+    HoughParametres houghParametres;
 
+    houghParametres.maxPeaksNum = atoi(maxLinesOption->answer);
+    houghParametres.threshold = 10;
+    houghParametres.angleWidth = atoi(angleWidthOption->answer);
+    extractParametres.gapSize = atoi(minGapOption->answer);
+    extractParametres.maxGap = atoi(maxGapOption->answer);
+    extractParametres.maxNumOfGaps = atoi(maxNumberOfGapsOption->answer);
+    extractParametres.lineLength = atoi(minSegmentLengthOption->answer);
+    extractParametres.lineWidth = atoi(lineWidthOption->answer);
+    houghParametres.sizeOfNeighbourhood = 1;
+
     /* returns NULL if the map was not found in any mapset,
      * mapset name otherwise */
     mapset = (char *)G_find_raster2(name, "");
@@ -178,7 +191,7 @@
 
     /* **** */
 
-    hough_peaks(maxPeaks, threshold, angleWidth, sizeOfNeighbourhood, gapSize, maxNumOfGaps, gap, minSegmentLength, name, mapset, nrows, ncols, anglesOption->answer, houghImageNameOption->answer, result);
+    hough_peaks(houghParametres, extractParametres, name, mapset, nrows, ncols, anglesOption->answer, houghImageNameOption->answer, result);
 
     /* **** */
 



More information about the grass-commit mailing list