* green/mln/clustering/kmean_rgb.hh: New library source code. * green/demo/clustering/kmean_rgb/Makefile: New Makefile. * green/demo/clustering/kmean_rgb/kmean_rgb.cc: New demo. code. --- milena/sandbox/ChangeLog | 8 + .../histo2d => clustering/kmean_rgb}/Makefile.am | 0 .../green/demo/clustering/kmean_rgb/kmean_rgb.cc | 91 ++ milena/sandbox/green/mln/clustering/kmean_rgb.hh | 973 ++++++++++++++++++++ 4 files changed, 1072 insertions(+), 0 deletions(-) copy milena/sandbox/green/demo/{accu/stat/histo2d => clustering/kmean_rgb}/Makefile.am (100%) create mode 100644 milena/sandbox/green/demo/clustering/kmean_rgb/kmean_rgb.cc create mode 100644 milena/sandbox/green/mln/clustering/kmean_rgb.hh diff --git a/milena/sandbox/ChangeLog b/milena/sandbox/ChangeLog index a2cc751..b3029e2 100644 --- a/milena/sandbox/ChangeLog +++ b/milena/sandbox/ChangeLog @@ -1,5 +1,13 @@ 2009-12-02 Yann Jacquelet <jacquelet@lrde.epita.fr> + Transform kmean object in a big function and then a canvas. + + * green/mln/clustering/kmean_rgb.hh: New library source code. + * green/demo/clustering/kmean_rgb/Makefile: New Makefile. + * green/demo/clustering/kmean_rgb/kmean_rgb.cc: New demo. code. + +2009-12-02 Yann Jacquelet <jacquelet@lrde.epita.fr> + Benchmark on distance, preliminary work before optimizing kmean. * green/bench/transform/distance/Makefile: New Makefile. diff --git a/milena/sandbox/green/demo/accu/stat/histo2d/Makefile.am b/milena/sandbox/green/demo/clustering/kmean_rgb/Makefile.am similarity index 100% copy from milena/sandbox/green/demo/accu/stat/histo2d/Makefile.am copy to milena/sandbox/green/demo/clustering/kmean_rgb/Makefile.am diff --git a/milena/sandbox/green/demo/clustering/kmean_rgb/kmean_rgb.cc b/milena/sandbox/green/demo/clustering/kmean_rgb/kmean_rgb.cc new file mode 100644 index 0000000..36641b6 --- /dev/null +++ b/milena/sandbox/green/demo/clustering/kmean_rgb/kmean_rgb.cc @@ -0,0 +1,91 @@ +// DEMO ON KMEAN_RGB + +#include <mln/clustering/kmean_rgb.hh> + +#include <iostream> +#include <sstream> + +#include <mln/img_path.hh> + +#include <mln/core/macros.hh> +#include <mln/core/image/image2d.hh> +//#include <mln/core/image/dmorph/image_if.hh> + +#include <mln/data/transform.hh> +#include <mln/fun/v2v/rgb8_to_rgbn.hh> + +#include <mln/io/ppm/load.hh> +//#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> +//#include <mln/io/ppm/save.hh> +//#include <mln/io/plot/save_image_sh.hh> + +//#include <mln/pw/value.hh> + +//#include <mln/value/int_u8.hh> +#include <mln/value/label_8.hh> +#include <mln/value/rgb8.hh> + +int main() +{ + const std::string& image = OLENA_IMG_PATH"/house.ppm"; + const unsigned k_center = 3; + const unsigned n_times = 10; + const unsigned watch_dog = 10; + + typedef mln::value::label_8 t_lbl8; + typedef mln::value::rgb8 t_rgb8; + typedef mln::value::rgb<5> t_rgb5; + typedef mln::image2d<t_rgb8> t_image2d_rgb8; + typedef mln::image2d<t_rgb5> t_image2d_rgb5; + typedef mln::image2d<t_lbl8> t_image2d_lbl8; + typedef mln::fun::v2v::rgb8_to_rgbn<5> t_rgb8_to_rgb5; + + t_image2d_rgb8 house_rgb8; + t_image2d_rgb5 house_rgb5; + t_image2d_lbl8 house_lbl8; + + mln::io::ppm::load(house_rgb8, image.c_str()); + house_rgb5 = mln::data::transform(house_rgb8, t_rgb8_to_rgb5()); + + mln::trace::quiet = false; + + house_lbl8 = mln::clustering::kmean_rgb<double,5>(house_rgb5, + k_center, + watch_dog, + n_times); + mln::trace::quiet = true; + + mln::io::pgm::save(house_lbl8, "label.pgm"); + /* + + + + t_kmean kmean(house_rgb5, k_center, watch_dog, n_times); + + + + //kmean.launch_one_time(); + kmean.launch_n_times(); + + // Not safe because we don't test kmean.is_valid() + + t_kmean::t_color_dbg color_img = kmean.get_color_dbg(); + t_kmean::t_mean_dbg mean_img = kmean.get_mean_dbg(); + t_kmean::t_label_dbg label_img = kmean.get_label_dbg(); + t_kmean::t_variance_cnv variance_cnv = kmean.get_variance_cnv(); + t_kmean::t_mean_cnv mean_cnv = kmean.get_mean_cnv(); + + mln::io::ppm::save(mean_img, "mean.ppm"); + mln::io::ppm::save(color_img, "color.ppm"); + mln::io::pgm::save(label_img, "label.pgm"); + + mln::io::plot::save_image_sh(mean_img, "mean.sh"); + mln::io::plot::save_image_sh(mean_cnv, "mean_cnv.sh"); + mln::io::plot::save_image_sh(variance_cnv, "variance_cnv.sh"); + + */ + + return 0; +} + diff --git a/milena/sandbox/green/mln/clustering/kmean_rgb.hh b/milena/sandbox/green/mln/clustering/kmean_rgb.hh new file mode 100644 index 0000000..253745b --- /dev/null +++ b/milena/sandbox/green/mln/clustering/kmean_rgb.hh @@ -0,0 +1,973 @@ +// Copyright (C) 2008,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 MLN_CLUSTERING_KMEAN_RGB_HH +# define MLN_CLUSTERING_KMEAN_RGB_HH + +/// \file +/// +/// \brief Implements the optimized kmean algorithm. +/// +/// This algorithm is optimized in the way it proceeds directly with +/// the rgb values inspite of the pixel attribute. The +/// algorithm is independant from the image dimension. But, we have to +/// compute one time the histogram. In fact, we move a recurrent cost +/// to a fix cost in the complexity. This version is adapted to +/// image with small quantification. + +/// APLATISSEMENT DES KMEAN3D + +# include <limits.h> +# include <iostream> + +# include <mln/accu/stat/histo3d_rgb.hh> + +# include <mln/algebra/vec.hh> + +# include <mln/core/concept/image.hh> +# include <mln/core/contract.hh> +# include <mln/core/image/image2d.hh> +# include <mln/core/macros.hh> + +# include <mln/data/compute.hh> +# include <mln/data/fill.hh> +# include <mln/data/transform.hh> + +# include <mln/debug/println.hh> + +# include <mln/io/ppm/save.hh> +# include <mln/io/pgm/save.hh> + +# include <mln/labeling/colorize.hh> +# include <mln/labeling/mean_values.hh> + +# include <mln/literal/zero.hh> +# include <mln/literal/one.hh> + +# include <mln/math/min.hh> +# include <mln/math/sqr.hh> + +# include <mln/norm/l2.hh> + +# include <mln/opt/at.hh> + +# include <mln/trace/entering.hh> +# include <mln/trace/exiting.hh> + +# include <mln/trait/value_.hh> + +# include <mln/util/array.hh> + +# include <mln/value/int_u.hh> +# include <mln/value/rgb8.hh> +# include <mln/value/label_8.hh> + + +//-------------------------------------------------------------------------- +// CODE APLATI +//-------------------------------------------------------------------------- + + +namespace mln +{ + + namespace clustering + { + + template <typename T, unsigned n, typename I> + inline + image2d<value::label_8> + kmean_rgb(const Image<I>& point, + const unsigned k_center, + const unsigned watch_dog, + const unsigned n_times); + + } // end of namespace mln::clustering + + namespace clustering + { + + +# ifndef MLN_INCLUDE_ONLY + + //-------------------------------------------------------------------------- + // Internal. + //-------------------------------------------------------------------------- + + namespace internal + { + + //------------------------------------------------------------------------ + // Debugging tools + //------------------------------------------------------------------------ + + /* + template <typename T, unsigned n> + inline + void kmean3d_a<T,n>::build_label_dbg() + { + trace::entering("mln::clustering::kmean3d_a::build_label_dbg"); + + mln_piter(t_point_img) pi(_point.domain()); + mln_piter(t_label_dbg) po(_label_dbg.domain()); + + for_all_2(pi, po) + { + t_value val = _point(pi); + t_label grp = _group(point3d(val.blue(), val.red(), val.green())); + + // As label zero has got a particular semantic, the first label is one + _label_dbg(po) = ++grp; + } + + trace::exiting("mln::clustering::kmean3d_a::build_label_dbg"); + } + + template <typename T, unsigned n> + inline + void kmean3d_a<T,n>::build_mean_dbg() + { + trace::entering("mln::clustering::kmean3d_a::build_mean_dbg"); + + mln_piter(t_mean_dbg) p(_mean_dbg.domain()); + + for_all(p) + { + _mean_dbg(p).red() = static_cast<unsigned>(_mean[_label_dbg(p)][0]); + _mean_dbg(p).green() = static_cast<unsigned>(_mean[_label_dbg(p)][1]); + _mean_dbg(p).blue() = static_cast<unsigned>(_mean[_label_dbg(p)][2]); + } + + trace::exiting("mln::clustering::kmean3d_a::build_mean_dbg"); + } + + + template <typename T, unsigned n> + inline + void kmean3d_a<T,n>::build_all_dbg() + { + trace::entering("mln::clustering::kmean3d_a::build_all_dbg"); + build_label_dbg(); + //build_mean_dbg(); + _mean_dbg = labeling::mean_values(_point, _label_dbg, _k_center); + _color_dbg = labeling::colorize(value::rgb8(), _label_dbg); + + trace::exiting("mln::clustering::kmean3d_a::build_all_dbg"); + } + + template <typename T, unsigned n> + inline + void kmean3d_a<T,n>::update_cnv() + { + trace::entering("mln::clustering::kmean3d_a::update_cnv"); + + _variance_cnv[_current_launching](point1d(_current_step)) + = _within_variance; + + mln_eiter(t_mean_img) l(_mean); + + for_all(l) + { + _mean_cnv[l.index_()][_current_launching](point1d(_current_step)) + = _mean[l.index_()]; + } + + trace::exiting("mln::clustering::kmean3d_a::update_cnv"); + } + + template <typename T, unsigned n> + inline + void kmean3d_a<T,n>::finalize_cnv() + { + trace::entering("mln::clustering::kmean3d_a::finalize_cnv"); + + // saturate the curv with the within variance + for (unsigned i = _current_step; i < _watch_dog; ++i) + _variance_cnv[_current_launching](point1d(i)) = _within_variance; + + for (unsigned i = _current_step; i < _watch_dog; ++i) + { + mln_eiter(t_mean_img) l(_mean); + + for_all(l) + { + _mean_cnv[l.index_()][_current_launching](point1d(i)) + = _mean[l.index_()]; + } + } + + trace::exiting("mln::clustering::kmean3d_a::finalize_cnv"); + } + + + + + //-------------------------------------------------------------------------- + // Printing temporary results + //-------------------------------------------------------------------------- + + template <typename T, unsigned n> + inline + void kmean3d_a<T,n>::print_mean() + { + trace::entering("mln::clustering::kmean3d_a::print_mean"); + + mln_eiter(t_mean_img) l(_mean); + + for_all(l) + { + std::cout << "mean(" << l.index_(); + std::cout << ") = [r=" << _mean[l.index_()][0]; + std::cout << ", g=" << _mean[l.index_()][1]; + std::cout << ", b=" << _mean[l.index_()][2]; + std::cout << "]" << std::endl; + } + + trace::exiting("mln::clustering::kmean3d_a::print_mean"); + } + + template <typename T, unsigned n> + inline + void kmean3d_a<T,n>::print_number() + { + trace::entering("mln::clustering::kmean3d_a::print_number"); + + mln_eiter(t_number_img) l(_number); + + for_all(l) + { + std::cout << "number(" << l.index_(); + std::cout << ") = " << _number[l.index_()]; + std::cout << std::endl; + } + + trace::exiting("mln::clustering::kmean3d_a::print_number"); + } + + template <typename T, unsigned n> + inline + void kmean3d_a<T,n>::print_variance() + { + trace::entering("mln::clustering::kmean3d_a::print_variance"); + + mln_eiter(t_variance_img) l(_number); + + for_all(l) + { + std::cout << "variance(" << l.index_(); + std::cout << ") = " << _variance[l.index_()]; + std::cout << std::endl; + } + + trace::exiting("mln::clustering::kmean3d_a::print_variance"); + } + + template <typename T, unsigned n> + inline + void kmean3d_a<T,n>::print_histo() + { + trace::entering("mln::clustering::kmean3d_a::print_histo"); + + mln_piter(t_histo_img) rgb(_histo.domain()); + + for_all(rgb) + { + if (0 < _histo(rgb)) + { + std::cout << "histo(r=" << rgb.row(); + std::cout << ", g=" << rgb.col(); + std::cout << ", b=" << rgb.sli(); + std::cout << ")= " << _histo(rgb); + std::cout << std::endl; + } + } + + trace::exiting("mln::clustering::kmean3d_a::print_histo"); + } + + template <typename T, unsigned n> + inline + void kmean3d_a<T,n>::print_group() + { + trace::entering("mln::clustering::kmean3d_a::print_group"); + + mln_piter(t_group_img) rgb(_group.domain()); + + for_all(rgb) + { + if (0 < _histo(rgb)) + { + std::cout << "group(r=" << rgb.row(); + std::cout << ", g=" << rgb.col(); + std::cout << ", b=" << rgb.sli(); + std::cout << ")= " << _group(rgb); + std::cout << std::endl; + } + } + + trace::exiting("mln::clustering::kmean3d_a::print_group"); + } + + template <typename T, unsigned n> + inline + void kmean3d_a<T,n>::print_distance() + { + trace::entering("mln::clustering::kmean3d_a::print_distance"); + + mln_eiter(t_distance_img) l(_distance); + + for_all(l) + { + mln_piter(t_distance_val) rgb(_distance[l.index_()].domain()); + + for_all(rgb) + { + if (0 < _histo(rgb)) + { + std::cout << "distance(l=" << l.index_(); + std::cout << ",r=" << rgb.row(); + std::cout << ", g=" << rgb.col(); + std::cout << ", b=" << rgb.sli(); + std::cout << ")= " << _distance[l.index_()](rgb); + std::cout << std::endl; + } + } + } + + trace::exiting("mln::clustering::kmean3d_a::print_distance"); + } + + template <typename T, unsigned n> + inline + void kmean3d_a<T,n>::print_point() + { + trace::entering("mln::clustering::kmean3d_a::print_point"); + + mln_piter(t_point_img) p(_point.domain()); + + for_all(p) + { + std::cout << "point(r=" << p.row(); + std::cout << ", c=" << p.col(); + std::cout << ")= " << _point(p); + std::cout << std::endl; + } + + trace::exiting("mln::clustering::kmean3d_a::print_point"); + } + + + + template <typename T, unsigned n> + inline + void rgb_rand_init(t_mean_img mean) + { + typedef value::rgb<n> t_value; + typedef mln_trait_value_comp(t_value,0) t_value_comp0; + typedef mln_trait_value_comp(t_value,1) t_value_comp1; + typedef mln_trait_value_comp(t_value,2) t_value_comp2; + typedef algebra::vec<3,T> t_result3d; + typedef util::array<t_result3d> t_mean_img; + + t_value_comp0 min_comp0 = mln_min(t_value_comp0); + t_value_comp0 max_comp0 = mln_max(t_value_comp0); + t_value_comp1 min_comp1 = mln_min(t_value_comp1); + t_value_comp1 max_comp1 = mln_max(t_value_comp1); + t_value_comp2 min_comp2 = mln_min(t_value_comp2); + t_value_comp2 max_comp2 = mln_max(t_value_comp2); + mln_eiter(t_mean_img) l(mean); + + for_all(l) + { + mean[l.index_()][0] = (rand() % (max_comp0 - min_comp0)) + min_comp0; + mean[l.index_()][1] = (rand() % (max_comp1 - min_comp1)) + min_comp1; + mean[l.index_()][2] = (rand() % (max_comp2 - min_comp2)) + min_comp2; + } + + return mean; + } + + */ + + } // end of namespace mln::clustering::internal + + + //-------------------------------------------------------------------------- + // Impl. + //-------------------------------------------------------------------------- + + namespace impl + { + + //------------------------------------------------------------------------ + // kmean_image2d_rgb(const t_point_img& point, + // const unsigned k_center, + // const unsigned watch_dog = 10, + // const unsigned n_times = 10) + //------------------------------------------------------------------------ + + template <unsigned n> + struct rgbn_to_lbl8 : Function_v2v< rgbn_to_lbl8<n> > + { + typedef value::rgb<n> argument; + typedef value::label_8 result; + typedef value::label_8 t_label; + typedef image3d<t_label> t_group_img; + + t_group_img _group; + + rgbn_to_lbl8(t_group_img group) : _group(group) {} + + result operator()(const argument& c) const + { + value::label_8 tmp = opt::at(_group, c.blue(), c.red(), c.green()); + + // FIXME WHY DO WE NOT USE +1 + return ++tmp; + } + }; + + template <typename T, unsigned n> + struct rgb_to_dist : Function_v2v< rgb_to_dist<T,n> > + { + typedef value::rgb<n> argument; + typedef T result; + typedef T t_result1d; + typedef algebra::vec<3,T> t_result3d; + typedef image3d<unsigned> t_histo_img; + + t_result3d _mean; + t_histo_img _histo; + + rgb_to_dist(t_result3d mean, t_histo_img histo) : _mean(mean), + _histo(histo) {} + + result operator()(const argument& c) const + { + t_result1d diff2_row = math::sqr(c.row() - _mean[0]); + t_result1d diff2_col = math::sqr(c.col() - _mean[1]); + t_result1d diff2_sli = math::sqr(c.sli() - _mean[2]); + t_result1d tmp = _histo(c)*(diff2_row + diff2_col + diff2_sli); + + return tmp; + } + }; + + template <typename T, unsigned n> + inline + image2d<value::label_8> + kmean_image2d_rgb(const image2d< value::rgb<n> >& point, + const unsigned k_center, + const unsigned watch_dog = 10, + const unsigned n_times = 10) + { + trace::entering("mln::clustering::impl::kmean_image2d_rgb"); + trace::quiet = true; + // BEGIN TYPEDEF + typedef value::rgb<8> t_rgb; + typedef value::label<8> t_label; + typedef value::rgb<n> t_value; + typedef mln_trait_value_comp(t_value,0) t_value_comp0; + typedef mln_trait_value_comp(t_value,1) t_value_comp1; + typedef mln_trait_value_comp(t_value,2) t_value_comp2; + typedef T t_result1d; + typedef algebra::vec<3,T> t_result3d; + + typedef image2d<t_value> t_point_img; + typedef image3d<unsigned> t_histo_img; + typedef util::array<t_result1d> t_number_img; + typedef util::array<t_result3d> t_mean_img; + typedef util::array<t_result1d> t_variance_img; + + typedef image3d<t_label> t_group_img; + typedef image3d<t_result1d> t_distance_val; + typedef util::array<t_distance_val> t_distance_img; + + typedef image2d<t_label> t_label_dbg; + typedef image2d<t_rgb> t_color_dbg; + typedef image2d<t_value> t_mean_dbg; + + typedef image1d<t_result3d> t_mean_val; + typedef util::array<t_mean_val> t_mean_set; + typedef util::array<t_mean_set> t_mean_cnv; + typedef image1d<t_result1d> t_variance_val; + typedef util::array<t_variance_val> t_variance_cnv; + // END TYPEDEF + + // BEGIN INITIALISATION + mln_precondition(point.is_valid()); + + static const unsigned _N_TRIES = 3; + + typedef accu::meta::stat::histo3d_rgb t_histo3d_rgb; + + t_result1d _within_variance; + + unsigned _k_center = k_center; + unsigned _watch_dog = watch_dog; + unsigned _n_times = n_times; + t_point_img _point = point; + + // HISTOGRAM INIT + t_histo_img _histo = data::compute(t_histo3d_rgb(), + _point); + + // CENTER STATS INIT + t_number_img _number; + t_mean_img _mean; + t_variance_img _variance; + + for (unsigned i = 0; i < _k_center; ++i) + { + _number.append(literal::zero); + _mean.append(literal::zero); + _variance.append(literal::zero); + } + + + unsigned _current_step = 0; + unsigned _current_launching = 0; + bool _is_number_valid = false; + + unsigned _launching_min; + t_result1d _variance_min; + t_mean_img _mean_min; + + + + t_group_img _group; + t_distance_img _distance; + + + t_label_dbg _label_dbg; + t_color_dbg _color_dbg; + t_mean_dbg _mean_dbg; + + + t_mean_cnv _mean_cnv; + t_variance_cnv _variance_cnv; + + + + + _group.init_(box3d(point3d(mln_min(t_value_comp2), + mln_min(t_value_comp0), + mln_min(t_value_comp1)), + point3d(mln_max(t_value_comp2), + mln_max(t_value_comp0), + mln_max(t_value_comp1)))); + + for (unsigned i = 0; i < _k_center; ++i) + { + t_distance_val img(box3d(point3d(mln_min(t_value_comp2), + mln_min(t_value_comp0), + mln_min(t_value_comp1)), + point3d(mln_max(t_value_comp2), + mln_max(t_value_comp0), + mln_max(t_value_comp1)))); + + _distance.append(img); + } + + // Debugging, calibrating and testing + initialize(_label_dbg, _point); + initialize(_color_dbg, _point); + initialize(_mean_dbg, _point); + + // Observing the convergence + + for (unsigned i = 0; i < _n_times; ++i) + { + t_variance_val img(box1d(point1d(0), point1d(_watch_dog-1))); + + data::fill(img, literal::zero); + + _variance_cnv.append(img); + } + + for (unsigned i = 0; i < _k_center; ++i) + { + t_mean_set mean_set; + + for (unsigned j = 0; j < _n_times; ++j) + { + t_mean_val img(box1d(point1d(0), point1d(_watch_dog-1))); + + data::fill(img, literal::zero); + + mean_set.append(img); + } + + _mean_cnv.append(mean_set); + } + // END INITIALISATION + + // BEGIN LOOP N TIMES + { + unsigned tries = 0; + _variance_min = mln_max(t_result1d); + _current_launching = 0; + + while (_current_launching < _n_times) + { + // BEGIN LAUNCH ONE TIME + trace::quiet = false; + trace::entering("Launch one time"); + trace::quiet = true; + { + t_result1d old_variance = mln_max(t_result1d); + _within_variance = mln_max(t_result1d); + _current_step = 0; + + // BEGIN INIT_MEAN + trace::quiet = false; + trace::entering("init mean"); + trace::quiet = true; + { + t_value_comp0 min_comp0 = mln_min(t_value_comp0); + t_value_comp0 max_comp0 = mln_max(t_value_comp0); + t_value_comp1 min_comp1 = mln_min(t_value_comp1); + t_value_comp1 max_comp1 = mln_max(t_value_comp1); + t_value_comp2 min_comp2 = mln_min(t_value_comp2); + t_value_comp2 max_comp2 = mln_max(t_value_comp2); + mln_eiter(t_mean_img) l(_mean); + + for_all(l) + { + _mean[l.index_()][0]=(rand()%(max_comp0-min_comp0))+min_comp0; + _mean[l.index_()][1]=(rand()%(max_comp1-min_comp1))+min_comp1; + _mean[l.index_()][2]=(rand()%(max_comp2-min_comp2))+min_comp2; + } + } + trace::quiet = false; + trace::exiting("init mean"); + trace::quiet = true; + // END INIT MEAN + + + // UPDATE DISTANCE + trace::quiet = false; + trace::entering("update distance"); + trace::quiet = true; + + for (unsigned i = 0; i < _k_center; ++i) + { + + // _distance[i] = data::transform(_histo, + // rgb_to_dist<T,n>(_mean[i], + // _histo)); + + mln_piter(t_distance_val) d(_distance[i].domain()); + + for_all(d) + { + t_result1d diff2_row = math::sqr(d.row() - _mean[i][0]); + t_result1d diff2_col = math::sqr(d.col() - _mean[i][1]); + t_result1d diff2_sli = math::sqr(d.sli() - _mean[i][2]); + _distance[i](d) = _histo(d)* + (diff2_row + diff2_col + diff2_sli); + } + } + + trace::quiet = false; + trace::exiting("update distance"); + trace::quiet = true; + // END UPDATE DISTANCE + + do + { + old_variance = _within_variance; + + // BEGIN UPDATE GROUP + trace::quiet = false; + trace::entering("update group"); + trace::quiet = true; + { + mln_piter(t_group_img) rgb(_group.domain()); + + for_all(rgb) + { + mln_eiter(t_distance_img) l(_distance); + t_result1d min = mln_max(t_result1d); + t_label label = mln_max(t_label); + + for_all(l) + { + if (min > _distance[l.index_()](rgb)) + { + min = _distance[l.index_()](rgb); + label = l.index_(); + } + } + + _group(rgb) = label; + } + + } + trace::quiet = false; + trace::exiting("update group"); + trace::quiet = true; + // END UPDATE GROUP + + // BEGIN UPDATE MEAN + trace::quiet = false; + trace::entering("update mean"); + trace::quiet = true; + { + mln_eiter(t_number_img) en(_number); + mln_eiter(t_mean_img) em(_mean); + + for_all_2(en,em) + { + _number[en.index_()] = literal::zero; + _mean[em.index_()] = literal::zero; + } + + mln_piter(t_group_img) rgb(_group.domain()); + + for_all(rgb) + { + _mean[_group(rgb)][0] += rgb.row() * _histo(rgb); + _mean[_group(rgb)][1] += rgb.col() * _histo(rgb); + _mean[_group(rgb)][2] += rgb.sli() * _histo(rgb); + _number(_group(rgb)) += _histo(rgb); + } + + mln_eiter(t_mean_img) l(_mean); + + for_all(l) + { + _is_number_valid = (0 != _number[l.index_()]); + + if (!_is_number_valid) + break; + + _mean[l.index_()] /= _number[l.index_()]; + } + } + trace::quiet = false; + trace::exiting("update mean"); + trace::quiet = true; + // END UPDATE MEAN + + + // Stopping Nan propagation + if (!_is_number_valid) + break; + + // UPDATE DISTANCE + trace::quiet = false; + trace::entering("update distance"); + trace::quiet = true; + + for (unsigned i = 0; i < _k_center; ++i) + { + mln_piter(t_distance_val) d(_distance[i].domain()); + + for_all(d) + { + // the square distance + t_result1d diff2_row = math::sqr(d.row() - _mean[i][0]); + t_result1d diff2_col = math::sqr(d.col() - _mean[i][1]); + t_result1d diff2_sli = math::sqr(d.sli() - _mean[i][2]); + _distance[i](d) = _histo(d)* + (diff2_row + diff2_col + diff2_sli); + } + } + trace::quiet = false; + trace::exiting("update distance"); + trace::quiet = true; + // END UPDATE DISTANCE + + // BEGIN UPDATE VARIANCE + trace::quiet = false; + trace::entering("update variance"); + trace::quiet = true; + { + _within_variance = literal::zero; + mln_eiter(t_variance_img) l(_variance); + + for_all(l) + { + _variance[l.index_()] = literal::zero; + + mln_piter(t_group_img) rgb(_group.domain()); + + for_all(rgb) + { + if (l.index_() == _group(rgb)) + _variance[l.index_()] += _distance[l.index_()](rgb); + } + + _within_variance += _variance[l.index_()]; + } + + } + trace::quiet = false; + trace::exiting("update variance"); + trace::quiet = true; + // END UPDATE VARIANCE + + //update_cnv(); + + ++_current_step; + } + while (_current_step < _watch_dog && + _within_variance < old_variance); + + //finalize_cnv(); + //build_all_dbg(); + } + trace::quiet = false; + trace::exiting("Launch one time"); + trace::quiet = true; + // END LAUNCH ONE TIME + + if ((_is_number_valid && (_current_step < _watch_dog))|| + _N_TRIES < tries) + { + if (_within_variance < _variance_min) + { + _variance_min = _within_variance; + _mean_min = _mean; + _launching_min = _current_launching; + } + + // Reinitialize the number of echecs possible + tries = 0; + + //std::cout << "_current_launching : " << _current_launching + // << std::endl; + + //std::cout << "within_variance[" << _current_launching << "] = " + // << _within_variance << std::endl; + + ++_current_launching; + } + else + ++tries; + } + + //Debugging code + //build_all_dbg(); + + } + // END LOOP N TIMES + + // BEGIN BUILD LABEL IMAGE + _label_dbg = data::transform(_point, rgbn_to_lbl8<n>(_group)); + +// { +// mln_piter(t_point_img) pi(_point.domain()); +// mln_piter(t_label_dbg) po(_label_dbg.domain()); + +// for_all_2(pi, po) +// { +// t_value val = _point(pi); +// t_label grp = _group(point3d(val.blue(),val.red(),val.green())); + +// _label_dbg(po) = ++grp; +// } +// } + + // END BUILD LABEL IMAGE + trace::quiet = false; + trace::exiting("mln::clustering::impl::kmean_image2d_rgb"); + + return _label_dbg; + + } + + } // end of namespace mln::clustering::impl + + + + + + //-------------------------------------------------------------------------- + // Internal. + //-------------------------------------------------------------------------- + + namespace internal + { + + template <typename T, unsigned n> + inline + image2d<value::label_8> + kmean_rgb_dispatch(const image2d< value::rgb<n> >& img, + const unsigned k_center, + const unsigned watch_dog, + const unsigned n_times) + { + return impl::kmean_image2d_rgb<T,n>(img, k_center, watch_dog, n_times); + } + + + template <typename T, unsigned n, typename I> + inline + image2d< value::label_8> + kmean_rgb_dispatch(const Image<I>& img, + const unsigned k_center, + const unsigned watch_dog, + const unsigned n_times) + { + return kmean_rgb_dispatch<T,n>(exact(img),k_center,watch_dog,n_times); + } + + + } // end of namespace mln::clustering::internal + + + //-------------------------------------------------------------------------- + // Facade. + //-------------------------------------------------------------------------- + + template <typename T, unsigned n, typename I> + inline + image2d<value::label_8> + kmean_rgb(const Image<I>& point, + const unsigned k_center, + const unsigned watch_dog, + const unsigned n_times) + { + trace::entering("mln::clustering::kmean_rgb"); + + image2d<value::label_8> tmp = internal::kmean_rgb_dispatch<T,n>(point, + k_center, + watch_dog, + n_times); + trace::exiting("mln::clustering::kmean_rgb"); + + return tmp; + } + + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::clustering + +} // end of namespace mln + +#endif // ! MLN_CLUSTERING_KMEAN_RGB_HH -- 1.5.6.5