migrate image processing filters

This commit is contained in:
Li-Ta Lo 2022-02-07 10:03:08 -07:00
parent f3c65d3404
commit c161ea03ce
24 changed files with 576 additions and 425 deletions

@ -14,6 +14,7 @@ set(deprecated_headers
CleanGrid.h
ClipWithField.h
ClipWithImplicitFunction.h
ComputeMoments.h
Contour.h
CoordinateSystemTransform.h
CrossProduct.h
@ -29,6 +30,8 @@ set(deprecated_headers
Gradient.h
Histogram.h
ImageConnectivity.h
ImageDifference.h
ImageMedian.h
Mask.h
MaskPoints.h
NDEntropy.h
@ -82,15 +85,12 @@ vtkm_declare_headers(${common_header_template_sources})
set(extra_headers
AmrArrays.h
ComputeMoments.h
ContourTreeUniformAugmented.h
ContourTreeUniformDistributed.h
ContourTreeUniform.h
CreateResult.h
FieldSelection.h
GhostCellClassify.h
ImageDifference.h
ImageMedian.h
Lagrangian.h
LagrangianStructures.h
MeshQuality.h
@ -116,13 +116,10 @@ set(extra_headers
set(extra_header_template_sources
AmrArrays.hxx
ComputeMoments.hxx
ContourTreeUniformAugmented.hxx
ContourTreeUniformDistributed.hxx
ContourTreeUniform.hxx
GhostCellClassify.hxx
ImageDifference.hxx
ImageMedian.hxx
Lagrangian.hxx
LagrangianStructures.hxx
MeshQuality.hxx
@ -210,6 +207,7 @@ add_subdirectory(connected_components)
add_subdirectory(contour)
add_subdirectory(density_estimate)
add_subdirectory(entity_extraction)
add_subdirectory(image_processing)
add_subdirectory(internal)
add_subdirectory(particleadvection)
add_subdirectory(field_conversion)

@ -6,55 +6,35 @@
// 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_ComputeMoments_h
#define vtk_m_filter_ComputeMoments_h
#include <vtkm/filter/FilterField.h>
#include <vtkm/Deprecated.h>
#include <vtkm/filter/image_processing/ComputeMoments.h>
namespace vtkm
{
namespace filter
{
class ComputeMoments : public vtkm::filter::FilterField<ComputeMoments>
VTKM_DEPRECATED(
1.8,
"Use vtkm/filter/image_processing/ComputeMoments.h instead of vtkm/filter/ComputeMoments.h.")
inline void ComputeMoments_deprecated() {}
inline void ComputeMoments_deprecated_warning()
{
public:
using SupportedTypes = vtkm::List<vtkm::Float32,
vtkm::Float64,
vtkm::Vec<vtkm::Float32, 2>,
vtkm::Vec<vtkm::Float64, 2>,
vtkm::Vec<vtkm::Float32, 3>,
vtkm::Vec<vtkm::Float64, 3>,
vtkm::Vec<vtkm::Float32, 4>,
vtkm::Vec<vtkm::Float64, 4>,
vtkm::Vec<vtkm::Float32, 6>,
vtkm::Vec<vtkm::Float64, 6>,
vtkm::Vec<vtkm::Float32, 9>,
vtkm::Vec<vtkm::Float64, 9>>;
ComputeMoments_deprecated();
}
VTKM_CONT ComputeMoments();
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>&);
VTKM_CONT void SetRadius(double _radius) { this->Radius = _radius; }
VTKM_CONT void SetSpacing(vtkm::Vec3f _spacing) { this->Spacing = _spacing; }
VTKM_CONT void SetOrder(vtkm::Int32 _order) { this->Order = _order; }
private:
double Radius = 1;
vtkm::Vec3f Spacing = { 1.0f, 1.0f, 1.0f };
vtkm::Int32 Order = 0;
class VTKM_DEPRECATED(1.8, "Use vtkm::filter::image_processing::ComputeMoments.") ComputeMoments
: public vtkm::filter::image_processing::ComputeMoments
{
using image_processing::ComputeMoments::ComputeMoments;
};
}
} // namespace vtkm::filter
#include <vtkm/filter/ComputeMoments.hxx>
#endif //vtk_m_filter_ComputeMoments_h

@ -1,50 +0,0 @@
//============================================================================
// 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_ComputeMoments_hxx
#define vtk_m_filter_ComputeMoments_hxx
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/worklet/moments/ComputeMoments.h>
namespace vtkm
{
namespace filter
{
VTKM_CONT ComputeMoments::ComputeMoments()
{
this->SetOutputFieldName("moments_");
}
template <typename T, typename StorageType, typename DerivedPolicy>
inline VTKM_CONT vtkm::cont::DataSet ComputeMoments::DoExecute(
const vtkm::cont::DataSet& input,
const vtkm::cont::ArrayHandle<T, StorageType>& field,
const vtkm::filter::FieldMetadata& fieldMetadata,
const vtkm::filter::PolicyBase<DerivedPolicy>&)
{
if (fieldMetadata.GetAssociation() != vtkm::cont::Field::Association::POINTS)
{
throw vtkm::cont::ErrorBadValue("Active field for ComputeMoments must be a point field.");
}
vtkm::cont::DataSet output;
output.CopyStructure(input);
auto worklet = vtkm::worklet::moments::ComputeMoments(this->Spacing, this->Radius);
worklet.Run(input.GetCellSet(), field, this->Order, output);
return output;
}
}
}
#endif //vtk_m_filter_ComputeMoments_hxx

@ -7,118 +7,34 @@
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_filter_ImageDifference_h
#define vtk_m_filter_ImageDifference_h
#include <vtkm/TypeList.h>
#include <vtkm/filter/FilterField.h>
#include <vtkm/Deprecated.h>
#include <vtkm/filter/image_processing/ImageDifference.h>
namespace vtkm
{
namespace filter
{
/// \brief Construct an ImageDifference of a given DataSet
///
/// The dataset generated by executing this filter is a Dataset with two Fields:
/// - "image-diff": Uniform Structured Dataset, difference values of image A - B
/// - "threshold-output": Uniform Structured Dataset, the magnitudes of the pixel differences
///
/// The threshold-output is calculated for each pixel using the `vtkm::Magnitude` vector function
/// on the individual pixel difference.
///
class ImageDifference : public vtkm::filter::FilterField<ImageDifference>
VTKM_DEPRECATED(
1.8,
"Use vtkm/filter/image_processing/ImageDifference.h instead of vtkm/filter/ImageDifference.h.")
inline void ImageDifference_deprecated() {}
inline void ImageDifference_deprecated_warning()
{
public:
using SupportedTypes = vtkm::TypeListFieldVec4;
ImageDifference_deprecated();
}
VTKM_CONT ImageDifference();
VTKM_CONT vtkm::IdComponent GetAverageRadius() const { return this->AverageRadius; }
VTKM_CONT void SetAverageRadius(const vtkm::IdComponent& averageRadius)
{
this->AverageRadius = averageRadius;
}
VTKM_CONT vtkm::IdComponent GetPixelShiftRadius() const { return this->PixelShiftRadius; }
VTKM_CONT void SetPixelShiftRadius(const vtkm::IdComponent& pixelShiftRadius)
{
this->PixelShiftRadius = pixelShiftRadius;
}
VTKM_CONT vtkm::FloatDefault GetAllowedPixelErrorRatio() const
{
return this->AllowedPixelErrorRatio;
}
VTKM_CONT void SetAllowedPixelErrorRatio(const vtkm::FloatDefault& pixelErrorRatio)
{
this->AllowedPixelErrorRatio = pixelErrorRatio;
}
VTKM_CONT vtkm::FloatDefault GetPixelDiffThreshold() const { return this->PixelDiffThreshold; }
VTKM_CONT void SetPixelDiffThreshold(const vtkm::FloatDefault& threshold)
{
this->PixelDiffThreshold = threshold;
}
VTKM_CONT bool GetImageDiffWithinThreshold() const { return this->ImageDiffWithinThreshold; }
VTKM_CONT void SetThresholdFieldName(const std::string& name) { this->ThresholdFieldName = name; }
VTKM_CONT std::string GetThresholdFieldName() const { return this->ThresholdFieldName; }
/// Choose the primary field to operate on. For Image difference A - B, A is the
/// primary field.
VTKM_CONT
void SetPrimaryField(
const std::string& name,
vtkm::cont::Field::Association association = vtkm::cont::Field::Association::ANY)
{
this->SetActiveField(name, association);
}
VTKM_CONT std::string GetPrimaryFieldName() const { return this->GetActiveFieldName(); }
VTKM_CONT vtkm::cont::Field::Association GetPrimaryFieldAssociation() const
{
return this->GetActiveFieldAssociation();
}
/// Choose the secondary field to operate on. For Image difference A - B, B is the
/// secondary field.
VTKM_CONT
void SetSecondaryField(
const std::string& name,
vtkm::cont::Field::Association association = vtkm::cont::Field::Association::ANY)
{
this->SecondaryFieldName = name;
this->SecondaryFieldAssociation = association;
}
VTKM_CONT std::string GetSecondaryFieldName() const { return this->SecondaryFieldName; }
VTKM_CONT vtkm::cont::Field::Association GetSecondaryFieldAssociation() const
{
return this->SecondaryFieldAssociation;
}
template <typename T, typename StorageType, typename DerivedPolicy>
VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input,
const vtkm::cont::ArrayHandle<T, StorageType>& primary,
const vtkm::filter::FieldMetadata& fieldMetadata,
vtkm::filter::PolicyBase<DerivedPolicy> policy);
private:
vtkm::IdComponent AverageRadius;
vtkm::IdComponent PixelShiftRadius;
vtkm::FloatDefault AllowedPixelErrorRatio;
vtkm::FloatDefault PixelDiffThreshold;
bool ImageDiffWithinThreshold;
std::string SecondaryFieldName;
vtkm::cont::Field::Association SecondaryFieldAssociation;
std::string ThresholdFieldName;
class VTKM_DEPRECATED(1.8, "Use vtkm::filter::image_processing::ImageDifference.") ImageDifference
: public vtkm::filter::image_processing::ImageDifference
{
using image_processing::ImageDifference::ImageDifference;
};
} // namespace filter
} // namespace vtkm
}
} // namespace vtkm::filter
#include <vtkm/filter/ImageDifference.hxx>
#endif // vtk_m_filter_ImageDifference_h
#endif //vtk_m_filter_ImageDifference_h

@ -1,145 +0,0 @@
//============================================================================
// 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_ImageDifference_hxx
#define vtk_m_filter_ImageDifference_hxx
#include <vtkm/filter/ImageDifference.h>
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/ArrayPortalToIterators.h>
#include <vtkm/cont/Logging.h>
#include <vtkm/worklet/AveragePointNeighborhood.h>
#include <vtkm/worklet/ImageDifference.h>
namespace vtkm
{
namespace filter
{
namespace detail
{
struct GreaterThanThreshold
{
GreaterThanThreshold(const vtkm::FloatDefault& thresholdError)
: ThresholdError(thresholdError)
{
}
VTKM_EXEC_CONT bool operator()(const vtkm::FloatDefault& x) const { return x > ThresholdError; }
vtkm::FloatDefault ThresholdError;
};
} // namespace detail
inline VTKM_CONT ImageDifference::ImageDifference()
: vtkm::filter::FilterField<ImageDifference>()
, AverageRadius(0)
, PixelShiftRadius(0)
, AllowedPixelErrorRatio(0.00025f)
, PixelDiffThreshold(0.05f)
, ImageDiffWithinThreshold(true)
, SecondaryFieldName("image-2")
, SecondaryFieldAssociation(vtkm::cont::Field::Association::ANY)
, ThresholdFieldName("threshold-output")
{
this->SetPrimaryField("image-1");
this->SetOutputFieldName("image-diff");
}
template <typename T, typename StorageType, typename DerivedPolicy>
inline VTKM_CONT vtkm::cont::DataSet ImageDifference::DoExecute(
const vtkm::cont::DataSet& input,
const vtkm::cont::ArrayHandle<T, StorageType>& primary,
const vtkm::filter::FieldMetadata& fieldMetadata,
vtkm::filter::PolicyBase<DerivedPolicy> policy)
{
this->ImageDiffWithinThreshold = true;
if (!fieldMetadata.IsPointField())
{
throw vtkm::cont::ErrorFilterExecution("Point field expected.");
}
VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Performing Image Difference");
vtkm::cont::Field secondaryField;
secondaryField = input.GetField(this->SecondaryFieldName, this->SecondaryFieldAssociation);
auto secondary = vtkm::filter::ApplyPolicyFieldOfType<T>(secondaryField, policy, *this);
auto cellSet = vtkm::filter::ApplyPolicyCellSetStructured(input.GetCellSet(), policy, *this);
vtkm::cont::ArrayHandle<T, StorageType> diffOutput;
vtkm::cont::ArrayHandle<T, StorageType> primaryOutput;
vtkm::cont::ArrayHandle<T, StorageType> secondaryOutput;
vtkm::cont::ArrayHandle<vtkm::FloatDefault> thresholdOutput;
if (this->AverageRadius > 0)
{
VTKM_LOG_S(vtkm::cont::LogLevel::Info,
"Performing Average with radius: " << this->AverageRadius);
auto averageWorklet = vtkm::worklet::AveragePointNeighborhood(this->AverageRadius);
this->Invoke(averageWorklet, cellSet, primary, primaryOutput);
this->Invoke(averageWorklet, cellSet, secondary, secondaryOutput);
}
else
{
VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Not performing average");
primaryOutput = primary;
vtkm::cont::ArrayCopyShallowIfPossible(secondaryField.GetData(), secondaryOutput);
}
if (this->PixelShiftRadius > 0)
{
VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Diffing image in Neighborhood");
auto diffWorklet =
vtkm::worklet::ImageDifferenceNeighborhood(this->PixelShiftRadius, this->PixelDiffThreshold);
this->Invoke(diffWorklet, cellSet, primaryOutput, secondaryOutput, diffOutput, thresholdOutput);
}
else
{
VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Diffing image directly");
auto diffWorklet = vtkm::worklet::ImageDifference();
this->Invoke(diffWorklet, primaryOutput, secondaryOutput, diffOutput, thresholdOutput);
}
vtkm::cont::ArrayHandle<vtkm::FloatDefault> errorPixels;
vtkm::cont::Algorithm::CopyIf(thresholdOutput,
thresholdOutput,
errorPixels,
detail::GreaterThanThreshold(this->PixelDiffThreshold));
if (errorPixels.GetNumberOfValues() >
thresholdOutput.GetNumberOfValues() * this->AllowedPixelErrorRatio)
{
this->ImageDiffWithinThreshold = false;
}
VTKM_LOG_S(vtkm::cont::LogLevel::Info,
"Difference within threshold: "
<< this->ImageDiffWithinThreshold
<< ", for pixels outside threshold: " << errorPixels.GetNumberOfValues()
<< ", with a total number of pixesl: " << thresholdOutput.GetNumberOfValues()
<< ", and an allowable pixel error ratio: " << this->AllowedPixelErrorRatio
<< ", with a total summed threshold error: "
<< vtkm::cont::Algorithm::Reduce(errorPixels, static_cast<FloatDefault>(0)));
vtkm::cont::DataSet clone;
clone.CopyStructure(input);
clone.AddField(fieldMetadata.AsField(this->GetOutputFieldName(), diffOutput));
clone.AddField(fieldMetadata.AsField(this->GetThresholdFieldName(), thresholdOutput));
VTKM_ASSERT(clone.HasField(this->GetOutputFieldName(), fieldMetadata.GetAssociation()));
VTKM_ASSERT(clone.HasField(this->GetThresholdFieldName(), fieldMetadata.GetAssociation()));
return clone;
}
} // namespace filter
} // namespace vtkm
#endif // vtk_m_filter_ImageDifference_hxx

@ -7,48 +7,34 @@
// 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/FilterField.h>
#include <vtkm/Deprecated.h>
#include <vtkm/filter/image_processing/ImageMedian.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>
VTKM_DEPRECATED(
1.8,
"Use vtkm/filter/image_processing/ImageMedian.h instead of vtkm/filter/ImageMedian.h.")
inline void ImageMedian_deprecated() {}
inline void ImageMedian_deprecated_warning()
{
int Neighborhood = 1;
ImageMedian_deprecated();
}
public:
using SupportedTypes = vtkm::TypeListScalarAll;
VTKM_CONT ImageMedian() { this->SetOutputFieldName("median"); }
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>&);
class VTKM_DEPRECATED(1.8, "Use vtkm::filter::image_processing::ImageMedian.") ImageMedian
: public vtkm::filter::image_processing::ImageMedian
{
using image_processing::ImageMedian::ImageMedian;
};
}
} // namespace vtkm::filter
#include <vtkm/filter/ImageMedian.hxx>
#endif //vtk_m_filter_ImageMedian_h

