Added hough_transform

dlib/image_transforms.h

@@ -12,6 +12,7 @@
 #include "image_transforms/draw.h"
 #include "image_transforms/integral_image.h"
 #include "image_transforms/image_pyramid.h"
+#include "image_transforms/hough_transform.h"
 #include "image_transforms/label_connected_blobs.h"
 #include "image_transforms/colormaps.h"
 #include "image_transforms/segment_image.h"

dlib/image_transforms/hough_transform.h

+// Copyright (C) 2014  Davis E. King (davis@dlib.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#ifndef DLIB_HOUGH_tRANSFORM_Hh_
+#define DLIB_HOUGH_tRANSFORM_Hh_
+
+#include "hough_transform_abstract.h"
+#include "../image_processing/generic_image.h"
+#include "../geometry.h"
+#include "../algs.h"
+#include "assign_image.h"
+
+namespace dlib
+{
+
+// ----------------------------------------------------------------------------------------
+
+    class hough_transform
+    {
+
+    public:
+        explicit hough_transform (
+            unsigned long size_
+        ) : _size(size_) 
+        {
+            DLIB_CASSERT(size_ > 0,
+                "\t hough_transform::hough_transform(size_)"
+                << "\n\t Invalid arguments given to this function."
+                );
+
+            even_size = _size - (_size%2);
+
+            const point cent = center(rectangle(0,0,size_-1,size_-1));
+            xcos_theta.set_size(size_, size_);
+            ysin_theta.set_size(size_, size_);
+
+            std::vector<double> cos_theta(size_), sin_theta(size_);
+            const double scale = 1<<16;
+            for (unsigned long t = 0; t < size_; ++t)
+            {
+                double theta = t*pi/even_size;
+
+                cos_theta[t] = scale*std::cos(theta)/sqrt_2;
+                sin_theta[t] = scale*std::sin(theta)/sqrt_2;
+            }
+            const double offset = scale*even_size/4.0 + 0.5;
+
+            for (unsigned long c = 0; c < size_; ++c)
+            {
+                const long x = c - cent.x();
+                for (unsigned long t = 0; t < size_; ++t)
+                    xcos_theta(c,t) = static_cast<int32>(x*cos_theta[t] + offset);
+            }
+            for (unsigned long r = 0; r < size_; ++r)
+            {
+                const long y = r - cent.y();
+                for (unsigned long t = 0; t < size_; ++t)
+                    ysin_theta(r,t) = static_cast<int32>(y*sin_theta[t] + offset);
+            }
+        }
+
+        unsigned long size(
+        ) const { return _size; }
+
+        long nr(
+        ) const { return _size; }
+
+        long nc(
+        ) const { return _size; }
+
+        std::pair<point, point> get_line (
+            const point& p
+        ) const
+        {
+            DLIB_ASSERT(rectangle(0,0,size()-1,size()-1).contains(p) == true,
+                "\t pair<point,point> hough_transform::get_line(point)"
+                << "\n\t Invalid arguments given to this function."
+                << "\n\t p:      " << p 
+                << "\n\t size(): " << size()
+                );
+
+            // First we compute the radius measured in pixels from the center and the theta
+            // angle in radians.
+            typedef dlib::vector<double,2> vect;
+            const rectangle box(0,0,size()-1,size()-1);
+            const vect cent = center(box);
+            double theta = p.x()-cent.x();
+            double radius = p.y()-cent.y();
+            theta = theta*pi/even_size;
+            radius = radius*sqrt_2 + 0.5;
+
+            // now make a line segment on the line.
+            vect v1 = cent + vect(size()+1000,0) + vect(0,radius);
+            vect v2 = cent - vect(size()+1000,0) + vect(0,radius);
+            point p1 = rotate_point(cent, v1, theta);
+            point p2 = rotate_point(cent, v2, theta);
+
+            clip_line_to_rectangle(box, p1, p2);
+
+            return std::make_pair(p1,p2);
+        }
+
+        template <
+            typename in_image_type,
+            typename out_image_type
+            >
+        void operator() (
+            const in_image_type& img_,
+            const rectangle& box,
+            out_image_type& himg_
+        ) const
+        {
+            typedef typename image_traits<in_image_type>::pixel_type in_pixel_type;
+            typedef typename image_traits<out_image_type>::pixel_type out_pixel_type;
+
+            DLIB_CASSERT(box.width() == size() && box.height() == size(),
+                "\t hough_transform::hough_transform(size_)"
+                << "\n\t Invalid arguments given to this function."
+                << "\n\t box.width():  " << box.width()
+                << "\n\t box.height(): " << box.height()
+                << "\n\t size():       " << size()
+                );
+
+            COMPILE_TIME_ASSERT(pixel_traits<in_pixel_type>::grayscale == true);
+            COMPILE_TIME_ASSERT(pixel_traits<out_pixel_type>::grayscale == true);
+
+            const_image_view<in_image_type> img(img_);
+            image_view<out_image_type> himg(himg_);
+
+            himg.set_size(size(), size());
+            assign_all_pixels(himg, 0);
+
+            const rectangle area = box.intersect(get_rect(img));
+
+            const long max_n8 = (himg.nc()/8)*8;
+            const long max_n4 = (himg.nc()/4)*4;
+            for (long r = area.top(); r <= area.bottom(); ++r)
+            {
+                const int32* ysin_base = &ysin_theta(r-box.top(),0);
+                for (long c = area.left(); c <= area.right(); ++c)
+                {
+                    const out_pixel_type val = static_cast<out_pixel_type>(img[r][c]);
+                    if (val != 0)
+                    {
+                        /*
+                        // The code in this comment is equivalent to the more complex but
+                        // faster code below.  We keep this simple version of the Hough
+                        // transform implementation here just to document what it's doing
+                        // more clearly.
+                        const point cent = center(box);
+                        const long x = c - cent.x();
+                        const long y = r - cent.y();
+                        for (long t = 0; t < himg.nc(); ++t)
+                        {
+                            double theta = t*pi/even_size;
+                            double radius = (x*std::cos(theta) + y*std::sin(theta))/sqrt_2 + even_size/2 + 0.5;
+                            long rr = static_cast<long>(radius);
+                            himg[rr][t] += val;
+                        }
+                        continue;
+                        */
+
+                        // Run the speed optimized version of the code in the above
+                        // comment.
+                        const int32* ysin = ysin_base;
+                        const int32* xcos = &xcos_theta(c-box.left(),0);
+                        long t = 0;
+                        while(t < max_n8)
+                        {
+                            long rr0 = (*xcos++ + *ysin++)>>16;
+                            long rr1 = (*xcos++ + *ysin++)>>16;
+                            long rr2 = (*xcos++ + *ysin++)>>16;
+                            long rr3 = (*xcos++ + *ysin++)>>16;
+                            long rr4 = (*xcos++ + *ysin++)>>16;
+                            long rr5 = (*xcos++ + *ysin++)>>16;
+                            long rr6 = (*xcos++ + *ysin++)>>16;
+                            long rr7 = (*xcos++ + *ysin++)>>16;
+
+                            himg[rr0][t++] += val;
+                            himg[rr1][t++] += val;
+                            himg[rr2][t++] += val;
+                            himg[rr3][t++] += val;
+                            himg[rr4][t++] += val;
+                            himg[rr5][t++] += val;
+                            himg[rr6][t++] += val;
+                            himg[rr7][t++] += val;
+                        }
+                        while(t < max_n4)
+                        {
+                            long rr0 = (*xcos++ + *ysin++)>>16;
+                            long rr1 = (*xcos++ + *ysin++)>>16;
+                            long rr2 = (*xcos++ + *ysin++)>>16;
+                            long rr3 = (*xcos++ + *ysin++)>>16;
+                            himg[rr0][t++] += val;
+                            himg[rr1][t++] += val;
+                            himg[rr2][t++] += val;
+                            himg[rr3][t++] += val;
+                        }
+                        while(t < himg.nc())
+                        {
+                            long rr0 = (*xcos++ + *ysin++)>>16;
+                            himg[rr0][t++] += val;
+                        }
+                    }
+                }
+            }
+        }
+
+    private:
+
+        unsigned long _size;
+        unsigned long even_size; // equal to _size if _size is even, otherwise equal to _size-1.
+        matrix<int32> xcos_theta, ysin_theta;
+    };
+}
+
+#endif // DLIB_HOUGH_tRANSFORM_Hh_
+

