VTK-m  2.0
FlyingEdgesPass4X.h
Go to the documentation of this file.
1 
2 //============================================================================
3 // Copyright (c) Kitware, Inc.
4 // All rights reserved.
5 // See LICENSE.txt for details.
6 //
7 // This software is distributed WITHOUT ANY WARRANTY; without even
8 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
9 // PURPOSE. See the above copyright notice for more information.
10 //============================================================================
11 
12 
13 #ifndef vtk_m_worklet_contour_flyingedges_pass4x_h
14 #define vtk_m_worklet_contour_flyingedges_pass4x_h
15 
16 
20 
22 
23 namespace vtkm
24 {
25 namespace worklet
26 {
27 namespace flying_edges
28 {
29 
30 template <typename T>
32 {
33 
37 
39 
42 
44  ComputePass4X(T value,
45  const vtkm::Id3& pdims,
46  const vtkm::Vec3f& origin,
47  const vtkm::Vec3f& spacing,
48  vtkm::Id multiContourCellOffset,
49  vtkm::Id multiContourPointOffset)
50  : PointDims(pdims)
51  , Origin(origin)
52  , Spacing(spacing)
53  , IsoValue(value)
54  , CellWriteOffset(multiContourCellOffset)
55  , PointWriteOffset(multiContourPointOffset)
56  {
57  }
58 
59  using ControlSignature = void(CellSetIn,
60  FieldInPoint axis_sums,
61  FieldInPoint axis_mins,
62  FieldInPoint axis_maxs,
63  WholeArrayIn cell_tri_count,
64  WholeArrayIn edgeData,
65  WholeArrayIn data,
66  WholeArrayOut connectivity,
67  WholeArrayOut edgeIds,
68  WholeArrayOut weights,
69  WholeArrayOut inputCellIds,
70  WholeArrayOut points);
71  using ExecutionSignature =
72  void(ThreadIndices, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12, WorkIndex);
73 
74  template <typename ThreadIndices,
75  typename FieldInPointId3,
76  typename FieldInPointId,
77  typename WholeTriField,
78  typename WholeEdgeField,
79  typename WholeDataField,
80  typename WholeConnField,
81  typename WholeEdgeIdField,
82  typename WholeWeightField,
83  typename WholeCellIdField,
84  typename WholePointField>
85  VTKM_EXEC void operator()(const ThreadIndices& threadIndices,
86  const FieldInPointId3& axis_sums,
87  const FieldInPointId& axis_mins,
88  const FieldInPointId& axis_maxs,
89  const WholeTriField& cellTriCount,
90  const WholeEdgeField& edges,
91  const WholeDataField& field,
92  const WholeConnField& conn,
93  const WholeEdgeIdField& interpolatedEdgeIds,
94  const WholeWeightField& weights,
95  const WholeCellIdField& inputCellIds,
96  const WholePointField& points,
97  vtkm::Id oidx) const
98  {
99  using AxisToSum = SumXAxis;
100 
101  //This works as cellTriCount was computed with ScanExtended
102  //and therefore has one more entry than the number of cells
103  vtkm::Id cell_tri_offset = cellTriCount.Get(oidx);
104  vtkm::Id next_tri_offset = cellTriCount.Get(oidx + 1);
105  if (cell_tri_offset == next_tri_offset)
106  { //we produce nothing
107  return;
108  }
109  cell_tri_offset += this->CellWriteOffset;
110 
111  Pass4TrimState state(
112  AxisToSum{}, this->PointDims, threadIndices, axis_sums, axis_mins, axis_maxs, edges);
113  if (!state.hasWork)
114  {
115  return;
116  }
117 
118  const vtkm::Id3 pdims = this->PointDims;
119  const vtkm::Id3 increments = compute_incs3d(pdims);
120  vtkm::Id edgeIds[12];
121 
122  auto edgeCase = getEdgeCase(edges, state.startPos, (state.axis_inc * state.left));
123  init_voxelIds(AxisToSum{}, this->PointWriteOffset, edgeCase, axis_sums, edgeIds);
124  for (vtkm::Id i = state.left; i < state.right; ++i) // run along the trimmed voxels
125  {
126  edgeCase = getEdgeCase(edges, state.startPos, (state.axis_inc * i));
127  vtkm::UInt8 numTris = data::GetNumberOfPrimitives(edgeCase);
128  if (numTris > 0)
129  {
130  // Start by generating triangles for this case
132  state.cellId, edgeCase, numTris, edgeIds, cell_tri_offset, conn, inputCellIds);
133 
134  // Now generate edgeIds and weights along voxel axes if needed. Remember to take
135  // boundary into account.
136  auto* edgeUses = data::GetEdgeUses(edgeCase);
137  if (!fully_interior(state.boundaryStatus) || case_includes_axes(edgeUses))
138  {
139  this->Generate(state.boundaryStatus,
140  state.ijk,
141  field,
142  interpolatedEdgeIds,
143  weights,
144  points,
145  state.startPos,
146  increments,
147  (state.axis_inc * i),
148  edgeUses,
149  edgeIds);
150  }
151  advance_voxelIds(edgeUses, edgeIds);
152  }
153  state.increment(AxisToSum{}, pdims);
154  }
155  }
156 
157  //----------------------------------------------------------------------------
158  template <typename WholeDataField,
159  typename WholeIEdgeField,
160  typename WholeWeightField,
161  typename WholePointField>
162  VTKM_EXEC inline void Generate(const vtkm::Vec<vtkm::UInt8, 3>& boundaryStatus,
163  const vtkm::Id3& ijk,
164  const WholeDataField& field,
165  const WholeIEdgeField& interpolatedEdgeIds,
166  const WholeWeightField& weights,
167  const WholePointField& points,
168  const vtkm::Id4& startPos,
169  const vtkm::Id3& incs,
170  vtkm::Id offset,
171  vtkm::UInt8 const* const edgeUses,
172  vtkm::Id* edgeIds) const
173  {
174  using AxisToSum = SumXAxis;
175 
176  vtkm::Id2 pos(startPos[0] + offset, 0);
177  {
178  auto s0 = field.Get(pos[0]);
179 
180  //EdgesUses 0,4,8 work for Y axis
181  if (edgeUses[0])
182  { // edgesUses[0] == i axes edge
183  auto writeIndex = edgeIds[0];
184  pos[1] = startPos[0] + offset + incs[AxisToSum::xindex];
185  auto s1 = field.Get(pos[1]);
186  T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
187 
188  interpolatedEdgeIds.Set(writeIndex, pos);
189  weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
190 
191  auto coord = this->InterpolateCoordinate(t, ijk, ijk + vtkm::Id3{ 1, 0, 0 });
192  points.Set(writeIndex, coord);
193  }
194  if (edgeUses[4])
195  { // edgesUses[4] == j axes edge
196  auto writeIndex = edgeIds[4];
197  pos[1] = startPos[1] + offset;
198  auto s1 = field.Get(pos[1]);
199  T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
200 
201  interpolatedEdgeIds.Set(writeIndex, pos);
202  weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
203 
204  auto coord = this->InterpolateCoordinate(t, ijk, ijk + vtkm::Id3{ 0, 1, 0 });
205  points.Set(writeIndex, coord);
206  }
207  if (edgeUses[8])
208  { // edgesUses[8] == k axes edge
209  auto writeIndex = edgeIds[8];
210  pos[1] = startPos[2] + offset;
211  auto s1 = field.Get(pos[1]);
212  T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
213 
214  interpolatedEdgeIds.Set(writeIndex, pos);
215  weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
216 
217  auto coord = this->InterpolateCoordinate(t, ijk, ijk + vtkm::Id3{ 0, 0, 1 });
218  points.Set(writeIndex, coord);
219  }
220  }
221 
222  // On the boundary cells special work has to be done to cover the partial
223  // cell axes. These are boundary situations where the voxel axes is not
224  // fully formed. These situations occur on the +x,+y,+z volume
225  // boundaries. The other cases such as interior, or -x,-y,-z boundaries fall through
226  // which is expected
227  //
228  // clang-format off
229  const bool onX = boundaryStatus[AxisToSum::xindex] & FlyingEdges3D::MaxBoundary;
230  const bool onY = boundaryStatus[AxisToSum::yindex] & FlyingEdges3D::MaxBoundary;
231  const bool onZ = boundaryStatus[AxisToSum::zindex] & FlyingEdges3D::MaxBoundary;
232  if (onX) //+x boundary
233  {
234  this->InterpolateEdge(ijk, pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
235  this->InterpolateEdge(ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
236  if (onY) //+x +y
237  {
238  this->InterpolateEdge(ijk, pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
239  }
240  if (onZ) //+x +z
241  {
242  this->InterpolateEdge(ijk, pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
243  }
244  }
245  if (onY) //+y boundary
246  {
247  this->InterpolateEdge(ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
248  this->InterpolateEdge(ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
249  if (onZ) //+y +z boundary
250  {
251  this->InterpolateEdge(ijk, pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
252  }
253  }
254  if (onZ) //+z boundary
255  {
256  this->InterpolateEdge(ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
257  this->InterpolateEdge(ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
258  }
259  // clang-format on
260  }
261 
262  // Indicate whether voxel axes need processing for this case.
263  //----------------------------------------------------------------------------
264  template <typename WholeField,
265  typename WholeIEdgeField,
266  typename WholeWeightField,
267  typename WholePointField>
268  VTKM_EXEC inline void InterpolateEdge(const vtkm::Id3& ijk,
269  vtkm::Id currentIdx,
270  const vtkm::Id3& incs,
271  vtkm::Id edgeNum,
272  vtkm::UInt8 const* const edgeUses,
273  vtkm::Id* edgeIds,
274  const WholeField& field,
275  const WholeIEdgeField& interpolatedEdgeIds,
276  const WholeWeightField& weights,
277  const WholePointField& points) const
278  {
279  using AxisToSum = SumXAxis;
280 
281  // if this edge is not used then get out
282  if (!edgeUses[edgeNum])
283  {
284  return;
285  }
286  const vtkm::Id writeIndex = edgeIds[edgeNum];
287 
288  // build the edge information
290 
291  vtkm::Id3 offsets1 = data::GetVertOffsets(AxisToSum{}, verts[0]);
292  vtkm::Id3 offsets2 = data::GetVertOffsets(AxisToSum{}, verts[1]);
293 
294  vtkm::Id2 iEdge(currentIdx + vtkm::Dot(offsets1, incs), currentIdx + vtkm::Dot(offsets2, incs));
295 
296  interpolatedEdgeIds.Set(writeIndex, iEdge);
297 
298  auto s0 = field.Get(iEdge[0]);
299  auto s1 = field.Get(iEdge[1]);
300  T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
301  weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
302 
303  auto coord = this->InterpolateCoordinate(t, ijk + offsets1, ijk + offsets2);
304  points.Set(writeIndex, coord);
305  }
306 
307  //----------------------------------------------------------------------------
309  const vtkm::Id3& ijk0,
310  const vtkm::Id3& ijk1) const
311  {
312  return vtkm::Vec3f(
313  this->Origin[0] +
314  this->Spacing[0] *
315  (static_cast<vtkm::FloatDefault>(ijk0[0]) +
316  static_cast<vtkm::FloatDefault>(t) * static_cast<vtkm::FloatDefault>(ijk1[0] - ijk0[0])),
317  this->Origin[1] +
318  this->Spacing[1] *
319  (static_cast<vtkm::FloatDefault>(ijk0[1]) +
320  static_cast<vtkm::FloatDefault>(t) * static_cast<vtkm::FloatDefault>(ijk1[1] - ijk0[1])),
321  this->Origin[2] +
322  this->Spacing[2] *
323  (static_cast<vtkm::FloatDefault>(ijk0[2]) +
324  static_cast<vtkm::FloatDefault>(t) *
325  static_cast<vtkm::FloatDefault>(ijk1[2] - ijk0[2])));
326  }
327 };
328 }
329 }
330 }
331 #endif
vtkm::worklet::flying_edges::generate_tris
VTKM_EXEC void generate_tris(vtkm::Id inputCellId, vtkm::UInt8 edgeCase, vtkm::UInt8 numTris, vtkm::Id *edgeIds, vtkm::Id &triId, const WholeConnField &conn, const WholeCellIdField &cellIds)
Definition: FlyingEdgesPass4Common.h:48
vtkm::worklet::flying_edges::case_includes_axes
VTKM_EXEC bool case_includes_axes(vtkm::UInt8 const *const edgeUses)
Definition: FlyingEdgesPass4Common.h:42
VTKM_EXEC
#define VTKM_EXEC
Definition: ExportMacros.h:51
vtkm::worklet::flying_edges::data::GetEdgeUses
VTKM_EXEC vtkm::UInt8 const * GetEdgeUses(vtkm::UInt8 edgecase)
Definition: FlyingEdgesTables.h:42
vtkm
Groups connected points that have the same field value.
Definition: Atomic.h:19
vtkm::worklet::flying_edges::FlyingEdges3D::MaxBoundary
@ MaxBoundary
Definition: FlyingEdgesHelpers.h:47
vtkm::worklet::flying_edges::compute_incs3d
VTKM_EXEC vtkm::Id3 compute_incs3d(const vtkm::Id3 &dims)
Definition: FlyingEdgesPass4Common.h:26
vtkm::worklet::flying_edges::data::GetVertMap
VTKM_EXEC vtkm::Vec< vtkm::UInt8, 2 > GetVertMap(vtkm::Id index)
Definition: FlyingEdgesTables.h:381
vtkm::worklet::flying_edges::ComputePass4X::CellWriteOffset
vtkm::Id CellWriteOffset
Definition: FlyingEdgesPass4X.h:40
FlyingEdgesHelpers.h
vtkm::worklet::flying_edges::ComputePass4X
Definition: FlyingEdgesPass4X.h:31
vtkm::worklet::flying_edges::ComputePass4X::IsoValue
T IsoValue
Definition: FlyingEdgesPass4X.h:38
vtkm::worklet::flying_edges::ComputePass4X::ComputePass4X
ComputePass4X(T value, const vtkm::Id3 &pdims, const vtkm::Vec3f &origin, const vtkm::Vec3f &spacing, vtkm::Id multiContourCellOffset, vtkm::Id multiContourPointOffset)
Definition: FlyingEdgesPass4X.h:44
vtkm::worklet::flying_edges::ComputePass4X::Generate
VTKM_EXEC void Generate(const vtkm::Vec< vtkm::UInt8, 3 > &boundaryStatus, const vtkm::Id3 &ijk, const WholeDataField &field, const WholeIEdgeField &interpolatedEdgeIds, const WholeWeightField &weights, const WholePointField &points, const vtkm::Id4 &startPos, const vtkm::Id3 &incs, vtkm::Id offset, vtkm::UInt8 const *const edgeUses, vtkm::Id *edgeIds) const
Definition: FlyingEdgesPass4X.h:162
vtkm::worklet::flying_edges::Pass4TrimState
Definition: FlyingEdgesPass4Common.h:120
vtkm::Id
vtkm::Int32 Id
Represents an ID (index into arrays).
Definition: Types.h:191
FlyingEdgesPass4Common.h
vtkm::worklet::flying_edges::SumXAxis
Definition: FlyingEdgesHelpers.h:51
vtkm::worklet::flying_edges::data::GetNumberOfPrimitives
VTKM_EXEC vtkm::UInt8 GetNumberOfPrimitives(vtkm::UInt8 edgecase)
Definition: FlyingEdgesTables.h:26
vtkm::worklet::flying_edges::fully_interior
VTKM_EXEC bool fully_interior(const vtkm::Vec< vtkm::UInt8, 3 > &boundaryStatus)
Definition: FlyingEdgesPass4Common.h:219
vtkm::worklet::flying_edges::data::GetVertOffsets
VTKM_EXEC vtkm::Id3 GetVertOffsets(SumXAxis, vtkm::UInt8 index)
Definition: FlyingEdgesTables.h:391
vtkm::worklet::flying_edges::ComputePass4X::ComputePass4X
ComputePass4X()
Definition: FlyingEdgesPass4X.h:43
vtkm::worklet::flying_edges::ComputePass4X::ExecutionSignature
void(ThreadIndices, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12, WorkIndex) ExecutionSignature
Definition: FlyingEdgesPass4X.h:72
vtkm::worklet::WorkletVisitCellsWithPoints
Base class for worklets that map from Points to Cells.
Definition: WorkletMapTopology.h:255
vtkm::exec::arg::ThreadIndices
The ExecutionSignature tag to use to get the thread indices.
Definition: ThreadIndices.h:41
vtkm::worklet::flying_edges::getEdgeCase
VTKM_EXEC vtkm::UInt8 getEdgeCase(const WholeEdgeField &edges, const vtkm::Id4 &startPos, vtkm::Id inc)
Definition: FlyingEdgesHelpers.h:185
vtkm::worklet::WorkletVisitCellsWithPoints::FieldInPoint
FieldInIncident FieldInPoint
Definition: WorkletMapTopology.h:259
vtkm::worklet::flying_edges::ComputePass4X::PointDims
vtkm::Id3 PointDims
Definition: FlyingEdgesPass4X.h:34
vtkm::UInt8
uint8_t UInt8
Definition: Types.h:157
vtkm::worklet::flying_edges::ComputePass4X::InterpolateCoordinate
VTKM_EXEC vtkm::Vec3f InterpolateCoordinate(T t, const vtkm::Id3 &ijk0, const vtkm::Id3 &ijk1) const
Definition: FlyingEdgesPass4X.h:308
vtkm::Vec3f
vtkm::Vec< vtkm::FloatDefault, 3 > Vec3f
Vec3f corresponds to a 3-dimensional vector of floating point values.
Definition: Types.h:1014
vtkm::worklet::flying_edges::advance_voxelIds
VTKM_EXEC void advance_voxelIds(vtkm::UInt8 const *const edgeUses, vtkm::Id *edgeIds)
Definition: FlyingEdgesPass4Common.h:103
vtkm::Vec< vtkm::Id, 3 >
vtkm::FloatDefault
vtkm::Float32 FloatDefault
The floating point type to use when no other precision is specified.
Definition: Types.h:198
FlyingEdgesTables.h
vtkm::worklet::flying_edges::init_voxelIds
VTKM_EXEC void init_voxelIds(AxisToSum, vtkm::Id writeOffset, vtkm::UInt8 edgeCase, const FieldInPointId3 &axis_sums, vtkm::Id *edgeIds)
Definition: FlyingEdgesPass4Common.h:80
vtkm::worklet::flying_edges::ComputePass4X::Origin
vtkm::Vec3f Origin
Definition: FlyingEdgesPass4X.h:35
WorkletMapTopology.h
vtkm::worklet::flying_edges::ComputePass4X::PointWriteOffset
vtkm::Id PointWriteOffset
Definition: FlyingEdgesPass4X.h:41
vtkm::worklet::flying_edges::ComputePass4X::InterpolateEdge
VTKM_EXEC void InterpolateEdge(const vtkm::Id3 &ijk, vtkm::Id currentIdx, const vtkm::Id3 &incs, vtkm::Id edgeNum, vtkm::UInt8 const *const edgeUses, vtkm::Id *edgeIds, const WholeField &field, const WholeIEdgeField &interpolatedEdgeIds, const WholeWeightField &weights, const WholePointField &points) const
Definition: FlyingEdgesPass4X.h:268
vtkm::worklet::flying_edges::ComputePass4X::ControlSignature
void(CellSetIn, FieldInPoint axis_sums, FieldInPoint axis_mins, FieldInPoint axis_maxs, WholeArrayIn cell_tri_count, WholeArrayIn edgeData, WholeArrayIn data, WholeArrayOut connectivity, WholeArrayOut edgeIds, WholeArrayOut weights, WholeArrayOut inputCellIds, WholeArrayOut points) ControlSignature
Definition: FlyingEdgesPass4X.h:70
vtkm::worklet::flying_edges::ComputePass4X::Spacing
vtkm::Vec3f Spacing
Definition: FlyingEdgesPass4X.h:36
vtkm::worklet::flying_edges::ComputePass4X::operator()
VTKM_EXEC void operator()(const ThreadIndices &threadIndices, const FieldInPointId3 &axis_sums, const FieldInPointId &axis_mins, const FieldInPointId &axis_maxs, const WholeTriField &cellTriCount, const WholeEdgeField &edges, const WholeDataField &field, const WholeConnField &conn, const WholeEdgeIdField &interpolatedEdgeIds, const WholeWeightField &weights, const WholeCellIdField &inputCellIds, const WholePointField &points, vtkm::Id oidx) const
Definition: FlyingEdgesPass4X.h:85
vtkm::exec::arg::WorkIndex
The ExecutionSignature tag to use to get the work index.
Definition: WorkIndex.h:39