Commit 32e461a6 authored by Guillaume Lazzara's avatar Guillaume Lazzara
Browse files

Optimize Sauvola's multiscale binarization.

	* core/init_integral_image.hh: New.

	* src/binarization/sauvola_ms.cc: Optimize and make it more
	robust.

	* canvas/integral_browsing.hh: New canvas to browse and compute
	data in an integral image.

	* binarization/internal/first_pass_functor.hh: New functor to be
	used in the integral browsing.

	* binarization/sauvola_threshold.hh: Add new overloads.

	* subsampling/integral_single_image.hh: Subsample an image and
	compute integral images at the same time.
parent 4bfd663c
2009-12-04 Guillaume Lazzara <z@lrde.epita.fr>
Optimize Sauvola's multiscale binarization.
* core/init_integral_image.hh: New.
* src/binarization/sauvola_ms.cc: Optimize and make it more
robust.
* canvas/integral_browsing.hh: New canvas to browse and compute
data in an integral image.
* binarization/internal/first_pass_functor.hh: New functor to be
used in the integral browsing.
* binarization/sauvola_threshold.hh: Add new overloads.
* subsampling/integral_single_image.hh: Subsample an image and
compute integral images at the same time.
2009-11-03 Guillaume Lazzara <z@lrde.epita.fr>
Add a new example in Scribo.
......
// Copyright (C) 2009 EPITA Research and Development Laboratory (LRDE)
//
// This file is part of Olena.
//
// Olena is free software: you can redistribute it and/or modify it under
// the terms of the GNU General Public License as published by the Free
// Software Foundation, version 2 of the License.
//
// Olena is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Olena. If not, see <http://www.gnu.org/licenses/>.
//
// As a special exception, you may use this file as part of a free
// software project without restriction. Specifically, if other files
// instantiate templates or use macros or inline functions from this
// file, or you compile this file and link it with other files to produce
// an executable, this file does not by itself cause the resulting
// executable to be covered by the GNU General Public License. This
// exception does not however invalidate any other reasons why the
// executable file might be covered by the GNU General Public License.
#ifndef SCRIBO_BINARIZATION_INTERNAL_FIRST_PASS_FUNCTOR_HH
# define SCRIBO_BINARIZATION_INTERNAL_FIRST_PASS_FUNCTOR_HH
# include <mln/value/int_u8.hh>
# include <scribo/binarization/sauvola_threshold.hh>
// #include <mln/border/adjust.hh>
namespace scribo
{
namespace binarization
{
namespace internal
{
using namespace mln;
unsigned my_find_root(image2d<unsigned>& parent, unsigned x)
{
if (parent.element(x) == x)
return x;
return parent.element(x) = my_find_root(parent,
parent.element(x));
}
template <typename I>
struct first_pass_functor
{
const I& input;
mln_fwd_pixter(const I) pxl;
double res;
image2d<unsigned> parent;
image2d<unsigned> card;
image2d<bool> msk;
image2d<value::int_u8> t_sub;
unsigned n_nbhs;
util::array<int> dp;
static const double one_k = 1 - 0.34;
static const double k_R = 0.34 / 128.0;
first_pass_functor(const I& input)
: input(input),
pxl(input)
{
res = 0;
pxl.start();
initialize(t_sub, input);
initialize(parent, input);
initialize(msk, input);
extension::fill(msk, false);
initialize(card, input);
data::fill(card, 1);
dp = negative_offsets_wrt(input, c4());
n_nbhs = dp.nelements();
}
void exec(double mean, double stddev)
{
mln_precondition(pxl.is_valid());
unsigned p = pxl.offset();
// Use an inlined and developed version of sauvola's
// threshold formula.
value::int_u8 t_p = mean * (one_k + k_R * stddev);
// value::int_u8 t_p = sauvola_threshold_formula(mean, stddev);
msk.element(p) = input.element(p) < t_p;
t_sub.element(p) = t_p;
if (! msk.element(p))
{
pxl.next();
return;
}
parent.element(p) = p;
for (unsigned i = 0; i < n_nbhs; ++i)
{
unsigned n = p + dp[i];
if (! msk.element(n))
continue;
unsigned r = my_find_root(parent, n);
if (r != p)
{
parent.element(r) = p;
card.element(p) += card.element(r);
}
}
pxl.next(); // next pixel
}
void finalize()
{
mln_assertion(! pxl.is_valid());
}
};
} // end of namespace scribo::binarization::internal
} // end of namespace scribo::binarization
} // end of namespace scribo
#endif // SCRIBO_BINARIZATION_INTERNAL_FIRST_PASS_FUNCTOR_HH
......@@ -44,10 +44,7 @@
# include <mln/pw/all.hh>
# include <mln/core/routine/duplicate.hh>
// Forward declaration.
namespace mln { template <typename T> class integral_image; }
# include <scribo/core/init_integral_image.hh>
namespace scribo
{
......@@ -60,7 +57,7 @@ namespace scribo
/*! \brief Compute a threshold image with Sauvola's algorithm.
\input[in] input An image.
\input[in] window_size The window size.x
\input[in] window_size The window size.
\input[out] simple The sum of all intensities of \p input.
\input[out] squared The sum of all squared intensities of \p
input.
......@@ -68,11 +65,11 @@ namespace scribo
\return An image of local thresholds.
*/
template <typename I, typename T>
template <typename I, typename J>
mln_ch_value(I, value::int_u8)
sauvola_threshold(const Image<I>& input, unsigned window_size,
integral_image<T>& simple,
integral_image<T>& squared);
Image<J>& simple,
Image<J>& squared);
/// \overload
......@@ -106,36 +103,63 @@ namespace scribo
};
unsigned long long square_(const value::int_u8& val)
/*! \brief compute Sauvola's threshold applying directly the formula.
\param[in] m_x_y Mean value.
\param[in] s_x_y Standard deviation.
\param[in] k Control the threshold value in the local
window. The higher, the lower the threshold
form the local mean m(x, y).
\param[in] R Maximum value of the standard deviation (128
for grayscale documents).
\return A threshold.
*/
inline
double
sauvola_threshold_formula(const double m_x_y, const double s_x_y,
const double k, const double R)
{
unsigned long long v = static_cast<unsigned long long>(val);
return v * v;
return m_x_y * (1.0 + k * ((s_x_y / R) - 1.0));
}
unsigned long long identity_(const value::int_u8& val)
/// \overload
//
inline
double
sauvola_threshold_formula(double m_x_y, double s_x_y)
{
return static_cast<unsigned long long>(val);
// Badekas et al. said 0.34 was best.
const double k = 0.34;
// 128 is best for grayscale documents.
const double R = 128;
return sauvola_threshold_formula(m_x_y, s_x_y, k, R);
}
/// \brief Compute a point wise threshold according Sauvola's
/// binarization.
///
/// \param[in] p A site.
/// \param[in] simple An integral image of mean values.
/// \param[in] squared An integral image of squared mean values.
/// \param[in] win_width Window width.
/// \param[in] k Control the threshold value in the local
/// window. The higher, the lower the threshold
/// form the local mean m(x, y).
/// \param[in] R Maximum value of the standard deviation (128
/// for grayscale documents).
template <typename P, typename T>
/*! \brief Compute a point wise threshold according Sauvola's
binarization.
\param[in] p A site.
\param[in] simple An integral image of mean values.
\param[in] squared An integral image of squared mean values.
\param[in] win_width Window width.
\param[in] k Control the threshold value in the local
window. The higher, the lower the threshold
form the local mean m(x, y).
\param[in] R Maximum value of the standard deviation (128
for grayscale documents).
\return A threshold.
*/
template <typename P, typename J>
double
compute_sauvola_threshold(const P& p,
const integral_image<T>& simple,
const integral_image<T>& squared,
const J& simple,
const J& squared,
int win_width, double k, double R)
{
mln_precondition(simple.nrows() == squared.nrows());
......@@ -146,142 +170,116 @@ namespace scribo
int row_min = std::max(0, p.row() - w_2);
int col_min = std::max(0, p.col() - w_2);
int row_max = std::min(simple.nrows() - 1, p.row() + w_2);
int col_max = std::min(simple.ncols() - 1, p.col() + w_2);
int row_max = std::min(static_cast<int>(simple.nrows()) - 1,
p.row() + w_2);
int col_max = std::min(static_cast<int>(simple.ncols()) - 1,
p.col() + w_2);
// std::cout << "sauvola threshold : "
// << simple.domain() << " - "
// << row_max << " - "
// << col_max << std::endl;
double wh = (row_max - row_min + 1) * (col_max - col_min + 1);
// Mean.
double m_x_y_tmp = (simple(row_max, col_max)
+ simple(row_min, col_min)
- simple(row_max, col_min)
- simple(row_min, col_max));
double m_x_y_tmp = (simple.at_(row_max, col_max)
+ simple.at_(row_min, col_min)
- simple.at_(row_max, col_min)
- simple.at_(row_min, col_max));
double m_x_y = m_x_y_tmp / wh;
// Standard deviation.
double s_x_y_tmp = (squared(row_max, col_max)
+ squared(row_min, col_min)
- squared(row_max, col_min)
- squared(row_min, col_max));
double s_x_y_tmp = (squared.at_(row_max, col_max)
+ squared.at_(row_min, col_min)
- squared.at_(row_max, col_min)
- squared.at_(row_min, col_max));
double s_x_y = std::sqrt((s_x_y_tmp - (m_x_y_tmp * m_x_y_tmp) / wh) / (wh - 1.f));
// Thresholding.
double t_x_y = m_x_y * (1.0 + k * ((s_x_y / R) - 1.0));
double t_x_y = sauvola_threshold_formula(m_x_y, s_x_y, k, R);
return t_x_y;
}
template <typename P, typename T>
template <typename P, typename J>
double
compute_sauvola_threshold(const P& p,
const integral_image<T>& simple,
const integral_image<T>& squared,
int win_width)
compute_sauvola_threshold_single_image(const P& p,
const J& integral,
int win_width)
{
// Badekas et al. said 0.34 was best.
const double k = 0.34;
// 128 is best for grayscale documents.
const double R = 128;
return compute_sauvola_threshold(p, simple, squared, win_width, k, R);
}
} // end of namespace scribo::binarization::internal
} // end of namespace scribo::binarization
// Window half width.
int w_2 = win_width >> 1;
} // end of namespace scribo
int row_min = std::max(0, p.row() - w_2);
int col_min = std::max(0, p.col() - w_2);
//FIXME: offset (-4) should be replaced by the number of
//padding pixels.
int row_max = std::min(static_cast<int>(integral.nrows()) - 1,
p.row() + w_2);
int col_max = std::min(static_cast<int>(integral.ncols()) - 1,
p.col() + w_2);
namespace mln
{
template <typename T>
class integral_image
{
public:
// std::cout << "sauvola threshold : "
// << simple.domain() << " - "
// << row_max << " - "
// << col_max << std::endl;
integral_image()
: img_(0), nrows_(0), ncols_(0)
{}
template<class F>
integral_image(const image2d<T>& i, F func)
: img_(0),
nrows_(0),
ncols_(0)
{
init_(i, func);
}
double wh = (row_max - row_min + 1) * (col_max - col_min + 1);
template<class F>
void init_(const image2d<T>& i, F func)
{
nrows_ = i.nrows();
ncols_ = i.ncols();
// Mean.
double m_x_y_tmp = (integral.at_(row_max, col_max).first()
+ integral.at_(row_min, col_min).first()
- integral.at_(row_max, col_min).first()
- integral.at_(row_min, col_max).first());
img_ = (unsigned long long**)malloc(sizeof(unsigned long long*) * nrows_);
for (int n = 0; n < nrows_; ++n)
img_[n] = (unsigned long long*)malloc(sizeof(unsigned long long) * ncols_);
double m_x_y = m_x_y_tmp / wh;
img_[0][0] = func(i.at_(0, 0));
// Standard deviation.
double s_x_y_tmp = (integral.at_(row_max, col_max).second()
+ integral.at_(row_min, col_min).second()
- integral.at_(row_max, col_min).second()
- integral.at_(row_min, col_max).second());
for (int row = 1; row < nrows_; ++row)
img_[row][0] = (*this)(row - 1, 0) + func(i.at_(row, 0));
double s_x_y = std::sqrt((s_x_y_tmp - (m_x_y_tmp * m_x_y_tmp) / wh) / (wh - 1.f));
for (int col = 1; col < ncols_; ++col)
img_[0][col] = (*this)(0, col - 1)
+ func(i.at_(0, col));
// Thresholding.
double t_x_y = m_x_y * (1.0 + 0.14 * ((s_x_y / 128) - 1.0));
for (int row = 1; row < nrows_; ++row)
for (int col = 1; col < ncols_; ++col)
img_[row][col] = (*this)(row - 1, col)
+ (*this)(row, col - 1)
- (*this)(row - 1, col - 1)
+ func(i.at_(row, col));
}
return t_x_y;
}
~integral_image()
{
for (int n = 0; n < nrows_; ++n)
free(img_[n]);
free(img_);
}
bool is_valid() const
{
return img_ != 0;
}
template <typename P, typename J>
double
compute_sauvola_threshold(const P& p,
const J& simple,
const J& squared,
int win_width)
{
// Badekas et al. said 0.34 was best.
const double k = 0.34;
unsigned long long operator()(int row, int col) const
{
return img_[row][col];
}
// 128 is best for grayscale documents.
const double R = 128;
int nrows() const
{
return nrows_;
}
return compute_sauvola_threshold(p, simple, squared, win_width, k, R);
}
int ncols() const
{
return ncols_;
}
} // end of namespace scribo::binarization::internal
private:
unsigned long long** img_;
int nrows_;
int ncols_;
};
} // end of namespace scribo::binarization
} // end of namespace mln
} // end of namespace scribo
......@@ -300,23 +298,25 @@ namespace scribo
namespace generic
{
template <typename I, typename T>
template <typename I, typename J>
inline
mln_concrete(I)
sauvola_threshold(const I& input, unsigned window_size,
integral_image<T>& simple,
integral_image<T>& squared)
sauvola_threshold(const Image<I>& input_, unsigned window_size,
Image<J>& simple_,
Image<J>& squared_)
{
trace::entering("scribo::binarization::impl::generic::sauvola_threshold");
typedef mln_value(I) V;
typedef mln_site(I) P;
const I& input = exact(input_);
J& simple = exact(simple_);
J& squared = exact(squared_);
// Compute the sum of all intensities of input
simple.init_(input, internal::identity_);
mln_assertion(input.is_valid());
mln_assertion(simple.is_valid());
mln_assertion(squared.is_valid());
// Compute the sum of all squared intensities of input
squared.init_(input, internal::square_);
typedef mln_value(I) V;
typedef mln_site(I) P;
// Savaula Algorithm with I.I.
......@@ -340,24 +340,24 @@ namespace scribo
template <typename I, typename T>
template <typename I, typename J>
inline
mln_concrete(I)
sauvola_threshold_gl(const I& input, unsigned window_size,
integral_image<T>& simple,
integral_image<T>& squared)
Image<J>& simple,
Image<J>& squared)
{
return impl::generic::sauvola_threshold(input, window_size,
simple, squared);
}
template <typename I, typename T>
template <typename I, typename J>
inline
mln_ch_value(I, value::int_u8)
sauvola_threshold_rgb8(const I& input, unsigned window_size,
integral_image<T>& simple,
integral_image<T>& squared)
Image<J>& simple,
Image<J>& squared)
{
trace::entering("scribo::binarization::impl::sauvola_threshold_rgb8");
......@@ -384,36 +384,36 @@ namespace scribo
namespace internal
{
template <unsigned n, typename I, typename T>
template <unsigned n, typename I, typename J>
inline
mln_ch_value(I, value::int_u<n>)
sauvola_threshold_dispatch(const value::int_u<n>&, const I& input,
unsigned window_size,
integral_image<T>& simple,
integral_image<T>& squared)
J& simple,
J& squared)
{
return impl::sauvola_threshold_gl(input, window_size, simple, squared);
}
template <typename I, typename T>
template <typename I, typename J>
inline
mln_ch_value(I, value::int_u8)
sauvola_threshold_dispatch(const value::rgb8&, const I& input,
unsigned window_size,
integral_image<T>& simple,
integral_image<T>& squared)
J& simple,
J& squared)
{
return impl::sauvola_threshold_rgb8(input, window_size,
simple, squared);
}
template <typename I, typename T>
template <typename I, typename J>
inline
mln_ch_value(I, value::int_u8)
sauvola_threshold_dispatch(