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 
15 template <typename FieldType = SIM_RawIndexField>
17 {
18 private:
19 
21  using RootPairs = std::pair<exint, exint>;
22  using ScalarType = typename FieldType::ScalarType;
23 
24  static constexpr exint UNVISITED_CELL = -2;
25  static constexpr exint VISITED_CELL = -1;
26  static constexpr bool USING_VOXEL_ARRAY = std::is_same_v<UT_VoxelArray<ScalarType>, FieldType>;
27 
28 public:
29  static constexpr exint INACTIVE_REGION = -1;
30 
31  // The provided labelled cell field should have active cells labelled as zero or greater.
32  // If FieldType is SIM_RawField or SIM_RawIndexField, then connected_regions will be aligned
33  // with label_cells. If FieldType is a UT_VoxelArray, then connected_regions will be merely
34  // of the same dimension.
35  // Faces are considered "open" wherever value of faceWeights exceeds min_weight.
37  const FieldType &label_cells,
38  const SIM_RawField **faceWeights,
39  float min_weight = 0)
40  : m_connected_regions(connected_regions)
41  , m_face_weights(faceWeights)
42  , m_min_weight(min_weight)
43  {
44  if constexpr (USING_VOXEL_ARRAY)
45  {
46  m_labelled_cells = &label_cells;
47  m_connected_regions.init(SIM_SAMPLE_CENTER, UT_Vector3(0.0f), UT_Vector3(1.0f),
48  label_cells.getXRes(), label_cells.getYRes(), label_cells.getZRes());
49  }
50  else
51  {
52  m_labelled_cells = label_cells.field();
53  m_connected_regions.match(label_cells);
54  }
55 
56  m_connected_regions.makeConstant(INACTIVE_REGION);
57  }
58 
59  template<typename IsLabelActiveFunctor>
60  exint buildConnectedComponents(const IsLabelActiveFunctor &isCellLabelActiveFunctor);
61 
62 private:
63 
65  exint
66  flattenIndex(const UT_Vector3I& index, const UT_Vector3I& size) const
67  {
68  return (index[2] * size[1] + index[1]) * size[0] + index[0];
69  }
70 
71  template<typename IsLabelActiveFunctor>
72  void buildLocalRoots(SIM_RawIndexField &visited_cells, const IsLabelActiveFunctor &isCellLabelActiveFunctor);
73 
74  THREADED_METHOD1_CONST(SIM_VolumetricConnectedComponentBuilder, m_labelled_cells->numTiles() > 20,
75  linkLocalRoots,
76  UT_Array<UT_Array<RootPairs>> &, parallel_root_pair_lists)
77 
78  void linkLocalRootsPartial(UT_Array<UT_Array<RootPairs>> &parallel_root_pair_lists,
79  const UT_JobInfo &info) const;
80 
81  THREADED_METHOD2(SIM_VolumetricConnectedComponentBuilder, m_connected_regions.shouldMultiThread(),
82  assignRegionLabels,
83  const IntegerMap &, root_map,
84  const IntegerMap &, region_index_map)
85 
86  void assignRegionLabelsPartial(const IntegerMap &root_map,
87  const IntegerMap &region_index_map,
88  const UT_JobInfo &info);
89 
90  const UT_VoxelArray<ScalarType> *m_labelled_cells;
91  SIM_RawIndexField &m_connected_regions;
92 
93  const SIM_RawField **m_face_weights;
94 
95  float m_min_weight;
96 };
97 
102 
103 template<typename FieldType>
104 template<typename IsLabelActiveFunctor>
105 exint SIM_VolumetricConnectedComponentBuilder<FieldType>::buildConnectedComponents(const IsLabelActiveFunctor &isCellLabelActiveFunctor)
106 {
107  // We want to replicate a flood fill throughout the entire contiguous
108  // volume without actually flooding in a BFS sense.
109 
110  // First we perform a flood locally within each tile and point every voxel
111  // in a connected region to a single representative cell that acts as a local root.
112  {
113  SIM_RawIndexField visited_cells;
114  visited_cells.match(m_connected_regions);
115  visited_cells.makeConstant(UNVISITED_CELL);
116 
117  buildLocalRoots(visited_cells, isCellLabelActiveFunctor);
118  }
119 
120  const int thread_count = UT_Thread::getNumProcessors();
121  UT_Array<UT_Array<RootPairs>> parallel_root_pair_lists;
122  parallel_root_pair_lists.setSize(thread_count);
123 
124  // Once we have a local root for each connected component in each tile, we can
125  // link the roots together across tile boundaries. Because there are far fewer tiles than voxels,
126  // this process is extremely fast.
127 
128  linkLocalRoots(parallel_root_pair_lists);
129 
130  // Compile the pairs of roots into a single list
131  UT_Array<RootPairs> root_pair_list;
132  exint list_size = 0;
133  for (int thread = 0; thread < thread_count; ++thread)
134  list_size += parallel_root_pair_lists[thread].size();
135 
136  root_pair_list.bumpCapacity(list_size);
137 
138  for (int thread = 0; thread < thread_count; ++thread)
139  root_pair_list.concat(parallel_root_pair_lists[thread]);
140 
141  // Build a hash map with the local root values as the keys and a list of all the roots they point to
142  // as the value. Since this is done at the tile level, there won't be many entries and it doesn't need
143  // to be parallelized.
144 
145  UT_Map<exint, UT_Array<exint>> root_to_adjacent_root_map;
146  root_to_adjacent_root_map.reserve(root_pair_list.size());
147 
148  // Sort the root pair list so we can pull from the hash table once per local root and fill its list.
149  UTparallelSort(root_pair_list.begin(), root_pair_list.end(), [](const RootPairs &pair0, const RootPairs &pair1)
150  {
151  if (pair0.first < pair1.first) return true;
152  return false;
153  });
154 
155  const exint rpl_size = root_pair_list.size();
156  for (int i = 0; i < rpl_size; )
157  {
158  exint root_key = root_pair_list[i].first;
159 
160  auto &local_root_list = root_to_adjacent_root_map[root_key];
161 
162  // Iterate over entries in the list that correspond to the same root key
163  // and append the adjacent roots to the list
164  while (i < rpl_size && root_pair_list[i].first == root_key)
165  {
166  exint adjacent_root_key = root_pair_list[i].second;
167 
168  // If the list of pairs was built correctly there shouldn't be any repeated pairs.
169  UT_ASSERT(local_root_list.find(adjacent_root_key) == -1);
170 
171  local_root_list.append(adjacent_root_key);
172  ++i;
173  }
174  }
175 
176 #if UT_ASSERT_LEVEL > 0
177  // Debug check: make sure that every link has a reciprocating link back
178  for (const auto& local_root_map : root_to_adjacent_root_map)
179  {
180  exint local_root_key = local_root_map.first;
181  const auto &local_root_list = local_root_map.second;
182 
183  for (const exint adjacent_root_key : local_root_list)
184  {
185  // Make sure the link exists
186  UT_ASSERT(root_to_adjacent_root_map.contains(adjacent_root_key));
187 
188  const auto &adjacent_root_list = root_to_adjacent_root_map[adjacent_root_key];
189 
190  // Make sure that the root we are linking to does indeed link back.
191  UT_ASSERT(adjacent_root_list.find(local_root_key) >= 0);
192  }
193  }
194 #endif
195 
196  // Now that the root connections are compiled, we can walk through the implied graph of root connections
197  // per region. Each root walks through the graph and stores the largest root they come across. After walking
198  // through the graph, the final root is stored in a separate hash map.
199 
200  IntegerMap final_root_map;
201  final_root_map.reserve(root_pair_list.size());
202 
203  IntegerMap visited_roots_map;
204  visited_roots_map.reserve(root_pair_list.size());
205 
206  {
207  UT_Array<exint> to_visit_list;
208  to_visit_list.bumpCapacity(root_pair_list.size());
209 
210  // This list stores all the nodes in the graph within
211  // a connected component so we only have to walk through
212  // the graph once to map all the local nodes to a global root
213  // within the connected component.
214  UT_Array<exint> roots_to_finish_list;
215  roots_to_finish_list.bumpCapacity(root_pair_list.size());
216 
217  for (const auto &local_root_map : root_to_adjacent_root_map)
218  {
219  exint final_root_index = local_root_map.first;
220 
221  // If this current root has been visited already, we don't need to walk through the graph.
222  if (visited_roots_map.contains(final_root_index))
223  {
224  UT_ASSERT(visited_roots_map[final_root_index] == VISITED_CELL);
225  continue;
226  }
227 
228  to_visit_list.clear();
229  to_visit_list.append(final_root_index);
230 
231  visited_roots_map[final_root_index] = VISITED_CELL;
232 
233  roots_to_finish_list.clear();
234  roots_to_finish_list.append(final_root_index);
235 
236  // Walk through the graph without retracing any steps. Record the largest
237  // root index along the way.
238  while (!to_visit_list.isEmpty())
239  {
240  exint local_root_index = to_visit_list.last();
241  to_visit_list.removeLast();
242 
243  // Store the largest index as the global root in the connected component chain.
244  final_root_index = SYSmax(final_root_index, local_root_index);
245 
246  // Queue up the new list. Only insert unvisited roots.
247  const auto &adjacent_root_list = root_to_adjacent_root_map[local_root_index];
248  for (const exint adjacent_root_index : adjacent_root_list)
249  {
250  // We should only add a local root to the list if it has not yet been added to the list
251  // and has not been visited previously.
252  if (!visited_roots_map.contains(adjacent_root_index))
253  {
254  visited_roots_map[adjacent_root_index] = VISITED_CELL;
255  roots_to_finish_list.append(adjacent_root_index);
256  to_visit_list.append(adjacent_root_index);
257  }
258  else
259  {
260  UT_ASSERT(visited_roots_map[adjacent_root_index] == VISITED_CELL);
261  }
262  }
263  }
264 
265  // Update all nodes in connected component
266  for (const exint root : roots_to_finish_list)
267  final_root_map[root] = final_root_index;
268  }
269  }
270 
271  // With each local root mapping to a global root per region, we need to give those global
272  // roots a unique indentifier.
273  IntegerMap region_index_map;
274 
275  exint region_index = 0;
276 
277  for (const auto &root : final_root_map)
278  {
279  // If the global root for a region has not yet been created,
280  // make a new entry in the hash and assign it a unique index.
281  if (!region_index_map.contains(root.second))
282  region_index_map[root.second] = region_index++;
283  }
284 
285  // Assign region number to all valid voxels.
286  assignRegionLabels(final_root_map, region_index_map);
287 
288  m_connected_regions.fieldNC()->collapseAllTiles();
289 
290  return region_index;
291 }
292 
293 template<typename FieldType>
294 template<typename IsCellLabelActiveFunctor>
295 void
297  const IsCellLabelActiveFunctor &isCellLabelActiveFunctor)
298 {
299  using namespace SIM::FieldUtils;
300 
301  UT_Interrupt *boss = UTgetInterrupt();
302 
303  UT_Vector3I voxel_res = m_labelled_cells->getVoxelRes();
304 
305  const exint tiles = m_labelled_cells->numTiles();
306 
308  {
310  vit.setConstArray(m_labelled_cells);
311 
313 
314  // TODO: change from for each number to something that uses a larger batch.
315  // Allocating this list for every time could be slow
316  UT_Array<UT_Vector3I> to_visit_cell_list;
317 
318  // TILESIZE is a Houdini defined value for the number of
319  // voxels in a dimension of a voxel tile.
320  to_visit_cell_list.bumpCapacity(TILESIZE * TILESIZE * TILESIZE);
321 
322  for (exint i = range.begin(); i != range.end(); ++i)
323  {
324  vit.myTileStart = i;
325  vit.myTileEnd = i + 1;
326  vit.rewind();
327 
328  if (boss->opInterrupt())
329  return;
330 
331  if (!vit.isTileConstant() || isCellLabelActiveFunctor(vit.getValue()))
332  {
333  UT_Vector3I tile_start, tile_end;
334  vit.getTileVoxels(tile_start, tile_end);
335 
336  vitt.setTile(vit);
337  for (vitt.rewind(); !vitt.atEnd(); vitt.advance())
338  {
339  UT_Vector3I cell(vitt.x(), vitt.y(), vitt.z());
340 
341  UT_ASSERT(cell[0] >= tile_start[0] && cell[0] < tile_end[0] &&
342  cell[1] >= tile_start[1] && cell[1] < tile_end[1] &&
343  cell[2] >= tile_start[2] && cell[2] < tile_end[2]);
344 
345  if (isCellLabelActiveFunctor(vitt.getValue()) && getFieldValue(visited_cells, cell) == UNVISITED_CELL)
346  {
347  const exint flat_root = flattenIndex(cell, voxel_res);
348 
349  setFieldValue(visited_cells, cell, VISITED_CELL);
350  setFieldValue(m_connected_regions, cell, flat_root);
351 
352  to_visit_cell_list.clear();
353  to_visit_cell_list.append(cell);
354 
355  // Begin the flood fill within the tile.
356  while (!to_visit_cell_list.isEmpty())
357  {
358  UT_Vector3I local_cell = to_visit_cell_list.last();
359  to_visit_cell_list.removeLast();
360 
361  // Load up adjacent cells that have not yet been labelled.
362  for (int axis : {0,1,2})
363  for (int direction : {0,1})
364  {
365  UT_Vector3I adjacent_cell = cellToCellMap(local_cell, axis, direction);
366 
367  // Check that the adjacent cell is within the same tile.
368  if (adjacent_cell[axis] < tile_start[axis] || adjacent_cell[axis] >= tile_end[axis])
369  continue;
370 
371  // We don't want to touch a cell that has already been visited.
372  // We ignore any cells that do not have a valid label.
373  if (getFieldValue(visited_cells, adjacent_cell) == UNVISITED_CELL &&
374  isCellLabelActiveFunctor((*m_labelled_cells)(adjacent_cell[0],
375  adjacent_cell[1],
376  adjacent_cell[2])))
377  {
378  UT_Vector3I face = cellToFaceMap(local_cell, axis, direction);
379 
380  // We only want to connect cells across non-zero face weights.
381  if (getFieldValue(*m_face_weights[axis], face) > m_min_weight)
382  {
383  setFieldValue(visited_cells, adjacent_cell, VISITED_CELL);
384  setFieldValue(m_connected_regions, adjacent_cell, flat_root);
385 
386  to_visit_cell_list.append(adjacent_cell);
387  }
388  }
389  }
390  }
391  }
392  }
393  }
394  }
395  });
396 
397  // Compress down the connected region tiles
398  m_connected_regions.fieldNC()->collapseAllTiles();
399 }
400 
401 #endif
void UTparallelSort(RandomAccessIterator begin, RandomAccessIterator end, const Compare &compare)
#define SYSmax(a, b)
Definition: SYS_Math.h:1538
T & last()
Definition: UT_Array.h:796
GLint first
Definition: glcorearb.h:405
GLenum GLint * range
Definition: glcorearb.h:1925
Unsorted map container.
Definition: UT_Map.h:107
void bumpCapacity(exint min_capacity)
Definition: UT_Array.h:612
bool contains(const Key &key) const
Definition: UT_ArrayMap.h:334
SYS_FORCE_INLINE void removeLast()
Definition: UT_Array.h:379
void UTparallelForEachNumber(IntType nitems, const Body &body, const bool force_use_task_scope=true)
IMF_EXPORT IMATH_NAMESPACE::V3f direction(const IMATH_NAMESPACE::Box2i &dataWindow, const IMATH_NAMESPACE::V2f &pixelPosition)
UT_Vector3T< float > UT_Vector3
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:982
bool atEnd() const
Returns true if we have iterated over all of the voxels.
exint size() const
Definition: UT_Array.h:646
void setSize(exint newsize)
Definition: UT_Array.h:666
void rewind()
Resets the iterator to point to the first voxel.
SYS_FORCE_INLINE void setFieldValue(SIM_RawField &field, const UT_Vector3I &cell, fpreal32 value)
GLfloat f
Definition: glcorearb.h:1926
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.
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:986
exint append()
Definition: UT_Array.h:142
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)
GLsizeiptr size
Definition: glcorearb.h:664
**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.
SIM_EXTERN_TEMPLATE(SIM_VolumetricConnectedComponentBuilder< SIM_RawField >)
void setConstArray(const UT_VoxelArray< T > *vox)
UT_API UT_Interrupt * UTgetInterrupt()
Obtain global UT_Interrupt singleton.
GLuint index
Definition: glcorearb.h:786
SYS_FORCE_INLINE UT_Vector3I cellToFaceMap(const UT_Vector3I &cell, const int axis, const int direction)
#define SIM_API
Definition: SIM_API.h:12
exint buildConnectedComponents(const IsLabelActiveFunctor &isCellLabelActiveFunctor)
#define UT_ASSERT(ZZ)
Definition: UT_Assert.h:156
void match(const SIM_RawField &src)
void clear()
Resets list to an empty list.
Definition: UT_Array.h:716
Declare prior to use.
iterator end()
End iterator.
Definition: UT_Array.h:991
bool isEmpty() const
Returns true iff there are no occupied elements in the array.
Definition: UT_Array.h:650
SIM_VolumetricConnectedComponentBuilder(SIM_RawIndexField &connected_regions, const FieldType &label_cells, const SIM_RawField **faceWeights, float min_weight=0)
void reserve(size_type num_items)
Definition: UT_ArraySet.h:464