Expose vtkm coordinate tranform filter

Accelerators/Vtkm/CMakeLists.txt

@@ -36,6 +36,7 @@ set(headers
   vtkmCleanGrid.h
   vtkmClip.h
   vtkmContour.h
+  vtkmCoordinateSystemTransform.h
   vtkmExternalFaces.h
   vtkmExtractVOI.h
   vtkmGradient.h
@@ -61,6 +62,7 @@ set(cpu_accelerator_srcs
   vtkmClip.cxx
   vtkmConnectivityExec.cxx
   vtkmContour.cxx
+  vtkmCoordinateSystemTransform.cxx
   vtkmExternalFaces.cxx
   vtkmExtractVOI.cxx
   vtkmGradient.cxx
@@ -89,6 +91,7 @@ set(cuda_accelerator_srcs
   vtkmClip.cu
   vtkmConnectivityExec.cu
   vtkmContour.cu
+  vtkmCoordinateSystemTransform.cu
   vtkmExternalFaces.cu
   vtkmExtractVOI.cu
   vtkmGradient.cu

Accelerators/Vtkm/Testing/Cxx/CMakeLists.txt

 vtk_add_test_cxx(vtkAcceleratorsVtkmCxxTests tests
   TestVTKMCleanGrid.cxx
+  TestVTKMCoordinateSystemTransform.cxx,NO_VALID
   TestVTKMClip.cxx
   TestVTKMClipWithImplicitFunction.cxx
   TestVTKMExternalFaces.cxx

Accelerators/Vtkm/Testing/Cxx/TestVTKMCoordinateSystemTransform.cxx

+/*=========================================================================
+
+  Program:   Visualization Toolkit
+  Module:    TestVTKMClip.cxx
+
+  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
+  All rights reserved.
+  See Copyright.txt or http://www.kitware.com/Copyright.htm 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.
+
+=========================================================================*/
+#include "vtkmCoordinateSystemTransform.h"
+
+#include "vtkDoubleArray.h"
+#include "vtkPoints.h"
+#include "vtkPolyData.h"
+#include "vtkSmartPointer.h"
+
+namespace {
+const double pi = 3.14159265358979323846264338327950288;
+const double twoPi = 6.28318530717958647692528676655900576;
+const int dim = 5;
+const double eps = 0.00001;
+const double tolerance = 0.0001;
+
+bool ArePointsWithinTolerance(double v1, double v2)
+{
+  if(v1 == v2 || fabs(v1)+fabs(v2) < tolerance)
+  {
+    return true;
+  }
+
+  if(v1 == 0.0)
+  {
+    if(fabs(v2) < tolerance)
+    {
+      return true;
+    }
+    return false;
+  }
+  if(fabs(fabs(v1) - fabs(v2)) < tolerance)
+  {
+    return true;
+  }
+  return false;
+}
+
+enum struct CoordinateType
+{
+  CART = 0,
+  CYL,
+  SPH
+};
+
+void MakeTestDataSet(vtkPolyData* pd, const CoordinateType& coordType)
+{
+  vtkSmartPointer<vtkDoubleArray> pcoords = vtkSmartPointer<vtkDoubleArray>::New();
+  pcoords->SetNumberOfComponents(3);
+  pcoords->SetNumberOfTuples(dim * dim);
+
+  if (coordType == CoordinateType::CART)
+  {
+    for (vtkIdType i = 0; i < dim; ++i)
+    {
+      double z =
+        static_cast<double>(i) / static_cast<double>(dim - 1);
+      for (vtkIdType j = 0; j < dim; ++j)
+      {
+        double x =
+          static_cast<double>(j) / static_cast<double>(dim - 1);
+        double y = (x * x + z * z) / 2.0f;
+        pcoords->SetTuple3(i*dim + j,x + 0, y + 0, z + 0);
+      }
+    }
+  }
+  else if (coordType == CoordinateType::CYL)
+  {
+    double R = 1.0f;
+    for (vtkIdType i = 0; i < dim; i++)
+    {
+      double Z =
+        static_cast<double>(i) / static_cast<double>(dim - 1);
+      for (vtkIdType j = 0; j < dim; j++)
+      {
+        double Theta = twoPi *
+          (static_cast<double>(j) / static_cast<double>(dim - 1));
+        pcoords->SetTuple3(i*dim + j, R, Theta, Z);
+      }
+    }
+  }
+  else if (coordType == CoordinateType::SPH)
+  {
+    //Spherical coordinates have some degenerate cases, so provide some good cases.
+    double R = 1.0f;
+    std::vector<double> Thetas = {
+      eps, pi / 4, pi / 3, pi / 2, pi - eps
+    };
+    std::vector<double> Phis = {
+      eps, twoPi / 4, twoPi / 3, twoPi / 2, twoPi - eps
+    };
+    for (std::size_t i = 0; i < Thetas.size(); i++)
+    {
+      for (std::size_t j = 0; j < Phis.size(); j++)
+      {
+        pcoords->SetTuple3(static_cast<vtkIdType>(i*dim + j), R, Thetas[i], Phis[j]);
+      }
+    }
+  }
+  pd->GetPoints()->SetData(pcoords);
+}
+
+void ValidateCoordTransform(vtkPolyData* pd, vtkPolyData* pdTrans,
+                            const std::vector<bool>& isAngle)
+{
+  vtkPoints* pdPoints = pd->GetPoints();
+  vtkPoints* pdTransPoints = pdTrans->GetPoints();
+  assert(pdPoints->GetNumberOfPoints() == pdTransPoints->GetNumberOfPoints());
+  for (vtkIdType i=0; i < pdPoints->GetNumberOfPoints(); i++)
+  {
+    double* point = pdPoints->GetPoint(i);
+    double* pointTrans = pdTransPoints->GetPoint(i);
+    bool isEqual = true;
+    for (size_t j=0; j <3; j++)
+    {
+      if (isAngle[j])
+      {
+        isEqual &= (ArePointsWithinTolerance(point[j], pointTrans[j]) ||
+                   ArePointsWithinTolerance(point[j] + static_cast<double>(twoPi), pointTrans[j]) ||
+                   ArePointsWithinTolerance(point[j], pointTrans[j] + static_cast<double>(twoPi)));
+      }
+      else
+      {
+        isEqual &= ArePointsWithinTolerance(point[j], pointTrans[j]);
+      }
+      if (isEqual==false)
+      {
+        std::cerr << "i=" << i << " is wrong! result value=" <<  pointTrans[j] <<
+                  " target value=" << point[j] <<std::endl;
+      }
+    }
+    assert (isEqual == true);
+  }
+}
+}
+
+int TestVTKMCoordinateSystemTransform(int, char*[])
+{
+  // Test cartesian to cylindrical
+  vtkSmartPointer<vtkPolyData> pdCart = vtkSmartPointer<vtkPolyData>::New();
+  vtkSmartPointer<vtkPoints> pdCartPoints = vtkSmartPointer<vtkPoints>::New();
+  pdCart->SetPoints(pdCartPoints);
+  MakeTestDataSet(pdCart, CoordinateType::CART);
+
+  vtkSmartPointer<vtkmCoordinateSystemTransform> cstFilter =
+      vtkSmartPointer<vtkmCoordinateSystemTransform>::New();
+  {
+    cstFilter->SetInputData(pdCart);
+    cstFilter->SetCartesianToCylindrical();
+    cstFilter->Update();
+    vtkPolyData* pdCarToCyl = vtkPolyData::SafeDownCast(cstFilter->GetOutput());
+
+    vtkSmartPointer<vtkPolyData> pdCarToCylCopy = vtkSmartPointer<vtkPolyData>::New();
+    pdCarToCylCopy->ShallowCopy(pdCarToCyl);
+    cstFilter->SetInputData(pdCarToCylCopy);
+    cstFilter->SetCylindricalToCartesian();
+    cstFilter->Update();
+    vtkPolyData* pdCylToCar = vtkPolyData::SafeDownCast(cstFilter->GetOutput());
+    ValidateCoordTransform(pdCart, pdCylToCar, {false, false, false});
+  }
+
+  // Test cylindrical to cartesian
+  vtkSmartPointer<vtkPolyData> pdCyl = vtkSmartPointer<vtkPolyData>::New();
+  vtkSmartPointer<vtkPoints> pdCylPoints = vtkSmartPointer<vtkPoints>::New();
+  pdCyl->SetPoints(pdCylPoints);
+  MakeTestDataSet(pdCyl, CoordinateType::CYL);
+
+  {
+    cstFilter->SetInputData(pdCyl);
+    cstFilter->SetCylindricalToCartesian();
+    cstFilter->Update();
+    vtkPolyData* pdCylToCal = vtkPolyData::SafeDownCast(cstFilter->GetOutput());
+
+    vtkSmartPointer<vtkPolyData> pdCylToCarCopy = vtkSmartPointer<vtkPolyData>::New();
+    pdCylToCarCopy->ShallowCopy(pdCylToCal);
+    cstFilter->SetInputData(pdCylToCarCopy);
+    cstFilter->SetCartesianToCylindrical();
+    cstFilter->Update();
+    vtkPolyData* pdCarToCyl = vtkPolyData::SafeDownCast(cstFilter->GetOutput());
+    ValidateCoordTransform(pdCyl, pdCarToCyl, {true, true, false});
+  }
+
+  // Test cartesian to spherical
+  {
+    cstFilter->SetInputData(pdCart);
+    cstFilter->SetCartesianToSpherical();
+    cstFilter->Update();
+    vtkPolyData* pdCarToSph = vtkPolyData::SafeDownCast(cstFilter->GetOutput());
+
+    vtkSmartPointer<vtkPolyData> pdCarToSphCopy = vtkSmartPointer<vtkPolyData>::New();
+    pdCarToSphCopy->ShallowCopy(pdCarToSph);
+    cstFilter->SetInputData(pdCarToSphCopy);
+    cstFilter->SetSphericalToCartesian();
+    cstFilter->Update();
+    vtkPolyData* pdSphToCar = vtkPolyData::SafeDownCast(cstFilter->GetOutput());
+    ValidateCoordTransform(pdCart, pdSphToCar, {false, false, false});
+  }
+
+  // Test spherical to cartesian
+  vtkSmartPointer<vtkPolyData> pdSph = vtkSmartPointer<vtkPolyData>::New();
+  vtkSmartPointer<vtkPoints> pdSphPoints = vtkSmartPointer<vtkPoints>::New();
+  pdSph->SetPoints(pdSphPoints);
+  MakeTestDataSet(pdSph, CoordinateType::SPH);
+  {
+    cstFilter->SetInputData(pdSph);
+    cstFilter->SetSphericalToCartesian();
+    cstFilter->Update();
+    vtkPolyData*  pdSphToCar = vtkPolyData::SafeDownCast(cstFilter->GetOutput());
+
+    vtkSmartPointer<vtkPolyData> pdSphToCarCopy = vtkSmartPointer<vtkPolyData>::New();
+    pdSphToCarCopy->ShallowCopy(pdSphToCar);
+    cstFilter->SetInputData(pdSphToCarCopy);
+    cstFilter->SetCartesianToSpherical();
+    cstFilter->Update();
+    vtkPolyData* pdCarToSph = vtkPolyData::SafeDownCast(cstFilter->GetOutput());
+    ValidateCoordTransform(pdSph, pdCarToSph, {false, true, true});
+  }
+  return 0;
+}

Accelerators/Vtkm/vtkmCoordinateSystemTransform.cu

+//=============================================================================
+//
+//  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 2012 Sandia Corporation.
+//  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
+//  the U.S. Government retains certain rights in this software.
+//
+//=============================================================================
+#include "vtkmCoordinateSystemTransform.cxx"

Accelerators/Vtkm/vtkmCoordinateSystemTransform.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 2012 Sandia Corporation.
+//  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
+//  the U.S. Government retains certain rights in this software.
+//
+//=============================================================================
+#include "vtkmCoordinateSystemTransform.h"
+#include "vtkmConfig.h"
+
+#include "vtkCellData.h"
+#include "vtkDataSet.h"
+#include "vtkImageData.h"
+#include "vtkImageDataToPointSet.h"
+#include "vtkInformation.h"
+#include "vtkInformationVector.h"
+#include "vtkObjectFactory.h"
+#include "vtkPointData.h"
+#include "vtkRectilinearGrid.h"
+#include "vtkRectilinearGridToPointSet.h"
+#include "vtkUnstructuredGrid.h"
+
+#include "vtkmlib/ArrayConverters.h"
+#include "vtkmlib/DataSetConverters.h"
+#include "vtkmlib/Storage.h"
+
+#include "vtkmFilterPolicy.h"
+
+#include <vtkm/filter/CoordinateSystemTransform.h>
+
+vtkStandardNewMacro(vtkmCoordinateSystemTransform)
+
+//------------------------------------------------------------------------------
+vtkmCoordinateSystemTransform::vtkmCoordinateSystemTransform()
+{
+  this->TransformType = TransformTypes::None;
+}
+
+//------------------------------------------------------------------------------
+vtkmCoordinateSystemTransform::~vtkmCoordinateSystemTransform()
+{
+}
+
+//------------------------------------------------------------------------------
+void vtkmCoordinateSystemTransform::SetCartesianToCylindrical()
+{
+  this->TransformType = TransformTypes::CarToCyl;
+}
+
+//------------------------------------------------------------------------------
+void vtkmCoordinateSystemTransform::SetCylindricalToCartesian()
+{
+  this->TransformType = TransformTypes::CylToCar;
+}
+
+//------------------------------------------------------------------------------
+void vtkmCoordinateSystemTransform::SetCartesianToSpherical()
+{
+  this->TransformType = TransformTypes::CarToSph;
+}
+
+//------------------------------------------------------------------------------
+void vtkmCoordinateSystemTransform::SetSphericalToCartesian()
+{
+  this->TransformType = TransformTypes::SphToCar;
+}
+
+//------------------------------------------------------------------------------
+int  vtkmCoordinateSystemTransform::FillInputPortInformation(int vtkNotUsed(port),
+                                                             vtkInformation* info)
+{
+  info->Remove(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE());
+  info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkPointSet");
+  info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkImageData");
+  info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkRectilinearGrid");
+  return 1;
+}
+
+//------------------------------------------------------------------------------
+int vtkmCoordinateSystemTransform::RequestDataObject(vtkInformation *request,
+                                          vtkInformationVector** inputVector,
+                                          vtkInformationVector* outputVector)
+{
+  vtkImageData* inImage = vtkImageData::GetData(inputVector[0]);
+  vtkRectilinearGrid* inRect = vtkRectilinearGrid::GetData(inputVector[0]);
+
+  if (inImage || inRect)
+  {
+    vtkStructuredGrid* output = vtkStructuredGrid::GetData(outputVector);
+    if (!output)
+    {
+      vtkNew<vtkStructuredGrid> newOutput;
+      outputVector->GetInformationObject(0)->Set(
+        vtkDataObject::DATA_OBJECT(), newOutput);
+    }
+    return 1;
+  }
+  else
+  {
+    return this->Superclass::RequestDataObject(request,
+                                               inputVector,
+                                               outputVector);
+  }
+}
+
+//------------------------------------------------------------------------------
+int vtkmCoordinateSystemTransform::RequestData(vtkInformation*,
+                                               vtkInformationVector** inputVector,
+                                               vtkInformationVector* outputVector)
+{
+  vtkSmartPointer<vtkPointSet> input = vtkPointSet::GetData(inputVector[0]);
+  vtkPointSet *output = vtkPointSet::GetData(outputVector);
+
+  if (!input)
+  {
+    // Try converting image data.
+    vtkImageData *inImage = vtkImageData::GetData(inputVector[0]);
+    if (inImage)
+    {
+      vtkNew<vtkImageDataToPointSet> image2points;
+      image2points->SetInputData(inImage);
+      image2points->Update();
+      input = image2points->GetOutput();
+    }
+  }
+
+  if (!input)
+  {
+    // Try converting rectilinear grid.
+    vtkRectilinearGrid *inRect = vtkRectilinearGrid::GetData(inputVector[0]);
+    if (inRect)
+    {
+      vtkNew<vtkRectilinearGridToPointSet> rect2points;
+      rect2points->SetInputData(inRect);
+      rect2points->Update();
+      input = rect2points->GetOutput();
+    }
+  }
+
+  if (!input)
+  {
+    vtkErrorMacro(<< "Invalid or missing input");
+    return 0;
+  }
+
+  output->CopyStructure(input);
+
+  vtkPoints* inPts = input->GetPoints();
+
+  if (!inPts || this->TransformType == TransformTypes::None)
+  {
+    vtkErrorMacro(<< "Miss input points or transform type has not been specified");
+    return 0;
+  }
+
+  try
+  {
+    vtkm::cont::DataSet in = tovtkm::Convert(input, tovtkm::FieldsFlag::Points);
+    vtkmInputFilterPolicy policy;
+    vtkDataArray* transformResult;
+    if (this->TransformType == TransformTypes::CarToCyl ||
+        this->TransformType == TransformTypes::CylToCar)
+    { // Cylindrical coordinate transform
+      vtkm::filter::CylindricalCoordinateTransform cylindricalCT;
+      cylindricalCT.SetUseCoordinateSystemAsField(true);
+      (this->TransformType == TransformTypes::CarToCyl) ?
+                                cylindricalCT.SetCartesianToCylindrical() :
+                                cylindricalCT.SetCylindricalToCartesian();
+      auto result = cylindricalCT.Execute(in, policy);
+      transformResult = fromvtkm::Convert(result.GetField(
+                                 "cylindricalCoordinateSystemTransform",
+                                 vtkm::cont::Field::Association::POINTS));
+    }
+    else
+    { // Spherical coordinate system
+      vtkm::filter::SphericalCoordinateTransform sphericalCT;
+      sphericalCT.SetUseCoordinateSystemAsField(true);
+      (this->TransformType == TransformTypes::CarToSph) ?
+                                sphericalCT.SetCartesianToSpherical() :
+                                sphericalCT.SetSphericalToCartesian();
+      auto result = sphericalCT.Execute(in, policy);
+      transformResult = fromvtkm::Convert(result.GetField(
+                                 "sphericalCoordinateSystemTransform",
+                                 vtkm::cont::Field::Association::POINTS));
+    }
+    vtkPoints* newPts = vtkPoints::New();
+    // Update points
+    newPts->SetNumberOfPoints(transformResult->GetNumberOfTuples());
+    newPts->SetData(transformResult);
+    output->SetPoints(newPts);
+    newPts->Delete();
+    transformResult->FastDelete();
+  }
+  catch(const vtkm::cont::Error& e)
+  {
+    vtkErrorMacro(<< "VTK-m error: " << e.GetMessage());
+  }
+
+  // Update ourselves and release memory
+  output->GetPointData()->CopyNormalsOff(); // distorted geometry
+  output->GetPointData()->PassData(input->GetPointData());
+  output->GetCellData()->CopyNormalsOff(); // distorted geometry
+  output->GetCellData()->PassData(input->GetCellData());
+
+  return 1;
+}
+
+//------------------------------------------------------------------------------
+void vtkmCoordinateSystemTransform::PrintSelf(std::ostream &os, vtkIndent indent)
+{
+  this->Superclass::PrintSelf(os, indent);
+}

Accelerators/Vtkm/vtkmCoordinateSystemTransform.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 2012 Sandia Corporation.
+//  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
+//  the U.S. Government retains certain rights in this software.
+//
+//=============================================================================
+/**
+ * @class   vtkmCoordinateSystemTransform
+ * @brief   transform a coordinate system between Cartesian&Cylindrical and
+ *          Cartesian&Spherical
+ *
+ * vtkmCoordinateSystemTransform is a filter that transforms a coordinate system
+ * between Cartesian&Cylindrical and Cartesian&Spherical.
+*/
+
+#ifndef vtkmCoordinateSystemTransform_h
+#define vtkmCoordinateSystemTransform_h
+
+#include "vtkPointSetAlgorithm.h"
+#include "vtkAcceleratorsVTKmModule.h" // required for correct export
+
+class VTKACCELERATORSVTKM_EXPORT vtkmCoordinateSystemTransform : public vtkPointSetAlgorithm
+{
+  enum struct TransformTypes {None, CarToCyl, CylToCar, CarToSph, SphToCar};
+public:
+  vtkTypeMacro(vtkmCoordinateSystemTransform, vtkPointSetAlgorithm)
+  void PrintSelf(ostream& os, vtkIndent indent) override;
+
+  static vtkmCoordinateSystemTransform* New();
+
+  void SetCartesianToCylindrical();
+  void SetCylindricalToCartesian();
+
+  void SetCartesianToSpherical();
+  void SetSphericalToCartesian();
+
+  int FillInputPortInformation(int port, vtkInformation* info) override;
+protected:
+  vtkmCoordinateSystemTransform();
+  ~vtkmCoordinateSystemTransform();
+
+  int RequestDataObject(vtkInformation* request,
+                        vtkInformationVector** inputVector,
+                        vtkInformationVector* outputVector) override;
+  int RequestData(vtkInformation* , vtkInformationVector**,
+                          vtkInformationVector*) override;
+
+private:
+  vtkmCoordinateSystemTransform(const vtkmCoordinateSystemTransform&) = delete;
+  void operator=(const vtkmCoordinateSystemTransform&) = delete;
+
+  TransformTypes TransformType;
+};
+
+#endif // vtkmCoordinateSystemTransform_h
+
+// VTK-HeaderTest-Exclude: vtkmCoordinateSystemTransform.h