April 2008 - Olena-patches

1835: Fix a bug in basic_najman.
by Alexandre Abraham 03 Apr '08

03 Apr '08

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Alexandre Abraham <abraham(a)lrde.epita.fr> Fix a bug in basic_najman. * sandbox/abraham/morpho/basic_najman.hh: (mln::morpho::basic_najman<I,N>::BuildComponentTree): Fix a bug, the algorithm is still bugged but now gives a nice tree. * sandbox/abraham/morpho/test.cc: Add picture inversion for test purpose. basic_najman.hh | 103 +++++++++++++++++++++++++++++++++++--------------------- test.cc | 10 ++++- 2 files changed, 75 insertions(+), 38 deletions(-) Index: sandbox/abraham/morpho/basic_najman.hh --- sandbox/abraham/morpho/basic_najman.hh (revision 1834) +++ sandbox/abraham/morpho/basic_najman.hh (working copy) @@ -37,13 +37,13 @@ const Image& ima; const Neighborhood<N>& nbh; - mln_ch_value(I, psite) par; + mln_ch_value(I, psite) Par_node; mln_ch_value(I, psite) Par_tree; // Par_node is handled by par (from ap_maxtree) mln_ch_value(I, int) Rnk_tree; mln_ch_value(I, int) Rnk_node; - mln_ch_value(I, psite) lowestNode; + mln_ch_value(I, psite) subtreeRoot; mln_ch_value(I, node) nodes; mln_ch_value(I, bool) isproc; psite Root; @@ -54,11 +54,11 @@ const Neighborhood<N>& n) : ima(i), nbh(n), - par(exact(i).domain(), exact(i).border()), + Par_node(exact(i).domain(), exact(i).border()), Par_tree(exact(i).domain(), exact(i).border()), Rnk_tree(exact(i).domain(), exact(i).border()), Rnk_node(exact(i).domain(), exact(i).border()), - lowestNode(exact(i).domain(), exact(i).border()), + subtreeRoot(exact(i).domain(), exact(i).border()), nodes(exact(i).domain(), exact(i).border()), isproc(exact(i).domain(), exact(i).border()) { @@ -73,7 +73,7 @@ void MakeSet_node(psite x) { - par(x) = x; // was Par_node(x) = x; + Par_node(x) = x; // was Par_node(x) = x; Rnk_node(x) = 0; } @@ -91,29 +91,29 @@ y = memo; } - bool is_root(const psite& x) const - { - return par(x) == x; - } +// bool is_root(const psite& x) const +// { +// return Par_node(x) == x; +// } - bool is_level_root(const psite& x) const - { - return is_root(x) || exact(ima)(par(x)) < exact(ima)(x); - } +// bool is_level_root(const psite& x) const +// { +// return is_root(x) || exact(ima)(Par_node(x)) < exact(ima)(x); +// } - psite find_level_root(const psite& x) - { - if (is_level_root(x)) - return x; - else - return par(x) = find_level_root(par(x)); - } +// psite find_level_root(const psite& x) +// { +// if (is_level_root(x)) +// return x; +// else +// return Par_node(x) = find_level_root(Par_node(x)); +// } psite Find_node(psite x) { - if (par(x) != x) - par(x) = Find_node(par(x)); - return par(x); + if (Par_node(x) != x) + Par_node(x) = Find_node(Par_node(x)); + return Par_node(x); } void go() @@ -135,8 +135,8 @@ psite p = S[ip]; MakeSet_node(p); MakeSet_tree(p); - // if (lowestNode.hold(p)) - lowestNode(p) = p; + // if (subtreeRoot.hold(p)) + subtreeRoot(p) = p; // if (nodes.hold(p)) nodes(p) = MakeNode(exact(ima)(p)); } @@ -148,22 +148,32 @@ for (int ip = 0; ip < int(S.npoints()); ++ip) { psite p = S[ip]; - isproc(p) = true; + // isproc(p) = true; psite curTree = Find_tree(p); - psite curNode = Find_node(lowestNode(curTree)); + psite curNode = Find_node(subtreeRoot(curTree)); mln_niter(N) q(nbh, p); for_all(q) - if (exact(ima).has(q) and isproc(q) and exact(ima)(q) >= exact(ima)(p)) // f[q] >= f[p] parce qu'on prend l'ordre de l'histogramme + if (exact(ima).has(q) and isproc(q) and exact(ima)(q) <= exact(ima)(p)) { psite adjTree = Find_tree(q); - psite adjNode = Find_node(lowestNode(adjTree)); + psite adjNode = Find_node(subtreeRoot(adjTree)); if (curNode != adjNode) { + std::cout << "curNode != adjNode; " << nodes(curNode).level << " vs. " << nodes(adjNode).level << std::endl; if (nodes(curNode).level == nodes(adjNode).level) - curNode = MergeNode(adjNode, curNode); + { + // curNode = MergeNode(adjNode, curNode); + psite tmpNode = MergeNode(adjNode, curNode); + /*if (tmpNode == curNode) + nodes(curNode).addChildren(nodes(adjNode)); + else + nodes(adjNode).addChildren(nodes(curNode)); + nodes(adjNode) = nodes(curNode);*/ + curNode = tmpNode; + } else { // we have nodes[curNode].level < nodes[adjNode].level @@ -172,12 +182,16 @@ nodes(curNode).area += nodes(adjNode).area; nodes(curNode).highest += nodes(adjNode).highest; } + // curTree = Link_tree(adjTree, curTree); + // subtreeRoot(curTree) = curNode; + } curTree = Link_tree(adjTree, curTree); - lowestNode(curTree) = curNode; + subtreeRoot(curTree) = curNode; } + isproc(p) = true; } - Root = lowestNode(Find_tree(Find_node(psite(0, 0)))); + Root = subtreeRoot(Find_tree(Find_node(psite(0, 0)))); // Pour garder une map de correspondance point <-> local_root // for (int ip = 0; ip < int(S.size()); ++ip) // { @@ -185,9 +199,10 @@ // M(p) = Find_node(p); // } - mln_piter(I) r(par.domain()); + mln_piter(I) r(Par_node.domain()); for_all(r) - par(r) = Find_node(r); + Par_node(r) = Find_node(r); + } @@ -228,7 +243,7 @@ else if (Rnk_node(x) == Rnk_node(y)) Rnk_node(y) += 1; - par(x) = y; // was Par_node(x) = y; + Par_node(x) = y; // was Par_node(x) = y; return y; } @@ -260,7 +275,7 @@ void print_tree(psite p) { node& n = nodes(p); - std::cout << "psite " << p << "=" << (p.row() * exact(lowestNode).domain().len(1) + p.col()) << " : "; + std::cout << "psite " <::fwd_piter it(n.children); for (it.start(); @@ -268,7 +283,7 @@ it.next()) { psite q = it.to_psite(); - std::cout << q << "=" << (q.row() * lowestNode.domain().len(2) + q.col()) << " "; + std::cout << q << "=" << (q.row() * subtreeRoot.domain().len(1) + q.col()) << " "; } std::cout << std::endl; @@ -280,6 +295,20 @@ } } + void print_sup_tree() + { + psite max = psite(0, 0); + + for(unsigned int i = 0; i < nodes.domain().len(0); ++i) + for(unsigned int j = 0; j < nodes.domain().len(1); ++j) + { + if (nodes(psite(i, j)).area > nodes(max).area) + max = psite(i, j); + } + + print_tree(max); + } + }; // struct basic_najman }; // namespace morpho Index: sandbox/abraham/morpho/test.cc --- sandbox/abraham/morpho/test.cc (revision 1834) +++ sandbox/abraham/morpho/test.cc (working copy) @@ -25,16 +25,24 @@ image2d<int_u8> input; io::pgm::load(input, "./images/test_component_tree.pgm"); + + for(unsigned int i = 0; i < input.domain().len(0); ++i) + for(unsigned int j = 0; j < input.domain().len(1); ++j) + input.at(i, j) = 255 - input.at(i, j); + + morpho::basic_najman< image2d<int_u8>, neighb2d> n(input, c8()); n.image_output(input); n.go(); n.image_output(n.Par_tree); n.print_tree(n.Root); - n.image_output(n.par); + n.image_output(n.Par_node); std::cout << " --------------------------- " << std::endl; n.image_output(n.Rnk_tree); std::cout << " --------------------------- " << std::endl; n.image_output(n.Rnk_node); + std::cout << " --------------------------- " << std::endl; + n.print_sup_tree(); // image2d<int_u8> output = morpho::topo_wst(input, c4()); // io::pgm::save(output, "out.pgm");

