Added image_gradients as well as find_bright_lines() and find_dark_lines().

dlib/image_transforms/edge_detector.h

@@ -8,6 +8,7 @@
 #include "../array2d.h"
 #include "../geometry.h"
 #include <vector>
+#include "../image_keypoint/build_separable_poly_filters.h"
 
 namespace dlib
 {
@@ -393,6 +394,233 @@ namespace dlib
         return newpixels;
     }
 
+// ----------------------------------------------------------------------------------------
+
+    class image_gradients
+    {
+
+    public:
+
+        image_gradients (
+        ) : image_gradients(1) {}
+
+        image_gradients (
+            long scale
+        ) : the_scale(scale)
+        {
+            DLIB_CASSERT(scale >= 1);
+
+            scale = 2*scale+1;
+
+            auto dfilters = build_separable_poly_filters(2,scale); 
+            DLIB_CASSERT(dfilters[1].size() == 1);
+            DLIB_CASSERT(dfilters[2].size() == 1);
+            DLIB_CASSERT(dfilters[3].size() == 1);
+            DLIB_CASSERT(dfilters[4].size() == 1);
+            DLIB_CASSERT(dfilters[5].size() == 1);
+
+            filter_x.first = matrix_cast<float>(dfilters[1][0].first);
+            filter_x.second = matrix_cast<float>(dfilters[1][0].second);
+
+            filter_y.first = matrix_cast<float>(dfilters[2][0].first);
+            filter_y.second = matrix_cast<float>(dfilters[2][0].second);
+
+            // We multiply by 2 so that the filter gives the gradient rather than the x^2
+            // polynomial coefficient.
+            filter_xx.first = 2*matrix_cast<float>(dfilters[3][0].first);
+            filter_xx.second = matrix_cast<float>(dfilters[3][0].second);
+
+            filter_xy.first = matrix_cast<float>(dfilters[4][0].first);
+            filter_xy.second = matrix_cast<float>(dfilters[4][0].second);
+
+            // We multiply by 2 so that the filter gives the gradient rather than the y^2
+            // polynomial coefficient.
+            filter_yy.first = 2*matrix_cast<float>(dfilters[5][0].first);
+            filter_yy.second = matrix_cast<float>(dfilters[5][0].second);
+        }
+
+
+        long get_scale() const { return the_scale; }
+
+        template <
+            typename in_image_type,
+            typename out_image_type
+            >
+        rectangle gradient_x(
+            const in_image_type& img,
+            out_image_type& out
+        ) const 
+        { 
+            return spatially_filter_image_separable(img, out, filter_x.second, filter_x.first); 
+        }
+
+        template <
+            typename in_image_type,
+            typename out_image_type
+            >
+        rectangle gradient_y(
+            const in_image_type& img,
+            out_image_type& out
+        ) const 
+        { 
+            return spatially_filter_image_separable(img, out, filter_y.second, filter_y.first); 
+        }
+
+
+        template <
+            typename in_image_type,
+            typename out_image_type
+            >
+        rectangle gradient_xx(
+            const in_image_type& img,
+            out_image_type& out
+        )  const
+        { 
+            return spatially_filter_image_separable(img, out, filter_xx.second, filter_xx.first); 
+        }
+
+        template <
+            typename in_image_type,
+            typename out_image_type
+            >
+        rectangle gradient_xy(
+            const in_image_type& img,
+            out_image_type& out
+        ) const 
+        { 
+            return spatially_filter_image_separable(img, out, filter_xy.second, filter_xy.first); 
+        }
+
+        template <
+            typename in_image_type,
+            typename out_image_type
+            >
+        rectangle gradient_yy(
+            const in_image_type& img,
+            out_image_type& out
+        ) const
+        { 
+            return spatially_filter_image_separable(img, out, filter_yy.second, filter_yy.first); 
+        }
+
+        matrix<float> get_x_filter()  const { return filter_x.first*trans(filter_x.second); }
+        matrix<float> get_y_filter()  const { return filter_y.first*trans(filter_y.second); }
+        matrix<float> get_xx_filter() const { return filter_xx.first*trans(filter_xx.second); }
+        matrix<float> get_xy_filter() const { return filter_xy.first*trans(filter_xy.second); }
+        matrix<float> get_yy_filter() const { return filter_yy.first*trans(filter_yy.second); }
+
+    private:
+        std::pair<matrix<float,0,1>,matrix<float,0,1>> filter_x;
+        std::pair<matrix<float,0,1>,matrix<float,0,1>> filter_y;
+        std::pair<matrix<float,0,1>,matrix<float,0,1>> filter_xx;
+        std::pair<matrix<float,0,1>,matrix<float,0,1>> filter_xy;
+        std::pair<matrix<float,0,1>,matrix<float,0,1>> filter_yy;
+
+        long the_scale;
+    };
+
+// ----------------------------------------------------------------------------------------
+
+    namespace impl
+    {
+        template <
+            typename in_image_type,
+            typename out_image_type
+            >
+        void find_lines(
+            const in_image_type& xx_,
+            const in_image_type& xy_,
+            const in_image_type& yy_,
+            out_image_type& horz_,
+            out_image_type& vert_,
+            double positive_if_should_find_dark_lines
+        )
+        {
+            typedef typename image_traits<out_image_type>::pixel_type out_pixel_type;
+            static_assert(std::is_same<float,out_pixel_type>::value || std::is_same<double,out_pixel_type>::value,
+                "Output images must contain either float or double valued pixels");
+
+            const_image_view<in_image_type> xx(xx_);
+            const_image_view<in_image_type> xy(xy_);
+            const_image_view<in_image_type> yy(yy_);
+
+            DLIB_CASSERT(xx.nr() == xy.nr());
+            DLIB_CASSERT(xx.nr() == yy.nr());
+            DLIB_CASSERT(xx.nc() == xy.nc());
+            DLIB_CASSERT(xx.nc() == yy.nc());
+
+
+            image_view<out_image_type> x(horz_);
+            image_view<out_image_type> y(vert_);
+
+            x.set_size(xx.nr(), xx.nc());
+            y.set_size(xx.nr(), xx.nc());
+
+
+            // store the max eigenvalue into xy and then the associated eigen vector into [xx,yy]
+            for (long r = 0; r < xx.nr(); ++r)
+            {
+                for (long c = 0; c < xx.nc(); ++c)
+                {
+                    // negate to that lambda will be the *minimum* eigenvalue
+                    double w1 = positive_if_should_find_dark_lines*xx[r][c]/2.0;
+                    double w2 = positive_if_should_find_dark_lines*yy[r][c]/2.0;
+                    double w3 = positive_if_should_find_dark_lines*xy[r][c];
+
+
+                    auto lambda = w1 + w2 + std::sqrt((w1-w2)*(w1-w2) + w3*w3);
+                    if (lambda < 0)
+                        lambda = 0;
+
+                    if (2*w1!=lambda)
+                    {
+                        x[r][c] = -w3/(2*w1-lambda);
+                        y[r][c] = 1;
+
+                        double norm = std::sqrt(x[r][c]*x[r][c] + y[r][c]*y[r][c]);
+                        x[r][c] *= lambda/norm;
+                        y[r][c] *= lambda/norm;
+                    }
+                    else
+                    {
+                        x[r][c] = lambda;
+                        y[r][c] = 0;
+                    }
+                }
+            }
+        }
+    }
+
+    template <
+        typename in_image_type,
+        typename out_image_type
+        >
+    void find_bright_lines(
+        const in_image_type& xx,
+        const in_image_type& xy,
+        const in_image_type& yy,
+        out_image_type& horz,
+        out_image_type& vert
+    )
+    {
+        impl::find_lines(xx,xy,yy,horz,vert,-1);
+    }
+
+    template <
+        typename in_image_type,
+        typename out_image_type
+        >
+    void find_dark_lines(
+        const in_image_type& xx,
+        const in_image_type& xy,
+        const in_image_type& yy,
+        out_image_type& horz,
+        out_image_type& vert
+    )
+    {
+        impl::find_lines(xx,xy,yy,horz,vert,+1);
+    }
+
 // ----------------------------------------------------------------------------------------
 
 }

