Commit 542fcf95 authored by Guillaume Lazzara's avatar Guillaume Lazzara
Browse files

Improve stats computation in Sauvola Multi-scale.

	* binarization/internal/first_pass_functor.hh: Use
	sauvola_threshold routine.

	* binarization/sauvola_threshold.hh: Remove debug and fix invalid
	read in integral image.

	* canvas/integral_browsing.hh,
	* subsampling/integral_single_image.hh: Fix stats computation.

	* src/binarization/sauvola_ms.cc: Fix window parameter and make
	debug output optional.
parent 32e461a6
2009-12-11 Guillaume Lazzara <z@lrde.epita.fr>
Improve Sauvola Multi-scale.
* binarization/internal/first_pass_functor.hh: Use
sauvola_threshold routine.
* binarization/sauvola_threshold.hh: Remove debug and fix invalid
read in integral image.
* canvas/integral_browsing.hh,
* subsampling/integral_single_image.hh: Fix stats computation.
* src/binarization/sauvola_ms.cc: Fix window parameter and make
debug output optional.
2009-12-04 Guillaume Lazzara <z@lrde.epita.fr>
Optimize Sauvola's multiscale binarization.
......
......@@ -98,8 +98,16 @@ namespace scribo
// 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);
// value::int_u8 t_p = mean * (one_k + k_R * stddev);
// std::cout << t_p << ", ";
// std::cout << input.element(p) << " - " << t_p << std::endl;
value::int_u8 t_p = sauvola_threshold_formula(mean, stddev);
// std::cout << input.point_at_index(p)
// << " - " << sauvola_threshold_formula(mean, stddev);
msk.element(p) = input.element(p) < t_p;
t_sub.element(p) = t_p;
......@@ -128,6 +136,8 @@ namespace scribo
void finalize()
{
mln_assertion(! pxl.is_valid());
// std::cout << std::endl << " ------- " << std::endl;
}
};
......
......@@ -46,6 +46,9 @@
# include <scribo/core/init_integral_image.hh>
#include <mln/io/pgm/save.hh>
namespace scribo
{
......@@ -168,8 +171,8 @@ namespace scribo
// Window half width.
int w_2 = win_width >> 1;
int row_min = std::max(0, p.row() - w_2);
int col_min = std::max(0, p.col() - w_2);
int row_min = std::max(0, p.row() - w_2 - 1);
int col_min = std::max(0, p.col() - w_2 - 1);
int row_max = std::min(static_cast<int>(simple.nrows()) - 1,
p.row() + w_2);
......@@ -177,12 +180,7 @@ namespace scribo
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);
double wh = (row_max - row_min) * (col_max - col_min);
// Mean.
double m_x_y_tmp = (simple.at_(row_max, col_max)
......@@ -200,9 +198,38 @@ namespace scribo
double s_x_y = std::sqrt((s_x_y_tmp - (m_x_y_tmp * m_x_y_tmp) / wh) / (wh - 1.f));
// if (p == point2d(3,3))// || p == point2d(4,4) || p == point2d(1,1))
// {
// // std::cout << "p" << p << " - A(" << row_min << ", " << col_min
// // << ") = " << simple.at_(row_min, col_min)
// << " - B(" << row_min << ", " << col_max
// << ") = " << simple.at_(row_min, col_max)
// << " - C(" << row_max << ", " << col_min
// << ") = " << simple.at_(row_max, col_min)
// << " - D(" << row_max << ", " << col_max
// << ") = " << simple.at_(row_max, col_max)
// << " - n = " << wh
// << std::endl;
// << std::endl;
// }
// Thresholding.
double t_x_y = sauvola_threshold_formula(m_x_y, s_x_y, k, R);
// std::cout << p
// << " - m = " << m_x_y
// << " - s = " << s_x_y
// << " - t = " << t_x_y
// << " - sum = " << m_x_y_tmp
// << " - sum_2 = " << s_x_y_tmp
// << std::endl;
return t_x_y;
}
......@@ -219,19 +246,12 @@ 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);
// 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.
......@@ -251,7 +271,7 @@ namespace scribo
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 + 0.14 * ((s_x_y / 128) - 1.0));
double t_x_y = m_x_y * (1.0 + 0.34 * ((s_x_y / 128) - 1.0));
return t_x_y;
}
......@@ -438,6 +458,25 @@ namespace scribo
mln_precondition(mln_site_(I)::dim == 2);
mln_precondition(exact(input).is_valid());
// {
// J& simple_ = exact(simple);
// J& squared_ = exact(squared);
// mln_piter(J) p(simple_.domain());
// for_all(p)
// {
// std::cout << simple_(p) << ", ";
// }
// std::cout << std::endl << " ------- " << std::endl;
// for_all(p)
// {
// std::cout << squared_(p) << ", ";
// }
// std::cout << std::endl << " ------- " << std::endl;
// }
typedef mln_value(I) value_t;
mln_ch_value(I, value::int_u8)
output = internal::sauvola_threshold_dispatch(value_t(), exact(input),
......@@ -445,6 +484,8 @@ namespace scribo
exact(simple),
exact(squared));
// std::cout << std::endl << " ------- " << std::endl;
io::pgm::save(output, "ref_2_t.pgm");
trace::exiting("scribo::text::ppm2pbm");
return output;
}
......
......@@ -58,9 +58,13 @@ namespace scribo
double& mean, double& stddev)
{
mean = sum / n;
stddev = std::sqrt(sum_2 / n - mean * mean);
// stddev = std::sqrt(sum_2 / n - mean * mean);
// std::cout << "(" << mean << " - " << stddev << " - " << n << "),";
// unbias version:
// stddev = std::sqrt((sum_2 - n * mean * mean) / (n - 1));
stddev = std::sqrt((sum_2 - sum * sum / n) / (n - 1));
}
} // end of namespace scribo::canvas::internal
......@@ -72,6 +76,7 @@ namespace scribo
void integral_browsing(const image2d<util::couple<double, double> >& ima,
unsigned step,
unsigned w, unsigned h,
unsigned s,
F& functor)
{
typedef util::couple<double, double> V;
......@@ -112,6 +117,8 @@ namespace scribo
double mean, stddev;
unsigned s_2 = s * s;
// -------------------------------
// T (top)
......@@ -146,7 +153,7 @@ namespace scribo
// D
internal::compute_stats(d_ima->first(),
d_ima->second(),
size_tl,
size_tl * s_2,
mean, stddev);
functor.exec(mean, stddev);
d_ima += step;
......@@ -166,7 +173,7 @@ namespace scribo
// D - C
internal::compute_stats(d_ima->first() - c_ima->first(),
d_ima->second() - c_ima->second(),
size_tc,
size_tc * s_2,
mean, stddev);
functor.exec(mean, stddev);
c_ima += step;
......@@ -188,7 +195,7 @@ namespace scribo
// D* - C
internal::compute_stats(d_sum - c_ima->first(),
d_sum_2 - c_ima->second(),
size_tr,
size_tr * s_2,
mean, stddev);
functor.exec(mean, stddev);
c_ima += step;
......@@ -239,7 +246,7 @@ namespace scribo
// D - B
internal::compute_stats(d_ima->first() - b_ima->first(),
d_ima->second() - b_ima->second(),
size_ml,
size_ml * s_2,
mean, stddev);
functor.exec(mean, stddev);
b_ima += step;
......@@ -258,11 +265,40 @@ namespace scribo
for (; col <= max_col_mid; col += step)
{
// D + A - B - C
// if (row == 3 && col == 3)
// std::cout << "p(" << row << "," << col << ") - "
// << "A" << ima.point_at_index(a_ima - ima.buffer())
// << "=" << a_ima->first() << " - "
// << "B" << ima.point_at_index(b_ima - ima.buffer())
// << "=" << b_ima->first() << " - "
// << "C" << ima.point_at_index(c_ima - ima.buffer())
// << "=" << c_ima->first() << " - "
// << "D" << ima.point_at_index(d_ima - ima.buffer())
// << "=" << d_ima->first() << " - "
// << "n =" << size_mc << " - "
// << "n*s_2 =" << size_mc * s_2
// << std::endl;
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,
size_mc * s_2,
mean, stddev);
functor.exec(mean, stddev);
// std::cout << " - " << mean
// << " - " << stddev
// << " - " << (d_ima->first() - b_ima->first()) + (a_ima->first() - c_ima->first())
// << " - " << (d_ima->second() - b_ima->second()) + (a_ima->second() - c_ima->second())
// << std::endl;
a_ima += step;
b_ima += step;
c_ima += step;
......@@ -283,7 +319,7 @@ namespace scribo
// 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,
size_mr * s_2,
mean, stddev);
functor.exec(mean, stddev);
a_ima += step;
......@@ -332,7 +368,7 @@ namespace scribo
// D* - B
internal::compute_stats(d_ima->first() - b_ima->first(),
d_ima->second() - b_ima->second(),
size_bl,
size_bl * s_2,
mean, stddev);
functor.exec(mean, stddev);
b_ima += step;
......@@ -354,7 +390,7 @@ namespace scribo
// 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,
size_bc * s_2,
mean, stddev);
// std::cout << (d_ima->second() - b_ima->second()) + (a_ima->second() - c_ima->second()) << std::endl;
......@@ -385,7 +421,7 @@ namespace scribo
// 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,
size_br * s_2,
mean, stddev);
functor.exec(mean, stddev);
a_ima += step;
......
......@@ -26,10 +26,10 @@
#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/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>
......@@ -96,6 +96,7 @@ namespace mln
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 s,
unsigned q, unsigned i, unsigned w,
const image2d<util::couple<double,double> >& integral_sum_sum_2)
// lambdas: limits of component cardinality at this scale
......@@ -111,22 +112,38 @@ namespace mln
tt.restart();
unsigned w_2 = w / 2 * ratio;
unsigned
w_local = w * ratio,
w_local_h = w_local,
w_local_w = w_local;
if (! (w_local % 2))
{
--w_local_w;
++w_local_h;
}
// std::cout << "Echelle " << i
// << " - w_local_h = " << w_local_h
// << " - w_local_w = " << w_local_w
// << " - w_1 = " << (w_local - 1) * s + 1
// << std::endl;
// std::cout << "Ratio = " << ratio << std::endl;
// std::cout << " -----------" << std::endl;
// 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,
w_local_w, w_local_h,
s,
f);
// debug::println("mask", f.msk);
// debug::println("parent", f.parent);
// debug::println("card", f.card);
t_ = tt;
if (mln::debug::quiet)
if (! mln::debug::quiet)
std::cout << "1st pass - " << t_ << std::endl;
tt.restart();
......@@ -175,7 +192,6 @@ namespace mln
{
for (unsigned l = 0; l < ptr.size(); ++l)
std::memset(ptr(l), i, ratio * sizeof(mln_value_(I)));
// debug(sq) = i;
}
}
......@@ -188,7 +204,6 @@ namespace mln
{
for (unsigned l = 0; l < ptr.size(); ++l)
std::memset(ptr(l), i, ratio * sizeof(mln_value_(I)));
// debug(sq) = i;
}
}
......@@ -204,12 +219,10 @@ namespace mln
}
t_ = tt;
if (mln::debug::quiet)
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
......@@ -245,9 +258,9 @@ namespace mln
// 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));
int more_offset = - ((4 * s) - in.ncols() % (4 * s));
if (more_offset == (4 * s))
if (more_offset == - (static_cast<int>(4*s)))
more_offset = 0; // No offset needed.
const int
......@@ -267,7 +280,7 @@ namespace mln
delta3 = t_ima[3].delta_index(dpoint2d(+1, -1)),
eor1 = in.delta_index(dpoint2d(+4 * s, - in.ncols() - in.border())) + more_offset,
eor1 = in.delta_index(dpoint2d(+4 * s, - in.ncols())) + 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()));
......@@ -700,6 +713,499 @@ namespace mln
// template <typename I, typename J, typename K>
// mln_ch_value(I, unsigned)
// binarize_generic_debug(const I& in, const J& e2, const util::array<K>& t_ima,
// unsigned s)
// {
// mln_ch_value(I,unsigned) out;
// initialize(out, in);
// data::fill(out, 0);
// typedef const mln_value(K)* ptr_type;
// 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);
// const mln_value(J)* ptr_e2 = & e2.at_(0, 0);
// const mln_value(I)* ptr__in = & in.at_(0, 0);
// unsigned* 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.
// std::cout << in.ncols() << std::endl;
// std::cout << in.ncols() % (4 * s) << std::endl;
// int more_offset = - ((4 * s) - in.ncols() % (4 * s));
// if (more_offset == - (4*s))
// more_offset = 0; // No offset needed.
// std::cout << "more_offset == " << more_offset << std::endl;
// std::cout << "- b1 = " << in.border()
// << "- b2 = " << t_ima[2].border()
// << "- b3 = " << t_ima[3].border()
// << "- b4 = " << t_ima[4].border()
// << std::endl;
// 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())) + 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()));
// unsigned pid = 0;
// 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 = pid++;
// ++ptr__out; ++ptr__in;
// }
// *ptr__out = pid++;
// ptr__out += delta1; ptr__in += delta1;
// }
// for (unsigned j = 1; j < s; ++j)
// {
// *ptr__out = pid++;
// ++ptr__out; ++ptr__in;
// }
// *ptr__out = pid++;
// 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 = pid++;
// ++ptr__out; ++ptr__in;
// }
// *ptr__out = pid++;
// ptr__out += delta1; ptr__in += delta1;
// }
// for (unsigned j = 1; j < s; ++j)
// {
// *ptr__out = pid++;
// ++ptr__out; ++ptr__in;
// }
// *ptr__out = pid++;
// 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 = pid++;
// ++ptr__out; ++ptr__in;
// }
// *ptr__out = pid++;
// ptr__out += delta1; ptr__in += delta1;
// }
// for (unsigned j = 1; j < s; ++j)
// {
// *ptr__out = pid++;
// ++ptr__out; ++ptr__in;
// }
// *ptr__out = pid++;
// 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 = pid++;
// ++ptr__out; ++ptr__in;
// }
// *ptr__out = pid++;
// ptr__out += delta1; ptr__in += delta1;
// }
// for (unsigned j = 1; j < s; ++j)
// {
// *ptr__out = pid++;
// ++ptr__out; ++ptr__in;
// }
// *ptr__out = pid++;
// 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 = pid++;
// ++ptr__out; ++ptr__in;