1 0

1834: Sandbox: ICP: Add icp_subsampled.
by Ugo Jardonnet 03 Apr '08

03 Apr '08

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Ugo Jardonnet <ugo.jardonnet(a)lrde.epita.fr> Sandbox: ICP: Add icp_subsampled. * sandbox/jardonnet/test/bench: New info file. * sandbox/jardonnet/test/icp_check.sh: Update. * sandbox/jardonnet/test/Makefile: Add rules icp_sub. * sandbox/jardonnet/test/icp_subsampled.cc: New icp_sub. * sandbox/jardonnet/registration/icp_subsampled.hh: New. Registers, at first, a subset of the cloud. * sandbox/jardonnet/TODO: Update. * sandbox/jardonnet/registration/quat7.hh: Fix vec initialization. * sandbox/jardonnet/registration/cloud.hh: Fix vec initialization. svn wrapper unexpectedly failed to send the news. => Index: ChangeLog =================================================================== --- ChangeLog (revision 1833) +++ ChangeLog (working copy) @@ -1,3 +1,19 @@ +2008-04-03 Ugo Jardonnet <jardonnet(a)lrde.epita.fr> + + Sandbox: ICP: Add icp_subsampled. + + * sandbox/jardonnet/test/bench: New info file. + + * sandbox/jardonnet/test/icp_check.sh: Update. + * sandbox/jardonnet/test/Makefile: Add rules icp_sub. + * sandbox/jardonnet/test/icp_subsampled.cc: New icp_sub. + * sandbox/jardonnet/registration/icp_subsampled.hh: New. Registers, at + first, a subset of the cloud. + + * sandbox/jardonnet/TODO: Update. + * sandbox/jardonnet/registration/quat7.hh: Fix vec initialisation. + * sandbox/jardonnet/registration/cloud.hh: Fix vec initialisation. + 2008-04-02 Ugo Jardonnet <ugo.jardonnet(a)lrde.epita.fr> Sandbox: ICP: Fix quaternion computation. Index: sandbox/jardonnet/test/bench =================================================================== --- sandbox/jardonnet/test/bench (revision 0) +++ sandbox/jardonnet/test/bench (revision 1834) @@ -0,0 +1,32 @@ + +icp_0 +icp_lazy_0 + 0m27.662s 0m8.157s + + + +icp_0D +icp_lazy_0D + 0m7.496s 0m4.960s + + + +icp_0Df +icp_lazy_0Df + 0m8.973s 0m5.568s + + + +icp_0f +icp_lazy_0f + 0m28.490s 0m8.273s + + + +icp_3 +icp_lazy_3 + 0m5.800s 0m1.256s + + + +icp_3D +icp_lazy_3D + 0m0.752s 0m0.316s + + + +icp_3df +icp_lazy_3df + 0m0.780s 0m0.532s + + + +icp_3f +icp_lazy_3f + 0m6.468s 0m1.212s + + Index: sandbox/jardonnet/test/icp_subsampled.cc =================================================================== --- sandbox/jardonnet/test/icp_subsampled.cc (revision 0) +++ sandbox/jardonnet/test/icp_subsampled.cc (revision 1834) @@ -0,0 +1,20 @@ +#include <mln/core/image3d.hh> + +#include <mln/io/pbm/load.hh> +#include <mln/io/pbm/save.hh> + +#include <sandbox/jardonnet/registration/icp_subsampled.hh> + +int main(int, char* argv[]) +{ + using namespace mln; + + image2d< bool > img1; + image2d< bool > img2; + + io::pbm::load(img1, argv[1]); + io::pbm::load(img2, argv[2]); + + io::pbm::save(registration::icp(img1,img2), "./+registred.pbm"); +} + Index: sandbox/jardonnet/test/icp_check.sh =================================================================== --- sandbox/jardonnet/test/icp_check.sh (revision 1833) +++ sandbox/jardonnet/test/icp_check.sh (working copy) @@ -3,6 +3,6 @@ for i in `\ls bin/` do echo execute $i 01.pbm 02.pbm - ./bin/$i 01.pbm 02.pbm > ./bin/log_$i + time ./bin/$i 01.pbm 02.pbm > ./bin/log_$i echo ./bin/log_$i done Index: sandbox/jardonnet/test/Makefile =================================================================== --- sandbox/jardonnet/test/Makefile (revision 1833) +++ sandbox/jardonnet/test/Makefile (working copy) @@ -13,18 +13,23 @@ gau: g++ gaussian.cc $(FLAG) -o '+gau.exe' +run: + time ./+sub.exe . . ; time ./+gsub.exe . . + + + icp: - g++ icp.cc -I../../.. -O1 -DNDEBUG -g -o '+icp.exe' + g++ icp.cc -I../../.. -g -o '+icp_map.exe' icp++: - g++ icp.cc -I../../.. -O3 -DNDEBUG -o '+icp.exe' + g++ icp.cc -I../../.. -O3 -DNDEBUG -o '+icp_map.exe' +icp_sub: + g++ icp_subsampled.cc -I../../.. -O3 -DNDEBUG -g -o '+icp_sub.exe' + icpsafe: g++ icp.cc -fsignaling-nans -ffloat-store -I../../.. -O1 -o '+icp.exe' -run: - time ./+sub.exe . . ; time ./+gsub.exe . . - plot: icp ./+icp.exe 01.pbm 02.pbm > log gnuplot plotscript @@ -46,8 +51,9 @@ g++ icp_lazy.cc -I../../.. -O3 -ffloat-store -o './bin/+icp_lazy_3f' g++ icp_lazy.cc -I../../.. -O3 -DNDEBUG -o './bin/+icp_lazy_3D' g++ icp_lazy.cc -I../../.. -O3 -DNDEBUG -ffloat-store -o './bin/+icp_lazy_3df' + ./icp_check.sh clean: - rm -- ./bin/* + rm -f -- ./bin/* .PHONY: check \ No newline at end of file Index: sandbox/jardonnet/TODO =================================================================== --- sandbox/jardonnet/TODO (revision 1833) +++ sandbox/jardonnet/TODO (working copy) @@ -5,18 +5,9 @@ gaussian.cc: In function 'int main(int, char*)': gaussian.cc:22: error: no match for 'operator==' in 'img == out' - - +Check gaussian -adapt test for threshold (old thresholding) - - -- Enlever static dans projection::memo - -Note: -+01.bmp: 943 pts -30 itérations : 30*943 = 28290 calcul de ppp -Domaine +02.bmp: 777x480 = 362082 -nb de ppp calculé memo = 11749 - -map 1:40 -de_base 10.5 -memo 9.5 \ No newline at end of file +adapt test for threshold (old thresholding) +- - Index: sandbox/jardonnet/registration/icp_subsampled.hh =================================================================== --- sandbox/jardonnet/registration/icp_subsampled.hh (revision 0) +++ sandbox/jardonnet/registration/icp_subsampled.hh (revision 1834) @@ -0,0 +1,188 @@ +// Copyright (C) 2008 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_REGISTRATION_ICP_HH +# define MLN_REGISTRATION_ICP_HH + +/*! \file mln/registration/icp.hh + * + * \brief image registration + */ + +# include <mln/algebra/quat.hh> +# include <mln/algebra/vec.hh> +# include <mln/make/w_window.hh> +# include <mln/make/w_window3d.hh> + +# include "tools.hh" + +# include "cloud.hh" +# include "quat7.hh" +# include "projection.hh" +# include "chamfer.hh" + +namespace mln +{ + + namespace registration + { + + /*! Registration FIXME : doxy + * + * + */ + template <typename I, typename J> + inline + mln_concrete(I) + icp(const Image& cloud, + const Image<J>& surface); + +# ifndef MLN_INCLUDE_ONLY + + namespace impl + { + + template <typename P, typename M> + inline + p_array + icp_(p_array& C, + const p_array&, + M& map, + quat7<P::dim>& qk) + { + trace::entering("registration::impl::icp_"); + + quat7<P::dim> old_qk; + float err; + //float err_bis; + + p_array Ck(C), Xk(C); //FIXME: Xk sould not copy C + + algebra::vec<P::dim,float> mu_C = center(C), mu_Xk; + + const float epsilon = 1; //e-3; + + qk.apply_on(C, Ck); + + unsigned int k = 0; + do { + + //project Ck over Xk + projection::fill_Xk(Ck, map, Xk); + + mu_Xk = center(Xk); + + //compute qk + old_qk = qk; + qk = match(C, mu_C, Xk, mu_Xk); + + //Ck+1 = qk(C) + qk.apply_on(C, Ck); + + //err = d(Ck+1,Xk) + err = rms(Ck, Xk); + + //generate plot file + std::cout << k << ' ' << err << ' ' << (qk - old_qk).sqr_norm() << std::endl; + + ++k; + } while (k < 3 || (qk - old_qk).sqr_norm() > epsilon); + + trace::exiting("registration::impl::icp_"); + return Ck; + } + + } // end of namespace mln::registration::impl + + + //Only for 2d and 3d image + template <typename I, typename J> + inline + mln_concrete(I) //FIXME: should return something else ? qk ? + icp(const Image& cloud_, + const Image<J>& surface_) + { + trace::entering("registration::icp"); + mln_precondition(exact(cloud_).has_data()); + mln_precondition(exact(surface_).has_data()); + + //convert to image: time consuming + typedef image3d<mln_value(I)> I3d; + I3d cloud = convert::to_image_3d(exact(cloud_)); + const I3d surface = convert::to_image_3d(exact(surface_)); + + + //create a pair (distance map, closest point) + float w[27] = {1.4142, 1, 1.4142, 1.4142, 1, 1.4142, 0, 0, 0, + 1, 1, 1, 1, 1, 0, 0, 0, 0, + 1.4142, 1, 1.4142, 0, 0, 0, 0, 0, 0}; + w_window<mln_dpoint(I3d), float> chamfer = make::w_window3d(w); + std::pair<mln_ch_value(I3d,float), mln_ch_value(I3d,mln_point(I3d)) > map = + dt::chamfer(surface, chamfer); + + + //build p_arrays. + p_array<mln_point(I3d)> c = convert::to_p_array(cloud); + p_array<mln_point(I3d)> x = convert::to_p_array(surface); + + //init qk + quat7<3> qk; + + //create subsampled p_array + p_array<mln_point(I3d)> cprime; + for (size_t i = 0; i < c.npoints(); i += 9) //FIXME: randomize + cprime.append(c[i]); + + //apply icp 1 + impl::icp_(cprime, x, map, qk); + + //apply icp 2 + p_array<mln_point(I3d)> res = impl::icp_(c, x, map, qk); + + //to 2d : projection (FIXME:if 3d) + //mln_concrete(I) output = convert::to_image(res)? + mln_concrete(I) output(exact(cloud_).domain()); + for (size_t i = 0; i < res.npoints(); i++) + { + point2d p(res[i][0], res[i][1]); + //FIXME: not necessary if output(res.bbox()) + if (output.has(p)) + output(p) = true; + } + + trace::exiting("registration::icp"); + return output; + } + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::registration + +} // end of namespace mln + + +#endif // ! MLN_REGISTRATION_ICP_HH Index: sandbox/jardonnet/registration/quat7.hh =================================================================== --- sandbox/jardonnet/registration/quat7.hh (revision 1833) +++ sandbox/jardonnet/registration/quat7.hh (working copy) @@ -29,6 +29,7 @@ algebra::vec<n,float> _qT; quat7() + : _qR(1,0,0,0), _qT(literal::zero) { } @@ -59,6 +60,9 @@ void apply_on(const p_array& input, p_array& output) const { assert(input.npoints() == output.npoints()); + + std::cout << _qR << std::endl; + assert(_qR.is_unit()); for (size_t i = 0; i < input.npoints(); i++) Index: sandbox/jardonnet/registration/cloud.hh =================================================================== --- sandbox/jardonnet/registration/cloud.hh (revision 1833) +++ sandbox/jardonnet/registration/cloud.hh (working copy) @@ -20,7 +20,7 @@ template <typename P> P center(const p_array& a) { - algebra::vec<P::dim,float> c; + algebra::vec<P::dim,float> c(literal::zero); for (size_t i = 0; i < a.npoints(); ++i) { algebra::vec<P::dim,float> ai = a[i]; Index: sandbox/jardonnet/registration/jacobi.hh =================================================================== --- sandbox/jardonnet/registration/jacobi.hh (revision 1833) +++ sandbox/jardonnet/registration/jacobi.hh (working copy) @@ -20,11 +20,12 @@ void jacobi(algebra::mat<4,4,float> a, algebra::quat& q) { float dd, d[4]; - algebra::mat<4,4,float> v; + algebra::mat<4,4,float> v(literal::zero); int j,iq,ip,i = 0; float tresh,theta,tau,t,sm,s,h,g,c,b[4],z[4]; for (ip=0;ip<4;ip++) { - for (iq=0;iq<4;iq++) v(ip,iq)=0.0; + for (iq=0;iq<4;iq++) + v(ip,iq)=0.0; v(ip,ip)=1.0; } for (ip=0;ip<4;ip++) { Index: sandbox/jardonnet/registration/rotation.hh =================================================================== --- sandbox/jardonnet/registration/rotation.hh (revision 1833) +++ sandbox/jardonnet/registration/rotation.hh (working copy) @@ -23,7 +23,6 @@ template <unsigned n> algebra::vec<n,float> rotate(const algebra::quat& q, const algebra::vec<n,float>& p) { - return (q * algebra::quat(0. ,p) * q.inv()).v(); }

