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