- Olena-patches

3166: Improve soft heap interface.
by Guillaume Lazzara 19 Jan '09

19 Jan '09

* mln/util/soft_heap.hh: rename take() as push(). * tests/util/soft_heap.cc: update test. --- milena/ChangeLog | 8 ++++++++ milena/mln/util/soft_heap.hh | 10 ++++------ milena/tests/util/soft_heap.cc | 2 +- 3 files changed, 13 insertions(+), 7 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index daed173..8201593 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,13 @@ 2009-01-19 Guillaume Lazzara <z(a)lrde.epita.fr> + Improve soft heap interface. + + * mln/util/soft_heap.hh: rename take() as push(). + + * tests/util/soft_heap.cc: update test. + +2009-01-19 Guillaume Lazzara <z(a)lrde.epita.fr> + Add implementation of Fibonacci heap. * headers.mk: add a new header to distribution. diff --git a/milena/mln/util/soft_heap.hh b/milena/mln/util/soft_heap.hh index 144110f..bff9d11 100644 --- a/milena/mln/util/soft_heap.hh +++ b/milena/mln/util/soft_heap.hh @@ -197,13 +197,11 @@ namespace mln /// Add a new element \p element. void push(const T& element); - /// Insert a site \p p (equivalent as 'push'). - void insert(const T& p); /// Merge \p sh with this heap. /// Be ware that after this call, \p sh will be empty. This heap will /// hold the elements which were part of \p sh. - void take(soft_heap<T,R>& sh); + void push(soft_heap<T,R>& sh); /// Returns the element with the lowest priority and remove it from the /// heap. @@ -227,14 +225,14 @@ namespace mln /// Reset the heap to an empty heap. Do *NOT* delete element which may /// have been inserted. - /// \sa take + /// \sa push void soft_clear_(); private: /// Merge a node \p q to this heap. - /// \sa take + /// \sa push void meld(node<T,R> *q); /// Update suffix_min pointer according to the new values inserted in @@ -655,7 +653,7 @@ namespace mln template <typename T, typename R> inline void - soft_heap<T,R>::take(soft_heap<T,R>& psh) + soft_heap<T,R>::push(soft_heap<T,R>& psh) { head<T,R> *head = psh.head_hook_(); while (head != 0) diff --git a/milena/tests/util/soft_heap.cc b/milena/tests/util/soft_heap.cc index d2bffcc..9ab1d8b 100644 --- a/milena/tests/util/soft_heap.cc +++ b/milena/tests/util/soft_heap.cc @@ -60,7 +60,7 @@ int main() fh2.push(p[i]); // Merge fh in fh2. - fh2.take(fh); + fh2.push(fh); // fh2 now holds both its elements and fh's. unsigned i = 0; -- 1.5.6.5

1 0

3165: Add implementation of Fibonacci heap.
by Guillaume Lazzara 19 Jan '09

19 Jan '09

* headers.mk: add a new header to distribution. * mln/util/fibonacci_heap.hh: new file. The implementation. * tests/unit_test/Makefile.am, * tests/unit_test/mln_util_fibonacci_heap.cc: new unit test. * tests/util/Makefile.am, * tests/util/fibonacci_heap.cc: new test. --- milena/ChangeLog | 14 + milena/headers.mk | 1 + milena/mln/util/fibonacci_heap.hh | 1001 +++++++++++++++++++++ milena/tests/unit_test/Makefile.am | 2 + milena/tests/unit_test/mln_util_fibonacci_heap.cc | 11 + milena/tests/util/Makefile.am | 12 +- milena/tests/util/fibonacci_heap.cc | 129 +++ 7 files changed, 1165 insertions(+), 5 deletions(-) create mode 100644 milena/mln/util/fibonacci_heap.hh create mode 100644 milena/tests/unit_test/mln_util_fibonacci_heap.cc create mode 100644 milena/tests/util/fibonacci_heap.cc diff --git a/milena/ChangeLog b/milena/ChangeLog index 135ea79..daed173 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,3 +1,17 @@ +2009-01-19 Guillaume Lazzara <z(a)lrde.epita.fr> + + Add implementation of Fibonacci heap. + + * headers.mk: add a new header to distribution. + + * mln/util/fibonacci_heap.hh: new file. The implementation. + + * tests/unit_test/Makefile.am, + * tests/unit_test/mln_util_fibonacci_heap.cc: new unit test. + + * tests/util/Makefile.am, + * tests/util/fibonacci_heap.cc: new test. + 2009-01-19 Thierry Geraud <thierry.geraud(a)lrde.epita.fr> Cleanup extension with image, forcing it to be const. diff --git a/milena/headers.mk b/milena/headers.mk index a526da5..a730437 100644 --- a/milena/headers.mk +++ b/milena/headers.mk @@ -62,6 +62,7 @@ mln/registration/registration.hh \ mln/registration/essential.hh \ mln/registration/icp.hh \ mln/util/graph.hh \ +mln/util/fibonacci_heap.hh \ mln/util/max.hh \ mln/util/lazy_set.hh \ mln/util/soft_heap.hh \ diff --git a/milena/mln/util/fibonacci_heap.hh b/milena/mln/util/fibonacci_heap.hh new file mode 100644 index 0000000..8a7ac2c --- /dev/null +++ b/milena/mln/util/fibonacci_heap.hh @@ -0,0 +1,1001 @@ +// Copyright (C) 2009 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. +// +//*************************************************************************** +// The Fibonacci heap implementation is Copyright (c) 1996 by John Boyer +// +// Once this Fibonacci heap implementation (the software) has been published +// by Dr. Dobb's Journal, permission to use and distribute the software is +// granted provided that this copyright notice remains in the source and +// and the author (John Boyer) is acknowledged in works that use this program. +// +// Every effort has been made to ensure that this implementation is free of +// errors. Nonetheless, the author (John Boyer) assumes no liability regarding +// your use of this software. +// +// The author would also be very glad to hear from anyone who uses the +// software or has any feedback about it. +// Email: jboyer(a)gulf.csc.uvic.ca +//*************************************************************************** + + +#ifndef MLN_UTIL_FIBONACCI_HEAP_HH +# define MLN_UTIL_FIBONACCI_HEAP_HH + + +# include <iostream> +# include <mln/core/concept/object.hh> +# include <mln/util/ord.hh> + + + +namespace mln +{ + + namespace util + { + + + namespace internal + { + + /*--------------------------\ + | Fibonacci Heap node Class | + `--------------------------*/ + + template <typename T> + class fibonacci_heap_node + { + + public: + /// Default constructor. + fibonacci_heap_node(); + + /// Construct a new node with \p value as value. + fibonacci_heap_node(const T& value); + + ~fibonacci_heap_node(); + + /// Return the node value. + const T& value() const; + + fibonacci_heap_node<T> *left() const; + fibonacci_heap_node<T> *right() const; + fibonacci_heap_node<T> *parent() const; + fibonacci_heap_node<T> *child() const; + + short degree() const; + + short mark() const; + + + void set_value(const T&); + void set_left(fibonacci_heap_node<T> *node); + void set_right(fibonacci_heap_node<T> *node); + void set_parent(fibonacci_heap_node<T> *node); + void set_child(fibonacci_heap_node<T> *node); + void set_degree(short degree); + void set_mark(short mark); + + void operator=(const T& val); + + void operator=(fibonacci_heap_node<T>& rhs); + int operator==(fibonacci_heap_node<T>& rhs); + int operator<(fibonacci_heap_node<T>& rhs); + + void print_(); + + + private: + + fibonacci_heap_node<T> *left_; + fibonacci_heap_node<T> *right_; + fibonacci_heap_node<T> *parent_; + fibonacci_heap_node<T> *child_; + short degree_; + short mark_; + T value_; + }; + + } // end of namespace mln::util::internal + + + + /*---------------------\ + | Fibonacci Heap Class | + `---------------------*/ + + template <typename T> + class fibonacci_heap : public Object< fibonacci_heap<T> > + { + public: + + /// Default constructor + fibonacci_heap(); + + ~fibonacci_heap(); + + /// Push a new element in the heap. + /// \sa insert + void push(const T& value); + + /// Take \p other_heap's elements and insert them in this heap. + /// After this call \p other_heap is cleared. + void push(fibonacci_heap<T>& other_heap); + + /// Return the minimum value in the heap. + const T& front() const; + + /// Return and remove the minimum value in the heap. + T pop_front(); + + /// Is it empty? + bool is_empty() const; + + /// return false if it is empty. + bool is_valid() const; + + /// Return the number of elements. + unsigned nelements() const; + + /// Clear all elements in the heap and make the heap empty. + void clear(); + + + void print_(internal::fibonacci_heap_node<T> *tree = 0, + internal::fibonacci_heap_node<T> *parent = 0); + + + private: + + // Internal functions that help to implement the Standard Operations + + + /// Allow to change a node value. + /// FIXME: cannot use that function efficiently except by passing the + /// node pointer. Any idea? + /// FIXME: may be part of the public interface. + int decrease_key(internal::fibonacci_heap_node<T> *node, + internal::fibonacci_heap_node<T>& key); + + /// Remove node \p node from the child list of its parent node. + /// FIXME: cannot use that function efficiently except by passing the + /// node pointer. Any idea? + /// FIXME: may be part of the public interface. + int remove(internal::fibonacci_heap_node<T> *node); + + /// Insert a new node \p node in the heap. + void insert(internal::fibonacci_heap_node<T> *node); + + void exchange(internal::fibonacci_heap_node<T>*&, + internal::fibonacci_heap_node<T>*&); + + /// Internal function that reorganizes heap as part of an pop_front(). + /// We must find new minimum in heap, which could be anywhere along the + /// root list. The search could be O(n) since all nodes could be on + /// the root list. So, we reorganize the tree into more of a binomial + /// forest structure, and then find the new minimum on the consolidated + /// O(lg n) sized root list, and in the process set each Parent pointer + /// to 0, and count the number of resulting subtrees. + /// + /// Note that after a list of n inserts, there will be n nodes on the + /// root list, so the first pop_front() will be O(n) regardless of + /// whether or not we consolidate. However, the consolidation causes + /// subsequent pop_front() operations to be O(lg n). Furthermore, + /// the extra cost of the first pop_front() is covered by the higher + /// amortized cost of insert(), which is the real governing factor in + /// how costly the first pop_front() will be. + void consolidate(); + + /// The node \p lhs is removed from the root list and becomes a subtree + /// of node \p rhs. + void link(internal::fibonacci_heap_node<T> *lhs, + internal::fibonacci_heap_node<T> *rhs); + + void add_to_root_list(internal::fibonacci_heap_node<T> *); + + /// Remove node \p lhs from the child list of its parent node \p rhs. + void cut(internal::fibonacci_heap_node<T> *lhs, + internal::fibonacci_heap_node<T> *rhs); + + /// Cuts each node in parent list, putting successive ancestor nodes on + /// the root list until we either arrive at the root list or until we + /// find an ancestor that is unmarked. When a mark is set (which only + /// happens during a cascading cut), it means that one child subtree has + /// been lost; if a second subtree is lost later during another + /// cascading cut, then we move the node to the root list so that it + /// can be re-balanced on the next consolidate. + void cascading_cut(internal::fibonacci_heap_node<T> *); + + /// Clear the heap but do *NOT* delete elements. + void soft_clear(); + + internal::fibonacci_heap_node<T> *min_root; + long num_nodes; + long num_trees; + long num_marked_nodes; + + }; + + + +# ifndef MLN_INCLUDE_ONLY + + + namespace internal + { + + + /*---------------\ + | fibonacci_heap_node | + `---------------*/ + + + template <typename T> + inline + fibonacci_heap_node<T>::fibonacci_heap_node() + : left_(0), right_(0), parent_(0), child_(0), + degree_(0), mark_(0) + { + } + + + template <typename T> + inline + fibonacci_heap_node<T>::fibonacci_heap_node(const T& new_value) + : left_(0), right_(0), parent_(0), child_(0), + degree_(0), mark_(0), value_(new_value) + { + } + + + template <typename T> + inline + fibonacci_heap_node<T>::~fibonacci_heap_node() + { + } + + + template <typename T> + inline + const T& + fibonacci_heap_node<T>::value() const + { + return value_; + } + + + template <typename T> + inline + fibonacci_heap_node<T> * + fibonacci_heap_node<T>::left() const + { + return left_; + } + + + template <typename T> + inline + fibonacci_heap_node<T> * + fibonacci_heap_node<T>::right() const + { + return right_; + } + + + template <typename T> + inline + fibonacci_heap_node<T> * + fibonacci_heap_node<T>::parent() const + { + return parent_; + } + + + template <typename T> + inline + fibonacci_heap_node<T> * + fibonacci_heap_node<T>::child() const + { + return child_; + } + + + template <typename T> + inline + short + fibonacci_heap_node<T>::degree() const + { + return degree_; + } + + + template <typename T> + inline + short + fibonacci_heap_node<T>::mark() const + { + return mark_; + } + + + template <typename T> + inline + void + fibonacci_heap_node<T>::set_value(const T& value) + { + value_ = value; + } + + + template <typename T> + inline + void + fibonacci_heap_node<T>::set_left(fibonacci_heap_node<T> *node) + { + left_ = node; + } + + + template <typename T> + inline + void + fibonacci_heap_node<T>::set_right(fibonacci_heap_node<T> *node) + { + right_ = node; + } + + + template <typename T> + inline + void + fibonacci_heap_node<T>::set_parent(fibonacci_heap_node<T> *node) + { + parent_ = node; + } + + + template <typename T> + inline + void + fibonacci_heap_node<T>::set_child(fibonacci_heap_node<T> *node) + { + child_ = node; + } + + + template <typename T> + inline + void + fibonacci_heap_node<T>::set_degree(short degree) + { + degree_ = degree; + } + + + template <typename T> + inline + void + fibonacci_heap_node<T>::set_mark(short mark) + { + mark_ = mark; + } + + + template <typename T> + inline + void + fibonacci_heap_node<T>::operator=(const T& value) + { + value_ = value; + } + + + template <typename T> + inline + void fibonacci_heap_node<T>::operator=(fibonacci_heap_node<T>& rhs) + { + value_ = rhs.value(); + } + + + template <typename T> + inline + int + fibonacci_heap_node<T>::operator==(fibonacci_heap_node<T>& rhs) + { + return value_ == rhs.value() ? 1 : 0; + } + + + template <typename T> + inline + int + fibonacci_heap_node<T>::operator<(fibonacci_heap_node<T>& rhs) + { + return util::ord_strict(value_, rhs.value()) ? 1 : 0; + } + + + template <typename T> + inline + void fibonacci_heap_node<T>::print_() + { + std::cout << value_; + } + + + } // end of namespace mln::util::internal + + + + /*----------\ + | fibonacci_heap | + `----------*/ + + template <typename T> + inline + fibonacci_heap<T>::fibonacci_heap() + { + min_root = 0; + num_nodes = 0; + num_trees = 0; + num_marked_nodes = 0; + } + + + template <typename T> + inline + fibonacci_heap<T>::~fibonacci_heap() + { + clear(); + } + + + template <typename T> + inline + void + fibonacci_heap<T>::push(const T& value) + { + typedef internal::fibonacci_heap_node<T> node_t; + node_t *new_node = new node_t(value); + + insert(new_node); + } + + + template <typename T> + inline + void + fibonacci_heap<T>::push(fibonacci_heap<T>& other_heap) + { + internal::fibonacci_heap_node<T> *min1, *min2, *next1, *next2; + + if (other_heap.is_empty()) + return; + + if (min_root != 0) + { + // We join the two circular lists by cutting each list between its + // min node and the node after the min. This code just pulls those + // nodes into temporary variables so we don't get lost as changes + // are made. + min1 = min_root; + min2 = other_heap.min_root; + next1 = min1->right(); + next2 = min2->right(); + + // To join the two circles, we join the minimum nodes to the next + // nodes on the opposite chains. Conceptually, it looks like the way + // two bubbles join to form one larger bubble. They meet at one point + // of contact, then expand out to make the bigger circle. + min1->set_right(next2); + next2->set_left(min1); + min2->set_right(next1); + next1->set_left(min2); + + // Choose the new minimum for the heap + if (*min2 < *min1) + min_root = min2; + } else + min_root = other_heap.min_root; + + // Set the amortized analysis statistics and size of the new heap + num_nodes += other_heap.num_nodes; + num_marked_nodes += other_heap.num_marked_nodes; + num_trees += other_heap.num_trees; + + // Complete the union by setting the other heap to emptiness + other_heap.min_root = 0; + other_heap.num_nodes = 0; + other_heap.num_trees = 0; + other_heap.num_marked_nodes = 0; + } + + + template <typename T> + inline + const T& + fibonacci_heap<T>::front() const + { + return min_root->value(); + } + + + template <typename T> + inline + T + fibonacci_heap<T>::pop_front() + { + mln_precondition(is_valid()); + mln_precondition(!is_empty()); + + internal::fibonacci_heap_node<T> *result = min_root; + fibonacci_heap<T> *child_heap = 0; + + // Remove minimum node and set min_root to next node + min_root = result->right(); + result->right()->set_left(result->left()); + result->left()->set_right(result->right()); + result->set_left(0); + result->set_right(0); + + num_nodes --; + if (result->mark()) + { + --num_marked_nodes; + result->set_mark(0); + } + result->set_degree(0); + + // Attach child list of minimum node to the root list of the heap + // If there is no child list, then do no work + if (result->child() == 0) + { + if (min_root == result) + min_root = 0; + } + + // If min_root==result then there was only one root tree, so the + // root list is simply the child list of that node (which is + // 0 if this is the last node in the list) + else if (min_root == result) + min_root = result->child(); + + // If min_root is different, then the child list is pushed into a + // new temporary heap, which is then merged by union() onto the + // root list of this heap. + else + { + child_heap = new fibonacci_heap<T>(); + child_heap->min_root = result->child(); + } + + // Complete the disassociation of the result node from the heap + if (result->child() != 0) + result->child()->set_parent(0); + result->set_child(0); + result->set_parent(0); + + // If there was a child list, then we now merge it with the + // rest of the root list + if (child_heap) + { + push(*child_heap); + delete child_heap; + } + + // Consolidate heap to find new minimum and do reorganize work + if (min_root != 0) + consolidate(); + + // Return the minimum node, which is now disassociated with the heap + // It has left, right, parent, child, mark and degree cleared. + T val = result->value(); + delete result; + + return val; + } + + + template <typename T> + inline + int + fibonacci_heap<T>::decrease_key(internal::fibonacci_heap_node<T> *node, + internal::fibonacci_heap_node<T>& key) + { + internal::fibonacci_heap_node<T> *parent; + + if (node == 0 || *node < key) + return -1; + + *node = key; + + parent = node->parent(); + if (parent != 0 && *node < *parent) + { + cut(node, parent); + cascading_cut(parent); + } + + if (*node < *min_root) + min_root = node; + + return 0; + } + + + template <typename T> + inline + int + fibonacci_heap<T>::remove(internal::fibonacci_heap_node<T> *node) + { + internal::fibonacci_heap_node<T> temp; + int result; + + if (node == 0) + return -1; + + result = decrease_key(node, temp); + + if (result == 0) + if (pop_front() == 0) + result = -1; + + if (result == 0) + delete node; + + return result; + } + + + template <typename T> + inline + bool + fibonacci_heap<T>::is_empty() const + { + return min_root == 0; + } + + + template <typename T> + inline + bool + fibonacci_heap<T>::is_valid() const + { + return min_root != 0; + } + + + template <typename T> + inline + void + fibonacci_heap<T>::clear() + { + internal::fibonacci_heap_node<T> *temp; + + for (int i = 0; i < num_nodes; ++i) + { + temp = min_root->right(); + delete min_root; + min_root = temp; + } + + soft_clear(); + } + + + template <typename T> + inline + void + fibonacci_heap<T>::soft_clear() + { + min_root = 0; + num_nodes = 0; + num_trees = 0; + num_marked_nodes = 0; + } + + + template <typename T> + inline + unsigned + fibonacci_heap<T>::nelements() const + { + return num_nodes; + }; + + + template <typename T> + inline + void + fibonacci_heap<T>::print_(internal::fibonacci_heap_node<T> *tree, + internal::fibonacci_heap_node<T> *parent) + { + internal::fibonacci_heap_node<T>* temp = 0; + + if (tree == 0) + tree = min_root; + + temp = tree; + if (temp != 0) + { + do { + if (temp->left() == 0) + std::cout << "(left is 0)"; + temp->print_(); + if (temp->parent() != parent) + std::cout << "(parent is incorrect)"; + if (temp->right() == 0) + std::cout << "(right is 0)"; + else if (temp->right()->left() != temp) + std::cout << "(Error in left link left) ->"; + else + std::cout << " <-> "; + + temp = temp->right(); + + } while (temp != 0 && temp != tree); + } + else + std::cout << " <empty>" << std::endl; + std::cout << std::endl; + + temp = tree; + if (temp != 0) + { + do { + std::cout << "children of "; + temp->print_(); + std::cout << ": "; + if (temp->child() == 0) + std::cout << "NONE" << std::endl; + else print_(temp->child(), temp); + temp = temp->right(); + } while (temp!=0 && temp != tree); + } + } + + + template <typename T> + inline + void fibonacci_heap<T>::consolidate() + { + internal::fibonacci_heap_node<T> *x, *y, *w; + internal::fibonacci_heap_node<T> *a[1 + 8 * sizeof (long)]; // 1+lg(n) + short dn = 1 + 8 * sizeof (long); + + // Initialize the consolidation detection array + for (int i = 0; i < dn; i++) + a[i] = 0; + + // We need to loop through all elements on root list. + // When a collision of degree is found, the two trees + // are consolidated in favor of the one with the lesser + // element key value. We first need to break the circle + // so that we can have a stopping condition (we can't go + // around until we reach the tree we started with + // because all root trees are subject to becoming a + // child during the consolidation). + min_root->left()->set_right(0); + min_root->set_left(0); + w = min_root; + + short d; + do { + x = w; + d = x->degree(); + w = w->right(); + + // We need another loop here because the consolidated result + // may collide with another large tree on the root list. + while (a[d] != 0) + { + y = a[d]; + if (*y < *x) + exchange(x, y); + if (w == y) w = y->right(); + link(y, x); + a[d] = 0; + d++; + } + a[d] = x; + + } while (w != 0); + + // Now we rebuild the root list, find the new minimum, + // set all root list nodes' parent pointers to 0 and + // count the number of subtrees. + min_root = 0; + num_trees = 0; + for (int i = 0; i < dn; i++) + if (a[i] != 0) + add_to_root_list(a[i]); + } + + + template <typename T> + inline + void + fibonacci_heap<T>::link(internal::fibonacci_heap_node<T> *y, + internal::fibonacci_heap_node<T> *x) + { + // Remove node y from root list + if (y->right() != 0) + y->right()->set_left(y->left()); + if (y->left() != 0) + y->left()->set_right(y->right()); + num_trees--; + + // Make node y a singleton circular list with a parent of x + y->set_left(y); + y->set_right(y); + y->set_parent(x); + + // If node x has no children, then list y is its new child list + if (x->child() == 0) + x->set_child(y); + + // Otherwise, node y must be added to node x's child list + else + { + y->set_left(x->child()); + y->set_right(x->child()->right()); + x->child()->set_right(y); + y->right()->set_left(y); + } + + // Increase the degree of node x because it's now a bigger tree + x->set_degree(x->degree() + 1); + + // node y has just been made a child, so clear its mark + if (y->mark()) num_marked_nodes--; + y->set_mark(0); + } + + + template <typename T> + inline + void + fibonacci_heap<T>::add_to_root_list(internal::fibonacci_heap_node<T> *x) + { + if (x->mark()) + --num_marked_nodes; + x->set_mark(0); + + num_nodes--; + insert(x); + } + + + template <typename T> + inline + void + fibonacci_heap<T>::cut(internal::fibonacci_heap_node<T> *x, + internal::fibonacci_heap_node<T> *y) + { + if (y->child() == x) + y->child() = x->right(); + if (y->child() == x) + y->child() = 0; + + y->degree --; + + x->left()->right() = x->right(); + x->right()->left() = x->left(); + + add_to_root_list(x); + } + + + template <typename T> + inline + void + fibonacci_heap<T>::cascading_cut(internal::fibonacci_heap_node<T> *y) + { + internal::fibonacci_heap_node<T> *z = y->parent(); + + while (z != 0) + { + if (y->mark() == 0) + { + y->mark() = 1; + num_marked_nodes++; + z = 0; + } + else + { + cut(y, z); + y = z; + z = y->parent(); + } + } + } + + + template <typename T> + inline + void + fibonacci_heap<T>::insert(internal::fibonacci_heap_node<T> *node) + { + if (node == 0) + return; + + // If the heap is currently empty, then new node becomes singleton + // circular root list + if (min_root == 0) + { + min_root = node; + node->set_left(node); + node->set_right(node); + } + else + { + // Pointers from node set to insert between min_root and + // min_root->right() + node->set_right(min_root->right()); + node->set_left(min_root); + + // Set Pointers to node + node->left()->set_right(node); + node->right()->set_left(node); + + // The new node becomes new min_root if it is less than current + // min_root + if (*node < *min_root) + min_root = node; + } + + // We have one more node in the heap, and it is a tree on the root list + num_nodes++; + num_trees++; + node->set_parent(0); + } + + + template <typename T> + inline + void fibonacci_heap<T>::exchange(internal::fibonacci_heap_node<T>*& n1, + internal::fibonacci_heap_node<T>*& n2) + { + internal::fibonacci_heap_node<T> *temp; + + temp = n1; + n1 = n2; + n2 = temp; + } + +# endif // ! MLN_INCLUDE_ONLY + + + + } // end of namespace mln::util + +} // end of namespace mln + +#endif // ! MLN_UTIL_FIBONACCI_HEAP_HH diff --git a/milena/tests/unit_test/Makefile.am b/milena/tests/unit_test/Makefile.am index 1d3d7cd..b5dd047 100644 --- a/milena/tests/unit_test/Makefile.am +++ b/milena/tests/unit_test/Makefile.am @@ -62,6 +62,7 @@ mln_registration_registration \ mln_registration_essential \ mln_registration_icp \ mln_util_graph \ +mln_util_fibonacci_heap \ mln_util_max \ mln_util_lazy_set \ mln_util_soft_heap \ @@ -1076,6 +1077,7 @@ mln_registration_registration_SOURCES = mln_registration_registration.cc mln_registration_essential_SOURCES = mln_registration_essential.cc mln_registration_icp_SOURCES = mln_registration_icp.cc mln_util_graph_SOURCES = mln_util_graph.cc +mln_util_fibonacci_heap_SOURCES = mln_util_fibonacci_heap.cc mln_util_max_SOURCES = mln_util_max.cc mln_util_lazy_set_SOURCES = mln_util_lazy_set.cc mln_util_soft_heap_SOURCES = mln_util_soft_heap.cc diff --git a/milena/tests/unit_test/mln_util_fibonacci_heap.cc b/milena/tests/unit_test/mln_util_fibonacci_heap.cc new file mode 100644 index 0000000..d143ea9 --- /dev/null +++ b/milena/tests/unit_test/mln_util_fibonacci_heap.cc @@ -0,0 +1,11 @@ +// Unit test for mln/util/fibonacci_heap.hh. +// Generated by ./build_unit_test.sh, do not modify. + +// Include the file twice, so we detect missing inclusion guards. +#include <mln/util/fibonacci_heap.hh> +#include <mln/util/fibonacci_heap.hh> + +int main() +{ + // Nothing. +} diff --git a/milena/tests/util/Makefile.am b/milena/tests/util/Makefile.am index 1dcd5c5..c095994 100644 --- a/milena/tests/util/Makefile.am +++ b/milena/tests/util/Makefile.am @@ -6,32 +6,34 @@ include $(top_srcdir)/milena/tests/tests.mk check_PROGRAMS = \ branch_iter \ branch_iter_ind \ - line_graph \ eat \ + fibonacci_heap \ graph \ lazy_set \ lemmings \ + line_graph \ ord_pair \ soft_heap \ tree \ tree_fast \ tree_fast_to_image \ - tree_to_image \ - tree_to_fast + tree_to_fast \ + tree_to_image branch_iter_SOURCES = branch_iter.cc branch_iter_ind_SOURCES = branch_iter_ind.cc -line_graph_SOURCES = line_graph.cc eat_SOURCES = eat.cc +fibonacci_heap_SOURCES = fibonacci_heap.cc graph_SOURCES = graph.cc lazy_set_SOURCES = lazy_set.cc lemmings_SOURCES = lemmings.cc +line_graph_SOURCES = line_graph.cc ord_pair_SOURCES = ord_pair.cc soft_heap_SOURCES = soft_heap.cc tree_SOURCES = tree.cc tree_fast_SOURCES = tree_fast.cc tree_fast_to_image_SOURCES = tree_to_image.cc -tree_to_image_SOURCES = tree_to_image.cc tree_to_fast_SOURCES = tree_to_fast.cc +tree_to_image_SOURCES = tree_to_image.cc TESTS = $(check_PROGRAMS) diff --git a/milena/tests/util/fibonacci_heap.cc b/milena/tests/util/fibonacci_heap.cc new file mode 100644 index 0000000..a7ce896 --- /dev/null +++ b/milena/tests/util/fibonacci_heap.cc @@ -0,0 +1,129 @@ +// Copyright (C) 2009 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/util/fibonacci_heap.cc +/// +/// Definition of a generic vector class. + +#include <mln/util/fibonacci_heap.hh> +#include <mln/core/alias/point2d.hh> + +using mln::point2d; + +point2d p[] = { point2d(4,5), point2d(3,4), point2d(3,2), + point2d(1,6), point2d(2,3), point2d(3,5), + point2d(2,4), point2d(7,2), point2d(9,6), + point2d(7,3) }; + +point2d res_1[] = { point2d(1,6), point2d(2,3), point2d(2,4) }; + +point2d res_2[] = { point2d(1,4), point2d(3,2), point2d(3,4), + point2d(3,5), point2d(4,5), point2d(5,4), + point2d(7,2), point2d(7,3), point2d(9,6), + point2d(53,4) }; + + +int main() +{ + using namespace mln; + + /// Init heap + util::fibonacci_heap<point2d> heap; + for (unsigned i = 0; i < 5; ++i) + heap.push(p[i]); + + /// Init heap2 + util::fibonacci_heap<point2d> heap2; + for (unsigned i = 5; i < 10; ++i) + heap2.push(p[i]); + + + /*! + ** + ** | | | | |x| | | + ** | | | | |x| | | + ** |x| + |x| => |x| + | | + ** |x| |x| |x| | | + ** - - - - + ** heap heap2 heap heap2 + ** + */ + /// Merge heap and heap2. + heap.push(heap2); + + /// Heap2 is now empty + mln_assertion(heap2.is_empty()); + mln_assertion(heap2.nelements() == 0); + /// Check if the front() is correct in heap + mln_assertion(heap.front() == res_1[0]); + + + /*! + ** + ** |x| | | | | |x| + ** |x| | | | | |x| + ** |x| + | | => | | + |x| + ** |x| | | | | |x| + ** - - - - + ** heap heap2 heap heap2 + ** + */ + /// Heap2 is empty and heap is full of elements. + /// Move the content of heap into heap2. + heap2.push(heap); + + /// Heap is now empty + mln_assertion(heap.is_empty()); + mln_assertion(heap.nelements() == 0); + + /// Check if the front() is correct + mln_assertion(heap2.front() == res_1[0]); + + + /// Extract and delete few front elements. + for (unsigned i = 0; i < 3; ++i) + mln_assertion(heap2.pop_front() == res_1[i]); + + + /// Re-insert data after deletion... + heap2.push(point2d(53,4)); + heap2.push(point2d(1,4)); + heap2.push(point2d(5,4)); + + /// ... and try to extract and delete data again. + unsigned i = 0; + while (heap2.is_valid()) + { + mln_assertion(heap2.pop_front() == res_2[i++]); + mln_assertion(heap2.nelements() == (10 - i)); + } + + /// The heap must be empty. + mln_assertion(heap2.is_empty()); + mln_assertion(heap2.nelements() == 0); +} -- 1.5.6.5

1 0

[Fwd: [LrdeTools] 3158: Commit test.]
by Thierry GERAUD 19 Jan '09

19 Jan '09

1 0

[Fwd: [LrdeTools] 3157: Add images for Dalila and create IGR directories .]
by Thierry GERAUD 19 Jan '09

19 Jan '09

1 0

3164: Cleanup extension with image, forcing it to be const.
by Thierry Geraud 19 Jan '09

19 Jan '09

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Thierry Geraud <thierry.geraud(a)lrde.epita.fr> Cleanup extension with image, forcing it to be const. * mln/core/image/extension_ima.hh (todo): Remove; obsolete. (ctor): Add const for extension. (change_extension): Remove. (extension): Remove the mutable version. (init_): Remove the mutable version. (init_): Make it robust to constness. * mln/core/routine/extend.hh: Upgrade file doc style. (extend): Remove overload for mutable extension. Misc. * mln/core/site_set/p_priority.hh: Upgrade file doc style. (insert): New overload for a priority queue. * mln/make/double_neighb2d.hh: Upgrade file doc style. * mln/make/box2d.hh (todo): New. * mln/labeling/compute.hh: Fix commentary. core/image/extension_ima.hh | 71 ++++++++++---------------------------------- core/routine/extend.hh | 31 ++++--------------- core/site_set/p_priority.hh | 32 ++++++++++++++++--- labeling/compute.hh | 2 - make/box2d.hh | 3 + make/double_neighb2d.hh | 13 +++----- 6 files changed, 61 insertions(+), 91 deletions(-) Index: mln/core/site_set/p_priority.hh --- mln/core/site_set/p_priority.hh (revision 3163) +++ mln/core/site_set/p_priority.hh (working copy) @@ -1,4 +1,5 @@ // Copyright (C) 2007, 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 @@ -28,12 +29,9 @@ #ifndef MLN_CORE_P_PRIORITY_HH # define MLN_CORE_P_PRIORITY_HH -/*! \file mln/core/site_set/p_priority.hh - * - * \brief Definition of a priority queue of sites. - * - * \todo Rename file as p_priority.hh - */ +/// \file mln/core/site_set/p_priority.hh +/// +/// Definition of a priority queue of sites. # include <map> # include <mln/core/site_set/p_double.hh> @@ -126,6 +124,9 @@ /// Insert a pair \p p_e (priority p, element e). void insert(const i_element& p_e); + /// Insert elements from another priority queue. + void insert(const p_priority<P,Q>& other); + /// Pop (remove) from the queue an element with highest priority. /// If several elements have this priority, the least recently @@ -175,6 +176,7 @@ /// Return the size of this site set in memory. std::size_t memory_size() const; + protected: typedef std::map<P, Q, util::ord<P> > q_type_; @@ -255,6 +257,24 @@ template <typename P, typename Q> inline void + p_priority<P,Q>::insert(const p_priority<P,Q>& other) + { + mln_invariant(run_()); + typename q_type_::const_iterator i; + for (i = other.q_.begin(); i != other.q_.end(); ++i) + { + P priority = i->first; + p_.insert(priority); + const Q& q_p = i->second; + q_[priority] += q_p; + } + n_ += other.n_; + mln_invariant(run_()); + } + + template <typename P, typename Q> + inline + void p_priority<P,Q>::pop() { mln_precondition(! this->is_empty()); // Also test invariants. Index: mln/core/image/extension_ima.hh --- mln/core/image/extension_ima.hh (revision 3163) +++ mln/core/image/extension_ima.hh (working copy) @@ -31,15 +31,9 @@ /// \file mln/core/image/extension_ima.hh /// /// Definition of a morpher that extends the domain of an image -/// with a function. -/// -/// \todo Use an envelop as lvalue to test extension writing. -/// -/// \todo Handle the couple of cases: either J is value_io::read_write -/// or value_io::read_only; then ext_io can be read_write... +/// with an image. # include <mln/core/internal/image_identity.hh> -# include <mln/data/fill_with_value.hh> @@ -57,7 +51,7 @@ template <typename I, typename J> struct data< extension_ima<I, J> > { - data(I& ima, J& ext); + data(I& ima, const J& ext); I ima_; J ext_; @@ -80,7 +74,7 @@ // extended domain typedef trait::image::ext_domain::extendable ext_domain; typedef trait::image::ext_value::multiple ext_value; - typedef trait::image::ext_io::read_only ext_io; // FIXME: Too restrictive? + typedef trait::image::ext_io::read_only ext_io; }; template <typename I, typename J, typename V> @@ -118,11 +112,11 @@ extension_ima(); /// Constructor from an image \p ima and a function \p ext. - extension_ima(I& ima, J& ext); + extension_ima(I& ima, const J& ext); /// Deferred initialization from an image \p ima and a function \p /// ext. - void init_(I& ima, J& ext); + void init_(I& ima, const J& ext); /// Test if \p p is valid. @@ -139,14 +133,7 @@ /// Read-only access to the extension domain (image). - mlc_const(J)& extension() const; - - /// Mutable access to the extension domain (image). This domain - /// can be modified if J a read-write image type. - J& extension(); - - /// Change the value in the extension domain (image). - void change_extension(const mln_value(I)& v); + const J& extension() const; }; @@ -158,10 +145,6 @@ template <typename J, typename I> void init_(tag::extension_t, J& target, const extension_ima<I,J>& model); - template <typename J, typename I> - void init_(tag::extension_t, J& target, const extension_ima<I,const J>& model); - - # ifndef MLN_INCLUDE_ONLY @@ -173,9 +156,9 @@ template <typename I, typename J> inline - data< extension_ima<I, J> >::data(I& ima, J& ext) + data< extension_ima<I, J> >::data(I& ima, const J& ext) : ima_(ima), - ext_(ext) + ext_(const_cast<J&>(ext)) { } @@ -191,7 +174,7 @@ template <typename I, typename J> inline - extension_ima<I, J>::extension_ima(I& ima, J& ext) + extension_ima<I, J>::extension_ima(I& ima, const J& ext) { init_(ima, ext); } @@ -199,7 +182,7 @@ template <typename I, typename J> inline void - extension_ima<I, J>::init_(I& ima, J& ext) + extension_ima<I, J>::init_(I& ima, const J& ext) { this->data_ = new internal::data< extension_ima<I, J> >(ima, ext); } @@ -255,35 +238,17 @@ template <typename I, typename J> inline - mlc_const(J)& + const J& extension_ima<I, J>::extension() const { mln_precondition(this->is_valid()); return this->data_->ext_; } - template <typename I, typename J> - inline - J& - extension_ima<I, J>::extension() - { - mln_precondition(this->is_valid()); - return this->data_->ext_; - } - - template <typename I, typename J> - inline - void - extension_ima<I, J>::change_extension(const mln_value(I)& v) - { - mlc_equal(mln_trait_image_value_io(J), - trait::image::value_io::read_write)::check(); - data::fill_with_value(v); - } - // init_ template <typename I, typename J, typename M> + inline void init_(tag::image_t, extension_ima<I,J>& target, const M& model) { I ima; @@ -294,15 +259,13 @@ } template <typename J, typename I> + inline void init_(tag::extension_t, J& target, const extension_ima<I,J>& model) { - target = model.extension(); - } - - template <typename J, typename I> - void init_(tag::extension_t, J& target, const extension_ima<I,const J>& model) - { - target = model.extension(); + typedef mlc_unconst(J) J_; + J_& ext_ = const_cast<J_&>(model.extension()); + J_& target_ = const_cast<J_&>(target); + target_ = ext_; } # endif // ! MLN_INCLUDE_ONLY Index: mln/core/routine/extend.hh --- mln/core/routine/extend.hh (revision 3163) +++ mln/core/routine/extend.hh (working copy) @@ -1,4 +1,4 @@ -// Copyright (C) 2008 EPITA Research and Development Laboratory +// 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 @@ -28,15 +28,13 @@ #ifndef MLN_CORE_ROUTINE_EXTEND_HH # define MLN_CORE_ROUTINE_EXTEND_HH -/*! - * \file mln/core/routine/extend.hh - * - * \brief Definition of a morpher that extends the extended domain of an image. - * - * \todo Use the 'instant' mechanism. - * \todo Deal with ambiguities. - * \todo Check that there is no extension yet (except "extendable"). - */ +/// \file mln/core/routine/extend.hh +/// +/// Definition of a morpher that extends the extended domain of an image. +/// +/// \todo Use the 'instant' mechanism. +/// \todo Deal with ambiguities. +/// \todo Check that there is no extension yet (except "extendable"). # include <mln/core/image/extension_ima.hh> # include <mln/core/image/extension_fun.hh> @@ -65,10 +63,6 @@ extension_ima<const I, const J> extend(const Image<I>& ima, const Image<J>& ext); - template <typename I, typename J> - extension_ima<I, J> - extend(Image<I>& ima, Image<J>& ext); - /// Routines for domain extension with a value. @@ -119,15 +113,6 @@ return tmp; } - template <typename I, typename J> - extension_ima<I, J> - extend(Image<I>& ima, Image<J>& ext) - { - mlc_converts_to(mln_value(J), mln_value(I))::check(); - extension_ima<I, J> tmp(exact(ima), exact(ext)); - return tmp; - } - // With a value. Index: mln/make/double_neighb2d.hh --- mln/make/double_neighb2d.hh (revision 3163) +++ mln/make/double_neighb2d.hh (working copy) @@ -1,4 +1,4 @@ -// Copyright (C) 2008 EPITA Research and Development Laboratory +// 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 @@ -28,12 +28,11 @@ #ifndef MLN_MAKE_DOUBLE_NEIGHB2D_HH # define MLN_MAKE_DOUBLE_NEIGHB2D_HH -/*! \file mln/make/double_neighb2d.hh - * - * \brief Routine to create a double neighborhood. - * - * \todo Add overload with 'when_*' being Neighborhood<N>... - */ +/// \file mln/make/double_neighb2d.hh +/// +/// Routine to create a double neighborhood. +/// +/// \todo Add overload with 'when_*' being Neighborhood<N>... # include <mln/convert/to.hh> # include <mln/core/alias/window2d.hh> Index: mln/make/box2d.hh --- mln/make/box2d.hh (revision 3163) +++ mln/make/box2d.hh (working copy) @@ -32,6 +32,9 @@ /// \file mln/make/box2d.hh /// /// Routines to construct an mln::box2d. +/// +/// \todo Accept long int as args and test overflow (out of def::coord +/// bounds). # include <mln/core/alias/box2d.hh> Index: mln/labeling/compute.hh --- mln/labeling/compute.hh (revision 3163) +++ mln/labeling/compute.hh (working copy) @@ -210,7 +210,7 @@ return convert::to< util::array<mln_result(A)> >(accus); } - } // end of namespace mln::labeling::impl::internal + } // end of namespace mln::labeling::impl::generic } // end of namespace mln::labeling::impl

1 0

3163: Update Laurent's ISMM code.
by Thierry Geraud 19 Jan '09

19 Jan '09

https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox Index: ChangeLog from Thierry Geraud <thierry.geraud(a)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. ismm2009.cc | 687 +++++++++++++++++++++++++++++++++++++++++++++++++--- ismm2009.hh | 33 -- ismm2009.v1.cc | 745 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 1403 insertions(+), 62 deletions(-) Index: laurent/ismm2009.cc --- laurent/ismm2009.cc (revision 3162) +++ laurent/ismm2009.cc (working copy) @@ -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 @@ + 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 @@ 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 @@ 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 @@ std::vector<pair_t> v = sort_attributes(a, n_basins); // Sort regions. - std::cout << "attributes = "; + // Display attributes. + { + std::cout << "(attribute, label) = { "; for (unsigned i = 1; i <= n_basins; ++i) - std::cout << v[i].first // Attribute (increasing). - << ':' + std::cout << '(' << v[i].first // Attribute (increasing). + << ',' << v[i].second // Region label. - << " - "; - std::cout << std::endl; + << ") "; + 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; + // 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) + { + // 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, // Region label. - lr = find_root(lpar, l); + 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; + } - point2d e = edge[l]; // FIXME: Use the heap! + 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); - 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! + 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::endl; + 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. + + std::cout << "region:(edge) = "; + for (L l = 1; l <= n_basins; ++l) + { + mln_invariant(edge[l] != null); + leaves.insert(edge[l]); + std::cout << l << ':' << edge[l] << " "; } + std::cout << std::endl + << std::endl; + + mln_piter_(p_set<E>) e(leaves); + + + // Tree "e -> epar(e)". + + std::cout << "edge tree: " << std::endl; + for_all(e) + { + E e_ = e; + do + { + std::cout << e_ << " -> "; + mln_invariant(aa(e_) != 0 && aa(epar(e_)) != 0); // aa is valid + mln_invariant(aa(epar(e_)) >= aa(e_)); + e_ = epar(e_); + } + while (epar(e_) != e_); + std::cout << e_ << std::endl; + } + std::cout << std::endl; + + + // Edge parenthood goes with 'a' increasing: + + std::cout << "edge tree new attributes: " << std::endl; + for_all(e) + { + E e_ = e; + do + { + std::cout << aa(e_) << " -> "; + // Edge parenthood goes with a increasing: + mln_invariant(aa(epar(e_)) >= aa(e_)); + e_ = epar(e_); + } + while (epar(e_) != e_); + std::cout << aa(e_) << std::endl; + } + std::cout << std::endl; + + } // end of tree display. + + + + // Testing both that all regions and all "smallest" edges have + // merged (assumption: connected domain!). + + { + L l_1 = v[1].second, + l_root = find_root_l(lpar, l_1); + mln_invariant(edge[l_1] != null); + E e_root = find_root_e(z_epar, edge[l_1]); + for (unsigned i = 2; i <= n_basins; ++i) + { + L l = v[i].second; + mln_invariant(find_root_l(lpar, l) == l_root); + mln_invariant(edge[l] != null); + mln_invariant(find_root_e(z_epar, edge[l]) == e_root); + } + } + + + // Reminders: + debug::println("wst(g) fully valuated:", wst_g_full); + debug::println("g_line:", g_line); + + + + // +---------------------------------------------------------------+ + // | We want an "image" of the "merging value of a" over all edges | + // | (not only the set of 'e' used to merge). | + // +---------------------------------------------------------------+ + + + { + + // Dealing with the watershed line. + + mln_VAR(aa_line, aa | is_line); + + { + // First, processing the 1D-faces (edges) of the watershed line. + + std::vector< std::vector<E> > lca; + lca = compute_lca(n_basins, edge, epar); // lca is auxiliary data. + + mln_piter_(g_line_t) e(g_line.domain()); + for_all(e) + { + L l1 = adj_line(e).first(), + l2 = adj_line(e).second(); + mln_invariant(l1 != 0 && l2 > l1); + E e_ = lca[l1][l2]; + mln_invariant(is_line(e_)); + mln_invariant(aa(e_) != 0); // aa is valid at e_. + + // The attribute value propagates from the lca to the current edge + // of the line: + aa(e) = aa(e_); + } + + debug::println("aa | wst(g) line:", (aa | is_edge) | is_line); + + // Second, processing the 0D-faces (points) of the watershed line. + + data::paste(morpho::dilation(extend(aa_line | (pw::value(aa_line) == pw::cst(0)), + aa_line), + c4().win()), + aa_line); + + /* + + // Version with "aa(x) = 0" to separate sheds (a shed = + // adjacency between a couple of basins). + + mln_piter_(aa_line_t) pxl(aa_line.domain()); + mln_niter_(neighb2d) ne(c4(), pxl); // FIXME: "pxl -> nbhoring edges" is explictly c4()... + for_all(pxl) + { + if (aa(pxl) != 0) + continue; + unsigned count = 0; + A na, na2; + for_all(ne) + if (aa_line.has(ne)) + { + mln_invariant(aa(ne) != 0); + if (count == 0) + na = aa(ne); // First attribute value encountered. + else + na2 = aa(ne); // Second (or more) attr value encountered. + ++count; + } + if (count == 2) + { + mln_invariant(na == na2); // Both attribute values have to be the same. + aa(pxl) = na; + } + } + */ + + + } + + debug::println("aa | line:", aa_line); + + + // Now dealing with all elements of basins: + + // - 2D-faces, i.e., original pixels; + // - 1D-faces, i.e., edges within regions; + // - 0D-faces, i.e., points within regions. + + mln_VAR(aa_basins, aa | (pw::value(wst_g_full) != pw::cst(0))); + + { + mln_piter_(aa_basins_t) p(aa_basins.domain()); + for_all(p) + aa(p) = a[wst_g_full(p)]; + } + + debug::println("aa | basins", aa_basins); + + } + + + + debug::println("aa:", aa); + } // end of main Index: laurent/ismm2009.v1.cc --- laurent/ismm2009.v1.cc (revision 0) +++ laurent/ismm2009.v1.cc (revision 0) @@ -0,0 +1,745 @@ +#include "ismm2009.hh" + +#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> + + + +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 Q, typename L, typename W> + mln_element(Q) + get_smallest_edge(Q& q, L l, const W& wst, std::vector<L>& lpar, bool& found) + { + L lr = find_root_l(lpar, l); + typedef mln_element(Q) E; + + mln_fwd_piter(Q) e(q); + for_all(e) + { + 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); + + if (l1r == l2r) + // 'e' is an internal edge => forget it. + continue; + + if (l1r == lr || l2r == lr) + { + found = true; + return e; + } + } + + found = false; + return E(); + } + + + // Compute an LCA; reference(?) code. + + template <typename L, typename E, typename Epar> + image2d<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()); + } + + + box2d b = make::box2d(1, 1, l_max, l_max); + image2d<E> lca(b); + data::fill(lca, E(0,0)); + // Superior right. + for (int l1 = 1; l1 <= l_max; ++l1) + for (int l2 = l1 + 1; l2 <= l_max; ++l2) + { + E e = edge[l1]; + while (chl(e).find(l2) == chl(e).end()) + e = epar(e); + lca.at_(l1, l2) = e; + } + // Diagonal. + for (int l = 1; l <= l_max; ++l) + lca.at_(l, l) = edge[l]; // Itself. + + // debug::println(lca); + + return lca; + } + + +} // mln + + + + +void usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " input.pgm" << std::endl; + std::cerr << "Laurent ISMM 2009 scheme." << std::endl; + abort(); +} + + + + +int main(int argc, char* argv[]) +{ + using namespace mln; + using value::int_u8; + using value::label_16; + + if (argc != 2) + usage(argv); + + image2d<int_u8> raw_f; + io::pgm::load(raw_f, argv[1]); + + 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 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) ); + + 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("adj:", adj); + mln_VAR(adj_line, adj | (pw::value(wst_g) == pw::cst(0))); + debug::println("adjacency:", adj_line); + + + + image2d<L> wst_g_full = wst_g.unmorph_(); + { + // edges (1D-faces) -> pixels (2D-faces) + mln_VAR(w_pixels, wst_g_full | is_pixel); + data::paste(morpho::dilation(extend(w_pixels, pw::value(wst_g_full)), + c4().win()), + wst_g_full); + // edges (1D-faces) -> points (0D-faces) + mln_VAR(w_points, wst_g_full | is_point); + data::paste(morpho::erosion(extend(w_points, pw::value(wst_g_full)), + c4().win()), + 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); + + mln_VAR(g_line, g | is_line); + debug::println("g | wst(g) line:", g_line); + + + + + + // Priority queue -> edge. + + typedef p_priority< T, p_queue<E> > Q; + Q q; + + { + mln_piter_(g_t) e(g.domain()); + for_all(e) + if (wst_g(e) == 0) + q.push(mln_max(T) - g(e), e); + // --- std::cout << "q = " << q << 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. +// { +// bool found; +// for (L l = 1; l <= n_basins; ++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(), + g, + wst_g, + n_basins); + + typedef std::pair<A,L> pair_t; + std::vector<pair_t> v = sort_attributes(a, n_basins); // Sort regions. + + + // 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; + } + + + + + // 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) + { + // 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; // 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; + + /* + // 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. + + std::cout << "region:(edge) = "; + for (L l = 1; l <= n_basins; ++l) + { + mln_invariant(edge[l] != null); + leaves.insert(edge[l]); + std::cout << l << ':' << edge[l] << " "; + } + std::cout << std::endl + << std::endl; + + mln_piter_(p_set<E>) e(leaves); + + + // Tree "e -> epar(e)". + + std::cout << "edge tree: " << std::endl; + for_all(e) + { + E e_ = e; + do + { + std::cout << e_ << " -> "; + mln_invariant(aa(e_) != 0 && aa(epar(e_)) != 0); // aa is valid + mln_invariant(aa(epar(e_)) >= aa(e_)); + e_ = epar(e_); + } + while (epar(e_) != e_); + std::cout << e_ << std::endl; + } + std::cout << std::endl; + + + // Edge parenthood goes with 'a' increasing: + + std::cout << "edge tree new attributes: " << std::endl; + for_all(e) + { + E e_ = e; + do + { + std::cout << aa(e_) << " -> "; + // Edge parenthood goes with a increasing: + mln_invariant(aa(epar(e_)) >= aa(e_)); + e_ = epar(e_); + } + while (epar(e_) != e_); + std::cout << aa(e_) << std::endl; + } + std::cout << std::endl; + + } // end of tree display. + + + + // Testing both that all regions and all "smallest" edges have + // merged (assumption: connected domain!). + + { + L l_1 = v[1].second, + l_root = find_root_l(lpar, l_1); + mln_invariant(edge[l_1] != null); + E e_root = find_root_e(z_epar, edge[l_1]); + for (unsigned i = 2; i <= n_basins; ++i) + { + L l = v[i].second; + mln_invariant(find_root_l(lpar, l) == l_root); + mln_invariant(edge[l] != null); + mln_invariant(find_root_e(z_epar, edge[l]) == e_root); + } + } + + + // Reminders: + debug::println("wst(g) fully valuated:", wst_g_full); + debug::println("g_line:", g_line); + + + + // +---------------------------------------------------------------+ + // | We want an "image" of the "merging value of a" over all edges | + // | (not only the set of 'e' used to merge). | + // +---------------------------------------------------------------+ + + + { + + // Dealing with the watershed line. + + mln_VAR(aa_line, aa | is_line); + + { + // First, processing the 1D-faces (edges) of the watershed line. + + image2d<E> lca = compute_lca(n_basins, edge, epar); // lca is auxiliary data. + + mln_piter_(g_line_t) e(g_line.domain()); + for_all(e) + { + E e_ = lca.at_(adj_line(e).first(), adj_line(e).second()); + mln_invariant(is_line(e_)); + mln_invariant(aa(e_) != 0); // aa is valid at e_. + + // The attribute value propagates from the lca to the current edge + // of the line: + aa(e) = aa(e_); + } + + debug::println("aa | wst(g) line:", (aa | is_edge) | is_line); + + // Second, processing the 0D-faces (points) of the watershed line. + + mln_piter_(aa_line_t) pxl(aa_line.domain()); + mln_niter_(neighb2d) ne(c4(), pxl); // FIXME: "pxl -> nbhoring edges" is explictly c4()... + for_all(pxl) + { + if (aa(pxl) != 0) + continue; + unsigned count = 0; + A na, na2; + for_all(ne) + if (aa_line.has(ne)) + { + mln_invariant(aa(ne) != 0); + if (count == 0) + na = aa(ne); // First attribute value encountered. + else + na2 = aa(ne); // Second (or more) attr value encountered. + ++count; + } + if (count == 2) + { + mln_invariant(na == na2); // Both attribute values have to be the same. + aa(pxl) = na; + } + } + + } + + debug::println("aa | line:", aa_line); + + + // Now dealing with all elements of basins: + + // - 2D-faces, i.e., original pixels; + // - 1D-faces, i.e., edges within regions; + // - 0D-faces, i.e., points within regions. + + mln_VAR(aa_basins, aa | (pw::value(wst_g_full) != pw::cst(0))); + + { + mln_piter_(aa_basins_t) p(aa_basins.domain()); + for_all(p) + aa(p) = a[wst_g_full(p)]; + } + + debug::println("aa | basins", aa_basins); + + } + + + + debug::println("aa:", aa); + + + +} // end of main + + + + +/* + + +for (i=0; i<V; i++) // Initialiser les heaps de fibonacci +{ + fh_setcmp(G->hp[i], cmp); // Chaque region a une heap de ses edges +} + +forall edges z that separates two regions v and w +{ + fh_insert(G->hp[v], (void *)(z)); // Ajouter les edges dans les heaps + fh_insert(G->hp[w], (void *)(z)); +} + +UFinit(G->V); // Initialiser l'union-find + +// Parcourir les regions par attribut croissant +for (j=0; j<V-1; j++) +{ + i = find(j); + + do + { // trouver l'edge minimum qui sorte de la region + e = fh_extractmin(G->hp[i]); + } + while ((UFfind(e->v, e->w)) && (e !=NULL)); + + if (e != NULL) + { // Normalement, e != NULL, sinon c'est un BIG pb!!! + int ui, uj, uk; + ui = find(e->v); + uj = find(e->w); + uk = UFunion(e->v,e->w); // Merger les regions + if (uk == ui) + { // et merger les edges + G->hp[ui] = fh_union(G->hp[ui], G->hp[uj]); + } + else + { + G->hp[uj] = fh_union(G->hp[uj], G->hp[ui]); + } + mst[k] = *e; // Garder l'arete + SaliencyWeight[k] = attribut[ui];// Poids dans la nouvelle hierarchie + OldWeight[k] = e->weight; // Poids dans l'ancienne hierarchie (inutile) + k++; + } + + // Calcul de la sortie + Pour toutes les edges separantes z=(x,y) + { + S[z] = max {SaliencyWeight[k] | sur le chemin reliant x a y dans mst} + } +} + + + */ Index: laurent/ismm2009.hh --- laurent/ismm2009.hh (revision 3162) +++ laurent/ismm2009.hh (working copy) @@ -49,6 +49,12 @@ return p.row() % 2 && p.col() % 2; } + inline + bool is_not_point(const point2d& p) + { + return ! is_point(p); + } + // Neighborhoods. @@ -268,19 +274,6 @@ } - // Find root. - - template <typename L> - inline - L find_root(std::vector<L>& par, L l) - { - if (par[l] == l) - return l; - else - return par[l] = find_root(par, par[l]); - } - - // Display. @@ -334,13 +327,13 @@ mln_VAR( wst_g, morpho::meyer_wst(g, e2e, n_basins) ); - // Test the consistency with regional minima. - { - L n_regmins; - mln_VAR( regmin_g, - labeling::regional_minima(g, e2e, n_regmins) ); - mln_invariant(n_basins == n_regmins); - } +// // Test the consistency with regional minima. +// { +// L n_regmins; +// mln_VAR( regmin_g, +// labeling::regional_minima(g, e2e, n_regmins) ); +// mln_invariant(n_basins == n_regmins); +// } return wst_g; }

