[Olena-patches] 1005: Better milena median.

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Thierry Geraud &lt;thierry.geraud(a)lrde.epita.fr&gt; Better milena median. * tests/median.cc: Update. * mln/convert/to_dpoint.hh: New. * mln/debug/println.hh (println): Add endl. * mln/level/median.hh: New. * mln/io/save_pgm.hh: New. * mln/io/load_pgm.hh: New. * mln/core/box.hh (len): New. * mln/core/concept/box.hh: Likewise. * mln/core/point.hh (zero): New. * mln/core/dpoint.hh: Likewise. * mln/core/concept/window.hh (dpoint, point): New. * mln/core/window.hh, * mln/core/rectangle2d.hh: Update. * mln/core/concept/genpoint.hh: Remove dead line. * mln/core/concept/dpoint.hh (operator+): New. * mln/core/image2d_b.hh: . * mln/core/internal/force_exact.hh: Avoid obj creation. * mln/accu: New. * mln/value/histo.hh: Rename as... * mln/value/median.hh: ...this. * mln/accu/histo.hh: Rename as... * mln/accu/median.hh: ...this. * mln/accu/median_alt.hh: New. mln/accu/histo.hh | 20 +- mln/accu/median.hh | 147 +++++++++----------- mln/accu/median_alt.hh | 279 +++++++++++++++++++++++++++++++++++++++ mln/convert/to_dpoint.hh | 66 +++++++++ mln/core/box.hh | 13 + mln/core/concept/box.hh | 6 mln/core/concept/dpoint.hh | 22 +++ mln/core/concept/genpoint.hh | 1 mln/core/concept/window.hh | 6 mln/core/dpoint.hh | 13 + mln/core/image2d_b.hh | 32 +--- mln/core/internal/force_exact.hh | 6 mln/core/point.hh | 11 + mln/core/rectangle2d.hh | 6 mln/core/window.hh | 61 ++++++++ mln/debug/println.hh | 1 mln/io/load_pgm.hh | 177 ++++++++++++++++++++++++ mln/io/save_pgm.hh | 74 ++++++++++ mln/level/median.hh | 136 +++++++++++++++++++ tests/median.cc | 230 ++++++++++++++++++-------------- 20 files changed, 1092 insertions(+), 215 deletions(-) Index: tests/median.cc --- tests/median.cc (revision 1004) +++ tests/median.cc (working copy) @@ -30,155 +30,185 @@ * \brief Tests on mln::value::median<S>. */ +#include <cmath> + +#include <mln/core/image2d_b.hh> #include <mln/value/int_u.hh> -#include <mln/value/median.hh> +#include <mln/level/fill.hh> +#include <mln/level/median.hh> +#include <mln/debug/println.hh> + +#include <mln/core/rectangle2d.hh> + +#include <mln/io/load_pgm.hh> +#include <mln/io/save_pgm.hh> -int main() -{ using namespace mln; using namespace mln::value; - typedef set_<int_u8> S; - median<S> m(S::the()); +int_u8 f(const point2d& p) +{ + return unsigned((2 + + std::cos(float(p.row())) + + std::sin(float(p.col()))) * 63.9); +} +int main() { - unsigned vals[] = { 42, 69, 51, 12, 51, 12, 42 }; - unsigned n = sizeof(vals)/sizeof(unsigned); - for (unsigned i = 0; i < n; ++i) - { - std::cout << "taking " << vals[i] << ':' << std::endl; - m.take(vals[i]); - std::cout << m << std::endl; - } +// typedef set_<int_u8> S; - for (int i = int(n) - 1; i >= 0; --i) - { - std::cout << "untaking " << vals[i] << ':' << std::endl; - m.untake(vals[i]); - std::cout << m << std::endl; - } +// median<S> m(S::the()); - } +// { +// unsigned vals[] = { 42, 69, 51, 12, 51, 12, 42 }; +// unsigned n = sizeof(vals)/sizeof(unsigned); - { +// for (unsigned i = 0; i < n; ++i) +// { +// std::cout << "taking " << vals[i] << ':' << std::endl; +// m.take(vals[i]); +// std::cout << m << std::endl; +// } - unsigned vals[] = { 42, 69, 51, 12, 51, 12, 42 }; - unsigned n = sizeof(vals)/sizeof(unsigned); +// for (int i = int(n) - 1; i >= 0; --i) +// { +// std::cout << "untaking " << vals[i] << ':' << std::endl; +// m.untake(vals[i]); +// std::cout << m << std::endl; +// } - for (unsigned i = 0; i < n; ++i) - { - std::cout << "taking " << vals[i] << ':' << std::endl; - m.take(vals[i]); - std::cout << m << std::endl; - } +// } - for (unsigned i = 0; i < n; ++i) - { - std::cout << "untaking " << vals[i] << ':' << std::endl; - m.untake(vals[i]); - std::cout << m << std::endl; - } + image2d_b<int_u8> + lena = io::load_pgm("lena.pgm"), + out(lena.domain()); + + rectangle2d rec(64, 64); + level::median(lena, rec, out); + io::save_pgm(out, "out.pgm"); } - { - unsigned vals[] = { 42, 42, 69, 69, 51, 51, 12, 12, 51, 51, 12, 12, 42, 42 }; - unsigned n = sizeof(vals)/sizeof(unsigned); +// { - for (unsigned i = 0; i < n; ++i) - { - std::cout << "taking " << vals[i] << ':' << std::endl; - m.take(vals[i]); - std::cout << m << std::endl; - } +// unsigned vals[] = { 42, 69, 51, 12, 51, 12, 42 }; +// unsigned n = sizeof(vals)/sizeof(unsigned); - for (int i = int(n) - 1; i >= 0; --i) - { - std::cout << "untaking " << vals[i] << ':' << std::endl; - m.untake(vals[i]); - std::cout << m << std::endl; - } +// for (unsigned i = 0; i < n; ++i) +// { +// std::cout << "taking " << vals[i] << ':' << std::endl; +// m.take(vals[i]); +// std::cout << m << std::endl; +// } - } +// for (unsigned i = 0; i < n; ++i) +// { +// std::cout << "untaking " << vals[i] << ':' << std::endl; +// m.untake(vals[i]); +// std::cout << m << std::endl; +// } +// } - { - unsigned vals[] = { 42, 42, 69, 69, 51, 51, 12, 12, 51, 51, 12, 12, 42, 42 }; - unsigned n = sizeof(vals)/sizeof(unsigned); - for (unsigned i = 0; i < n; ++i) - { - std::cout << "taking " << vals[i] << ':' << std::endl; - m.take(vals[i]); - std::cout << m << std::endl; - } +// { - for (unsigned i = 0; i < n; ++i) - { - std::cout << "untaking " << vals[i] << ':' << std::endl; - m.untake(vals[i]); - std::cout << m << std::endl; - } +// unsigned vals[] = { 42, 42, 69, 69, 51, 51, 12, 12, 51, 51, 12, 12, 42, 42 }; +// unsigned n = sizeof(vals)/sizeof(unsigned); - } +// for (unsigned i = 0; i < n; ++i) +// { +// std::cout << "taking " << vals[i] << ':' << std::endl; +// m.take(vals[i]); +// std::cout << m << std::endl; +// } +// for (int i = int(n) - 1; i >= 0; --i) +// { +// std::cout << "untaking " << vals[i] << ':' << std::endl; +// m.untake(vals[i]); +// std::cout << m << std::endl; +// } +// } - { +// { - unsigned vals[] = { 42, 69, 51, 12, 51, 12, 42, 69, 51, 12, 51, 12, 42 }; - unsigned n = sizeof(vals)/sizeof(unsigned); +// unsigned vals[] = { 42, 42, 69, 69, 51, 51, 12, 12, 51, 51, 12, 12, 42, 42 }; +// unsigned n = sizeof(vals)/sizeof(unsigned); - for (unsigned i = 0; i < n; ++i) - { - std::cout << "taking " << vals[i] << ':' << std::endl; - m.take(vals[i]); - std::cout << m << std::endl; - } +// for (unsigned i = 0; i < n; ++i) +// { +// std::cout << "taking " << vals[i] << ':' << std::endl; +// m.take(vals[i]); +// std::cout << m << std::endl; +// } - for (int i = int(n) - 1; i >= 0; --i) - { - std::cout << "untaking " << vals[i] << ':' << std::endl; - m.untake(vals[i]); - std::cout << m << std::endl; - } +// for (unsigned i = 0; i < n; ++i) +// { +// std::cout << "untaking " << vals[i] << ':' << std::endl; +// m.untake(vals[i]); +// std::cout << m << std::endl; +// } - } +// } - { - unsigned vals[] = { 42, 69, 51, 12, 51, 12, 42, 69, 51, 12, 51, 12, 42 }; - unsigned n = sizeof(vals)/sizeof(unsigned); - for (unsigned i = 0; i < n; ++i) - { - std::cout << "taking " << vals[i] << ':' << std::endl; - m.take(vals[i]); - std::cout << m << std::endl; - } +// { - for (unsigned i = 0; i < n; ++i) - { - std::cout << "untaking " << vals[i] << ':' << std::endl; - m.untake(vals[i]); - std::cout << m << std::endl; - } +// unsigned vals[] = { 42, 69, 51, 12, 51, 12, 42, 69, 51, 12, 51, 12, 42 }; +// unsigned n = sizeof(vals)/sizeof(unsigned); - } +// for (unsigned i = 0; i < n; ++i) +// { +// std::cout << "taking " << vals[i] << ':' << std::endl; +// m.take(vals[i]); +// std::cout << m << std::endl; +// } +// for (int i = int(n) - 1; i >= 0; --i) +// { +// std::cout << "untaking " << vals[i] << ':' << std::endl; +// m.untake(vals[i]); +// std::cout << m << std::endl; +// } + +// } + + +// { + +// unsigned vals[] = { 42, 69, 51, 12, 51, 12, 42, 69, 51, 12, 51, 12, 42 }; +// unsigned n = sizeof(vals)/sizeof(unsigned); + +// for (unsigned i = 0; i < n; ++i) +// { +// std::cout << "taking " << vals[i] << ':' << std::endl; +// m.take(vals[i]); +// std::cout << m << std::endl; +// } + +// for (unsigned i = 0; i < n; ++i) +// { +// std::cout << "untaking " << vals[i] << ':' << std::endl; +// m.untake(vals[i]); +// std::cout << m << std::endl; +// } + +// } -} Index: mln/convert/to_dpoint.hh --- mln/convert/to_dpoint.hh (revision 0) +++ mln/convert/to_dpoint.hh (revision 0) @@ -0,0 +1,66 @@ +// Copyright (C) 2007 EPITA Research and Development Laboratory +// +// This file is part of the Olena Library. This library is free +// software; you can redistribute it and/or modify it under the terms +// of the GNU General Public License version 2 as published by the +// Free Software Foundation. +// +// This library 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 this library; see the file COPYING. If not, write to +// the Free Software Foundation, 51 Franklin Street, Fifth Floor, +// Boston, MA 02111-1307, USA. +// +// As a special exception, you may use this file as part of a free +// software library 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 MLN_CONVERT_TO_DPOINT_HH +# define MLN_CONVERT_TO_DPOINT_HH + +/*! \file mln/convert/to_dpoint.hh + * + * \brief Convertions to mln::Dpoint. + */ + + +namespace mln +{ + + namespace convert + { + + template <typename P> + mln_dpoint(P) to_dpoint(const GenPoint<P>& p); + + +# ifndef MLN_INCLUDE_ONLY + + template <typename P> + mln_dpoint(P) to_dpoint(const GenPoint<P>& p_) + { + const P& p = p_.force_exact_(); + mln_dpoint(P) tmp; + for (unsigned i = 0; i < P::dim; ++i) + tmp[i] = p[i]; + return tmp; + } + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::convert + +} // end of namespace mln + + +#endif // ! MLN_CONVERT_TO_DPOINT_HH Index: mln/debug/println.hh --- mln/debug/println.hh (revision 1004) +++ mln/debug/println.hh (working copy) @@ -74,6 +74,7 @@ std::cout << input(mk_point2d(row, col)) << ' '; std::cout << std::endl; } + std::cout << std::endl; } } // end of namespace mln::debug::impl Index: mln/level/median.hh --- mln/level/median.hh (revision 0) +++ mln/level/median.hh (revision 0) @@ -0,0 +1,136 @@ +// Copyright (C) 2007 EPITA Research and Development Laboratory +// +// This file is part of the Olena Library. This library is free +// software; you can redistribute it and/or modify it under the terms +// of the GNU General Public License version 2 as published by the +// Free Software Foundation. +// +// This library 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 this library; see the file COPYING. If not, write to +// the Free Software Foundation, 51 Franklin Street, Fifth Floor, +// Boston, MA 02111-1307, USA. +// +// As a special exception, you may use this file as part of a free +// software library 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 MLN_LEVEL_MEDIAN_HH +# define MLN_LEVEL_MEDIAN_HH + +/*! \file mln/level/median.hh + * + * \brief Median filter on an image. + */ + +# include <mln/core/concept/image.hh> +# include <mln/core/window2d.hh> +# include <mln/accu/median.hh> + + +namespace mln +{ + + namespace level + { + + /*! FIXME: Median the image \p ima with the values of the image \p data. + * + * \param[in] input The image to be filtered. + * \param[in] win The window. + * \param[in,out] output The output image. + * + * \warning The definition domain of \p ima has to be included in + * the one of \p data. + * + * \pre \p ima has to be initialized. + * + * \todo Test domain inclusion. + */ + template <typename I, typename W, typename O> + void median(const Image<I>& input, const Window<W>& win, + Image<O>& output); + + +# ifndef MLN_INCLUDE_ONLY + + + namespace impl + { + + template <typename I, typename W, typename O> + void median(const I& input, + const W& win, + O& output) + { + dpoint2d dp = mk_dpoint2d(0, 1); + + window2d + w_plus = win - (win - dp), + w_minus = (win - dp) - win; + + accu::median_on<mln_value(I)> med; + + point2d p; + mln_qiter(W) q(win, p); + mln_qiter(W) q_plus(w_plus, p); + mln_qiter(W) q_minus(w_minus, p); + + for (p.row() = input.domain().pmin().row(); + p.row() <= input.domain().pmax().row(); + ++p.row()) + { + // init + med.init(); + p.col() = input.domain().pmin().col(); + for_all(q) + if (input.has(q)) + med.take(input(q)); + output(p) = med; + + // other col + for (++p.col(); + p.col() <= input.domain().pmax().col(); + ++p.col()) + { + for_all(q_plus) + if (input.has(q_plus)) + med.take(input(q_plus)); + for_all(q_minus) + if (input.has(q_minus)) + med.untake(input(q_minus)); + output(p) = med; + } + } + } + + } // end of namespace mln::level::impl + + + // facade + + template <typename I, typename W, typename O> + void median(const Image<I>& input, const Window<W>& win, + Image<O>& output) + { + impl::median(exact(input), exact(win), exact(output)); + } + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::level + +} // end of namespace mln + + +#endif // ! MLN_LEVEL_MEDIAN_HH Index: mln/io/save_pgm.hh --- mln/io/save_pgm.hh (revision 0) +++ mln/io/save_pgm.hh (revision 0) @@ -0,0 +1,74 @@ +// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007 EPITA +// Research and Development Laboratory +// +// This file is part of the Olena Library. This library is free +// software; you can redistribute it and/or modify it under the terms +// of the GNU General Public License version 2 as published by the +// Free Software Foundation. +// +// This library 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 this library; see the file COPYING. If not, write to +// the Free Software Foundation, 51 Franklin Street, Fifth Floor, +// Boston, MA 02111-1307, USA. +// +// As a special exception, you may use this file as part of a free +// software library 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 MLN_IO_SAVE_PGM_HH +# define MLN_IO_SAVE_PGM_HH + +# include <iostream> +# include <fstream> + +# include <mln/core/image2d_b.hh> +# include <mln/value/int_u.hh> + + +namespace mln +{ + + namespace io + { + + void save_pgm(const image2d_b<value::int_u8>& ima, const std::string& filename) + { + std::ofstream file(filename.c_str()); + if (not file) + { + std::cerr << "error: cannot open file '" << filename + << "'!"; + abort(); + } + file << "P5" << std::endl; + file << "# olena" << std::endl; + file << ima.ncols() << ' ' << ima.nrows() << std::endl; + file << "255" << std::endl; + point2d p = mk_point2d(ima.domain().pmin().row(), + ima.domain().pmin().col()); + size_t len = ima.ncols() * sizeof(unsigned char); + for (; + p.row() <= ima.domain().pmax().row(); + ++p.row()) + { + file.write((char*)(& ima(p)), len); + } + } + + } // end of namespace mln::io + +} // end of namespace mln + + +#endif // ! MLN_IO_SAVE_PGM_HH Index: mln/io/load_pgm.hh --- mln/io/load_pgm.hh (revision 0) +++ mln/io/load_pgm.hh (revision 0) @@ -0,0 +1,177 @@ +// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007 EPITA +// Research and Development Laboratory +// +// This file is part of the Olena Library. This library is free +// software; you can redistribute it and/or modify it under the terms +// of the GNU General Public License version 2 as published by the +// Free Software Foundation. +// +// This library 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 this library; see the file COPYING. If not, write to +// the Free Software Foundation, 51 Franklin Street, Fifth Floor, +// Boston, MA 02111-1307, USA. +// +// As a special exception, you may use this file as part of a free +// software library 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 MLN_IO_LOAD_PGM_HH +# define MLN_IO_LOAD_PGM_HH + +# include <iostream> +# include <fstream> +# include <string> + +# include <mln/core/image2d_b.hh> +# include <mln/value/int_u.hh> + + +namespace mln +{ + + namespace io + { + + namespace internal + { + + void abort() + { + std::cerr << " aborting." << std::endl; + exit(0); + } + + bool read_pnm_header(std::istream& istr, + char& type, + int& nrows, int& ncols, + bool test = false) + { + // check magic + if (istr.get() != 'P' ) + goto err; + type = istr.get(); + if (type < '1' or type > '6') + goto err; + if (istr.get() != '\n') + goto err; + + // skip comments + while (istr.peek() = '#') + { + std::string line; + std::getline(istr, line); + } + + // get size + istr >> ncols >> nrows; + if (nrows <= 0 or ncols <= 0) + goto err; + + // skip maxvalue + if (istr.get() != '\n') + goto err; + if (type != '1' and type != '4') + { + std::string line; + std::getline(istr, line); + } + return true; + + err: + if (not test) + { + std::cerr << "error: badly formed header!"; + abort(); + } + return false; + } + + void read_pnm_header(char ascii, char raw, + std::istream& istr, + char& type, + int& nrows, int& ncols) + { + read_pnm_header(istr, type, nrows, ncols); + if (not (type = ascii or type = raw)) + { + std::cerr << "error: bad pnm type; " + << "expected P" << ascii + << " or P" << raw + << ", get P" << type << "!"; + abort(); + } + } + + + /// load_ascii. + template <typename I> + void load_pgm_ascii(std::ifstream& file, I& ima) + { + mln_fwd_piter(I) p(ima.domain()); + for_all(p) + { + unsigned value; + file >> value; + ima(p) = value; + // FIXME: Test alt code below. + // file >> ima(p); + } + } + + + /// load_raw_2d. + template <typename I> + void load_pgm_raw_2d(std::ifstream& file, I& ima) + { + point2d p = mk_point2d(0, ima.domain().pmin().col()); + size_t len = ima.ncols() * sizeof(mln_value(I)); + for (p.row() = ima.domain().pmin().row(); + p.row() <= ima.domain().pmax().row(); + ++p.row()) + { + file.read((char*)(& ima(p)), len); + } + } + + + } // end of namespace mln::io::internal + + + image2d_b<value::int_u8> load_pgm(const std::string& filename) + { + std::ifstream file(filename.c_str()); + if (not file) + { + std::cerr << "error: file '" << filename + << "' not found!"; + abort(); + } + char type; + int nrows, ncols; + internal::read_pnm_header('2', '5', file, type, nrows, ncols); + image2d_b<value::int_u8> ima(nrows, ncols); + if (type = '5') + internal::load_pgm_raw_2d(file, ima); + else + if (type = '2') + internal::load_pgm_ascii(file, ima); + return ima; + } + + } // end of namespace mln::io + +} // end of namespace mln + + +#endif // ! MLN_IO_LOAD_PGM_HH Index: mln/core/dpoints_piter.hh Index: mln/core/box.hh --- mln/core/box.hh (revision 1004) +++ mln/core/box.hh (working copy) @@ -88,6 +88,11 @@ */ P& pmax(); + /*! \brief Give the length of the \p i-th side. + * \pre i < dim + */ + unsigned len(unsigned i) const; + /*! \brief Constructor without argument. */ box_(); @@ -151,6 +156,14 @@ } template <typename P> + unsigned + box_<P>::len(unsigned i) const + { + mln_precondition(i < P::dim); + return 1 + pmax_[i] - pmin_[i]; + } + + template <typename P> box_<P>::box_() { } Index: mln/core/point.hh --- mln/core/point.hh (revision 1004) +++ mln/core/point.hh (working copy) @@ -35,6 +35,7 @@ # include <mln/core/concept/point.hh> # include <mln/core/internal/coord_impl.hh> +# include <mln/fun/all.hh> namespace mln @@ -92,6 +93,9 @@ */ void set_all(C c); + /// Give the origin (all coordinates are 0). + static const point_<n,C>& zero(); + protected: C coord_[n]; }; @@ -133,6 +137,13 @@ coord_[i] = c; } + template <unsigned n, typename C> + const point_<n,C>& point_<n,C>::zero() + { + static const point_<n,C> zero_(all(0)); + return zero_; + } + # endif // ! MLN_INCLUDE_ONLY } // end of namespace mln Index: mln/core/rectangle2d.hh --- mln/core/rectangle2d.hh (revision 1004) +++ mln/core/rectangle2d.hh (working copy) @@ -56,6 +56,12 @@ struct rectangle2d : public Window< rectangle2d >, public internal::set_of_<dpoint2d> { + /// Point associated type. + typedef point2d point; + + /// Dpoint associated type. + typedef dpoint2d dpoint; + /*! \brief Piter type to browse a rectangle such as: "for each row * (increasing), for each column (increasing)." */ Index: mln/core/concept/genpoint.hh --- mln/core/concept/genpoint.hh (revision 1004) +++ mln/core/concept/genpoint.hh (working copy) @@ -60,7 +60,6 @@ */ template <typename E> struct GenPoint - // : virtual internal::object_<E> { /* Index: mln/core/concept/box.hh --- mln/core/concept/box.hh (revision 1004) +++ mln/core/concept/box.hh (working copy) @@ -52,6 +52,7 @@ /* const point& pmin() const; const point& pmax() const; + unsigned len(unsigned i) const; // FIXME: Doc! */ /*! \brief Return the bounding box of this point set. @@ -102,10 +103,7 @@ std::size_t count = 1; typedef typename E::point P; // helps g++-3.3.5 for (unsigned i = 0; i < P::dim; ++i) - count *- exact(this)->pmax()[i] - + 1 - - exact(this)->pmin()[i]; + count *= exact(this)->len(i); return count; } Index: mln/core/concept/dpoint.hh --- mln/core/concept/dpoint.hh (revision 1004) +++ mln/core/concept/dpoint.hh (working copy) @@ -78,6 +78,19 @@ D operator-(const Dpoint<D>& rhs); + /*! \brief Add the couple of delta-points \p lhs and \p rhs. + * + * \param[in] lhs A delta-point. + * \param[in] rhs Another delta-point. + * + * \return A delta-point (temporary object). + * + * \relates mln::Dpoint + */ + template <typename D> + D operator+(const Dpoint<D>& lhs, const Dpoint<D>& rhs); + + /*! \brief Equality comparison between a couple of delta-point \p lhs * and \p rhs. * @@ -149,6 +162,15 @@ return tmp; } + template <typename D> + D operator+(const Dpoint<D>& lhs, const Dpoint<D>& rhs) + { + D tmp; + for (unsigned i = 0; i < D::dim; ++i) + tmp[i] = exact(lhs)[i] + exact(rhs)[i]; + return tmp; + } + template <typename Dl, typename Dr> bool operator=(const Dpoint<Dl>& lhs, const Dpoint<Dr>& rhs) { Index: mln/core/concept/window.hh --- mln/core/concept/window.hh (revision 1004) +++ mln/core/concept/window.hh (working copy) @@ -47,6 +47,9 @@ struct Window : public Object<E> { /* + typedef point; + typedef dpoint; + typedef qiter; typedef fwd_qiter; typedef bkd_qiter; @@ -66,6 +69,9 @@ template <typename E> Window<E>::Window() { + typedef mln_point(E) point; + typedef mln_dpoint(E) dpoint; + typedef mln_qiter(E) qiter; typedef mln_fwd_qiter(E) fwd_qiter; typedef mln_bkd_qiter(E) bkd_qiter; Index: mln/core/dpoint.hh --- mln/core/dpoint.hh (revision 1004) +++ mln/core/dpoint.hh (working copy) @@ -35,6 +35,7 @@ # include <mln/core/concept/dpoint.hh> # include <mln/core/internal/coord_impl.hh> +# include <mln/fun/all.hh> namespace mln @@ -86,12 +87,15 @@ /*! \brief Constructor; coordinates are set by function \p f. */ template <typename F> - dpoint_(F f); + dpoint_(F f); // FIXME: Bound parameter! /*! \brief Set all coordinates to the value \p c. */ void set_all(C c); + /// Give the null delta-point (all coordinates are 0). + static const dpoint_<n,C>& zero(); + protected: C coord_[n]; }; @@ -133,6 +137,13 @@ coord_[i] = c; } + template <unsigned n, typename C> + const dpoint_<n,C>& dpoint_<n,C>::zero() + { + static const dpoint_<n,C> zero_(all(0)); + return zero_; + } + # endif // ! MLN_INCLUDE_ONLY } // end of namespace mln Index: mln/core/window.hh --- mln/core/window.hh (revision 1004) +++ mln/core/window.hh (working copy) @@ -34,8 +34,11 @@ */ # include <mln/core/concept/window.hh> +# include <mln/core/concept/genpoint.hh> # include <mln/core/internal/set_of.hh> # include <mln/core/dpoint.hh> + +# include <mln/convert/to_dpoint.hh> # include <mln/fun/all.hh> @@ -56,6 +59,12 @@ struct window_ : public Window< window_<D> >, public internal::set_of_<D> { + /// Point associated type. + typedef mln_point(D) point; + + /// Dpoint associated type. + typedef D dpoint; + /*! \brief Piter type to browse the points of a generic window * w.r.t. the ordering of delta-points. */ @@ -90,6 +99,25 @@ }; + // FIXME: Move both ops below to mln/core/concept/window.hh + + /// Shift a window \p win with a delta-point \p dp. + template <typename W> + window_<mln_dpoint(W)> operator+(const Window<W>& win, + const mln_dpoint(W)& dp); + + /// Shift a window \p win with the delta-point (-\p dp). + template <typename W> + window_<mln_dpoint(W)> operator-(const Window<W>& win, + const mln_dpoint(W)& dp); + + /// Substract \p rhs from \p lhs. + // FIXME: Give details! + template <typename Wl, typename Wr> + window_<mln_dpoint(Wl)> operator-(const Window<Wl>& lhs, + const Window<Wr>& rhs); + + # ifndef MLN_INCLUDE_ONLY @@ -112,6 +140,39 @@ return false; } + template <typename W> + window_<mln_dpoint(W)> operator+(const Window<W>& win, + const mln_dpoint(W)& dp) + { + window_<mln_dpoint(W)> tmp; + mln_qiter(W) q(win, W::point::zero()); + for_all(q) + tmp.insert(convert::to_dpoint(q) + dp); + return tmp; + } + + template <typename W> + window_<mln_dpoint(W)> operator-(const Window<W>& win, + const mln_dpoint(W)& dp) + { + return win + (-dp); + } + + template <typename W, typename Wr> + window_<mln_dpoint(W)> operator-(const Window<W>& lhs, + const Window<Wr>& rhs) + { + window_<mln_dpoint(W)> tmp; + mln_qiter(W) q(lhs, W::point::zero()); + for_all(q) + { + mln_dpoint(W) dp = convert::to_dpoint(q); + if (not exact(rhs).has(dp)) + tmp.insert(dp); + } + return tmp; + } + # endif // ! MLN_INCLUDE_ONLY } // end of namespace mln Index: mln/core/image2d_b.hh --- mln/core/image2d_b.hh (revision 1004) +++ mln/core/image2d_b.hh (working copy) @@ -180,7 +180,7 @@ allocate_(); std::memcpy(this->buffer_, rhs.buffer_, - ncells() * sizeof(value)); + ncells() * sizeof(T)); } // assignment @@ -199,7 +199,7 @@ allocate_(); std::memcpy(this->buffer_, rhs.buffer_, - ncells() * sizeof(value)); + ncells() * sizeof(T)); return *this; } @@ -233,7 +233,7 @@ image2d_b<T>::nrows() const { mln_precondition(this->has_data()); - return 1 + b_.pmax().row() - b_.pmin().row(); + return b_.len(0); } template <typename T> @@ -241,7 +241,7 @@ image2d_b<T>::ncols() const { mln_precondition(this->has_data()); - return 1 + b_.pmax().col() - b_.pmin().col(); + return b_.len(1); } template <typename T> @@ -249,10 +249,7 @@ image2d_b<T>::ncells() const { mln_precondition(this->has_data()); - std::size_t s = 1; - s *= nrows() + 2 * bdr_; - s *= ncols() + 2 * bdr_; - return s; + return vb_.npoints(); } template <typename T> @@ -260,10 +257,7 @@ image2d_b<T>::owns_(const point2d& p) const { mln_precondition(this->has_data()); - return p.row() >= vb_.pmin().row() - && p.row() <= vb_.pmax().row() - && p.col() >= vb_.pmin().col() - && p.col() <= vb_.pmax().col(); + return vb_.has(p); } template <typename T> @@ -304,17 +298,19 @@ { update_vb_(); unsigned - nr = nrows() + 2 * bdr_, - nc = ncols() + 2 * bdr_; - buffer_ = new T[ncells()]; + nr = vb_.len(0), + nc = vb_.len(1); + buffer_ = new T[nr * nc]; array_ = new T*[nr]; - T* buf = buffer_ - (b_.pmin().col() - bdr_); + T* buf = buffer_ - vb_.pmin().col(); for (unsigned i = 0; i < nr; ++i) { array_[i] = buf; buf += nc; } - array_ -= b_.pmin().row() - bdr_; + array_ -= vb_.pmin().row(); + mln_postcondition(vb_.len(0) = b_.nrows() + 2 * bdr_); + mln_postcondition(vb_.len(1) = b_.ncols() + 2 * bdr_); } template <typename T> @@ -328,7 +324,7 @@ } if (array_) { - array_ += b_.pmin().row() - bdr_; + array_ += vb_.pmin().row(); delete[] array_; array_ = 0; } Index: mln/core/internal/force_exact.hh --- mln/core/internal/force_exact.hh (revision 1004) +++ mln/core/internal/force_exact.hh (working copy) @@ -50,11 +50,11 @@ \ E& force_exact_() const \ { \ - static const E exact_obj; \ - static const Type& exact_obj_ref = exact_obj; \ + static const E* exact_obj; \ + static const Type& exact_obj_ref = *exact_obj; \ static const int exact_offset = \ (const char*)(void*)(&exact_obj_ref) \ - - (const char*)(void*)(&exact_obj); \ + - (const char*)(void*)( exact_obj); \ return *(E*)((char*)(this) - exact_offset); \ } Index: mln/accu/histo.hh --- mln/accu/histo.hh (revision 0) +++ mln/accu/histo.hh (working copy) @@ -25,8 +25,8 @@ // reasons why the executable file might be covered by the GNU General // Public License. -#ifndef MLN_VALUE_HISTO_HH -# define MLN_VALUE_HISTO_HH +#ifndef MLN_ACCU_HISTO_HH +# define MLN_ACCU_HISTO_HH /*! \file mln/value/histo.hh * @@ -43,7 +43,7 @@ namespace mln { - namespace value + namespace accu { @@ -58,7 +58,7 @@ void take(const value& v); void untake(const value& v); - void clear(); + void init(); std::size_t operator()(const value& v) const; std::size_t operator[](std::size_t i) const; @@ -86,7 +86,7 @@ /*! Generic histogram class over the set of values of type \c T. */ template <typename T> - struct histo_on_type : public histo_on_set< set_<T> > + struct histo_on_type : public histo_on_set< value::set_<T> > { histo_on_type(); }; @@ -127,9 +127,10 @@ template <typename S> void - histo_on_set<S>::clear() + histo_on_set<S>::init() { std::fill(h_.begin(), h_.end(), 0); + sum_ = 0; } template <typename S> @@ -180,6 +181,7 @@ { mln_viter(S) v(h.vset()); for_all(v) + if (h(v) != 0) ostr << v << ':' << h(v) << ' '; ostr << std::endl; return ostr; @@ -190,7 +192,7 @@ template <typename T> histo_on_type<T>::histo_on_type() - : histo_on_set< set_<T> >(set_<T>::the()) + : histo_on_set< value::set_<T> >(value::set_<T>::the()) { } @@ -198,9 +200,9 @@ # endif // ! MLN_INCLUDE_ONLY - } // end of namespace mln::value + } // end of namespace mln::accu } // end of namespace mln -#endif // ! MLN_VALUE_HISTO_HH +#endif // ! MLN_ACCU_HISTO_HH Index: mln/accu/median.hh --- mln/accu/median.hh (revision 0) +++ mln/accu/median.hh (working copy) @@ -25,21 +25,21 @@ // reasons why the executable file might be covered by the GNU General // Public License. -#ifndef MLN_VALUE_MEDIAN_HH -# define MLN_VALUE_MEDIAN_HH +#ifndef MLN_ACCU_MEDIAN_HH +# define MLN_ACCU_MEDIAN_HH -/*! \file mln/value/median.hh +/*! \file mln/accu/median.hh * * \brief Define FIXME */ -# include <mln/value/histo.hh> +# include <mln/accu/histo.hh> namespace mln { - namespace value + namespace accu { @@ -55,11 +55,15 @@ void take(const value& v); void untake(const value& v); - void clear(); + void init(); + + unsigned card() const { return h_.sum(); } operator value() const; value to_value() const; + const histo_on_set<S>& histo() const; + // FIXME: remove void debug__() const { @@ -71,20 +75,31 @@ protected: - histo_on_set<S> h_; + mutable histo_on_set<S> h_; const S& s_; // derived from h_ - std::size_t sum_minus_, sum_plus_; + mutable std::size_t sum_minus_, sum_plus_; - std::size_t i_; // the median index - value v_; // the median value + mutable bool valid_; + mutable std::size_t i_; // the median index + mutable value v_; // the median value // Auxiliary methods - void go_minus_(); - void go_plus_(); + void update_() const; + void go_minus_() const; + void go_plus_() const; }; + template <typename T> + struct median_on : public median< value::set_<T> > + { + median_on() + : median< value::set_<T> >(value::set_<T>::the()) + { + } + }; + # ifndef MLN_INCLUDE_ONLY @@ -94,7 +109,7 @@ : h_(s), s_(h_.vset()) { - clear(); + init(); } @@ -102,42 +117,15 @@ void median<S>::take(const value& v) { - // update h_ h_.take(v); - // particular case: - // current state was initialization - if (h_[i_] = 0) - { - // std::cout << "init!" << std::endl; - i_ = s_.index_of(v); - v_ = v; - return; - } - - // particular case: - // the median does not change - if (v = v_) - { - // std::cout << "no change!" << std::endl; - return; - } - - // general case: - if (v < v_) - { ++sum_minus_; - if (2 * sum_minus_ > h_.sum()) - go_minus_(); - } - else - // v > v_ - { + else if (v > v_) ++sum_plus_; - if (2 * sum_plus_ > h_.sum()) - go_plus_(); - } + + if (valid_) + valid_ = false; } @@ -146,34 +134,33 @@ median<S>::untake(const value& v) { mln_precondition(h_(v) != 0); - - // update h_ h_.untake(v); - // particular case: - // the only value has been removed - if (h_.sum() = 0) - { - clear(); - return; - } - - // general case: if (v < v_) - { --sum_minus_; - if (2 * sum_plus_ > h_.sum()) - go_plus_(); - } else if (v > v_) - { --sum_plus_; + + if (valid_) + valid_ = false; + } + + + template <typename S> + void + median<S>::update_() const + { + valid_ = true; + + if (h_.sum() = 0) + return; + if (2 * sum_minus_ > h_.sum()) go_minus_(); - } else - // v = v_ - { + if (2 * sum_plus_ > h_.sum()) + go_plus_(); + else if (h_[i_] = 0) { // go to the heaviest side @@ -182,21 +169,12 @@ else go_minus_(); // default when both sides are balanced } - else - { - if (2 * sum_plus_ > h_.sum()) - go_plus_(); - else if (2 * sum_minus_ > h_.sum()) - go_minus_(); - // else no change - } - } } template <typename S> void - median<S>::go_minus_() + median<S>::go_minus_() const { do { @@ -213,7 +191,7 @@ template <typename S> void - median<S>::go_plus_() + median<S>::go_plus_() const { do { @@ -230,29 +208,40 @@ template <typename S> void - median<S>::clear() + median<S>::init() { - h_.clear(); + h_.init(); sum_minus_ = 0; sum_plus_ = 0; i_ = (mln_max(value) - mln_min(value)) / 2; v_ = s_[i_]; + valid_ = true; } + template <typename S> median<S>::operator typename median<S>::value () const { - return v_; + return to_value(); } template <typename S> typename median<S>::value median<S>::to_value() const { + if (not valid_) + update_(); return v_; } template <typename S> + const histo_on_set<S>& + median<S>::histo() const + { + return h_; + } + + template <typename S> std::ostream& operator<<(std::ostream& ostr, const median<S>& m) { m.debug__(); @@ -262,9 +251,9 @@ # endif // ! MLN_INCLUDE_ONLY - } // end of namespace mln::value + } // end of namespace mln::accu } // end of namespace mln -#endif // ! MLN_VALUE_MEDIAN_HH +#endif // ! MLN_ACCU_MEDIAN_HH Index: mln/accu/median_alt.hh --- mln/accu/median_alt.hh (revision 0) +++ mln/accu/median_alt.hh (revision 0) @@ -0,0 +1,279 @@ +// Copyright (C) 2007 EPITA Research and Development Laboratory +// +// This file is part of the Olena Library. This library is free +// software; you can redistribute it and/or modify it under the terms +// of the GNU General Public License version 2 as published by the +// Free Software Foundation. +// +// This library 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 this library; see the file COPYING. If not, write to +// the Free Software Foundation, 51 Franklin Street, Fifth Floor, +// Boston, MA 02111-1307, USA. +// +// As a special exception, you may use this file as part of a free +// software library 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 MLN_ACCU_MEDIAN_ALT_HH +# define MLN_ACCU_MEDIAN_ALT_HH + +/*! \file mln/accu/median_alt.hh + * + * \brief Define FIXME + */ + +# include <mln/accu/histo.hh> + + +namespace mln +{ + + namespace accu + { + + + /*! Generic median_alt function based on histogram over a value set + * with type \c S. + */ + template <typename S> + struct median_alt + { + typedef mln_value(S) value; + + median_alt(const Value_Set<S>& s); + + void take(const value& v); + void untake(const value& v); + void init(); + + operator value() const; + value to_value() const; + + // FIXME: remove + void debug__() const + { + std::cout << " i = " << i_ + << " v = " << v_ + << " s = " << sum_minus_ << " ; " << h_[i_] << " ; " << sum_plus_ << " = " << h_.sum() + << std::endl; + } + + protected: + + histo_on_set<S> h_; + const S& s_; // derived from h_ + + std::size_t sum_minus_, sum_plus_; + + std::size_t i_; // the median index + value v_; // the median value + + // Auxiliary methods + void go_minus_(); + void go_plus_(); + }; + + + template <typename T> + struct median_alt_on : public median_alt< value::set_<T> > + { + median_alt_on() + : median_alt< value::set_<T> >(value::set_<T>::the()) + { + } + }; + + +# ifndef MLN_INCLUDE_ONLY + + + template <typename S> + median_alt<S>::median_alt(const Value_Set<S>& s) + : h_(s), + s_(h_.vset()) + { + init(); + } + + + template <typename S> + void + median_alt<S>::take(const value& v) + { + // update h_ + h_.take(v); + + // particular case: + // current state was initialization + if (h_[i_] = 0) + { + // std::cout << "init!" << std::endl; + i_ = s_.index_of(v); + v_ = v; + return; + } + + // particular case: + // the median does not change + if (v = v_) + { + // std::cout << "no change!" << std::endl; + return; + } + + // general case: + + if (v < v_) + { + ++sum_minus_; + if (2 * sum_minus_ > h_.sum()) + go_minus_(); + } + else + // v > v_ + { + ++sum_plus_; + if (2 * sum_plus_ > h_.sum()) + go_plus_(); + } + } + + + template <typename S> + void + median_alt<S>::untake(const value& v) + { + mln_precondition(h_(v) != 0); + + // update h_ + h_.untake(v); + + // particular case: + // the only value has been removed + if (h_.sum() = 0) + { + init(); + return; + } + + // general case: + if (v < v_) + { + --sum_minus_; + if (2 * sum_plus_ > h_.sum()) + go_plus_(); + } + else if (v > v_) + { + --sum_plus_; + if (2 * sum_minus_ > h_.sum()) + go_minus_(); + } + else + // v = v_ + { + if (h_[i_] = 0) + { + // go to the heaviest side + if (sum_plus_ > sum_minus_) + go_plus_(); + else + go_minus_(); // default when both sides are balanced + } + else + { + if (2 * sum_plus_ > h_.sum()) + go_plus_(); + else if (2 * sum_minus_ > h_.sum()) + go_minus_(); + // else no change + } + } + } + + + template <typename S> + void + median_alt<S>::go_minus_() + { + do + { + sum_plus_ += h_[i_]; + do + --i_; + while (h_[i_] = 0); + sum_minus_ -= h_[i_]; + } + while (2 * sum_minus_ > h_.sum()); + v_ = s_[i_]; + } + + + template <typename S> + void + median_alt<S>::go_plus_() + { + do + { + sum_minus_ += h_[i_]; + do + ++i_; + while (h_[i_] = 0); + sum_plus_ -= h_[i_]; + } + while (2 * sum_plus_ > h_.sum()); + v_ = s_[i_]; + } + + + template <typename S> + void + median_alt<S>::init() + { + h_.init(); + sum_minus_ = 0; + sum_plus_ = 0; + i_ = (mln_max(value) - mln_min(value)) / 2; + v_ = s_[i_]; + } + + template <typename S> + median_alt<S>::operator typename median_alt<S>::value () const + { + return v_; + } + + template <typename S> + typename median_alt<S>::value + median_alt<S>::to_value() const + { + return v_; + } + + template <typename S> + std::ostream& operator<<(std::ostream& ostr, const median_alt<S>& m) + { + m.debug__(); + return ostr << m.to_value(); + } + + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::accu + +} // end of namespace mln + + +#endif // ! MLN_ACCU_MEDIAN_ALT_HH