From 1147edb192645a9b8727461343c7e49ace9013df Mon Sep 17 00:00:00 2001 From: Robert Maynard Date: Tue, 8 Aug 2017 16:31:18 -0400 Subject: [PATCH] MarchingCubes now uses Gradient fast paths when possible. When running on image data we now use central differences to compute the point gradient --- vtkm/exec/ConnectivityStructured.h | 8 ++ .../testing/UnitTestMarchingCubesFilter.cxx | 20 ++++- vtkm/worklet/CMakeLists.txt | 2 +- vtkm/worklet/MarchingCubes.h | 78 ++++++++++++++++++- vtkm/worklet/contour/CMakeLists.txt | 26 +++++++ .../DataTables.h} | 0 .../gradient/StructuredPointGradient.h | 15 ++-- 7 files changed, 135 insertions(+), 14 deletions(-) create mode 100644 vtkm/worklet/contour/CMakeLists.txt rename vtkm/worklet/{MarchingCubesDataTables.h => contour/DataTables.h} (100%) diff --git a/vtkm/exec/ConnectivityStructured.h b/vtkm/exec/ConnectivityStructured.h index 31e186258..04dd74db7 100644 --- a/vtkm/exec/ConnectivityStructured.h +++ b/vtkm/exec/ConnectivityStructured.h @@ -62,6 +62,12 @@ public: { } + VTKM_EXEC_CONT + ConnectivityStructured(const ConnectivityStructured& src) + : Internals(src.Internals) + { + } + VTKM_EXEC vtkm::Id GetNumberOfElements() const { return Helper::GetNumberOfElements(this->Internals); } @@ -113,6 +119,8 @@ public: return this->Internals.GetPointDimensions(); } + friend class ConnectivityStructured; + private: InternalsType Internals; }; diff --git a/vtkm/filter/testing/UnitTestMarchingCubesFilter.cxx b/vtkm/filter/testing/UnitTestMarchingCubesFilter.cxx index fe41e694b..56ec5cc8a 100644 --- a/vtkm/filter/testing/UnitTestMarchingCubesFilter.cxx +++ b/vtkm/filter/testing/UnitTestMarchingCubesFilter.cxx @@ -429,7 +429,8 @@ void TestNormals(const vtkm::cont::DataSet& dataset, bool structured) { const vtkm::Id numVerts = 16; - const vtkm::Vec hq[numVerts] = { + //Calculated using PointGradient + const vtkm::Vec hq_ug[numVerts] = { { 0.1510f, 0.6268f, 0.7644f }, { 0.1333f, -0.3974f, 0.9079f }, { 0.1626f, 0.7642f, 0.6242f }, { 0.3853f, 0.6643f, 0.6405f }, { -0.1337f, 0.7136f, 0.6876f }, { 0.7705f, -0.4212f, 0.4784f }, @@ -440,6 +441,19 @@ void TestNormals(const vtkm::cont::DataSet& dataset, bool structured) { 0.1234f, -0.8871f, -0.4448f }, { 0.1333f, -0.3974f, -0.9079f } }; + //Calculated using StructuredPointGradient + const vtkm::Vec hq_sg[numVerts] = { + { 0.165519f, 0.687006f, 0.707549f }, { 0.188441f, -0.561729f, 0.805574f }, + { 0.179543f, 0.702158f, 0.689012f }, { 0.271085f, 0.692957f, 0.668074f }, + { 0.00313049f, 0.720109f, 0.693854f }, { 0.549947f, -0.551974f, 0.626804f }, + { -0.447526f, -0.588187f, 0.673614f }, { 0.167553f, -0.779396f, 0.603711f }, + { 0.179543f, 0.702158f, -0.689012f }, { 0.271085f, 0.692957f, -0.668074f }, + { 0.00313049f, 0.720109f, -0.693854f }, { 0.165519f, 0.687006f, -0.707549f }, + { 0.549947f, -0.551974f, -0.626804f }, { -0.447526f, -0.588187f, -0.673614f }, + { 0.167553f, -0.779396f, -0.603711f }, { 0.188441f, -0.561729f, -0.805574f } + }; + + //Calculated using normals of the output triangles const vtkm::Vec fast[numVerts] = { { -0.1351f, 0.4377f, 0.8889f }, { 0.2863f, -0.1721f, 0.9426f }, { 0.3629f, 0.8155f, 0.4509f }, { 0.8486f, 0.3560f, 0.3914f }, @@ -458,7 +472,7 @@ void TestNormals(const vtkm::cont::DataSet& dataset, bool structured) mc.SetGenerateNormals(true); // Test default normals generation: high quality for structured, fast for unstructured. - auto expected = structured ? hq : fast; + auto expected = structured ? hq_sg : fast; auto result = mc.Execute(dataset, dataset.GetField("pointvar")); result.GetDataSet().GetField("normals").GetData().CopyTo(normals); @@ -479,7 +493,7 @@ void TestNormals(const vtkm::cont::DataSet& dataset, bool structured) else { mc.SetComputeFastNormalsForUnstructured(false); - expected = hq; + expected = hq_ug; } result = mc.Execute(dataset, dataset.GetField("pointvar")); diff --git a/vtkm/worklet/CMakeLists.txt b/vtkm/worklet/CMakeLists.txt index 54a47a93b..6f062cb54 100644 --- a/vtkm/worklet/CMakeLists.txt +++ b/vtkm/worklet/CMakeLists.txt @@ -43,7 +43,6 @@ set(headers Keys.h Magnitude.h MarchingCubes.h - MarchingCubesDataTables.h Mask.h MaskPoints.h NDimsEntropy.h @@ -73,6 +72,7 @@ set(headers #----------------------------------------------------------------------------- add_subdirectory(internal) +add_subdirectory(contour) add_subdirectory(contourtree) add_subdirectory(cosmotools) add_subdirectory(gradient) diff --git a/vtkm/worklet/MarchingCubes.h b/vtkm/worklet/MarchingCubes.h index c1249ff7a..d4796ccc8 100644 --- a/vtkm/worklet/MarchingCubes.h +++ b/vtkm/worklet/MarchingCubes.h @@ -41,15 +41,18 @@ #include #include +#include #include #include #include #include #include +#include #include -#include -#include +#include +#include +#include namespace vtkm { @@ -526,10 +529,40 @@ public: const WholeFieldIn& inputField, NormalType& normal) const { - vtkm::worklet::gradient::PointGradient gradient; + using T = typename WholeFieldIn::ValueType; + vtkm::worklet::gradient::PointGradient gradient; gradient(numCells, cellIds, pointId, geometry, pointCoordinates, inputField, normal); } + template + VTKM_EXEC void operator()(const vtkm::IdComponent& vtkmNotUsed(numCells), + const FromIndexType& vtkmNotUsed(cellIds), + vtkm::Id pointId, + vtkm::exec::ConnectivityStructured& geometry, + const WholeCoordinatesIn& pointCoordinates, + const WholeFieldIn& inputField, + NormalType& normal) const + { + using T = typename WholeFieldIn::ValueType; + + //Optimization for structured cellsets so we can call StructuredPointGradient + //and have way faster gradients + vtkm::exec::ConnectivityStructured pointGeom(geometry); + vtkm::exec::arg::ThreadIndicesPointNeighborhood<3> tpn(pointId, pointId, 0, pointGeom, 0); + + const auto& boundary = tpn.GetBoundaryState(); + auto pointPortal = pointCoordinates.GetPortal(); + auto fieldPortal = inputField.GetPortal(); + vtkm::exec::arg::Neighborhood<1, decltype(pointPortal)> points(pointPortal, boundary); + vtkm::exec::arg::Neighborhood<1, decltype(fieldPortal)> field(fieldPortal, boundary); + + vtkm::worklet::gradient::StructuredPointGradient gradient; + gradient(boundary, points, field, normal); + } + ScatterType GetScatter() const { return this->Scatter; } private: @@ -577,7 +610,8 @@ public: const WholeWeightsIn& weights, NormalType& normal) const { - vtkm::worklet::gradient::PointGradient gradient; + using T = typename WholeFieldIn::ValueType; + vtkm::worklet::gradient::PointGradient gradient; NormalType grad1; gradient(numCells, cellIds, pointId, geometry, pointCoordinates, inputField, grad1); @@ -586,6 +620,42 @@ public: normal = vtkm::Normal(vtkm::Lerp(grad0, grad1, weight)); } + template + VTKM_EXEC void operator()(const vtkm::IdComponent& vtkmNotUsed(numCells), + const FromIndexType& vtkmNotUsed(cellIds), + vtkm::Id pointId, + vtkm::exec::ConnectivityStructured& geometry, + const WholeCoordinatesIn& pointCoordinates, + const WholeFieldIn& inputField, + vtkm::Id edgeId, + const WholeWeightsIn& weights, + NormalType& normal) const + { + using T = typename WholeFieldIn::ValueType; + //Optimization for structured cellsets so we can call StructuredPointGradient + //and have way faster gradients + vtkm::exec::ConnectivityStructured pointGeom(geometry); + vtkm::exec::arg::ThreadIndicesPointNeighborhood<3> tpn(pointId, pointId, 0, pointGeom, 0); + + const auto& boundary = tpn.GetBoundaryState(); + auto pointPortal = pointCoordinates.GetPortal(); + auto fieldPortal = inputField.GetPortal(); + vtkm::exec::arg::Neighborhood<1, decltype(pointPortal)> points(pointPortal, boundary); + vtkm::exec::arg::Neighborhood<1, decltype(fieldPortal)> field(fieldPortal, boundary); + + vtkm::worklet::gradient::StructuredPointGradient gradient; + NormalType grad1; + gradient(boundary, points, field, grad1); + + NormalType grad0 = normal; + auto weight = weights.Get(edgeId); + normal = vtkm::Normal(vtkm::Lerp(grad0, grad1, weight)); + } + ScatterType GetScatter() const { return this->Scatter; } private: diff --git a/vtkm/worklet/contour/CMakeLists.txt b/vtkm/worklet/contour/CMakeLists.txt new file mode 100644 index 000000000..e54d31188 --- /dev/null +++ b/vtkm/worklet/contour/CMakeLists.txt @@ -0,0 +1,26 @@ +##============================================================================ +## 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 2016 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +## Copyright 2016 UT-Battelle, LLC. +## Copyright 2016 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. +##============================================================================ + +set(headers + DataTables.h + ) + +#----------------------------------------------------------------------------- +vtkm_declare_headers(${headers}) diff --git a/vtkm/worklet/MarchingCubesDataTables.h b/vtkm/worklet/contour/DataTables.h similarity index 100% rename from vtkm/worklet/MarchingCubesDataTables.h rename to vtkm/worklet/contour/DataTables.h diff --git a/vtkm/worklet/gradient/StructuredPointGradient.h b/vtkm/worklet/gradient/StructuredPointGradient.h index 22ed8ef6e..fade0ba06 100644 --- a/vtkm/worklet/gradient/StructuredPointGradient.h +++ b/vtkm/worklet/gradient/StructuredPointGradient.h @@ -58,6 +58,7 @@ struct StructuredPointGradient : public vtkm::worklet::WorkletPointNeighborhood3 { using CoordType = typename PointsIn::ValueType; using CT = typename vtkm::BaseComponent::Type; + using OT = typename GradientOutType::ComponentType; vtkm::Vec xi, eta, zeta; this->Jacobian(inputPoints, boundary, xi, eta, zeta); //store the metrics in xi,eta,zeta @@ -70,9 +71,9 @@ struct StructuredPointGradient : public vtkm::worklet::WorkletPointNeighborhood3 deta = (boundary.OnY() ? deta : deta * 0.5f); dzeta = (boundary.OnZ() ? dzeta : dzeta * 0.5f); - outputGradient[0] = static_cast(xi[0] * dxi + eta[0] * deta + zeta[0] * dzeta); - outputGradient[1] = static_cast(xi[1] * dxi + eta[1] * deta + zeta[1] * dzeta); - outputGradient[2] = static_cast(xi[2] * dxi + eta[2] * deta + zeta[2] * dzeta); + outputGradient[0] = static_cast(xi[0] * dxi + eta[0] * deta + zeta[0] * dzeta); + outputGradient[1] = static_cast(xi[1] * dxi + eta[1] * deta + zeta[1] * dzeta); + outputGradient[2] = static_cast(xi[2] * dxi + eta[2] * deta + zeta[2] * dzeta); } template @@ -89,6 +90,8 @@ struct StructuredPointGradient : public vtkm::worklet::WorkletPointNeighborhood3 using PointsIn = vtkm::exec::arg::Neighborhood<1, vtkm::internal::ArrayPortalUniformPointCoordinates>; using CoordType = typename PointsIn::ValueType; + using OT = typename GradientOutType::ComponentType; + CoordType r = inputPoints.Portal.GetSpacing(); @@ -100,9 +103,9 @@ struct StructuredPointGradient : public vtkm::worklet::WorkletPointNeighborhood3 const T dy = inputField.Get(0, 1, 0) - inputField.Get(0, -1, 0); const T dz = inputField.Get(0, 0, 1) - inputField.Get(0, 0, -1); - outputGradient[0] = dx * r[0]; - outputGradient[1] = dy * r[1]; - outputGradient[2] = dz * r[2]; + outputGradient[0] = static_cast(dx * r[0]); + outputGradient[1] = static_cast(dy * r[1]); + outputGradient[2] = static_cast(dz * r[2]); } //we need to pass the coordinates into this function, and instead