@ -36,6 +36,6 @@ endif()
vtkm_unit_tests(
SOURCES ${unit_tests}
LIBRARIES ${libraries}
ALL_BACKENDS # ArrayCopy/Sort called, needs device compiler
ALL_BACKENDS # Algorithm::Sort called, needs device compiler
USE_VTKM_JOB_POOL
)

@ -0,0 +1,37 @@
##============================================================================
## 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.
##============================================================================
set(image_processing_headers
ComputeMoments.h
ImageDifference.h
ImageMedian.h
)
set(image_processing_sources
ComputeMoments.cxx
ImageDifference.cxx
ImageMedian.cxx
)
vtkm_library(
NAME vtkm_filter_image_processing
HEADERS ${image_processing_headers}
DEVICE_SOURCES ${image_processing_sources}
USE_VTKM_JOB_POOL
)
target_link_libraries(vtkm_filter_image_processing PUBLIC vtkm_worklet vtkm_filter_core)
target_link_libraries(vtkm_filter PUBLIC INTERFACE vtkm_filter_image_processing)
add_subdirectory(worklet)
#-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -
if (VTKm_ENABLE_TESTING)
add_subdirectory(testing)
endif ()

@ -0,0 +1,64 @@
//============================================================================
// 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_ComputeMoments_hxx
#define vtk_m_filter_ComputeMoments_hxx
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/filter/image_processing/ComputeMoments.h>
#include <vtkm/filter/image_processing/worklet/ComputeMoments.h>
namespace vtkm
{
namespace filter
{
namespace image_processing
{
using SupportedTypes = vtkm::List<vtkm::Float32,
vtkm::Float64,
vtkm::Vec<vtkm::Float32, 2>,
vtkm::Vec<vtkm::Float64, 2>,
vtkm::Vec<vtkm::Float32, 3>,
vtkm::Vec<vtkm::Float64, 3>,
vtkm::Vec<vtkm::Float32, 4>,
vtkm::Vec<vtkm::Float64, 4>,
vtkm::Vec<vtkm::Float32, 6>,
vtkm::Vec<vtkm::Float64, 6>,
vtkm::Vec<vtkm::Float32, 9>,
vtkm::Vec<vtkm::Float64, 9>>;
VTKM_CONT ComputeMoments::ComputeMoments()
{
this->SetOutputFieldName("moments_");
}
VTKM_CONT vtkm::cont::DataSet ComputeMoments::DoExecute(const vtkm::cont::DataSet& input)
{
const auto& field = this->GetFieldFromDataSet(input);
if (!field.IsFieldPoint())
{
throw vtkm::cont::ErrorBadValue("Active field for ComputeMoments must be a point field.");
}
vtkm::cont::DataSet output = this->CreateResult(input);
auto worklet = vtkm::worklet::moments::ComputeMoments(this->Radius, this->Spacing);
auto resolveType = [&](const auto& concrete) {
worklet.Run(input.GetCellSet(), concrete, this->Order, output);
};
field.GetData().CastAndCallForTypesWithFloatFallback<SupportedTypes, VTKM_DEFAULT_STORAGE_LIST>(
resolveType);
return output;
}
} // namespace image_processing
} // namespace filter
} // namespace vtkm
#endif //vtk_m_filter_ComputeMoments_hxx