dlib/image_transforms/hough_transform_abstract.h

+// Copyright (C) 2014  Davis E. King (davis@dlib.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#undef DLIB_HOUGH_tRANSFORM_ABSTRACT_Hh_
+#ifdef DLIB_HOUGH_tRANSFORM_ABSTRACT_Hh_
+
+#include "../geometry.h"
+#include "../image_processing/generic_image.h"
+
+namespace dlib
+{
+
+// ----------------------------------------------------------------------------------------
+
+    class hough_transform
+    {
+        /*!
+            WHAT THIS OBJECT REPRESENTS
+                This object is a tool for computing the line finding version of the Hough
+                transform given some kind of edge detection image as input.  It also allows
+                the edge pixels to be weighted such that higher weighted edge pixels
+                contribute correspondingly more to the output of the Hough transform,
+                allowing stronger edges to create correspondingly stronger line detections
+                in the final Hough transform.
+        !*/
+
+    public:
+
+        explicit hough_transform (
+            unsigned long size_
+        ); 
+        /*!
+            requires
+                - size_ > 0
+            ensures
+                - This object will compute Hough transforms that are size_ by size_ pixels.  
+                  This is in terms of both the Hough accumulator array size as well as the
+                  input image size.
+                - #size() == size_
+        !*/
+
+        unsigned long size(
+        ) const;
+        /*!
+            ensures
+                - returns the size of the Hough transforms generated by this object.  In
+                  particular, this object creates Hough transform images that are size() by
+                  size() pixels in size.
+        !*/
+
+        long nr(
+        ) const;
+        /*!
+            ensures
+                - returns size()
+        !*/
+
+        long nc(
+        ) const;
+        /*!
+            ensures
+                - returns size()
+        !*/
+
+        std::pair<point, point> get_line (
+            const point& p
+        ) const;
+        /*!
+            requires
+                - rectangle(0,0,size()-1,size()-1).contains(p) == true
+                  (i.e. p must be a point inside the Hough accumulator array)
+            ensures
+                - returns the line segment in the original image space corresponding
+                  to Hough transform point p. 
+                - The returned points are inside rectangle(0,0,size()-1,size()-1).
+        !*/
+
+        template <
+            typename in_image_type,
+            typename out_image_type
+            >
+        void operator() (
+            const in_image_type& img,
+            const rectangle& box,
+            out_image_type& himg
+        ) const;
+        /*!
+            requires
+                - in_image_type == an image object that implements the interface defined in
+                  dlib/image_processing/generic_image.h and it must contain grayscale pixels.
+                - out_image_type == an image object that implements the interface defined in
+                  dlib/image_processing/generic_image.h and it must contain grayscale pixels.
+                - box.width() == size()
+                - box.height() == size()
+            ensures
+                - Computes the Hough transform of the part of img contained within box.
+                  In particular, we do a grayscale version of the Hough transform where any
+                  non-zero pixel in img is treated as a potential component of a line and
+                  accumulated into the Hough accumulator #himg.  However, rather than
+                  adding 1 to each relevant accumulator bin we add the value of the pixel
+                  in img to each Hough accumulator bin.  This means that, if all the
+                  pixels in img are 0 or 1 then this routine performs a normal Hough
+                  transform.  However, if some pixels have larger values then they will be
+                  weighted correspondingly more in the resulting Hough transform.
+                - #himg.nr() == size()
+                - #himg.nc() == size()
+                - #himg is the Hough transform of the part of img contained in box.  Each
+                  point in #himg corresponds to a line in the input box.  In particular,
+                  the line for #himg[y][x] is given by get_line(point(x,y)).  Also, when
+                  viewing the #himg image, the x-axis gives the angle of the line and the
+                  y-axis the distance of the line from the center of the box.
+        !*/
+
+    };
+}
+
+#endif // DLIB_HOUGH_tRANSFORM_ABSTRACT_Hh_
+
+