[Olena-patches] last-svn-commit-19-gc6851fd 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. --- scribo/ChangeLog | 20 + scribo/binarization/internal/first_pass_functor.hh | 141 +++ scribo/binarization/sauvola_threshold.hh | 318 +++--- scribo/canvas/integral_browsing.hh | 412 ++++++++ scribo/core/init_integral_image.hh | 96 ++ scribo/src/binarization/sauvola_ms.cc | 1072 ++++++++++++++++---- scribo/subsampling/integral_single_image.hh | 403 ++++++++ 7 files changed, 2103 insertions(+), 359 deletions(-) create mode 100644 scribo/binarization/internal/first_pass_functor.hh create mode 100644 scribo/canvas/integral_browsing.hh create mode 100644 scribo/core/init_integral_image.hh create mode 100644 scribo/subsampling/integral_single_image.hh diff --git a/scribo/ChangeLog b/scribo/ChangeLog index 4fed85c..9eac19c 100644 --- a/scribo/ChangeLog +++ b/scribo/ChangeLog @@ -1,3 +1,23 @@ +2009-12-04 Guillaume Lazzara &lt;z(a)lrde.epita.fr&gt; + + 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 &lt;z(a)lrde.epita.fr&gt; Add a new example in Scribo. diff --git a/scribo/binarization/internal/first_pass_functor.hh b/scribo/binarization/internal/first_pass_functor.hh new file mode 100644 index 0000000..13fc62a --- /dev/null +++ b/scribo/binarization/internal/first_pass_functor.hh @@ -0,0 +1,141 @@ +// 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 diff --git a/scribo/binarization/sauvola_threshold.hh b/scribo/binarization/sauvola_threshold.hh index 9a9e439..857de1f 100644 --- a/scribo/binarization/sauvola_threshold.hh +++ b/scribo/binarization/sauvola_threshold.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(const mln_value(I)&, const I& input, unsigned window_size, - integral_image<T>& simple, - integral_image<T>& squared) + J& simple, + J& squared) { // No dispatch for this kind of value type. mlc_abort(I)::check(); @@ -427,11 +427,11 @@ namespace scribo - 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) { trace::entering("scribo::binarization::sauvola_threshold"); @@ -441,8 +441,9 @@ namespace scribo typedef mln_value(I) value_t; mln_ch_value(I, value::int_u8) output = internal::sauvola_threshold_dispatch(value_t(), exact(input), - window_size, simple, - squared); + window_size, + exact(simple), + exact(squared)); trace::exiting("scribo::text::ppm2pbm"); return output; @@ -454,7 +455,10 @@ namespace scribo mln_ch_value(I, value::int_u8) sauvola_threshold(const Image<I>& input, unsigned window_size) { - mln::integral_image<value::int_u8> simple, squared; + mln_ch_value(I, double) + simple = init_integral_image(input, scribo::internal::identity_), + squared = init_integral_image(input, scribo::internal::square_); + return sauvola_threshold(input, window_size, simple, squared); } diff --git a/scribo/canvas/integral_browsing.hh b/scribo/canvas/integral_browsing.hh new file mode 100644 index 0000000..814733f --- /dev/null +++ b/scribo/canvas/integral_browsing.hh @@ -0,0 +1,412 @@ +// 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_CANVAS_INTEGRAL_BROWSING_HH +# define SCRIBO_CANVAS_INTEGRAL_BROWSING_HH + +# include <mln/core/image/image2d.hh> +# include <mln/util/couple.hh> + +namespace scribo +{ + + namespace canvas + { + using namespace mln; + + + template <typename F> + void integral_browsing(const image2d<util::couple<double, double> >& ima, + unsigned step, + unsigned w, unsigned h, + F& functor); + + +# ifndef MLN_INCLUDE_ONLY + + + namespace internal + { + + inline + void compute_stats(// in + double sum, double sum_2, unsigned n, + // out + double& mean, double& stddev) + { + mean = sum / n; + stddev = std::sqrt(sum_2 / n - mean * mean); + // unbias version: + // stddev = std::sqrt((sum_2 - n * mean * mean) / (n - 1)); + } + + } // end of namespace scribo::canvas::internal + + + + + template <typename F> + void integral_browsing(const image2d<util::couple<double, double> >& ima, + unsigned step, + unsigned w, unsigned h, + F& functor) + { + typedef util::couple<double, double> V; + typedef const V* Ptr; + Ptr a_ima, b_ima, c_ima, d_ima; + + const unsigned + nrows = ima.nrows(), + ncols = ima.ncols(), + row_0 = step / 2, + col_0 = step / 2; + + const unsigned + offset_down = ima.delta_index(dpoint2d(step, 0)), + offset_ante = ima.delta_index(dpoint2d(0, -w)), + offset_below = ima.delta_index(dpoint2d(+h, 0)); + + const unsigned + max_row_top = h/2, + max_row_mid = nrows - 1 - h/2, + max_col_left = w/2, + max_col_mid = ncols - 1 - w/2, + + step_2 = step * step, + + h_top = row_0 + h/2 + 1, + w_left = col_0 + w/2 + 1; + + unsigned row, col; + + for (col = col_0; col <= max_col_mid; col += step) ; + unsigned w_right = ncols - col + w/2; + + Ptr + d_tl_start, d_tr_start, + b_ml_start = 0, d_ml_start = 0, b_mr_start = 0, d_mr_start = 0, + b_bl_start = 0, d_bl_start = 0, b_br_start = 0, d_br = 0; + + double mean, stddev; + + + // ------------------------------- + // T (top) + // ------------------------------- + + const unsigned + delta_start_left = step * w_left, + delta_start_right = step * w_right, + step_w = step * w; + unsigned + size_tl_start = h_top * w_left, + size_tl, + delta_size_tl = h_top * step, + size_tc = h_top * w, + delta_size_tr = h_top * step, + size_tr_start = h_top * w_right, + size_tr; + + d_tl_start = & ima.at_(row_0 + h/2, col_0 + w/2); + d_tr_start = & ima.at_(row_0 + h/2, ncols - 1); + + for (row = row_0; row <= max_row_top; row += step) + { + + // TL (top left) + + d_ima = d_tl_start; + size_tl = size_tl_start; + + for (col = col_0; col <= max_col_left; col += step) + { + // D + internal::compute_stats(d_ima->first(), + d_ima->second(), + size_tl, + mean, stddev); + functor.exec(mean, stddev); + d_ima += step; + size_tl += delta_size_tl; + } + + delta_size_tl += step_2; + size_tl_start += delta_start_left; + d_tl_start += offset_down; + + // TC (top center) + + c_ima = d_ima + offset_ante; + + for (; col <= max_col_mid; col += step) + { + // D - C + internal::compute_stats(d_ima->first() - c_ima->first(), + d_ima->second() - c_ima->second(), + size_tc, + mean, stddev); + functor.exec(mean, stddev); + c_ima += step; + d_ima += step; + } + + size_tc += step_w; + + // TR (top right) + + d_ima = d_tr_start; + double + d_sum = d_ima->first(), + d_sum_2 = d_ima->second(); + size_tr = size_tr_start; + + for (; col < ncols; col += step) + { + // D* - C + internal::compute_stats(d_sum - c_ima->first(), + d_sum_2 - c_ima->second(), + size_tr, + mean, stddev); + functor.exec(mean, stddev); + c_ima += step; + size_tr -= delta_size_tr; + } + + delta_size_tr += step_2; + size_tr_start += delta_start_right; + d_tr_start += offset_down; + + } + + + + // ------------------------------- + // (M) middle + // ------------------------------- + + + const unsigned + size_ml_start = h * w_left, + h_step = h * step, + size_mc = w * h; + unsigned + size_ml, + size_mr_start = h * w_right, + size_mr; + + if (row <= max_row_mid) + { + b_ml_start = & ima.at_(row - h/2 - 1, col_0 + w/2); + d_ml_start = b_ml_start + offset_below; + b_mr_start = & ima.at_(row - h/2 - 1, ncols - 1); + d_mr_start = b_mr_start + offset_below; + } + + for (; row <= max_row_mid; row += step) + { + + // ML (middle left) + + size_ml = size_ml_start; + b_ima = b_ml_start; + d_ima = d_ml_start; + + for (col = col_0; col <= max_col_left; col += step) + { + // D - B + internal::compute_stats(d_ima->first() - b_ima->first(), + d_ima->second() - b_ima->second(), + size_ml, + mean, stddev); + functor.exec(mean, stddev); + b_ima += step; + d_ima += step; + size_ml += h_step; + } + + b_ml_start += offset_down; + d_ml_start += offset_down; + + // MC (middle center) + + a_ima = b_ima + offset_ante; + c_ima = d_ima + offset_ante; + + for (; col <= max_col_mid; col += step) + { + // D + A - B - C + internal::compute_stats((d_ima->first() - b_ima->first()) + (a_ima->first() - c_ima->first()), + (d_ima->second() - b_ima->second()) + (a_ima->second() - c_ima->second()), + size_mc, + mean, stddev); + functor.exec(mean, stddev); + a_ima += step; + b_ima += step; + c_ima += step; + d_ima += step; + } + + // MR (middle right) + + size_mr = size_mr_start; + b_ima = b_mr_start; + d_ima = d_mr_start; + double + d_b_sum = d_ima->first() - b_ima->first(), + d_b_sum_2 = d_ima->second() - b_ima->second(); + + for (; col < ncols; col += step) + { + // D* + A - B* - C + internal::compute_stats(d_b_sum + (a_ima->first() - c_ima->first()), + d_b_sum_2 + (a_ima->second() - c_ima->second()), + size_mr, + mean, stddev); + functor.exec(mean, stddev); + a_ima += step; + c_ima += step; + size_mr -= h_step; + } + + b_mr_start += offset_down; + d_mr_start += offset_down; + } + + + + // ------------------------------- + // B (bottom) + // ------------------------------- + + unsigned + size_bl_start = (nrows - row + h/2) * w_left, + size_bl, + delta_size_bl = (nrows - row + h/2) * step, + size_bc = (nrows - row + h/2) * w, + size_br_start = (nrows - row + h/2) * w_right, + delta_size_br = (nrows - row + h/2) * step, + size_br; + + if (row < nrows) + { + b_bl_start = & ima.at_(row - h/2 - 1, col_0 + w/2); + d_bl_start = & ima.at_(nrows - 1, col_0 + w/2); + b_br_start = & ima.at_(row - h/2 - 1, ncols - 1); + d_br = & ima.at_(nrows - 1, ncols - 1); + } + + for (; row < nrows; row += step) + { + + // BL (bottom left) + + size_bl = size_bl_start; + b_ima = b_bl_start; + d_ima = d_bl_start; + + for (col = col_0; col <= max_col_left; col += step) + { + // D* - B + internal::compute_stats(d_ima->first() - b_ima->first(), + d_ima->second() - b_ima->second(), + size_bl, + mean, stddev); + functor.exec(mean, stddev); + b_ima += step; + d_ima += step; + size_bl += delta_size_bl; + } + + delta_size_bl -= step_2; + size_bl_start -= delta_start_left; + b_bl_start += offset_down; + + // BC (bottom center) + + a_ima = b_ima + offset_ante; + c_ima = d_ima + offset_ante; + + for (; col <= max_col_mid; col += step) + { + // D* + A - B - C* + internal::compute_stats((d_ima->first() - b_ima->first()) + (a_ima->first() - c_ima->first()), + (d_ima->second() - b_ima->second()) + (a_ima->second() - c_ima->second()), + size_bc, + mean, stddev); +// std::cout << (d_ima->second() - b_ima->second()) + (a_ima->second() - c_ima->second()) << std::endl; + +// std::cout << d_ima->second() << " - " << b_ima->second() << " - " +// << a_ima->second() << " - " << c_ima->second() << std::endl; +// std::cout << d_ima->first() << " - " << b_ima->first() << " - " +// << a_ima->first() << " - " << c_ima->first() << std::endl; + functor.exec(mean, stddev); + a_ima += step; + b_ima += step; + c_ima += step; + d_ima += step; + } + + size_bc -= step_w; + + // BR (bottom right) + + size_br = size_br_start; + b_ima = b_br_start; + d_ima = d_br; + double + d_b_sum = d_ima->first() - b_ima->first(), + d_b_sum_2 = d_ima->second() - b_ima->second(); + + for (; col < ncols; col += step) + { + // D* + A - B* - C* + internal::compute_stats(d_b_sum + (a_ima->first() - c_ima->first()), + d_b_sum_2 + (a_ima->second() - c_ima->second()), + size_br, + mean, stddev); + functor.exec(mean, stddev); + a_ima += step; + c_ima += step; + size_br -= delta_size_br; + } + + delta_size_br -= step_2; + size_br_start -= delta_start_right; + b_br_start += offset_down; + } + + functor.finalize(); + } + +# endif // ! MLN_INCLUDE_ONLY + + + } // end of namespace scribo::canvas + +} // end of namespace mln + + +#endif // ! SCRIBO_CANVAS_INTEGRAL_BROWSING_HH diff --git a/scribo/core/init_integral_image.hh b/scribo/core/init_integral_image.hh new file mode 100644 index 0000000..ec19070 --- /dev/null +++ b/scribo/core/init_integral_image.hh @@ -0,0 +1,96 @@ +// 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_CORE_INIT_INTEGRAL_IMAGE_HH +# define SCRIBO_CORE_INIT_INTEGRAL_IMAGE_HH + +/// \file +/// +/// + +# include <mln/core/image/image2d.hh> + +namespace scribo +{ + using namespace mln; + + + namespace internal + { + + double square_(const double& val) + { + double v = static_cast<double>(val); + return v * v; + } + + double identity_(const double& val) + { + return static_cast<double>(val); + } + + } // end of namespace scribo::internal + + + template <typename I, typename F> + mln_ch_value(I,double) + init_integral_image(const Image<I>& input_, F& func) + { + trace::entering("scribo::init_integral_image"); + + const I& input = exact(input_); + mln_precondition(input.is_valid()); + mln_precondition(input.domain().pmin() == literal::origin); + + mln_ch_value(I,double) output; + initialize(output, input); + + unsigned + nrows_ = input.nrows(), + ncols_ = input.ncols(); + + output.at_(0,0) = func(input.at_(0, 0)); + + for (unsigned row = 1; row < nrows_; ++row) + output.at_(row,0) = output.at_(row - 1, 0) + func(input.at_(row, 0)); + + for (unsigned col = 1; col < ncols_; ++col) + output.at_(0,col) = output.at_(0, col - 1) + + func(input.at_(0, col)); + + for (unsigned row = 1; row < nrows_; ++row) + for (unsigned col = 1; col < ncols_; ++col) + output.at_(row, col) = output.at_(row - 1, col) + + output.at_(row, col - 1) + - output.at_(row - 1, col - 1) + + func(input.at_(row, col)); + + trace::exiting("scribo::init_integral_image"); + return output; + } + +} // end of namespace scribo + +#endif // ! SCRIBO_CORE_INIT_INTEGRAL_IMAGE_HH diff --git a/scribo/src/binarization/sauvola_ms.cc b/scribo/src/binarization/sauvola_ms.cc index 747bc11..8a07dcc 100644 --- a/scribo/src/binarization/sauvola_ms.cc +++ b/scribo/src/binarization/sauvola_ms.cc @@ -25,6 +25,11 @@ #include <mln/core/alias/neighb2d.hh> #include <mln/data/stretch.hh> +#include <mln/data/paste.hh> +// #include <mln/debug/iota.hh> +// #include <mln/debug/quiet.hh> +// #include <mln/debug/println.hh> +// #include <mln/debug/println_with_border.hh> #include <mln/debug/filename.hh> #include <mln/fun/i2v/array.hh> #include <mln/io/pbm/all.hh> @@ -32,17 +37,30 @@ #include <mln/io/ppm/all.hh> #include <mln/literal/colors.hh> #include <mln/math/sqr.hh> -#include <mln/subsampling/subsampling.hh> +#include <mln/math/abs.hh> + +#include <mln/subsampling/antialiased.hh> + #include <mln/transform/influence_zone_geodesic.hh> #include <mln/util/timer.hh> #include <mln/value/int_u16.hh> #include <mln/value/int_u8.hh> #include <mln/value/label_16.hh> #include <mln/value/rgb8.hh> +#include <mln/border/equalize.hh> +#include <mln/border/mirror.hh> +#include <mln/border/adjust.hh> #include <mln/debug/filename.hh> +#include <mln/core/box_runend_piter.hh> +#include <mln/core/box_runstart_piter.hh> + +#include <scribo/subsampling/integral_single_image.hh> +//#include <scribo/subsampling/integral.hh> + #include <scribo/core/macros.hh> #include <scribo/core/object_image.hh> + #include <scribo/filter/objects_small.hh> #include <scribo/filter/objects_large.hh> @@ -50,155 +68,737 @@ #include <scribo/filter/objects_thick.hh> #include <scribo/primitive/extract/objects.hh> + #include <scribo/binarization/sauvola_threshold.hh> +#include <scribo/binarization/internal/first_pass_functor.hh> + +#include <scribo/canvas/integral_browsing.hh> + #include <scribo/debug/usage.hh> #include <scribo/debug/save_object_diff.hh> + namespace mln { + using value::int_u8; -// FIXME: do not use scribo's filters since they use object images. -// Object images computes object bounding boxes which are not needed here. -// -// template <typename L> -// mln_concrete(L) -// filter_small_objects(const Image<L>& objects, -// const mln_value(L)& nobjects, unsigned min_size) -// { -// card = labeling::compute(card_t(), objects, objects.nlabels()); - -// mln_concrete(L) output; -// initialize(output, objects); -// fun::i2v::array<bool> f(nobjects.next(), 0); -// for (unsigned i = 0; i < card.size(); ++i) + 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, typename T> - mln_ch_value(I, bool) - apply_sauvola(const I& intensity, const T& threshold, unsigned s) + image2d<int_u8> + compute_t_n_and_e_2(const image2d<int_u8>& sub, image2d<int_u8>& e_2, + unsigned lambda_min, unsigned lambda_max, + unsigned q, unsigned i, unsigned w, + const image2d<util::couple<double,double> >& integral_sum_sum_2) + // lambdas: limits of component cardinality at this scale { - mln_precondition(intensity.domain() == threshold.domain()); + typedef image2d<int_u8> I; + typedef point2d P; - mln_ch_value(I, bool) output; - initialize(output, intensity); + util::timer tt; + float t_; - mln_piter(I) p(intensity.domain()); - for_all(p) - output(p) = (intensity(p) <= threshold(p / s)); + unsigned ratio = std::pow(q, i - 2); // Ratio in comparison to e_2 + + + tt.restart(); + + unsigned w_2 = w / 2 * ratio; + + // 1st pass + scribo::binarization::internal::first_pass_functor< image2d<int_u8> > + f(sub); + scribo::canvas::integral_browsing(integral_sum_sum_2, + ratio, + w_2, w_2, + f); + +// debug::println("mask", f.msk); +// debug::println("parent", f.parent); +// debug::println("card", f.card); + + t_ = tt; + if (mln::debug::quiet) + std::cout << "1st pass - " << t_ << std::endl; + + tt.restart(); + + { + util::array<mln_value_(I) *> ptr(ratio); + unsigned nrows = geom::nrows(e_2); + + mln_box_runend_piter_(I) sp(sub.domain()); // Backward. + unsigned ncols = sp.run_length(); + for_all(sp) + { + unsigned p = &sub(sp) - sub.buffer(); // Offset + P site = sp; + + { + P tmp = site * ratio; + + // FIXME: to be removed! + if (tmp.row() + ratio >= nrows) + ptr.resize(nrows - tmp.row()); + + ptr(0) = &e_2(tmp); + // FIXME: pointers could just be updated with an offset. + for (unsigned j = 1; j < ptr.size(); ++j) + { + tmp[0] += 1; + ptr(j) = & e_2(tmp); + } + } + + for (unsigned j = 0; j < ncols; ++j) + { + if (f.msk.element(p)) + { + + mln_site_(I) sq = site * ratio; + + if (f.parent.element(p) == p) + { + // test over the component cardinality + f.msk.element(p) = f.card.element(p) > lambda_min + && f.card.element(p) < lambda_max; + + if (f.msk.element(p) && e_2(sq) == 0u) + { + for (unsigned l = 0; l < ptr.size(); ++l) + std::memset(ptr(l), i, ratio * sizeof(mln_value_(I))); +// debug(sq) = i; + } + + } + else + { + // Propagation + f.msk.element(p) = f.msk.element(f.parent.element(p)); + + if (f.msk.element(p) && e_2(sq) == 0u) + { + for (unsigned l = 0; l < ptr.size(); ++l) + std::memset(ptr(l), i, ratio * sizeof(mln_value_(I))); +// debug(sq) = i; + } + + } + } + + for (unsigned l = 0; l < ptr.size(); ++l) + ptr(l) -= ratio; + + --site[1]; + --p; + } + + } - return output; + t_ = tt; + if (mln::debug::quiet) + std::cout << "2nd pass - " << t_ << std::endl; + +// io::pgm::save(e_2, mln::debug::filename("e.pgm", i)); +// io::pgm::save(debug, mln::debug::filename("debug.pgm", i)); +// debug::println(msk); +// io::pbm::save(f.msk, mln::debug::filename("mask.pbm", i)); +// io::pgm::save(data::stretch(int_u8(), card), mln::debug::filename("card.pgm")); + } // end of 2nd pass + +// io::pgm::save(f.t_sub, mln::debug::filename("t.pgm", i)); + return f.t_sub; } - template <typename I> - mln_concrete(I) - enlarge(const Image<I>& input_, unsigned ratio) + + template <typename I, typename J, typename K> + mln_ch_value(I, bool) + binarize_generic(const I& in, const J& e2, const util::array<K>& t_ima, + unsigned s) { - const I& input = exact(input_); - mln_precondition(input.is_valid()); + mln_ch_value(I,bool) out; + initialize(out, in); - mln_domain(I) bbox(input.domain().pmin() * ratio, - input.domain().pmax() * ratio - + (ratio - 1) * mln::down_right); + typedef const mln_value(K)* ptr_type; - mln_concrete(I) output(bbox); + ptr_type ptr_t[5]; + ptr_t[2] = & t_ima[2].at_(0, 0); + ptr_t[3] = & t_ima[3].at_(0, 0); + ptr_t[4] = & t_ima[4].at_(0, 0); - mln_site(I) first_p = input.domain().pmin(); - for (def::coord j = geom::min_col(input); - j <= geom::max_col(input); ++j) - for (def::coord i = geom::min_row(input); - i <= geom::max_row(input); ++i) - { - point2d p1(i, j); - point2d p2(first_p[0] + i * ratio, first_p[1] + j * ratio); - output(p2) = input(p1); - output(p2 + mln::right) = input(p1); - output(p2 + mln::down_right) = input(p1); - output(p2 + mln::down) = input(p1); + const mln_value(J)* ptr_e2 = & e2.at_(0, 0); + const mln_value(I)* ptr__in = & in.at_(0, 0); + bool* ptr__out = & out.at_(0, 0); + + + // Since we iterate from a smaller image in the largest ones and + // image at scale 1 does not always have a size which can be + // divided by (4*s), some sites in the border may not be processed + // and we must skip them. + unsigned more_offset = in.border() - ((4 * s) - in.ncols() % (4 * s)); + + if (more_offset == (4 * s)) + more_offset = 0; // No offset needed. + + const int + nrows4 = t_ima[4].nrows(), ncols4 = t_ima[4].ncols(), + + + delta1 = in.delta_index(dpoint2d(+1, -(s - 1))), + delta1b = in.delta_index(dpoint2d(+1, -(s + s - 1))), + delta1c = in.delta_index(dpoint2d(-(s + s - 1), +1)), + delta1d = in.delta_index(dpoint2d(+1, -(s * 4 - 1))), + delta1e = in.delta_index(dpoint2d(-(s * 4 - 1), +1)), + delta1f = in.delta_index(dpoint2d(-(s - 1), +1)), + + delta2 = t_ima[2].delta_index(dpoint2d(+1, -1)), + delta2b = t_ima[2].delta_index(dpoint2d(+1, -3)), + delta2c = t_ima[2].delta_index(dpoint2d(-3, +1)), + + delta3 = t_ima[3].delta_index(dpoint2d(+1, -1)), + + eor1 = in.delta_index(dpoint2d(+4 * s, - in.ncols() - in.border())) + more_offset, + eor2 = t_ima[2].delta_index(dpoint2d(+4,- t_ima[2].ncols())), + eor3 = t_ima[3].delta_index(dpoint2d(+2,- t_ima[3].ncols())), + eor4 = t_ima[4].delta_index(dpoint2d(+1,- t_ima[4].ncols())); + + mln_value(J) threshold; + for (int row4 = 0; row4 < nrows4; ++row4) + { + for (int col4 = 0; col4 < ncols4; ++col4) + { + // top left 1 + { + threshold = *ptr_t[*ptr_e2]; + { + for (unsigned i = 1; i < s; ++i) + { + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + + *ptr__out = *ptr__in < threshold; + ptr__out += delta1; ptr__in += delta1; + } + + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + *ptr__out = *ptr__in < threshold; + ptr__out += delta1f; ptr__in += delta1f; + } + + ++ptr_t[2]; ++ptr_e2; + threshold = *ptr_t[*ptr_e2]; + { + for (unsigned i = 1; i < s; ++i) + { + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + + *ptr__out = *ptr__in < threshold; + ptr__out += delta1; ptr__in += delta1; + } + + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + *ptr__out = *ptr__in < threshold; + ptr__out += delta1b; ptr__in += delta1b; + } + + ptr_t[2] += delta2; ptr_e2 += delta2; + threshold = *ptr_t[*ptr_e2]; + { + for (unsigned i = 1; i < s; ++i) + { + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + + *ptr__out = *ptr__in < threshold; + ptr__out += delta1; ptr__in += delta1; + } + + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + *ptr__out = *ptr__in < threshold; + ptr__out += delta1f; ptr__in += delta1f; + + } + + ++ptr_t[2]; ++ptr_e2; + threshold = *ptr_t[*ptr_e2]; + { + for (unsigned i = 1; i < s; ++i) + { + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + + *ptr__out = *ptr__in < threshold; + ptr__out += delta1; ptr__in += delta1; + } + + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + *ptr__out = *ptr__in < threshold; + ptr__out += delta1c; ptr__in += delta1c; + } + + ptr_t[2] -= delta2; ptr_e2 -= delta2; + } + + // top right 1 + ptr_t[3] += 1; + { + threshold = *ptr_t[*ptr_e2]; + { + for (unsigned i = 1; i < s; ++i) + { + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + + *ptr__out = *ptr__in < threshold; + ptr__out += delta1; ptr__in += delta1; + } + + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + *ptr__out = *ptr__in < threshold; + ptr__out += delta1f; ptr__in += delta1f; + } + + ++ptr_t[2]; ++ptr_e2; + threshold = *ptr_t[*ptr_e2]; + { + for (unsigned i = 1; i < s; ++i) + { + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + + *ptr__out = *ptr__in < threshold; + ptr__out += delta1; ptr__in += delta1; + } + + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + *ptr__out = *ptr__in < threshold; + ptr__out += delta1b; ptr__in += delta1b; + } + + ptr_t[2] += delta2; ptr_e2 += delta2; + threshold = *ptr_t[*ptr_e2]; + { + for (unsigned i = 1; i < s; ++i) + { + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + + *ptr__out = *ptr__in < threshold; + ptr__out += delta1; ptr__in += delta1; + } + + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + *ptr__out = *ptr__in < threshold; + ptr__out += delta1f; ptr__in += delta1f; + } + + ++ptr_t[2]; ++ptr_e2; + threshold = *ptr_t[*ptr_e2]; + { + for (unsigned i = 1; i < s; ++i) + { + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + + *ptr__out = *ptr__in < threshold; + ptr__out += delta1; ptr__in += delta1; + } + + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + *ptr__out = *ptr__in < threshold; + ptr__out += delta1d; ptr__in += delta1d; + } + + ptr_t[2] += delta2b; ptr_e2 += delta2b; + } + + // bot left 1 + ptr_t[3] += delta3; + { + threshold = *ptr_t[*ptr_e2]; + { + for (unsigned i = 1; i < s; ++i) + { + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + + *ptr__out = *ptr__in < threshold; + ptr__out += delta1; ptr__in += delta1; + } + + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + *ptr__out = *ptr__in < threshold; + ptr__out += delta1f; ptr__in += delta1f; + } + + ++ptr_t[2]; ++ptr_e2; + threshold = *ptr_t[*ptr_e2]; + { + for (unsigned i = 1; i < s; ++i) + { + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + + *ptr__out = *ptr__in < threshold; + ptr__out += delta1; ptr__in += delta1; + } + + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + *ptr__out = *ptr__in < threshold; + ptr__out += delta1b; ptr__in += delta1b; + } + + ptr_t[2] += delta2; ptr_e2 += delta2; + threshold = *ptr_t[*ptr_e2]; + { + for (unsigned i = 1; i < s; ++i) + { + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + + *ptr__out = *ptr__in < threshold; + ptr__out += delta1; ptr__in += delta1; + } + + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + *ptr__out = *ptr__in < threshold; + ptr__out += delta1f; ptr__in += delta1f; + } + + ++ptr_t[2]; ++ptr_e2; + threshold = *ptr_t[*ptr_e2]; + { + for (unsigned i = 1; i < s; ++i) + { + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + + *ptr__out = *ptr__in < threshold; + ptr__out += delta1; ptr__in += delta1; + } + + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + *ptr__out = *ptr__in < threshold; + ptr__out += delta1c; ptr__in += delta1c; + } + + ptr_t[2] -= delta2; ptr_e2 -= delta2; + } + + // bot right 1 + ptr_t[3] += 1; + { + threshold = *ptr_t[*ptr_e2]; + { + for (unsigned i = 1; i < s; ++i) + { + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + + *ptr__out = *ptr__in < threshold; + ptr__out += delta1; ptr__in += delta1; + } + + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + *ptr__out = *ptr__in < threshold; + ptr__out += delta1f; ptr__in += delta1f; + } + + ++ptr_t[2]; ++ptr_e2; + threshold = *ptr_t[*ptr_e2]; + { + for (unsigned i = 1; i < s; ++i) + { + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + + *ptr__out = *ptr__in < threshold; + ptr__out += delta1; ptr__in += delta1; + } + + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + *ptr__out = *ptr__in < threshold; + ptr__out += delta1b; ptr__in += delta1b; + } + + ptr_t[2] += delta2; ptr_e2 += delta2; + threshold = *ptr_t[*ptr_e2]; + { + for (unsigned i = 1; i < s; ++i) + { + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + + *ptr__out = *ptr__in < threshold; + ptr__out += delta1; ptr__in += delta1; + } + + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + *ptr__out = *ptr__in < threshold; + ptr__out += delta1f; ptr__in += delta1f; + } + + ++ptr_t[2]; ++ptr_e2; + threshold = *ptr_t[*ptr_e2]; + { + for (unsigned i = 1; i < s; ++i) + { + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + + *ptr__out = *ptr__in < threshold; + ptr__out += delta1; ptr__in += delta1; + } + + for (unsigned j = 1; j < s; ++j) + { + *ptr__out = *ptr__in < threshold; + ++ptr__out; ++ptr__in; + } + *ptr__out = *ptr__in < threshold; + ptr__out += delta1e; ptr__in += delta1e; + } + } + + // bot right -> next top left + ptr_t[2] += delta2c; ptr_e2 += delta2c; + ptr_t[3] = ptr_t[3] - delta3; + ptr_t[4] += 1; } - return output; + // eof -> next bof + ptr__out += eor1; ptr__in += eor1; + ptr_t[2] += eor2; ptr_e2 += eor2; + ptr_t[3] += eor3; + ptr_t[4] += eor4; + } + + +// mln::debug::println(out); + + return out; + } + + + + + unsigned sub(unsigned nbr, unsigned down_scaling) + { + return (nbr + down_scaling - 1) / down_scaling; } template <typename I> - object_image(mln_ch_value(I, value::label_16)) - get_objects(const Image<I>& input_, - unsigned min_size, unsigned max_size, - unsigned ratio, unsigned w, unsigned scale) + util::array<util::couple<mln_domain(I), unsigned> > + compute_sub_domains(const I& ima, unsigned n_scales, unsigned s) { - const I& input = exact(input_); + util::array<util::couple<unsigned, unsigned> > n(n_scales + 2); - mln_precondition(input.is_valid()); + n(1) = make::couple(ima.nrows(), ima.ncols()); + n(2) = make::couple(sub(n(1).first(), s), + sub(n(1).second(), s)); + for (unsigned i = 3; i <= n_scales + 1; ++i) + n(i) = make::couple(sub(n(i - 1).first(), 2), + sub(n(i - 1).second(), 2)); - typedef mln_ch_value(I, value::label_16) L; - dpoint2d none(0, 0); - I input_sub = mln::subsampling::subsampling(input, none, ratio); + util::array<util::couple<mln_domain(I), unsigned> > out(n.size()); + out(0) = make::couple(make::box2d(1,1), 1u); + out(1) = make::couple(make::box2d(ima.nrows(), ima.ncols()), 2u); + out(n_scales + 1) = make::couple(make::box2d(n(n_scales + 1).first(), + n(n_scales + 1).second()), + 1u); - image2d<value::int_u8> - t = scribo::binarization::sauvola_threshold(input_sub, w); + for (unsigned i = n_scales; i > 1; --i) + out(i) = make::couple(make::box2d(2 * out(i + 1).first().nrows(), + 2 * out(i + 1).first().ncols()), + 2 * out(i + 1).second()); - image2d<bool> b = apply_sauvola(input_sub, t, 1); + out(1).second() = std::max(out(2).first().ncols() * s - ima.ncols(), + out(2).first().nrows() * s - ima.nrows()); - value::label_16 nb; - object_image(L) lbl = scribo::primitive::extract::objects(b, c8(), nb); - object_image(L) lbl_raw(lbl); +// out(1).second() = std::max(ima.ncols() % (4 * s), +// ima.nrows() % (4 * s)); - if (min_size > 0) - lbl = scribo::filter::objects_small(lbl, min_size); - if (max_size > 0) - lbl = scribo::filter::objects_large(lbl, max_size); - scribo::debug::save_object_diff(lbl_raw, lbl, - mln::debug::filename("filter_diff.ppm", scale + 2)); +// out(2).first().ncols() * s - ima.ncols(), +// out(2).first().nrows() * s - ima.nrows() ); - return lbl; + return out; } + bool check_args(int argc, char * argv[]) { - if (argc < 7) + if (argc < 5 || argc > 7) return false; - int nb_scale = atoi(argv[3]); - int s = atoi(argv[4]); - int q = atoi(argv[5]); +// int nb_scale = atoi(argv[3]); + int s = atoi(argv[3]); +// int q = atoi(argv[5]); - if (q < 2) - { - std::cout << "q must be greater than 2." << std::endl; - return false; - } - if (s < 1 || s < q) +// if (q < 2) +// { +// std::cout << "q must be greater than 2." << std::endl; +// return false; +// } + if (s < 1 || s > 3)// || s < q) { - std::cout << "s must be greater or equal to 1 and greater than q." + std::cout << "s must be set to 2 or 3." << std::endl; return false; } - if (nb_scale < 1) - { - std::cout << "Not enough scales." << std::endl; - return false; - } - if ((argc - 7) != (nb_scale - 1)) - { - std::cout << "Not enough area threshold." - << "There must be nb_scale - 1 thresholds." - << std::endl; - return false; - } +// if (nb_scale < 1) +// { +// std::cout << "Not enough scales." << std::endl; +// return false; +// } + +// if ((argc - 7) != (nb_scale - 1)) +// { +// std::cout << "Not enough area threshold." +// << "There must be nb_scale - 1 thresholds." +// << std::endl; +// return false; +// } return true; } + void data_rand(image2d<unsigned>& e2) + { + unsigned v = 2; + mln_piter_(box2d) p(e2.domain()); + for_all(p) + { + e2(p) = v++; + if (v == 5) v = 2; + } + } + + } // end of namespace mln; @@ -207,11 +807,10 @@ namespace mln const char *args_desc[][2] = { { "input.pgm", "A graylevel image." }, - { "w", "Window size." }, - { "nb_scale", "Number of scales (Common value: 3)." }, + { "w", "Window size. (Common value: 51)" }, { "s", "First subsampling ratio (Common value: 2)." }, - { "q", "Next subsampling ratio (Common value: 2)." }, - { "<area threshold> (Common values: 200 800)", "Area threshold" }, + { "min_area", "Minimum object area (relative to scale 2) (Common value: 200)" }, + { "debug", "Display debug/bench data if set to 1" }, {0, 0} }; @@ -229,11 +828,10 @@ int main(int argc, char *argv[]) typedef image2d<label_16> L; - if (!check_args(argc, argv)) return scribo::debug::usage(argv, - "Binarization of a color image based on Sauvola's algorithm.", - "input.pgm w nb_scale s q <area thresholds>", + "Multi-Scale Binarization of a color image based on Sauvola's algorithm.", + "input.pgm w s area_thresholds output.pbm [debug]", args_desc, "A binary image."); trace::entering("main"); @@ -244,163 +842,233 @@ int main(int argc, char *argv[]) unsigned w = atoi(argv[2]); // First subsampling scale. - unsigned s = atoi(argv[4]); + unsigned s = atoi(argv[3]); // Number of subscales. - unsigned nb_subscale = atoi(argv[3]); + unsigned nb_subscale = 3;//atoi(argv[3]); // Subscale step. - unsigned q = atoi(argv[5]); + unsigned q = 2;//atoi(argv[5]); + + mln::debug::quiet = true; - std::cout << "Running Sauvola_ms with w = " << w - << ", s = " << s - << ", nb_subscale = " << nb_subscale - << ", q = " << q - << std::endl; + if (argc == 7) + mln::debug::quiet = ! atoi(argv[6]); + + + if (mln::debug::quiet) + std::cout << "Running Sauvola_ms with w = " << w + << ", s = " << s + << ", nb_subscale = " << nb_subscale + << ", q = " << q + << std::endl; typedef image2d<value::int_u8> I; dpoint2d none(0, 0); - mln::util::timer timer_; + mln::util::timer + timer_, + sauvola_timer_;; + + // Tmp variable used for timer results; + float t_; - timer_.start(); I input_full; io::pgm::load(input_full, argv[1]); - std::cout << "Image loaded - " << timer_ << std::endl; - timer_.restart(); + { + unsigned max_dim = math::max(input_full.ncols() / 2, + input_full.nrows() / 2); + if (w > max_dim) + { + std::cout << "------------------" << std::endl; + std::cout << "The window is too large! Image size is only " + << input_full.nrows() << "x" << input_full.ncols() + << std::endl + << "Window size must not exceed " << max_dim + << std::endl; + return 1; + } + } - util::array<object_image(L)> lbls; +// I input_full(9,9); +// mln::debug::iota(input_full); + sauvola_timer_.start(); + unsigned lambda_min = atoi(argv[4]); + unsigned lambda_max = lambda_min * q; // * atoi(argv[7]) - // First subsampling 1/s -// std::cout << "1/" << s << std::endl; - timer_.start(); - I input = mln::subsampling::subsampling(input_full, none, s); + util::array<I> t_ima; + + // Make sure t_ima indexes start from 2. + { + I dummy(1,1); + for (unsigned i = 0; i < nb_subscale + 2; ++i) + t_ima.append(dummy); + } - integral_image<int_u8> simple, squared; - image2d<int_u8> - t_1 = scribo::binarization::sauvola_threshold(input, w, simple, squared); - label_16 nb1; - image2d<bool> b_1 = apply_sauvola(input, t_1, 1); + image2d<int_u8> e_2; + util::array<I> sub_ima; + + // Make sure sub_ima indexes start from 2. { - lbls.append(primitive::extract::objects(b_1, c8(), nb1)); - object_image(L) lbl_raw(lbls[0]); - lbls[0] = filter::objects_large(lbls[0], atoi(argv[6])); - scribo::debug::save_object_diff(lbl_raw, lbls[0], - mln::debug::filename("filter_diff.ppm", 2)); + I dummy(1,1); + sub_ima.append(dummy); + sub_ima.append(dummy); } - std::cout << "Scale 2 - 1/" << s << " Done - " << timer_ << std::endl; + timer_.restart(); + util::array<util::couple<box2d, unsigned> > + sub_domains = compute_sub_domains(input_full, nb_subscale, s); + if (mln::debug::quiet) + std::cout << "adjusting input border to " << sub_domains(1).second() + << std::endl; - // Additional subscales. - for (unsigned i = 1; i < nb_subscale; ++i) - { - unsigned ratio = std::pow(q, i); -// std::cout << "Scale " << 2 + i << " - 1/" << s * ratio << std::endl; - timer_.restart(); - unsigned - min_size = atoi(argv[5 + i]) / ratio, - max_size; - if ((3 + i) == static_cast<unsigned>(argc)) - max_size = 0; // Last subscale, so not max size limit. - else - max_size = atoi(argv[6 + i]) / ratio; - - lbls.append(get_objects(input, min_size, max_size, ratio, w, i)); - std::cout << "Scale " << 2 + i - << " - 1/" << s * ratio - << " Done - " << timer_ << std::endl; - } + border::adjust(input_full, sub_domains(1).second()); + border::mirror(input_full); - std::cout << "--------" << std::endl; +// mln::debug::println_with_border(input_full); - // Constructing "scale image". - image2d<int_u8> e; - initialize(e, input); + t_ = timer_; + if (mln::debug::quiet) + std::cout << "sub domains computed and adjust input border size - " + << t_ << std::endl; - data::fill(e, 0); + // Resize input and compute integral images. + timer_.restart(); + image2d<util::couple<double,double> > integral_sum_sum_2; - typedef object_image(L) obj_t; + if (mln::debug::quiet) + std::cout << "sub_domain(2).domain() == " << sub_domains(2).first() << std::endl; - for (int i = nb_subscale - 1; i >= 0; --i) - { - unsigned ratio = std::pow(q, i); + sub_ima.append(scribo::subsampling::integral(input_full, s, + integral_sum_sum_2, + sub_domains[2].first(), + sub_domains[2].second())); - std::cout << "Scale " << 2 + i << " - 1/" << s * ratio << " merged" << std::endl; +// mln::debug::println(integral_sum_sum_2); - mln_piter_(obj_t) p(lbls[i].domain()); - for_all(p) - if (lbls[i](p) != 0 && e(p * ratio) == 0) - { - box2d b(p * ratio, p * ratio + mln::down_right * (ratio - 1)); - data::fill((e | b).rw(), i + 1); - } - } + t_ = timer_; + if (mln::debug::quiet) + std::cout << "subsampling 1 -> 2 And integral images - " << t_ + << " - nsites = " + << input_full.domain().nsites() << " -> " + << sub_ima[2].domain().nsites() << " - " + << input_full.domain() << " -> " + << sub_ima[2].domain() << std::endl; - std::cout << "--------" << std::endl; - /// Saving "scale image". + for (unsigned i = 3; i <= nb_subscale + 1; ++i) { - mln::fun::i2v::array<value::int_u8> f(nb_subscale + 1, 0); - for (unsigned i = 1; i <= nb_subscale; ++i) - f(i) = 255 / nb_subscale * i; - - io::pgm::save(e, "e_raw.pgm"); - io::pgm::save(data::transform(e, f), "e.pgm"); + timer_.restart(); + sub_ima.append(mln::subsampling::antialiased(sub_ima[i - 1], q, none, + sub_domains[i].first(), + sub_domains[i].second())); + t_ = timer_; + if (mln::debug::quiet) + std::cout << "subsampling " << (i - 1) << " -> " << i + << " - " << t_ + << " - nsites = " + << sub_ima[i].domain().nsites() << " - " + << sub_ima[i].domain() + << std::endl; } - /// Saving influence zone scale image. - image2d<int_u8> - e_ext = transform::influence_zone_geodesic(e, c8(), mln_max(unsigned)); - io::pgm::save(e_ext, "e_ext.pgm"); + initialize(e_2, sub_ima[2]); + data::fill(e_2, 0u); + // Compute threshold image. - /// Creating window size image - std::cout << "Creating window size image..." << std::endl; - mln::fun::i2v::array<value::int_u16> f(nb_subscale, 0); - for (unsigned i = 0; i < nb_subscale; ++i) + // Highest scale -> no maximum component size. { - unsigned ratio = std::pow(q, i); - f(i) = ratio * w; + timer_.restart(); + int i = sub_ima.size() - 1; + unsigned ratio = std::pow(q, i - 2); // Ratio compared to e_2 + t_ima[i] = compute_t_n_and_e_2(sub_ima[i], e_2, + lambda_min / ratio, + mln_max(unsigned), + q, i, w, integral_sum_sum_2); + + t_ = timer_; + if (mln::debug::quiet) + std::cout << "Scale " << i + << " - 1/" << s * ratio + << " compute t_n and update e - " << t_ << std::endl; } - image2d<int_u16> wsize = data::transform(e_ext, f); - io::pgm::save(data::stretch(int_u8(), wsize), "wsize.pgm"); + // Other scales -> maximum and minimum component size. + { + for (int i = sub_ima.size() - 2; i > 2; --i) + { + timer_.restart(); + unsigned ratio = std::pow(q, i - 2); // Ratio compared to e_2 + t_ima[i] = compute_t_n_and_e_2(sub_ima[i], e_2, + lambda_min / ratio, + lambda_max / ratio, + q, i, w, integral_sum_sum_2); + + t_ = timer_; + if (mln::debug::quiet) + std::cout << "Scale " << i + << " - 1/" << s * ratio + << " compute t_n and update e - " << t_ << std::endl; + } + } + // Lowest scale -> no minimum component size. + { + timer_.restart(); + t_ima[2] = compute_t_n_and_e_2(sub_ima[2], e_2, 0, lambda_max, + 1, 2, w, integral_sum_sum_2); + t_ = timer_; + if (mln::debug::quiet) + std::cout << "Scale " << 2 + << " - 1/" << s + << " compute t_n and update e - " << t_ << std::endl; + + } + + if (mln::debug::quiet) + std::cout << "--------" << std::endl; +// io::pgm::save(e_2, mln::debug::filename("e.pgm")); - /// Constructing threshold image - image2d<int_u8> t; - initialize(t, input); - std::cout << "Computing threshold image..." << std::endl; timer_.restart(); - mln_piter_(L) p(input.domain()); - for_all(p) - t(p) = binarization::internal::compute_sauvola_threshold(p, - simple, - squared, - wsize(p)); - std::cout << "Compute t Done - " << timer_ << std::endl; + e_2 = transform::influence_zone_geodesic(e_2, c8()); + t_ = timer_; + if (mln::debug::quiet) + std::cout << "influence zone - " << t_ << std::endl; - io::pgm::save(t, "t.pgm"); - timer_.restart(); +// Saving influence zone scale image. +// io::pgm::save(e_2, mln::debug::filename("e_ext.pgm")); +// io::pbm::save(bin, argv[8]); +// io::pgm::save(t, mln::debug::filename("t.pgm")); - /// Applying threshold image and save. - io::pbm::save(apply_sauvola(input_full, t, s), argv[argc - 1]); - std::cout << "sauvola applied and saved Done - " << timer_ << std::endl; +// for (unsigned i = 2; i < t_ima.size(); ++i) +// io::pgm::save(t_ima[i], mln::debug::filename("t.pgm", i)); - trace::exiting("main"); + timer_.restart(); + image2d<bool> out_new = binarize_generic(input_full, e_2, t_ima, s); +// image2d<bool> out_new = binarize(input_full, e_2, t_ima); + t_ = timer_; + if (mln::debug::quiet) + std::cout << "Compute bin - " << t_ << std::endl; + + t_ = sauvola_timer_; + if (mln::debug::quiet) + std::cout << "Sauvola : " << t_ << std::endl; + io::pbm::save(out_new, argv[5]); } + diff --git a/scribo/subsampling/integral_single_image.hh b/scribo/subsampling/integral_single_image.hh new file mode 100644 index 0000000..07b391d --- /dev/null +++ b/scribo/subsampling/integral_single_image.hh @@ -0,0 +1,403 @@ +// 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_SUBSAMPLING_INTEGRAL_HH +# define SCRIBO_SUBSAMPLING_INTEGRAL_HH + +/// \file +/// +/// Both subsampling and integral image computation. + +#include <mln/core/concept/image.hh> +#include <mln/metal/equal.hh> +#include <mln/extension/fill.hh> +#include <mln/debug/println.hh> +#include <mln/debug/println_with_border.hh> + + + +namespace scribo +{ + + namespace subsampling + { + using namespace mln; + + + /*! \brief Subsample an image and compute tow integral images. + + \param[in] input An image of Scalar. + \param[in] scale The scale factor. + \param[in] integral_sum Integral image of mean values. + \param[in] integral_sum_2 Integral image of squared mean values. + \param[in] output_domain The domain of the subscaled image. + \param[in] border_thickness Border of the integral and subscaled images. + + \p integral_sum, \p integral_sum_2 and output image have the same domain. + The output domain is set with \p output_domain. + + */ + template <typename I, typename J> + inline + mln_concrete(I) + integral(const Image<I>& input, unsigned scale, + Image<J>& integral_sum, Image<J>& integral_sum_2, + const mln_domain(I)& output_domain, unsigned border_thickness); + + /*! \overload + + The output domain is defined as follows : + + for_each(dimension) + input.domain()[dimension] = (N + scale - 1) / scale + + where N is the number of elements in one run in a given dimension. + + The border thickness is set to mln::border::thickness. + */ + template <typename I, typename J> + inline + mln_concrete(I) + integral(const Image<I>& input, unsigned scale, + Image<J>& integral_sum, Image<J>& integral_sum_2); + + + +# ifndef MLN_INCLUDE_ONLY + + + // Implementation. + + namespace impl + { + + template <typename I, typename J> + inline + mln_concrete(I) + integral_3(const Image<I>& input_, + Image<J>& integral_sum_sum_2_, + const mln_domain(I)& output_domain, + unsigned border_thickness) + { + trace::entering("subsampling::impl::integral_3"); + + const unsigned scale = 3; + + const I& input = exact(input_); + J& integral_sum_sum_2 = exact(integral_sum_sum_2_); + const unsigned area = scale * scale; + + mln_precondition(input.is_valid()); + mln_precondition(input.domain().pmin() == literal::origin); + mln_precondition(scale > 1); + + typedef mln_value(I) V; + typedef mln_sum(V) S; + typedef mln_value(J) V2; + typedef mln_site(I) P; + + mln_concrete(I) sub(output_domain, border_thickness); + V* p_sub = sub.buffer(); + + integral_sum_sum_2.init_(output_domain, border_thickness); + V2* p_integ = integral_sum_sum_2.buffer(); + + const unsigned up = sub.delta_index(dpoint2d(-1, 0)); + + const unsigned nrows = 3 * output_domain.nrows(); + const unsigned ncols = 3 * output_domain.ncols(); + + unsigned row = 0; + + unsigned b_offset = sub.delta_index(dpoint2d(border_thickness, + border_thickness)); + p_sub += b_offset; + p_integ += b_offset; + { + S h_sum = 0, h_sum_2 = 0; + const V* ptr1 = & input.at_(row, 0); + const V* ptr2 = & input.at_(row + 1, 0); + const V* ptr3 = & input.at_(row + 2, 0); + for (unsigned col = 0; col < ncols; col += scale) + { + S sum; + sum = *ptr1 + *(ptr1 + 1) + *(ptr1 + 2); + sum += *ptr2 + *(ptr2 + 1) + *(ptr2 + 2); + sum += *ptr3 + *(ptr3 + 1) + *(ptr3 + 2); + ptr1 += 3; + ptr2 += 3; + ptr3 += 3; + + S val = sum / area; + *p_sub++ = val; + + h_sum += val; + h_sum_2 += val * val; + + // exception + p_integ->first() = h_sum; + p_integ->second() = h_sum_2; + + p_integ += 1; + } + } + + unsigned b_next = 2 * border_thickness; + + p_sub += b_next; + p_integ += b_next; + + for (row += scale; row < nrows; row += scale) + { + S h_sum = 0, h_sum_2 = 0; + const V* ptr1 = & input.at_(row, 0); + const V* ptr2 = & input.at_(row + 1, 0); + const V* ptr3 = & input.at_(row + 2, 0); + for (unsigned col = 0; col < ncols; col += scale) + { + S sum; + sum = *ptr1 + *(ptr1 + 1) + *(ptr1 + 2); + sum += *ptr2 + *(ptr2 + 1) + *(ptr2 + 2); + sum += *ptr3 + *(ptr3 + 1) + *(ptr3 + 2); + ptr1 += 3; + ptr2 += 3; + ptr3 += 3; + + S val = sum / area; + *p_sub++ = val; + + h_sum += val; + h_sum_2 += val * val; + + p_integ->first() = h_sum + (p_integ + up)->first(); + p_integ->second() = h_sum_2 + (p_integ + up)->second(); + + p_integ += 1; + } + + p_sub += b_next; + p_integ += b_next; + } + + trace::exiting("subsampling::impl::integral_3"); + return sub; + } + + + template <typename I, typename J> + inline + mln_concrete(I) + integral_2(const Image<I>& input_, + Image<J>& integral_sum_sum_2_, + const mln_domain(I)& output_domain, + unsigned border_thickness) + { + trace::entering("subsampling::impl::integral_2"); + + const unsigned scale = 2; + + const I& input = exact(input_); + J& integral_sum_sum_2 = exact(integral_sum_sum_2_); + const unsigned area = scale * scale; + + typedef mln_value(I) V; + typedef mln_sum(V) S; + typedef mln_site(I) P; + typedef mln_value(J) V2; + + mlc_bool(P::dim == 2)::check(); + mln_precondition(input.is_valid()); + mln_precondition(input.domain().pmin() == literal::origin); + mln_precondition(scale > 1); + + mln_concrete(I) sub(output_domain, border_thickness); + V* p_sub = sub.buffer(); + + integral_sum_sum_2.init_(output_domain, border_thickness); + V2* p_integ = integral_sum_sum_2.buffer(); + + const unsigned up = sub.delta_index(dpoint2d(-1, 0)); + + const unsigned nrows = 2 * output_domain.nrows(); + const unsigned ncols = 2 * output_domain.ncols(); + + extension::fill(sub, 0); + + unsigned b_offset = sub.delta_index(dpoint2d(border_thickness, + border_thickness)); + p_sub += b_offset; + p_integ += b_offset; + + unsigned row = 0; + { + S h_sum = 0, h_sum_2 = 0; + const V* ptr1 = & input.at_(row, 0); + const V* ptr2 = & input.at_(row + 1, 0); + for (unsigned col = 0; col < ncols; col += scale) + { + S sum; + sum = *ptr1 + *(ptr1 + 1); + sum += *ptr2 + *(ptr2 + 1); + ptr1 += 2; + ptr2 += 2; + + S val = sum / area; + *p_sub++ = val; + + h_sum += val; + h_sum_2 += val * val; + + // exception + p_integ->first() = h_sum; + p_integ->second() = h_sum_2; + + p_integ += 1; + } + } + + unsigned b_next = 2 * border_thickness; + + p_sub += b_next; + p_integ += b_next; + + for (row += scale; row < nrows; row += scale) + { + S h_sum = 0, h_sum_2 = 0; + const V* ptr1 = & input.at_(row, 0); + const V* ptr2 = & input.at_(row + 1, 0); + for (unsigned col = 0; col < ncols; col += scale) + { + S sum; + sum = *ptr1 + *(ptr1 + 1); + sum += *ptr2 + *(ptr2 + 1); + ptr1 += 2; + ptr2 += 2; + + S val = sum / area; + *p_sub++ = val; + + h_sum += val; + h_sum_2 += val * val; + + p_integ->first() = h_sum + (p_integ + up)->first(); + p_integ->second() = h_sum_2 + (p_integ + up)->second(); + + p_integ += 1; + } + + p_sub += b_next; + p_integ += b_next; + } + + trace::exiting("subsampling::impl::integral_2"); + return sub; + } + + + template <typename I, typename J> + inline + mln_concrete(I) + integral(const Image<I>& input, unsigned scale, + Image<J>& integral_sum_sum_2, + const mln_domain(I)& output_domain, unsigned border_thickness) + { + // mln_precondition(input.nrows() % scale == 0); + // mln_precondition(input.ncols() % scale == 0); + if (scale == 3) + return integral_3(input, integral_sum_sum_2, + output_domain, border_thickness); + else if (scale == 2) + return integral_2(input, integral_sum_sum_2, + output_domain, border_thickness); + else + std::cerr << "NYI!" << std::endl; + + typedef mln_concrete(I) output_t; + return output_t(); + } + + + } // end of namespace mln::subsampling::impl + + + + + // Facades. + + template <typename I, typename J> + inline + mln_concrete(I) + integral(const Image<I>& input_, unsigned scale, + Image<J>& integral_sum_sum_2, + const mln_domain(I)& output_domain, unsigned border_thickness) + { + trace::entering("subsampling::integral"); + + const I& input = exact(input_); + + mln_precondition(input.is_valid()); + mln_precondition(input.domain().pmin() == literal::origin); + mln_precondition(scale > 1); + + mln_concrete(I) + output = impl::integral(input, scale, integral_sum_sum_2, + output_domain, border_thickness); + + trace::exiting("subsampling::integral"); + return output; + } + + template <typename I, typename J> + inline + mln_concrete(I) + integral(const Image<I>& input_, unsigned scale, + Image<J>& integral_sum_sum_2) + { + trace::entering("subsampling::integral"); + + const I& input = exact(input_); + + mln_precondition(input.is_valid()); + mln_precondition(input.domain().pmin() == literal::origin); + mln_precondition(scale > 1); + + box<mln_site(I)> b = make::box2d((input.nrows() + scale - 1) / scale, + (input.ncols() + scale - 1) / scale); + mln_concrete(I) output; + output = integral(input_, scale, integral_sum_sum_2, + b, mln::border::thickness); + + trace::exiting("subsampling::integral"); + return output; + } + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace scribo::subsampling + +} // end of namespace scribo + + +#endif // ! SCRIBO_SUBSAMPLING_INTEGRAL_HH -- 1.5.6.5