From 961f6a585124b7d75ae12c566369ecc341307caa Mon Sep 17 00:00:00 2001 From: Dave Pugmire Date: Wed, 20 Jun 2018 16:40:47 -0400 Subject: [PATCH] Add coordinate system transformation. --- vtkm/worklet/CMakeLists.txt | 1 + vtkm/worklet/CoordinateSystemTransform.h | 224 ++++++++++++++++++ vtkm/worklet/testing/CMakeLists.txt | 1 + .../UnitTestCoordinateSystemTransform.cxx | 214 +++++++++++++++++ 4 files changed, 440 insertions(+) create mode 100644 vtkm/worklet/CoordinateSystemTransform.h create mode 100644 vtkm/worklet/testing/UnitTestCoordinateSystemTransform.cxx diff --git a/vtkm/worklet/CMakeLists.txt b/vtkm/worklet/CMakeLists.txt index d88af1cfd..ac3120821 100644 --- a/vtkm/worklet/CMakeLists.txt +++ b/vtkm/worklet/CMakeLists.txt @@ -25,6 +25,7 @@ set(headers CellMeasure.h Clip.h ContourTreeUniform.h + CoordinateSystemTransform.h CosmoTools.h CrossProduct.h DispatcherMapField.h diff --git a/vtkm/worklet/CoordinateSystemTransform.h b/vtkm/worklet/CoordinateSystemTransform.h new file mode 100644 index 000000000..f8c178c37 --- /dev/null +++ b/vtkm/worklet/CoordinateSystemTransform.h @@ -0,0 +1,224 @@ +//============================================================================ +// 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 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2014 UT-Battelle, LLC. +// Copyright 2014 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 vtk_m_worklet_CoordinateSystemTransform_h +#define vtk_m_worklet_CoordinateSystemTransform_h + +#include +#include +#include +#include +#include + +namespace vtkm +{ +namespace worklet +{ + +template +class CylindricalCoordinateTransform +{ +public: + struct CylToCar : public vtkm::worklet::WorkletMapField + { + using ControlSignature = void(FieldIn, FieldOut); + using ExecutionSignature = _2(_1); + + + //Functor + VTKM_EXEC vtkm::Vec operator()(const vtkm::Vec& vec) const + { + vtkm::Vec res(vec[0] * static_cast(vtkm::Cos(vec[1])), + vec[0] * static_cast(vtkm::Sin(vec[1])), + vec[2]); + return res; + } + }; + + struct CarToCyl : public vtkm::worklet::WorkletMapField + { + using ControlSignature = void(FieldIn, FieldOut); + using ExecutionSignature = _2(_1); + + //Functor + VTKM_EXEC vtkm::Vec operator()(const vtkm::Vec& vec) const + { + T R = vtkm::Sqrt(vec[0] * vec[0] + vec[1] * vec[1]); + T Theta = 0; + + if (vec[0] == 0 && vec[1] == 0) + Theta = 0; + else if (vec[0] < 0) + Theta = -vtkm::ASin(vec[1] / R) + static_cast(vtkm::Pi()); + else + Theta = vtkm::ASin(vec[1] / R); + + vtkm::Vec res(R, Theta, vec[2]); + return res; + } + }; + + VTKM_CONT + CylindricalCoordinateTransform() + : cartesianToCylindrical(true) + { + } + + VTKM_CONT void SetCartesianToCylindrical() { cartesianToCylindrical = true; } + VTKM_CONT void SetCylindricalToCartesian() { cartesianToCylindrical = false; } + + template + void Run(const vtkm::cont::ArrayHandle, CoordsStorageType>& inPoints, + vtkm::cont::ArrayHandle, CoordsStorageType>& outPoints) + { + if (cartesianToCylindrical) + { + vtkm::worklet::DispatcherMapField dispatcher; + dispatcher.Invoke(inPoints, outPoints); + } + else + { + vtkm::worklet::DispatcherMapField dispatcher; + dispatcher.Invoke(inPoints, outPoints); + } + } + + template + void Run(const vtkm::cont::CoordinateSystem& inPoints, + vtkm::cont::ArrayHandle, CoordsStorageType>& outPoints) + { + if (cartesianToCylindrical) + { + vtkm::worklet::DispatcherMapField dispatcher; + dispatcher.Invoke(inPoints, outPoints); + } + else + { + vtkm::worklet::DispatcherMapField dispatcher; + dispatcher.Invoke(inPoints, outPoints); + } + } + +private: + bool cartesianToCylindrical; +}; + +template +class SphericalCoordinateTransform +{ +public: + struct SphereToCar : public vtkm::worklet::WorkletMapField + { + using ControlSignature = void(FieldIn, FieldOut); + using ExecutionSignature = _2(_1); + + //Functor + VTKM_EXEC vtkm::Vec operator()(const vtkm::Vec& vec) const + { + T R = vec[0]; + T Theta = vec[1]; + T Phi = vec[2]; + + T sinTheta = static_cast(vtkm::Sin(Theta)); + T cosTheta = static_cast(vtkm::Cos(Theta)); + T sinPhi = static_cast(vtkm::Sin(Phi)); + T cosPhi = static_cast(vtkm::Cos(Phi)); + + T x = R * sinTheta * cosPhi; + T y = R * sinTheta * sinPhi; + T z = R * cosTheta; + + vtkm::Vec r(x, y, z); + return r; + } + }; + + struct CarToSphere : public vtkm::worklet::WorkletMapField + { + using ControlSignature = void(FieldIn, FieldOut); + using ExecutionSignature = _2(_1); + + //Functor + VTKM_EXEC vtkm::Vec operator()(const vtkm::Vec& vec) const + { + T x = vec[0]; + T y = vec[1]; + T z = vec[2]; + + T R = vtkm::Sqrt(x * x + y * y + z * z); + T Theta = 0; + if (R > 0) + Theta = vtkm::ACos(z / R); + T Phi = vtkm::ATan2(y, x); + if (Phi < 0) + Phi += static_cast(vtkm::TwoPi()); + + return vtkm::Vec(R, Theta, Phi); + } + }; + + VTKM_CONT + SphericalCoordinateTransform() + : CartesianToSpherical(true) + { + } + + VTKM_CONT void SetCartesianToSpherical() { CartesianToSpherical = true; } + VTKM_CONT void SetSphericalToCartesian() { CartesianToSpherical = false; } + + template + void Run(const vtkm::cont::ArrayHandle, CoordsStorageType>& inPoints, + vtkm::cont::ArrayHandle, CoordsStorageType>& outPoints) + { + if (CartesianToSpherical) + { + vtkm::worklet::DispatcherMapField dispatcher; + dispatcher.Invoke(inPoints, outPoints); + } + else + { + vtkm::worklet::DispatcherMapField dispatcher; + dispatcher.Invoke(inPoints, outPoints); + } + } + + template + void Run(const vtkm::cont::CoordinateSystem& inPoints, + vtkm::cont::ArrayHandle, CoordsStorageType>& outPoints) + { + if (CartesianToSpherical) + { + vtkm::worklet::DispatcherMapField dispatcher; + dispatcher.Invoke(inPoints, outPoints); + } + else + { + vtkm::worklet::DispatcherMapField dispatcher; + dispatcher.Invoke(inPoints, outPoints); + } + } + +private: + bool CartesianToSpherical; +}; +} +} // namespace vtkm::worklet + +#endif // vtk_m_worklet_CoordinateSystemTransform_h diff --git a/vtkm/worklet/testing/CMakeLists.txt b/vtkm/worklet/testing/CMakeLists.txt index e63d7ca55..e1a9f88c7 100644 --- a/vtkm/worklet/testing/CMakeLists.txt +++ b/vtkm/worklet/testing/CMakeLists.txt @@ -28,6 +28,7 @@ set(unit_tests UnitTestCellMeasure.cxx UnitTestClipping.cxx UnitTestContourTreeUniform.cxx + UnitTestCoordinateSystemTransform.cxx UnitTestCosmoTools.cxx UnitTestCrossProduct.cxx UnitTestDotProduct.cxx diff --git a/vtkm/worklet/testing/UnitTestCoordinateSystemTransform.cxx b/vtkm/worklet/testing/UnitTestCoordinateSystemTransform.cxx new file mode 100644 index 000000000..db8c102c8 --- /dev/null +++ b/vtkm/worklet/testing/UnitTestCoordinateSystemTransform.cxx @@ -0,0 +1,214 @@ +//============================================================================ +// 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 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2014 UT-Battelle, LLC. +// Copyright 2014 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 + +namespace +{ +std::mt19937 randGenerator; + +enum CoordinateType +{ + CART = 0, + CYL, + SPH +}; + +vtkm::cont::DataSet MakeTestDataSet(const CoordinateType& cType) +{ + vtkm::cont::DataSet dataSet; + + std::vector> coordinates; + const vtkm::Id dim = 5; + if (cType == CART) + { + for (vtkm::Id j = 0; j < dim; ++j) + { + vtkm::FloatDefault z = + static_cast(j) / static_cast(dim - 1); + for (vtkm::Id i = 0; i < dim; ++i) + { + vtkm::FloatDefault x = + static_cast(i) / static_cast(dim - 1); + vtkm::FloatDefault y = (x * x + z * z) / 2.0f; + coordinates.push_back(vtkm::make_Vec(x + 0, y + 0, z + 0)); + } + } + } + else if (cType == CYL) + { + vtkm::FloatDefault R = 1.0f; + for (vtkm::Id j = 0; j < dim; j++) + { + vtkm::FloatDefault Z = + static_cast(j) / static_cast(dim - 1); + for (vtkm::Id i = 0; i < dim; i++) + { + vtkm::FloatDefault Theta = vtkm::TwoPif() * + (static_cast(i) / static_cast(dim - 1)); + coordinates.push_back(vtkm::make_Vec(R, Theta, Z)); + } + } + } + else if (cType == SPH) + { + //Spherical coordinates have some degenerate cases, so provide some good cases. + vtkm::FloatDefault R = 1.0f; + vtkm::FloatDefault eps = vtkm::Epsilon(); + std::vector Thetas = { + eps, vtkm::Pif() / 4, vtkm::Pif() / 3, vtkm::Pif() / 2, vtkm::Pif() - eps + }; + std::vector Phis = { + eps, vtkm::TwoPif() / 4, vtkm::TwoPif() / 3, vtkm::TwoPif() / 2, vtkm::TwoPif() - eps + }; + for (std::size_t i = 0; i < Thetas.size(); i++) + for (std::size_t j = 0; j < Phis.size(); j++) + coordinates.push_back(vtkm::make_Vec(R, Thetas[i], Phis[j])); + } + + vtkm::Id numCells = (dim - 1) * (dim - 1); + dataSet.AddCoordinateSystem( + vtkm::cont::make_CoordinateSystem("coordinates", coordinates, vtkm::CopyFlag::On)); + + vtkm::cont::CellSetExplicit<> cellSet("cells"); + cellSet.PrepareToAddCells(numCells, numCells * 4); + for (vtkm::Id j = 0; j < dim - 1; ++j) + { + for (vtkm::Id i = 0; i < dim - 1; ++i) + { + cellSet.AddCell(vtkm::CELL_SHAPE_QUAD, + 4, + vtkm::make_Vec( + j * dim + i, j * dim + i + 1, (j + 1) * dim + i + 1, (j + 1) * dim + i)); + } + } + cellSet.CompleteAddingCells(vtkm::Id(coordinates.size())); + + dataSet.AddCellSet(cellSet); + return dataSet; +} + +void ValidateCoordTransform( + const vtkm::cont::CoordinateSystem& coords, + const vtkm::cont::ArrayHandle>& transform, + const vtkm::cont::ArrayHandle>& doubleTransform, + const std::vector& isAngle) +{ + auto points = coords.GetData(); + VTKM_TEST_ASSERT(points.GetNumberOfValues() == transform.GetNumberOfValues() && + points.GetNumberOfValues() == doubleTransform.GetNumberOfValues(), + "Incorrect number of points in point transform"); + + //The double transform should produce the same result. + auto pointsPortal = points.GetPortalConstControl(); + auto resultsPortal = doubleTransform.GetPortalConstControl(); + + for (vtkm::Id i = 0; i < points.GetNumberOfValues(); i++) + { + vtkm::Vec p = pointsPortal.Get(i); + vtkm::Vec r = resultsPortal.Get(i); + bool isEqual = true; + for (vtkm::IdComponent j = 0; j < 3; j++) + { + if (isAngle[static_cast(j)]) + isEqual &= (test_equal(p[j], r[j]) || test_equal(p[j] + vtkm::TwoPif(), r[j]) || + test_equal(p[j], r[j] + vtkm::TwoPif())); + else + isEqual &= test_equal(p[j], r[j]); + } + VTKM_TEST_ASSERT(isEqual, "Wrong result for PointTransform worklet"); + } +} +} + +void TestCoordinateSystemTransform() +{ + std::cout << "Testing CylindricalCoordinateTransform Worklet" << std::endl; + + //Test cartesian to cyl + vtkm::cont::DataSet dsCart = MakeTestDataSet(CART); + vtkm::worklet::CylindricalCoordinateTransform cylTrn; + + vtkm::cont::ArrayHandle> carToCylPts; + vtkm::cont::ArrayHandle> revResult; + + cylTrn.SetCartesianToCylindrical(); + cylTrn.Run(dsCart.GetCoordinateSystem(), carToCylPts); + + cylTrn.SetCylindricalToCartesian(); + cylTrn.Run(carToCylPts, revResult); + ValidateCoordTransform( + dsCart.GetCoordinateSystem(), carToCylPts, revResult, { false, false, false }); + + //Test cylindrical to cartesian + vtkm::cont::DataSet dsCyl = MakeTestDataSet(CYL); + vtkm::cont::ArrayHandle> cylToCarPts; + cylTrn.SetCylindricalToCartesian(); + cylTrn.Run(dsCyl.GetCoordinateSystem(), cylToCarPts); + + cylTrn.SetCartesianToCylindrical(); + cylTrn.Run(cylToCarPts, revResult); + ValidateCoordTransform( + dsCyl.GetCoordinateSystem(), cylToCarPts, revResult, { false, true, false }); + + //Spherical transform + //Test cartesian to sph + vtkm::worklet::SphericalCoordinateTransform sphTrn; + vtkm::cont::ArrayHandle> carToSphPts; + + sphTrn.SetCartesianToSpherical(); + sphTrn.Run(dsCart.GetCoordinateSystem(), carToSphPts); + + sphTrn.SetSphericalToCartesian(); + sphTrn.Run(carToSphPts, revResult); + ValidateCoordTransform( + dsCart.GetCoordinateSystem(), carToSphPts, revResult, { false, true, true }); + + //Test spherical to cartesian + vtkm::cont::ArrayHandle> sphToCarPts; + vtkm::cont::DataSet dsSph = MakeTestDataSet(SPH); + + sphTrn.SetSphericalToCartesian(); + sphTrn.Run(dsSph.GetCoordinateSystem(), sphToCarPts); + + sphTrn.SetCartesianToSpherical(); + sphTrn.Run(sphToCarPts, revResult); + + ValidateCoordTransform( + dsSph.GetCoordinateSystem(), sphToCarPts, revResult, { false, true, true }); + sphTrn.SetSphericalToCartesian(); + sphTrn.Run(dsSph.GetCoordinateSystem(), sphToCarPts); + sphTrn.SetCartesianToSpherical(); + sphTrn.Run(sphToCarPts, revResult); + ValidateCoordTransform( + dsSph.GetCoordinateSystem(), sphToCarPts, revResult, { false, true, true }); +} + +int UnitTestCoordinateSystemTransform(int, char* []) +{ + return vtkm::cont::testing::Testing::Run(TestCoordinateSystemTransform); +}