Commit 94996ead authored by Thierry Geraud's avatar Thierry Geraud
Browse files

Update Laurent's ISMM code.

	* laurent/ismm2009.cc: Update.
	* laurent/ismm2009.v1.cc: New.
	* laurent/ismm2009.hh (is_not_point): New.
	(find_root): Remove.


git-svn-id: https://svn.lrde.epita.fr/svn/oln/trunk@3163 4aad255d-cdde-0310-9447-f3009e2ae8c0
parent beddc3b8
2009-01-19 Thierry Geraud <thierry.geraud@lrde.epita.fr>
Update Laurent's ISMM code.
* laurent/ismm2009.cc: Update.
* laurent/ismm2009.v1.cc: New.
* laurent/ismm2009.hh (is_not_point): New.
(find_root): Remove.
2009-01-16 Matthieu Garrigues <matthieu.garrigues@gmail.com>
Work on reconstruction.
......
......@@ -2,16 +2,216 @@
#include <mln/value/int_u8.hh>
#include <mln/value/label_8.hh>
#include <mln/value/label_16.hh>
#include <mln/io/pgm/load.hh>
#include <mln/util/ord_pair.hh>
#include <mln/debug/println.hh>
#include <mln/core/routine/duplicate.hh>
# include <mln/core/site_set/p_set.hh>
# include <mln/core/site_set/p_queue.hh>
# include <mln/core/site_set/p_priority.hh>
#include <mln/labeling/compute.hh>
#include <mln/accu/count.hh>
/*
TO-DO list:
-----------
- having a heap for every lr (support merge)
- handling adjacency on the fly (instead of having an image)
*/
namespace mln
{
// Region label parenthood routines.
template <typename L>
inline
void make_sets(std::vector<L>& lpar, L l_max)
{
for (L l = 1; l <= l_max; ++l)
lpar[l] = l; // Make-set.
}
template <typename L>
inline
L find_root_l(std::vector<L>& lpar, L l)
{
if (lpar[l] == l)
return l;
else
return lpar[l] = find_root_l(lpar, lpar[l]);
}
// Edge tree.
template <typename I>
inline
mln_psite(I)
find_root_e(I& z_epar, mln_psite(I) e)
{
if (z_epar(e) == e)
return e;
else
return z_epar(e) = find_root_e(z_epar, z_epar(e));
}
// Get smallest edge.
template <typename Qs, typename L, typename W>
mln_deduce(Qs, element, element)
get_smallest_edge(Qs& qs, L l, const W& wst, std::vector<L>& lpar, bool& found)
{
typedef mln_element(Qs) Q; // qs is an array of queues with type Q.
{
// Test that, when a region has merged, its edge queue has been
// emptied.
L l_ = l;
while (lpar[l_] != l_)
{
mln_invariant(qs[l_].is_empty());
l_ = lpar[l_];
}
}
L lr = find_root_l(lpar, l);
Q& q = qs[lr];
typedef mln_element(Q) E;
while (! q.is_empty())
{
const E& e = q.front();
L l1 = wst(p1_from_e(e)),
l2 = wst(p2_from_e(e));
mln_invariant(l1 != l2);
L l1r = find_root_l(lpar, l1),
l2r = find_root_l(lpar, l2);
mln_invariant(l1r == lr || l2r == lr);
if (l1r == l2r)
// 'e' is an internal edge => forget it.
{
q.pop();
continue;
}
found = true;
return e;
}
found = false;
return E(0,0); // FIXME: null edge.
}
// Compute an LCA; reference(?) code.
template <typename L, typename E, typename Epar>
std::vector< std::vector<E> >
compute_lca(const L& l_max,
const std::vector<E>& edge,
const Epar& epar)
{
mln_ch_value(Epar, std::vector<L>) chl_; // Children (known as array).
{
initialize(chl_, epar);
E e;
for (L l = 1; l <= l_max; ++l)
{
e = edge[l];
do
{
chl_(e).push_back(l);
e = epar(e);
}
while (e != epar(e));
chl_(e).push_back(l);
}
// e is the root edge.
mln_invariant(chl_(e).size() == l_max); // All labels are children.
}
/*
// Display children tree.
{
std::cout << "chl_ tree: " << std::endl;
for (L l = 1; l <= l_max; ++l)
{
E e_ = edge[l];
std::cout << "l = " << l << " => ";
do
{
std::cout << e_ << " : ";
for (unsigned i = 0; i < chl_(e_).size(); ++i)
std::cout << chl_(e_)[i] << ' ';
std::cout << " -> ";
e_ = epar(e_);
}
while (epar(e_) != e_);
std::cout << e_ << " : ";
for (unsigned i = 0; i < chl_(e_).size(); ++i)
std::cout << chl_(e_)[i] << ' ';
std::cout << std::endl;
}
}
*/
mln_ch_value(Epar, std::set<L>) chl; // Children (as a set).
{
initialize(chl, epar);
mln_piter(Epar) e(chl.domain());
for_all(e)
if (chl_(e).size() != 0)
chl(e).insert(chl_(e).begin(), chl_(e).end());
}
unsigned size = l_max.next();
std::vector< std::vector<E> > lca(size);
for (L l = 1; l <= l_max; ++l)
{
lca[l].resize(size);
// Diagonal.
lca[l][l] = edge[l]; // Itself.
// Superior right.
for (L l2 = l.next(); l2 <= l_max; ++l2)
{
E e = edge[l];
while (chl(e).find(l2) == chl(e).end())
e = epar(e);
lca[l][l2] = e;
}
}
return lca;
}
} // mln
void usage(char* argv[])
{
......@@ -22,11 +222,12 @@ void usage(char* argv[])
int main(int argc, char* argv[])
{
using namespace mln;
using value::int_u8;
using value::label_8;
using value::label_16;
if (argc != 2)
usage(argv);
......@@ -34,27 +235,37 @@ int main(int argc, char* argv[])
image2d<int_u8> raw_f;
io::pgm::load(raw_f, argv[1]);
image2d<int_u8> f_ = add_edges(raw_f);
typedef image2d<int_u8> Complex;
typedef mln_pset_(Complex) full_D;
Complex f_ = add_edges(raw_f);
mln_VAR(f, f_ | is_pixel);
// debug::println("f:", f);
mln_VAR(g, f_ | is_edge);
typedef label_8 L; // Type of labels.
typedef mln_value_(g_t) T; // Type of edge values.
typedef mln_psite_(g_t) E; // Type of edges.
typedef label_16 L; // Type of labels.
L n_basins;
mln_VAR( wst_g,
f_to_wst_g(f, g, e2p(), e2e(), n_basins) );
// debug::println("g:", g);
std::cout << "n basins = " << n_basins << std::endl;
debug::println("g:", g);
debug::println("wst(g):", wst_g);
// Just to see things.
//
mln_VAR(adj, adjacency(wst_g, e2e()));
debug::println("adjacency:", adj | (pw::value(wst_g) == pw::cst(0)));
debug::println("adj:", adj);
mln_VAR(adj_line, adj | (pw::value(wst_g) == pw::cst(0)));
debug::println("adjacency:", adj_line);
......@@ -71,33 +282,89 @@ int main(int argc, char* argv[])
c4().win()),
wst_g_full);
}
debug::println("wst(g) fully valuated:", wst_g_full);
// Depict the watershed line.
mln_VAR(is_line, pw::value(wst_g_full) == pw::cst(0));
debug::println("wst(g) line:", wst_g | is_line);
// --- debug::println("wst(g) line:", wst_g | is_line);
mln_VAR(g_line, g | is_line);
debug::println("g | wst(g) line:", g_line);
// Get the smallest edge out of every basin.
std::vector<point2d> edge = smallest_edges(extend(wst_g | is_line, wst_g_full),
n_basins, e2p(), g);
for (L l = 1; l <= n_basins; ++l)
std::cout << int(l) << ": " << edge[l] << std::endl;
// Priority queue -> edge.
typedef p_priority< T, p_queue<E> > Q;
util::array<Q> q(n_basins.next());
{
mln_piter_(g_t) e(g.domain());
for_all(e)
if (wst_g(e) == 0)
{
L l1 = wst_g_full(p1_from_e(e)),
l2 = wst_g_full(p2_from_e(e));
if (l1 > l2)
std::swap(l1, l2);
mln_invariant(l1 < l2);
q[l1].push(mln_max(T) - g(e), e);
q[l2].push(mln_max(T) - g(e), e);
}
// --- for (L l = 1; l <= n_basins; ++l)
// --- std::cout << "q["<< l << "] = " << q[l] << std::endl;
}
// To know if an edge is out of a given region (label l), we
// maintain the information about region merging using
// an union-find structure.
std::vector<L> lpar(n_basins.next());
make_sets(lpar, n_basins);
// Given a region R with label l and an edge e = (l1, l2) from the
// priority queue, we get know if that edge is out of the set of
// regions containing l: we shall have l1r = lr or l2r = lr.
// Please note that an edge (l1, l2) is internal to a set of regions
// if l1r = l2r.
// // Test the 'get_smallest_edge' routine.
// {
// // Test the result with an alternate code.
// std::vector<point2d> edge_alt = smallest_edges_alt(wst_g, n_basins, e2e(), g);
// bool found;
// for (L l = 1; l <= n_basins; ++l)
// mln_invariant(edge_alt[l] == edge[l]);
// {
// E e = get_smallest_edge(q, l, wst_g_full, lpar, found);
// std::cout << l << ": e = " << e
// << " -- g(e) = " << g(e) << std::endl;
// }
// }
E null = E(0,0); // Impossible value.
// Information "label l -> edge e".
std::vector<E> edge(n_basins.next());
for (L l = 1; l <= n_basins; ++l)
edge[l] = null;
// Compute an attribute per region.
// --------------------------------
typedef unsigned A;
util::array<A> a = labeling::compute(accu::meta::count(),
......@@ -109,49 +376,385 @@ int main(int argc, char* argv[])
std::vector<pair_t> v = sort_attributes(a, n_basins); // Sort regions.
std::cout << "attributes = ";
for (unsigned i = 1; i <= n_basins; ++i)
std::cout << v[i].first // Attribute (increasing).
<< ':'
<< v[i].second // Region label.
<< " - ";
std::cout << std::endl;
// Display attributes.
{
std::cout << "(attribute, label) = { ";
for (unsigned i = 1; i <= n_basins; ++i)
std::cout << '(' << v[i].first // Attribute (increasing).
<< ','
<< v[i].second // Region label.
<< ") ";
std::cout << '}' << std::endl
<< std::endl;
}
std::vector<L> lpar(n_basins.next());
for (L l = 1; l <= n_basins; ++l)
lpar[l] = l; // Make-set.
util::array<A> a_merged = a;
for (unsigned i = 1; i <= n_basins; ++i)
// Go!
// --------------------------------
mln_ch_value_(g_line_t, E)
epar, // Edge forest.
z_epar; // Auxiliary data: edge forest with compression and balancing.
{
initialize(epar, g_line);
initialize(z_epar, g_line);
mln_piter_(g_line_t) e(g_line.domain());
for_all(e)
{
L l = v[i].second, // Region label.
lr = find_root(lpar, l);
point2d e = edge[l]; // FIXME: Use the heap!
mln_invariant(wst_g_full(p1_from_e(e)) == l ||
wst_g_full(p2_from_e(e)) == l);
L l2 = ( wst_g_full(p1_from_e(e)) == l ?
wst_g_full(p2_from_e(e)) :
wst_g_full(p1_from_e(e)) ),
l2r = find_root(lpar, l2);
if (lr == l2r)
continue; // Already merged.
if (l2r < lr)
std::swap(lr, l2r);
mln_invariant(l2r > lr);
lpar[lr] = l2r;
a_merged[l2r] += lr; // FIXME: We want accumulators here!
// Make-Set.
epar(e) = e;
z_epar(e) = e;
}
debug::println("all edges:", epar); // epar(e) == e so we depict the edges!
}
// 'aa' is the result attribute image; it is defined on the complex
// (so it has edges, pixels, and 0-face-points).
image2d<A> aa;
initialize(aa, f_);
data::fill(aa, 0); // Safety iff 0 is an invalidate attribute value!
for (unsigned i = 1; i <= n_basins; ++i)
{
L l = v[i].second;
std::cout << l << " -> " << lpar[l] << std::endl;
L l = v[i].second; // Region label.
bool found;
E e = get_smallest_edge(q, l, wst_g_full, lpar, found);
if (! found)
{
// The last iteration is useless: all regions have already
// merged. We could have written: "for i = 1 to n_basins - 1".
mln_invariant(i == n_basins);
break;
}
L l1 = wst_g_full(p1_from_e(e)),
l2 = wst_g_full(p2_from_e(e)),
l1r = find_root_l(lpar, l1),
l2r = find_root_l(lpar, l2);
aa(e) = a[l];
// Consistency tests.
{
mln_invariant(a[l] == v[i].first);
mln_invariant(epar(e) == e); // 'e' has not been processed yet.
mln_invariant(l1 != l2); // It is a valid edge, i.e., between 2 different regions...
mln_invariant(l1r != l2r); // ...that are not already merged.
L lr = find_root_l(lpar, l);
mln_invariant(lr == l // Either l is not deja-vu
|| ((lr == l1r && lr != l2r) || // or l belongs to l1r xor l2r.
(lr == l2r && lr != l1r)));
}
/*
std::cout << "l = " << l
<< " e = " << e
<< " (l1, l2) = (" << l1 << ", " << l2 << ") => "
<< " merging R" << l1r << " and R" << l2r
<< std::endl;
*/
// Merging both regions.
{
if (l2r < l1r)
std::swap(l1r, l2r);
mln_invariant(l2r > l1r);
lpar[l1r] = l2r;
/*
std::cout << "q[l1r] = " << q[l1r] << std::endl;
std::cout << "q[l2r] = " << q[l2r] << std::endl;
*/
q[l2r].insert(q[l1r]);
q[l1r].clear();
/*
std::cout << "q[l2r] = " << q[l2r] << std::endl
<< std::endl;
*/
/*
// Displaying lpar.
{
std::cout << "lpar = ";
for (unsigned i = 1; i <= n_basins; ++i)
{
L l = v[i].second;
std::cout << l << "->" << lpar[l] << " ";
}
std::cout << std::endl;
}
*/
}
E e1 = edge[l1],
e2 = edge[l2];
if (e1 == null && e2 == null)
{
// New edge-component (singleton)
// Put differently: new edge-node!
edge[l1] = e;
edge[l2] = e;
// after:
// e
// / \
// l1 l2
}
else if (e1 != null && e2 != null)
{
// Two trees shall merge.
E e1r = find_root_e(z_epar, e1),
e2r = find_root_e(z_epar, e2);
mln_invariant(e1r != e2r); // Otherwise, there's a bug!
// before:
// e1r e2r
// |... |...
// e1 e2
// / \ / \
// l1 . l2 .
epar(e1r) = e; z_epar(e1r) = e;
epar(e2r) = e; z_epar(e2r) = e;
// after:
// e
// / \
// e1r e2r
// |... |...
// e1 e2
// / \ / \
// l1 . l2 .
}
else if (e1 != null && e2 == null)
{
E e1r = find_root_e(z_epar, e1);
// before:
// e1r
// |...
// e1 null
// / \ /
// l1 . l2
epar(e1r) = e; z_epar(e1r) = e;
edge[l2] = e;
// after:
// e
// / \
// e1r l2
// |...
// e1
// / \
// l1 .
}
else
{
mln_invariant(e1 == null && e2 != null);
E e2r = find_root_e(z_epar, e2);
epar(e2r) = e; z_epar(e2r) = e;
edge[l1] = e;
}
} // end of "for every region with increasing attribute"
// Displaying the "edge tree".
{
p_set<E> leaves;
// Region l -> edge e.