1 0

1834: Sandbox: ICP: Add icp_subsampled.
by Ugo Jardonnet 03 Apr '08

03 Apr '08

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Ugo Jardonnet <ugo.jardonnet(a)lrde.epita.fr> Sandbox: ICP: Add icp_subsampled. * sandbox/jardonnet/test/bench: New info file. * sandbox/jardonnet/test/icp_check.sh: Update. * sandbox/jardonnet/test/Makefile: Add rules icp_sub. * sandbox/jardonnet/test/icp_subsampled.cc: New icp_sub. * sandbox/jardonnet/registration/icp_subsampled.hh: New. Registers, at first, a subset of the cloud. * sandbox/jardonnet/TODO: Update. * sandbox/jardonnet/registration/quat7.hh: Fix vec initialization. * sandbox/jardonnet/registration/cloud.hh: Fix vec initialization. svn wrapper unexpectedly failed to send the the news. => Index: ChangeLog =================================================================== --- ChangeLog (revision 1833) +++ ChangeLog (working copy) @@ -1,3 +1,19 @@ +2008-04-03 Ugo Jardonnet <jardonnet(a)lrde.epita.fr> + + Sandbox: ICP: Add icp_subsampled. + + * sandbox/jardonnet/test/bench: New info file. + + * sandbox/jardonnet/test/icp_check.sh: Update. + * sandbox/jardonnet/test/Makefile: Add rules icp_sub. + * sandbox/jardonnet/test/icp_subsampled.cc: New icp_sub. + * sandbox/jardonnet/registration/icp_subsampled.hh: New. Registers, at + first, a subset of the cloud. + + * sandbox/jardonnet/TODO: Update. + * sandbox/jardonnet/registration/quat7.hh: Fix vec initialisation. + * sandbox/jardonnet/registration/cloud.hh: Fix vec initialisation. + 2008-04-02 Ugo Jardonnet <ugo.jardonnet(a)lrde.epita.fr> Sandbox: ICP: Fix quaternion computation. Index: sandbox/jardonnet/test/bench =================================================================== --- sandbox/jardonnet/test/bench (revision 0) +++ sandbox/jardonnet/test/bench (revision 1834) @@ -0,0 +1,32 @@ + +icp_0 +icp_lazy_0 + 0m27.662s 0m8.157s + + + +icp_0D +icp_lazy_0D + 0m7.496s 0m4.960s + + + +icp_0Df +icp_lazy_0Df + 0m8.973s 0m5.568s + + + +icp_0f +icp_lazy_0f + 0m28.490s 0m8.273s + + + +icp_3 +icp_lazy_3 + 0m5.800s 0m1.256s + + + +icp_3D +icp_lazy_3D + 0m0.752s 0m0.316s + + + +icp_3df +icp_lazy_3df + 0m0.780s 0m0.532s + + + +icp_3f +icp_lazy_3f + 0m6.468s 0m1.212s + + Index: sandbox/jardonnet/test/icp_subsampled.cc =================================================================== --- sandbox/jardonnet/test/icp_subsampled.cc (revision 0) +++ sandbox/jardonnet/test/icp_subsampled.cc (revision 1834) @@ -0,0 +1,20 @@ +#include <mln/core/image3d.hh> + +#include <mln/io/pbm/load.hh> +#include <mln/io/pbm/save.hh> + +#include <sandbox/jardonnet/registration/icp_subsampled.hh> + +int main(int, char* argv[]) +{ + using namespace mln; + + image2d< bool > img1; + image2d< bool > img2; + + io::pbm::load(img1, argv[1]); + io::pbm::load(img2, argv[2]); + + io::pbm::save(registration::icp(img1,img2), "./+registred.pbm"); +} + Index: sandbox/jardonnet/test/icp_check.sh =================================================================== --- sandbox/jardonnet/test/icp_check.sh (revision 1833) +++ sandbox/jardonnet/test/icp_check.sh (working copy) @@ -3,6 +3,6 @@ for i in `\ls bin/` do echo execute $i 01.pbm 02.pbm - ./bin/$i 01.pbm 02.pbm > ./bin/log_$i + time ./bin/$i 01.pbm 02.pbm > ./bin/log_$i echo ./bin/log_$i done Index: sandbox/jardonnet/test/Makefile =================================================================== --- sandbox/jardonnet/test/Makefile (revision 1833) +++ sandbox/jardonnet/test/Makefile (working copy) @@ -13,18 +13,23 @@ gau: g++ gaussian.cc $(FLAG) -o '+gau.exe' +run: + time ./+sub.exe . . ; time ./+gsub.exe . . + + + icp: - g++ icp.cc -I../../.. -O1 -DNDEBUG -g -o '+icp.exe' + g++ icp.cc -I../../.. -g -o '+icp_map.exe' icp++: - g++ icp.cc -I../../.. -O3 -DNDEBUG -o '+icp.exe' + g++ icp.cc -I../../.. -O3 -DNDEBUG -o '+icp_map.exe' +icp_sub: + g++ icp_subsampled.cc -I../../.. -O3 -DNDEBUG -g -o '+icp_sub.exe' + icpsafe: g++ icp.cc -fsignaling-nans -ffloat-store -I../../.. -O1 -o '+icp.exe' -run: - time ./+sub.exe . . ; time ./+gsub.exe . . - plot: icp ./+icp.exe 01.pbm 02.pbm > log gnuplot plotscript @@ -46,8 +51,9 @@ g++ icp_lazy.cc -I../../.. -O3 -ffloat-store -o './bin/+icp_lazy_3f' g++ icp_lazy.cc -I../../.. -O3 -DNDEBUG -o './bin/+icp_lazy_3D' g++ icp_lazy.cc -I../../.. -O3 -DNDEBUG -ffloat-store -o './bin/+icp_lazy_3df' + ./icp_check.sh clean: - rm -- ./bin/* + rm -f -- ./bin/* .PHONY: check \ No newline at end of file Index: sandbox/jardonnet/TODO =================================================================== --- sandbox/jardonnet/TODO (revision 1833) +++ sandbox/jardonnet/TODO (working copy) @@ -5,18 +5,9 @@ gaussian.cc: In function 'int main(int, char*)': gaussian.cc:22: error: no match for 'operator==' in 'img == out' - - +Check gaussian -adapt test for threshold (old thresholding) - - -- Enlever static dans projection::memo - -Note: -+01.bmp: 943 pts -30 itérations : 30*943 = 28290 calcul de ppp -Domaine +02.bmp: 777x480 = 362082 -nb de ppp calculé memo = 11749 - -map 1:40 -de_base 10.5 -memo 9.5 \ No newline at end of file +adapt test for threshold (old thresholding) +- - Index: sandbox/jardonnet/registration/icp_subsampled.hh =================================================================== --- sandbox/jardonnet/registration/icp_subsampled.hh (revision 0) +++ sandbox/jardonnet/registration/icp_subsampled.hh (revision 1834) @@ -0,0 +1,188 @@ +// Copyright (C) 2008 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_REGISTRATION_ICP_HH +# define MLN_REGISTRATION_ICP_HH + +/*! \file mln/registration/icp.hh + * + * \brief image registration + */ + +# include <mln/algebra/quat.hh> +# include <mln/algebra/vec.hh> +# include <mln/make/w_window.hh> +# include <mln/make/w_window3d.hh> + +# include "tools.hh" + +# include "cloud.hh" +# include "quat7.hh" +# include "projection.hh" +# include "chamfer.hh" + +namespace mln +{ + + namespace registration + { + + /*! Registration FIXME : doxy + * + * + */ + template <typename I, typename J> + inline + mln_concrete(I) + icp(const Image& cloud, + const Image<J>& surface); + +# ifndef MLN_INCLUDE_ONLY + + namespace impl + { + + template <typename P, typename M> + inline + p_array + icp_(p_array& C, + const p_array&, + M& map, + quat7<P::dim>& qk) + { + trace::entering("registration::impl::icp_"); + + quat7<P::dim> old_qk; + float err; + //float err_bis; + + p_array Ck(C), Xk(C); //FIXME: Xk sould not copy C + + algebra::vec<P::dim,float> mu_C = center(C), mu_Xk; + + const float epsilon = 1; //e-3; + + qk.apply_on(C, Ck); + + unsigned int k = 0; + do { + + //project Ck over Xk + projection::fill_Xk(Ck, map, Xk); + + mu_Xk = center(Xk); + + //compute qk + old_qk = qk; + qk = match(C, mu_C, Xk, mu_Xk); + + //Ck+1 = qk(C) + qk.apply_on(C, Ck); + + //err = d(Ck+1,Xk) + err = rms(Ck, Xk); + + //generate plot file + std::cout << k << ' ' << err << ' ' << (qk - old_qk).sqr_norm() << std::endl; + + ++k; + } while (k < 3 || (qk - old_qk).sqr_norm() > epsilon); + + trace::exiting("registration::impl::icp_"); + return Ck; + } + + } // end of namespace mln::registration::impl + + + //Only for 2d and 3d image + template <typename I, typename J> + inline + mln_concrete(I) //FIXME: should return something else ? qk ? + icp(const Image& cloud_, + const Image<J>& surface_) + { + trace::entering("registration::icp"); + mln_precondition(exact(cloud_).has_data()); + mln_precondition(exact(surface_).has_data()); + + //convert to image: time consuming + typedef image3d<mln_value(I)> I3d; + I3d cloud = convert::to_image_3d(exact(cloud_)); + const I3d surface = convert::to_image_3d(exact(surface_)); + + + //create a pair (distance map, closest point) + float w[27] = {1.4142, 1, 1.4142, 1.4142, 1, 1.4142, 0, 0, 0, + 1, 1, 1, 1, 1, 0, 0, 0, 0, + 1.4142, 1, 1.4142, 0, 0, 0, 0, 0, 0}; + w_window<mln_dpoint(I3d), float> chamfer = make::w_window3d(w); + std::pair<mln_ch_value(I3d,float), mln_ch_value(I3d,mln_point(I3d)) > map = + dt::chamfer(surface, chamfer); + + + //build p_arrays. + p_array<mln_point(I3d)> c = convert::to_p_array(cloud); + p_array<mln_point(I3d)> x = convert::to_p_array(surface); + + //init qk + quat7<3> qk; + + //create subsampled p_array + p_array<mln_point(I3d)> cprime; + for (size_t i = 0; i < c.npoints(); i += 9) //FIXME: randomize + cprime.append(c[i]); + + //apply icp 1 + impl::icp_(cprime, x, map, qk); + + //apply icp 2 + p_array<mln_point(I3d)> res = impl::icp_(c, x, map, qk); + + //to 2d : projection (FIXME:if 3d) + //mln_concrete(I) output = convert::to_image(res)? + mln_concrete(I) output(exact(cloud_).domain()); + for (size_t i = 0; i < res.npoints(); i++) + { + point2d p(res[i][0], res[i][1]); + //FIXME: not necessary if output(res.bbox()) + if (output.has(p)) + output(p) = true; + } + + trace::exiting("registration::icp"); + return output; + } + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::registration + +} // end of namespace mln + + +#endif // ! MLN_REGISTRATION_ICP_HH Index: sandbox/jardonnet/registration/quat7.hh =================================================================== --- sandbox/jardonnet/registration/quat7.hh (revision 1833) +++ sandbox/jardonnet/registration/quat7.hh (working copy) @@ -29,6 +29,7 @@ algebra::vec<n,float> _qT; quat7() + : _qR(1,0,0,0), _qT(literal::zero) { } @@ -59,6 +60,9 @@ void apply_on(const p_array& input, p_array& output) const { assert(input.npoints() == output.npoints()); + + std::cout << _qR << std::endl; + assert(_qR.is_unit()); for (size_t i = 0; i < input.npoints(); i++) Index: sandbox/jardonnet/registration/cloud.hh =================================================================== --- sandbox/jardonnet/registration/cloud.hh (revision 1833) +++ sandbox/jardonnet/registration/cloud.hh (working copy) @@ -20,7 +20,7 @@ template <typename P> P center(const p_array& a) { - algebra::vec<P::dim,float> c; + algebra::vec<P::dim,float> c(literal::zero); for (size_t i = 0; i < a.npoints(); ++i) { algebra::vec<P::dim,float> ai = a[i]; Index: sandbox/jardonnet/registration/jacobi.hh =================================================================== --- sandbox/jardonnet/registration/jacobi.hh (revision 1833) +++ sandbox/jardonnet/registration/jacobi.hh (working copy) @@ -20,11 +20,12 @@ void jacobi(algebra::mat<4,4,float> a, algebra::quat& q) { float dd, d[4]; - algebra::mat<4,4,float> v; + algebra::mat<4,4,float> v(literal::zero); int j,iq,ip,i = 0; float tresh,theta,tau,t,sm,s,h,g,c,b[4],z[4]; for (ip=0;ip<4;ip++) { - for (iq=0;iq<4;iq++) v(ip,iq)=0.0; + for (iq=0;iq<4;iq++) + v(ip,iq)=0.0; v(ip,ip)=1.0; } for (ip=0;ip<4;ip++) { Index: sandbox/jardonnet/registration/rotation.hh =================================================================== --- sandbox/jardonnet/registration/rotation.hh (revision 1833) +++ sandbox/jardonnet/registration/rotation.hh (working copy) @@ -23,7 +23,6 @@ template <unsigned n> algebra::vec<n,float> rotate(const algebra::quat& q, const algebra::vec<n,float>& p) { - return (q * algebra::quat(0. ,p) * q.inv()).v(); }

