HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SIM_VolumetricConnectedComponentBuilder.h
Go to the documentation of this file.
1 #ifndef __SIM__VOLUMETRICCONNECTEDCOMPONENTBUILDER_H__
2 #define __SIM__VOLUMETRICCONNECTEDCOMPONENTBUILDER_H__
3 
4 #include <SIM/SIM_FieldUtils.h>
5 
6 #include <UT/UT_ArrayMap.h>
7 #include <UT/UT_ParallelUtil.h>
8 #include <SYS/SYS_Math.h>
9 
10 #include <utility>
11 
12 // Helper class that links together regions of an index field that
13 // consist of non-negative values.
14 
16 {
17 private:
18 
20  using RootPairs = std::pair<exint, exint>;
21 
22  static constexpr exint UNVISITED_CELL = -2;
23  static constexpr exint VISITED_CELL = -1;
24 
25 public:
26  static constexpr exint INACTIVE_REGION = -1;
27 
28  // The provided labelled cell field should have active cells labelled as zero or greater.
29  // TODO: figure out a better way to get weights[3] in here.
31  const SIM_RawIndexField &label_cells,
32  const SIM_RawField **faceWeights)
33  : m_connected_regions(connected_regions)
34  , m_labelled_cells(label_cells)
35  , m_face_weights(faceWeights)
36  {
37  m_connected_regions.match(label_cells);
38  m_connected_regions.makeConstant(INACTIVE_REGION);
39  }
40 
41  template<typename IsLabelActiveFunctor>
42  exint buildConnectedComponents(const IsLabelActiveFunctor &isCellLabelActiveFunctor);
43 
44 private:
45 
47  exint
48  flattenIndex(const UT_Vector3I& index, const UT_Vector3I& size) const
49  {
50  return (index[2] * size[1] + index[1]) * size[0] + index[0];
51  }
52 
53  template<typename IsLabelActiveFunctor>
54  void buildLocalRoots(SIM_RawIndexField &visited_cells, const IsLabelActiveFunctor &isCellLabelActiveFunctor);
55 
56  THREADED_METHOD1_CONST(SIM_VolumetricConnectedComponentBuilder, m_labelled_cells.shouldMultiThread(),
57  linkLocalRoots,
58  UT_Array<UT_Array<RootPairs>> &, parallel_root_pair_lists)
59 
60  void linkLocalRootsPartial(UT_Array<UT_Array<RootPairs>> &parallel_root_pair_lists,
61  const UT_JobInfo &info) const;
62 
63  THREADED_METHOD2(SIM_VolumetricConnectedComponentBuilder, m_connected_regions.shouldMultiThread(),
64  assignRegionLabels,
65  const IntegerMap &, root_map,
66  const IntegerMap &, region_index_map)
67 
68  void assignRegionLabelsPartial(const IntegerMap &root_map,
69  const IntegerMap &region_index_map,
70  const UT_JobInfo &info);
71 
72  const SIM_RawIndexField &m_labelled_cells;
73  SIM_RawIndexField &m_connected_regions;
74 
75  const SIM_RawField **m_face_weights;
76 };
77 
78 template<typename IsLabelActiveFunctor>
79 exint SIM_VolumetricConnectedComponentBuilder::buildConnectedComponents(const IsLabelActiveFunctor &isCellLabelActiveFunctor)
80 {
81  // We want to replicate a flood fill throughout the entire contiguous
82  // volume without actually flooding in a BFS sense.
83 
84  // First we perform a flood locally within each tile and point every voxel
85  // in a connected region to a single representative cell that acts as a local root.
86  {
87  SIM_RawIndexField visited_cells;
88  visited_cells.match(m_labelled_cells);
89  visited_cells.makeConstant(UNVISITED_CELL);
90 
91  buildLocalRoots(visited_cells, isCellLabelActiveFunctor);
92  }
93 
94  const int thread_count = UT_Thread::getNumProcessors();
95  UT_Array<UT_Array<RootPairs>> parallel_root_pair_lists;
96  parallel_root_pair_lists.setSize(thread_count);
97 
98  // Once we have a local root for each connected component in each tile, we can
99  // link the roots together across tile boundaries. Because there are far fewer tiles than voxels,
100  // this process is extremely fast.
101 
102  linkLocalRoots(parallel_root_pair_lists);
103 
104  // Compile the pairs of roots into a single list
105  UT_Array<RootPairs> root_pair_list;
106  exint list_size = 0;
107  for (int thread = 0; thread < thread_count; ++thread)
108  list_size += parallel_root_pair_lists[thread].size();
109 
110  root_pair_list.bumpCapacity(list_size);
111 
112  for (int thread = 0; thread < thread_count; ++thread)
113  root_pair_list.concat(parallel_root_pair_lists[thread]);
114 
115  // Build a hash map with the local root values as the keys and a list of all the roots they point to
116  // as the value. Since this is done at the tile level, there won't be many entries and it doesn't need
117  // to be parallelized.
118 
119  UT_Map<exint, UT_Array<exint>> root_to_adjacent_root_map;
120  root_to_adjacent_root_map.reserve(root_pair_list.size());
121 
122  // Sort the root pair list so we can pull from the hash table once per local root and fill its list.
123  UTparallelSort(root_pair_list.begin(), root_pair_list.end(), [](const RootPairs &pair0, const RootPairs &pair1)
124  {
125  if (pair0.first < pair1.first) return true;
126  return false;
127  });
128 
129  const exint rpl_size = root_pair_list.size();
130  for (int i = 0; i < rpl_size; )
131  {
132  exint root_key = root_pair_list[i].first;
133 
134  auto &local_root_list = root_to_adjacent_root_map[root_key];
135 
136  // Iterate over entries in the list that correspond to the same root key
137  // and append the adjacent roots to the list
138  while (i < rpl_size && root_pair_list[i].first == root_key)
139  {
140  exint adjacent_root_key = root_pair_list[i].second;
141 
142  // If the list of pairs was built correctly there shouldn't be any repeated pairs.
143  UT_ASSERT(local_root_list.find(adjacent_root_key) == -1);
144 
145  local_root_list.append(adjacent_root_key);
146  ++i;
147  }
148  }
149 
150 #if UT_ASSERT_LEVEL > 0
151  // Debug check: make sure that every link has a reciprocating link back
152  for (const auto& local_root_map : root_to_adjacent_root_map)
153  {
154  exint local_root_key = local_root_map.first;
155  const auto &local_root_list = local_root_map.second;
156 
157  for (const exint adjacent_root_key : local_root_list)
158  {
159  // Make sure the link exists
160  UT_ASSERT(root_to_adjacent_root_map.contains(adjacent_root_key));
161 
162  const auto &adjacent_root_list = root_to_adjacent_root_map[adjacent_root_key];
163 
164  // Make sure that the root we are linking to does indeed link back.
165  UT_ASSERT(adjacent_root_list.find(local_root_key) >= 0);
166  }
167  }
168 #endif
169 
170  // Now that the root connections are compiled, we can walk through the implied graph of root connections
171  // per region. Each root walks through the graph and stores the largest root they come across. After walking
172  // through the graph, the final root is stored in a separate hash map.
173 
174  IntegerMap final_root_map;
175  final_root_map.reserve(root_pair_list.size());
176 
177  IntegerMap visited_roots_map;
178  visited_roots_map.reserve(root_pair_list.size());
179 
180  {
181  UT_Array<exint> to_visit_list;
182  to_visit_list.bumpCapacity(root_pair_list.size());
183 
184  // This list stores all the nodes in the graph within
185  // a connected component so we only have to walk through
186  // the graph once to map all the local nodes to a global root
187  // within the connected component.
188  UT_Array<exint> roots_to_finish_list;
189  roots_to_finish_list.bumpCapacity(root_pair_list.size());
190 
191  for (const auto &local_root_map : root_to_adjacent_root_map)
192  {
193  exint final_root_index = local_root_map.first;
194 
195  // If this current root has been visited already, we don't need to walk through the graph.
196  if (visited_roots_map.contains(final_root_index))
197  {
198  UT_ASSERT(visited_roots_map[final_root_index] == VISITED_CELL);
199  continue;
200  }
201 
202  to_visit_list.clear();
203  to_visit_list.append(final_root_index);
204 
205  visited_roots_map[final_root_index] = VISITED_CELL;
206 
207  roots_to_finish_list.clear();
208  roots_to_finish_list.append(final_root_index);
209 
210  // Walk through the graph without retracing any steps. Record the largest
211  // root index along the way.
212  while (!to_visit_list.isEmpty())
213  {
214  exint local_root_index = to_visit_list.last();
215  to_visit_list.removeLast();
216 
217  // Store the largest index as the global root in the connected component chain.
218  final_root_index = SYSmax(final_root_index, local_root_index);
219 
220  // Queue up the new list. Only insert unvisited roots.
221  const auto &adjacent_root_list = root_to_adjacent_root_map[local_root_index];
222  for (const exint adjacent_root_index : adjacent_root_list)
223  {
224  // We should only add a local root to the list if it has not yet been added to the list
225  // and has not been visited previously.
226  if (!visited_roots_map.contains(adjacent_root_index))
227  {
228  visited_roots_map[adjacent_root_index] = VISITED_CELL;
229  roots_to_finish_list.append(adjacent_root_index);
230  to_visit_list.append(adjacent_root_index);
231  }
232  else
233  {
234  UT_ASSERT(visited_roots_map[adjacent_root_index] == VISITED_CELL);
235  }
236  }
237  }
238 
239  // Update all nodes in connected component
240  for (const exint root : roots_to_finish_list)
241  final_root_map[root] = final_root_index;
242  }
243  }
244 
245  // With each local root mapping to a global root per region, we need to give those global
246  // roots a unique indentifier.
247  IntegerMap region_index_map;
248 
249  exint region_index = 0;
250 
251  for (const auto &root : final_root_map)
252  {
253  // If the global root for a region has not yet been created,
254  // make a new entry in the hash and assign it a unique index.
255  if (!region_index_map.contains(root.second))
256  region_index_map[root.second] = region_index++;
257  }
258 
259  // Assign region number to all valid voxels.
260  assignRegionLabels(final_root_map, region_index_map);
261 
262  m_connected_regions.fieldNC()->collapseAllTiles();
263 
264  return region_index;
265 }
266 
267 template<typename IsCellLabelActiveFunctor>
268 void
269 SIM_VolumetricConnectedComponentBuilder::buildLocalRoots(SIM_RawIndexField &visited_cells,
270  const IsCellLabelActiveFunctor &isCellLabelActiveFunctor)
271 {
272  using namespace SIM::FieldUtils;
273 
274  UT_Interrupt *boss = UTgetInterrupt();
275 
276  UT_Vector3I voxel_res = m_labelled_cells.getVoxelRes();
277 
278  const exint tiles = m_labelled_cells.field()->numTiles();
279 
281  {
283  vit.setConstArray(m_labelled_cells.field());
284 
286 
287  // TODO: change from for each number to something that uses a larger batch.
288  // Allocating this list for every time could be slow
289  UT_Array<UT_Vector3I> to_visit_cell_list;
290 
291  // TILESIZE is a Houdini defined value for the number of
292  // voxels in a dimension of a voxel tile.
293  to_visit_cell_list.bumpCapacity(TILESIZE * TILESIZE * TILESIZE);
294 
295  for (exint i = range.begin(); i != range.end(); ++i)
296  {
297  vit.myTileStart = i;
298  vit.myTileEnd = i + 1;
299  vit.rewind();
300 
301  if (boss->opInterrupt())
302  return;
303 
304  if (!vit.isTileConstant() || isCellLabelActiveFunctor(vit.getValue()))
305  {
306  UT_Vector3I tile_start, tile_end;
307  vit.getTileVoxels(tile_start, tile_end);
308 
309  vitt.setTile(vit);
310  for (vitt.rewind(); !vitt.atEnd(); vitt.advance())
311  {
312  UT_Vector3I cell(vitt.x(), vitt.y(), vitt.z());
313 
314  UT_ASSERT(cell[0] >= tile_start[0] && cell[0] < tile_end[0] &&
315  cell[1] >= tile_start[1] && cell[1] < tile_end[1] &&
316  cell[2] >= tile_start[2] && cell[2] < tile_end[2]);
317 
318  if (isCellLabelActiveFunctor(vitt.getValue()) && getFieldValue(visited_cells, cell) == UNVISITED_CELL)
319  {
320  const exint flat_root = flattenIndex(cell, voxel_res);
321 
322  setFieldValue(visited_cells, cell, VISITED_CELL);
323  setFieldValue(m_connected_regions, cell, flat_root);
324 
325  to_visit_cell_list.clear();
326  to_visit_cell_list.append(cell);
327 
328  // Begin the flood fill within the tile.
329  while (!to_visit_cell_list.isEmpty())
330  {
331  UT_Vector3I local_cell = to_visit_cell_list.last();
332  to_visit_cell_list.removeLast();
333 
334  // Load up adjacent cells that have not yet been labelled.
335  for (int axis : {0,1,2})
336  for (int direction : {0,1})
337  {
338  UT_Vector3I adjacent_cell = cellToCellMap(local_cell, axis, direction);
339 
340  // Check that the adjacent cell is within the same tile.
341  if (adjacent_cell[axis] < tile_start[axis] || adjacent_cell[axis] >= tile_end[axis])
342  continue;
343 
344  // We don't want to touch a cell that has already been visited.
345  // We ignore any cells that do not have a valid label.
346  if (getFieldValue(visited_cells, adjacent_cell) == UNVISITED_CELL &&
347  isCellLabelActiveFunctor(getFieldValue(m_labelled_cells, adjacent_cell)))
348  {
349  UT_Vector3I face = cellToFaceMap(local_cell, axis, direction);
350 
351  // We only want to connect cells across non-zero face weights.
352  if (getFieldValue(*m_face_weights[axis], face) > 0)
353  {
354  setFieldValue(visited_cells, adjacent_cell, VISITED_CELL);
355  setFieldValue(m_connected_regions, adjacent_cell, flat_root);
356 
357  to_visit_cell_list.append(adjacent_cell);
358  }
359  }
360  }
361  }
362  }
363  }
364  }
365  }
366  });
367 
368  // Compress down the connected region tiles
369  m_connected_regions.fieldNC()->collapseAllTiles();
370 }
371 
372 #endif
void UTparallelSort(RandomAccessIterator begin, RandomAccessIterator end, const Compare &compare)
#define SYSmax(a, b)
Definition: SYS_Math.h:1535
T & last()
Definition: UT_Array.h:629
void UTparallelForEachNumber(IntType nitems, const Body &body)
Unsorted map container.
Definition: UT_Map.h:83
bool contains(const Key &key) const
Definition: UT_ArrayMap.h:334
GLint first
Definition: glcorearb.h:404
SYS_FORCE_INLINE void removeLast()
Definition: UT_Array.h:230
IMF_EXPORT IMATH_NAMESPACE::V3f direction(const IMATH_NAMESPACE::Box2i &dataWindow, const IMATH_NAMESPACE::V2f &pixelPosition)
int64 exint
Definition: SYS_Types.h:125
exint concat(const UT_Array< T > &a)
Takes another T array and concatenate it onto my end.
Definition: UT_ArrayImpl.h:423
void bumpCapacity(exint mincapacity)
Definition: UT_Array.h:445
SIM_VolumetricConnectedComponentBuilder(SIM_RawIndexField &connected_regions, const SIM_RawIndexField &label_cells, const SIM_RawField **faceWeights)
bool atEnd() const
Returns true if we have iterated over all of the voxels.
exint size() const
Definition: UT_Array.h:479
void setSize(exint newsize)
Definition: UT_Array.h:499
void rewind()
Resets the iterator to point to the first voxel.
GLsizeiptr size
Definition: glcorearb.h:663
SYS_FORCE_INLINE void setFieldValue(SIM_RawField &field, const UT_Vector3I &cell, fpreal32 value)
void advance()
Advances the iterator to point to the next voxel.
SYS_FORCE_INLINE UT_Vector3I cellToCellMap(const UT_Vector3I &cell, const int axis, const int direction)
void rewind()
Resets the iterator to point to the first voxel.
exint buildConnectedComponents(const IsLabelActiveFunctor &isCellLabelActiveFunctor)
static int getNumProcessors()
#define SYS_FORCE_INLINE
Definition: SYS_Inline.h:45
#define THREADED_METHOD2(CLASSNAME, DOMULTI, METHOD, PARMTYPE1, PARMNAME1, PARMTYPE2, PARMNAME2)
void setTile(const UT_VoxelArrayIterator< S > &vit, UT_VoxelArray< T > *array)
SYS_FORCE_INLINE fpreal32 getFieldValue(const SIM_RawField &field, const UT_Vector3I &cell)
iterator begin()
Definition: UT_Array.h:817
exint append()
Definition: UT_Array.h:95
bool isTileConstant() const
Returns true if the tile we are currently in is a constant tile.
void getTileVoxels(UT_Vector3I &start, UT_Vector3I &end) const
This tile will iterate over the voxels indexed [start,end).
#define THREADED_METHOD1_CONST(CLASSNAME, DOMULTI, METHOD, PARMTYPE1, PARMNAME1)
**Note that the tasks the is the thread number *for the or if it s being executed by a non pool thread(this *can happen in cases where the whole pool is occupied and the calling *thread contributes to running the work load).**Thread pool.Have fun
void makeConstant(exint cval)
Sets this field to a constant value.
void setConstArray(const UT_VoxelArray< T > *vox)
UT_API UT_Interrupt * UTgetInterrupt()
Obtain global UT_Interrupt singleton.
GLuint index
Definition: glcorearb.h:785
SYS_FORCE_INLINE UT_Vector3I cellToFaceMap(const UT_Vector3I &cell, const int axis, const int direction)
#define SIM_API
Definition: SIM_API.h:10
GLenum GLint * range
Definition: glcorearb.h:1924
#define UT_ASSERT(ZZ)
Definition: UT_Assert.h:171
UT_VoxelArrayI * fieldNC() const
void match(const SIM_RawField &src)
void clear()
Resets list to an empty list.
Definition: UT_Array.h:549
Declare prior to use.
GLenum GLuint GLint GLenum face
Definition: glew.h:4630
iterator end()
End iterator.
Definition: UT_Array.h:822
bool isEmpty() const
Returns true iff there are no occupied elements in the array.
Definition: UT_Array.h:483
void reserve(size_type num_items)
Definition: UT_ArraySet.h:464