From 112024dae217381c273ebfe263b58c4a9cada225 Mon Sep 17 00:00:00 2001 From: Allison Vacanti Date: Thu, 11 Jul 2019 12:41:43 -0400 Subject: [PATCH 1/4] Fix CUDA shfl usage. There was a bug in the implementations of CountSetBits and BitFieldToUnorderedSet. --- .../internal/DeviceAdapterAlgorithmCuda.h | 12 ++++++++-- vtkm/cont/testing/TestingDeviceAdapter.h | 22 +++++++++++++++++++ 2 files changed, 32 insertions(+), 2 deletions(-) diff --git a/vtkm/cont/cuda/internal/DeviceAdapterAlgorithmCuda.h b/vtkm/cont/cuda/internal/DeviceAdapterAlgorithmCuda.h index e6bdb32bf..3a82df3fa 100644 --- a/vtkm/cont/cuda/internal/DeviceAdapterAlgorithmCuda.h +++ b/vtkm/cont/cuda/internal/DeviceAdapterAlgorithmCuda.h @@ -320,7 +320,11 @@ private: vtkm::Int32 rVal = this->LocalPopCount; for (int delta = 1; delta < activeSize; delta *= 2) { - rVal += activeLanes.shfl_down(rVal, delta); + const vtkm::Int32 shflVal = activeLanes.shfl_down(rVal, delta); + if (activeRank + delta < activeSize) + { + rVal += shflVal; + } } if (activeRank == 0) @@ -511,7 +515,11 @@ private: vtkm::Int32 rVal = this->LocalPopCount; for (int delta = 1; delta < activeSize; delta *= 2) { - rVal += activeLanes.shfl_down(rVal, delta); + const vtkm::Int32 shflVal = activeLanes.shfl_down(rVal, delta); + if (activeRank + delta < activeSize) + { + rVal += shflVal; + } } if (activeRank == 0) diff --git a/vtkm/cont/testing/TestingDeviceAdapter.h b/vtkm/cont/testing/TestingDeviceAdapter.h index e6d12ebc5..171fa909e 100644 --- a/vtkm/cont/testing/TestingDeviceAdapter.h +++ b/vtkm/cont/testing/TestingDeviceAdapter.h @@ -2483,6 +2483,17 @@ private: testRandomMask(0xffffffff); testRandomMask(0x1c0fd395); testRandomMask(0xdeadbeef); + + // This case was causing issues on CUDA: + { + BitField bits; + Algorithm::Fill(bits, false, 32 * 32); + auto portal = bits.GetPortalControl(); + portal.SetWord(2, 0x00100000ul); + portal.SetWord(8, 0x00100010ul); + portal.SetWord(11, 0x10000000ul); + testIndexArray(bits); + } } static VTKM_CONT void TestCountSetBits() @@ -2562,6 +2573,17 @@ private: testRandomMask(0xffffffff); testRandomMask(0x1c0fd395); testRandomMask(0xdeadbeef); + + // This case was causing issues on CUDA: + { + BitField bits; + Algorithm::Fill(bits, false, 32 * 32); + auto portal = bits.GetPortalControl(); + portal.SetWord(2, 0x00100000ul); + portal.SetWord(8, 0x00100010ul); + portal.SetWord(11, 0x10000000ul); + verifyPopCount(bits); + } } template From ab627b61d6061b81dad352759b20c557399779e9 Mon Sep 17 00:00:00 2001 From: Allison Vacanti Date: Thu, 24 Jan 2019 11:09:21 -0500 Subject: [PATCH 2/4] Add OrientNormals worklet. --- vtkm/worklet/CMakeLists.txt | 4 + vtkm/worklet/OrientCellNormals.h | 451 ++++++++++++++++++ vtkm/worklet/OrientNormals.h | 114 +++++ vtkm/worklet/OrientPointAndCellNormals.h | 421 ++++++++++++++++ vtkm/worklet/OrientPointNormals.h | 369 ++++++++++++++ vtkm/worklet/testing/CMakeLists.txt | 1 + .../worklet/testing/UnitTestOrientNormals.cxx | 408 ++++++++++++++++ 7 files changed, 1768 insertions(+) create mode 100644 vtkm/worklet/OrientCellNormals.h create mode 100644 vtkm/worklet/OrientNormals.h create mode 100644 vtkm/worklet/OrientPointAndCellNormals.h create mode 100644 vtkm/worklet/OrientPointNormals.h create mode 100644 vtkm/worklet/testing/UnitTestOrientNormals.cxx diff --git a/vtkm/worklet/CMakeLists.txt b/vtkm/worklet/CMakeLists.txt index 500252a03..a23e9ea76 100644 --- a/vtkm/worklet/CMakeLists.txt +++ b/vtkm/worklet/CMakeLists.txt @@ -49,6 +49,10 @@ set(headers NDimsHistMarginalization.h NDimsHistogram.h Normalize.h + OrientCellNormals.h + OrientNormals.h + OrientPointNormals.h + OrientPointAndCellNormals.h OscillatorSource.h ParticleAdvection.h PointAverage.h diff --git a/vtkm/worklet/OrientCellNormals.h b/vtkm/worklet/OrientCellNormals.h new file mode 100644 index 000000000..deac6f10b --- /dev/null +++ b/vtkm/worklet/OrientCellNormals.h @@ -0,0 +1,451 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +// +// Copyright 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2019 UT-Battelle, LLC. +// Copyright 2019 Los Alamos National Security. +// +// Under the terms of Contract DE-NA0003525 with NTESS, +// the U.S. Government retains certain rights in this software. +// +// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National +// Laboratory (LANL), the U.S. Government retains certain rights in +// this software. +//============================================================================ +#ifndef vtkm_m_worklet_OrientCellNormals_h +#define vtkm_m_worklet_OrientCellNormals_h + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +namespace vtkm +{ +namespace worklet +{ + +/// +/// Orients normals to point outside of the dataset. This requires a closed +/// manifold surface or else the behavior is undefined. This requires an +/// unstructured cellset as input. +/// +class OrientCellNormals +{ + static constexpr vtkm::Id INVALID_ID = -1; + + // Returns true if v1 and v2 are pointing in the same hemisphere. + template + VTKM_EXEC static bool SameDirection(const vtkm::Vec& v1, const vtkm::Vec& v2) + { + return vtkm::Dot(v1, v2) >= 0; + } + + // Ensure that the normal is pointing in the same hemisphere as ref. + // Returns true if normal is modified. + template + VTKM_EXEC static bool Align(vtkm::Vec& normal, const vtkm::Vec& ref) + { + if (!SameDirection(normal, ref)) + { + normal = -normal; + return true; + } + return false; + } + +public: + // Locates starting points for BFS traversal of dataset by finding points + // on the dataset boundaries. These points are marked as active. + // Initializes the ActivePoints array. + class WorkletMarkSourcePoints : public vtkm::worklet::WorkletMapField + { + public: + using ControlSignature = void(FieldIn coords, WholeArrayIn ranges, FieldOut activePoints); + using ExecutionSignature = _3(_1 coord, _2 ranges); + + template + VTKM_EXEC bool operator()(const vtkm::Vec& point, const RangePortal& ranges) const + { + bool isActive = false; + for (vtkm::IdComponent dim = 0; dim < 3; ++dim) + { + const auto& range = ranges.Get(dim); + const auto val = static_cast(point[dim]); + if (val <= range.Min || val >= range.Max) + { + isActive = true; + } + } + return isActive; + } + }; + + // For each of the source points, determine the boundaries it lies on. Align + // each incident cell's normal to point out of the boundary, marking each cell + // as both visited and active. + // Clears the active flags for points, and marks the current point as visited. + class WorkletProcessSourceCells : public vtkm::worklet::WorkletMapCellToPoint + { + public: + using ControlSignature = void(CellSetIn cells, + FieldInPoint coords, + WholeArrayIn ranges, + WholeArrayInOut cellNormals, + // InOut for preserve data on masked indices + BitFieldInOut activeCells, + BitFieldInOut visitedCells, + FieldInOutPoint activePoints, + FieldInOutPoint visitedPoints); + using ExecutionSignature = void(CellIndices cellIds, + _2 coords, + _3 ranges, + _4 cellNormals, + _5 activeCells, + _6 visitedCells, + _7 activePoints, + _8 visitedPoints); + + using MaskType = vtkm::worklet::MaskIndices; + + template + VTKM_EXEC void operator()(const CellList& cellIds, + const vtkm::Vec& coord, + const RangePortal& ranges, + CellNormalPortal& cellNormals, + ActiveCellsBitPortal& activeCells, + VisitedCellsBitPortal& visitedCells, + bool& pointIsActive, + bool& pointIsVisited) const + { + using NormalType = typename CellNormalPortal::ValueType; + VTKM_STATIC_ASSERT_MSG(vtkm::VecTraits::NUM_COMPONENTS == 3, + "Cell Normals expected to have 3 components."); + using NormalCompType = typename NormalType::ComponentType; + + // Find the vector that points out of the dataset from the current point: + const NormalType refNormal = [&]() -> NormalType { + NormalType normal{ NormalCompType{ 0 } }; + NormalCompType numNormals{ 0 }; + for (vtkm::IdComponent dim = 0; dim < 3; ++dim) + { + const auto range = ranges.Get(dim); + const auto val = static_cast(coord[dim]); + if (val <= range.Min) + { + NormalType compNormal{ NormalCompType{ 0 } }; + compNormal[dim] = NormalCompType{ -1 }; + normal += compNormal; + ++numNormals; + } + else if (val >= range.Max) + { + NormalType compNormal{ NormalCompType{ 0 } }; + compNormal[dim] = NormalCompType{ 1 }; + normal += compNormal; + ++numNormals; + } + } + + VTKM_ASSERT("Source point is not on a boundary?" && numNormals > 0.5); + return normal / numNormals; + }(); + + // Align all cell normals to the reference, marking them as active and + // visited. + const vtkm::IdComponent numCells = cellIds.GetNumberOfComponents(); + for (vtkm::IdComponent c = 0; c < numCells; ++c) + { + const vtkm::Id cellId = cellIds[c]; + + if (!visitedCells.OrBitAtomic(cellId, true)) + { // This thread is the first to touch this cell. + activeCells.SetBitAtomic(cellId, true); + + NormalType cellNormal = cellNormals.Get(cellId); + if (Align(cellNormal, refNormal)) + { + cellNormals.Set(cellId, cellNormal); + } + } + } + + // Mark current point as inactive but visited: + pointIsActive = false; + pointIsVisited = true; + } + }; + + // Mark each incident point as active and visited. + // Marks the current cell as inactive. + class WorkletMarkActivePoints : public vtkm::worklet::WorkletMapPointToCell + { + public: + using ControlSignature = void(CellSetIn cell, + BitFieldInOut activePoints, + BitFieldInOut visitedPoints, + FieldInOutCell activeCells); + using ExecutionSignature = _4(PointIndices pointIds, _2 activePoints, _3 visitedPoints); + + using MaskType = vtkm::worklet::MaskIndices; + + template + VTKM_EXEC bool operator()(const PointList& pointIds, + ActivePointsBitPortal& activePoints, + VisitedPointsBitPortal& visitedPoints) const + { + const vtkm::IdComponent numPoints = pointIds.GetNumberOfComponents(); + for (vtkm::IdComponent p = 0; p < numPoints; ++p) + { + const vtkm::Id pointId = pointIds[p]; + if (!visitedPoints.OrBitAtomic(pointId, true)) + { // This thread owns this point. + activePoints.SetBitAtomic(pointId, true); + } + } + + // Mark current cell as inactive: + return false; + } + }; + + // Mark each incident cell as active, setting a visited neighbor + // cell as its reference for alignment. + // Marks the current point as inactive. + class WorkletMarkActiveCells : public vtkm::worklet::WorkletMapCellToPoint + { + public: + using ControlSignature = void(CellSetIn cells, + WholeArrayOut refCells, + BitFieldInOut activeCells, + BitFieldIn visitedCells, + FieldInOutPoint activePoints); + using ExecutionSignature = _5(CellIndices cellIds, + _2 refCells, + _3 activeCells, + _4 visitedCells); + + using MaskType = vtkm::worklet::MaskIndices; + + template + VTKM_EXEC bool operator()(const CellList& cellIds, + RefCellPortal& refCells, + ActiveCellBitPortal& activeCells, + const VisitedCellBitPortal& visitedCells) const + { + // One of the cells must be marked visited already. Find it and use it as + // an alignment reference for the others: + const vtkm::IdComponent numCells = cellIds.GetNumberOfComponents(); + const vtkm::Id refCellId = [&]() -> vtkm::Id { + for (vtkm::IdComponent c = 0; c < numCells; ++c) + { + const vtkm::Id cellId = cellIds[c]; + if (visitedCells.GetBit(cellId)) + { + return cellId; + } + } + return INVALID_ID; + }(); + + VTKM_ASSERT("No reference cell found." && refCellId != INVALID_ID); + + for (vtkm::IdComponent c = 0; c < numCells; ++c) + { + const vtkm::Id cellId = cellIds[c]; + if (!visitedCells.GetBit(cellId)) + { + if (!activeCells.OrBitAtomic(cellId, true)) + { // This thread owns this cell. + refCells.Set(cellId, refCellId); + } + } + } + + // Mark current point as inactive: + return false; + } + }; + + // Align the normal of each active cell, to its reference cell normal. + // The cell is marked visited. + class WorkletProcessCellNormals : public vtkm::worklet::WorkletMapField + { + public: + using ControlSignature = void(FieldIn refCells, + WholeArrayInOut cellNormals, + FieldInOut visitedCells); + using ExecutionSignature = _3(InputIndex cellId, _1 refCellId, _2 cellNormals); + + using MaskType = vtkm::worklet::MaskIndices; + + template + VTKM_EXEC bool operator()(const vtkm::Id cellId, + const vtkm::Id refCellId, + CellNormalsPortal& cellNormals) const + { + const auto refNormal = cellNormals.Get(refCellId); + auto normal = cellNormals.Get(cellId); + if (Align(normal, refNormal)) + { + cellNormals.Set(cellId, normal); + } + + // Mark cell as visited: + return true; + } + }; + + template + VTKM_CONT static void Run( + const CellSetType& cells, + const vtkm::cont::ArrayHandle, CoordsStorageType>& coords, + vtkm::cont::ArrayHandle, CellNormalStorageType>& cellNormals) + { + using RangeType = vtkm::cont::ArrayHandle; + + using MarkSourcePoints = vtkm::worklet::DispatcherMapField; + using ProcessSourceCells = vtkm::worklet::DispatcherMapTopology; + using MarkActivePoints = vtkm::worklet::DispatcherMapTopology; + using MarkActiveCells = vtkm::worklet::DispatcherMapTopology; + using ProcessCellNormals = vtkm::worklet::DispatcherMapField; + + const vtkm::Id numPoints = coords.GetNumberOfValues(); + const vtkm::Id numCells = cells.GetNumberOfCells(); + + VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf, + "OrientCellNormals worklet (%lld points, %lld cells)", + static_cast(numPoints), + static_cast(numCells)); + + // active = cells / point to be used in the next worklet invocation mask. + vtkm::cont::BitField activePointBits; // Initialized by MarkSourcePoints + auto activePoints = vtkm::cont::make_ArrayHandleBitField(activePointBits); + + vtkm::cont::BitField activeCellBits; + vtkm::cont::Algorithm::Fill(activeCellBits, false, numCells); + auto activeCells = vtkm::cont::make_ArrayHandleBitField(activeCellBits); + + // visited = cells / points that have been corrected. + vtkm::cont::BitField visitedPointBits; + vtkm::cont::Algorithm::Fill(visitedPointBits, false, numPoints); + auto visitedPoints = vtkm::cont::make_ArrayHandleBitField(visitedPointBits); + + vtkm::cont::BitField visitedCellBits; + vtkm::cont::Algorithm::Fill(visitedCellBits, false, numCells); + auto visitedCells = vtkm::cont::make_ArrayHandleBitField(visitedCellBits); + + vtkm::cont::ArrayHandle mask; // Allocated as needed + + // For each cell, store a reference alignment cell. + vtkm::cont::ArrayHandle refCells; + { + vtkm::cont::Algorithm::Copy( + vtkm::cont::make_ArrayHandleConstant(INVALID_ID, numCells), refCells); + } + + // 1) Compute range of coords. + const RangeType ranges = vtkm::cont::ArrayRangeCompute(coords); + + // 2) Locate points on a boundary, since their normal alignment direction + // is known. + { + MarkSourcePoints dispatcher; + dispatcher.Invoke(coords, ranges, activePoints); + } + + // 3) For each source point, align the normals of the adjacent cells. + { + vtkm::Id numActive = vtkm::cont::Algorithm::BitFieldToUnorderedSet(activePointBits, mask); + (void)numActive; + VTKM_LOG_S(vtkm::cont::LogLevel::Perf, + "ProcessSourceCells from " << numActive << " source points."); + ProcessSourceCells dispatcher{ vtkm::worklet::MaskIndices{ mask } }; + dispatcher.Invoke(cells, + coords, + ranges, + cellNormals, + activeCellBits, + visitedCellBits, + activePoints, + visitedPoints); + } + + for (size_t iter = 1;; ++iter) + { + // 4) Mark unvisited points adjacent to active cells + { + vtkm::Id numActive = vtkm::cont::Algorithm::BitFieldToUnorderedSet(activeCellBits, mask); + (void)numActive; + VTKM_LOG_S(vtkm::cont::LogLevel::Perf, + "MarkActivePoints from " << numActive << " active cells."); + MarkActivePoints dispatcher{ vtkm::worklet::MaskIndices{ mask } }; + dispatcher.Invoke(cells, activePointBits, visitedPointBits, activeCells); + } + + // 5) Mark unvisited cells adjacent to active points + { + vtkm::Id numActive = vtkm::cont::Algorithm::BitFieldToUnorderedSet(activePointBits, mask); + (void)numActive; + VTKM_LOG_S(vtkm::cont::LogLevel::Perf, + "MarkActiveCells from " << numActive << " active points."); + MarkActiveCells dispatcher{ vtkm::worklet::MaskIndices{ mask } }; + dispatcher.Invoke(cells, refCells, activeCellBits, visitedCellBits, activePoints); + } + + vtkm::Id numActiveCells = vtkm::cont::Algorithm::BitFieldToUnorderedSet(activeCellBits, mask); + + if (numActiveCells == 0) + { // Done! + VTKM_LOG_S(vtkm::cont::LogLevel::Perf, "Iteration " << iter << ": Traversal complete."); + break; + } + + VTKM_LOG_S(vtkm::cont::LogLevel::Perf, + "Iteration " << iter << ": Processing " << numActiveCells << " normals."); + + // 5) Correct normals for active cells. + { + ProcessCellNormals dispatcher{ vtkm::worklet::MaskIndices{ mask } }; + dispatcher.Invoke(refCells, cellNormals, visitedCells); + } + } + } +}; +} +} // end namespace vtkm::worklet + + +#endif // vtkm_m_worklet_OrientCellNormals_h diff --git a/vtkm/worklet/OrientNormals.h b/vtkm/worklet/OrientNormals.h new file mode 100644 index 000000000..f15f19b1f --- /dev/null +++ b/vtkm/worklet/OrientNormals.h @@ -0,0 +1,114 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +// +// Copyright 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2019 UT-Battelle, LLC. +// Copyright 2019 Los Alamos National Security. +// +// Under the terms of Contract DE-NA0003525 with NTESS, +// the U.S. Government retains certain rights in this software. +// +// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National +// Laboratory (LANL), the U.S. Government retains certain rights in +// this software. +//============================================================================ +#ifndef vtkm_m_worklet_OrientNormals_h +#define vtkm_m_worklet_OrientNormals_h + +#include + +#include +#include +#include + +#include +#include +#include + +namespace vtkm +{ +namespace worklet +{ + +/// +/// Orients normals to point outside of the dataset. This requires a closed +/// manifold surface or else the behavior is undefined. This requires an +/// unstructured cellset as input. +/// +class OrientNormals +{ +public: + template + VTKM_CONT static void RunCellNormals( + const CellSetType& cells, + const vtkm::cont::ArrayHandle, CoordsStorageType>& coords, + vtkm::cont::ArrayHandle, CellNormalStorageType>& cellNormals) + { + OrientCellNormals::Run(cells, coords, cellNormals); + } + + template + VTKM_CONT static void RunPointNormals( + const CellSetType& cells, + const vtkm::cont::ArrayHandle, CoordsStorageType>& coords, + vtkm::cont::ArrayHandle, PointNormalStorageType>& + pointNormals) + { + OrientPointNormals::Run(cells, coords, pointNormals); + } + + template + VTKM_CONT static void RunPointAndCellNormals( + const CellSetType& cells, + const vtkm::cont::ArrayHandle, CoordsStorageType>& coords, + vtkm::cont::ArrayHandle, PointNormalStorageType>& + pointNormals, + vtkm::cont::ArrayHandle, CellNormalStorageType>& cellNormals) + { + OrientPointAndCellNormals::Run(cells, coords, pointNormals, cellNormals); + } + + struct NegateFunctor + { + template + VTKM_EXEC_CONT T operator()(const T& val) const + { + return -val; + } + }; + + /// + /// Reverse the normals to point in the opposite direction. + /// + template + VTKM_CONT static void RunFlipNormals( + vtkm::cont::ArrayHandle, NormalStorageType>& normals) + { + const auto flippedAlias = vtkm::cont::make_ArrayHandleTransform(normals, NegateFunctor{}); + vtkm::cont::Algorithm::Copy(flippedAlias, normals); + } +}; +} +} // end namespace vtkm::worklet + + +#endif // vtkm_m_worklet_OrientNormals_h diff --git a/vtkm/worklet/OrientPointAndCellNormals.h b/vtkm/worklet/OrientPointAndCellNormals.h new file mode 100644 index 000000000..613166a2f --- /dev/null +++ b/vtkm/worklet/OrientPointAndCellNormals.h @@ -0,0 +1,421 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +// +// Copyright 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2019 UT-Battelle, LLC. +// Copyright 2019 Los Alamos National Security. +// +// Under the terms of Contract DE-NA0003525 with NTESS, +// the U.S. Government retains certain rights in this software. +// +// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National +// Laboratory (LANL), the U.S. Government retains certain rights in +// this software. +//============================================================================ +#ifndef vtkm_m_worklet_OrientPointAndCellNormals_h +#define vtkm_m_worklet_OrientPointAndCellNormals_h + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +namespace vtkm +{ +namespace worklet +{ + +/// +/// Orients normals to point outside of the dataset. This requires a closed +/// manifold surface or else the behavior is undefined. This requires an +/// unstructured cellset as input. +/// +class OrientPointAndCellNormals +{ + static constexpr vtkm::Id INVALID_ID = -1; + + // Returns true if v1 and v2 are pointing in the same hemisphere. + template + VTKM_EXEC static bool SameDirection(const vtkm::Vec& v1, const vtkm::Vec& v2) + { + return vtkm::Dot(v1, v2) >= 0; + } + + // Ensure that the normal is pointing in the same hemisphere as ref. + // Returns true if normal is modified. + template + VTKM_EXEC static bool Align(vtkm::Vec& normal, const vtkm::Vec& ref) + { + if (!SameDirection(normal, ref)) + { + normal = -normal; + return true; + } + return false; + } + +public: + // Locates starting points for BFS traversal of dataset by finding points + // on the dataset boundaries. The normals for these points are corrected by + // making them point outside of the dataset, and they are marked as both + // active and visited. + class WorkletMarkSourcePoints : public vtkm::worklet::WorkletMapField + { + public: + using ControlSignature = void(FieldIn coords, + FieldInOut pointNormals, + WholeArrayIn ranges, + FieldOut activePoints, + FieldOut visitedPoints); + using ExecutionSignature = + void(_1 coord, _2 pointNormal, _3 ranges, _4 activePoints, _5 visitedPoints); + + template + VTKM_EXEC void operator()(const vtkm::Vec& point, + vtkm::Vec& pointNormal, + const RangePortal& ranges, + bool& isActive, + bool& isVisited) const + { + for (vtkm::IdComponent dim = 0; dim < 3; ++dim) + { + const auto& range = ranges.Get(dim); + const auto val = static_cast(point[dim]); + if (val <= range.Min) + { + vtkm::Vec ref{ static_cast(0) }; + ref[dim] = static_cast(-1); + Align(pointNormal, ref); + isActive = true; + isVisited = true; + return; + } + else if (val >= range.Max) + { + vtkm::Vec ref{ static_cast(0) }; + ref[dim] = static_cast(1); + Align(pointNormal, ref); + isActive = true; + isVisited = true; + return; + } + } + + isActive = false; + isVisited = false; + } + }; + + // Mark each incident cell as active and visited. + // Marks the current point as inactive. + class WorkletMarkActiveCells : public vtkm::worklet::WorkletMapCellToPoint + { + public: + using ControlSignature = void(CellSetIn cell, + BitFieldInOut activeCells, + BitFieldInOut visitedCells, + FieldInOutPoint activePoints); + using ExecutionSignature = _4(CellIndices cellIds, _2 activeCells, _3 visitedCells); + + using MaskType = vtkm::worklet::MaskIndices; + + template + VTKM_EXEC bool operator()(const CellList& cellIds, + ActiveCellsBitPortal& activeCells, + VisitedCellsBitPortal& visitedCells) const + { + const vtkm::IdComponent numCells = cellIds.GetNumberOfComponents(); + for (vtkm::IdComponent c = 0; c < numCells; ++c) + { + const vtkm::Id cellId = cellIds[c]; + if (!visitedCells.OrBitAtomic(cellId, true)) + { // This thread owns this cell. + activeCells.SetBitAtomic(cellId, true); + } + } + + // Mark current point as inactive: + return false; + } + }; + + // Align the current cell's normals to an adjacent visited point's normal. + class WorkletProcessCellNormals : public vtkm::worklet::WorkletMapPointToCell + { + public: + using ControlSignature = void(CellSetIn cells, + WholeArrayIn pointNormals, + WholeArrayInOut cellNormals, + BitFieldIn visitedPoints); + using ExecutionSignature = void(PointIndices pointIds, + InputIndex cellId, + _2 pointNormals, + _3 cellNormals, + _4 visitedPoints); + + using MaskType = vtkm::worklet::MaskIndices; + + template + VTKM_EXEC void operator()(const PointList& pointIds, + const vtkm::Id cellId, + const PointNormalsPortal& pointNormals, + CellNormalsPortal& cellNormals, + const VisitedPointsBitPortal& visitedPoints) const + { + // Use the normal of a visited point as a reference: + const vtkm::Id refPointId = [&]() -> vtkm::Id { + const vtkm::IdComponent numPoints = pointIds.GetNumberOfComponents(); + for (vtkm::IdComponent p = 0; p < numPoints; ++p) + { + const vtkm::Id pointId = pointIds[p]; + if (visitedPoints.GetBit(pointId)) + { + return pointId; + } + } + + return INVALID_ID; + }(); + + VTKM_ASSERT("No reference point found." && refPointId != INVALID_ID); + + const auto refNormal = pointNormals.Get(refPointId); + auto normal = cellNormals.Get(cellId); + if (Align(normal, refNormal)) + { + cellNormals.Set(cellId, normal); + } + } + }; + + // Mark each incident point as active and visited. + // Marks the current cell as inactive. + class WorkletMarkActivePoints : public vtkm::worklet::WorkletMapPointToCell + { + public: + using ControlSignature = void(CellSetIn cell, + BitFieldInOut activePoints, + BitFieldInOut visitedPoints, + FieldInOutCell activeCells); + using ExecutionSignature = _4(PointIndices pointIds, _2 activePoint, _3 visitedPoint); + + using MaskType = vtkm::worklet::MaskIndices; + + template + VTKM_EXEC bool operator()(const PointList& pointIds, + ActivePointsBitPortal& activePoints, + VisitedPointsBitPortal& visitedPoints) const + { + const vtkm::IdComponent numPoints = pointIds.GetNumberOfComponents(); + for (vtkm::IdComponent p = 0; p < numPoints; ++p) + { + const vtkm::Id pointId = pointIds[p]; + if (!visitedPoints.OrBitAtomic(pointId, true)) + { // This thread owns this point. + activePoints.SetBitAtomic(pointId, true); + } + } + + // Mark current cell as inactive: + return false; + } + }; + + // Align the current point's normals to an adjacent visited cell's normal. + class WorkletProcessPointNormals : public vtkm::worklet::WorkletMapCellToPoint + { + public: + using ControlSignature = void(CellSetIn cells, + WholeArrayInOut pointNormals, + WholeArrayIn cellNormals, + BitFieldIn visitedCells); + using ExecutionSignature = void(CellIndices cellIds, + InputIndex pointId, + _2 pointNormals, + _3 cellNormals, + _4 visitedCells); + + using MaskType = vtkm::worklet::MaskIndices; + + template + VTKM_EXEC void operator()(const CellList& cellIds, + const vtkm::Id pointId, + PointNormalsPortal& pointNormals, + const CellNormalsPortal& cellNormals, + const VisitedCellsBitPortal& visitedCells) const + { + // Use the normal of a visited cell as a reference: + const vtkm::Id refCellId = [&]() -> vtkm::Id { + const vtkm::IdComponent numCells = cellIds.GetNumberOfComponents(); + for (vtkm::IdComponent c = 0; c < numCells; ++c) + { + const vtkm::Id cellId = cellIds[c]; + if (visitedCells.GetBit(cellId)) + { + return cellId; + } + } + + return INVALID_ID; + }(); + + VTKM_ASSERT("No reference cell found." && refCellId != INVALID_ID); + + const auto refNormal = cellNormals.Get(refCellId); + auto normal = pointNormals.Get(pointId); + if (Align(normal, refNormal)) + { + pointNormals.Set(pointId, normal); + } + } + }; + + template + VTKM_CONT static void Run( + const CellSetType& cells, + const vtkm::cont::ArrayHandle, CoordsStorageType>& coords, + vtkm::cont::ArrayHandle, PointNormalStorageType>& + pointNormals, + vtkm::cont::ArrayHandle, CellNormalStorageType>& cellNormals) + { + using RangeType = vtkm::cont::ArrayHandle; + + using MarkSourcePoints = vtkm::worklet::DispatcherMapField; + using MarkActiveCells = vtkm::worklet::DispatcherMapTopology; + using ProcessCellNormals = vtkm::worklet::DispatcherMapTopology; + using MarkActivePoints = vtkm::worklet::DispatcherMapTopology; + using ProcessPointNormals = vtkm::worklet::DispatcherMapTopology; + + const vtkm::Id numCells = cells.GetNumberOfCells(); + + VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf, + "OrientPointAndCellNormals worklet (%lld points, %lld cells)", + static_cast(coords.GetNumberOfValues()), + static_cast(numCells)); + + // active = cells / point to be used in the next worklet invocation mask. + vtkm::cont::BitField activePointBits; // Initialized by MarkSourcePoints + auto activePoints = vtkm::cont::make_ArrayHandleBitField(activePointBits); + + vtkm::cont::BitField activeCellBits; + vtkm::cont::Algorithm::Fill(activeCellBits, false, numCells); + auto activeCells = vtkm::cont::make_ArrayHandleBitField(activeCellBits); + + // visited = cells / points that have been corrected. + vtkm::cont::BitField visitedPointBits; // Initialized by MarkSourcePoints + auto visitedPoints = vtkm::cont::make_ArrayHandleBitField(visitedPointBits); + + vtkm::cont::BitField visitedCellBits; + vtkm::cont::Algorithm::Fill(visitedCellBits, false, numCells); + auto visitedCells = vtkm::cont::make_ArrayHandleBitField(visitedCellBits); + + vtkm::cont::ArrayHandle mask; // Allocated as needed + + // 1) Compute range of coords. + const RangeType ranges = vtkm::cont::ArrayRangeCompute(coords); + + // 2) Locate points on a boundary and align their normal to point out of the + // dataset: + { + MarkSourcePoints dispatcher; + dispatcher.Invoke(coords, pointNormals, ranges, activePoints, visitedPoints); + } + + for (size_t iter = 1;; ++iter) + { + // 3) Mark unvisited cells adjacent to active points + { + vtkm::Id numActive = vtkm::cont::Algorithm::BitFieldToUnorderedSet(activePointBits, mask); + (void)numActive; + VTKM_LOG_S(vtkm::cont::LogLevel::Perf, + "MarkActiveCells from " << numActive << " active points."); + MarkActiveCells dispatcher{ vtkm::worklet::MaskIndices{ mask } }; + dispatcher.Invoke(cells, activeCellBits, visitedCellBits, activePoints); + } + + vtkm::Id numActiveCells = vtkm::cont::Algorithm::BitFieldToUnorderedSet(activeCellBits, mask); + + if (numActiveCells == 0) + { // Done! + VTKM_LOG_S(vtkm::cont::LogLevel::Perf, + "Iteration " << iter << ": Traversal complete; no more cells"); + break; + } + + VTKM_LOG_S(vtkm::cont::LogLevel::Perf, + "Iteration " << iter << ": Processing " << numActiveCells << " cell normals."); + + // 4) Correct normals for active cells. + { + ProcessCellNormals dispatcher{ vtkm::worklet::MaskIndices{ mask } }; + dispatcher.Invoke(cells, pointNormals, cellNormals, visitedPointBits); + } + + // 5) Mark unvisited points adjacent to active cells + { + vtkm::Id numActive = vtkm::cont::Algorithm::BitFieldToUnorderedSet(activeCellBits, mask); + (void)numActive; + VTKM_LOG_S(vtkm::cont::LogLevel::Perf, + "MarkActivePoints from " << numActive << " active cells."); + MarkActivePoints dispatcher{ vtkm::worklet::MaskIndices{ mask } }; + dispatcher.Invoke(cells, activePointBits, visitedPointBits, activeCells); + } + + vtkm::Id numActivePoints = + vtkm::cont::Algorithm::BitFieldToUnorderedSet(activePointBits, mask); + + if (numActivePoints == 0) + { // Done! + VTKM_LOG_S(vtkm::cont::LogLevel::Perf, + "Iteration " << iter << ": Traversal complete; no more points"); + break; + } + + VTKM_LOG_S(vtkm::cont::LogLevel::Perf, + "Iteration " << iter << ": Processing " << numActivePoints << " point normals."); + + // 4) Correct normals for active points. + { + ProcessPointNormals dispatcher{ vtkm::worklet::MaskIndices{ mask } }; + dispatcher.Invoke(cells, pointNormals, cellNormals, visitedCellBits); + } + } + } +}; +} +} // end namespace vtkm::worklet + + +#endif // vtkm_m_worklet_OrientPointAndCellNormals_h diff --git a/vtkm/worklet/OrientPointNormals.h b/vtkm/worklet/OrientPointNormals.h new file mode 100644 index 000000000..3704f4f79 --- /dev/null +++ b/vtkm/worklet/OrientPointNormals.h @@ -0,0 +1,369 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +// +// Copyright 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2019 UT-Battelle, LLC. +// Copyright 2019 Los Alamos National Security. +// +// Under the terms of Contract DE-NA0003525 with NTESS, +// the U.S. Government retains certain rights in this software. +// +// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National +// Laboratory (LANL), the U.S. Government retains certain rights in +// this software. +//============================================================================ +#ifndef vtkm_m_worklet_OrientPointNormals_h +#define vtkm_m_worklet_OrientPointNormals_h + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace vtkm +{ +namespace worklet +{ + +/// +/// Orients normals to point outside of the dataset. This requires a closed +/// manifold surface or else the behavior is undefined. This requires an +/// unstructured cellset as input. +/// +class OrientPointNormals +{ + static constexpr vtkm::Id INVALID_ID = -1; + + // Returns true if v1 and v2 are pointing in the same hemisphere. + template + VTKM_EXEC static bool SameDirection(const vtkm::Vec& v1, const vtkm::Vec& v2) + { + return vtkm::Dot(v1, v2) >= 0; + } + + // Ensure that the normal is pointing in the same hemisphere as ref. + // Returns true if normal is modified. + template + VTKM_EXEC static bool Align(vtkm::Vec& normal, const vtkm::Vec& ref) + { + if (!SameDirection(normal, ref)) + { + normal = -normal; + return true; + } + return false; + } + +public: + // Locates starting points for BFS traversal of dataset by finding points + // on the dataset boundaries. The normals for these points are corrected by + // making them point outside of the dataset, and they are marked as both + // active and visited. + class WorkletMarkSourcePoints : public vtkm::worklet::WorkletMapField + { + public: + using ControlSignature = void(FieldIn coords, + FieldInOut normals, + WholeArrayIn ranges, + FieldOut activePoints, + FieldOut visitedPoints, + FieldOut refPoints); + using ExecutionSignature = + _6(InputIndex pointId, _1 coord, _2 normal, _3 ranges, _4 activePoints, _5 visitedPoints); + + template + VTKM_EXEC vtkm::Id operator()(const vtkm::Id pointId, + const vtkm::Vec& point, + vtkm::Vec& normal, + const RangePortal& ranges, + bool& isActive, + bool& isVisited) const + { + for (vtkm::IdComponent dim = 0; dim < 3; ++dim) + { + const auto& range = ranges.Get(dim); + const auto val = static_cast(point[dim]); + if (val <= range.Min) + { + vtkm::Vec ref{ static_cast(0) }; + ref[dim] = static_cast(-1); + Align(normal, ref); + isActive = true; + isVisited = true; + return pointId; + } + else if (val >= range.Max) + { + vtkm::Vec ref{ static_cast(0) }; + ref[dim] = static_cast(1); + Align(normal, ref); + isActive = true; + isVisited = true; + return pointId; + } + } + + isActive = false; + isVisited = false; + return INVALID_ID; + } + }; + + // Traverses the active points (via mask) and marks the connected cells as + // active. Set the reference point for all adjacent cells to the current + // point. + class WorkletMarkActiveCells : public vtkm::worklet::WorkletMapCellToPoint + { + public: + using ControlSignature = void(CellSetIn cellSet, + // InOut to preserve data on masked indices + BitFieldInOut activeCells, + BitFieldInOut visitedCells, + FieldInOutPoint activePoints); + using ExecutionSignature = _4(CellIndices cells, _2 activeCells, _3 visitedCells); + + using MaskType = vtkm::worklet::MaskIndices; + + // Mark all unvisited cells as active: + template + VTKM_EXEC bool operator()(const CellListT& cells, + ActiveCellsT& activeCells, + VisitedCellsT& visitedCells) const + { + for (vtkm::IdComponent c = 0; c < cells.GetNumberOfComponents(); ++c) + { + const vtkm::Id cellId = cells[c]; + const bool alreadyVisited = visitedCells.CompareAndSwapBitAtomic(cellId, true, false); + + if (!alreadyVisited) + { // This thread is first to visit cell + activeCells.SetBitAtomic(cellId, true); + } + } + + // Mark the current point as inactive: + return false; + } + }; + + // Traverses the active cells and mark the connected points as active, + // propogating the reference pointId. + class WorkletMarkActivePoints : public vtkm::worklet::WorkletMapPointToCell + { + public: + using ControlSignature = void(CellSetIn cellSet, + BitFieldInOut activePoints, + BitFieldIn visitedPoints, + WholeArrayInOut refPoints, + FieldInOutCell activeCells); + using ExecutionSignature = _5(PointIndices points, + _2 activePoints, + _3 visitedPoints, + _4 refPoints); + + using MaskType = vtkm::worklet::MaskIndices; + + template + VTKM_EXEC bool operator()(const PointListT& points, + ActivePointsT& activePoints, + VisitedPointsT& visitedPoints, + RefPointsT& refPoints) const + { + // Find any point in the cell that has already been visited, and take + // its id as the reference for this cell. + vtkm::Id refPtId = INVALID_ID; + for (vtkm::IdComponent p = 0; p < points.GetNumberOfComponents(); ++p) + { + const vtkm::Id pointId = points[p]; + const bool alreadyVisited = visitedPoints.GetBit(pointId); + if (alreadyVisited) + { + refPtId = pointId; + break; + } + } + + // There must be one valid point in each cell: + VTKM_ASSERT("Reference point not found." && refPtId != INVALID_ID); + + // Propogate the reference point to other cell members + for (vtkm::IdComponent p = 0; p < points.GetNumberOfComponents(); ++p) + { + const vtkm::Id pointId = points[p]; + + // Mark this point as active + const bool alreadyVisited = visitedPoints.GetBit(pointId); + if (!alreadyVisited) + { + const bool alreadyActive = activePoints.CompareAndSwapBitAtomic(pointId, true, false); + if (!alreadyActive) + { // If we're the first thread to mark point active, set ref point: + refPoints.Set(pointId, refPtId); + } + } + } + + // Mark current cell as inactive: + return false; + } + }; + + // For each point with a refPtId set, ensure that the associated normal is + // in the same hemisphere as the reference normal. + // This must be done in a separate step from MarkActivePoints since modifying + // visitedPoints in that worklet would create race conditions. + class WorkletProcessNormals : public vtkm::worklet::WorkletMapField + { + public: + using ControlSignature = void(FieldIn refIds, + WholeArrayInOut normals, + // InOut to preserve data on masked indices + BitFieldInOut visitedPoints); + using ExecutionSignature = void(InputIndex ptId, _1 refPtId, _2 normals, _3 visitedPoints); + + using MaskType = vtkm::worklet::MaskIndices; + + template + VTKM_EXEC void operator()(const vtkm::Id ptId, + const vtkm::Id refPtId, + NormalsPortal& normals, + VisitedPointsT& visitedPoints) const + { + visitedPoints.SetBitAtomic(ptId, true); + + using Normal = typename NormalsPortal::ValueType; + Normal normal = normals.Get(ptId); + const Normal ref = normals.Get(refPtId); + if (Align(normal, ref)) + { + normals.Set(ptId, normal); + } + } + }; + + template + VTKM_CONT static void Run( + const CellSetType& cells, + const vtkm::cont::ArrayHandle, CoordsStorageType>& coords, + vtkm::cont::ArrayHandle, PointNormalStorageType>& + pointNormals) + { + using RangeType = vtkm::cont::ArrayHandle; + + using MarkSourcePoints = vtkm::worklet::DispatcherMapField; + using MarkActiveCells = vtkm::worklet::DispatcherMapTopology; + using MarkActivePoints = vtkm::worklet::DispatcherMapTopology; + using ProcessNormals = vtkm::worklet::DispatcherMapField; + + const vtkm::Id numCells = cells.GetNumberOfCells(); + + VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf, + "OrientPointNormals worklet (%lld points, %lld cells)", + static_cast(coords.GetNumberOfValues()), + static_cast(numCells)); + + // active = cells / point to be used in the next worklet invocation mask. + vtkm::cont::BitField activePointBits; // Initialized by MarkSourcePoints + auto activePoints = vtkm::cont::make_ArrayHandleBitField(activePointBits); + + vtkm::cont::BitField activeCellBits; + vtkm::cont::Algorithm::Fill(activeCellBits, false, numCells); + auto activeCells = vtkm::cont::make_ArrayHandleBitField(activeCellBits); + + // visited = cells / points that have been corrected. + vtkm::cont::BitField visitedPointBits; // Initialized by MarkSourcePoints + auto visitedPoints = vtkm::cont::make_ArrayHandleBitField(visitedPointBits); + + vtkm::cont::BitField visitedCellBits; + vtkm::cont::Algorithm::Fill(visitedCellBits, false, numCells); + auto visitedCells = vtkm::cont::make_ArrayHandleBitField(visitedCellBits); + + vtkm::cont::ArrayHandle mask; // Allocated as needed + + // For each point, store a reference alignment point. Allocated by + // MarkSourcePoints. + vtkm::cont::ArrayHandle refPoints; + + // 1) Compute range of coords. + const RangeType ranges = vtkm::cont::ArrayRangeCompute(coords); + + // 2) Label source points for traversal (use those on a boundary). + // Correct the normals for these points by making them point towards the + // boundary. + { + MarkSourcePoints dispatcher; + dispatcher.Invoke(coords, pointNormals, ranges, activePoints, visitedPoints, refPoints); + } + + for (size_t iter = 1;; ++iter) + { + // 3) Mark unvisited cells adjacent to active points + { + vtkm::Id numActive = vtkm::cont::Algorithm::BitFieldToUnorderedSet(activePointBits, mask); + (void)numActive; + VTKM_LOG_S(vtkm::cont::LogLevel::Perf, + "MarkActiveCells from " << numActive << " active points."); + MarkActiveCells dispatcher{ vtkm::worklet::MaskIndices{ mask } }; + dispatcher.Invoke(cells, activeCellBits, visitedCellBits, activePoints); + } + + // 4) Mark unvisited points in active cells, using ref point from cell. + { + vtkm::Id numActive = vtkm::cont::Algorithm::BitFieldToUnorderedSet(activeCellBits, mask); + (void)numActive; + VTKM_LOG_S(vtkm::cont::LogLevel::Perf, + "MarkActivePoints from " << numActive << " active cells."); + MarkActivePoints dispatcher{ vtkm::worklet::MaskIndices{ mask } }; + dispatcher.Invoke(cells, activePointBits, visitedPointBits, refPoints, activeCells); + } + + vtkm::Id numActivePoints = + vtkm::cont::Algorithm::BitFieldToUnorderedSet(activePointBits, mask); + + if (numActivePoints == 0) + { // Done! + VTKM_LOG_S(vtkm::cont::LogLevel::Perf, "Iteration " << iter << ": Traversal complete."); + break; + } + + VTKM_LOG_S(vtkm::cont::LogLevel::Perf, + "Iteration " << iter << ": Processing " << numActivePoints << " normals."); + + // 5) Correct normals for active points. + { + ProcessNormals dispatcher{ vtkm::worklet::MaskIndices{ mask } }; + dispatcher.Invoke(refPoints, pointNormals, visitedPointBits); + } + } + } +}; +} +} // end namespace vtkm::worklet + + +#endif // vtkm_m_worklet_OrientPointNormals_h diff --git a/vtkm/worklet/testing/CMakeLists.txt b/vtkm/worklet/testing/CMakeLists.txt index fd3db8708..9d81bdf5d 100644 --- a/vtkm/worklet/testing/CMakeLists.txt +++ b/vtkm/worklet/testing/CMakeLists.txt @@ -45,6 +45,7 @@ set(unit_tests UnitTestNDimsEntropy.cxx UnitTestNDimsHistogram.cxx UnitTestNDimsHistMarginalization.cxx + UnitTestOrientNormals.cxx UnitTestParticleAdvection.cxx UnitTestPointElevation.cxx UnitTestPointGradient.cxx diff --git a/vtkm/worklet/testing/UnitTestOrientNormals.cxx b/vtkm/worklet/testing/UnitTestOrientNormals.cxx new file mode 100644 index 000000000..b77ebbaa7 --- /dev/null +++ b/vtkm/worklet/testing/UnitTestOrientNormals.cxx @@ -0,0 +1,408 @@ +//============================================================================= +// +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +// +// Copyright 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2019 UT-Battelle, LLC. +// Copyright 2019 Los Alamos National Security. +// +// Under the terms of Contract DE-NA0003525 with NTESS, +// the U.S. Government retains certain rights in this software. +// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National +// Laboratory (LANL), the U.S. Government retains certain rights in +// this software. +// +//============================================================================= + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include + +#include + + +namespace +{ + +struct TestPolicy : public vtkm::filter::PolicyBase +{ + using StructuredCellSetList = vtkm::ListTagBase>; + using UnstructuredCellSetList = vtkm::ListTagBase>; + using AllCellSetList = vtkm::ListTagJoin; + using FieldTypeList = vtkm::ListTagBase>; +}; + +VTKM_CONT +vtkm::cont::DataSet CreateDataSet(bool pointNormals, bool cellNormals) +{ + vtkm::worklet::WaveletGenerator wavelet; + wavelet.SetMinimumExtent({ -25 }); + wavelet.SetMaximumExtent({ 25 }); + wavelet.SetFrequency({ 20, 15, 25 }); + wavelet.SetMagnitude({ 5 }); + auto dataSet = wavelet.GenerateDataSet(); + + // Cut a contour + vtkm::filter::MarchingCubes contour; + contour.SetActiveField("scalars", vtkm::cont::Field::Association::POINTS); + contour.SetNumberOfIsoValues(1); + contour.SetIsoValue(192); + contour.SetMergeDuplicatePoints(true); + contour.SetGenerateNormals(false); + dataSet = contour.Execute(dataSet, TestPolicy{}); + + vtkm::filter::SurfaceNormals normals; + normals.SetGeneratePointNormals(pointNormals); + normals.SetGenerateCellNormals(cellNormals); + normals.SetPointNormalsName("normals"); + normals.SetCellNormalsName("normals"); + normals.SetAutoOrientNormals(false); + dataSet = normals.Execute(dataSet, TestPolicy{}); + + return dataSet; +} + +struct ValidateNormals +{ + using CellSetType = vtkm::cont::CellSetSingleType<>; + using NormalType = vtkm::Vec; + using NormalsArrayType = vtkm::cont::ArrayHandleVirtual; + using NormalsPortalType = decltype(std::declval().GetPortalConstControl()); + using ConnType = + decltype(std::declval().PrepareForInput(vtkm::cont::DeviceAdapterTagSerial{}, + vtkm::TopologyElementTagPoint{}, + vtkm::TopologyElementTagCell{})); + using RConnType = + decltype(std::declval().PrepareForInput(vtkm::cont::DeviceAdapterTagSerial{}, + vtkm::TopologyElementTagCell{}, + vtkm::TopologyElementTagPoint{})); + using PointsType = + decltype(std::declval().GetData().GetPortalConstControl()); + + vtkm::cont::CoordinateSystem Coords; + CellSetType Cells; + + PointsType Points; + ConnType Conn; + RConnType RConn; + + NormalsArrayType PointNormalsArray; + NormalsPortalType PointNormals; + NormalsArrayType CellNormalsArray; + NormalsPortalType CellNormals; + + vtkm::cont::BitField VisitedCellsField; + vtkm::cont::BitField VisitedPointsField; + vtkm::cont::BitField::PortalControl VisitedCells; + vtkm::cont::BitField::PortalControl VisitedPoints; + + bool CheckPoints; + bool CheckCells; + + VTKM_CONT + static void Run(vtkm::cont::DataSet& dataset, + bool checkPoints, + bool checkCells, + const std::string& normalsName = "normals") + { + vtkm::cont::Field pointNormals; + vtkm::cont::Field cellNormals; + + if (checkPoints) + { + pointNormals = dataset.GetField(normalsName, vtkm::cont::Field::Association::POINTS); + } + if (checkCells) + { + cellNormals = dataset.GetField(normalsName, vtkm::cont::Field::Association::CELL_SET); + } + + ValidateNormals obj{ dataset, checkPoints, checkCells, pointNormals, cellNormals }; + obj.Validate(); + } + + VTKM_CONT + ValidateNormals(const vtkm::cont::DataSet& dataset, + bool checkPoints, + bool checkCells, + const vtkm::cont::Field& pointNormalsField, + const vtkm::cont::Field& cellNormalsField) + : Coords{ dataset.GetCoordinateSystem() } + , Cells{ dataset.GetCellSet().Cast() } + , Points{ this->Coords.GetData().GetPortalConstControl() } + , CheckPoints(checkPoints) + , CheckCells(checkCells) + { + // FIXME This would be much simplier if we had a GetPointCells() method on + // cell sets.... + // Build the connectivity table on any device, then get a portal for serial + // so we can do lookups on the CPU. + this->Cells.GetConnectivityArray(vtkm::TopologyElementTagPoint{}, + vtkm::TopologyElementTagCell{}); + this->Cells.GetConnectivityArray(vtkm::TopologyElementTagCell{}, + vtkm::TopologyElementTagPoint{}); + this->Conn = this->Cells.PrepareForInput(vtkm::cont::DeviceAdapterTagSerial{}, + vtkm::TopologyElementTagPoint{}, + vtkm::TopologyElementTagCell{}); + this->RConn = this->Cells.PrepareForInput(vtkm::cont::DeviceAdapterTagSerial{}, + vtkm::TopologyElementTagCell{}, + vtkm::TopologyElementTagPoint{}); + + if (this->CheckPoints) + { + this->PointNormalsArray = pointNormalsField.GetData().AsVirtual(); + this->PointNormals = this->PointNormalsArray.GetPortalConstControl(); + } + if (this->CheckCells) + { + this->CellNormalsArray = cellNormalsField.GetData().AsVirtual(); + this->CellNormals = this->CellNormalsArray.GetPortalConstControl(); + } + } + + VTKM_CONT + void Validate() + { + // Locate a point with the minimum x coordinate: + const vtkm::Id startPoint = [&]() -> vtkm::Id { + const vtkm::Float64 xMin = this->Coords.GetBounds().X.Min; + const auto points = this->Coords.GetData().GetPortalConstControl(); + const vtkm::Id numPoints = points.GetNumberOfValues(); + vtkm::Id resultIdx = -1; + for (vtkm::Id pointIdx = 0; pointIdx < numPoints; ++pointIdx) + { + const auto point = points.Get(pointIdx); + if (static_cast(point[0]) <= xMin) + { + resultIdx = pointIdx; + break; + } + } + if (resultIdx < 0) + { + throw vtkm::cont::ErrorBadValue("Minimum point not found!"); + } + + return resultIdx; + }(); + + // Start recursive validation. + this->Prepare(); + this->ValidateImpl(startPoint, NormalType{ -1, 0, 0 }); + + vtkm::Id numPoints = this->Points.GetNumberOfValues(); + vtkm::Id numCells = this->Cells.GetNumberOfCells(); + vtkm::Id numVisitedPoints = vtkm::cont::Algorithm::CountSetBits(this->VisitedPointsField); + vtkm::Id numVisitedCells = vtkm::cont::Algorithm::CountSetBits(this->VisitedCellsField); + if (numPoints != numVisitedPoints) + { + throw vtkm::cont::ErrorBadValue("Unvisited point!"); + } + if (numCells != numVisitedCells) + { + throw vtkm::cont::ErrorBadValue("Unvisited cell!"); + } + } + +private: + static bool SameHemisphere(const NormalType& a, const NormalType& b) + { + return vtkm::Dot(a, b) >= 0; + } + + void Prepare() + { + vtkm::cont::Algorithm::Fill(this->VisitedPointsField, false, this->Coords.GetNumberOfPoints()); + vtkm::cont::Algorithm::Fill(this->VisitedCellsField, false, this->Cells.GetNumberOfCells()); + + this->VisitedPoints = this->VisitedPointsField.GetPortalControl(); + this->VisitedCells = this->VisitedCellsField.GetPortalControl(); + } + + void ValidateImpl(vtkm::Id startPtIdx, const NormalType& startRefNormal) + { + using Entry = vtkm::Pair; + std::vector queue; + queue.emplace_back(startPtIdx, startRefNormal); + this->VisitedPoints.SetBit(startPtIdx, true); + + while (!queue.empty()) + { + const vtkm::Id curPtIdx = queue.back().first; + NormalType refNormal = queue.back().second; + queue.pop_back(); + + if (this->CheckPoints) + { + const NormalType curNormal = this->PointNormals.Get(curPtIdx); + if (!this->SameHemisphere(curNormal, refNormal)) + { + const auto coord = this->Points.Get(curPtIdx); + std::cerr << "BadPointNormal PtId: " << curPtIdx << "\n\t" + << "- Normal: {" << curNormal[0] << ", " << curNormal[1] << ", " << curNormal[2] + << "}\n\t" + << "- Reference: {" << refNormal[0] << ", " << refNormal[1] << ", " + << refNormal[2] << "}\n\t" + << "- Coord: {" << coord[0] << ", " << coord[1] << ", " << coord[2] << "}\n"; + throw vtkm::cont::ErrorBadValue("Bad point normal found!"); + } + refNormal = curNormal; + } + + // Lookup and visit neighbor cells: + const auto neighborCells = this->RConn.GetIndices(curPtIdx); + const auto numNeighborCells = neighborCells.GetNumberOfComponents(); + for (vtkm::IdComponent nCellIdx = 0; nCellIdx < numNeighborCells; ++nCellIdx) + { + const vtkm::Id curCellIdx = neighborCells[nCellIdx]; + + // Skip this cell if already visited: + if (this->VisitedCells.GetBit(curCellIdx)) + { + continue; + } + this->VisitedCells.SetBit(curCellIdx, true); + + if (this->CheckCells) + { + const NormalType curNormal = this->CellNormals.Get(curCellIdx); + if (!this->SameHemisphere(curNormal, refNormal)) + { + std::cerr << "BadCellNormal: CellId: " << curCellIdx << "\n\t" + << "- Normal: {" << curNormal[0] << ", " << curNormal[1] << ", " + << curNormal[2] << "}\n\t" + << "- Reference: {" << refNormal[0] << ", " << refNormal[1] << ", " + << refNormal[2] << "}\n"; + throw vtkm::cont::ErrorBadValue("Bad cell normal found!"); + } + refNormal = curNormal; + } + + // Lookup and visit points in this cell: + const auto neighborPoints = this->Conn.GetIndices(curCellIdx); + const auto numNeighborPoints = neighborPoints.GetNumberOfComponents(); + for (vtkm::IdComponent nPtIdx = 0; nPtIdx < numNeighborPoints; ++nPtIdx) + { + const vtkm::Id nextPtIdx = neighborPoints[nPtIdx]; + + // Skip if already visited: + if (this->VisitedPoints.GetBit(nextPtIdx)) + { + continue; + } + + // Otherwise, queue next point using current normal as reference: + queue.emplace_back(nextPtIdx, refNormal); + this->VisitedPoints.SetBit(nextPtIdx, true); + } + } + } + } +}; + +VTKM_CONT +void TestOrientNormals(bool testPoints, bool testCells) +{ + using NormalArrayT = vtkm::cont::ArrayHandle>; + + auto dataset = CreateDataSet(testPoints, testCells); + + // Check that the input actually has bad normals: + const bool inputValid = [&]() -> bool { + try + { + std::cerr << "Expecting to throw a BadNormal exception...\n"; + ValidateNormals::Run(dataset, testPoints, testCells); + return true; // Dataset is already oriented + } + catch (vtkm::cont::ErrorBadValue&) + { + std::cerr << "Expected exception caught! All is well. " + "Future exceptions are errors.\n"; + return false; // Dataset is unoriented + } + }(); + + if (inputValid) + { + throw vtkm::cont::ErrorBadValue("Error: Input doesn't have bad normals."); + } + + // modify normals in place + const auto coords = dataset.GetCoordinateSystem().GetData(); + const auto cells = dataset.GetCellSet(); + if (testPoints && testCells) + { + std::cerr << "Testing point and cell normals...\n"; + + const auto pointNormalField = + dataset.GetField("normals", vtkm::cont::Field::Association::POINTS); + const auto cellNormalField = + dataset.GetField("normals", vtkm::cont::Field::Association::CELL_SET); + auto pointNormals = pointNormalField.GetData().Cast(); + auto cellNormals = cellNormalField.GetData().Cast(); + + vtkm::worklet::OrientNormals::RunPointAndCellNormals(cells, coords, pointNormals, cellNormals); + } + else if (testPoints) + { + std::cerr << "Testing point normals...\n"; + const auto pointNormalField = + dataset.GetField("normals", vtkm::cont::Field::Association::POINTS); + auto pointNormals = pointNormalField.GetData().Cast(); + + vtkm::worklet::OrientNormals::RunPointNormals(cells, coords, pointNormals); + } + else if (testCells) + { + std::cerr << "Testing cell normals...\n"; + const auto cellNormalField = + dataset.GetField("normals", vtkm::cont::Field::Association::CELL_SET); + auto cellNormals = cellNormalField.GetData().Cast(); + + vtkm::worklet::OrientNormals::RunCellNormals(cells, coords, cellNormals); + } + else + { + throw "Nothing tested..."; + } + + ValidateNormals::Run(dataset, testPoints, testCells); +} + +void DoTest() +{ + TestOrientNormals(true, false); + TestOrientNormals(false, true); + TestOrientNormals(true, true); +} + +} // end anon namespace + +int UnitTestOrientNormals(int argc, char* argv[]) +{ + return vtkm::cont::testing::Testing::Run(DoTest, argc, argv); +} From 320a5257ad56b051b26fcd87bca8bf2088c9db9f Mon Sep 17 00:00:00 2001 From: Allison Vacanti Date: Wed, 26 Jun 2019 12:53:15 -0400 Subject: [PATCH 3/4] Add TriangleWinding worklet. --- vtkm/worklet/CMakeLists.txt | 1 + vtkm/worklet/TriangleWinding.h | 185 ++++++++++++++++++ vtkm/worklet/testing/CMakeLists.txt | 1 + .../testing/UnitTestTriangleWinding.cxx | 139 +++++++++++++ 4 files changed, 326 insertions(+) create mode 100644 vtkm/worklet/TriangleWinding.h create mode 100644 vtkm/worklet/testing/UnitTestTriangleWinding.cxx diff --git a/vtkm/worklet/CMakeLists.txt b/vtkm/worklet/CMakeLists.txt index a23e9ea76..4c05d377f 100644 --- a/vtkm/worklet/CMakeLists.txt +++ b/vtkm/worklet/CMakeLists.txt @@ -75,6 +75,7 @@ set(headers Tetrahedralize.h Threshold.h ThresholdPoints.h + TriangleWinding.h Triangulate.h Tube.h VertexClustering.h diff --git a/vtkm/worklet/TriangleWinding.h b/vtkm/worklet/TriangleWinding.h new file mode 100644 index 000000000..28e444c02 --- /dev/null +++ b/vtkm/worklet/TriangleWinding.h @@ -0,0 +1,185 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +// +// Copyright 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2019 UT-Battelle, LLC. +// Copyright 2019 Los Alamos National Security. +// +// Under the terms of Contract DE-NA0003525 with NTESS, +// the U.S. Government retains certain rights in this software. +// +// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National +// Laboratory (LANL), the U.S. Government retains certain rights in +// this software. +//============================================================================ + +#ifndef vtkm_m_worklet_TriangleWinding_h +#define vtkm_m_worklet_TriangleWinding_h + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +namespace vtkm +{ +namespace worklet +{ + +/** + * This worklet ensures that triangle windings are consistent with provided + * cell normals. The triangles are wound CCW around the cell normals, and + * all other cells are ignored. + * + * The input cellset must be unstructured. + */ +class TriangleWinding +{ +public: + struct WorkletWindToCellNormals : public WorkletMapField + { + using ControlSignature = void(FieldIn cellNormals, FieldInOut cellPoints, WholeArrayIn coords); + using ExecutionSignature = void(_1 cellNormal, _2 cellPoints, _3 coords); + + template + VTKM_EXEC void operator()(const vtkm::Vec& cellNormal, + CellPointsType& cellPoints, + const CoordsPortal& coords) const + { + // We only care about triangles: + if (cellPoints.GetNumberOfComponents() != 3) + { + return; + } + + using NormalType = vtkm::Vec; + + const NormalType p0 = coords.Get(cellPoints[0]); + const NormalType p1 = coords.Get(cellPoints[1]); + const NormalType p2 = coords.Get(cellPoints[2]); + const NormalType v01 = p1 - p0; + const NormalType v02 = p2 - p0; + const NormalType triangleNormal = vtkm::Cross(v01, v02); + if (vtkm::Dot(cellNormal, triangleNormal) < 0) + { + // Can't just use std::swap from exec function: + const vtkm::Id tmp = cellPoints[1]; + cellPoints[1] = cellPoints[2]; + cellPoints[2] = tmp; + } + } + }; + + struct Launcher + { + vtkm::cont::DynamicCellSet Result; + + template + VTKM_CONT void operator()( + const vtkm::cont::CellSetExplicit& cellSet, + const vtkm::cont::ArrayHandle, PointStorageType>& coords, + const vtkm::cont::ArrayHandle, CellNormalStorageType>& + cellNormals) + { + using WindToCellNormals = vtkm::worklet::DispatcherMapField; + + vtkm::cont::ArrayHandle conn; + { + const auto& connIn = cellSet.GetConnectivityArray(vtkm::TopologyElementTagPoint{}, + vtkm::TopologyElementTagCell{}); + vtkm::cont::Algorithm::Copy(connIn, conn); + } + + const auto& offsets = cellSet.GetIndexOffsetArray(vtkm::TopologyElementTagPoint{}, + vtkm::TopologyElementTagCell{}); + auto cells = vtkm::cont::make_ArrayHandleGroupVecVariable(conn, offsets); + + WindToCellNormals dispatcher; + dispatcher.Invoke(cellNormals, cells, coords); + + const auto& shapes = + cellSet.GetShapesArray(vtkm::TopologyElementTagPoint{}, vtkm::TopologyElementTagCell{}); + const auto& numIndices = + cellSet.GetNumIndicesArray(vtkm::TopologyElementTagPoint{}, vtkm::TopologyElementTagCell{}); + vtkm::cont::CellSetExplicit newCells; + newCells.Fill(cellSet.GetNumberOfPoints(), shapes, numIndices, conn, offsets); + + this->Result = newCells; + } + + template + VTKM_CONT void operator()( + const vtkm::cont::CellSetSingleType& cellSet, + const vtkm::cont::ArrayHandle, PointStorageType>& coords, + const vtkm::cont::ArrayHandle, CellNormalStorageType>& + cellNormals) + { + using WindToCellNormals = vtkm::worklet::DispatcherMapField; + + vtkm::cont::ArrayHandle conn; + { + const auto& connIn = cellSet.GetConnectivityArray(vtkm::TopologyElementTagPoint{}, + vtkm::TopologyElementTagCell{}); + vtkm::cont::Algorithm::Copy(connIn, conn); + } + + const auto& offsets = cellSet.GetIndexOffsetArray(vtkm::TopologyElementTagPoint{}, + vtkm::TopologyElementTagCell{}); + auto cells = vtkm::cont::make_ArrayHandleGroupVecVariable(conn, offsets); + + WindToCellNormals dispatcher; + dispatcher.Invoke(cellNormals, cells, coords); + + vtkm::cont::CellSetSingleType newCells; + newCells.Fill(cellSet.GetNumberOfPoints(), + cellSet.GetCellShape(0), + cellSet.GetNumberOfPointsInCell(0), + conn); + + this->Result = newCells; + } + }; + + template + VTKM_CONT static vtkm::cont::DynamicCellSet Run( + const CellSetType& cellSet, + const vtkm::cont::ArrayHandle, PointStorageType>& coords, + const vtkm::cont::ArrayHandle, CellNormalStorageType>& + cellNormals) + { + Launcher launcher; + vtkm::cont::CastAndCall(cellSet, launcher, coords, cellNormals); + return launcher.Result; + } +}; +} +} // end namespace vtkm::worklet + +#endif // vtkm_m_worklet_TriangleWinding_h diff --git a/vtkm/worklet/testing/CMakeLists.txt b/vtkm/worklet/testing/CMakeLists.txt index 9d81bdf5d..6514b2e67 100644 --- a/vtkm/worklet/testing/CMakeLists.txt +++ b/vtkm/worklet/testing/CMakeLists.txt @@ -66,6 +66,7 @@ set(unit_tests UnitTestTetrahedralize.cxx UnitTestThreshold.cxx UnitTestThresholdPoints.cxx + UnitTestTriangleWinding.cxx UnitTestTriangulate.cxx UnitTestTube.cxx UnitTestWholeCellSetIn.cxx diff --git a/vtkm/worklet/testing/UnitTestTriangleWinding.cxx b/vtkm/worklet/testing/UnitTestTriangleWinding.cxx new file mode 100644 index 000000000..86222a2cb --- /dev/null +++ b/vtkm/worklet/testing/UnitTestTriangleWinding.cxx @@ -0,0 +1,139 @@ +//============================================================================= +// +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +// +// Copyright 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2019 UT-Battelle, LLC. +// Copyright 2019 Los Alamos National Security. +// +// Under the terms of Contract DE-NA0003525 with NTESS, +// the U.S. Government retains certain rights in this software. +// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National +// Laboratory (LANL), the U.S. Government retains certain rights in +// this software. +// +//============================================================================= + +#include + +#include + +#include +#include +#include +#include +#include + +#include +#include + +using NormalType = vtkm::Vec; + +namespace +{ + +vtkm::cont::DataSet GenerateDataSet() +{ + auto ds = vtkm::cont::testing::MakeTestDataSet{}.Make3DExplicitDataSetPolygonal(); + const auto numCells = ds.GetCellSet().GetNumberOfCells(); + + vtkm::cont::ArrayHandle cellNormals; + vtkm::cont::Algorithm::Fill(cellNormals, NormalType{ 1., 0., 0. }, numCells); + + ds.AddField(vtkm::cont::Field{ + "normals", vtkm::cont::Field::Association::CELL_SET, ds.GetCellSet().GetName(), cellNormals }); + return ds; +} + +void Validate(vtkm::cont::DataSet dataSet) +{ + const auto cellSet = dataSet.GetCellSet().Cast>(); + const auto coordsArray = dataSet.GetCoordinateSystem().GetData(); + const auto conn = + cellSet.GetConnectivityArray(vtkm::TopologyElementTagPoint{}, vtkm::TopologyElementTagCell{}); + const auto offsets = + cellSet.GetIndexOffsetArray(vtkm::TopologyElementTagPoint{}, vtkm::TopologyElementTagCell{}); + const auto cellArray = vtkm::cont::make_ArrayHandleGroupVecVariable(conn, offsets); + const auto cellNormalsVar = + dataSet.GetField("normals", vtkm::cont::Field::Association::CELL_SET).GetData(); + const auto cellNormalsArray = cellNormalsVar.Cast>(); + + const auto cellPortal = cellArray.GetPortalConstControl(); + const auto cellNormals = cellNormalsArray.GetPortalConstControl(); + const auto coords = coordsArray.GetPortalConstControl(); + + const auto numCells = cellPortal.GetNumberOfValues(); + VTKM_TEST_ASSERT(numCells == cellNormals.GetNumberOfValues()); + + for (vtkm::Id cellId = 0; cellId < numCells; ++cellId) + { + const auto cell = cellPortal.Get(cellId); + if (cell.GetNumberOfComponents() != 3) + { // Triangles only! + continue; + } + + const NormalType cellNormal = cellNormals.Get(cellId); + const NormalType p0 = coords.Get(cell[0]); + const NormalType p1 = coords.Get(cell[1]); + const NormalType p2 = coords.Get(cell[2]); + const NormalType v01 = p1 - p0; + const NormalType v02 = p2 - p0; + const NormalType triangleNormal = vtkm::Cross(v01, v02); + VTKM_TEST_ASSERT(vtkm::Dot(triangleNormal, cellNormal) > 0, + "Triangle at index ", + cellId, + " incorrectly wound."); + } +} + +void DoTest() +{ + auto ds = GenerateDataSet(); + + // Ensure that the test dataset needs to be rewound: + bool threw = false; + try + { + std::cerr << "Expecting an exception...\n"; + Validate(ds); + } + catch (...) + { + threw = true; + } + + VTKM_TEST_ASSERT(threw, "Test dataset is already wound consistently wrt normals."); + + auto cellSet = ds.GetCellSet().Cast>(); + const auto coords = ds.GetCoordinateSystem().GetData(); + const auto cellNormalsVar = + ds.GetField("normals", vtkm::cont::Field::Association::CELL_SET).GetData(); + const auto cellNormals = cellNormalsVar.Cast>(); + + + auto newCells = vtkm::worklet::TriangleWinding::Run(cellSet, coords, cellNormals); + + vtkm::cont::DataSet result; + result.AddCoordinateSystem(ds.GetCoordinateSystem()); + result.AddCellSet(newCells); + for (vtkm::Id i = 0; i < ds.GetNumberOfFields(); ++i) + { + result.AddField(ds.GetField(i)); + } + + Validate(result); +} + +} // end anon namespace + +int UnitTestTriangleWinding(int argc, char* argv[]) +{ + return vtkm::cont::testing::Testing::Run(DoTest, argc, argv); +} From 91e7c4d0b9b43120f2f6dc12c01c9fd5c1c9baa5 Mon Sep 17 00:00:00 2001 From: Allison Vacanti Date: Wed, 26 Jun 2019 12:53:46 -0400 Subject: [PATCH 4/4] Add orientation and winding options to SurfaceNormals. --- docs/changelog/orient-normals.md | 18 +++++ vtkm/filter/SurfaceNormals.h | 39 ++++++++++ vtkm/filter/SurfaceNormals.hxx | 78 ++++++++++++++++++- .../testing/UnitTestSurfaceNormalsFilter.cxx | 1 + 4 files changed, 132 insertions(+), 4 deletions(-) create mode 100644 docs/changelog/orient-normals.md diff --git a/docs/changelog/orient-normals.md b/docs/changelog/orient-normals.md new file mode 100644 index 000000000..c84a3daf6 --- /dev/null +++ b/docs/changelog/orient-normals.md @@ -0,0 +1,18 @@ +# `SurfaceNormals` filter can now orient normals. + +The `OrientNormals` worklet has been added to the `SurfaceNormals` filter, and +is enabled by turning on the `AutoOrientNormals` option. This feature ensures +that all normals generated by the filter will point out of the dataset (or +inward if the `FlipNormals` option is true). In addition, +`SurfaceNormals` now has a `Consistency` option that forces all triangle +windings to be consistent with the cell normal direction (the cell points are +specified in counter-clockwise order around the normal). + +This functionality is provided by the following new worklets: + +* OrientNormals + * RunOrientCellNormals + * RunOrientPointNormals + * RunOrientPointAndCellNormals + * RunFlipNormals +* TriangleWinding diff --git a/vtkm/filter/SurfaceNormals.h b/vtkm/filter/SurfaceNormals.h index 82e62aee8..cc3447a6c 100644 --- a/vtkm/filter/SurfaceNormals.h +++ b/vtkm/filter/SurfaceNormals.h @@ -32,28 +32,64 @@ public: SurfaceNormals(); /// Set/Get if cell normals should be generated. Default is off. + /// @{ void SetGenerateCellNormals(bool value) { this->GenerateCellNormals = value; } bool GetGenerateCellNormals() const { return this->GenerateCellNormals; } + /// @} /// Set/Get if the cell normals should be normalized. Default value is true. /// The intended use case of this flag is for faster, approximate point /// normals generation by skipping the normalization of the face normals. /// Note that when set to false, the result cell normals will not be unit /// length normals and the point normals will be different. + /// @{ void SetNormalizeCellNormals(bool value) { this->NormalizeCellNormals = value; } bool GetNormalizeCellNormals() const { return this->NormalizeCellNormals; } + /// @} /// Set/Get if the point normals should be generated. Default is on. + /// @{ void SetGeneratePointNormals(bool value) { this->GeneratePointNormals = value; } bool GetGeneratePointNormals() const { return this->GeneratePointNormals; } + /// @} /// Set/Get the name of the cell normals field. Default is "Normals". + /// @{ void SetCellNormalsName(const std::string& name) { this->CellNormalsName = name; } const std::string& GetCellNormalsName() const { return this->CellNormalsName; } + /// @} /// Set/Get the name of the point normals field. Default is "Normals". + /// @{ void SetPointNormalsName(const std::string& name) { this->PointNormalsName = name; } const std::string& GetPointNormalsName() const { return this->PointNormalsName; } + /// @} + + /// If true, the normals will be oriented to face outwards from the surface. + /// This requires a closed manifold surface or the behavior is undefined. + /// This option is expensive but necessary for rendering. + /// To make the normals point inward, set FlipNormals to true. + /// @{ + void SetAutoOrientNormals(bool v) { this->AutoOrientNormals = v; } + bool GetAutoOrientNormals() const { return this->AutoOrientNormals; } + /// @} + + /// Reverse the normals to point inward when AutoOrientNormals is true. + /// Default is false. + /// @{ + void SetFlipNormals(bool v) { this->FlipNormals = v; } + bool GetFlipNormals() const { return this->FlipNormals; } + /// @} + + /// Ensure that polygon winding is consistent with normal orientation. + /// Triangles are wound such that their points are counter-clockwise around + /// the generated cell normal. Default is true. + /// @note This currently only affects triangles. + /// @note This is only applied when cell normals are generated. + /// @{ + void SetConsistency(bool v) { this->Consistency = v; } + bool GetConsistency() const { return this->Consistency; } + /// @} template vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input, @@ -65,6 +101,9 @@ private: bool GenerateCellNormals; bool NormalizeCellNormals; bool GeneratePointNormals; + bool AutoOrientNormals; + bool FlipNormals; + bool Consistency; std::string CellNormalsName; std::string PointNormalsName; diff --git a/vtkm/filter/SurfaceNormals.hxx b/vtkm/filter/SurfaceNormals.hxx index 4d703247d..27fd42f4f 100644 --- a/vtkm/filter/SurfaceNormals.hxx +++ b/vtkm/filter/SurfaceNormals.hxx @@ -8,9 +8,13 @@ // PURPOSE. See the above copyright notice for more information. //============================================================================ +#include + #include #include +#include #include +#include namespace vtkm { @@ -58,6 +62,9 @@ inline SurfaceNormals::SurfaceNormals() : GenerateCellNormals(false) , NormalizeCellNormals(true) , GeneratePointNormals(true) + , AutoOrientNormals(false) + , FlipNormals(false) + , Consistency(true) { this->SetUseCoordinateSystemAsField(true); } @@ -76,19 +83,22 @@ inline vtkm::cont::DataSet SurfaceNormals::DoExecute( throw vtkm::cont::ErrorFilterExecution("No normals selected."); } - const auto& cellset = input.GetCellSet(this->GetActiveCellSetIndex()); + const auto cellset = + vtkm::filter::ApplyPolicyUnstructured(input.GetCellSet(this->GetActiveCellSetIndex()), policy); + const auto& coords = input.GetCoordinateSystem(this->GetActiveCoordinateSystemIndex()).GetData(); vtkm::cont::ArrayHandle faceNormals; vtkm::worklet::FacetedSurfaceNormals faceted; faceted.SetNormalize(this->NormalizeCellNormals); - faceted.Run(vtkm::filter::ApplyPolicy(cellset, policy), points, faceNormals); + faceted.Run(cellset, points, faceNormals); vtkm::cont::DataSet result; + vtkm::cont::ArrayHandle pointNormals; if (this->GeneratePointNormals) { - vtkm::cont::ArrayHandle pointNormals; vtkm::worklet::SmoothSurfaceNormals smooth; - smooth.Run(vtkm::filter::ApplyPolicy(cellset, policy), faceNormals, pointNormals); + smooth.Run(cellset, faceNormals, pointNormals); + result = internal::CreateResult(input, pointNormals, @@ -111,6 +121,66 @@ inline vtkm::cont::DataSet SurfaceNormals::DoExecute( cellset.GetName()); } + if (this->AutoOrientNormals) + { + using Orient = vtkm::worklet::OrientNormals; + + if (this->GenerateCellNormals && this->GeneratePointNormals) + { + Orient::RunPointAndCellNormals(cellset, coords, pointNormals, faceNormals); + } + else if (this->GenerateCellNormals) + { + Orient::RunCellNormals(cellset, coords, faceNormals); + } + else if (this->GeneratePointNormals) + { + Orient::RunPointNormals(cellset, coords, pointNormals); + } + + if (this->FlipNormals) + { + if (this->GenerateCellNormals) + { + Orient::RunFlipNormals(faceNormals); + } + if (this->GeneratePointNormals) + { + Orient::RunFlipNormals(pointNormals); + } + } + } + + if (this->Consistency && this->GenerateCellNormals) + { + auto newCells = vtkm::worklet::TriangleWinding::Run(cellset, coords, faceNormals); + + // Add the new cells into the result: + vtkm::cont::DataSet newResult; + for (vtkm::Id i = 0; i < result.GetNumberOfCoordinateSystems(); ++i) + { + newResult.AddCoordinateSystem(result.GetCoordinateSystem(i)); + } + const vtkm::Id activeCells = this->GetActiveCellSetIndex(); + for (vtkm::Id i = 0; i < result.GetNumberOfCellSets(); ++i) + { + if (i != activeCells) + { + newResult.AddCellSet(result.GetCellSet(i)); + } + else + { + newResult.AddCellSet(newCells); + } + } + for (vtkm::Id i = 0; i < result.GetNumberOfFields(); ++i) + { + newResult.AddField(result.GetField(i)); + } + + result = newResult; + } + return result; } } diff --git a/vtkm/filter/testing/UnitTestSurfaceNormalsFilter.cxx b/vtkm/filter/testing/UnitTestSurfaceNormalsFilter.cxx index 9a9c62057..2cd4e7d42 100644 --- a/vtkm/filter/testing/UnitTestSurfaceNormalsFilter.cxx +++ b/vtkm/filter/testing/UnitTestSurfaceNormalsFilter.cxx @@ -74,6 +74,7 @@ void TestSurfaceNormals() std::cout << "generate both cell and point normals:\n"; filter.SetGeneratePointNormals(true); + filter.SetAutoOrientNormals(true); result = filter.Execute(ds); VTKM_TEST_ASSERT(result.HasField("Normals", vtkm::cont::Field::Association::POINTS), "Point normals missing.");