Added an initial cut of the SURF and image keypoint finding code

--HG--
extra : convert_revision : svn%3Afdd8eb12-d10e-0410-9acb-85c331704f74/trunk%402989

dlib/image_keypoint.h

+// Copyright (C) 2009  Davis E. King (davisking@users.sourceforge.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#ifndef DLIB_IMAGE_KEYPOINt_H_
+#define DLIB_IMAGE_KEYPOINt_H_ 
+
+#include "image_keypoint/surf.h"
+#include "image_keypoint/hessian_pyramid.h"
+
+#endif // DLIB_IMAGE_KEYPOINt_H_
+

dlib/image_keypoint/hessian_pyramid.h

+// Copyright (C) 2009  Davis E. King (davisking@users.sourceforge.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#ifndef DLIB_HESSIAN_PYRAMId_H__
+#define DLIB_HESSIAN_PYRAMId_H__
+
+#include "hessian_pyramid_abstract.h"
+#include "../algs.h"
+#include "../image_transforms/integral_image.h"
+#include "../array.h"
+#include "../array2d.h"
+#include "../noncopyable.h"
+#include "../matrix.h"
+#include "../stl_checked.h"
+#include <algorithm>
+#include <vector>
+
+namespace dlib
+{
+
+// ----------------------------------------------------------------------------------------
+
+    struct interest_point
+    {
+        interest_point() : scale(0), score(0), laplacian(0) {}
+
+        dlib::vector<double,2> center;
+        double scale;
+        double score;
+        double laplacian;
+
+        bool operator < (const interest_point& p) const { return score < p.score; }
+    };
+
+// ----------------------------------------------------------------------------------------
+
+    class hessian_pyramid : noncopyable
+    {
+    public:
+        hessian_pyramid()
+        {
+            num_octaves = 0;
+            num_intervals = 0;
+            initial_step_size = 0;
+        }
+
+        template <typename integral_image_type>
+        void build_pyramid (
+            const integral_image_type& img,
+            long num_octaves,
+            long num_intervals,
+            long initial_step_size
+        )
+        {
+            DLIB_ASSERT(num_octaves > 0 && num_intervals > 0 && initial_step_size > 0,
+                "\tvoid build_pyramid()"
+                << "\n\tAll arguments to this function must be > 0"
+                << "\n\t this:              " << this
+                << "\n\t num_octaves:       " << num_octaves 
+                << "\n\t num_intervals:     " << num_intervals 
+                << "\n\t initial_step_size: " << initial_step_size 
+            );
+
+            this->num_octaves = num_octaves;
+            this->num_intervals = num_intervals;
+            this->initial_step_size = initial_step_size;
+
+            // allocate space for the pyramid
+            pyramid.resize(num_octaves*num_intervals);
+            for (long o = 0; o < num_octaves; ++o)
+            {
+                const long step_size = get_step_size(o);
+                for (long i = 0; i < num_intervals; ++i)
+                {
+                    pyramid[num_intervals*o + i].set_size(img.nr()/step_size, img.nc()/step_size);
+                }
+            }
+
+            // now fill out the pyramid with data
+            for (long o = 0; o < num_octaves; ++o)
+            {
+                const long step_size = get_step_size(o);
+                const long border_size = get_border_size(o)*step_size;
+
+                for (long i = 0; i < num_intervals; ++i)
+                {
+                    const long lobe_size = static_cast<long>(std::pow(2.0, o+1.0)+0.5)*(i+1) + 1;
+                    const double area_inv = 1.0/std::pow(3.0*lobe_size, 2.0);
+
+                    const long lobe_offset = lobe_size/2+1;
+                    const point tl(-lobe_offset,-lobe_offset);
+                    const point tr(lobe_offset,-lobe_offset);
+                    const point bl(-lobe_offset,lobe_offset);
+                    const point br(lobe_offset,lobe_offset);
+
+                    for (long r = border_size; r < img.nr() - border_size; r += step_size)
+                    {
+                        for (long c = border_size; c < img.nc() - border_size; c += step_size)
+                        {
+                            const point p(c,r);
+
+                            double Dxx = img.get_sum_of_area(centered_rect(p, lobe_size*3, 2*lobe_size-1)) - 
+                                         img.get_sum_of_area(centered_rect(p, lobe_size,   2*lobe_size-1))*3.0;
+
+                            double Dyy = img.get_sum_of_area(centered_rect(p, 2*lobe_size-1, lobe_size*3)) - 
+                                         img.get_sum_of_area(centered_rect(p, 2*lobe_size-1, lobe_size))*3.0;
+
+                            double Dxy = img.get_sum_of_area(centered_rect(p+bl, lobe_size, lobe_size)) + 
+                                         img.get_sum_of_area(centered_rect(p+tr, lobe_size, lobe_size)) -
+                                         img.get_sum_of_area(centered_rect(p+tl, lobe_size, lobe_size)) -
+                                         img.get_sum_of_area(centered_rect(p+br, lobe_size, lobe_size));
+
+                            //std::cout << "(((((: " << Dxx << "  " << Dyy << "  " << Dxy << std::endl;
+                            // now we normalize the filter responses
+                            Dxx *= area_inv;
+                            Dyy *= area_inv;
+                            Dxy *= area_inv;
+
+
+                            double sign_of_laplacian = +1;
+                            if (Dxx + Dyy < 0)
+                                sign_of_laplacian = -1;
+
+                            double determinant = Dxx*Dyy - 0.81*Dxy*Dxy;
+
+                            // If the determinant is negative then just blank it out by setting
+                            // it to zero.
+                            if (determinant < 0)
+                                determinant = 0;
+
+                            //std::cout << "****: " << Dxx << "  " << Dyy << "  " << Dxy << std::endl;
+                            // Save the determinant of the Hessian into our image pyramid.  Also
+                            // pack the laplacian sign into the value so we can get it out later.
+                            pyramid[o*num_intervals + i][r/step_size][c/step_size] = sign_of_laplacian*determinant;
+
+                        }
+                    }
+
+                }
+            }
+        }
+
+        long get_border_size (
+            long octave
+        ) const
+        {
+            DLIB_ASSERT(0 <= octave && octave < octaves(),
+                "\tlong get_border_size(octave)"
+                << "\n\tInvalid octave value"
+                << "\n\t this:   " << this
+                << "\n\t octave: " << octave 
+            );
+
+            const double lobe_size = std::pow(2.0, octave+1.0)*(num_intervals+1) + 1;
+            const double filter_size = 3*lobe_size;
+
+            const long bs = static_cast<long>(std::ceil(filter_size/2.0));
+            return bs;
+        }
+
+        long get_step_size (
+            long octave
+        ) const
+        {
+            DLIB_ASSERT(0 <= octave && octave < octaves(),
+                "\tlong get_step_size(octave)"
+                << "\n\tInvalid octave value"
+                << "\n\t this:   " << this
+                << "\n\t octave: " << octave 
+            );
+
+            return initial_step_size*static_cast<long>(std::pow(2.0, (double)octave)+0.5);
+        }
+
+        long nr (
+            long octave
+        ) const
+        {
+            DLIB_ASSERT(0 <= octave && octave < octaves(),
+                "\tlong nr(octave)"
+                << "\n\tInvalid octave value"
+                << "\n\t this:   " << this
+                << "\n\t octave: " << octave 
+            );
+
+            return pyramid[num_intervals*octave].nr();
+        }
+
+        long nc (
+            long octave
+        ) const
+        {
+            DLIB_ASSERT(0 <= octave && octave < octaves(),
+                "\tlong nc(octave)"
+                << "\n\tInvalid octave value"
+                << "\n\t this:   " << this
+                << "\n\t octave: " << octave 
+            );
+
+            return pyramid[num_intervals*octave].nc();
+        }
+
+        double get_value (
+            long octave,
+            long interval,
+            long r,
+            long c
+        ) const
+        {
+            DLIB_ASSERT(0 <= octave && octave < octaves() &&
+                        0 <= interval && interval < intervals() &&
+                        get_border_size(octave) <= r && r < nr(octave)-get_border_size(octave) &&
+                        get_border_size(octave) <= c && c < nc(octave)-get_border_size(octave),
+                "\tdouble get_value(octave, interval, r, c)"
+                << "\n\tInvalid inputs to this function"
+                << "\n\t this:      " << this
+                << "\n\t octave:    " << octave 
+                << "\n\t interval:  " << interval 
+                << "\n\t octaves:   " << octaves() 
+                << "\n\t intervals: " << intervals()
+                << "\n\t r:         " << r  
+                << "\n\t c:         " << c 
+                << "\n\t nr(octave): " << nr(octave)  
+                << "\n\t nc(octave): " << nc(octave) 
+                << "\n\t get_border_size(octave): " << get_border_size(octave) 
+            );
+
+            return std::abs(pyramid[num_intervals*octave + interval][r][c]);
+        }
+
+        double get_laplacian (
+            long octave,
+            long interval,
+            long r,
+            long c
+        ) const
+        {
+            DLIB_ASSERT(0 <= octave && octave < octaves() &&
+                        0 <= interval && interval < intervals() &&
+                        get_border_size(octave) <= r && r < nr(octave)-get_border_size(octave) &&
+                        get_border_size(octave) <= c && c < nc(octave)-get_border_size(octave),
+                "\tdouble get_laplacian(octave, interval, r, c)"
+                << "\n\tInvalid inputs to this function"
+                << "\n\t this:      " << this
+                << "\n\t octave:    " << octave 
+                << "\n\t interval:  " << interval 
+                << "\n\t octaves:   " << octaves() 
+                << "\n\t intervals: " << intervals()
+                << "\n\t r:         " << r  
+                << "\n\t c:         " << c 
+                << "\n\t nr(octave): " << nr(octave)  
+                << "\n\t nc(octave): " << nc(octave) 
+                << "\n\t get_border_size(octave): " << get_border_size(octave) 
+            );
+
+            // return the sign of the laplacian
+            if (pyramid[num_intervals*octave + interval][r][c] > 0)
+                return +1;
+            else
+                return -1;
+        }
+
+        long octaves (
+        ) const { return num_octaves; }
+
+        long intervals (
+        ) const { return num_intervals; }
+
+    private:
+
+        long num_octaves;
+        long num_intervals;
+        long initial_step_size;
+
+        typedef array2d<double>::kernel_1a image_type;
+        typedef array<image_type>::expand_1d pyramid_type;
+
+        pyramid_type pyramid;
+    };
+
+// ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+// ----------------------------------------------------------------------------------------
+
+    namespace hessian_pyramid_helpers
+    {
+        inline bool is_maximum_in_region(
+            const hessian_pyramid& pyr,
+            long o, 
+            long i, 
+            long r, 
+            long c
+        )
+        {
+            // First check if this point is near the edge of the octave 
+            // If it is then we say it isn't a maximum as these points are
+            // not as reliable.
+            if (i <= 0 || i+1 >= pyr.intervals())
+            {
+                return false;
+            }
+
+            const double val = pyr.get_value(o,i,r,c);
+
+            // now check if there are any bigger values around this guy
+            for (long ii = i-1; ii <= i+1; ++ii)
+            {
+                for (long rr = r-1; rr <= r+1; ++rr)
+                {
+                    for (long cc = c-1; cc <= c+1; ++cc)
+                    {
+                        if (pyr.get_value(o,ii,rr,cc) > val)
+                            return false;
+                    }
+                }
+            }
+
+            return true;
+        }
+
+    // ------------------------------------------------------------------------------------
+
+        inline const matrix<double,3,1> get_hessian_gradient (
+            const hessian_pyramid& pyr,
+            long o, 
+            long i, 
+            long r, 
+            long c
+        )
+        {
+            matrix<double,3,1> grad;
+            grad(0) = (pyr.get_value(o,i,r,c+1) - pyr.get_value(o,i,r,c-1))/2.0;
+            grad(1) = (pyr.get_value(o,i,r+1,c) - pyr.get_value(o,i,r-1,c))/2.0;
+            grad(2) = (pyr.get_value(o,i+1,r,c) - pyr.get_value(o,i-1,r,c))/2.0;
+            return grad;
+        }
+
+    // ------------------------------------------------------------------------------------
+
+        inline const matrix<double,3,3> get_hessian_hessian (
+            const hessian_pyramid& pyr,
+            long o, 
+            long i, 
+            long r, 
+            long c
+        )
+        {
+            matrix<double,3,3> hess;
+            const double val = pyr.get_value(o,i,r,c);
+
+            double Dxx = (pyr.get_value(o,i,r,c+1) + pyr.get_value(o,i,r,c-1)) - 2*val;
+            double Dyy = (pyr.get_value(o,i,r+1,c) + pyr.get_value(o,i,r-1,c)) - 2*val;
+            double Dss = (pyr.get_value(o,i+1,r,c) + pyr.get_value(o,i-1,r,c)) - 2*val;
+
+            double Dxy = (pyr.get_value(o,i,r+1,c+1) + pyr.get_value(o,i,r-1,c-1) -
+                          pyr.get_value(o,i,r-1,c+1) - pyr.get_value(o,i,r+1,c-1)) / 4.0;
+
+            double Dxs = (pyr.get_value(o,i+1,r,c+1) + pyr.get_value(o,i-1,r,c-1) -
+                          pyr.get_value(o,i-1,r,c+1) - pyr.get_value(o,i+1,r,c-1)) / 4.0;
+
+            double Dys = (pyr.get_value(o,i+1,r+1,c) + pyr.get_value(o,i-1,r-1,c) -
+                          pyr.get_value(o,i-1,r+1,c) - pyr.get_value(o,i+1,r-1,c)) / 4.0;
+
+
+            hess = Dxx, Dxy, Dxs,
+            Dxy, Dyy, Dys,
+            Dxs, Dys, Dss;
+
+            return hess;
+        }
+
+    // ------------------------------------------------------------------------------------
+
+        inline const interest_point interpolate_point (
+            const hessian_pyramid& pyr, 
+            long o, 
+            long i, 
+            long r, 
+            long c
+        )
+        {
+            dlib::vector<double,2> p(c,r);
+
+            dlib::vector<double,3> start_point(c,r,i);
+            dlib::vector<double,3> interpolated_point = -inv(get_hessian_hessian(pyr,o,i,r,c))*get_hessian_gradient(pyr,o,i,r,c);
+
+            //cout << "inter: " <<  trans(interpolated_point);
+
+            interest_point temp;
+            if (max(abs(interpolated_point)) < 0.5)
+            {
+                p = (start_point+interpolated_point)*pyr.get_step_size(o);
+                const double lobe_size = std::pow(2.0, o+1.0)*(i+interpolated_point.z()+1) + 1;
+                const double filter_size = 3*lobe_size;
+                const double scale = 1.2/9.0 * filter_size;
+
+                temp.center = p;
+                temp.scale = scale;
+                temp.score = pyr.get_value(o,i,r,c);
+                temp.laplacian = pyr.get_laplacian(o,i,r,c);
+            }
+
+            return temp;
+        }
+
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename Alloc>
+    void get_interest_points (
+        const hessian_pyramid& pyr,
+        double threshold,
+        std::vector<interest_point,Alloc>& result_points
+    )
+    {
+        using namespace std;
+        using namespace hessian_pyramid_helpers;
+
+        result_points.clear();
+
+        for (long o = 0; o < pyr.octaves(); ++o)
+        {
+            const long border_size = pyr.get_border_size(o);
+            const long nr = pyr.nr(o);
+            const long nc = pyr.nc(o);
+
+            // do non-maximum suppression on all the intervals in the current octave and 
+            // accumulate the results in result_points
+            for (long i = 1; i < pyr.intervals()-1;  i += 3)
+            {
+                for (long r = border_size+1; r < nr - border_size-1; r += 3)
+                {
+                    for (long c = border_size+1; c < nc - border_size-1; c += 3)
+                    {
+                        double max_val = pyr.get_value(o,i,r,c);
+                        long max_i = i;
+                        long max_r = r;
+                        long max_c = c;
+
+                        // loop over this 3x3x3 block and find the largest element
+                        for (long ii = i; ii < std::min(i + 3, pyr.intervals()-1); ++ii)
+                        {
+                            for (long rr = r; rr < std::min(r + 3, nr - border_size - 1); ++rr)
+                            {
+                                for (long cc = c; cc < std::min(c + 3, nc - border_size - 1); ++cc)
+                                {
+                                    double temp = pyr.get_value(o,ii,rr,cc);
+                                    if (temp > max_val)
+                                    {
+                                        max_val = temp;
+                                        max_i = ii;
+                                        max_r = rr;
+                                        max_c = cc;
+                                    }
+                                }
+                            }
+                        }
+
+                        // If the max point we found is really a maximum in its own region and
+                        // is big enough then add it to the results.
+                        if (max_val > threshold && is_maximum_in_region(pyr, o, max_i, max_r, max_c))
+                        {
+                            //cout << max_val << endl;
+                            interest_point sp = interpolate_point (pyr, o, max_i, max_r, max_c);
+                            if (sp.score > threshold)
+                            {
+                                result_points.push_back(sp);
+                            }
+                        }
+
+                    }
+                }
+            }
+        }
+
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename Alloc>
+    void get_interest_points (
+        const hessian_pyramid& pyr,
+        double threshold,
+        std_vector_c<interest_point,Alloc>& result_points
+    )
+    /*!
+        This function is just an overload that automatically casts std_vector_c objects
+        into std::vector objects.  (Usually this is automatic but the template argument
+        there messes up the conversion so we have to do it explicitly)
+    !*/
+    {
+        std::vector<interest_point,Alloc>& v = result_points;
+        get_interest_points(pyr, threshold, v);
+    }
+
+// ----------------------------------------------------------------------------------------
+
+}
+
+#endif  // DLIB_HESSIAN_PYRAMID_H__
+

dlib/image_keypoint/hessian_pyramid_abstract.h

+// Copyright (C) 2009  Davis E. King (davisking@users.sourceforge.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#undef DLIB_HESSIAN_PYRAMId_ABSTRACT_H__
+#ifdef DLIB_HESSIAN_PYRAMId_ABSTRACT_H__
+
+#include "../image_transforms/integral_image_abstract.h"
+#include "../noncopyable.h"
+#include <vector>
+
+namespace dlib
+{
+
+    class hessian_pyramid : noncopyable
+    {
+        /*!
+            INITIAL VALUE
+                - octaves() == 0
+                - intervals() == 0
+
+            WHAT THIS OBJECT REPRESENTS
+                This object represents an image pyramid where each level in the
+                pyramid holds determinants of Hessian matrices for the original 
+                input image.  This object can be used to find stable interest
+                points in an image.  For further details consult the following
+                papers.
+
+                This object is an implementation of the fast Hessian pyramid 
+                as described in the paper: 
+                   SURF: Speeded Up Robust Features
+                   By Herbert Bay1 , Tinne Tuytelaars2 , and Luc Van Gool12
+
+                This implementation was also influenced by the very well documented
+                OpenSURF library and its corresponding description of how the fast
+                Hessian algorithm functions:  
+                    Notes on the OpenSURF Library
+                    Christopher Evans
+        !*/
+    public:
+
+        template <typename integral_image_type>
+        void build_pyramid (
+            const integral_image_type& img,
+            long num_octaves,
+            long num_intervals,
+            long initial_step_size
+        );
+        /*!
+            requires
+                - num_octaves > 0
+                - num_intervals > 0
+                - initial_step_size > 0
+                - integral_image_type == an object such as dlib::integral_image or another
+                  type that implements the interface defined in image_transforms/integral_image_abstract.h
+            ensures
+                - #get_step_size(0) == initial_step_size
+                - #octaves() == num_octaves
+                - #intervals() == num_intervals
+                - creates a Hessian pyramid from the given input image.  
+        !*/
+
+        long octaves (
+        ) const;
+        /*!
+            ensures
+                - returns the number of octaves in this pyramid
+        !*/
+
+        long intervals (
+        ) const; 
+        /*!
+            ensures
+                - returns the number of intervals in this pyramid
+        !*/
+
+        long get_border_size (
+            long octave
+        ) const;
+        /*!
+            requires
+                - 0 <= octave < octaves()
+            ensures
+                - Each octave of the pyramid has a certain sized border region where we
+                  can't compute the Hessian values since they are too close to the edge
+                  of the input image.  This function returns the size of that border.
+        !*/
+
+        long get_step_size (
+            long octave
+        ) const;
+        /*!
+            requires
+                - 0 <= octave < octaves()
+            ensures
+                - Each octave has a step size value.  This value determines how many
+                  input image pixels separate each pixel in the given pyramid octave.
+                  As the octave gets larger (i.e. as it goes to the top of the pyramid) the
+                  step size gets bigger and thus the pyramid narrows.
+        !*/
+
+        long nr (
+            long octave
+        ) const;
+        /*!
+            requires
+                - 0 <= octave < octaves()
+            ensures
+                - returns the number of rows there are per layer in the given 
+                  octave of pyramid
+        !*/
+
+        long nc (
+            long octave
+        ) const;
+        /*!
+            requires
+                - 0 <= octave < octaves()
+            ensures
+                - returns the number of columns there are per layer in the given 
+                  octave of pyramid
+        !*/
+
+        double get_value (
+            long octave,
+            long interval,
+            long r,
+            long c
+        ) const;
+        /*!
+            requires
+                - 0 <= octave < octaves()
+                - 0 <= interval < intervals()
+                - Let BS == get_border_size(octave): then
+                    - BS <= r < nr(octave)-BS
+                    - BS <= c < nc(octave)-BS
+            ensures
+                - returns the determinant of the Hessian from the given octave and interval
+                  of the pyramid.  The specific point sampled at this pyramid level is
+                  the one that corresponds to the input image point (point(r,c)*get_step_size(octave)).
+        !*/
+
+        double get_laplacian (
+            long octave,
+            long interval,
+            long r,
+            long c
+        ) const;
+        /*!
+            requires
+                - 0 <= octave < octaves()
+                - 0 <= interval < intervals()
+                - Let BS == get_border_size(octave): then
+                    - BS <= r < nr(octave)-BS
+                    - BS <= c < nc(octave)-BS
+            ensures
+                - returns the sign of the laplacian for the given octave and interval
+                  of the pyramid.  The specific point sampled at this pyramid level is
+                  the one that corresponds to the input image point (point(r,c)*get_step_size(octave)).
+                - The laplacian is the trace of the Hessian at the given point.  So this 
+                  function returns either +1 or -1 depending on this number's sign.  This
+                  value can be used to distinguish bright blobs on dark backgrounds from
+                  the reverse.
+        !*/
+
+    };
+
+// ----------------------------------------------------------------------------------------
+
+    struct interest_point
+    {
+        /*!
+            WHAT THIS OBJECT REPRESENTS
+                This object contains the interest points found using the 
+                hessian_pyramid object.  Its fields have the following
+                meanings:
+                    - center == the x/y location of the center of the interest point
+                    - scale == the scale at which the point was detected
+                    - score == the determinant of the Hessian for the interest point
+                    - laplacian == the sign of the laplacian for the interest point
+        !*/
+
+        interest_point() : scale(0), score(0), laplacian(0) {}
+
+        dlib::vector<double,2> center;
+        double scale;
+        double score;
+        double laplacian;
+
+        bool operator < (const interest_point& p) const { return score < p.score; }
+        /*!
+            This function is here so you can sort interest points according to 
+            their scores
+        !*/
+    };
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename Alloc>
+    void get_interest_points (
+        const hessian_pyramid& pyr,
+        double threshold,
+        std::vector<interest_point,Alloc>& result_points
+    )
+    /*!
+        requires
+            - threshold >= 0
+        ensures
+            - extracts interest points from the pyramid pyr and stores them into
+              result_points (note that result_points is cleared before these new interest
+              points are added to it).
+            - Only interest points with determinant values in the pyramid larger than
+              threshold are output.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+}
+
+#endif  // DLIB_HESSIAN_PYRAMId_ABSTRACT_H__
+

dlib/image_keypoint/surf.h

+// Copyright (C) 2009  Davis E. King (davisking@users.sourceforge.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#ifndef DLIB_SURf_H_
+#define DLIB_SURf_H_
+
+#include "surf_abstract.h"
+#include "hessian_pyramid.h"
+#include "../matrix.h"
+
+namespace dlib
+{
+
+// ----------------------------------------------------------------------------------------
+    
+    struct surf_point
+    {
+        interest_point p;
+        matrix<double,64,1> des;
+        double angle;
+
+        double match_score;
+
+        bool operator < (const surf_point& p) const { return match_score < p.match_score; }
+    };
+
+// ----------------------------------------------------------------------------------------
+
+    inline double gaussian (double x, double y, double sig)
+    {
+        const double pi = 3.1415926535898;
+        return 1.0/(sig*std::sqrt(2*pi)) * std::exp( -(x*x + y*y)/(2*sig*sig));
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename integral_image_type, typename T>
+    double compute_dominant_angle (
+        const integral_image_type& img,
+        const dlib::vector<T,2>& center,
+        const double& scale
+    )
+    {
+        DLIB_ASSERT(get_rect(img).contains(centered_rect(center, 17*scale,17*scale)) == true,
+            "\tdouble compute_dominant_angle(img, center, scale)"
+            << "\n\tAll arguments to this function must be > 0"
+            << "\n\t get_rect(img): " << get_rect(img) 
+            << "\n\t center:        " << center 
+            << "\n\t scale:         " << scale 
+        );
+
+        const double pi = 3.1415926535898;
+
+        std::vector<double> ang;
+        std::vector<dlib::vector<double,2> > samples;
+
+        // accumulate a bunch of angle and vector samples
+        dlib::vector<double,2> vect;
+        for (long r = -6; r <= 6; ++r)
+        {
+            for (long c = -6; c <= 6; ++c)
+            {
+                if (r*r + c*c < 36)
+                {
+                    // compute a Gaussian weighted gradient and the gradient's angle.
+                    const double gauss = gaussian(c,r, 2.5);
+                    vect.x() = gauss*haar_x(img, scale*point(c,r)+center, static_cast<long>(4*scale+0.5));
+                    vect.y() = gauss*haar_y(img, scale*point(c,r)+center, static_cast<long>(4*scale+0.5));
+                    samples.push_back(vect);
+                    ang.push_back(atan2(vect.y(), vect.x()));
+                }
+            }
+        }
+
+        //cout << "ang size: " << ang.size() << endl;
+
+        // now find the dominant direction
+        double max_length = 0;
+        double best_ang = 0;
+        // look at a bunch of pie shaped slices of a circle 
+        const long slices = 45;
+        const double ang_step = (2*pi)/slices;
+        for (long ang_i = 0; ang_i < slices; ++ang_i)
+        {
+            // compute the bounding angles
+            double ang1 = ang_step*ang_i - pi;
+            double ang2 = ang1 + pi/3;
+
+
+            // compute sum of all vectors that are within the above two angles
+            vect.x() = 0;
+            vect.y() = 0;
+            for (unsigned long i = 0; i < ang.size(); ++i)
+            {
+                if (ang1 <= ang[i] && ang[i] <= ang2)
+                {
+                    vect += samples[i];
+                    //cout << ".";
+                }
+                else if (ang2 > pi && (ang[i] >= ang1 || ang[i] <= (-2*pi+ang2)))
+                {
+                    vect += samples[i];
+                    //cout << ".";
+                }
+            }
+            //cout << "$";
+
+
+            // record the angle of the best vectors
+            if (length_squared(vect) > max_length)
+            {
+                max_length = length_squared(vect);
+                best_ang = atan2(vect.y(), vect.x());
+            }
+        }
+
+        return best_ang;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename integral_image_type, typename T, typename MM, typename L>
+    void compute_surf_descriptor (
+        const integral_image_type& img,
+        const dlib::vector<T,2>& center,
+        const double scale,
+        const double angle,
+        matrix<double,64,1,MM,L>& des
+    )
+    {
+        DLIB_ASSERT(get_rect(img).contains(centered_rect(center, 31*scale,31*scale)) == true,
+            "\tvoid compute_surf_descriptor(img, center, scale, angle)"
+            << "\n\tAll arguments to this function must be > 0"
+            << "\n\t get_rect(img): " << get_rect(img) 
+            << "\n\t center:        " << center 
+            << "\n\t scale:         " << scale 
+        );
+
+        point_rotator rot(angle);
+        point_rotator inv_rot(-angle);
+
+        long count = 0;
+
+        // loop over the 4x4 grid of histogram buckets 
+        for (long r = -10; r < 10; r += 5)
+        {
+            for (long c = -10; c < 10; c += 5)
+            {
+                dlib::vector<double,2> vect, abs_vect, temp;
+
+                // now loop over 25 points in this bucket and sum their features 
+                for (long y = r; y < r+5; ++y)
+                {
+                    for (long x = c; x < c+5; ++x)
+                    {
+                        // get the rotated point for this extraction point
+                        point p(rot(point(x,y)*scale) + center); 
+
+                        const double gauss = gaussian(x,y, 3.3);
+                        temp.x() = gauss*haar_x(img, p, static_cast<long>(2*scale+0.5));
+                        temp.y() = gauss*haar_y(img, p, static_cast<long>(2*scale+0.5));
+
+                        // rotate this vector into alignment with the surf descriptor box 
+                        temp = inv_rot(temp);
+
+                        vect += temp;
+                        abs_vect += abs(temp);
+                    }
+                }
+
+                des(count++) = vect.x();
+                des(count++) = vect.y();
+                des(count++) = abs_vect.x();
+                des(count++) = abs_vect.y();
+            }
+        }
+
+        // return the length normalized descriptor
+        des = des/length(des);
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename image_type>
+    const std::vector<surf_point> get_surf_points (
+        const image_type& img,
+        long max_points
+    )
+    {
+        integral_image int_img;
+
+        int_img.load(img);
+        hessian_pyramid pyr;
+        pyr.build_pyramid(int_img, 4, 6, 2);
+
+        std::vector<interest_point> points; 
+        get_interest_points(pyr, 0.10, points);
+        std::vector<surf_point> spoints;
+
+        // sort all the points by how strong their detect is
+        std::sort(points.rbegin(), points.rend());
+
+
+        // now extract SURF descriptors for the points
+        for (unsigned long i = 0; i < std::min((size_t)max_points,points.size()); ++i)
+        {
+            // ignore points that are close to the edge of the image
+            const double border = 31;
+            //const double border = std::sqrt(22.0*22 + 22*22);
+            if (get_rect(int_img).contains(centered_rect(points[i].center, border*points[i].scale, border*points[i].scale)))
+            {
+                surf_point sp;
+                sp.angle = compute_dominant_angle(int_img, points[i].center, points[i].scale);
+                compute_surf_descriptor(int_img, points[i].center, points[i].scale, sp.angle, sp.des);
+                sp.p = points[i];
+
+                spoints.push_back(sp);
+
+            }
+        }
+
+        return spoints;
+    }
+
+// ----------------------------------------------------------------------------------------
+
+}
+
+#endif // DLIB_SURf_H_
+

dlib/image_keypoint/surf_abstract.h

+// Copyright (C) 2009  Davis E. King (davisking@users.sourceforge.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#undef DLIB_SURf_ABSTRACT_H_
+#ifdef DLIB_SURf_ABSTRACT_H_
+
+#include "hessian_pyramid_abstract.h"
+#include "../geometry/vector_abstract.h"
+#include "../matrix/matrix_abstract.h"
+
+namespace dlib
+{
+
+// ----------------------------------------------------------------------------------------
+    
+    double gaussian (
+        double x, 
+        double y,
+        double sig
+    );
+    /*!
+    !*/
+    {
+        const double pi = 3.1415926535898;
+        return 1.0/(sig*std::sqrt(2*pi)) * std::exp( -(x*x + y*y)/(2*sig*sig));
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename integral_image_type, typename T>
+    double compute_dominant_angle (
+        const integral_image_type& img,
+        const dlib::vector<T,2>& center,
+        const double& scale
+    );
+    /*!
+        requires
+            - integral_image_type == an object such as dlib::integral_image or another
+              type that implements the interface defined in image_transforms/integral_image_abstract.h
+            - scale > 0
+            - get_rect(img).contains(centered_rect(center, 17*scale, 17*scale)) == true
+              (i.e. center can't be within 17*scale pixels of the edge of the image)
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename integral_image_type, typename T, typename MM, typename L>
+    void compute_surf_descriptor (
+        const integral_image_type& img,
+        const dlib::vector<T,2>& center,
+        const double scale,
+        const double angle,
+        matrix<double,64,1,MM,L>& des
+    )
+    /*!
+        requires
+            - integral_image_type == an object such as dlib::integral_image or another
+              type that implements the interface defined in image_transforms/integral_image_abstract.h
+            - scale > 0
+            - get_rect(img).contains(centered_rect(center, 31*scale, 31*scale)) == true
+              (i.e. center can't be within 31*scale pixels of the edge of the image)
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    struct surf_point
+    {
+        /*!
+            WHAT THIS OBJECT REPRESENTS
+        !*/
+
+        interest_point p;
+        matrix<double,64,1> des;
+        double angle;
+
+        double match_score;
+
+        bool operator < (const surf_point& p) const { return match_score < p.match_score; }
+    };
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename image_type>
+    const std::vector<surf_point> get_surf_points (
+        const image_type& img,
+        long max_points
+    );
+    /*!
+        requires
+            - max_points > 0
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+}
+
+#endif // DLIB_SURf_ABSTRACT_H_
+
+

dlib/image_transforms.h

@@ -10,6 +10,7 @@
 #include "image_transforms/thresholding.h"
 #include "image_transforms/edge_detector.h"
 #include "image_transforms/draw.h"
+#include "image_transforms/integral_image.h"
 
 #endif // DLIB_IMAGE_TRANSFORMs_
 

dlib/image_transforms/integral_image.h

+// Copyright (C) 2009  Davis E. King (davisking@users.sourceforge.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#ifndef DLIB_INTEGRAL_IMAGE
+#define DLIB_INTEGRAL_IMAGE
+
+#include "integral_image_abstract.h"
+
+#include "../algs.h"
+#include "../assert.h"
+#include "../geometry.h"
+#include "../array2d.h"
+#include "../matrix.h"
+#include "../pixel.h"
+#include "../noncopyable.h"
+
+namespace dlib
+{
+
+// ----------------------------------------------------------------------------------------
+
+    class integral_image;
+    inline const rectangle get_rect (
+        const integral_image& img
+    );
+
+// ----------------------------------------------------------------------------------------
+
+    class integral_image : noncopyable
+    {
+    public:
+        typedef long value_type;
+
+        const long nr() const { return int_img.nr(); }
+        const long nc() const { return int_img.nc(); }
+
+        template <typename image_type>
+        void load (
+            const image_type& img
+        )
+        {
+            unsigned long pixel;
+            int_img.set_size(img.nr(), img.nc());
+
+            // compute the first row of the integral image
+            unsigned long temp = 0;
+            for (unsigned long c = 0; c < img.nc(); ++c)
+            {
+                assign_pixel(pixel, img[0][c]);
+                temp += pixel;
+                int_img[0][c] = temp;
+            }
+
+            // now compute the rest of the integral image
+            for (unsigned long r = 1; r < img.nr(); ++r)
+            {
+                temp = 0;
+                for (unsigned long c = 0; c < img.nc(); ++c)
+                {
+                    assign_pixel(pixel, img[r][c]);
+                    temp += pixel;
+                    int_img[r][c] = temp + int_img[r-1][c];
+                }
+            }
+
+        }
+
+        long get_sum_of_area (
+            const rectangle& rect
+        ) const
+        {
+            DLIB_ASSERT(get_rect(*this).contains(rect) == true,
+                "\tlong get_sum_of_area(rect)"
+                << "\n\tYou have given a rectangle that goes outside the image"
+                << "\n\tthis:            " << this
+                << "\n\trect:            " << rect 
+                << "\n\tget_rect(*this): " << get_rect(*this) 
+            );
+
+            unsigned long top_left = 0, top_right = 0, bottom_left = 0, bottom_right = 0;
+
+            bottom_right = int_img[rect.bottom()][rect.right()];
+            if (rect.left()-1 >= 0 && rect.top()-1 >= 0)
+            {
+                top_left = int_img[rect.top()-1][rect.left()-1];
+                bottom_left = int_img[rect.bottom()][rect.left()-1];
+                top_right = int_img[rect.top()-1][rect.right()];
+            }
+            else if (rect.left()-1 >= 0)
+            {
+                bottom_left = int_img[rect.bottom()][rect.left()-1];
+            }
+            else if (rect.top()-1 >= 0)
+            {
+                top_right = int_img[rect.top()-1][rect.right()];
+            }
+
+            return bottom_right - bottom_left - top_right + top_left;
+        }
+
+    private:
+
+        array2d<unsigned long>::kernel_1a int_img;
+
+
+    };
+
+// ----------------------------------------------------------------------------------------
+
+    inline const rectangle get_rect (
+        const integral_image& img
+    )
+    {
+        return rectangle(0, 0, img.nc()-1, img.nr()-1);
+    }
+
+    //  ----------------------------------------------------------------------------
+
+    template <typename integral_image_type>
+    typename integral_image_type::value_type haar_x (
+        const integral_image_type& img,
+        const point& p,
+        long width
+    )
+    {
+        DLIB_ASSERT(get_rect(img).contains(centered_rect(p,width,width)) == true,
+            "\tlong haar_x(img,p,width)"
+            << "\n\tYou have given a point and with that goes outside the image"
+            << "\n\tget_rect(img):  " << get_rect(img) 
+            << "\n\tp:              " << p 
+            << "\n\twidth:          " << width 
+        );
+
+        rectangle left_rect;
+        left_rect.set_left ( p.x() - width / 2 );
+        left_rect.set_top ( p.y() - width / 2 );
+        left_rect.set_right ( p.x()-1 );
+        left_rect.set_bottom ( left_rect.top() + width - 1 );
+
+        rectangle right_rect;
+        right_rect.set_left ( p.x() );
+        right_rect.set_top ( left_rect.top() );
+        right_rect.set_right ( left_rect.left() + width -1 );
+        right_rect.set_bottom ( left_rect.bottom() );
+
+        return img.get_sum_of_area(right_rect) - img.get_sum_of_area(left_rect);
+    }
+
+    //  ----------------------------------------------------------------------------
+
+    template <typename integral_image_type>
+    typename integral_image_type::value_type haar_y (
+        const integral_image_type& img,
+        const point& p,
+        long width
+    )
+    {
+        DLIB_ASSERT(get_rect(img).contains(centered_rect(p,width,width)) == true,
+            "\tlong haar_y(img,p,width)"
+            << "\n\tYou have given a point and with that goes outside the image"
+            << "\n\tget_rect(img):  " << get_rect(img) 
+            << "\n\tp:              " << p 
+            << "\n\twidth:          " << width 
+        );
+
+        rectangle top_rect;
+        top_rect.set_left ( p.x() - width / 2 );
+        top_rect.set_top ( p.y() - width / 2 );
+        top_rect.set_right ( top_rect.left() + width - 1 );
+        top_rect.set_bottom ( p.y()-1 );
+
+        rectangle bottom_rect;
+        bottom_rect.set_left ( top_rect.left() );
+        bottom_rect.set_top ( p.y() );
+        bottom_rect.set_right ( top_rect.right() );
+        bottom_rect.set_bottom ( top_rect.top() + width - 1 );
+
+        return img.get_sum_of_area(bottom_rect) - img.get_sum_of_area(top_rect);
+    }
+
+// ----------------------------------------------------------------------------------------
+
+}
+
+#endif // DLIB_INTEGRAL_IMAGE
+

dlib/image_transforms/integral_image_abstract.h

+// Copyright (C) 2009  Davis E. King (davisking@users.sourceforge.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#undef DLIB_INTEGRAL_IMAGe_ABSTRACT_
+#ifdef DLIB_INTEGRAL_IMAGe_ABSTRACT_
+
+#include "../geometry/rectangle_abstract.h"
+#include "../array2d/array2d_kernel_abstract.h"
+#include "../pixel.h"
+#include "../noncopyable.h"
+
+namespace dlib
+{
+
+// ----------------------------------------------------------------------------------------
+
+    class integral_image : noncopyable
+    {
+        /*!
+            INITIAL VALUE
+                - nr() == 0
+                - nc() == 0
+
+            WHAT THIS OBJECT REPRESENTS
+                This object is an alternate way of representing image data
+                that allows for very fast computations of sums of pixels in 
+                rectangular regions.  To use this object you load it with a
+                normal image and then you can use the get_sum_of_area()
+                function to compute sums of pixels in a given area in
+                constant time.
+        !*/
+    public:
+        typedef long value_type;
+
+        const long nr(
+        ) const;
+        /*!
+            ensures
+                - returns the number of rows in this integral image object
+        !*/
+
+        const long nc(
+        ) const;
+        /*!
+            ensures
+                - returns the number of columns in this integral image object
+        !*/
+
+        template <typename image_type>
+        void load (
+            const image_type& img
+        );
+        /*!
+            requires
+                - image_type == a type that implements the array2d/array2d_kernel_abstract.h interface
+                - pixel_traits<image_type::type> must be defined
+            ensures
+                - #nr() == img.nr()
+                - #nc() == img.nc()
+                - #*this will now contain an "integral image" representation of the
+                  given input image.  
+        !*/
+
+        value_type get_sum_of_area (
+            const rectangle& rect
+        ) const;
+        /*!
+            requires
+                - get_rect(*this).contains(rect) == true
+                  (i.e. rect must not be outside the integral image)
+            ensures
+                - Let O denote the image this integral image was generated from.
+                  Then this function returns sum(subm(array_to_matrix(O),rect)).
+                  That is, this function returns the sum of the pixels in O that
+                  are contained within the given rectangle.
+        !*/
+
+    };
+
+// ----------------------------------------------------------------------------------------
+
+    const rectangle get_rect (
+        const integral_image& img
+    ); 
+    /*!
+        ensures
+            - returns rectangle(0, 0, img.nc()-1, img.nr()-1)
+              (i.e. returns a rectangle that has the same dimensions as
+              the integral_image img)
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename integral_image_type>
+    typename integral_image_type::value_type haar_x (
+        const integral_image_type& img,
+        const point& p,
+        long width
+    )
+    /*!
+        requires
+            - get_rect(img).contains(centered_rect(p,width,width)) == true
+            - integral_image_type == a type that implements the integral_image interface 
+              defined above
+        ensures
+            - returns the response of a Haar wavelet centered at the point p
+              with the given width.  The wavelet is oriented along the X axis
+              and has the following shape:
+                ----++++
+                ----++++
+                ----++++
+                ----++++
+              That is, the wavelet is square and computes the sum of pixels on the
+              right minus the sum of pixels on the left.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <typename integral_image_type>
+    typename integral_image_type::value_type haar_y (
+        const integral_image_type& img,
+        const point& p,
+        long width
+    )
+    /*!
+        requires
+            - get_rect(img).contains(centered_rect(p,width,width)) == true
+            - integral_image_type == a type that implements the integral_image interface 
+              defined above
+        ensures
+            - returns the response of a Haar wavelet centered at the point p
+              with the given width in the given image.  The wavelet is oriented 
+              along the Y axis and has the following shape:
+                --------
+                --------
+                ++++++++
+                ++++++++
+              That is, the wavelet is square and computes the sum of pixels on the
+              bottom minus the sum of pixels on the top.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+}
+
+#endif // DLIB_INTEGRAL_IMAGe_ABSTRACT_
+