Commit 05143f63 authored by Dave Pugmire's avatar Dave Pugmire Committed by Kitware Robot

Merge topic 'coordSysTransform'

74f8885d Move helper classes into detail namespace.
9fd821ed Template on DeviceAdapter.
961f6a58 Add coordinate system transformation.
Acked-by: Kitware Robot's avatarKitware Robot <kwrobot@kitware.com>
Merge-request: !1275
parents 55eafbaf 74f8885d
......@@ -25,6 +25,7 @@ set(headers
CellMeasure.h
Clip.h
ContourTreeUniform.h
CoordinateSystemTransform.h
CosmoTools.h
CrossProduct.h
DispatcherMapField.h
......
//============================================================================
// 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 <vtkm/Math.h>
#include <vtkm/Matrix.h>
#include <vtkm/Transform3D.h>
#include <vtkm/VectorAnalysis.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/WorkletMapField.h>
namespace vtkm
{
namespace worklet
{
namespace detail
{
template <typename T>
struct CylToCar : public vtkm::worklet::WorkletMapField
{
using ControlSignature = void(FieldIn<Vec3>, FieldOut<Vec3>);
using ExecutionSignature = _2(_1);
//Functor
VTKM_EXEC vtkm::Vec<T, 3> operator()(const vtkm::Vec<T, 3>& vec) const
{
vtkm::Vec<T, 3> res(vec[0] * static_cast<T>(vtkm::Cos(vec[1])),
vec[0] * static_cast<T>(vtkm::Sin(vec[1])),
vec[2]);
return res;
}
};
template <typename T>
struct CarToCyl : public vtkm::worklet::WorkletMapField
{
using ControlSignature = void(FieldIn<Vec3>, FieldOut<Vec3>);
using ExecutionSignature = _2(_1);
//Functor
VTKM_EXEC vtkm::Vec<T, 3> operator()(const vtkm::Vec<T, 3>& 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<T>(vtkm::Pi());
else
Theta = vtkm::ASin(vec[1] / R);
vtkm::Vec<T, 3> res(R, Theta, vec[2]);
return res;
}
};
template <typename T>
struct SphereToCar : public vtkm::worklet::WorkletMapField
{
using ControlSignature = void(FieldIn<Vec3>, FieldOut<Vec3>);
using ExecutionSignature = _2(_1);
//Functor
VTKM_EXEC vtkm::Vec<T, 3> operator()(const vtkm::Vec<T, 3>& vec) const
{
T R = vec[0];
T Theta = vec[1];
T Phi = vec[2];
T sinTheta = static_cast<T>(vtkm::Sin(Theta));
T cosTheta = static_cast<T>(vtkm::Cos(Theta));
T sinPhi = static_cast<T>(vtkm::Sin(Phi));
T cosPhi = static_cast<T>(vtkm::Cos(Phi));
T x = R * sinTheta * cosPhi;
T y = R * sinTheta * sinPhi;
T z = R * cosTheta;
vtkm::Vec<T, 3> r(x, y, z);
return r;
}
};
template <typename T>
struct CarToSphere : public vtkm::worklet::WorkletMapField
{
using ControlSignature = void(FieldIn<Vec3>, FieldOut<Vec3>);
using ExecutionSignature = _2(_1);
//Functor
VTKM_EXEC vtkm::Vec<T, 3> operator()(const vtkm::Vec<T, 3>& vec) const
{
T R = vtkm::Sqrt(vtkm::Dot(vec, vec));
T Theta = 0;
if (R > 0)
Theta = vtkm::ACos(vec[2] / R);
T Phi = vtkm::ATan2(vec[1], vec[0]);
if (Phi < 0)
Phi += static_cast<T>(vtkm::TwoPi());
return vtkm::Vec<T, 3>(R, Theta, Phi);
}
};
};
template <typename T>
class CylindricalCoordinateTransform
{
public:
VTKM_CONT
CylindricalCoordinateTransform()
: cartesianToCylindrical(true)
{
}
VTKM_CONT void SetCartesianToCylindrical() { cartesianToCylindrical = true; }
VTKM_CONT void SetCylindricalToCartesian() { cartesianToCylindrical = false; }
template <typename CoordsStorageType, typename DeviceAdapterTag>
void Run(const vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, CoordsStorageType>& inPoints,
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, CoordsStorageType>& outPoints,
DeviceAdapterTag) const
{
if (cartesianToCylindrical)
{
vtkm::worklet::DispatcherMapField<detail::CarToCyl<T>> dispatcher;
dispatcher.Invoke(inPoints, outPoints);
}
else
{
vtkm::worklet::DispatcherMapField<detail::CylToCar<T>> dispatcher;
dispatcher.Invoke(inPoints, outPoints);
}
}
template <typename CoordsStorageType, typename DeviceAdapterTag>
void Run(const vtkm::cont::CoordinateSystem& inPoints,
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, CoordsStorageType>& outPoints,
DeviceAdapterTag) const
{
if (cartesianToCylindrical)
{
vtkm::worklet::DispatcherMapField<detail::CarToCyl<T>> dispatcher;
dispatcher.Invoke(inPoints, outPoints);
}
else
{
vtkm::worklet::DispatcherMapField<detail::CylToCar<T>> dispatcher;
dispatcher.Invoke(inPoints, outPoints);
}
}
private:
bool cartesianToCylindrical;
};
template <typename T>
class SphericalCoordinateTransform
{
public:
VTKM_CONT
SphericalCoordinateTransform()
: CartesianToSpherical(true)
{
}
VTKM_CONT void SetCartesianToSpherical() { CartesianToSpherical = true; }
VTKM_CONT void SetSphericalToCartesian() { CartesianToSpherical = false; }
template <typename CoordsStorageType, typename DeviceAdapterTag>
void Run(const vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, CoordsStorageType>& inPoints,
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, CoordsStorageType>& outPoints,
DeviceAdapterTag) const
{
if (CartesianToSpherical)
{
vtkm::worklet::DispatcherMapField<detail::CarToSphere<T>, DeviceAdapterTag> dispatcher;
dispatcher.Invoke(inPoints, outPoints);
}
else
{
vtkm::worklet::DispatcherMapField<detail::SphereToCar<T>, DeviceAdapterTag> dispatcher;
dispatcher.Invoke(inPoints, outPoints);
}
}
template <typename CoordsStorageType, typename DeviceAdapterTag>
void Run(const vtkm::cont::CoordinateSystem& inPoints,
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, CoordsStorageType>& outPoints,
DeviceAdapterTag) const
{
if (CartesianToSpherical)
{
vtkm::worklet::DispatcherMapField<detail::CarToSphere<T>, DeviceAdapterTag> dispatcher;
dispatcher.Invoke(inPoints, outPoints);
}
else
{
vtkm::worklet::DispatcherMapField<detail::SphereToCar<T>, DeviceAdapterTag> dispatcher;
dispatcher.Invoke(inPoints, outPoints);
}
}
private:
bool CartesianToSpherical;
};
}
} // namespace vtkm::worklet
#endif // vtk_m_worklet_CoordinateSystemTransform_h
......@@ -28,6 +28,7 @@ set(unit_tests
UnitTestCellMeasure.cxx
UnitTestClipping.cxx
UnitTestContourTreeUniform.cxx
UnitTestCoordinateSystemTransform.cxx
UnitTestCosmoTools.cxx
UnitTestCrossProduct.cxx
UnitTestDotProduct.cxx
......
//============================================================================
// 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 <vtkm/cont/CellSetExplicit.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/worklet/CoordinateSystemTransform.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <random>
#include <vector>
namespace
{
std::mt19937 randGenerator;
enum CoordinateType
{
CART = 0,
CYL,
SPH
};
vtkm::cont::DataSet MakeTestDataSet(const CoordinateType& cType)
{
vtkm::cont::DataSet dataSet;
std::vector<vtkm::Vec<vtkm::FloatDefault, 3>> coordinates;
const vtkm::Id dim = 5;
if (cType == CART)
{
for (vtkm::Id j = 0; j < dim; ++j)
{
vtkm::FloatDefault z =
static_cast<vtkm::FloatDefault>(j) / static_cast<vtkm::FloatDefault>(dim - 1);
for (vtkm::Id i = 0; i < dim; ++i)
{
vtkm::FloatDefault x =
static_cast<vtkm::FloatDefault>(i) / static_cast<vtkm::FloatDefault>(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<vtkm::FloatDefault>(j) / static_cast<vtkm::FloatDefault>(dim - 1);
for (vtkm::Id i = 0; i < dim; i++)
{
vtkm::FloatDefault Theta = vtkm::TwoPif() *
(static_cast<vtkm::FloatDefault>(i) / static_cast<vtkm::FloatDefault>(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<float>();
std::vector<vtkm::FloatDefault> Thetas = {
eps, vtkm::Pif() / 4, vtkm::Pif() / 3, vtkm::Pif() / 2, vtkm::Pif() - eps
};
std::vector<vtkm::FloatDefault> 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<vtkm::Id>(
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<vtkm::Vec<vtkm::FloatDefault, 3>>& transform,
const vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>>& doubleTransform,
const std::vector<bool>& 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<vtkm::FloatDefault, 3> p = pointsPortal.Get(i);
vtkm::Vec<vtkm::FloatDefault, 3> r = resultsPortal.Get(i);
bool isEqual = true;
for (vtkm::IdComponent j = 0; j < 3; j++)
{
if (isAngle[static_cast<std::size_t>(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;
using DeviceAdapter = VTKM_DEFAULT_DEVICE_ADAPTER_TAG;
//Test cartesian to cyl
vtkm::cont::DataSet dsCart = MakeTestDataSet(CART);
vtkm::worklet::CylindricalCoordinateTransform<vtkm::FloatDefault> cylTrn;
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>> carToCylPts;
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>> revResult;
cylTrn.SetCartesianToCylindrical();
cylTrn.Run(dsCart.GetCoordinateSystem(), carToCylPts, DeviceAdapter());
cylTrn.SetCylindricalToCartesian();
cylTrn.Run(carToCylPts, revResult, DeviceAdapter());
ValidateCoordTransform(
dsCart.GetCoordinateSystem(), carToCylPts, revResult, { false, false, false });
//Test cylindrical to cartesian
vtkm::cont::DataSet dsCyl = MakeTestDataSet(CYL);
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>> cylToCarPts;
cylTrn.SetCylindricalToCartesian();
cylTrn.Run(dsCyl.GetCoordinateSystem(), cylToCarPts, DeviceAdapter());
cylTrn.SetCartesianToCylindrical();
cylTrn.Run(cylToCarPts, revResult, DeviceAdapter());
ValidateCoordTransform(
dsCyl.GetCoordinateSystem(), cylToCarPts, revResult, { false, true, false });
//Spherical transform
//Test cartesian to sph
vtkm::worklet::SphericalCoordinateTransform<vtkm::FloatDefault> sphTrn;
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>> carToSphPts;
sphTrn.SetCartesianToSpherical();
sphTrn.Run(dsCart.GetCoordinateSystem(), carToSphPts, DeviceAdapter());
sphTrn.SetSphericalToCartesian();
sphTrn.Run(carToSphPts, revResult, DeviceAdapter());
ValidateCoordTransform(
dsCart.GetCoordinateSystem(), carToSphPts, revResult, { false, true, true });
//Test spherical to cartesian
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>> sphToCarPts;
vtkm::cont::DataSet dsSph = MakeTestDataSet(SPH);
sphTrn.SetSphericalToCartesian();
sphTrn.Run(dsSph.GetCoordinateSystem(), sphToCarPts, DeviceAdapter());
sphTrn.SetCartesianToSpherical();
sphTrn.Run(sphToCarPts, revResult, DeviceAdapter());
ValidateCoordTransform(
dsSph.GetCoordinateSystem(), sphToCarPts, revResult, { false, true, true });
sphTrn.SetSphericalToCartesian();
sphTrn.Run(dsSph.GetCoordinateSystem(), sphToCarPts, DeviceAdapter());
sphTrn.SetCartesianToSpherical();
sphTrn.Run(sphToCarPts, revResult, DeviceAdapter());
ValidateCoordTransform(
dsSph.GetCoordinateSystem(), sphToCarPts, revResult, { false, true, true });
}
int UnitTestCoordinateSystemTransform(int, char* [])
{
return vtkm::cont::testing::Testing::Run(TestCoordinateSystemTransform);
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment