[Getfem-commits] r4689 - in /trunk/getfem/src: bgeot_rtree.cc getfem/bgeot_rtree.h

[Yves . Renard]

Author: renard
Date: Thu Jun 19 15:08:11 2014
New Revision: 4689

URL: http://svn.gna.org/viewcvs/getfem?rev=4689&view=rev
Log:
adding the search of boxes intersection a line

Modified:
    trunk/getfem/src/bgeot_rtree.cc
    trunk/getfem/src/getfem/bgeot_rtree.h

Modified: trunk/getfem/src/bgeot_rtree.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/bgeot_rtree.cc?rev=4689&r1=4688&r2=4689&view=diff
==============================================================================
--- trunk/getfem/src/bgeot_rtree.cc     (original)
+++ trunk/getfem/src/bgeot_rtree.cc     Thu Jun 19 15:08:11 2014
@@ -1,6 +1,6 @@
 /*===========================================================================
  
- Copyright (C) 2000-2012 Julien Pommier
+ Copyright (C) 2000-2014 Julien Pommier
  
  This file is a part of GETFEM++
  
@@ -46,7 +46,8 @@
   };
 
   /* enlarge box to hold [a..b] */
-  static void update_box(base_node& bmin, base_node& bmax, const base_node& a, 
const base_node& b) {
+  static void update_box(base_node& bmin, base_node& bmax,
+                         const base_node& a, const base_node& b) {
     base_node::iterator itmin=bmin.begin(), itmax=bmax.begin();
     for (size_type i=0; i < a.size(); ++i) {
       itmin[i] = std::min(itmin[i], a.at(i));
@@ -71,7 +72,8 @@
   /* some predicates for searches */
   struct intersection_p {
     const base_node min,max;
-    intersection_p(const base_node& min_, const base_node& max_) : min(min_), 
max(max_) {}
+    intersection_p(const base_node& min_, const base_node& max_)
+      : min(min_), max(max_) {}
     bool operator()(const base_node& min2, const base_node& max2)
     { return r1_inter_r2(min,max,min2,max2); }
     bool accept(const base_node& min2, const base_node& max2) 
@@ -81,7 +83,8 @@
   /* match boxes containing [min..max] */
   struct contains_p {
     const base_node min,max;
-    contains_p(const base_node& min_, const base_node& max_) : min(min_), 
max(max_) {}
+    contains_p(const base_node& min_, const base_node& max_)
+      : min(min_), max(max_) {}
     bool operator()(const base_node& min2, const base_node& max2)
     { return r1_ge_r2(min2,max2,min,max); }
     bool accept(const base_node& min2, const base_node& max2) 
@@ -91,13 +94,15 @@
   /* match boxes contained in [min..max] */
   struct contained_p {
     const base_node min,max;
-    contained_p(const base_node& min_, const base_node& max_) : min(min_), 
max(max_) {}
+    contained_p(const base_node& min_, const base_node& max_)
+      : min(min_), max(max_) {}
     bool accept(const base_node& min2, const base_node& max2)
     { return r1_inter_r2(min,max,min2,max2); }
     bool operator()(const base_node& min2, const base_node& max2) 
     { return r1_ge_r2(min,max,min2,max2); }
   };
 
+  /* match boxes containing P */
   struct has_point_p {
     const base_node P;
     has_point_p(const base_node& P_) : P(P_) {}
@@ -110,19 +115,49 @@
     { return operator()(min2,max2); }
   };
 
+  /* match boxes intersection with the line passing by org and of
+     direction vector dirv.*/
+  struct intersect_line {
+    const base_node org;
+    const base_small_vector dirv;
+    intersect_line(const base_node& org_, const base_small_vector &dirv_)
+      : org(org_), dirv(dirv_) {}
+    bool operator()(const base_node& min2, const base_node& max2) {
+      size_type N = org.size();
+      GMM_ASSERT1(N == min2.size(), "Dimensions mismatch");
+      for (size_type i = 0; i < N; ++i)
+        if (dirv[i] != scalar_type(0)) {
+          scalar_type a1=(min2[i]-org[i])/dirv[i], a2=(max2[i]-org[i])/dirv[i];
+          bool interf1 = true, interf2 = true;
+          for (size_type j = 0; j < N; ++j)
+            if (j != i) {
+              scalar_type y1 = org[j] + a1*dirv[j], y2 = org[j] + a2*dirv[j];
+              if (y1 < min2[j] || y1 > max2[j]) interf1 = false;
+              if (y2 < min2[j] || y2 > max2[j]) interf2 = false;
+            }
+          if (interf1 || interf2) return true;
+        }
+      return false;
+    }
+    bool accept(const base_node& min2, const base_node& max2) 
+    { return operator()(min2,max2); }
+  };
+
   template <typename Predicate>
   static void find_matching_boxes_(rtree_elt_base *n, rtree::pbox_set& boxlst,
                                   Predicate p) {
-    //cout << "find_matching_boxes_: " << n->rmin << ".." << n->rmax << "\n";
+    // cout << "find_matching_boxes_: " << n->rmin << ".." << n->rmax << "\n";
     if (n->isleaf()) {
-      //cout << "find_matching_boxes_ in leaf\n";
+      // cout << "find_matching_boxes_ in leaf\n";
       const rtree_leaf *rl = static_cast<rtree_leaf*>(n);
-      for (rtree::pbox_cont::const_iterator it = rl->lst.begin(); it != 
rl->lst.end(); ++it) {
-       //cout << "  ->match(" << (*it)->id << "=" << (*it)->min << "," << 
(*it)->max << " -> " << p((*it)->min, (*it)->max) << "\n";
+      for (rtree::pbox_cont::const_iterator it = rl->lst.begin();
+           it != rl->lst.end(); ++it) {
+       // cout << "  ->match(" << (*it)->id << "=" << (*it)->min << ","
+        //      << (*it)->max << " -> " << p((*it)->min, (*it)->max) << "\n";
        if (p((*it)->min, (*it)->max)) { boxlst.insert(*it); }
       }
     } else {
-      //cout << "find_matching_boxes_ in branch\n";
+      // cout << "find_matching_boxes_ in branch\n";
       const rtree_node *rn = static_cast<rtree_node*>(n);
       if (p.accept(rn->left->rmin,rn->left->rmax)) 
        bgeot::find_matching_boxes_(rn->left, boxlst, p);
@@ -131,26 +166,38 @@
     }
   }
 
-  void rtree::find_intersecting_boxes(const base_node& bmin, const base_node& 
bmax, pbox_set& boxlst) {
+  void rtree::find_intersecting_boxes(const base_node& bmin,
+                                      const base_node& bmax,
+                                      pbox_set& boxlst) {
     boxlst.clear(); if (!root) build_tree(); 
     if (root) find_matching_boxes_(root, boxlst, intersection_p(bmin,bmax));
-    //cout << "find_intersecting_boxes : found " << boxlst.size() << " 
matches\n";
-  }
-
-  void rtree::find_containing_boxes(const base_node& bmin, const base_node& 
bmax, pbox_set& boxlst) {
+    // cout << "find_intersecting_boxes : found " << boxlst.size()
+    //      << " matches\n";
+  }
+
+  void rtree::find_containing_boxes(const base_node& bmin,
+                                    const base_node& bmax, pbox_set& boxlst) {
     boxlst.clear(); if (!root) build_tree();
     if (root) find_matching_boxes_(root, boxlst, contains_p(bmin,bmax));
   }
 
-  void rtree::find_contained_boxes(const base_node& bmin, const base_node& 
bmax, pbox_set& boxlst) {
+  void rtree::find_contained_boxes(const base_node& bmin,
+                                   const base_node& bmax, pbox_set& boxlst) {
     boxlst.clear(); if (!root) build_tree();
     if (root) find_matching_boxes_(root, boxlst, contained_p(bmin,bmax));
-    //cout << "find_matching_boxes : found " << boxlst.size() << " matches\n";
+    // cout << "find_matching_boxes : found " << boxlst.size() << " matches\n";
   }
 
   void rtree::find_boxes_at_point(const base_node& P, pbox_set& boxlst) {
     boxlst.clear(); if (!root) build_tree();
     if (root) find_matching_boxes_(root, boxlst, has_point_p(P));
+  }
+
+  void rtree::find_line_intersecting_boxes(const base_node& org,
+                                           const base_small_vector& dirv,
+                                           pbox_set& boxlst) {
+    boxlst.clear(); if (!root) build_tree();
+    if (root) find_matching_boxes_(root, boxlst, intersect_line(org, dirv));
   }
 
   /* 
@@ -162,59 +209,71 @@
                         unsigned dir, scalar_type& split_v) {
     scalar_type v = bmin[dir] + (bmax[dir] - bmin[dir])/2; split_v = v;
     size_type cnt = 0;
-    //cout << "[enter]Split_test: dir=" << dir << ", split_v=" << v << ", 
bmin=" << bmin << ", bmax=" << bmax << "\n";
-    for (rtree::pbox_cont::const_iterator it = b.begin(); it != b.end(); ++it) 
{
+    // cout << "[enter]Split_test: dir=" << dir << ", split_v=" << v
+    //      << ", bmin=" << bmin << ", bmax=" << bmax << "\n";
+    for (rtree::pbox_cont::const_iterator it = b.begin(); it!=b.end(); ++it) {
       if ((*it)->max[dir] < v) {
        if (cnt == 0) split_v = (*it)->max[dir]; 
        else split_v = std::max((*it)->max[dir],split_v);
        cnt++; 
       }
     }
-    //cout << "[exit] Split_test cnt = " << cnt << ", b.size()=" << b.size() 
<< ", split_v=" << split_v << "\n";
+    // cout << "[exit] Split_test cnt = " << cnt << ", b.size()="
+    //      << b.size() << ", split_v=" << split_v << "\n";
     return (cnt > 0 && cnt < b.size());
   }
 
   /* 
      there are many flavors of rtree ... this one is more or less a quadtree
-     where splitting does not occurs at predefined positions (hence the 
split_test
-     function above). Regions of the tree do not overlap (box are splitted).
+     where splitting does not occurs at predefined positions (hence the
+     split_test function above).
+     Regions of the tree do not overlap (box are splitted).
   */
   static rtree_elt_base* build_tree_(rtree::pbox_cont b, 
-                                    const base_node& bmin, const base_node& 
bmax,
+                                    const base_node& bmin,
+                                     const base_node& bmax,
                                     unsigned last_dir) {
     size_type N=bmin.size();
     scalar_type split_v(0); 
     unsigned split_dir = unsigned((last_dir+1)%N);
-    //cout << " build_tree_ [b.size=" << b.size() << "], bmin=" << bmin << ", 
bmax=" << bmax << "\n";
+    // cout << " build_tree_ [b.size=" << b.size() << "], bmin=" << bmin
+    //      << ", bmax=" << bmax << "\n";
     bool split_ok = false;
     if (b.size() > rtree_elt_base::RECTS_PER_LEAF) {
       for (size_type itry=0; itry < N; ++itry) {
        //cout << "split_test: dir=" << split_dir << "\n";
-       if (split_test(b, bmin, bmax, split_dir, split_v)) { split_ok = true; 
break; }
+       if (split_test(b, bmin, bmax, split_dir, split_v))
+          { split_ok = true; break; }
        split_dir = unsigned((split_dir+1)%N);
       }
-      //if (!split_ok && b.size() > rtree_elt_base::RECTS_PER_LEAF*2) cout << 
"FAILED TO SPLIT ...\n";
+      // if (!split_ok && b.size() > rtree_elt_base::RECTS_PER_LEAF*2)
+      //   cout << "FAILED TO SPLIT ...\n";
     }
     if (split_ok) {
       size_type cnt1=0,cnt2=0;
-      //cout << "splitting with v=" << split_v << "\n";
-      for (rtree::pbox_cont::const_iterator it = b.begin(); it != b.end(); 
++it) {
-       //cout << " . test box" << (*it)->min[split_dir] << ".." << 
(*it)->max[split_dir] << "\n";
+      // cout << "splitting with v=" << split_v << "\n";
+      for (rtree::pbox_cont::const_iterator it = b.begin();
+           it != b.end(); ++it) {
+       // cout << " . test box" << (*it)->min[split_dir] << ".."
+        //      << (*it)->max[split_dir] << "\n";
        if ((*it)->min[split_dir] < split_v) cnt1++; 
        if ((*it)->max[split_dir] > split_v) cnt2++;
       }
-      //cout << "  -> left : " << cnt1 << " boxes, right : " << cnt2 << " 
boxes\n";
+      // cout << "  -> left : " << cnt1 << " boxes, right : "
+      //      << cnt2 << " boxes\n";
       assert(cnt1); assert(cnt2);
       GMM_ASSERT1(cnt1+cnt2 >= b.size(), "internal error");
       rtree::pbox_cont v1(cnt1), v2(cnt2);
       base_node bmin1(bmax), bmax1(bmin); 
       base_node bmin2(bmax), bmax2(bmin);
       cnt1 = cnt2 = 0;
-      for (rtree::pbox_cont::const_iterator it = b.begin(); it != b.end(); 
++it) {
+      for (rtree::pbox_cont::const_iterator it = b.begin();
+           it != b.end(); ++it) {
        if ((*it)->min[split_dir] < split_v) {
          v1[cnt1++] = *it; 
          update_box(bmin1,bmax1,(*it)->min,(*it)->max);
-         //cout << "update_box bmin1=" << bmin1 << ", bmax1=" << bmax1 << "\n";
+         // cout << "update_box bmin1=" << bmin1 << ", bmax1="
+          //      << bmax1 << "\n";
        }
        if ((*it)->max[split_dir] > split_v) {
          v2[cnt2++] = *it;
@@ -239,7 +298,6 @@
   }
 
  void rtree::build_tree() {
-   //cout << "build tree\n";
    if (boxes.size() == 0) return;
    getfem::local_guard lock = locks_.get_lock();
    assert(root == 0);
@@ -276,7 +334,8 @@
     if (p->isleaf()) {
       rtree_leaf *rl = static_cast<rtree_leaf*>(p);
       cout << "Leaf [" << rl->lst.size() << " elts] = ";
-      for (size_type i=0; i < rl->lst.size(); ++i) cout << " " << 
rl->lst[i]->id; cout << "\n";
+      for (size_type i=0; i < rl->lst.size(); ++i)
+        cout << " " << rl->lst[i]->id; cout << "\n";
       count += rl->lst.size();
     } else {
       cout << "Node\n";

Modified: trunk/getfem/src/getfem/bgeot_rtree.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem/bgeot_rtree.h?rev=4689&r1=4688&r2=4689&view=diff
==============================================================================
--- trunk/getfem/src/getfem/bgeot_rtree.h       (original)
+++ trunk/getfem/src/getfem/bgeot_rtree.h       Thu Jun 19 15:08:11 2014
@@ -1,7 +1,7 @@
 /* -*- c++ -*- (enables emacs c++ mode) */
 /*===========================================================================
  
- Copyright (C) 2000-2012 Julien Pommier
+ Copyright (C) 2000-2014 Julien Pommier
  
  This file is a part of GETFEM++
  
@@ -64,32 +64,61 @@
     rtree() : root(0) {}
     ~rtree() { destroy_tree(); }
     void add_box(base_node min, base_node max, size_type id=size_type(-1)) {
-      box_index bi; bi.min = min; bi.max = max; bi.id = (id + 1) ? id : 
boxes.size();
+      box_index bi; bi.min = min; bi.max = max;
+      bi.id = (id + 1) ? id : boxes.size();
       boxes.push_back(bi);
     }
     size_type nb_boxes() const { return boxes.size(); }
     void clear() { destroy_tree(); boxes.clear(); }
 
-    void find_intersecting_boxes(const base_node& bmin, const base_node& bmax, 
pbox_set& boxlst);
-    void find_containing_boxes(const base_node& bmin, const base_node& bmax, 
pbox_set& boxlst);
-    void find_contained_boxes(const base_node& bmin, const base_node& bmax, 
pbox_set& boxlst);
+    void find_intersecting_boxes(const base_node& bmin, const base_node& bmax,
+                                 pbox_set& boxlst);
+    void find_containing_boxes(const base_node& bmin, const base_node& bmax,
+                               pbox_set& boxlst);
+    void find_contained_boxes(const base_node& bmin, const base_node& bmax,
+                              pbox_set& boxlst);
     void find_boxes_at_point(const base_node& P, pbox_set& boxlst);
+    void find_line_intersecting_boxes(const base_node& org,
+                                      const base_small_vector& dirv,
+                                      pbox_set& boxlst);
 
-    void find_intersecting_boxes(const base_node& bmin, const base_node& bmax, 
std::vector<size_type>& idvec)
-    { pbox_set bs; find_intersecting_boxes(bmin, bmax, bs);  
pbox_set_to_idvec(bs, idvec); }
-    void find_containing_boxes(const base_node& bmin, const base_node& bmax, 
std::vector<size_type>& idvec)
-    { pbox_set bs; find_containing_boxes(bmin, bmax, bs);  
pbox_set_to_idvec(bs, idvec); }
-    void find_contained_boxes(const base_node& bmin, const base_node& bmax, 
std::vector<size_type>& idvec)
-    { pbox_set bs; find_contained_boxes(bmin, bmax, bs);  
pbox_set_to_idvec(bs, idvec); }
+    void find_intersecting_boxes(const base_node& bmin, const base_node& bmax,
+                                 std::vector<size_type>& idvec) {
+      pbox_set bs;
+      find_intersecting_boxes(bmin, bmax, bs);
+      pbox_set_to_idvec(bs, idvec);
+    }
+    void find_containing_boxes(const base_node& bmin, const base_node& bmax,
+                               std::vector<size_type>& idvec) {
+      pbox_set bs;
+      find_containing_boxes(bmin, bmax, bs);
+      pbox_set_to_idvec(bs, idvec);
+    }
+    void find_contained_boxes(const base_node& bmin,
+                              const base_node& bmax,
+                              std::vector<size_type>& idvec) {
+      pbox_set bs;
+      find_contained_boxes(bmin, bmax, bs);
+      pbox_set_to_idvec(bs, idvec);
+    }
     void find_boxes_at_point(const base_node& P, std::vector<size_type>& idvec)
     { pbox_set bs; find_boxes_at_point(P, bs);  pbox_set_to_idvec(bs, idvec); }
+    void find_line_intersecting_boxes(const base_node& org,
+                                      const base_small_vector& dirv,
+                                      std::vector<size_type>& idvec) {
+      pbox_set bs;
+      find_line_intersecting_boxes(org, dirv, bs);
+      pbox_set_to_idvec(bs, idvec);
+    }
+
     void dump();
     void build_tree();
   private:
     void destroy_tree();
     static void pbox_set_to_idvec(pbox_set bs, std::vector<size_type>& idvec) {
       idvec.reserve(bs.size()); idvec.resize(0);
-      for (pbox_set::const_iterator it=bs.begin(); it != bs.end(); ++it) 
idvec.push_back((*it)->id);
+      for (pbox_set::const_iterator it=bs.begin(); it != bs.end(); ++it)
+        idvec.push_back((*it)->id);
     }
 
     box_cont boxes;