1 0

1833: Sandbox: ICP: Fix quaternion computation.
by Ugo Jardonnet 02 Apr '08

02 Apr '08

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Ugo Jardonnet <ugo.jardonnet(a)lrde.epita.fr> Sandbox: ICP: Fix quaternion computation. * sandbox/jardonnet/test/Makefile: Add rules for ICP. * sandbox/jardonnet/registration/quat7.hh (match): Fix matrix instantiation. * sandbox/jardonnet/registration/cloud.hh (rms): Make it use norm::l2. * sandbox/jardonnet/registration/jacobi.hh: Add safer checks, on float values. registration/cloud.hh | 14 ++++++++++++-- registration/icp.hh | 4 ---- registration/jacobi.hh | 13 +++++++++---- registration/quat7.hh | 8 +++++--- test/Makefile | 5 ++++- 5 files changed, 30 insertions(+), 14 deletions(-) Index: sandbox/jardonnet/test/Makefile --- sandbox/jardonnet/test/Makefile (revision 1832) +++ sandbox/jardonnet/test/Makefile (working copy) @@ -14,7 +14,10 @@ g++ gaussian.cc $(FLAG) -o '+gau.exe' icp: - g++ icp.cc -ffloat-store -I../../.. -O3 -DNDEBUG -o '+icp.exe' + g++ icp.cc -I../../.. -O1 -DNDEBUG -g -o '+icp.exe' + +icp++: + g++ icp.cc -I../../.. -O3 -DNDEBUG -o '+icp.exe' icpsafe: g++ icp.cc -fsignaling-nans -ffloat-store -I../../.. -O1 -o '+icp.exe' Index: sandbox/jardonnet/registration/quat7.hh --- sandbox/jardonnet/registration/quat7.hh (revision 1832) +++ sandbox/jardonnet/registration/quat7.hh (working copy) @@ -98,7 +98,7 @@ // qR - algebra::mat<P::dim,P::dim,float> Mk; + algebra::mat<P::dim,P::dim,float> Mk(literal::zero); for (size_t i = 0; i < C.npoints(); ++i) { algebra::vec<P::dim,float> Ci = C[i]; @@ -113,14 +113,16 @@ const algebra::mat<3,1,float> D = make::mat<3,1,3,float>(v); - algebra::mat<4,4,float> Qk; + algebra::mat<4,4,float> Qk(literal::zero); Qk(0,0) = tr(Mk); put(trans(D), 0,1, Qk); put(D, 1,0, Qk); put(Mk + trans(Mk) - algebra::mat<P::dim,P::dim,float>::identity() * tr(Mk), 1,1, Qk); - algebra::quat qR; + algebra::quat qR(literal::zero); + + //std::cout << qR << std::endl; jacobi(Qk, qR); // qT Index: sandbox/jardonnet/registration/cloud.hh --- sandbox/jardonnet/registration/cloud.hh (revision 1832) +++ sandbox/jardonnet/registration/cloud.hh (working copy) @@ -9,6 +9,7 @@ # include <mln/algebra/vec.hh> # include <mln/core/p_array.hh> +# include <mln/norm/l2.hh> namespace mln { @@ -32,9 +33,9 @@ // FIXME : move //exist for P? template <typename P> - mln_coord(P) sqr_norm(const P& v) + float sqr_norm(const P& v) { - mln_coord(P) tmp = 0; + float tmp = 0; for (unsigned i = 0; i < P::dim; i++) tmp += v[i] * v[i]; return tmp; @@ -45,9 +46,18 @@ const p_array& a2) { assert(a1.npoints() == a2.npoints()); + /* float f = 0.f; for (size_t i = 0; i < a1.npoints(); ++i) f += sqr_norm(a1[i] - a2[i]); + return f / a1.npoints();*/ + float f = 0.f; + for (size_t i = 0; i < a1.npoints(); ++i) + { + algebra::vec<3,float> a1f = a1[i]; + algebra::vec<3,float> a2f = a2[i]; + f += norm::l2_distance(a1f,a2f); + } return f / a1.npoints(); } Index: sandbox/jardonnet/registration/jacobi.hh --- sandbox/jardonnet/registration/jacobi.hh (revision 1832) +++ sandbox/jardonnet/registration/jacobi.hh (working copy) @@ -4,6 +4,8 @@ #include <mln/algebra/mat.hh> +#include <cmath> +#include "misc.hh" // from num. rec. in C @@ -35,7 +37,7 @@ for (iq=ip+1;iq<4;iq++) sm += fabs(a(ip,iq)); } - if (sm < 1e-12) { + if (sm < 1e-6) { //1e-12 dd = d[0]; iq = 0; for (ip=1;ip<4;ip++) @@ -61,14 +63,17 @@ /* unusefull */ g=100.0*fabs(a(ip,iq)); - if (i > 4 && (float)(fabs(d[ip])+g) == (float)fabs(d[ip]) - && (float)(fabs(d[iq])+g) == (float)fabs(d[iq])) + //if (i > 4 && (float)(fabs(d[ip])+g) == (float)fabs(d[ip]) + // && (float)(fabs(d[iq])+g) == (float)fabs(d[iq])) + if (i > 4 && about_equal((float)(fabs(d[ip])+g), (float)fabs(d[ip])) + && about_equal((float)(fabs(d[iq])+g), (float)fabs(d[iq]))) a(ip,iq)=0.0; /* unusefull */ else if (fabs(a(ip,iq)) > tresh) { h=d[iq]-d[ip]; - if ((float)(fabs(h)+g) == (float)fabs(h)) + //if ((float)(fabs(h)+g) == (float)fabs(h)) // unsafe ? + if (about_equal((float)(fabs(h)+g), (float)fabs(h))) t=(a(ip,iq))/h; else { theta=0.5*h/(a(ip,iq)); Index: sandbox/jardonnet/registration/icp.hh --- sandbox/jardonnet/registration/icp.hh (revision 1832) +++ sandbox/jardonnet/registration/icp.hh (working copy) @@ -89,7 +89,6 @@ k = 0; do { - //std::cout << "step 2" << std::endl; //// step 2 projection::fill_Xk(Ck, map, Xk); @@ -98,16 +97,13 @@ mu_Xk = center(Xk); - //std::cout << "step 3" << std::endl; //// step 3 old_qk = qk; qk = match(C, mu_C, Xk, mu_Xk); - //std::cout << "step 4" << std::endl; //// step 4 qk.apply_on(C, Ck); // Ck+1 = qk(C) - //std::cout << "step err" << std::endl; //// err = d(Ck+1,Xk) err = rms(Ck, Xk); std::cout << k << ' ' << err << ' ' << (qk - old_qk).sqr_norm() << std::endl; //plot file