@ -0,0 +1,44 @@
//============================================================================
// 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_image_processing_ComputeMoments_h
#define vtk_m_filter_image_processing_ComputeMoments_h
#include <vtkm/filter/NewFilterField.h>
#include <vtkm/filter/image_processing/vtkm_filter_image_processing_export.h>
namespace vtkm
{
namespace filter
{
namespace image_processing
{
class VTKM_FILTER_IMAGE_PROCESSING_EXPORT ComputeMoments : public vtkm::filter::NewFilterField
{
public:
VTKM_CONT ComputeMoments();
VTKM_CONT void SetRadius(double _radius) { this->Radius = _radius; }
VTKM_CONT void SetSpacing(vtkm::Vec3f _spacing) { this->Spacing = _spacing; }
VTKM_CONT void SetOrder(vtkm::Int32 _order) { this->Order = _order; }
private:
VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input) override;
double Radius = 1;
vtkm::Vec3f Spacing = { 1.0f, 1.0f, 1.0f };
vtkm::Int32 Order = 0;
};
} // namespace image_processing
} // namespace filter
} // namespace vtkm
#endif //vtk_m_filter_image_processing_ComputeMoments_h

@ -0,0 +1,138 @@
//============================================================================
// 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/Algorithm.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/cont/Logging.h>
#include <vtkm/cont/UncertainCellSet.h>
#include <vtkm/filter/image_processing/ImageDifference.h>
#include <vtkm/filter/image_processing/worklet/ImageDifference.h>
#include <vtkm/worklet/AveragePointNeighborhood.h>
namespace vtkm
{
namespace filter
{
namespace image_processing
{
namespace
{
struct GreaterThanThreshold
{
GreaterThanThreshold(const vtkm::FloatDefault& thresholdError)
: ThresholdError(thresholdError)
{
}
VTKM_EXEC_CONT bool operator()(const vtkm::FloatDefault& x) const { return x > ThresholdError; }
vtkm::FloatDefault ThresholdError;
};
} // anonymous namespace
VTKM_CONT ImageDifference::ImageDifference()
{
this->SetPrimaryField("image-1");
this->SetSecondaryField("image-2");
this->SetOutputFieldName("image-diff");
}
VTKM_CONT vtkm::cont::DataSet ImageDifference::DoExecute(const vtkm::cont::DataSet& input)
{
this->ImageDiffWithinThreshold = true;
const auto& primaryField = this->GetFieldFromDataSet(input);
if (!primaryField.IsFieldPoint())
{
throw vtkm::cont::ErrorFilterExecution("Point field expected.");
}
VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Performing Image Difference");
auto inputCellSet = input.GetCellSet().ResetCellSetList<VTKM_DEFAULT_CELL_SET_LIST_STRUCTURED>();
const auto& secondaryField = this->GetFieldFromDataSet(1, input);
vtkm::cont::UnknownArrayHandle diffOutput;
vtkm::cont::ArrayHandle<vtkm::FloatDefault> thresholdOutput;
auto resolveType = [&](const auto& primaryArray) {
// use std::decay to remove const ref from the decltype of primaryArray.
using T = typename std::decay_t<decltype(primaryArray)>::ValueType;
vtkm::cont::ArrayHandle<T> secondaryArray;
vtkm::cont::ArrayCopyShallowIfPossible(secondaryField.GetData(), secondaryArray);
vtkm::cont::ArrayHandle<T> primaryOutput;
vtkm::cont::ArrayHandle<T> secondaryOutput;
if (this->AverageRadius > 0)
{
VTKM_LOG_S(vtkm::cont::LogLevel::Info,
"Performing Average with radius: " << this->AverageRadius);
auto averageWorklet = vtkm::worklet::AveragePointNeighborhood(this->AverageRadius);
this->Invoke(averageWorklet, inputCellSet, primaryArray, primaryOutput);
this->Invoke(averageWorklet, inputCellSet, secondaryArray, secondaryOutput);
}
else
{
VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Not performing average");
vtkm::cont::ArrayCopyShallowIfPossible(primaryArray, primaryOutput);
secondaryOutput = secondaryArray;
}
vtkm::cont::ArrayHandle<T> diffArray;
if (this->PixelShiftRadius > 0)
{
VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Diffing image in Neighborhood");
this->Invoke(vtkm::worklet::ImageDifferenceNeighborhood(this->PixelShiftRadius,
this->PixelDiffThreshold),
inputCellSet,
primaryOutput,
secondaryOutput,
diffArray,
thresholdOutput);
}
else
{
VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Diffing image directly");
this->Invoke(vtkm::worklet::ImageDifference(),
primaryOutput,
secondaryOutput,
diffArray,
thresholdOutput);
}
diffOutput = diffArray;
};
this->CastAndCallVecField<4>(primaryField, resolveType);
vtkm::cont::ArrayHandle<vtkm::FloatDefault> errorPixels;
vtkm::cont::Algorithm::CopyIf(
thresholdOutput, thresholdOutput, errorPixels, GreaterThanThreshold(this->PixelDiffThreshold));
if (errorPixels.GetNumberOfValues() >
thresholdOutput.GetNumberOfValues() * this->AllowedPixelErrorRatio)
{
this->ImageDiffWithinThreshold = false;
}
VTKM_LOG_S(vtkm::cont::LogLevel::Info,
"Difference within threshold: "
<< this->ImageDiffWithinThreshold
<< ", for pixels outside threshold: " << errorPixels.GetNumberOfValues()
<< ", with a total number of pixesl: " << thresholdOutput.GetNumberOfValues()
<< ", and an allowable pixel error ratio: " << this->AllowedPixelErrorRatio
<< ", with a total summed threshold error: "
<< vtkm::cont::Algorithm::Reduce(errorPixels, static_cast<FloatDefault>(0)));
auto outputDataSet = this->CreateResultFieldPoint(input, this->GetOutputFieldName(), diffOutput);
outputDataSet.AddPointField(this->GetThresholdFieldName(), thresholdOutput);
return outputDataSet;
}
}
} // namespace filter
} // namespace vtkm