dlib/image_transforms/edge_detector_abstract.h

@@ -165,6 +165,285 @@ namespace dlib
                 - PTS is just line with some elements removed.
     !*/
 
+// ----------------------------------------------------------------------------------------
+
+    class image_gradients
+    {
+        /*!
+            WHAT THIS OBJECT REPRESENTS
+                This class is a tool for computing first and second derivatives of an
+                image.  It does this by fitting a quadratic surface around each pixel and
+                then computing the gradients of that quadratic surface.  For the details
+                see the paper:
+                    Quadratic models for curved line detection in SAR CCD by Davis E. King
+                    and Rhonda D. Phillips
+
+
+                This technique gives very accurate gradient estimates and is also very fast
+                since the entire gradient estimation procedure, for each type of gradient,
+                is accomplished with a single separable filter.  This means you can compute
+                gradients at very large scales (e.g. by fitting the quadratic to a large
+                window like a 99x99 window) and it still runs very quickly.
+        !*/
+
+    public:
+
+        image_gradients (
+        );
+        /*!
+            ensures
+                - #get_scale() == 1
+        !*/
+
+        image_gradients (
+            long scale
+        ); 
+        /*!
+            requires
+                - scale >= 1
+            ensures
+                - #get_scale() == scale
+        !*/
+
+        long get_scale(
+        ) const;
+        /*!
+            ensures
+                - When we estimate a gradient we do so by fitting a quadratic filter so a
+                  window of size get_scale()*2+1 centered on each pixel.  Therefore, the
+                  scale parameter controls the size of gradients we will find.  For
+                  example, a very large scale will cause the gradient_xx() to be
+                  insensitive to high frequency noise in the image while smaller scales
+                  would be more sensitive to such fluctuations in the image.
+        !*/
+
+        template <
+            typename in_image_type,
+            typename out_image_type
+            >
+        rectangle gradient_x(
+            const in_image_type& img,
+            out_image_type& out
+        ) const;
+        /*!
+            requires
+                - in_image_type == an image object that implements the interface defined in
+                  dlib/image_processing/generic_image.h 
+                - out_image_type == an image object that implements the interface defined in
+                  dlib/image_processing/generic_image.h 
+                - img and out do not contain pixels with an alpha channel.  That is,
+                  pixel_traits::has_alpha is false for the pixels in these objects.
+                - is_same_object(img, out) == false 
+                - out_image_type must use signed grayscale pixels
+            ensures
+                - Let VALID_AREA = shrink_rect(get_rect(img),get_scale()).
+                - Computes the 1st order gradient of img in the x direction at each
+                  location in VALID_AREA.  The gradients are stored in out.  All pixels in
+                  #out that are outside VALID_AREA are set to 0.
+                - #num_rows(out) == num_rows(img)
+                - #num_columns(out) == num_columns(img)
+                - While not a requirement, it is a good idea if the output image contains
+                  float or double pixels.  If get_scale() is small then this is less of an
+                  issue, but at large scales the gradient can easily be a small number
+                  which is not well represented by integer pixel type such as short. Also,
+                  if you use float pixels in the input and output images then this routine
+                  will use SIMD instructions and is particularly fast.
+                - returns VALID_AREA.  That is, returns the part of the output image which
+                  contains actual valid gradient values.
+        !*/
+
+        template <
+            typename in_image_type,
+            typename out_image_type
+            >
+        rectangle gradient_y(
+            const in_image_type& img,
+            out_image_type& out
+        ) const;
+        /*!
+            This routine is identical to gradient_x() (defined above) except that it
+            computes the 1st order y gradient.
+        !*/
+
+        template <
+            typename in_image_type,
+            typename out_image_type
+            >
+        rectangle gradient_xx(
+            const in_image_type& img,
+            out_image_type& out
+        )  const;
+        /*!
+            This routine is identical to gradient_x() (defined above) except that it
+            computes the 2nd order x gradient. 
+        !*/
+
+        template <
+            typename in_image_type,
+            typename out_image_type
+            >
+        rectangle gradient_xy(
+            const in_image_type& img,
+            out_image_type& out
+        ) const;
+        /*!
+            This routine is identical to gradient_x() (defined above) except that it
+            computes the partial derivative with respect to x and y.
+        !*/
+
+        template <
+            typename in_image_type,
+            typename out_image_type
+            >
+        rectangle gradient_yy(
+            const in_image_type& img,
+            out_image_type& out
+        ) const;
+        /*!
+            This routine is identical to gradient_x() (defined above) except that it
+            computes the 2nd order y gradient. 
+        !*/
+
+        matrix<float> get_x_filter(
+        )  const; 
+        /*!
+            ensures
+                - Returns the filter used by gradient_x() to compute the image gradient.
+                  That is, the output of gradient_x() is found by cross correlating the
+                  filter get_x_filter() with the image.
+                - The returned filter has get_scale()*2+1 rows and columns.
+        !*/
+
+        matrix<float> get_y_filter(
+        )  const; 
+        /*!
+            ensures
+                - Returns the filter used by gradient_y() to compute the image gradient.
+                  That is, the output of gradient_y() is found by cross correlating the
+                  filter get_y_filter() with the image.
+                - The returned filter has get_scale()*2+1 rows and columns.
+        !*/
+
+        matrix<float> get_xx_filter(
+        ) const; 
+        /*!
+            ensures
+                - Returns the filter used by gradient_xx() to compute the image gradient.
+                  That is, the output of gradient_xx() is found by cross correlating the
+                  filter get_xx_filter() with the image.
+                - The returned filter has get_scale()*2+1 rows and columns.
+        !*/
+
+
+        matrix<float> get_xy_filter(
+        ) const; 
+        /*!
+            ensures
+                - Returns the filter used by gradient_xy() to compute the image gradient.
+                  That is, the output of gradient_xy() is found by cross correlating the
+                  filter get_xy_filter() with the image.
+                - The returned filter has get_scale()*2+1 rows and columns.
+        !*/
+
+
+        matrix<float> get_yy_filter(
+        ) const; 
+        /*!
+            ensures
+                - Returns the filter used by gradient_yy() to compute the image gradient.
+                  That is, the output of gradient_yy() is found by cross correlating the
+                  filter get_yy_filter() with the image.
+                - The returned filter has get_scale()*2+1 rows and columns.
+        !*/
+
+    };
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename in_image_type,
+        typename out_image_type
+        >
+    void find_bright_lines(
+        const in_image_type& xx,
+        const in_image_type& xy,
+        const in_image_type& yy,
+        out_image_type& horz,
+        out_image_type& vert
+    );
+    /*!
+        requires
+            - in_image_type == an image object that implements the interface defined in
+              dlib/image_processing/generic_image.h 
+            - out_image_type == an image object that implements the interface defined in
+              dlib/image_processing/generic_image.h 
+            - All images are grayscale and the horz and vert images must contain float or
+              double pixel types.
+            - num_rows(xx) == num_rows(xy) == num_rows(yy)
+            - num_columns(xx) == num_columns(xy) == num_columns(yy)
+        ensures
+            - This routine is similar to sobel_edge_detector(), except instead of finding
+              an edge it finds a bright/white line.  For example, the border between a
+              black piece of paper and a white table is an edge, but a curve drawn with a
+              pencil on a piece of paper makes a line.  Therefore, the output of this
+              routine is a vector field encoded in the horz and vert images.  The vector
+              obtains a large magnitude when centered on a bright line in an image and the
+              direction of the vector is perpendicular to the line.  To be very precise,
+              each vector points in the direction of greatest change in second derivative
+              and the magnitude of the vector encodes the derivative magnitude in that
+              direction.  Moreover, if the second derivative is positive then the output
+              vector is zero.  This zeroing if positive gradients causes the output to be
+              sensitive only to bright lines surrounded by darker pixels.
+            - We assume that xx, xy, and yy are the 3 different second order gradients of
+              the image in question.  You can obtain these gradients using the
+              image_gradients class.
+            - The output images will have the same sizes as the input images, that is:
+                - #num_rows(horz) == #num_rows(vert) == num_rows(xx)
+                - #num_columns(horz) == #num_columns(vert) == num_columns(xx)
+    !*/
+
+    template <
+        typename in_image_type,
+        typename out_image_type
+        >
+    void find_dark_lines(
+        const in_image_type& xx,
+        const in_image_type& xy,
+        const in_image_type& yy,
+        out_image_type& horz,
+        out_image_type& vert
+    );
+    /*!
+        requires
+            - in_image_type == an image object that implements the interface defined in
+              dlib/image_processing/generic_image.h 
+            - out_image_type == an image object that implements the interface defined in
+              dlib/image_processing/generic_image.h 
+            - All images are grayscale and the horz and vert images must contain float or
+              double pixel types.
+            - num_rows(xx) == num_rows(xy) == num_rows(yy)
+            - num_columns(xx) == num_columns(xy) == num_columns(yy)
+        ensures
+            - This routine is similar to sobel_edge_detector(), except instead of finding
+              an edge it finds a dark/black line.  For example, the border between a
+              black piece of paper and a white table is an edge, but a curve drawn with a
+              pencil on a piece of paper makes a line.  Therefore, the output of this
+              routine is a vector field encoded in the horz and vert images.  The vector
+              obtains a large magnitude when centered on a dark line in an image and the
+              direction of the vector is perpendicular to the line.  To be very precise,
+              each vector points in the direction of greatest change in second derivative
+              and the magnitude of the vector encodes the derivative magnitude in that
+              direction.  Moreover, if the second derivative is negative then the output
+              vector is zero.  This zeroing if negative gradients causes the output to be
+              sensitive only to dark lines surrounded by light pixels.
+            - We assume that xx, xy, and yy are the 3 different second order gradients of
+              the image in question.  You can obtain these gradients using the
+              image_gradients class.
+            - The output images will have the same sizes as the input images, that is:
+                - #num_rows(horz) == #num_rows(vert) == num_rows(xx)
+                - #num_columns(horz) == #num_columns(vert) == num_columns(xx)
+    !*/
+
 // ----------------------------------------------------------------------------------------
 
 }