// Copyright (C) 2009 Davis E. King (davis@dlib.net)
// License: Boost Software License See LICENSE.txt for the full license.
#ifndef DLIB_HESSIAN_PYRAMId_Hh_
#define DLIB_HESSIAN_PYRAMId_Hh_
#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; }
};
// ----------------------------------------------------------------------------------------
inline void serialize(
const interest_point& item,
std::ostream& out
)
{
try
{
serialize(item.center,out);
serialize(item.scale,out);
serialize(item.score,out);
serialize(item.laplacian,out);
}
catch (serialization_error& e)
{
throw serialization_error(e.info + "\n while serializing object of type interest_point");
}
}
// ----------------------------------------------------------------------------------------
inline void deserialize(
interest_point& item,
std::istream& in
)
{
try
{
deserialize(item.center,in);
deserialize(item.scale,in);
deserialize(item.score,in);
deserialize(item.laplacian,in);
}
catch (serialization_error& e)
{
throw serialization_error(e.info + "\n while deserializing object of type interest_point");
}
}
// ----------------------------------------------------------------------------------------
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);
for (long i = 0; i < num_intervals; ++i)
{
const long border_size = get_border_size(i)*step_size;
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));
// 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;
// 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 interval
) const
{
DLIB_ASSERT(0 <= interval && interval < intervals(),
"\tlong get_border_size(interval)"
<< "\n\tInvalid interval value"
<< "\n\t this: " << this
<< "\n\t interval: " << interval
);
const double lobe_size = 2.0*(interval+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(interval) <= r && r < nr(octave)-get_border_size(interval) &&
get_border_size(interval) <= c && c < nc(octave)-get_border_size(interval),
"\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(interval): " << get_border_size(interval)
);
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(interval) <= r && r < nr(octave)-get_border_size(interval) &&
get_border_size(interval) <= c && c < nc(octave)-get_border_size(interval),
"\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(interval): " << get_border_size(interval)
);
// 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> image_type;
typedef array<image_type> 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);
}
else
{
// this indicates to the caller that no interest point was found.
temp.score = -1;
}
return temp;
}
}
// ----------------------------------------------------------------------------------------
template <typename Alloc>
void get_interest_points (
const hessian_pyramid& pyr,
double threshold,
std::vector<interest_point,Alloc>& result_points
)
{
DLIB_ASSERT(threshold >= 0,
"\tvoid get_interest_points()"
<< "\n\t Invalid arguments to this function"
<< "\n\t threshold: " << threshold
);
using namespace std;
using namespace hessian_pyramid_helpers;
result_points.clear();
for (long o = 0; o < pyr.octaves(); ++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 += 1)
{
const long border_size = pyr.get_border_size(i+1);
for (long r = border_size+1; r < nr - border_size-1; r += 1)
{
for (long c = border_size+1; c < nc - border_size-1; c += 1)
{
double max_val = pyr.get_value(o,i,r,c);
long max_i = i;
long max_r = r;
long max_c = c;
// 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_Hh_