@ -0,0 +1,114 @@
//============================================================================
// 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_image_processing_ImageDifference_h
#define vtk_m_filter_image_processing_ImageDifference_h
#include <vtkm/filter/NewFilterField.h>
#include <vtkm/filter/image_processing/vtkm_filter_image_processing_export.h>
namespace vtkm
{
namespace filter
{
namespace image_processing
{
/// \brief Construct an ImageDifference of a given DataSet
///
/// The dataset generated by executing this filter is a Dataset with two Fields:
/// - "image-diff": Uniform Structured Dataset, difference values of image A - B
/// - "threshold-output": Uniform Structured Dataset, the magnitudes of the pixel differences
///
/// The threshold-output is calculated for each pixel using the `vtkm::Magnitude` vector function
/// on the individual pixel difference.
///
class VTKM_FILTER_IMAGE_PROCESSING_EXPORT ImageDifference : public vtkm::filter::NewFilterField
{
public:
VTKM_CONT ImageDifference();
VTKM_CONT vtkm::IdComponent GetAverageRadius() const { return this->AverageRadius; }
VTKM_CONT void SetAverageRadius(const vtkm::IdComponent& averageRadius)
{
this->AverageRadius = averageRadius;
}
VTKM_CONT vtkm::IdComponent GetPixelShiftRadius() const { return this->PixelShiftRadius; }
VTKM_CONT void SetPixelShiftRadius(const vtkm::IdComponent& pixelShiftRadius)
{
this->PixelShiftRadius = pixelShiftRadius;
}
VTKM_CONT vtkm::FloatDefault GetAllowedPixelErrorRatio() const
{
return this->AllowedPixelErrorRatio;
}
VTKM_CONT void SetAllowedPixelErrorRatio(const vtkm::FloatDefault& pixelErrorRatio)
{
this->AllowedPixelErrorRatio = pixelErrorRatio;
}
VTKM_CONT vtkm::FloatDefault GetPixelDiffThreshold() const { return this->PixelDiffThreshold; }
VTKM_CONT void SetPixelDiffThreshold(const vtkm::FloatDefault& threshold)
{
this->PixelDiffThreshold = threshold;
}
VTKM_CONT bool GetImageDiffWithinThreshold() const { return this->ImageDiffWithinThreshold; }
VTKM_CONT void SetThresholdFieldName(const std::string& name) { this->ThresholdFieldName = name; }
VTKM_CONT std::string GetThresholdFieldName() const { return this->ThresholdFieldName; }
/// Choose the primary field to operate on. For Image difference A - B, A is the
/// primary field.
VTKM_CONT
void SetPrimaryField(
const std::string& name,
vtkm::cont::Field::Association association = vtkm::cont::Field::Association::ANY)
{
this->SetActiveField(name, association);
}
VTKM_CONT std::string GetPrimaryFieldName() const { return this->GetActiveFieldName(); }
VTKM_CONT vtkm::cont::Field::Association GetPrimaryFieldAssociation() const
{
return this->GetActiveFieldAssociation();
}
/// Choose the secondary field to operate on. For Image difference A - B, B is the
/// secondary field.
VTKM_CONT
void SetSecondaryField(
const std::string& name,
vtkm::cont::Field::Association association = vtkm::cont::Field::Association::ANY)
{
this->SetActiveField(1, name, association);
}
VTKM_CONT std::string GetSecondaryFieldName() const { return this->GetActiveFieldName(1); }
VTKM_CONT vtkm::cont::Field::Association GetSecondaryFieldAssociation() const
{
return this->GetActiveFieldAssociation(1);
}
private:
VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& primaryArray) override;
vtkm::IdComponent AverageRadius = 0;
vtkm::IdComponent PixelShiftRadius = 0;
vtkm::FloatDefault AllowedPixelErrorRatio = 0.00025f;
vtkm::FloatDefault PixelDiffThreshold = 0.05f;
bool ImageDiffWithinThreshold = true;
std::string ThresholdFieldName = "threshold-output";
};
} // namespace image_processing
} // namespace filter
} // namespace vtkm
#endif // vtk_m_filter_image_processing_ImageDifference_h

