Merge topic 'image_median_filters'

d47f7aaa2 Add an ImageMedian filter

Acked-by: Kitware Robot <kwrobot@kitware.com>
Merge-request: !1890

vtkm/filter/CMakeLists.txt

@@ -42,6 +42,7 @@ set(headers
   Gradient.h
   Histogram.h
   ImageConnectivity.h
+  ImageMedian.h
   Lagrangian.h
   LagrangianStructures.h
   Mask.h
@@ -107,6 +108,7 @@ set(header_template_sources
   Gradient.hxx
   Histogram.hxx
   ImageConnectivity.hxx
+  ImageMedian.hxx
   Lagrangian.hxx
   LagrangianStructures.hxx
   Mask.hxx

vtkm/filter/ImageMedian.h

@@ -0,0 +1,54 @@
+//============================================================================
+//  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.
+//============================================================================
+
+#ifndef vtk_m_filter_ImageMedian_h
+#define vtk_m_filter_ImageMedian_h
+
+#include <vtkm/filter/FilterCell.h>
+
+/// \brief Median algorithm for general image blur
+///
+/// The ImageMedian filter finds the median value for each pixel in an image.
+/// Currently the algorithm has the following restrictions.
+///   - Only supports a neighborhood of 5x5x1 or 3x3x1
+///
+/// This means that volumes are basically treated as an image stack
+/// along the z axis
+///
+/// Default output field name is 'median'
+namespace vtkm
+{
+namespace filter
+{
+class ImageMedian : public vtkm::filter::FilterField<ImageMedian>
+{
+  int Neighborhood = 1;
+
+public:
+  using SupportedTypes = vtkm::TypeListTagScalarAll;
+
+  VTKM_CONT ImageMedian();
+
+  VTKM_CONT void Perform3x3() { this->Neighborhood = 1; };
+  VTKM_CONT void Perform5x5() { this->Neighborhood = 2; };
+
+
+  template <typename T, typename StorageType, typename DerivedPolicy>
+  VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input,
+                                          const vtkm::cont::ArrayHandle<T, StorageType>& field,
+                                          const vtkm::filter::FieldMetadata& fieldMetadata,
+                                          const vtkm::filter::PolicyBase<DerivedPolicy>&);
+};
+}
+} // namespace vtkm::filter
+
+#include <vtkm/filter/ImageMedian.hxx>
+
+#endif //vtk_m_filter_ImageMedian_h

vtkm/filter/ImageMedian.hxx

@@ -0,0 +1,131 @@
+//============================================================================
+//  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.
+//============================================================================
+
+#ifndef vtk_m_filter_ImageMedian_hxx
+#define vtk_m_filter_ImageMedian_hxx
+
+#include <vtkm/Swap.h>
+
+namespace vtkm
+{
+namespace worklet
+{
+
+// An implementation of the quickselect/Hoare's selection algorithm to find medians
+// inplace, generally fairly fast for reasonable sized data.
+//
+template <typename T>
+VTKM_EXEC inline T find_median(T* values, std::size_t mid, std::size_t size)
+{
+  std::size_t begin = 0;
+  std::size_t end = size - 1;
+  while (begin < end)
+  {
+    T x = values[mid];
+    std::size_t i = begin;
+    std::size_t j = end;
+    do
+    {
+      for (; values[i] < x; i++)
+      {
+      }
+      for (; x < values[j]; j--)
+      {
+      }
+      if (i <= j)
+      {
+        vtkm::Swap(values[i], values[j]);
+        i++;
+        j--;
+      }
+    } while (i <= j);
+
+    begin = (j < mid) ? i : begin;
+    end = (mid < i) ? j : end;
+  }
+  return values[mid];
+}
+struct ImageMedian : public vtkm::worklet::WorkletPointNeighborhood
+{
+  int Neighborhood;
+  ImageMedian(int neighborhoodSize)
+    : Neighborhood(neighborhoodSize)
+  {
+  }
+  using ControlSignature = void(CellSetIn, FieldInNeighborhood, FieldOut);
+  using ExecutionSignature = void(_2, _3);
+
+  template <typename InNeighborhoodT, typename T>
+  VTKM_EXEC void operator()(const InNeighborhoodT& input, T& out) const
+  {
+    vtkm::Vec<T, 25> values;
+    int index = 0;
+    for (int x = -this->Neighborhood; x <= this->Neighborhood; ++x)
+    {
+      for (int y = -this->Neighborhood; y <= this->Neighborhood; ++y)
+      {
+        values[index++] = input.Get(x, y, 0);
+      }
+    }
+
+    std::size_t len =
+      static_cast<std::size_t>((this->Neighborhood * 2 + 1) * (this->Neighborhood * 2 + 1));
+    std::size_t mid = len / 2;
+    out = find_median(&values[0], mid, len);
+  }
+};
+}
+
+namespace filter
+{
+
+VTKM_CONT ImageMedian::ImageMedian()
+{
+  this->SetOutputFieldName("median");
+}
+
+template <typename T, typename StorageType, typename DerivedPolicy>
+inline VTKM_CONT vtkm::cont::DataSet ImageMedian::DoExecute(
+  const vtkm::cont::DataSet& input,
+  const vtkm::cont::ArrayHandle<T, StorageType>& field,
+  const vtkm::filter::FieldMetadata& fieldMetadata,
+  const vtkm::filter::PolicyBase<DerivedPolicy>& policy)
+{
+  if (!fieldMetadata.IsPointField())
+  {
+    throw vtkm::cont::ErrorBadValue("Active field for ImageMedian must be a point field.");
+  }
+
+  const vtkm::cont::DynamicCellSet& cells = input.GetCellSet();
+  vtkm::cont::ArrayHandle<T> result;
+  if (this->Neighborhood == 1 || this->Neighborhood == 2)
+  {
+    this->Invoke(worklet::ImageMedian{ this->Neighborhood },
+                 vtkm::filter::ApplyPolicyCellSetStructured(cells, policy),
+                 field,
+                 result);
+  }
+  else
+  {
+    throw vtkm::cont::ErrorBadValue("ImageMedian only support a 3x3 or 5x5 stencil.");
+  }
+
+  std::string name = this->GetOutputFieldName();
+  if (name.empty())
+  {
+    name = fieldMetadata.GetName();
+  }
+
+  return CreateResult(input, fieldMetadata.AsField(name, result));
+}
+}
+}
+
+#endif

vtkm/filter/testing/CMakeLists.txt

@@ -34,6 +34,7 @@ set(unit_tests
   UnitTestGhostCellRemove.cxx
   UnitTestHistogramFilter.cxx
   UnitTestImageConnectivityFilter.cxx
+  UnitTestImageMedianFilter.cxx
   UnitTestLagrangianFilter.cxx
   UnitTestLagrangianStructuresFilter.cxx
   UnitTestMaskFilter.cxx

vtkm/filter/testing/UnitTestImageMedianFilter.cxx

@@ -0,0 +1,56 @@
+//============================================================================
+//  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.
+//============================================================================
+
+#include <vtkm/cont/testing/MakeTestDataSet.h>
+#include <vtkm/cont/testing/Testing.h>
+#include <vtkm/filter/ImageMedian.h>
+
+namespace
+{
+
+void TestImageMedian()
+{
+  std::cout << "Testing Image Median Filter on 3D structured data" << std::endl;
+
+  vtkm::cont::testing::MakeTestDataSet testDataSet;
+  vtkm::cont::DataSet dataSet = testDataSet.Make3DUniformDataSet2();
+
+  vtkm::filter::ImageMedian median;
+  median.Perform3x3();
+  median.SetActiveField("pointvar");
+  auto result = median.Execute(dataSet);
+
+  VTKM_TEST_ASSERT(result.HasPointField("median"), "Field missing.");
+  vtkm::cont::ArrayHandle<vtkm::Float32> resultArrayHandle;
+  result.GetPointField("median").GetData().CopyTo(resultArrayHandle);
+
+  auto cells = result.GetCellSet().Cast<vtkm::cont::CellSetStructured<3>>();
+  auto pdims = cells.GetPointDimensions();
+
+  //verified by hand
+  {
+    auto portal = resultArrayHandle.GetPortalConstControl();
+    std::cout << "spot to verify x = 1, y = 1, z = 0 is: ";
+    vtkm::Float32 temp = portal.Get(1 + pdims[0]);
+    std::cout << temp << std::endl << std::endl;
+    VTKM_TEST_ASSERT(test_equal(temp, 2), "incorrect median value");
+
+    std::cout << "spot to verify x = 1, y = 1, z = 2 is: ";
+    temp = portal.Get(1 + pdims[0] + (pdims[1] * pdims[0] * 2));
+    std::cout << temp << std::endl << std::endl;
+    VTKM_TEST_ASSERT(test_equal(temp, 2.82843), "incorrect median value");
+  }
+}
+}
+
+int UnitTestImageMedianFilter(int argc, char* argv[])
+{
+  return vtkm::cont::testing::Testing::Run(TestImageMedian, argc, argv);
+}