1 0

1831: Improve ctors of mln::algebra::mat.
by Ugo Jardonnet 02 Apr '08

02 Apr '08

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Ugo Jardonnet <ugo.jardonnet(a)lrde.epita.fr> Improve ctors of mln::algebra::mat. * mln/algebra/mat.hh (mat(const literal::zero_t&)): Add ctor. mat.hh | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) Index: mln/algebra/mat.hh --- mln/algebra/mat.hh (revision 1831) +++ mln/algebra/mat.hh (working copy) @@ -95,9 +95,9 @@ static const mat<n,m,T> Id; - mat() - { - } + mat(); + + mat(const literal::zero_t&); template <typename U> mat(const mat<n,m,U>& rhs); @@ -280,6 +280,19 @@ } template <unsigned n, unsigned m, typename T> + inline + mat<n,m,T>::mat() + { + } + + template <unsigned n, unsigned m, typename T> + inline + mat<n,m,T>::mat(const literal::zero_t&) + { + this->set_all(0); + } + + template <unsigned n, unsigned m, typename T> template <typename U> inline mat<n,m,T>::mat(const mat<n,m,U>& rhs)

1 0

1831: Take milena/tools/ into account in configure.ac.
by Roland Levillain 02 Apr '08

02 Apr '08

https://svn.lrde.epita.fr/svn/oln/trunk Index: ChangeLog from Roland Levillain <roland(a)lrde.epita.fr> Take milena/tools/ into account in configure.ac. * configure.ac: Configure milena/tools/Makefile (TESTS_CXXFLAGS): Adjust these flags. (TOOLS_CXXFLAGS): New. configure.ac | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) Index: configure.ac --- configure.ac (revision 1830) +++ configure.ac (working copy) @@ -42,7 +42,14 @@ AC_ARG_VAR([TESTS_CXXFLAGS]) # We want no optimization for the tests (too slow), and a lot of debugging. if test "$GXX" = yes && test -z "$TESTS_CXXFLAGS"; then - TESTS_CXXFLAGS="-O0 -ggdb -Wall -pipe" + TESTS_CXXFLAGS="-O0 -ggdb -Wall -W" +fi + +# C++ compiler flags for tools. +AC_ARG_VAR([TOOLS_CXXFLAGS]) +# On the contrary, we want fast binaries for tools. +if test "$GXX" = yes && test -z "$TOOLS_CXXFLAGS"; then + TOOLS_CXXFLAGS="-O3 -NDEBUG -ggdb" fi @@ -71,6 +78,7 @@ Makefile milena/Makefile milena/doc/Makefile + milena/tools/Makefile ])