1 0

3162: Small fixes.
by Guillaume Lazzara 16 Jan '09

16 Jan '09

* headers.mk: add new headers to distribution. * mln/core/internal/site_relative_iterator_base.hh: avoid a warning. * mln/labeling/blobs.hh, * mln/draw/line.hh: use trace::warning. * mln/trace/warning.hh: always print warnings. * mln/fun/x2x/rotation.hh: fix a precondition. * mln/morpho/tree/data.hh: add a missing method. * tests/level/median.cc: add a missing include. * tests/unit_test/Makefile.am, * tests/unit_test/mln_accu_transform_line.cc, * tests/unit_test/mln_convert_impl_from_double_to_value.cc, * tests/unit_test/mln_opt_value.cc: add missing unit_tests. * tests/unit_test/mln_histo_data.cc: rename as... * tests/unit_test/mln_histo_array.cc: ... this. --- milena/ChangeLog | 27 ++++++++++++++++++++ milena/headers.mk | 5 +++- .../core/internal/site_relative_iterator_base.hh | 2 +- milena/mln/draw/line.hh | 13 +++++---- milena/mln/fun/x2x/rotation.hh | 21 +++++++++------ milena/mln/labeling/blobs.hh | 6 ++-- milena/mln/morpho/tree/data.hh | 7 +++++ milena/mln/trace/warning.hh | 11 ++++--- milena/tests/level/median.cc | 1 + milena/tests/unit_test/Makefile.am | 10 ++++++- ...ln_histo_data.cc => mln_accu_transform_line.cc} | 6 ++-- .../mln_convert_impl_from_double_to_value.cc | 11 ++++++++ .../{mln_histo_data.cc => mln_histo_array.cc} | 6 ++-- .../{mln_histo_data.cc => mln_opt_value.cc} | 6 ++-- 14 files changed, 96 insertions(+), 36 deletions(-) copy milena/tests/unit_test/{mln_histo_data.cc => mln_accu_transform_line.cc} (55%) create mode 100644 milena/tests/unit_test/mln_convert_impl_from_double_to_value.cc copy milena/tests/unit_test/{mln_histo_data.cc => mln_histo_array.cc} (60%) rename milena/tests/unit_test/{mln_histo_data.cc => mln_opt_value.cc} (61%) diff --git a/milena/ChangeLog b/milena/ChangeLog index 583f2d7..5dee5c1 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,32 @@ 2009-01-16 Guillaume Lazzara <z(a)lrde.epita.fr> + Small fixes. + + * headers.mk: add new headers to distribution. + + * mln/core/internal/site_relative_iterator_base.hh: avoid a warning. + + * mln/labeling/blobs.hh, + * mln/draw/line.hh: use trace::warning. + + * mln/trace/warning.hh: always print warnings. + + * mln/fun/x2x/rotation.hh: fix a precondition. + + * mln/morpho/tree/data.hh: add a missing method. + + * tests/level/median.cc: add a missing include. + + * tests/unit_test/Makefile.am, + * tests/unit_test/mln_accu_transform_line.cc, + * tests/unit_test/mln_convert_impl_from_double_to_value.cc, + * tests/unit_test/mln_opt_value.cc: add missing unit_tests. + + * tests/unit_test/mln_histo_data.cc: rename as... + * tests/unit_test/mln_histo_array.cc: ... this. + +2009-01-16 Guillaume Lazzara <z(a)lrde.epita.fr> + Add soft heap implementation. * headers.mk: add new headers to distribution. diff --git a/milena/headers.mk b/milena/headers.mk index 9380631..a526da5 100644 --- a/milena/headers.mk +++ b/milena/headers.mk @@ -3,7 +3,7 @@ nobase_include_HEADERS = \ mln/version.hh \ mln/histo/all.hh \ -mln/histo/data.hh \ +mln/histo/array.hh \ mln/histo/compute.hh \ mln/histo/compute.spe.hh \ mln/histo/essential.hh \ @@ -263,6 +263,7 @@ mln/convert/to_window.hh \ mln/convert/from_to.hh \ mln/convert/impl/from_value_to_value.hh \ mln/convert/impl/all.hh \ +mln/convert/impl/from_double_to_value.hh \ mln/convert/impl/from_image_to_site_set.hh \ mln/convert/impl/from_float_to_value.hh \ mln/convert/impl/from_int_to_value.hh \ @@ -372,6 +373,7 @@ mln/accu/count_adjacent_vertices.hh \ mln/accu/convolve.hh \ mln/accu/v.hh \ mln/accu/rank_bool.hh \ +mln/accu/transform_line.hh \ mln/accu/min.hh \ mln/accu/transform_directional.hh \ mln/accu/compute.hh \ @@ -425,6 +427,7 @@ mln/math/essential.hh \ mln/math/acos.hh \ mln/math/round_sat.hh \ mln/opt/at.hh \ +mln/opt/value.hh \ mln/binarization/includes.hh \ mln/binarization/all.hh \ mln/binarization/binarization.hh \ diff --git a/milena/mln/core/internal/site_relative_iterator_base.hh b/milena/mln/core/internal/site_relative_iterator_base.hh index 78e42d3..d186b87 100644 --- a/milena/mln/core/internal/site_relative_iterator_base.hh +++ b/milena/mln/core/internal/site_relative_iterator_base.hh @@ -238,7 +238,7 @@ namespace mln template <typename P> inline void - site_relative_iterator_base<S,E>::center_at_(const P& c) + site_relative_iterator_base<S,E>::center_at_(const P&) { // Default is no-op, meaning "no extra code". } diff --git a/milena/mln/draw/line.hh b/milena/mln/draw/line.hh index 096b041..09d8003 100644 --- a/milena/mln/draw/line.hh +++ b/milena/mln/draw/line.hh @@ -1,4 +1,5 @@ -// Copyright (C) 2007, 2008 EPITA Research and Development Laboratory +// Copyright (C) 2007, 2008, 2009 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 @@ -28,10 +29,9 @@ #ifndef MLN_DRAW_LINE_HH # define MLN_DRAW_LINE_HH -/*! \file mln/draw/line.hh - * - * \brief Draw a line in an image. - */ +/// \file mln/draw/line.hh +/// +/// Draw a line in an image. # include <mln/core/concept/image.hh> # include <mln/core/site_set/p_line2d.hh> @@ -76,7 +76,8 @@ namespace mln { I& ima = exact(ima_); mln_precondition(ima.is_valid()); - // if (! ima.has(beg) || ! ima.has(end)) trace::warning("out"); + if (! ima.has(beg) || ! ima.has(end)) + trace::warning("Begin or end site is not part of the given image."); data::paste(pw::cst(v) | p_line2d(beg, end), safe(ima).rw()); } diff --git a/milena/mln/fun/x2x/rotation.hh b/milena/mln/fun/x2x/rotation.hh index ff6bbcf..039fe66 100644 --- a/milena/mln/fun/x2x/rotation.hh +++ b/milena/mln/fun/x2x/rotation.hh @@ -1,5 +1,5 @@ -// Copyright (C) 2007, 2008 EPITA Research and Development Laboratory -// (LRDE) +// Copyright (C) 2007, 2008, 2009 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 @@ -53,29 +53,32 @@ namespace mln { template < unsigned n, typename C > algebra::h_mat<n, C> - get_rot_h_mat(const float alpha_, const algebra::vec<3,C>& axis_) + get_rot_h_mat(const float alpha_, const algebra::vec<n,C>& axis_) { - assert(!"get_h_mat : n not implemented"); + mln_assertion(!"get_h_mat : n not implemented"); } template <typename C > algebra::h_mat<3, C> get_rot_h_mat(const float alpha_, const algebra::vec<3,C>& axis_) { - algebra::h_mat<3, C> m_; + //test axis is valid + typedef algebra::vec<3,C> vec_t; + //FIXME: cannot use '!=' operator. + mln_precondition(!(axis_ == vec_t(literal::zero))); const float cos_a = cos(alpha_); const float sin_a = sin(alpha_); const float u = axis_[0]; const float v = axis_[1]; const float w = axis_[2]; - //test axis is valid - assert(u != 0 && v != 0 && w != 0); const float u2 = u * u; const float v2 = v * v; const float w2 = w * w; const float uvw2 = u2 + v2 + w2; + algebra::h_mat<3, C> m_; + m_(0,0) = (u2 + (v2 + w2) * cos_a) / uvw2; m_(0,1) = (u*v * (1 - cos_a) - u * std::sqrt(uvw2) * sin_a) / uvw2; m_(0,2) = (u*w * (1 - cos_a) + v * std::sqrt(uvw2) * sin_a) / uvw2; @@ -103,11 +106,11 @@ namespace mln algebra::h_mat<2, C> get_rot_h_mat(const float alpha_, const algebra::vec<2,C>&) { - algebra::h_mat<2, C> m_; - const float cos_a = cos(alpha_); const float sin_a = sin(alpha_); + algebra::h_mat<2, C> m_; + m_(0,0) = cos_a; m_(0,1) = -sin_a; m_(0,2) = 0; m_(1,0) = sin_a; m_(1,1) = cos_a; m_(1,2) = 0; m_(2,0) = 0; m_(2,1) = 0; m_(2,2) = 1; diff --git a/milena/mln/labeling/blobs.hh b/milena/mln/labeling/blobs.hh index 789eec4..dc2dba0 100644 --- a/milena/mln/labeling/blobs.hh +++ b/milena/mln/labeling/blobs.hh @@ -1,5 +1,5 @@ -// Copyright (C) 2007, 2008 EPITA Research and Development Laboratory -// (LRDE) +// Copyright (C) 2007, 2008, 2009 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 @@ -101,7 +101,7 @@ namespace mln // Label this point component. if (nlabels == mln_max(L)) { - std::cerr << "FIXME: labeling aborted!" << std::endl; + trace::warning("labeling aborted!"); return output; } ++nlabels; diff --git a/milena/mln/morpho/tree/data.hh b/milena/mln/morpho/tree/data.hh index c7c7f86..d7871c1 100644 --- a/milena/mln/morpho/tree/data.hh +++ b/milena/mln/morpho/tree/data.hh @@ -105,6 +105,13 @@ namespace mln return parent_(p) == p || f_(parent_(p)) != f_(p); } + bool is_a_non_root_node(const mln_psite(I)& p) const + { + mln_precondition(is_valid()); + mln_precondition(parent_.domain().has(p)); + return f_(parent_(p)) != f_(p); + } + /// \} diff --git a/milena/mln/trace/warning.hh b/milena/mln/trace/warning.hh index c130da6..5e3ba29 100644 --- a/milena/mln/trace/warning.hh +++ b/milena/mln/trace/warning.hh @@ -1,4 +1,5 @@ -// Copyright (C) 2008 EPITA Research and Development Laboratory +// Copyright (C) 2008, 2009 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 @@ -29,6 +30,7 @@ # define MLN_TRACE_WARNING_HH /// \file mln/trace/warning.hh +/// /// Display warning message in trace output. # include <string> @@ -50,10 +52,9 @@ namespace mln inline void warning(const std::string& message) { - if (!quiet) - std::cout << "Warning: " - << message - << std::endl; + std::cout << "Warning: " + << message + << std::endl; } # endif // ! MLN_INCLUDE_ONLY diff --git a/milena/tests/level/median.cc b/milena/tests/level/median.cc index 0d749e2..40ca061 100644 --- a/milena/tests/level/median.cc +++ b/milena/tests/level/median.cc @@ -31,6 +31,7 @@ /// Test on mln::level::median. #include <mln/core/image/image2d.hh> +#include <mln/level/compare.hh> #include <mln/win/rectangle2d.hh> #include <mln/io/pgm/load.hh> diff --git a/milena/tests/unit_test/Makefile.am b/milena/tests/unit_test/Makefile.am index 2eef268..1d3d7cd 100644 --- a/milena/tests/unit_test/Makefile.am +++ b/milena/tests/unit_test/Makefile.am @@ -4,7 +4,7 @@ include $(top_srcdir)/milena/tests/tests.mk check_PROGRAMS = \ mln_histo_all \ -mln_histo_data \ +mln_histo_array \ mln_histo_compute \ mln_histo_essential \ mln_geom_pmin_pmax \ @@ -260,6 +260,7 @@ mln_convert_to_window \ mln_convert_from_to \ mln_convert_impl_from_value_to_value \ mln_convert_impl_all \ +mln_convert_impl_from_double_to_value \ mln_convert_impl_from_image_to_site_set \ mln_convert_impl_from_float_to_value \ mln_convert_impl_from_int_to_value \ @@ -368,6 +369,7 @@ mln_accu_count_adjacent_vertices \ mln_accu_convolve \ mln_accu_v \ mln_accu_rank_bool \ +mln_accu_transform_line \ mln_accu_min \ mln_accu_transform_directional \ mln_accu_compute \ @@ -419,6 +421,7 @@ mln_math_essential \ mln_math_acos \ mln_math_round_sat \ mln_opt_at \ +mln_opt_value \ mln_binarization_includes \ mln_binarization_all \ mln_binarization_binarization \ @@ -1015,7 +1018,7 @@ mln_subsampling_subsampling \ mln_subsampling_essential mln_histo_all_SOURCES = mln_histo_all.cc -mln_histo_data_SOURCES = mln_histo_data.cc +mln_histo_array_SOURCES = mln_histo_array.cc mln_histo_compute_SOURCES = mln_histo_compute.cc mln_histo_essential_SOURCES = mln_histo_essential.cc mln_geom_pmin_pmax_SOURCES = mln_geom_pmin_pmax.cc @@ -1271,6 +1274,7 @@ mln_convert_to_window_SOURCES = mln_convert_to_window.cc mln_convert_from_to_SOURCES = mln_convert_from_to.cc mln_convert_impl_from_value_to_value_SOURCES = mln_convert_impl_from_value_to_value.cc mln_convert_impl_all_SOURCES = mln_convert_impl_all.cc +mln_convert_impl_from_double_to_value_SOURCES = mln_convert_impl_from_double_to_value.cc mln_convert_impl_from_image_to_site_set_SOURCES = mln_convert_impl_from_image_to_site_set.cc mln_convert_impl_from_float_to_value_SOURCES = mln_convert_impl_from_float_to_value.cc mln_convert_impl_from_int_to_value_SOURCES = mln_convert_impl_from_int_to_value.cc @@ -1379,6 +1383,7 @@ mln_accu_count_adjacent_vertices_SOURCES = mln_accu_count_adjacent_vertices.cc mln_accu_convolve_SOURCES = mln_accu_convolve.cc mln_accu_v_SOURCES = mln_accu_v.cc mln_accu_rank_bool_SOURCES = mln_accu_rank_bool.cc +mln_accu_transform_line_SOURCES = mln_accu_transform_line.cc mln_accu_min_SOURCES = mln_accu_min.cc mln_accu_transform_directional_SOURCES = mln_accu_transform_directional.cc mln_accu_compute_SOURCES = mln_accu_compute.cc @@ -1430,6 +1435,7 @@ mln_math_essential_SOURCES = mln_math_essential.cc mln_math_acos_SOURCES = mln_math_acos.cc mln_math_round_sat_SOURCES = mln_math_round_sat.cc mln_opt_at_SOURCES = mln_opt_at.cc +mln_opt_value_SOURCES = mln_opt_value.cc mln_binarization_includes_SOURCES = mln_binarization_includes.cc mln_binarization_all_SOURCES = mln_binarization_all.cc mln_binarization_binarization_SOURCES = mln_binarization_binarization.cc diff --git a/milena/tests/unit_test/mln_histo_data.cc b/milena/tests/unit_test/mln_accu_transform_line.cc similarity index 55% copy from milena/tests/unit_test/mln_histo_data.cc copy to milena/tests/unit_test/mln_accu_transform_line.cc index 1921ad8..acabfd1 100644 --- a/milena/tests/unit_test/mln_histo_data.cc +++ b/milena/tests/unit_test/mln_accu_transform_line.cc @@ -1,9 +1,9 @@ -// Unit test for mln/histo/data.hh. +// Unit test for mln/accu/transform_line.hh. // Generated by ./build_unit_test.sh, do not modify. // Include the file twice, so we detect missing inclusion guards. -#include <mln/histo/data.hh> -#include <mln/histo/data.hh> +#include <mln/accu/transform_line.hh> +#include <mln/accu/transform_line.hh> int main() { diff --git a/milena/tests/unit_test/mln_convert_impl_from_double_to_value.cc b/milena/tests/unit_test/mln_convert_impl_from_double_to_value.cc new file mode 100644 index 0000000..c1270ba --- /dev/null +++ b/milena/tests/unit_test/mln_convert_impl_from_double_to_value.cc @@ -0,0 +1,11 @@ +// Unit test for mln/convert/impl/from_double_to_value.hh. +// Generated by ./build_unit_test.sh, do not modify. + +// Include the file twice, so we detect missing inclusion guards. +#include <mln/convert/impl/from_double_to_value.hh> +#include <mln/convert/impl/from_double_to_value.hh> + +int main() +{ + // Nothing. +} diff --git a/milena/tests/unit_test/mln_histo_data.cc b/milena/tests/unit_test/mln_histo_array.cc similarity index 60% copy from milena/tests/unit_test/mln_histo_data.cc copy to milena/tests/unit_test/mln_histo_array.cc index 1921ad8..31baac6 100644 --- a/milena/tests/unit_test/mln_histo_data.cc +++ b/milena/tests/unit_test/mln_histo_array.cc @@ -1,9 +1,9 @@ -// Unit test for mln/histo/data.hh. +// Unit test for mln/histo/array.hh. // Generated by ./build_unit_test.sh, do not modify. // Include the file twice, so we detect missing inclusion guards. -#include <mln/histo/data.hh> -#include <mln/histo/data.hh> +#include <mln/histo/array.hh> +#include <mln/histo/array.hh> int main() { diff --git a/milena/tests/unit_test/mln_histo_data.cc b/milena/tests/unit_test/mln_opt_value.cc similarity index 61% rename from milena/tests/unit_test/mln_histo_data.cc rename to milena/tests/unit_test/mln_opt_value.cc index 1921ad8..48469ea 100644 --- a/milena/tests/unit_test/mln_histo_data.cc +++ b/milena/tests/unit_test/mln_opt_value.cc @@ -1,9 +1,9 @@ -// Unit test for mln/histo/data.hh. +// Unit test for mln/opt/value.hh. // Generated by ./build_unit_test.sh, do not modify. // Include the file twice, so we detect missing inclusion guards. -#include <mln/histo/data.hh> -#include <mln/histo/data.hh> +#include <mln/opt/value.hh> +#include <mln/opt/value.hh> int main() { -- 1.5.6.5

1 0

3161: Add soft heap implementation.
by Guillaume Lazzara 16 Jan '09

16 Jan '09

* headers.mk: add new headers to distribution. * mln/util/soft_heap.hh: new file. New implementation. New members may be added later. * mln/util/tracked_ptr.hh: Fix an issue while assigning a null tracked_ptr to another tracked_ptr. An internal pointer was not initialized but dereferenced and led to a segfault. * tests/unit_test/Makefile.am, * tests/unit_test/mln_util_soft_heap.cc: add a new unit test. * tests/util/Makefile.am, * tests/util/soft_heap.cc: add a new test. --- milena/ChangeLog | 19 + milena/headers.mk | 1 + milena/mln/util/soft_heap.hh | 1137 ++++++++++++++++++++++++++ milena/mln/util/tracked_ptr.hh | 46 +- milena/tests/unit_test/Makefile.am | 2 + milena/tests/unit_test/mln_util_soft_heap.cc | 11 + milena/tests/util/Makefile.am | 2 + milena/tests/util/soft_heap.cc | 79 ++ 8 files changed, 1273 insertions(+), 24 deletions(-) create mode 100644 milena/mln/util/soft_heap.hh create mode 100644 milena/tests/unit_test/mln_util_soft_heap.cc create mode 100644 milena/tests/util/soft_heap.cc diff --git a/milena/ChangeLog b/milena/ChangeLog index 0166541..583f2d7 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,24 @@ 2009-01-16 Guillaume Lazzara <z(a)lrde.epita.fr> + Add soft heap implementation. + + * headers.mk: add new headers to distribution. + + * mln/util/soft_heap.hh: new file. New implementation. New members may + be added later. + + * mln/util/tracked_ptr.hh: Fix an issue while assigning a null + tracked_ptr to another tracked_ptr. An internal pointer was not + initialized but dereferenced and led to a segfault. + + * tests/unit_test/Makefile.am, + * tests/unit_test/mln_util_soft_heap.cc: add a new unit test. + + * tests/util/Makefile.am, + * tests/util/soft_heap.cc: add a new test. + +2009-01-16 Guillaume Lazzara <z(a)lrde.epita.fr> + Rename histo::data as histo::array. * mln/histo/data.hh: rename file as... diff --git a/milena/headers.mk b/milena/headers.mk index fc78f90..9380631 100644 --- a/milena/headers.mk +++ b/milena/headers.mk @@ -64,6 +64,7 @@ mln/registration/icp.hh \ mln/util/graph.hh \ mln/util/max.hh \ mln/util/lazy_set.hh \ +mln/util/soft_heap.hh \ mln/util/set.hh \ mln/util/tree_to_image.hh \ mln/util/lemmings.hh \ diff --git a/milena/mln/util/soft_heap.hh b/milena/mln/util/soft_heap.hh new file mode 100644 index 0000000..144110f --- /dev/null +++ b/milena/mln/util/soft_heap.hh @@ -0,0 +1,1137 @@ +// Copyright (C) 2009 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. + +#ifndef MLN_UTIL_FAST_HEAP_HH +# define MLN_UTIL_FAST_HEAP_HH + +/// \file mln/util/soft_heap.hh +/// +/// Define a generic soft heap. +/* + This implementation is not an usual heap, it allows to set an error + rate so that some nodes may be "corrupted". A "corrupted node" means + that its correct order is not totally preserved for performance reasons. + Of course, it will have an impact on the returned values. + As a result, be ware of not using this data structure if the element + order is relevant for to you. + + A corruption threshold can be passed to the constructor. + This threshold means that if nodes have a rank higher than this + threshold they can be "corrupted" and therefore their rank can be + reduced. Tuning this threshold may have an impact on the structure + entropy thus on the returned values order. It may also have an impact + on the performance. + + + More implementation details are available in: + "The soft heap: an approximate priority queue with optimal error rate", + Bernard Chazelle, JACM, 2000. + + URL: http://www.cs.princeton.edu/~chazelle/pubs/sheap.pdf + + \todo try not to use tracked_ptr. + \todo try not to call T::plus_infty(). This is not generic enough but + it works with mln::point<>. +*/ + +# include <mln/core/concept/object.hh> +# include <mln/trait/value_.hh> +# include <mln/util/tracked_ptr.hh> +# include <climits> +# include <stack> + + + +namespace mln +{ + + namespace util + { + + + /// Element of an item list. Store the data (key) used in soft_heap. + template <typename T> + struct ilcell + { + typedef util::tracked_ptr< ilcell<T> > ilcell_t; + + ilcell(); + ilcell(const T& key, ilcell_t next = 0); + + ilcell_t next() const; + const T& key() const; + + void set_next(ilcell_t next); + void set_key(const T& key); + + private: + T key_; + ilcell_t next_; + }; + + + /// Meta-data of an element in the heap. + template <typename T, typename R> + class node + { + + typedef util::tracked_ptr< ilcell<T> > ilcell_t; + + public: + node(); + node(const T& ckey, const R& rank, + node<T,R> *next = 0, node<T,R> *child = 0, + ilcell_t il = 0, ilcell_t il_tail = 0); + ~node(); + + const T& ckey() const; + const R& rank() const; + node<T,R> * next() const; + node<T,R> * child() const; + ilcell_t il() const; + ilcell_t il_tail() const; + + void set_il(ilcell_t il); + void set_il_tail(ilcell_t il_tail); + void set_ckey(const T& ckey); + void set_rank(const R& rank); + void set_next(node<T,R> * next); + void set_child(node<T,R> *child); + + private: + T ckey_; + R rank_; + node<T,R> *next_; + node<T,R> *child_; + ilcell_t il_; + ilcell_t il_tail_; + + }; + + + + + /// Top structure of the soft heap. + template <typename T, typename R> + class head + { + public: + + head(); + head(const R& rank, node<T,R> *queue = 0, head<T,R> *next = 0, + head<T,R> *prev = 0, head<T,R> *suffix_min = 0); + ~head(); + + node<T,R> *queue() const; + head<T,R> *next() const; + head<T,R> *prev() const; + head<T,R> *suffix_min() const; + const R& rank() const; + + void set_queue(node<T,R> *queue); + void set_next(head<T,R> *next); + void set_prev(head<T,R> *prev); + void set_suffix_min(head<T,R> *suffix_min); + void set_rank(const R& rank); + + private: + node<T,R> *queue_; + head<T,R> *next_; + head<T,R> *prev_; + head<T,R> *suffix_min_; + R rank_; + + }; + + + + + /// T key, the data to store in the heap. For instance a point 2d. + /// R rank, for instance int_u8 + template <typename T, typename R> + class soft_heap : public Object< soft_heap<T,R> > + { + typedef util::tracked_ptr< ilcell<T> > ilcell_t; + + public: + + /// Element associated type. + typedef T element; + + + /// Default constructor. + /// A corruption threshold \p r can be specified. + /// This threshold means that if nodes have a rank higher than this + /// threshold they can be "corrupted" and therefore their rank can be + /// reduced. + /// + // FIXME: is it a correct value for r? + soft_heap(unsigned r = 20); + /// Destructor + ~soft_heap(); + + + /// Add a new element \p element. + void push(const T& element); + /// Insert a site \p p (equivalent as 'push'). + void insert(const T& p); + + /// Merge \p sh with this heap. + /// Be ware that after this call, \p sh will be empty. This heap will + /// hold the elements which were part of \p sh. + void take(soft_heap<T,R>& sh); + + /// Returns the element with the lowest priority and remove it from the + /// heap. + T pop_front(); + + /// Return true if there is at least one element. + bool is_valid() const; + + /// Return true if there is at least one element. + bool is_empty() const; + + /// Return the number of element in the heap. + int nelements() const; + + /// Clear the heap. + void clear(); + + /// Internal use only. Return a pointer to the first header struct + /// of this heap. + head<T,R> *head_hook_() const; + + /// Reset the heap to an empty heap. Do *NOT* delete element which may + /// have been inserted. + /// \sa take + void soft_clear_(); + + + + private: + /// Merge a node \p q to this heap. + /// \sa take + void meld(node<T,R> *q); + + /// Update suffix_min pointer according to the new values inserted in + /// the item lists. + void fix_minlist(head<T,R> *h); + + /// After the deletion of an element, this function move/update item lists + /// in the proper nodes. + node<T,R> * sift(node<T,R> *v); + + /// Delete all element in a C-list from \p begin to \p end. \p end is not + /// deleted. + /// The type L must provide 'L *next()'. + /// + /// next + /// + /// |x| -> |x| -> |x| + /// + template <typename L> + void clear_list(L *begin, L *end = 0); + + /// Delete all element in a 2 dimensional C-list. + /// The type L must provide 'L *next()' and 'L *child()'. + /// + /// next + /// + /// |x| -> |x| -> |x| + /// | + /// child v + /// |x| -> |x| + /// + template <typename L> + void deep_clear_list(L *n); + + /// Print which node is part of which head. + void debug_head_list() const; + + /// Display the heap, preserving the hierarchy of the elements. + void println_() const; + + head<T,R> *header_; + head<T,R> *tail_; + + /// Corruption threshold. Must be set once through the constructor. + unsigned r_; + + int nelements_; + }; + + + +# ifndef MLN_INCLUDE_ONLY + + + /*-------` + | ilcell | + \-------*/ + + + template <typename T> + inline + ilcell<T>::ilcell() + { + } + + + template <typename T> + inline + ilcell<T>::ilcell(const T& key, ilcell_t next) + : key_(key), next_(next) + { + } + + + template <typename T> + inline + typename ilcell<T>::ilcell_t + ilcell<T>::next() const + { + return next_; + } + + + template <typename T> + inline + const T& + ilcell<T>::key() const + { + return key_; + } + + + template <typename T> + inline + void + ilcell<T>::set_next(ilcell_t next) + { + next_ = next; + } + + + template <typename T> + inline + void + ilcell<T>::set_key(const T& key) + { + key_ = key; + } + + + + + /*-----` + | node | + \-----*/ + + + template <typename T, typename R> + inline + node<T,R>::node() + : il_(0), il_tail_(0) + { + } + + + template <typename T, typename R> + inline + node<T,R>::node(const T& ckey, const R& rank, + node<T,R> *next, node<T,R> *child, + ilcell_t il, ilcell_t il_tail) + : ckey_(ckey), rank_(rank), next_(next), child_(child), + il_(il), il_tail_(il_tail) + { + } + + + template <typename T, typename R> + inline + node<T,R>::~node() + { + } + + + template <typename T, typename R> + inline + const T& + node<T,R>::ckey() const + { + return ckey_; + } + + + template <typename T, typename R> + inline + const R& + node<T,R>::rank() const + { + return rank_; + } + + + template <typename T, typename R> + inline + node<T,R> * + node<T,R>::next() const + { + return next_; + } + + + template <typename T, typename R> + inline + node<T,R> * + node<T,R>::child() const + { + return child_; + } + + + template <typename T, typename R> + inline + typename node<T,R>::ilcell_t + node<T,R>::il() const + { + return il_; + } + + + template <typename T, typename R> + inline + typename node<T,R>::ilcell_t + node<T,R>::il_tail() const + { + return il_tail_; + } + + + template <typename T, typename R> + inline + void + node<T,R>::set_il(ilcell_t il) + { + il_ = il; + } + + + template <typename T, typename R> + inline + void + node<T,R>::set_il_tail(ilcell_t il_tail) + { + il_tail_ = il_tail; + } + + + template <typename T, typename R> + inline + void + node<T,R>::set_ckey(const T& ckey) + { + ckey_ = ckey; + } + + + template <typename T, typename R> + inline + void + node<T,R>::set_rank(const R& rank) + { + rank_ = rank; + } + + + template <typename T, typename R> + inline + void + node<T,R>::set_next(node<T,R> * next) + { + next_ = next; + } + + + template <typename T, typename R> + inline + void + node<T,R>::set_child(node<T,R> *child) + { + child_ = child; + } + + + + + /*-----` + | head | + \-----*/ + + + template <typename T, typename R> + inline + head<T,R>::head() + { + } + + + template <typename T, typename R> + inline + head<T,R>::head(const R& rank, node<T,R> *queue, head<T,R> *next, + head<T,R> *prev, head<T,R> *suffix_min) + : queue_(queue), next_(next), prev_(prev), + suffix_min_(suffix_min), rank_(rank) + { + } + + + template <typename T, typename R> + inline + head<T,R>::~head() + { + } + + + template <typename T, typename R> + inline + node<T,R> * + head<T,R>::queue() const + { + return queue_; + } + + + template <typename T, typename R> + inline + head<T,R> * + head<T,R>::next() const + { + return next_; + } + + + template <typename T, typename R> + inline + head<T,R> * + head<T,R>::prev() const + { + return prev_; + } + + + template <typename T, typename R> + inline + head<T,R> * + head<T,R>::suffix_min() const + { + return suffix_min_; + } + + + template <typename T, typename R> + inline + const R& + head<T,R>::rank() const + { + return rank_; + } + + + template <typename T, typename R> + inline + void + head<T,R>::set_queue(node<T,R> *queue) + { + queue_ = queue; + } + + + template <typename T, typename R> + inline + void + head<T,R>::set_next(head<T,R> *next) + { + next_ = next; + } + + + template <typename T, typename R> + inline + void + head<T,R>::set_prev(head<T,R> *prev) + { + prev_ = prev; + } + + + template <typename T, typename R> + inline + void + head<T,R>::set_suffix_min(head<T,R> *suffix_min) + { + suffix_min_ = suffix_min; + } + + + template <typename T, typename R> + inline + void + head<T,R>::set_rank(const R& rank) + { + rank_ = rank; + } + + + + + /*------------` + | soft_heap | + \------------*/ + + + template <typename T, typename R> + inline + soft_heap<T,R>::soft_heap(unsigned r) + { + nelements_ = 0; + r_ = r; + header_ = new head<T,R>(mln_max(R)); + tail_ = new head<T,R>(mln_max(R), 0, 0, header_); + header_->set_next(tail_); + } + + + template <typename T, typename R> + inline + soft_heap<T,R>::~soft_heap() + { + head<T,R> *tmp = header_; + while (tmp != 0) + { + deep_clear_list(tmp->queue()); + tmp = tmp->next(); + } + clear_list(header_, tail_->next()); + } + + + template <typename T, typename R> + inline + void + soft_heap<T,R>::push(const T& element) + { + ilcell_t l(new ilcell<T>(element)); + node<T,R> *q = new node<T,R>(element, 0, 0, 0, l, l); + meld(q); + ++nelements_; + } + + + template <typename T, typename R> + inline + void + soft_heap<T,R>::take(soft_heap<T,R>& psh) + { + head<T,R> *head = psh.head_hook_(); + while (head != 0) + { + if (head->queue() != 0) + meld(head->queue()); + head = head->next(); + } + nelements_ += psh.nelements(); + psh.soft_clear_(); + } + + + template <typename T, typename R> + inline + T + soft_heap<T,R>::pop_front() + { + mln_precondition(is_valid()); + + node<T,R> *tmp; + + T min; + int childcount; + head<T,R> *h = header_->next()->suffix_min(); + while (h->queue()->il() == 0) + { + tmp = h->queue(); + childcount = 0; + while (tmp->next() != 0) + { + tmp = tmp->next(); + ++childcount; + } + + if (childcount < h->rank() / 2) + { + h->prev()->set_next(h->next()); + h->next()->set_prev(h->prev()); + fix_minlist(h->prev()); + tmp = h->queue(); + while (tmp->next() != 0) + { + meld(tmp->child()); + tmp = tmp->next(); + } + //FIXME: is it right? + deep_clear_list(h->queue()); + delete h; + } + else + { + h->set_queue(sift(h->queue())); + // FIXME: not generic enough. + if (h->queue()->ckey() == T::plus_infty()) + { + h->prev()->set_next(h->next()); + h->next()->set_prev(h->prev()); + + // Remove h and switch to h->prev. + head<T,R> *h_bak = h; + h = h->prev(); + deep_clear_list(h_bak->queue()); + delete h_bak; + } + fix_minlist(h); + } + h = header_->next()->suffix_min(); + } + + min = h->queue()->il()->key(); + + //Don't need to delete h->queue()->il(). May be used in another node. + h->queue()->set_il(h->queue()->il()->next()); + if (h->queue()->il() == 0) + h->queue()->set_il_tail(0); + + --nelements_; + return min; + } + + + template <typename T, typename R> + inline + bool + soft_heap<T,R>::is_valid() const + { + return header_ != 0 && tail_ != 0; + } + + + template <typename T, typename R> + inline + bool + soft_heap<T,R>::is_empty() const + { + return nelements_ == 0 ; + } + + + template <typename T, typename R> + inline + int + soft_heap<T,R>::nelements() const + { + return nelements_; + } + + + template <typename T, typename R> + inline + void + soft_heap<T,R>::clear() + { + if (header_->next() == tail_) + return; + + head<T,R> *tohead = header_->next(); + head<T,R> *prevtail = tail_->prev(); + prevtail->set_next(0); + tohead->set_prev(0); + header_->set_next(tail_); + tail_->set_prev(header_); + + head<T,R> *tmp = tohead; + while (tmp != 0) + { + deep_clear_list(tmp->queue()); + tmp = tmp->next(); + } + clear_list(tohead); + + nelements_ = 0; + } + + + template <typename T, typename R> + inline + void + soft_heap<T,R>::soft_clear_() + { + clear_list(header_->next(), tail_); + header_->set_next(tail_); + header_->set_prev(0); + tail_->set_next(0); + tail_->set_prev(header_); + nelements_ = 0; + } + + + template <typename T, typename R> + inline + head<T,R> * + soft_heap<T,R>::head_hook_() const + { + return header_; + } + + + template <typename T, typename R> + inline + void + soft_heap<T,R>::meld(node<T,R>* q) + { + head<T,R> *tohead = header_->next(); + while (q->rank() > tohead->rank()) + tohead = tohead->next(); + head<T,R> *prevhead = tohead->prev(); + + node<T,R> *top, *bottom; + + while (q->rank() == tohead->rank()) + { + if (util::ord_strict(q->ckey(), tohead->queue()->ckey())) + { + top = q; + bottom = tohead->queue(); + } + else + { + top = tohead->queue(); + bottom = q; + } + + q = new node<T,R>(top->ckey(), top->rank() + 1, top, bottom, + top->il(), top->il_tail()); + + tohead = tohead->next(); + } + + head<T,R> *h; + if (prevhead == tohead->prev()) + { + h = new head<T,R>(q->rank(), q, tohead, prevhead); + } + else + { + // Do not create a new head. + h = prevhead->next(); + + //Discard/delete everything between h and tohead. + head<T,R> *head_del = h->next(), *tmp_del; + while (head_del != tohead) + { + tmp_del = head_del->next(); + delete head_del; + head_del = tmp_del; + } + + } + h->set_queue(q); + h->set_rank(q->rank()); + h->set_prev(prevhead); + h->set_next(tohead); + + prevhead->set_next(h); + tohead->set_prev(h); + + fix_minlist(h); + } + + + template <typename T, typename R> + inline + void + soft_heap<T,R>::fix_minlist(head<T,R> *h) + { + head<T,R> *tmpmin; + if (h->next() == tail_) + tmpmin = h; + else + tmpmin = h->next()->suffix_min(); + + while (h != header_) + { + if (util::ord_strict(h->queue()->ckey(), tmpmin->queue()->ckey())) + tmpmin = h; + h->set_suffix_min(tmpmin); + h = h->prev(); + } + } + + + template <typename T, typename R> + inline + node<T,R> * + soft_heap<T,R>::sift(node<T,R> *v) + { + node<T,R> *tmp; + // We do not need to free the list since these object are also pointed by + // other structs. + v->set_il(0); + v->set_il_tail(0); + + if (v->next() == 0 && v->child() == 0) + { + //FIXME: not generic enough... + v->set_ckey(T::plus_infty()); + return v; + } + node<T,R> *v_next_bak = v->next(); + v->set_next(sift(v->next())); + + // Cleanup possibly removed nodes. + while (v_next_bak != v->next()) + { + deep_clear_list(v_next_bak->next()); + deep_clear_list(v_next_bak->child()); + tmp = v_next_bak->next(); + delete v_next_bak; + v_next_bak = tmp; + } + + if (util::ord_strict(v->child()->ckey(), v->next()->ckey())) + { + // Swap child and next. + tmp = v->child(); + v->set_child(v->next()); + v->set_next(tmp); + } + + v->set_il(v->next()->il()); + v->set_il_tail(v->next()->il_tail()); + v->set_ckey(v->next()->ckey()); + + if (v->rank() > r_ && (v->rank() % 2 == 1 + || v->child()->rank() < static_cast<unsigned>((v->rank() - 1)))) + { + node<T,R> *v_next_bak = v->next(); + v->set_next(sift(v->next())); + + // Cleanup possibly removed nodes. + while (v_next_bak != v->next()) + { + deep_clear_list(v_next_bak->next()); + deep_clear_list(v_next_bak->child()); + tmp = v_next_bak->next(); + delete v_next_bak; + v_next_bak = tmp; + } + + if (util::ord_strict(v->child()->ckey(), v->next()->ckey())) + { + tmp = v->child(); + v->set_child(v->next()); + v->set_next(tmp); + } + + if (v->next()->ckey() != T::plus_infty() && v->next()->il() != 0) + { + v->next()->il_tail()->set_next(v->il()); + v->set_il(v->next()->il()); + if (v->il_tail() == 0) + v->set_il_tail(v->next()->il_tail()); + v->set_ckey(v->next()->ckey()); + } + } + + if (v->child()->ckey() == T::plus_infty()) + { + if (v->next()->ckey() == T::plus_infty()) + { + // Clear child and next lists. + deep_clear_list(v->child()); + deep_clear_list(v->next()); + + v->set_child(0); + v->set_next(0); + } + else + { + node<T,R> *next_bak = v->next(); + + // There may be several children, we must clear all of them. + deep_clear_list(v->child()); + + v->set_child(v->next()->child()); + v->set_next(v->next()->next()); + delete next_bak; + } + } + + return v; + } + + + template <typename T, typename R> + template <typename L> + inline + void + soft_heap<T,R>::clear_list(L *begin, L *end) + { + L *tmp; + while (begin != end) + { + tmp = begin->next(); + delete begin; + begin = tmp; + } + } + + + template <typename T, typename R> + template <typename L> + inline + void + soft_heap<T,R>::deep_clear_list(L *n) + { + L *current_node, *tmp; + std::stack<L *> st; + + st.push(n); + while (!st.empty()) + { + current_node = st.top(); + st.pop(); + while (current_node != 0) + { + if (current_node->child() != 0) + st.push(current_node->child()); + + tmp = current_node->next(); + delete current_node; + + current_node = tmp; + } + } + } + + + template <typename T, typename R> + inline + void + soft_heap<T,R>::debug_head_list() const + { + head<T,R> *n = header_; + std::cout << "Head list = " << std::endl; + while (n != 0) + { + std::cout << n->id << "("; + + node<T,R> *current_node; + std::stack< node<T,R> *> st; + st.push(n->queue()); + while (!st.empty()) + { + current_node = st.top(); + st.pop(); + while (current_node != 0) + { + if (current_node->child() != 0) + st.push(current_node->child()); + std::cout << current_node->id << ","; + current_node = current_node->next(); + } + } + + std::cout << ") - "; + n = n->next(); + } + std::cout << std::endl; + } + + + template <typename T, typename R> + inline + void + soft_heap<T,R>::println_() const + { + + std::cout << "===============" << std::endl; + head<T,R> *head = header_; + while (head != 0) + { + std::cout << "<Head>" << std::endl; + node<T,R> *node = head->queue(), *child; + while (node != 0) + { + std::cout << " <node>" << std::endl; + + ilcell_t il(node->il()); + while (il != 0) + { + std::cout << il->item() << std::endl; + il = il->next(); + } + + child = node->child(); + while (child != 0) + { + std::cout << " <child>" << std::endl; + ilcell_t il(child->il()); + while (il != 0) + { + std::cout << il->item() << std::endl; + il = il->next(); + } + child = child->child(); + std::cout << " </child>" << std::endl; + } + node = node->next(); + + std::cout << " </node>" << std::endl; + } + std::cout << "</Head>" << std::endl; + head = head->next(); + } + } + + +# endif // ! MLN_INCLUDE_ONLY + + + } // end of namespace mln::util + +} // end of namespace mln + + +#endif // ! MLN_UTIL_FAST_HEAP_HH diff --git a/milena/mln/util/tracked_ptr.hh b/milena/mln/util/tracked_ptr.hh index 6d22451..0148880 100644 --- a/milena/mln/util/tracked_ptr.hh +++ b/milena/mln/util/tracked_ptr.hh @@ -1,4 +1,5 @@ -// Copyright (C) 2006 EPITA Research and Development Laboratory +// Copyright (C) 2006, 2009 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 @@ -28,12 +29,11 @@ #ifndef MLN_UTIL_TRACKED_PTR_HH # define MLN_UTIL_TRACKED_PTR_HH -/*! \file mln/util/tracked_ptr.hh - * - * \brief Definition of a smart pointer for shared data with tracking. - * - * \todo Split defs from decls. - */ +/// \file mln/util/tracked_ptr.hh +/// +/// Definition of a smart pointer for shared data with tracking. +/// +/// \todo Split defs from decls. # include <set> # include <iostream> @@ -63,16 +63,14 @@ namespace mln /// Negation (for arithmetical tests). bool operator !() const; - /*! \brief Mimics the behavior of op-> for a pointer in the const case. - ** - ** \invariant Pointer proxy exists. - */ + /// Mimics the behavior of op-> for a pointer in the const case. + /// + /// \invariant Pointer proxy exists. const T* operator->() const; - /*! \brief Mimics the behavior of op-> for a pointer in the mutable case. - ** - ** \invariant Pointer proxy exists. - */ + /// Mimics the behavior of op-> for a pointer in the mutable case. + /// + /// \invariant Pointer proxy exists. T* operator->(); /// Ctor. @@ -122,10 +120,9 @@ namespace mln } template <typename T> - /*! \brief Mimics the behavior of op-> for a pointer in the const case. - ** - ** \invariant Pointer proxy exists. - */ + /// Mimics the behavior of op-> for a pointer in the const case. + /// + /// \invariant Pointer proxy exists. inline const T* tracked_ptr<T>::operator->() const { @@ -135,10 +132,9 @@ namespace mln } template <typename T> - /*! \brief Mimics the behavior of op-> for a pointer in the mutable case. - ** - ** \invariant Pointer proxy exists. - */ + /// \brief Mimics the behavior of op-> for a pointer in the mutable case. + /// + /// \invariant Pointer proxy exists. inline T* tracked_ptr<T>::operator->() { @@ -201,7 +197,9 @@ namespace mln clean_(); ptr_ = rhs.ptr_; holders_ = rhs.holders_; - holders_->insert(this); + // If ptr == 0, holders_ == 0 so we cannot insert anything in it. + if (holders_ != 0) + holders_->insert(this); return *this; } diff --git a/milena/tests/unit_test/Makefile.am b/milena/tests/unit_test/Makefile.am index e055ee3..2eef268 100644 --- a/milena/tests/unit_test/Makefile.am +++ b/milena/tests/unit_test/Makefile.am @@ -64,6 +64,7 @@ mln_registration_icp \ mln_util_graph \ mln_util_max \ mln_util_lazy_set \ +mln_util_soft_heap \ mln_util_set \ mln_util_tree_to_image \ mln_util_lemmings \ @@ -1074,6 +1075,7 @@ mln_registration_icp_SOURCES = mln_registration_icp.cc mln_util_graph_SOURCES = mln_util_graph.cc mln_util_max_SOURCES = mln_util_max.cc mln_util_lazy_set_SOURCES = mln_util_lazy_set.cc +mln_util_soft_heap_SOURCES = mln_util_soft_heap.cc mln_util_set_SOURCES = mln_util_set.cc mln_util_tree_to_image_SOURCES = mln_util_tree_to_image.cc mln_util_lemmings_SOURCES = mln_util_lemmings.cc diff --git a/milena/tests/unit_test/mln_util_soft_heap.cc b/milena/tests/unit_test/mln_util_soft_heap.cc new file mode 100644 index 0000000..6a98bb0 --- /dev/null +++ b/milena/tests/unit_test/mln_util_soft_heap.cc @@ -0,0 +1,11 @@ +// Unit test for mln/util/soft_heap.hh. +// Generated by ./build_unit_test.sh, do not modify. + +// Include the file twice, so we detect missing inclusion guards. +#include <mln/util/soft_heap.hh> +#include <mln/util/soft_heap.hh> + +int main() +{ + // Nothing. +} diff --git a/milena/tests/util/Makefile.am b/milena/tests/util/Makefile.am index 0e9c5ab..1dcd5c5 100644 --- a/milena/tests/util/Makefile.am +++ b/milena/tests/util/Makefile.am @@ -12,6 +12,7 @@ check_PROGRAMS = \ lazy_set \ lemmings \ ord_pair \ + soft_heap \ tree \ tree_fast \ tree_fast_to_image \ @@ -26,6 +27,7 @@ graph_SOURCES = graph.cc lazy_set_SOURCES = lazy_set.cc lemmings_SOURCES = lemmings.cc ord_pair_SOURCES = ord_pair.cc +soft_heap_SOURCES = soft_heap.cc tree_SOURCES = tree.cc tree_fast_SOURCES = tree_fast.cc tree_fast_to_image_SOURCES = tree_to_image.cc diff --git a/milena/tests/util/soft_heap.cc b/milena/tests/util/soft_heap.cc new file mode 100644 index 0000000..d2bffcc --- /dev/null +++ b/milena/tests/util/soft_heap.cc @@ -0,0 +1,79 @@ +// Copyright (C) 2009 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 mln/core/site_set/p_soft_heap.hh +/// +/// Test for the soft heap (approximate priority queue). + +#include <mln/util/soft_heap.hh> +#include <mln/core/alias/point2d.hh> +#include <mln/value/int_u8.hh> + +using mln::point2d; + +point2d p[] = { point2d(4,5), point2d(3,4), point2d(3,2), + point2d(1,6), point2d(2,3), point2d(3,5), + point2d(2,4), point2d(7,2), point2d(9,6), + point2d(7,3) }; + +point2d res[] = { point2d(1,6), point2d(2,3), point2d(2,4), + point2d(3,2), point2d(3,4), point2d(3,5), + point2d(4,5), point2d(7,2), point2d(7,3), + point2d(9,6) }; + + +int main() +{ + using namespace mln; + + util::soft_heap<point2d, value::int_u8> fh; + + for (unsigned i = 0; i < 5; ++i) + fh.push(p[i]); + + util::soft_heap<point2d, value::int_u8> fh2; + for (unsigned i = 5; i < 10; ++i) + fh2.push(p[i]); + + // Merge fh in fh2. + fh2.take(fh); + + // fh2 now holds both its elements and fh's. + unsigned i = 0; + while (!fh2.is_empty()) + mln_assertion(fh2.pop_front() == res[i++]); + + // fh must be empty, fh2 now holds its elements. + mln_assertion(fh.is_empty()); + + for (unsigned i = 5; i < 10; ++i) + fh2.push(p[i]); + + fh2.clear(); + mln_assertion(fh2.nelements() == 0); + mln_assertion(fh2.is_empty()); +} -- 1.5.6.5

1 0

3160: Rename histo::data as histo::array.
by Guillaume Lazzara 16 Jan '09

16 Jan '09

* mln/histo/data.hh: rename file as... * mln/histo/array.hh: ... this. * mln/accu/histo.hh, * mln/convert/to_image.hh, * mln/histo/all.hh, * mln/histo/compute.hh, * mln/histo/compute.spe.hh, * mln/level/sort_psites.hh: update use of histo::array. --- milena/ChangeLog | 14 ++++++++++ milena/mln/accu/histo.hh | 6 ++-- milena/mln/convert/to_image.hh | 10 +++--- milena/mln/histo/all.hh | 12 ++++---- milena/mln/histo/{data.hh => array.hh} | 45 +++++++++++++++---------------- milena/mln/histo/compute.hh | 24 ++++++++-------- milena/mln/histo/compute.spe.hh | 8 +++--- milena/mln/level/sort_psites.hh | 4 +- 8 files changed, 68 insertions(+), 55 deletions(-) rename milena/mln/histo/{data.hh => array.hh} (80%) diff --git a/milena/ChangeLog b/milena/ChangeLog index 5ebfb83..0166541 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,3 +1,17 @@ +2009-01-16 Guillaume Lazzara <z(a)lrde.epita.fr> + + Rename histo::data as histo::array. + + * mln/histo/data.hh: rename file as... + * mln/histo/array.hh: ... this. + + * mln/accu/histo.hh, + * mln/convert/to_image.hh, + * mln/histo/all.hh, + * mln/histo/compute.hh, + * mln/histo/compute.spe.hh, + * mln/level/sort_psites.hh: update use of histo::array. + 2009-01-14 Alexandre Abraham <abraham(a)lrde.epita.fr> Fix for two ways functions. diff --git a/milena/mln/accu/histo.hh b/milena/mln/accu/histo.hh index 2bf9bdf..b66ed52 100644 --- a/milena/mln/accu/histo.hh +++ b/milena/mln/accu/histo.hh @@ -33,7 +33,7 @@ /// /// Define a generic histogram accumulator class. /// -/// \todo Use histo::data instead of std::vector! +/// \todo Use histo::array instead of std::vector! # include <vector> # include <algorithm> @@ -42,7 +42,7 @@ # include <mln/core/concept/meta_accumulator.hh> # include <mln/accu/internal/base.hh> # include <mln/value/set.hh> -# include <mln/histo/data.hh> +# include <mln/histo/array.hh> namespace mln @@ -91,7 +91,7 @@ namespace mln protected: - mln::histo::data<V> h_; + mln::histo::array<V> h_; unsigned sum_; }; diff --git a/milena/mln/convert/to_image.hh b/milena/mln/convert/to_image.hh index 809bb14..e71aeda 100644 --- a/milena/mln/convert/to_image.hh +++ b/milena/mln/convert/to_image.hh @@ -1,5 +1,5 @@ -// Copyright (C) 2007, 2008 EPITA Research and Development Laboratory -// (LRDE) +// Copyright (C) 2007, 2008, 2009 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 @@ -42,7 +42,7 @@ # include <mln/geom/bbox.hh> # include <mln/data/fill.hh> -# include <mln/histo/data.hh> +# include <mln/histo/array.hh> # include <mln/core/image/image1d.hh> # include <mln/core/image/image2d.hh> @@ -81,7 +81,7 @@ namespace mln /// Convert an histo \p h into an image1d<unsigned>. template <typename T> image1d<unsigned> - to_image(const histo::data<T>& h); + to_image(const histo::array<T>& h); @@ -138,7 +138,7 @@ namespace mln template <typename T> inline image1d<unsigned> - to_image(const histo::data<T>& h) + to_image(const histo::array<T>& h) { def::coord n = static_cast<def::coord>(h.vset().nvalues()), diff --git a/milena/mln/histo/all.hh b/milena/mln/histo/all.hh index 746df6b..225a9b5 100644 --- a/milena/mln/histo/all.hh +++ b/milena/mln/histo/all.hh @@ -1,4 +1,5 @@ -// Copyright (C) 2007 EPITA Research and Development Laboratory +// Copyright (C) 2007, 2009 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 @@ -28,10 +29,9 @@ #ifndef MLN_HISTO_ALL_HH # define MLN_HISTO_ALL_HH -/*! \file mln/histo/all.hh - * - * \brief File that includes histogram files. - */ +/// \file mln/histo/all.hh +/// +/// File that includes histogram files. namespace mln @@ -53,6 +53,6 @@ namespace mln } # include <mln/histo/compute.hh> -# include <mln/histo/data.hh> +# include <mln/histo/array.hh> #endif // ! MLN_HISTO_ALL_HH diff --git a/milena/mln/histo/data.hh b/milena/mln/histo/array.hh similarity index 80% rename from milena/mln/histo/data.hh rename to milena/mln/histo/array.hh index 000c727..ad58e1a 100644 --- a/milena/mln/histo/data.hh +++ b/milena/mln/histo/array.hh @@ -1,4 +1,5 @@ -// Copyright (C) 2007, 2008 EPITA Research and Development Laboratory +// Copyright (C) 2007, 2008, 2009 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 @@ -25,13 +26,12 @@ // reasons why the executable file might be covered by the GNU General // Public License. -#ifndef MLN_HISTO_DATA_HH -# define MLN_HISTO_DATA_HH +#ifndef MLN_HISTO_ARRAY_HH +# define MLN_HISTO_ARRAY_HH -/*! \file mln/histo/data.hh - * - * \brief Define a generic histogram class. - */ +/// \file mln/histo/array.hh +/// +/// Define a generic histogram class. # include <vector> # include <algorithm> @@ -46,14 +46,13 @@ namespace mln { - /*! Generic histogram class over a value set with type \c T. - */ + /// Generic histogram class over a value set with type \c T. template <typename T> - struct data + struct array { typedef T value; - data(); + array(); void clear(); @@ -75,7 +74,7 @@ namespace mln template <typename T> - std::ostream& operator<<(std::ostream& ostr, const data<T>& h); + std::ostream& operator<<(std::ostream& ostr, const array<T>& h); @@ -83,7 +82,7 @@ namespace mln template <typename T> inline - data<T>::data() + array<T>::array() : s_(mln::value::set<T>::the()), h_(s_.nvalues(), 0) { @@ -93,7 +92,7 @@ namespace mln template <typename T> inline void - data<T>::clear() + array<T>::clear() { std::fill(h_.begin(), h_.end(), 0); } @@ -101,7 +100,7 @@ namespace mln template <typename T> inline unsigned - data<T>::operator()(const T& v) const + array<T>::operator()(const T& v) const { return h_[s_.index_of(v)]; } @@ -109,7 +108,7 @@ namespace mln template <typename T> inline unsigned& - data<T>::operator()(const T& v) + array<T>::operator()(const T& v) { return h_[s_.index_of(v)]; } @@ -117,7 +116,7 @@ namespace mln template <typename T> inline const mln::value::set<T>& - data<T>::vset() const + array<T>::vset() const { return s_; } @@ -125,7 +124,7 @@ namespace mln template <typename T> inline unsigned - data<T>::operator[](unsigned i) const + array<T>::operator[](unsigned i) const { mln_precondition(i < s_.nvalues()); return h_[i]; @@ -134,7 +133,7 @@ namespace mln template <typename T> inline unsigned& - data<T>::operator[](unsigned i) + array<T>::operator[](unsigned i) { mln_precondition(i < s_.nvalues()); return h_[i]; @@ -143,21 +142,21 @@ namespace mln template <typename T> inline const std::vector<unsigned>& - data<T>::vect() const + array<T>::vect() const { return h_; } template <typename T> inline - unsigned data<T>::nvalues() const + unsigned array<T>::nvalues() const { return h_.size(); } template <typename T> inline - std::ostream& operator<<(std::ostream& ostr, const data<T>& h) + std::ostream& operator<<(std::ostream& ostr, const array<T>& h) { mln_viter(mln::value::set<T>) v(h.vset()); for_all(v) @@ -173,4 +172,4 @@ namespace mln } // end of namespace mln -#endif // ! MLN_HISTO_DATA_HH +#endif // ! MLN_HISTO_ARRAY_HH diff --git a/milena/mln/histo/compute.hh b/milena/mln/histo/compute.hh index 14bf3a3..99379c9 100644 --- a/milena/mln/histo/compute.hh +++ b/milena/mln/histo/compute.hh @@ -1,4 +1,5 @@ -// Copyright (C) 2007, 2008 EPITA Research and Development Laboratory +// Copyright (C) 2007, 2008, 2009 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 @@ -28,13 +29,12 @@ #ifndef MLN_HISTO_COMPUTE_HH # define MLN_HISTO_COMPUTE_HH -/*! \file mln/histo/compute.hh - * - * \brief Routine to compute an histogram. - */ +/// \file mln/histo/compute.hh +/// +/// Routine to compute an histogram. # include <mln/core/concept/image.hh> -# include <mln/histo/data.hh> +# include <mln/histo/array.hh> // Specializations are in: @@ -49,7 +49,7 @@ namespace mln /// Compute the histogram of image \p input. template <typename I> - data<mln_value(I)> compute(const Image<I>& input); + array<mln_value(I)> compute(const Image<I>& input); # ifndef MLN_INCLUDE_ONLY @@ -62,9 +62,9 @@ namespace mln template <typename I> inline - data<mln_value(I)> compute_(const I& input) + array<mln_value(I)> compute_(const I& input) { - data<mln_value(I)> h; + array<mln_value(I)> h; mln_piter(I) p(input.domain()); for_all(p) ++h(input(p)); @@ -78,14 +78,14 @@ namespace mln template <typename I> inline - data<mln_value(I)> compute(const Image<I>& input) + array<mln_value(I)> compute(const Image<I>& input) { trace::entering("histo::compute"); mlc_equal(mln_trait_image_quant(I), mln::trait::image::quant::low)::check(); mln_precondition(exact(input).is_valid()); - data<mln_value(I)> h = impl::compute_(mln_trait_image_speed(I)(), - exact(input)); + array<mln_value(I)> h = impl::compute_(mln_trait_image_speed(I)(), + exact(input)); trace::exiting("histo::compute"); return h; diff --git a/milena/mln/histo/compute.spe.hh b/milena/mln/histo/compute.spe.hh index 616804d..9b149b9 100644 --- a/milena/mln/histo/compute.spe.hh +++ b/milena/mln/histo/compute.spe.hh @@ -52,12 +52,12 @@ namespace mln namespace generic { template <typename I> - data<mln_value(I)> compute_(const I& input); + array<mln_value(I)> compute_(const I& input); } template <typename I> inline - data<mln_value(I)> + array<mln_value(I)> compute_(trait::image::speed::any, const I& input) { return generic::compute_(input); @@ -65,10 +65,10 @@ namespace mln template <typename I> inline - data<mln_value(I)> + array<mln_value(I)> compute_(trait::image::speed::fastest, const I& input) { - data<mln_value(I)> h; + array<mln_value(I)> h; mln_pixter(const I) p(input); for_all(p) ++h(p.val()); diff --git a/milena/mln/level/sort_psites.hh b/milena/mln/level/sort_psites.hh index ef920cd..a7bc8ba 100644 --- a/milena/mln/level/sort_psites.hh +++ b/milena/mln/level/sort_psites.hh @@ -147,7 +147,7 @@ namespace mln const unsigned n = vset.nvalues(); // h - histo::data<mln_value(I)> h = histo::compute(input); + histo::array<mln_value(I)> h = histo::compute(input); // preparing output data std::vector<unsigned> loc(vset.nvalues()); @@ -191,7 +191,7 @@ namespace mln const unsigned n = vset.nvalues(); // h - histo::data<mln_value(I)> h = histo::compute(input); + histo::array<mln_value(I)> h = histo::compute(input); // preparing output data std::vector<unsigned> loc(vset.nvalues()); -- 1.5.6.5

1 0

3156: Split ismm2009 into routines and main.
by Thierry Geraud 15 Jan '09

15 Jan '09

https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox Index: ChangeLog from Thierry Geraud <thierry.geraud(a)lrde.epita.fr> Split ismm2009 into routines and main. * laurent/ismm2009.cc: Copy to... * laurent/ismm2009.v0.cc: ...this new file (memorization). * laurent/ismm2009.cc: Copy routines to... * laurent/ismm2009.hh: ...this new file. * laurent/ismm2009.cc: Remove routines. ismm2009.cc | 379 +-------------------------------------------------------- ismm2009.hh | 275 +++-------------------------------------- ismm2009.v0.cc | 33 +--- 3 files changed, 35 insertions(+), 652 deletions(-) Index: laurent/ismm2009.cc --- laurent/ismm2009.cc (revision 3155) +++ laurent/ismm2009.cc (working copy) @@ -1,11 +1,4 @@ -#include <vector> -#include <algorithm> - -#include <mln/core/var.hh> -#include <mln/core/image/image2d.hh> -#include <mln/core/image/image_if.hh> -#include <mln/core/alias/neighb2d.hh> -#include <mln/make/double_neighb2d.hh> +#include "ismm2009.hh" #include <mln/value/int_u8.hh> #include <mln/value/label_8.hh> @@ -13,340 +6,11 @@ #include <mln/io/pgm/load.hh> #include <mln/util/ord_pair.hh> #include <mln/debug/println.hh> - -#include <mln/core/routine/extend.hh> #include <mln/core/routine/duplicate.hh> -#include <mln/data/fill.hh> -#include <mln/data/paste.hh> -#include <mln/labeling/regional_minima.hh> #include <mln/labeling/compute.hh> - #include <mln/accu/count.hh> -#include <mln/morpho/gradient.hh> -#include <mln/morpho/meyer_wst.hh> - - - -namespace mln -{ - - // Functions. - - inline - bool is_row_odd(const point2d& p) - { - return p.row() % 2; - } - - inline - bool is_edge(const point2d& p) - { - return p.row() % 2 + p.col() % 2 == 1; - } - - inline - bool is_pixel(const point2d& p) - { - return p.row() % 2 == 0 && p.col() % 2 == 0; - } - - inline - bool is_point(const point2d& p) - { - return p.row() % 2 && p.col() % 2; - } - - - // Neighborhoods. - - typedef neighb< win::multiple<window2d, bool(*)(const point2d&)> > dbl_neighb2d; - - const dbl_neighb2d& e2p() // Edge (1D face) to neighboring original pixels (2D faces). - { - static bool e2p_h[] = { 0, 1, 0, - 0, 0, 0, - 0, 1, 0 }; - static bool e2p_v[] = { 0, 0, 0, - 1, 0, 1, - 0, 0, 0 }; - static dbl_neighb2d nbh = make::double_neighb2d(is_row_odd, e2p_h, e2p_v); - return nbh; - } - - const dbl_neighb2d& e2e() // Edge to neighboring edges. - { - static bool e2e_h[] = { 0, 0, 1, 0, 0, - 0, 1, 0, 1, 0, - 0, 0, 0, 0, 0, - 0, 1, 0, 1, 0, - 0, 0, 1, 0, 0 }; - static bool e2e_v[] = { 0, 0, 0, 0, 0, - 0, 1, 0, 1, 0, - 1, 0, 0, 0, 1, - 0, 1, 0, 1, 0, - 0, 0, 0, 0, 0 }; - static dbl_neighb2d nbh = make::double_neighb2d(is_row_odd, e2e_h, e2e_v); - return nbh; - } - - - inline - point2d p1_from_e(const point2d& e) - { - return e + (is_row_odd(e) ? up : left); - } - - inline - point2d p2_from_e(const point2d& e) - { - return e + (is_row_odd(e) ? down : right); - } - - - // Transform to make room for edges. - - template <typename T> - image2d<T> - add_edges(const image2d<T>& input) - { - image2d<T> output(2 * input.nrows() - 1, - 2 * input.ncols() - 1); - data::fill(output, 0); // Useless but for display! - for (int row = 0; row < input.nrows(); ++row) - for (int col = 0; col < input.ncols(); ++col) - opt::at(output, 2 * row, 2 * col) = opt::at(input, row, col); - return output; - } - - - template <typename I, typename N> - mln_concrete(I) - magnitude(const I& input, const N& nbh) - { - mln_concrete(I) output; - initialize(output, input); - data::fill(output, 0); - - mln_piter(I) p(input.domain()); - mln_niter(N) n(nbh, p); - for_all(p) - { - n.start(); - mln_value(I) v1 = input(n); - n.next(); - mln_value(I) v2 = input(n); - output(p) = v2 > v1 ? v2 - v1 : v1 - v2; - } - - return output; - } - - - template <typename I, typename N> - mln_ch_value(I, util::ord_pair<mln_value(I)>) - adjacency(const I& input, const N& nbh) - { - typedef mln_value(I) L; - typedef util::ord_pair<L> LL; - - mln_ch_value(I, LL) output; - initialize(output, input); - - mln_piter(I) p(input.domain()); - mln_niter(N) n(nbh, p); - for_all(p) - { - if (input(p) == 0) // Watershed edge. - { - L l1 = 0, l2 = 0; - for_all(n) - if (input.has(n) && input(n) != 0) - { - if (l1 == 0) // First label to be stored. - l1 = input(n); - else - if (input(n) != l1 && l2 == 0) // Second label to be stored. - l2 = input(n); - else - mln_invariant(input(n) == l1 || input(n) == l2); - } - mln_invariant(l1 != 0 && l2 != 0); - output(p) = LL(l1, l2); - } - else - { - L l = input(p); - output(p) = LL(l, l); - // Tests: - for_all(n) - if (input.has(n)) - mln_invariant(input(n) == 0 || input(n) == l); - } - } - - return output; - } - - - - - // Get the smallest edge out of a basin. - // - // Version with the watershed extended to all faces. - - template <typename W, typename N, typename G> - std::vector<mln_site(W)> // FIXME: Use p_array! - smallest_edges(const W& wst, unsigned nlabels, - const N& nbh, // edge (1D-face) -> pixels (2D-faces) - const G& g) - { - typedef mln_value(W) L; - std::vector<mln_site(W)> edge(nlabels + 1); - std::vector<L> g_min(nlabels + 1); - std::fill(g_min.begin(), g_min.end(), mln_max(mln_value(G))); - mln_piter(W) e(wst.domain()); - mln_niter(N) n(nbh, e); - for_all(e) - { - mln_invariant(wst(e) == 0); // Watershed line only. - n.start(); - L l1 = wst(n); - n.next(); - L l2 = wst(n); - if (g(e) < g_min[l1]) - { - g_min[l1] = g(e); - edge[l1] = e; - } - if (g(e) < g_min[l2]) - { - g_min[l2] = g(e); - edge[l2] = e; - } - } - return edge; - } - - - - // Get the smallest edge out of a basin. - // - // Version with the watershed on edges only (not extended to other faces). - // This is an ALTERNATE version (just to test that we get the same result - // as the "regular" version given above). - - template <typename W, typename N, typename G> - std::vector<mln_site(W)> - smallest_edges_alt(const W& wst, unsigned nlabels, - const N& nbh, // edge -> edges - const G& g) - { - typedef mln_value(W) L; - std::vector<mln_site(W)> edge(nlabels + 1); - std::vector<L> g_min(nlabels + 1); - std::fill(g_min.begin(), g_min.end(), mln_max(mln_value(G))); - mln_piter(W) e(wst.domain()); - mln_niter(N) n(nbh, e); - for_all(e) if (wst(e) == 0) // Watershed edge. - { - L l1 = 0, l2 = 0; - for_all(n) - if (wst.has(n) && wst(n) != 0) - { - if (l1 == 0) // First label. - l1 = wst(n); - else - if (wst(n) != l1 && l2 == 0) // Second label. - l2 = wst(n); - else - mln_invariant(wst(n) == l1 || wst(n) == l2); - } - mln_invariant(l1 != 0 && l2 != 0); - if (g(e) < g_min[l1]) - { - g_min[l1] = g(e); - edge[l1] = e; - } - if (g(e) < g_min[l2]) - { - g_min[l2] = g(e); - edge[l2] = e; - } - } - return edge; - } - - - - // Sort attributes. - - template <typename A, typename L> - std::vector< std::pair<A,L> > - sort_attributes(const util::array<A>& a, L n_basins) - { - typedef std::pair<A,L> pair_t; - std::vector<pair_t> v(n_basins.next()); - - v[0] = pair_t(mln_min(A), 0); // First elt, even after sorting. - - for (L l = 1; l <= n_basins; ++l) - v[l] = pair_t(a[l], l); - - std::sort(v.begin(), v.end()); - - return v; - } - - - // Find root. - - template <typename L> - inline - L find_root(std::vector<L>& par, L l) - { - if (par[l] == l) - return l; - else - return par[l] = find_root(par, par[l]); - } - - - - // Display. - - template <typename I> - I display_edge(const I& ima, mln_value(I) bg, unsigned zoom) - { - unsigned nrows = ima.nrows() / 2 + 1; - unsigned ncols = ima.ncols() / 2 + 1; - I output(nrows * (zoom + 1) - 1, - ncols * (zoom + 1) - 1); - data::fill(output, bg); - mln_VAR( edge, ima | is_edge ); - mln_piter(edge_t) p(edge.domain()); - for_all(p) - if (p.row() % 2) // horizontal edge - { - unsigned row = (p.row() / 2 + 1) * (zoom + 1) - 1; - unsigned col = (p.col() / 2) * (zoom + 1); - for (unsigned i = 0; i < zoom; ++i) - opt::at(output, row, col + i) = ima(p); - } - else // vertical edge - { - unsigned row = (p.row() / 2) * (zoom + 1); - unsigned col = (p.col() / 2 + 1) * (zoom + 1) - 1; - for (unsigned i = 0; i < zoom; ++i) - opt::at(output, row + i, col) = ima(p); - } - return output; - } - - -} // end of namespace mln - void usage(char* argv[]) @@ -369,59 +33,30 @@ image2d<int_u8> raw_f; io::pgm::load(raw_f, argv[1]); - debug::println("raw_f:", raw_f); image2d<int_u8> f_ = add_edges(raw_f); - debug::println("f_:", f_); mln_VAR(f, f_ | is_pixel); - debug::println("f:", f); + // debug::println("f:", f); mln_VAR(g, f_ | is_edge); - data::fill(g, magnitude(extend(f_ | is_edge, pw::value(f_)), - e2p())); - debug::println("g:", g); - - // surprise: - debug::println("g without the 'edge' predicate:", g.unmorph_()); - typedef label_8 L; // Type of labels. - - L n_regmins; - mln_VAR( regmin_g, - labeling::regional_minima(g, e2e(), n_regmins) ); - debug::println("regmin(g):", regmin_g); - L n_basins; mln_VAR( wst_g, - morpho::meyer_wst(g, e2e(), n_basins) ); - mln_invariant(n_basins == n_regmins); + f_to_wst_g(f, g, e2p(), e2e(), n_basins) ); + + // 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))); - /* // Délire! - { - box2d b = make::box2d(1,1, n_basins, n_basins); - point2d null(0, 0); - image2d<point2d> adj_edge(b); - data::fill(adj_edge, null); - - mln_piter_(adj_t) e(adj.domain()); - for_all(e) - if (adj(e).first() != adj(e).second()) - adj_edge.at_(adj(e).first(), adj(e).second()) = e; - - debug::println(adj_edge); - } - */ - image2d<L> wst_g_full = wst_g.unmorph_(); { @@ -445,7 +80,7 @@ - // Get the smallest esge out of every basin. + // 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); Index: laurent/ismm2009.v0.cc --- laurent/ismm2009.v0.cc (revision 3155) +++ laurent/ismm2009.v0.cc (working copy) @@ -122,29 +122,6 @@ template <typename I, typename N> - mln_concrete(I) - magnitude(const I& input, const N& nbh) - { - mln_concrete(I) output; - initialize(output, input); - data::fill(output, 0); - - mln_piter(I) p(input.domain()); - mln_niter(N) n(nbh, p); - for_all(p) - { - n.start(); - mln_value(I) v1 = input(n); - n.next(); - mln_value(I) v2 = input(n); - output(p) = v2 > v1 ? v2 - v1 : v1 - v2; - } - - return output; - } - - - template <typename I, typename N> mln_ch_value(I, util::ord_pair<mln_value(I)>) adjacency(const I& input, const N& nbh) { @@ -378,12 +355,16 @@ debug::println("f:", f); mln_VAR(g, f_ | is_edge); - data::fill(g, magnitude(extend(f_ | is_edge, pw::value(f_)), - e2p())); + + data::paste(morpho::gradient(extend(g, f), + e2p().win()), + g); debug::println("g:", g); // surprise: debug::println("g without the 'edge' predicate:", g.unmorph_()); + std::cout << "great: we actually encode both f and g in the same image" + << "so we do save memory!" << std::endl; typedef label_8 L; // Type of labels. @@ -445,7 +426,7 @@ - // Get the smallest esge out of every basin. + // 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); Index: laurent/ismm2009.hh --- laurent/ismm2009.hh (revision 3155) +++ laurent/ismm2009.hh (working copy) @@ -7,22 +7,13 @@ #include <mln/core/alias/neighb2d.hh> #include <mln/make/double_neighb2d.hh> -#include <mln/value/int_u8.hh> -#include <mln/value/label_8.hh> - -#include <mln/io/pgm/load.hh> #include <mln/util/ord_pair.hh> -#include <mln/debug/println.hh> #include <mln/core/routine/extend.hh> -#include <mln/core/routine/duplicate.hh> #include <mln/data/fill.hh> #include <mln/data/paste.hh> #include <mln/labeling/regional_minima.hh> -#include <mln/labeling/compute.hh> - -#include <mln/accu/count.hh> #include <mln/morpho/gradient.hh> #include <mln/morpho/meyer_wst.hh> @@ -122,29 +113,6 @@ template <typename I, typename N> - mln_concrete(I) - magnitude(const I& input, const N& nbh) - { - mln_concrete(I) output; - initialize(output, input); - data::fill(output, 0); - - mln_piter(I) p(input.domain()); - mln_niter(N) n(nbh, p); - for_all(p) - { - n.start(); - mln_value(I) v1 = input(n); - n.next(); - mln_value(I) v2 = input(n); - output(p) = v2 > v1 ? v2 - v1 : v1 - v2; - } - - return output; - } - - - template <typename I, typename N> mln_ch_value(I, util::ord_pair<mln_value(I)>) adjacency(const I& input, const N& nbh) { @@ -345,238 +313,37 @@ } -} // end of namespace mln + template < typename F, + typename G, + typename N_e2p, + typename N_e2e, + typename L > + mln_ch_value(G, L) + f_to_wst_g(F& f, // On pixels. + G& g, // On edges. + const N_e2p& e2p, + const N_e2e& e2e, + L& n_basins) + { + data::paste(morpho::gradient(extend(g, f), + e2p.win()), + g); -void usage(char* argv[]) -{ - std::cerr << "usage: " << argv[0] << " input.pgm" << std::endl; - std::cerr << "Laurent ISMM 2009 scheme." << std::endl; - abort(); -} - - + mln_VAR( wst_g, + morpho::meyer_wst(g, e2e, n_basins) ); -int main(int argc, char* argv[]) + // Test the consistency with regional minima. { - using namespace mln; - using value::int_u8; - using value::label_8; - - if (argc != 2) - usage(argv); - - image2d<int_u8> raw_f; - io::pgm::load(raw_f, argv[1]); - debug::println("raw_f:", raw_f); - - image2d<int_u8> f_ = add_edges(raw_f); - debug::println("f_:", f_); - - mln_VAR(f, f_ | is_pixel); - debug::println("f:", f); - - mln_VAR(g, f_ | is_edge); - data::fill(g, magnitude(extend(f_ | is_edge, pw::value(f_)), - e2p())); - debug::println("g:", g); - - // surprise: - debug::println("g without the 'edge' predicate:", g.unmorph_()); - - - typedef label_8 L; // Type of labels. - - L n_regmins; mln_VAR( regmin_g, - labeling::regional_minima(g, e2e(), n_regmins) ); - debug::println("regmin(g):", regmin_g); - - L n_basins; - mln_VAR( wst_g, - morpho::meyer_wst(g, e2e(), n_basins) ); + labeling::regional_minima(g, e2e, n_regmins) ); mln_invariant(n_basins == n_regmins); - 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))); - - - /* // Délire! - { - box2d b = make::box2d(1,1, n_basins, n_basins); - point2d null(0, 0); - image2d<point2d> adj_edge(b); - data::fill(adj_edge, null); - - mln_piter_(adj_t) e(adj.domain()); - for_all(e) - if (adj(e).first() != adj(e).second()) - adj_edge.at_(adj(e).first(), adj(e).second()) = e; - - debug::println(adj_edge); } - */ - - image2d<L> wst_g_full = wst_g.unmorph_(); - { - // edges (1D-faces) -> pixels (2D-faces) - mln_VAR(w_pixels, wst_g_full | is_pixel); - data::paste(morpho::dilation(extend(w_pixels, pw::value(wst_g_full)), - c4().win()), - wst_g_full); - // edges (1D-faces) -> points (0D-faces) - mln_VAR(w_points, wst_g_full | is_point); - data::paste(morpho::erosion(extend(w_points, pw::value(wst_g_full)), - c4().win()), - wst_g_full); + return wst_g; } - 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); - - - - // Get the smallest esge 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; -// { -// // Test the result with an alternate code. -// std::vector<point2d> edge_alt = smallest_edges_alt(wst_g, n_basins, e2e(), g); -// for (L l = 1; l <= n_basins; ++l) -// mln_invariant(edge_alt[l] == edge[l]); -// } - - - - // Compute an attribute per region. - - typedef unsigned A; - util::array<A> a = labeling::compute(accu::meta::count(), - g, - wst_g, - n_basins); - - typedef std::pair<A,L> pair_t; - 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; - - - 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) - { - 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! - } - - for (unsigned i = 1; i <= n_basins; ++i) - { - L l = v[i].second; - std::cout << l << " -> " << lpar[l] << std::endl; - } - - -} // end of main - - - - -/* - - -for (i=0; i<V; i++) // Initialiser les heaps de fibonacci -{ - fh_setcmp(G->hp[i], cmp); // Chaque region a une heap de ses edges -} - -forall edges z that separates two regions v and w -{ - fh_insert(G->hp[v], (void *)(z)); // Ajouter les edges dans les heaps - fh_insert(G->hp[w], (void *)(z)); -} - -UFinit(G->V); // Initialiser l'union-find - -// Parcourir les regions par attribut croissant -for (j=0; j<V-1; j++) -{ - i = find(j); - - do - { // trouver l'edge minimum qui sorte de la region - e = fh_extractmin(G->hp[i]); - } - while ((UFfind(e->v, e->w)) && (e !=NULL)); - - if (e != NULL) - { // Normalement, e != NULL, sinon c'est un BIG pb!!! - int ui, uj, uk; - ui = find(e->v); - uj = find(e->w); - uk = UFunion(e->v,e->w); // Merger les regions - if (uk == ui) - { // et merger les edges - G->hp[ui] = fh_union(G->hp[ui], G->hp[uj]); - } - else - { - G->hp[uj] = fh_union(G->hp[uj], G->hp[ui]); - } - mst[k] = *e; // Garder l'arete - SaliencyWeight[k] = attribut[ui];// Poids dans la nouvelle hierarchie - OldWeight[k] = e->weight; // Poids dans l'ancienne hierarchie (inutile) - k++; - } - - // Calcul de la sortie - Pour toutes les edges separantes z=(x,y) - { - S[z] = max {SaliencyWeight[k] | sur le chemin reliant x a y dans mst} - } -} - - - */ +} // end of namespace mln

1 0