@ -8,17 +8,14 @@
// 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>
#include <vtkm/filter/image_processing/ImageMedian.h>
#include <vtkm/worklet/WorkletPointNeighborhood.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.
//
@ -86,42 +83,45 @@ struct ImageMedian : public vtkm::worklet::WorkletPointNeighborhood
namespace filter
{
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)
namespace image_processing
{
if (!fieldMetadata.IsPointField())
VTKM_CONT vtkm::cont::DataSet ImageMedian::DoExecute(const vtkm::cont::DataSet& input)
{
const auto& field = this->GetFieldFromDataSet(input);
if (!field.IsFieldPoint())
{
throw vtkm::cont::ErrorBadValue("Active field for ImageMedian must be a point field.");
}
const vtkm::cont::UnknownCellSet& 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, *this),
field,
result);
}
else
{
throw vtkm::cont::ErrorBadValue("ImageMedian only support a 3x3 or 5x5 stencil.");
}
const vtkm::cont::UnknownCellSet& inputCellSet = input.GetCellSet();
vtkm::cont::UnknownArrayHandle outArray;
auto resolveType = [&](const auto& concrete) {
// use std::decay to remove const ref from the decltype of concrete.
using T = typename std::decay_t<decltype(concrete)>::ValueType;
vtkm::cont::ArrayHandle<T> result;
if (this->Neighborhood == 1 || this->Neighborhood == 2) // TODO: unnecessary test, see below.
{
this->Invoke(worklet::ImageMedian{ this->Neighborhood }, inputCellSet, concrete, result);
}
else
{
// TODO: this->Neighborhood is already either 1 or 2 by construction. This line is
// unreachable!!! Remove it if Ken is O.K..
throw vtkm::cont::ErrorBadValue("ImageMedian only support a 3x3 or 5x5 stencil.");
}
outArray = result;
};
this->CastAndCallScalarField(field, resolveType);
std::string name = this->GetOutputFieldName();
// TODO: this test is also questionable, didn't we set the output name in the constructor?
if (name.empty())
{
name = fieldMetadata.GetName();
name = field.GetName();
}
return CreateResult(input, fieldMetadata.AsField(name, result));
return this->CreateResultFieldPoint(input, name, outArray);
}
}
}
#endif
} // namespace image_processing
} // namespace filter
} // namespace vtkm

