1 // Copyright (C) 2011 Jordi Gutiérrez Hermoso <jordigh@octave.org>
2 // Copyright (C) 2014 Carnë Draug <carandraug@octave.org>
3 //
4 // This program is free software; you can redistribute it and/or modify it under
5 // the terms of the GNU General Public License as published by the Free Software
6 // Foundation; either version 3 of the License, or (at your option) any later
7 // version.
8 //
9 // This program is distributed in the hope that it will be useful, but WITHOUT
10 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11 // FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12 // details.
13 //
14 // You should have received a copy of the GNU General Public License along with
15 // this program; if not, see <http://www.gnu.org/licenses/>.
16 
17 #include <vector>
18 
19 struct voxel
20 {
21   octave_idx_type rank;
22   octave_idx_type parent;
23 
24   voxel () = default;
25 };
26 
27 class union_find
28 {
29   // Union-find data structure, see e.g.
30   // http://en.wikipedia.org/wiki/Union-find
31 
32 private:
33   std::vector<voxel> voxels;
34 
35 public:
36 
union_find(octave_idx_type s)37   explicit union_find (octave_idx_type s) : voxels (s) {};
38 
39   // Use only when adding new elements for the first time
40   void
add(const octave_idx_type idx)41   add (const octave_idx_type idx)
42   {
43     voxels[idx].parent  = idx;
44     voxels[idx].rank    = 0;
45     return;
46   }
47 
48   // Give the root representative id for this object
49   octave_idx_type
find(const octave_idx_type idx)50   find (const octave_idx_type idx)
51   {
52     voxel* elt = &voxels[idx];
53     if (elt->parent != idx)
54       elt->parent = find (elt->parent);
55 
56     return elt->parent;
57   }
58 
59   //Given two objects, unite the sets to which they belong
60   void
unite(const octave_idx_type idx1,const octave_idx_type idx2)61   unite (const octave_idx_type idx1, const octave_idx_type idx2)
62   {
63     octave_idx_type root1 = find (idx1);
64     octave_idx_type root2 = find (idx2);
65 
66     //Check if any union needs to be done, maybe they already are
67     //in the same set.
68     if (root1 != root2)
69       {
70         voxel* v1 = &voxels[root1];
71         voxel* v2 = &voxels[root2];
72         if (v1->rank > v2->rank)
73           v1->parent = root2;
74         else if (v1->rank < v2->rank)
75           v2->parent = root1;
76         else
77           {
78             v2->parent = root1;
79             v1->rank++;
80           }
81       }
82   }
83 
84   std::vector<octave_idx_type>
get_ids(const NDArray & L)85   get_ids (const NDArray& L) const
86   {
87     std::vector<octave_idx_type> ids;
88     const double* v = L.fortran_vec ();
89 
90     for (size_t i = 0; i < voxels.size (); i++)
91       if (v[i])
92         ids.push_back (i);
93 
94     return ids;
95   };
96 
97 };
98