1 0

1830: Initiate a set of command-line tools.
by Roland Levillain 02 Apr '08

02 Apr '08

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Roland Levillain <roland(a)lrde.epita.fr> Initiate a set of command-line tools. * tools/: New directory. * tools/area_flooding.cc: New. * tools/Makefile.am, tools/tools.mk: New. * Makefile.am (SUBDIRS): Add tools. Makefile.am | 2 tools/Makefile.am | 7 + tools/area_flooding.cc | 189 +++++++++++++++++++++++++++++++++++++++++++++++++ tools/tools.mk | 15 +-- 4 files changed, 202 insertions(+), 11 deletions(-) Index: Makefile.am --- Makefile.am (revision 1829) +++ Makefile.am (working copy) @@ -1,6 +1,6 @@ ## Process this file through Automake to create Makefile.in -*- Makefile -*- -SUBDIRS = doc tests +SUBDIRS = doc tests tools .PHONY: doc doc: Index: tools/tools.mk --- tools/tools.mk (revision 0) +++ tools/tools.mk (working copy) @@ -1,18 +1,13 @@ ## FIXME: Revamp (see Olena 0.11's test suite). -# Look for Milena header in srcdir and for (generated) test headers -# (like data.hh) in builddir. -AM_CPPFLAGS = -I$(top_srcdir)/milena -I$(top_builddir)/milena +# Look for Milena header in srcdir. +AM_CPPFLAGS = -I$(top_srcdir)/milena # FIXME: Add # # AM_CXXFLAGS = $(CXXFLAGS_STRICT) $(CXXFLAGS_DEBUG) # # when oln.m4 (or just a part of it) is available in the distribution, -# instead of using the hand-made TESTS_CXXFLAGS. -TESTS_CXXFLAGS = @TESTS_CXXFLAGS@ -AM_CXXFLAGS = $(TESTS_CXXFLAGS) - -# FIXME: Likewise, we should compute these values at configure time. -# Hard-code them for the moment. -LONG_TESTS_CXXFLAGS = -O3 +# instead of using the hand-made TOOLS_CXXFLAGS. +TOOLS_CXXFLAGS = @TOOLS_CXXFLAGS@ +AM_CXXFLAGS = $(TOOLS_CXXFLAGS) Index: tools/Makefile.am --- tools/Makefile.am (revision 0) +++ tools/Makefile.am (revision 0) @@ -0,0 +1,7 @@ +## Process this file through Automake to create Makefile.in -*- Makefile -*- + +include $(top_srcdir)/milena/tools/tools.mk + +bin_PROGRAMS = area_flooding + +area_flooding_SOURCES = area_flooding.cc Index: tools/area_flooding.cc --- tools/area_flooding.cc (revision 0) +++ tools/area_flooding.cc (revision 0) @@ -0,0 +1,189 @@ +// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) +// +// 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. + +/** \file tests/morpho/lena_line_graph_image_wst2.cc + \brief More tests on the Watershed Transform (WST) on a + mln::line_graph_image. + + The scenario is as follows: + \li load a 2-D, gray-level image from a PGM file; + \li convert this 2-D image into a line graph-based one, where values + on edges are computed as the absolute value f the difference + between the values on the nodes adjacent to the edge, so as to + create a (norm of the gradient) ``between the pixels'' of the + input image; + \li perform a hierarchical flooding using an area opening, until a + maximum number of regions has been reached; + \li compute a watershed transform on this image; + \li create a 2D image of this watershed image, where each basin has + an average grey level of the corresponding region in the input + image. +*/ + +#include <cstdlib> + +#include <vector> + +#include <mln/value/int_u8.hh> +#include <mln/value/int_u16.hh> + +#include <mln/core/image2d.hh> +#include <mln/core/neighb2d.hh> + +#include <mln/core/line_graph_image.hh> +#include <mln/core/line_graph_elt_neighborhood.hh> +#include <mln/core/line_graph_neighborhood_piter.hh> + +#include <mln/morpho/line_gradient.hh> +#include <mln/morpho/closing_area_on_vertices.hh> +#include <mln/labeling/regional_minima.hh> +#include <mln/morpho/meyer_wst.hh> + +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> + +int main(int argc, char* argv[]) +{ + if (argc != 4) + { + std::cerr + << "usage: " << argv[0] << " max_nregions input.pgm output.pgm" + << std::endl; + std::exit(EXIT_FAILURE); + } + + /*--------. + | Input. | + `--------*/ + + using namespace mln; + using value::int_u8; + using value::int_u16; + + typedef int_u8 val_t; + typedef image2d<val_t> orig_ima_t; + + orig_ima_t input; + io::pgm::load(input, argv[2]); + if (!input.has_data()) + { + std::cerr << "Error reading input " << argv[2] << std::endl; + std::exit(2); + } + + /*----------------. + | Line gradient. | + `----------------*/ + + // Line graph image. + typedef line_graph_image<point2d, val_t> ima_t; + ima_t lg_ima = morpho::line_gradient(input); + + /*-----------. + | Flooding. | + `-----------*/ + + /* FIXME: I'm not sure this is the way it should be done. Anyway, + we should implemented this as a canvas. */ + + typedef line_graph_elt_neighborhood<point2d> nbh_t; + nbh_t nbh; + + unsigned area = 0; + unsigned max_area = input.npoints(); + unsigned nregions = mln_max(unsigned); + unsigned max_nregions = atoi(argv[1]); + + ima_t result = clone(lg_ima); + while (area < max_area && nregions > max_nregions) + { + ++area; + std::cerr << "area = " << area << " \t" + << "nregions = " << nregions << std::endl; + ima_t work = clone(result); + // Compute the closing. + morpho::closing_area_on_vertices(work, nbh, area, result); + // Compute the number of local minima (but get rid of the image, + // as we don't need it). + labeling::regional_minima(result, nbh, nregions); + } + + /*------. + | WST. | + `------*/ + + // Perform a Watershed Transform. + typedef int_u16 wst_val_t; + wst_val_t nbasins; + typedef line_graph_image<point2d, wst_val_t> wst_ima_t; + wst_ima_t wshed = morpho::meyer_wst(result, nbh, nbasins); + std::cout << "nbasins = " << nbasins << std::endl; + + /*---------. + | Output. | + `---------*/ + + const wst_val_t wshed_label = 0; + + // Create a 2D-equivalent of WSHED. + image2d<wst_val_t> wshed2d(input.domain()); + + /* FIXME: It'd better if we could iterate over the *vertices* of a + line graph image. We could avoid all this lengthy code. */ + // Iterate over each edge of the watershed image, and propagate the + // label of an edge to its adjacent vertices when this edge is not + // part of the watershed. + mln_piter_(wst_ima_t) p(wshed.domain()); + for_all(p) + if (wshed(p) != wshed_label) + { + wshed2d(p.to_psite().first()) = wshed(p); + wshed2d(p.to_psite().second()) = wshed(p); + } + + // For each basin, compute the average gray level. + std::vector<mln_sum_(val_t)> sum(nbasins + 1, 0); + std::vector<unsigned> nsites(nbasins + 1, 0); + mln_piter_(orig_ima_t) q(input.domain()); + for_all(q) + { + sum[wshed2d(q)] += input(q); + ++nsites[wshed2d(q)]; + } + std::vector<float> average(nbasins + 1); + for (unsigned i = 1; i <= nbasins; ++i) + average[i] = float (sum[i]) / float (nsites[i]); + + // Create an output image using the average gray levels of the basins. + orig_ima_t output(input.domain()); + for_all(q) + output(q) = average[wshed2d(q)]; + + std::cout << "area = " << area << " \t" + << "nregions = " << nregions << std::endl; + io::pgm::save(output, argv[3]); +}