@ -0,0 +1,50 @@
//============================================================================
// 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_image_processing_ImageMedian_h
#define vtk_m_filter_image_processing_ImageMedian_h
#include <vtkm/filter/NewFilterField.h>
#include <vtkm/filter/image_processing/vtkm_filter_image_processing_export.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
{
namespace image_processing
{
class VTKM_FILTER_IMAGE_PROCESSING_EXPORT ImageMedian : public vtkm::filter::NewFilterField
{
public:
VTKM_CONT ImageMedian() { this->SetOutputFieldName("median"); }
VTKM_CONT void Perform3x3() { this->Neighborhood = 1; };
VTKM_CONT void Perform5x5() { this->Neighborhood = 2; };
private:
VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input) override;
int Neighborhood = 1;
};
} // namespace image_processing
} // namespace filter
} // namespace vtkm
#endif //vtk_m_filter_image_processing_ImageMedian_h

@ -0,0 +1,24 @@
##============================================================================
## 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.
##============================================================================
set(unit_tests
UnitTestImageDifferenceFilter.cxx
UnitTestImageMedianFilter.cxx
)
set(libraries
vtkm_filter_image_processing
vtkm_source)
vtkm_unit_tests(
SOURCES ${unit_tests}
LIBRARIES ${libraries}
USE_VTKM_JOB_POOL
)

