coco
591 lines
22 KiB
C++
591 lines
22 KiB
C++
// 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"
|
|
#include <limits>
|
|
|
|
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);
|
|
}
|
|
}
|
|
|
|
inline unsigned long size(
|
|
) const { return _size; }
|
|
|
|
long nr(
|
|
) const { return _size; }
|
|
|
|
long nc(
|
|
) const { return _size; }
|
|
|
|
std::pair<dpoint, dpoint> get_line (
|
|
const dpoint& p
|
|
) const
|
|
{
|
|
DLIB_ASSERT(rectangle(0,0,size()-1,size()-1).contains(p) == true,
|
|
"\t pair<dpoint,dpoint> hough_transform::get_line(dpoint)"
|
|
<< "\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.
|
|
double theta, radius;
|
|
get_line_properties(p, theta, radius);
|
|
theta *= pi/180;
|
|
|
|
// now make a line segment on the line.
|
|
const rectangle box = get_rect(*this);
|
|
const dpoint cent = center(box);
|
|
dpoint v1 = cent + dpoint(size()+1000,0) + dpoint(0,radius);
|
|
dpoint v2 = cent - dpoint(size()+1000,0) + dpoint(0,radius);
|
|
dpoint p1 = rotate_point(cent, v1, theta);
|
|
dpoint p2 = rotate_point(cent, v2, theta);
|
|
|
|
clip_line_to_rectangle(box, p1, p2);
|
|
|
|
return std::make_pair(p1,p2);
|
|
}
|
|
|
|
double get_line_angle_in_degrees (
|
|
const dpoint& p
|
|
) const
|
|
{
|
|
double angle, radius;
|
|
get_line_properties(p, angle, radius);
|
|
return angle;
|
|
}
|
|
|
|
void get_line_properties (
|
|
const dpoint& p,
|
|
double& angle_in_degrees,
|
|
double& radius
|
|
) const
|
|
{
|
|
const dpoint cent = center(get_rect(*this));
|
|
double theta = p.x()-cent.x();
|
|
radius = p.y()-cent.y();
|
|
angle_in_degrees = 180*theta/even_size;
|
|
radius = radius*sqrt_2 + 0.5;
|
|
}
|
|
|
|
template <
|
|
typename image_type
|
|
>
|
|
point get_best_hough_point (
|
|
const point& p,
|
|
const image_type& himg_
|
|
)
|
|
{
|
|
const const_image_view<image_type> himg(himg_);
|
|
|
|
DLIB_ASSERT(himg.nr() == size() && himg.nc() == size() &&
|
|
rectangle(0,0,size()-1,size()-1).contains(p) == true,
|
|
"\t point hough_transform::get_best_hough_point()"
|
|
<< "\n\t Invalid arguments given to this function."
|
|
<< "\n\t himg.nr(): " << himg.nr()
|
|
<< "\n\t himg.nc(): " << himg.nc()
|
|
<< "\n\t size(): " << size()
|
|
<< "\n\t p: " << p
|
|
);
|
|
|
|
|
|
typedef typename image_traits<image_type>::pixel_type pixel_type;
|
|
COMPILE_TIME_ASSERT(pixel_traits<pixel_type>::grayscale == true);
|
|
pixel_type best_val = std::numeric_limits<pixel_type>::min();
|
|
point best_point;
|
|
|
|
|
|
const long max_n8 = (himg.nc()/8)*8;
|
|
const long max_n4 = (himg.nc()/4)*4;
|
|
const long r = p.y();
|
|
const long c = p.x();
|
|
|
|
const int32* ysin = &ysin_theta(r,0);
|
|
const int32* xcos = &xcos_theta(c,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;
|
|
|
|
if (himg[rr0][t++] > best_val)
|
|
{
|
|
best_val = himg[rr0][t-1];
|
|
best_point.x() = t-1;
|
|
best_point.y() = rr0;
|
|
}
|
|
if (himg[rr1][t++] > best_val)
|
|
{
|
|
best_val = himg[rr1][t-1];
|
|
best_point.x() = t-1;
|
|
best_point.y() = rr1;
|
|
}
|
|
if (himg[rr2][t++] > best_val)
|
|
{
|
|
best_val = himg[rr2][t-1];
|
|
best_point.x() = t-1;
|
|
best_point.y() = rr2;
|
|
}
|
|
if (himg[rr3][t++] > best_val)
|
|
{
|
|
best_val = himg[rr3][t-1];
|
|
best_point.x() = t-1;
|
|
best_point.y() = rr3;
|
|
}
|
|
if (himg[rr4][t++] > best_val)
|
|
{
|
|
best_val = himg[rr4][t-1];
|
|
best_point.x() = t-1;
|
|
best_point.y() = rr4;
|
|
}
|
|
if (himg[rr5][t++] > best_val)
|
|
{
|
|
best_val = himg[rr5][t-1];
|
|
best_point.x() = t-1;
|
|
best_point.y() = rr5;
|
|
}
|
|
if (himg[rr6][t++] > best_val)
|
|
{
|
|
best_val = himg[rr6][t-1];
|
|
best_point.x() = t-1;
|
|
best_point.y() = rr6;
|
|
}
|
|
if (himg[rr7][t++] > best_val)
|
|
{
|
|
best_val = himg[rr7][t-1];
|
|
best_point.x() = t-1;
|
|
best_point.y() = rr7;
|
|
}
|
|
}
|
|
while(t < max_n4)
|
|
{
|
|
long rr0 = (*xcos++ + *ysin++)>>16;
|
|
long rr1 = (*xcos++ + *ysin++)>>16;
|
|
long rr2 = (*xcos++ + *ysin++)>>16;
|
|
long rr3 = (*xcos++ + *ysin++)>>16;
|
|
if (himg[rr0][t++] > best_val)
|
|
{
|
|
best_val = himg[rr0][t-1];
|
|
best_point.x() = t-1;
|
|
best_point.y() = rr0;
|
|
}
|
|
if (himg[rr1][t++] > best_val)
|
|
{
|
|
best_val = himg[rr1][t-1];
|
|
best_point.x() = t-1;
|
|
best_point.y() = rr1;
|
|
}
|
|
if (himg[rr2][t++] > best_val)
|
|
{
|
|
best_val = himg[rr2][t-1];
|
|
best_point.x() = t-1;
|
|
best_point.y() = rr2;
|
|
}
|
|
if (himg[rr3][t++] > best_val)
|
|
{
|
|
best_val = himg[rr3][t-1];
|
|
best_point.x() = t-1;
|
|
best_point.y() = rr3;
|
|
}
|
|
}
|
|
while(t < himg.nc())
|
|
{
|
|
long rr0 = (*xcos++ + *ysin++)>>16;
|
|
if (himg[rr0][t++] > best_val)
|
|
{
|
|
best_val = himg[rr0][t-1];
|
|
best_point.x() = t-1;
|
|
best_point.y() = rr0;
|
|
}
|
|
}
|
|
|
|
return best_point;
|
|
}
|
|
|
|
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 void hough_transform::operator()"
|
|
<< "\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);
|
|
|
|
auto record_hit = [&](const point& hough_point, const point& /*img_point*/, const in_pixel_type& val)
|
|
{
|
|
himg[hough_point.y()][hough_point.x()] += val;
|
|
};
|
|
perform_generic_hough_transform(img_, box, record_hit);
|
|
}
|
|
|
|
template <
|
|
typename in_image_type,
|
|
typename out_image_type
|
|
>
|
|
void operator() (
|
|
const in_image_type& img_,
|
|
out_image_type& himg_
|
|
) const
|
|
{
|
|
rectangle box(0,0, num_columns(img_)-1, num_rows(img_)-1);
|
|
(*this)(img_, box, himg_);
|
|
}
|
|
|
|
template <
|
|
typename in_image_type
|
|
>
|
|
std::vector<std::vector<point>> find_pixels_voting_for_lines (
|
|
const in_image_type& img,
|
|
const rectangle& box,
|
|
const std::vector<point>& hough_points,
|
|
const unsigned long angle_window_size = 1,
|
|
const unsigned long radius_window_size = 1
|
|
) const
|
|
{
|
|
|
|
typedef typename image_traits<in_image_type>::pixel_type in_pixel_type;
|
|
|
|
DLIB_CASSERT(angle_window_size >= 1);
|
|
DLIB_CASSERT(radius_window_size >= 1);
|
|
DLIB_CASSERT(box.width() == size() && box.height() == size(),
|
|
"\t std::vector<std::vector<point>> hough_transform::find_pixels_voting_for_lines()"
|
|
<< "\n\t Invalid arguments given to this function."
|
|
<< "\n\t box.width(): " << box.width()
|
|
<< "\n\t box.height(): " << box.height()
|
|
<< "\n\t size(): " << size()
|
|
);
|
|
#ifdef ENABLE_ASSERTS
|
|
for (auto& p : hough_points)
|
|
DLIB_CASSERT(get_rect(*this).contains(p),
|
|
"You gave a hough_points that isn't actually in the Hough space of this object."
|
|
<< "\n\t get_rect(*this): "<< get_rect(*this)
|
|
<< "\n\t p: "<< p
|
|
);
|
|
#endif
|
|
|
|
std::vector<std::vector<point>> constituent_points(hough_points.size());
|
|
|
|
// make a map that lets us look up in constant time if a hough point is in the
|
|
// constituent_points output and if so where.
|
|
matrix<uint32> hmap(size(),size());
|
|
hmap = hough_points.size();
|
|
for (size_t i = 0; i < hough_points.size(); ++i)
|
|
{
|
|
rectangle area = centered_rect(hough_points[i],angle_window_size,radius_window_size).intersect(get_rect(hmap));
|
|
for (long r = area.top(); r <= area.bottom(); ++r)
|
|
{
|
|
for (long c = area.left(); c <= area.right(); ++c)
|
|
{
|
|
hmap(r,c) = i;
|
|
}
|
|
}
|
|
}
|
|
|
|
// record that this image point voted for this Hough point
|
|
auto record_hit = [&](const point& hough_point, const point& img_point, in_pixel_type)
|
|
{
|
|
auto idx = hmap(hough_point.y(), hough_point.x());
|
|
if (idx < constituent_points.size())
|
|
{
|
|
// don't add img_point if it's already in the list.
|
|
if (constituent_points[idx].size() == 0 || constituent_points[idx].back() != img_point)
|
|
constituent_points[idx].push_back(img_point);
|
|
}
|
|
};
|
|
|
|
perform_generic_hough_transform(img, box, record_hit);
|
|
|
|
return constituent_points;
|
|
}
|
|
|
|
template <
|
|
typename in_image_type
|
|
>
|
|
std::vector<std::vector<point>> find_pixels_voting_for_lines (
|
|
const in_image_type& img,
|
|
const std::vector<point>& hough_points,
|
|
const unsigned long angle_window_size = 1,
|
|
const unsigned long radius_window_size = 1
|
|
) const
|
|
{
|
|
rectangle box(0,0, num_columns(img)-1, num_rows(img)-1);
|
|
return find_pixels_voting_for_lines(img, box, hough_points, angle_window_size, radius_window_size);
|
|
}
|
|
|
|
template <
|
|
typename image_type,
|
|
typename thresh_type
|
|
>
|
|
std::vector<point> find_strong_hough_points(
|
|
const image_type& himg_,
|
|
const thresh_type hough_count_threshold,
|
|
const double angle_nms_thresh,
|
|
const double radius_nms_thresh
|
|
)
|
|
{
|
|
const_image_view<image_type> himg(himg_);
|
|
|
|
DLIB_CASSERT(himg.nr() == size());
|
|
DLIB_CASSERT(himg.nc() == size());
|
|
DLIB_CASSERT(angle_nms_thresh >= 0)
|
|
DLIB_CASSERT(radius_nms_thresh >= 0)
|
|
|
|
std::vector<std::pair<double,point>> initial_lines;
|
|
for (long r = 0; r < himg.nr(); ++r)
|
|
{
|
|
for (long c = 0; c < himg.nc(); ++c)
|
|
{
|
|
if (himg[r][c] >= hough_count_threshold)
|
|
initial_lines.emplace_back(himg[r][c], point(c,r));
|
|
}
|
|
}
|
|
|
|
|
|
std::vector<point> final_lines;
|
|
std::vector<std::pair<double,double>> final_angle_and_radius;
|
|
|
|
// Now do non-max suppression. First, sort the initial_lines so the best lines come first.
|
|
std::sort(initial_lines.rbegin(), initial_lines.rend(),
|
|
[](const std::pair<double,point>& a, const std::pair<double,point>& b){ return a.first<b.first;});
|
|
for (auto& r : initial_lines)
|
|
{
|
|
double angle, radius;
|
|
get_line_properties(r.second, angle, radius);
|
|
|
|
// check if anything in final_lines is too close to r.second. If
|
|
// something is found then discard r.second.
|
|
auto too_close = false;
|
|
for (auto& ref : final_angle_and_radius)
|
|
{
|
|
auto& ref_angle = ref.first;
|
|
auto& ref_radius = ref.second;
|
|
|
|
// We need to check for wrap around in angle since, for instance, a
|
|
// line with angle and radius of 90 and 10 is the same line as one with
|
|
// angle -90 and radius -10.
|
|
if ((std::abs(ref_angle - angle) < angle_nms_thresh && std::abs(ref_radius-radius) < radius_nms_thresh) ||
|
|
(180 - std::abs(ref_angle - angle) < angle_nms_thresh && std::abs(ref_radius+radius) < radius_nms_thresh))
|
|
{
|
|
too_close = true;
|
|
break;
|
|
}
|
|
}
|
|
|
|
if (!too_close)
|
|
{
|
|
final_lines.emplace_back(r.second);
|
|
final_angle_and_radius.emplace_back(angle, radius);
|
|
}
|
|
}
|
|
|
|
return final_lines;
|
|
}
|
|
|
|
|
|
template <
|
|
typename in_image_type,
|
|
typename record_hit_function_type
|
|
>
|
|
void perform_generic_hough_transform (
|
|
const in_image_type& img_,
|
|
const rectangle& box,
|
|
record_hit_function_type record_hit
|
|
) const
|
|
{
|
|
|
|
typedef typename image_traits<in_image_type>::pixel_type in_pixel_type;
|
|
|
|
DLIB_ASSERT(box.width() == size() && box.height() == size(),
|
|
"\t void hough_transform::perform_generic_hough_transform()"
|
|
<< "\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);
|
|
|
|
|
|
const_image_view<in_image_type> img(img_);
|
|
|
|
|
|
const rectangle area = box.intersect(get_rect(img));
|
|
|
|
const long max_n8 = (size()/8)*8;
|
|
const long max_n4 = (size()/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 auto val = 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 < size(); ++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);
|
|
|
|
record_hit(point(t,rr), point(c,r), 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;
|
|
|
|
record_hit(point(t++,rr0), point(c,r), val);
|
|
record_hit(point(t++,rr1), point(c,r), val);
|
|
record_hit(point(t++,rr2), point(c,r), val);
|
|
record_hit(point(t++,rr3), point(c,r), val);
|
|
record_hit(point(t++,rr4), point(c,r), val);
|
|
record_hit(point(t++,rr5), point(c,r), val);
|
|
record_hit(point(t++,rr6), point(c,r), val);
|
|
record_hit(point(t++,rr7), point(c,r), 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;
|
|
record_hit(point(t++,rr0), point(c,r), val);
|
|
record_hit(point(t++,rr1), point(c,r), val);
|
|
record_hit(point(t++,rr2), point(c,r), val);
|
|
record_hit(point(t++,rr3), point(c,r), val);
|
|
}
|
|
while(t < (long)size())
|
|
{
|
|
long rr0 = (*xcos++ + *ysin++)>>16;
|
|
record_hit(point(t++,rr0), point(c,r), val);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
template <
|
|
typename in_image_type,
|
|
typename record_hit_function_type
|
|
>
|
|
void perform_generic_hough_transform (
|
|
const in_image_type& img_,
|
|
record_hit_function_type record_hit
|
|
) const
|
|
{
|
|
rectangle box(0,0, num_columns(img_)-1, num_rows(img_)-1);
|
|
perform_generic_hough_transform(img_, box, record_hit);
|
|
}
|
|
|
|
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_
|
|
|