1 0

1829: Add inclusion operators for mln::p_graph and mln::p_line_graph.
by Roland Levillain 02 Apr '08

02 Apr '08

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Roland Levillain <roland(a)lrde.epita.fr> Add inclusion operators for mln::p_graph and mln::p_line_graph. * mln/core/p_graph.hh (operator<=(const p_graph&, const p_graph&)) * mln/core/p_line_graph.hh (operator<=(const p_line_graph&, const p_line_graph&)): New operators. p_graph.hh | 19 +++++++++++++++++++ p_line_graph.hh | 19 +++++++++++++++++++ 2 files changed, 38 insertions(+) Index: mln/core/p_graph.hh --- mln/core/p_graph.hh (revision 1828) +++ mln/core/p_graph.hh (working copy) @@ -126,6 +126,18 @@ bool operator==(const p_graph& lhs, const p_graph& rhs); + /// \brief Inclusion of a mln::p_graph in another one. + /// + /// This inclusion relation is very strict for the moment, since our + /// infrastrure for graphs i simple: a mln::p_graph is included + /// in another one if their are equal. + /// + /// \todo Refine this later, when we are able to express subgraph + /// relations. + template <typename P> + bool + operator<=(const p_graph& lhs, const p_graph& rhs); + # ifndef MLN_INCLUDE_ONLY @@ -273,6 +285,13 @@ return lhs.gr_.ptr_ == rhs.gr_.ptr_; } + template <typename P> + bool + operator<=(const p_graph& lhs, const p_graph& rhs) + { + return lhs == rhs; + } + # endif // ! MLN_INCLUDE_ONLY } // end of mln Index: mln/core/p_line_graph.hh --- mln/core/p_line_graph.hh (revision 1828) +++ mln/core/p_line_graph.hh (working copy) @@ -105,6 +105,18 @@ bool operator==(const p_line_graph& lhs, const p_line_graph& rhs); + /// \brief Inclusion of a mln::p_line_graph in another one. + /// + /// This inclusion relation is very strict for the moment, since our + /// infrastrure for graphs i simple: a mln::p_line_graph is included + /// in another one if their are equal. + /// + /// \todo Refine this later, when we are able to express subgraph + /// relations. + template <typename P> + bool + operator<=(const p_line_graph& lhs, const p_line_graph& rhs); + # ifndef MLN_INCLUDE_ONLY @@ -173,6 +185,13 @@ return lhs.gr_.ptr_ == rhs.gr_.ptr_; } + template <typename P> + bool + operator<=(const p_line_graph& lhs, const p_line_graph& rhs) + { + return lhs == rhs; + } + # endif // ! MLN_INCLUDE_ONLY } // end of mln