@ -11,7 +11,7 @@
#include <vtkm/cont/DataSetBuilderUniform.h>
#include <vtkm/cont/Logging.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/ImageDifference.h>
#include <vtkm/filter/image_processing/ImageDifference.h>
#include <random>
#include <vector>
@ -99,7 +99,7 @@ void TestImageDifference()
{
VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Matching Images");
auto dataSet = FillDataSet(static_cast<vtkm::FloatDefault>(1));
vtkm::filter::ImageDifference filter;
vtkm::filter::image_processing::ImageDifference filter;
filter.SetPrimaryField("primary");
filter.SetSecondaryField("secondary");
filter.SetPixelDiffThreshold(0.05f);
@ -122,7 +122,7 @@ void TestImageDifference()
{
VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Matching Images with Average");
auto dataSet = FillDataSet(static_cast<vtkm::FloatDefault>(1));
vtkm::filter::ImageDifference filter;
vtkm::filter::image_processing::ImageDifference filter;
filter.SetPrimaryField("primary");
filter.SetSecondaryField("secondary");
filter.SetPixelDiffThreshold(0.05f);
@ -146,7 +146,7 @@ void TestImageDifference()
{
VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Non Matching Images (Different R pixel)");
auto dataSet = FillDataSet(static_cast<vtkm::FloatDefault>(3));
vtkm::filter::ImageDifference filter;
vtkm::filter::image_processing::ImageDifference filter;
filter.SetPrimaryField("primary");
filter.SetSecondaryField("secondary");
filter.SetPixelDiffThreshold(0.05f);
@ -169,7 +169,7 @@ void TestImageDifference()
{
VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Non Matching Images (Different R pixel)");
auto dataSet = FillDataSet(static_cast<vtkm::FloatDefault>(3));
vtkm::filter::ImageDifference filter;
vtkm::filter::image_processing::ImageDifference filter;
filter.SetPrimaryField("primary");
filter.SetSecondaryField("secondary");
filter.SetPixelDiffThreshold(0.05f);
@ -194,7 +194,7 @@ void TestImageDifference()
VTKM_LOG_S(vtkm::cont::LogLevel::Info,
"Non Matching Images (Different R pixel), Large Threshold");
auto dataSet = FillDataSet(static_cast<vtkm::FloatDefault>(3));
vtkm::filter::ImageDifference filter;
vtkm::filter::image_processing::ImageDifference filter;
filter.SetPrimaryField("primary");
filter.SetSecondaryField("secondary");
filter.SetPixelDiffThreshold(3.0f);

@ -10,7 +10,7 @@
#include <vtkm/cont/testing/MakeTestDataSet.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/ImageMedian.h>
#include <vtkm/filter/image_processing/ImageMedian.h>
namespace
{
@ -22,7 +22,7 @@ void TestImageMedian()
vtkm::cont::testing::MakeTestDataSet testDataSet;
vtkm::cont::DataSet dataSet = testDataSet.Make3DUniformDataSet2();
vtkm::filter::ImageMedian median;
vtkm::filter::image_processing::ImageMedian median;
median.Perform3x3();
median.SetActiveField("pointvar");
auto result = median.Execute(dataSet);

@ -1,5 +1,4 @@
##============================================================================
##
## Copyright (c) Kitware, Inc.
## All rights reserved.
## See LICENSE.txt for details.
@ -7,12 +6,11 @@
## 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.
##
##============================================================================
set(headers
ComputeMoments.h
ImageDifference.h
)
#-----------------------------------------------------------------------------
vtkm_declare_headers(${headers})

@ -17,6 +17,7 @@
#include <vtkm/cont/Field.h>
#include <vtkm/cont/UncertainArrayHandle.h>
#include <vtkm/cont/UncertainCellSet.h>
#include <vtkm/exec/BoundaryState.h>
@ -187,9 +188,9 @@ private:
class ComputeMoments
{
public:
ComputeMoments(const vtkm::Vec3f& _spacing, const double _radius)
: Spacing(_spacing)
, Radius(_radius)
ComputeMoments(double _radius, const vtkm::Vec3f& _spacing)
: Radius(_radius)
, Spacing(_spacing)
{
}
@ -275,8 +276,8 @@ public:
}
private:
const vtkm::FloatDefault Radius = 1;
const vtkm::Vec3f Spacing = { 1, 1, 1 };
const vtkm::FloatDefault Radius;
const vtkm::Vec3f Spacing;
};
}
}

@ -20,8 +20,6 @@ set(unit_tests
UnitTestFieldMetadata.cxx
UnitTestFieldSelection.cxx
UnitTestGhostCellClassify.cxx
UnitTestImageDifferenceFilter.cxx
UnitTestImageMedianFilter.cxx
UnitTestLagrangianFilter.cxx
UnitTestLagrangianStructuresFilter.cxx
UnitTestMapFieldMergeAverage.cxx

@ -7,7 +7,7 @@
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#include <vtkm/cont/ErrorExecution.h>
#include <vtkm/rendering/testing/Testing.h>
TestEqualResult test_equal_images(vtkm::rendering::View& view,
@ -116,7 +116,7 @@ TestEqualResult test_equal_images(const vtkm::cont::DataSet& dataset,
dartXML << "</DartMeasurementFile>\n";
imageDataSet.AddPointField("generated-image", dataset.GetField(0).GetData());
vtkm::filter::ImageDifference filter;
vtkm::filter::image_processing::ImageDifference filter;
filter.SetPrimaryField("baseline-image");
filter.SetSecondaryField("generated-image");
filter.SetAverageRadius(averageRadius);

@ -17,7 +17,7 @@
#include <vtkm/cont/RuntimeDeviceTracker.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/ImageDifference.h>
#include <vtkm/filter/image_processing/ImageDifference.h>
#include <vtkm/internal/Configure.h>
#include <vtkm/io/FileUtils.h>
#include <vtkm/io/ImageUtils.h>

@ -23,7 +23,6 @@ set(headers
DispatcherPointNeighborhood.h
DispatcherReduceByKey.h
FieldStatistics.h
ImageDifference.h
KdTree3D.h # Deprecated
KernelSplatter.h
Keys.h
@ -99,7 +98,6 @@ add_subdirectory(contourtree_distributed)
add_subdirectory(cosmotools)
add_subdirectory(lcs)
add_subdirectory(mir)
add_subdirectory(moments)
add_subdirectory(splatkernels)
add_subdirectory(spatialstructure)
add_subdirectory(tetrahedralize)