1 0

1828: Improve ctors of mln::graph_psite and mln::line_graph_psite.
by Roland Levillain 02 Apr '08

02 Apr '08

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Roland Levillain <roland(a)lrde.epita.fr> Improve ctors of mln::graph_psite and mln::line_graph_psite. * mln/core/graph_psite.hh (super_): New typedef. (mln::graph_psite::graph_psite): Explicitly call the ctor of the super class. * mln/core/line_graph_psite.hh (super_): New typedef. (mln::line_graph_psite::line_graph_psite): Explicitly call the ctor of the super class. graph_psite.hh | 10 +++++++--- line_graph_psite.hh | 10 +++++++--- 2 files changed, 14 insertions(+), 6 deletions(-) Index: mln/core/graph_psite.hh --- mln/core/graph_psite.hh (revision 1827) +++ mln/core/graph_psite.hh (working copy) @@ -46,6 +46,7 @@ class graph_psite : public Point_Site< graph_psite > { typedef graph_psite self_; + typedef Point_Site<self_> super_; public: typedef mln_mesh(P) mesh; @@ -100,7 +101,8 @@ inline graph_psite::graph_psite() // Dummy initializations. - : pg_(0), + : super_(), + pg_(0), id_(-1) { } @@ -108,7 +110,8 @@ template<typename P> inline graph_psite::graph_psite(const p_graph& g, util::node_id id) - : pg_(&g), + : super_(), + pg_(&g), id_(id) { } @@ -116,7 +119,8 @@ template<typename P> inline graph_psite::graph_psite(const graph_psite& rhs) - : pg_(rhs.pg_), + : super_(), + pg_(rhs.pg_), id_(rhs.id_) { } Index: mln/core/line_graph_psite.hh --- mln/core/line_graph_psite.hh (revision 1827) +++ mln/core/line_graph_psite.hh (working copy) @@ -49,6 +49,7 @@ class line_graph_psite : public Point_Site< line_graph_psite > { typedef line_graph_psite self_; + typedef Point_Site<self_> super_; public: typedef P point; @@ -138,7 +139,8 @@ inline line_graph_psite::line_graph_psite() // Dummy initializations. - : plg_(0), + : super_(), + plg_(0), id_(-1), p_() { @@ -148,7 +150,8 @@ inline line_graph_psite::line_graph_psite(const p_line_graph& plg, util::edge_id id) - : plg_(&plg), + : super_(), + plg_(&plg), id_(id), p_() { @@ -157,7 +160,8 @@ template<typename P> inline line_graph_psite::line_graph_psite(const line_graph_psite& rhs) - : plg_(rhs.plg_), + : super_(), + plg_(rhs.plg_), id_(rhs.id_), p_() {

1 0

1827: Fix a few comments.
by Roland Levillain 02 Apr '08

02 Apr '08

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Roland Levillain <roland(a)lrde.epita.fr> * mln/accu/count_adjacent_vertices.hh, * mln/morpho/meyer_wst.hh, * tests/morpho/lena_line_graph_image_wst1.cc, * tests/morpho/lena_line_graph_image_wst2.cc: Fix comments. mln/accu/count_adjacent_vertices.hh | 3 +++ mln/morpho/meyer_wst.hh | 6 ------ tests/morpho/lena_line_graph_image_wst1.cc | 6 +++--- tests/morpho/lena_line_graph_image_wst2.cc | 7 +++---- 4 files changed, 9 insertions(+), 13 deletions(-) Index: mln/accu/count_adjacent_vertices.hh --- mln/accu/count_adjacent_vertices.hh (revision 1826) +++ mln/accu/count_adjacent_vertices.hh (working copy) @@ -49,6 +49,9 @@ /// /// The type to be count is mln::util::pix< mln::line_graph_image<P, V> > /// where \p P and \p V are the parameters of this class. + /// + /// This accumulator is used by mln::closing_area_on_vertices and + /// mln::opening_area_on_vertices. template <typename P, typename V> struct count_adjacent_vertices_ : public mln::accu::internal::base_< std::size_t, Index: mln/morpho/meyer_wst.hh --- mln/morpho/meyer_wst.hh (revision 1826) +++ mln/morpho/meyer_wst.hh (working copy) @@ -113,12 +113,6 @@ const marker unmarked = mln_min(marker); // Initialize the output with the markers (minima components). - /* FIXME: labeling::regional_minima doesn't work on non - point-accessible images! We temporarily reuse the old - morpho::minima_components routine to work-around this - limitation. As soon as labeling::regional_minima work, get rid of - - mln/morpho/level_components.hh, and - - mln/morpho/extrema_components.hh. */ mln_ch_value(I, marker) markers = labeling::regional_minima (input, nbh, nbasins); Index: tests/morpho/lena_line_graph_image_wst1.cc --- tests/morpho/lena_line_graph_image_wst1.cc (revision 1826) +++ tests/morpho/lena_line_graph_image_wst1.cc (working copy) @@ -26,11 +26,11 @@ // Public License. /* FIXME: We should factor as much things as possible between - tests/core/lena_line_graph_image_wst1.cc and - tests/core/lena_line_graph_image_wst2.cc, starting from conversion + tests/morpho/lena_line_graph_image_wst1.cc and + tests/morpho/lena_line_graph_image_wst2.cc, starting from conversion routines. */ -/** \file tests/core/lena_line_graph_image_wst1.cc +/** \file tests/morpho/lena_line_graph_image_wst1.cc \brief Tests on the Watershed Transform (WST) on a mln::line_graph_image. The scenario is as follows: Index: tests/morpho/lena_line_graph_image_wst2.cc --- tests/morpho/lena_line_graph_image_wst2.cc (revision 1826) +++ tests/morpho/lena_line_graph_image_wst2.cc (working copy) @@ -26,17 +26,16 @@ // Public License. /* FIXME: We should factor as much things as possible between - tests/core/lena_line_graph_image_wst1.cc and - tests/core/lena_line_graph_image_wst2.cc, starting from conversion + tests/morpho/lena_line_graph_image_wst1.cc and + tests/morpho/lena_line_graph_image_wst2.cc, starting from conversion routines. */ -/** \file tests/core/lena_line_graph_image_wst2.cc +/** \file tests/morpho/lena_line_graph_image_wst2.cc \brief More tests on the Watershed Transform (WST) on a mln::line_graph_image. The scenario is as follows: \li load a 2-D, gray-level image from a PGM file; - \li simplify the image to reduce the number of local minima; \li convert this 2-D image into a line graph-based one, where values on edges are computed as the absolute value f the difference between the values on the nodes adjacent to the edge, so as to

1 0