From a173b9d9d908c911850ca7780ec7f82f3ad8c286 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?D=C5=BEenan=20Zuki=C4=87?= Date: Thu, 9 Apr 2026 17:34:51 -0400 Subject: [PATCH] WIP: Add StructuralSimilarityImageFilter Salvaged from copilot session: https://github.com/InsightSoftwareConsortium/ITK/agents/pull/6031?session_id=100b799f-7c42-4c54-a76d-4b35dbe7ce91 --- .../itkStructuralSimilarityImageFilter.h | 177 +++++++++++++ .../itkStructuralSimilarityImageFilter.hxx | 238 ++++++++++++++++++ .../ImageCompare/test/CMakeLists.txt | 7 + ...tkStructuralSimilarityImageFilterGTest.cxx | 197 +++++++++++++++ .../itkStructuralSimilarityImageFilter.wrap | 3 + 5 files changed, 622 insertions(+) create mode 100644 Modules/Filtering/ImageCompare/include/itkStructuralSimilarityImageFilter.h create mode 100644 Modules/Filtering/ImageCompare/include/itkStructuralSimilarityImageFilter.hxx create mode 100644 Modules/Filtering/ImageCompare/test/itkStructuralSimilarityImageFilterGTest.cxx create mode 100644 Modules/Filtering/ImageCompare/wrapping/itkStructuralSimilarityImageFilter.wrap diff --git a/Modules/Filtering/ImageCompare/include/itkStructuralSimilarityImageFilter.h b/Modules/Filtering/ImageCompare/include/itkStructuralSimilarityImageFilter.h new file mode 100644 index 00000000000..cea9ae174ab --- /dev/null +++ b/Modules/Filtering/ImageCompare/include/itkStructuralSimilarityImageFilter.h @@ -0,0 +1,177 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkStructuralSimilarityImageFilter_h +#define itkStructuralSimilarityImageFilter_h +#include "itkImageToImageFilter.h" +#include "itkNumericTraits.h" +#include "itkArray.h" +namespace itk +{ +/** + * \class StructuralSimilarityImageFilter + * \brief Computes the Structural Similarity Index Measure (SSIM) between two images. + * + * This filter computes the local SSIM map and an overall scalar SSIM value + * between two input images. The SSIM is defined as: + * + * \f[ + * \text{SSIM}(x,y) = [l(x,y)]^{\alpha} \cdot [c(x,y)]^{\beta} \cdot [s(x,y)]^{\gamma} + * \f] + * + * where + * \f[ + * l(x,y) = \frac{2\mu_x\mu_y + C_1}{\mu_x^2 + \mu_y^2 + C_1}, \quad + * c(x,y) = \frac{2\sigma_x\sigma_y + C_2}{\sigma_x^2 + \sigma_y^2 + C_2}, \quad + * s(x,y) = \frac{\sigma_{xy} + C_3}{\sigma_x\sigma_y + C_3} + * \f] + * + * with \f$C_1 = (K_1 L)^2\f$, \f$C_2 = (K_2 L)^2\f$, \f$C_3 = C_2 / 2\f$. + * \f$L\f$ is the dynamic range (defaults to 255), and \f$K_1\f$ and \f$K_2\f$ + * are small stability constants (defaults 0.01 and 0.03 respectively). + * + * The exponents \f$\alpha\f$, \f$\beta\f$, and \f$\gamma\f$ control the + * relative importance of luminance, contrast, and structure. With the default + * values of 1.0 the formula simplifies to the standard SSIM: + * + * \f[ + * \text{SSIM}(x,y) = \frac{(2\mu_x\mu_y + C_1)(2\sigma_{xy} + C_2)}{(\mu_x^2 + \mu_y^2 + C_1)(\sigma_x^2 + \sigma_y^2 + C_2)} + * \f] + * + * The output image contains the per-pixel SSIM map. The scalar mean SSIM + * across all pixels can be retrieved with GetSSIM(). + * + * The filter is N-dimensional and multi-threaded. + * + * \par Parameters + * - Radius: neighborhood radius for local statistics (default: 1 in each dimension, giving a 3x3 window in 2D). + * - LuminanceWeight (\f$\alpha\f$), ContrastWeight (\f$\beta\f$), StructureWeight (\f$\gamma\f$): + * exponents for the three components (all default to 1.0). + * - DynamicRange (\f$L\f$): range of pixel values (default: 255.0). + * - K1, K2: stabilization constants (defaults: 0.01, 0.03). + * + * \sa SimilarityIndexImageFilter + * + * \ingroup ImageCompare + * \ingroup ITKImageCompare + */ +template +class ITK_TEMPLATE_EXPORT StructuralSimilarityImageFilter : public ImageToImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(StructuralSimilarityImageFilter); + /** Standard class type aliases. */ + using Self = StructuralSimilarityImageFilter; + using Superclass = ImageToImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + /** Method for creation through the object factory. */ + itkNewMacro(Self); + /** \see LightObject::GetNameOfClass() */ + itkOverrideGetNameOfClassMacro(StructuralSimilarityImageFilter); + /** Image type aliases. */ + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + using InputPixelType = typename InputImageType::PixelType; + using OutputPixelType = typename OutputImageType::PixelType; + using InputImageRegionType = typename InputImageType::RegionType; + using OutputImageRegionType = typename OutputImageType::RegionType; + using SizeType = typename InputImageType::SizeType; + using IndexType = typename InputImageType::IndexType; + static constexpr unsigned int ImageDimension = InputImageType::ImageDimension; + /** Floating point type for computations. */ + using RealType = typename NumericTraits::RealType; + /** Neighborhood radius type. */ + using RadiusType = SizeType; + using RadiusValueType = typename RadiusType::SizeValueType; + /** Set/Get the neighborhood radius for local statistics computation. */ + virtual void + SetRadius(const RadiusType & radius); + virtual void + SetRadius(RadiusValueType radius); + itkGetConstReferenceMacro(Radius, RadiusType); + /** Set the first input image. */ + void + SetInput1(const InputImageType * image) + { + this->SetInput(image); + } + /** Set the second input image. */ + void + SetInput2(const InputImageType * image); + /** Get the first input image. */ + const InputImageType * + GetInput1() const + { + return this->GetInput(0); + } + /** Get the second input image. */ + const InputImageType * + GetInput2() const; + /** Set/Get the luminance exponent (alpha). Default: 1.0. */ + itkSetMacro(LuminanceWeight, RealType); + itkGetConstMacro(LuminanceWeight, RealType); + /** Set/Get the contrast exponent (beta). Default: 1.0. */ + itkSetMacro(ContrastWeight, RealType); + itkGetConstMacro(ContrastWeight, RealType); + /** Set/Get the structure exponent (gamma). Default: 1.0. */ + itkSetMacro(StructureWeight, RealType); + itkGetConstMacro(StructureWeight, RealType); + /** Set/Get the dynamic range L of the pixel values. Default: 255.0. */ + itkSetMacro(DynamicRange, RealType); + itkGetConstMacro(DynamicRange, RealType); + /** Set/Get the K1 stabilization constant. Default: 0.01. */ + itkSetMacro(K1, RealType); + itkGetConstMacro(K1, RealType); + /** Set/Get the K2 stabilization constant. Default: 0.03. */ + itkSetMacro(K2, RealType); + itkGetConstMacro(K2, RealType); + /** Get the computed mean SSIM value (available after Update()). */ + itkGetConstMacro(SSIM, RealType); + itkConceptMacro(InputHasNumericTraitsCheck, (Concept::HasNumericTraits)); +protected: + StructuralSimilarityImageFilter(); + ~StructuralSimilarityImageFilter() override = default; + void + PrintSelf(std::ostream & os, Indent indent) const override; + void + GenerateInputRequestedRegion() override; + void + BeforeThreadedGenerateData() override; + void + AfterThreadedGenerateData() override; + void + DynamicThreadedGenerateData(const OutputImageRegionType & outputRegionForThread) override; +private: + RadiusType m_Radius{}; + RealType m_LuminanceWeight{ 1.0 }; + RealType m_ContrastWeight{ 1.0 }; + RealType m_StructureWeight{ 1.0 }; + RealType m_DynamicRange{ 255.0 }; + RealType m_K1{ 0.01 }; + RealType m_K2{ 0.03 }; + RealType m_SSIM{}; + /** Thread-local accumulators for computing the global mean SSIM. */ + std::mutex m_Mutex{}; + RealType m_AccumulatedSSIM{}; + SizeValueType m_PixelCount{}; +}; +} // end namespace itk +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkStructuralSimilarityImageFilter.hxx" +#endif +#endif diff --git a/Modules/Filtering/ImageCompare/include/itkStructuralSimilarityImageFilter.hxx b/Modules/Filtering/ImageCompare/include/itkStructuralSimilarityImageFilter.hxx new file mode 100644 index 00000000000..30d4c6ef80c --- /dev/null +++ b/Modules/Filtering/ImageCompare/include/itkStructuralSimilarityImageFilter.hxx @@ -0,0 +1,238 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkStructuralSimilarityImageFilter_hxx +#define itkStructuralSimilarityImageFilter_hxx +#include "itkBufferedImageNeighborhoodPixelAccessPolicy.h" +#include "itkImageNeighborhoodOffsets.h" +#include "itkImageRegionRange.h" +#include "itkIndexRange.h" +#include "itkNeighborhoodAlgorithm.h" +#include "itkShapedImageNeighborhoodRange.h" +#include "itkZeroFluxNeumannImageNeighborhoodPixelAccessPolicy.h" +#include +#include +namespace itk +{ +template +StructuralSimilarityImageFilter::StructuralSimilarityImageFilter() +{ + this->SetNumberOfRequiredInputs(2); + this->DynamicMultiThreadingOn(); + m_Radius.Fill(1); +} +template +void +StructuralSimilarityImageFilter::SetRadius(const RadiusType & radius) +{ + if (m_Radius != radius) + { + m_Radius = radius; + this->Modified(); + } +} +template +void +StructuralSimilarityImageFilter::SetRadius(RadiusValueType radius) +{ + RadiusType r; + r.Fill(radius); + this->SetRadius(r); +} +template +void +StructuralSimilarityImageFilter::SetInput2(const InputImageType * image) +{ + this->SetNthInput(1, const_cast(image)); +} +template +auto +StructuralSimilarityImageFilter::GetInput2() const -> const InputImageType * +{ + return itkDynamicCastInDebugMode(this->ProcessObject::GetInput(1)); +} +template +void +StructuralSimilarityImageFilter::GenerateInputRequestedRegion() +{ + // Call the superclass' implementation + Superclass::GenerateInputRequestedRegion(); + // Pad the requested region by the neighborhood radius + for (unsigned int i = 0; i < 2; ++i) + { + auto * input = const_cast(this->GetInput(i)); + if (!input) + { + continue; + } + InputImageRegionType requestedRegion = input->GetRequestedRegion(); + requestedRegion.PadByRadius(m_Radius); + requestedRegion.Crop(input->GetLargestPossibleRegion()); + input->SetRequestedRegion(requestedRegion); + } +} +template +void +StructuralSimilarityImageFilter::BeforeThreadedGenerateData() +{ + m_AccumulatedSSIM = RealType{}; + m_PixelCount = 0; +} +template +void +StructuralSimilarityImageFilter::AfterThreadedGenerateData() +{ + if (m_PixelCount > 0) + { + m_SSIM = m_AccumulatedSSIM / static_cast(m_PixelCount); + } + else + { + m_SSIM = RealType{}; + } +} +template +void +StructuralSimilarityImageFilter::DynamicThreadedGenerateData( + const OutputImageRegionType & outputRegionForThread) +{ + const auto * input1 = this->GetInput(0); + const auto * input2 = this->GetInput(1); + auto * output = this->GetOutput(); + const RealType C1 = (m_K1 * m_DynamicRange) * (m_K1 * m_DynamicRange); + const RealType C2 = (m_K2 * m_DynamicRange) * (m_K2 * m_DynamicRange); + const RealType C3 = C2 / RealType{ 2.0 }; + const bool useSimplifiedFormula = + (itk::Math::FloatAlmostEqual(m_LuminanceWeight, RealType{ 1.0 }) && + itk::Math::FloatAlmostEqual(m_ContrastWeight, RealType{ 1.0 }) && + itk::Math::FloatAlmostEqual(m_StructureWeight, RealType{ 1.0 })); + const auto neighborhoodOffsets = GenerateRectangularImageNeighborhoodOffsets(m_Radius); + const auto neighborhoodSize = static_cast(neighborhoodOffsets.size()); + // Use ImageBoundaryFacesCalculator to split into interior and boundary regions + const auto calculatorResult = + NeighborhoodAlgorithm::ImageBoundaryFacesCalculator::Compute(*input1, outputRegionForThread, m_Radius); + // Lambda to process a subregion with a given pixel access policy + auto processSubregion = [&](const auto & subregion, auto policyTag1, auto policyTag2) + { + using Policy1 = typename decltype(policyTag1)::type; + using Policy2 = typename decltype(policyTag2)::type; + auto range1 = ShapedImageNeighborhoodRange( + *input1, IndexType(), neighborhoodOffsets); + auto range2 = ShapedImageNeighborhoodRange( + *input2, IndexType(), neighborhoodOffsets); + auto outputIt = ImageRegionRange(*output, subregion).begin(); + RealType localSum{}; + SizeValueType localCount{}; + for (const auto & index : MakeIndexRange(subregion)) + { + range1.SetLocation(index); + range2.SetLocation(index); + // Compute local means + RealType sumX{}; + RealType sumY{}; + RealType sumXX{}; + RealType sumYY{}; + RealType sumXY{}; + auto it1 = range1.begin(); + auto it2 = range2.begin(); + for (; it1 != range1.end(); ++it1, ++it2) + { + const RealType x = static_cast(*it1); + const RealType y = static_cast(*it2); + sumX += x; + sumY += y; + sumXX += x * x; + sumYY += y * y; + sumXY += x * y; + } + const RealType muX = sumX / neighborhoodSize; + const RealType muY = sumY / neighborhoodSize; + // Variance and covariance (using population formula consistent with SSIM literature) + const RealType sigmaXX = sumXX / neighborhoodSize - muX * muX; + const RealType sigmaYY = sumYY / neighborhoodSize - muY * muY; + const RealType sigmaXY = sumXY / neighborhoodSize - muX * muY; + RealType ssimValue; + if (useSimplifiedFormula) + { + // Standard SSIM: (2*muX*muY + C1)*(2*sigmaXY + C2) / + // ((muX^2 + muY^2 + C1)*(sigmaXX + sigmaYY + C2)) + const RealType numerator = (RealType{ 2.0 } * muX * muY + C1) * (RealType{ 2.0 } * sigmaXY + C2); + const RealType denominator = (muX * muX + muY * muY + C1) * (sigmaXX + sigmaYY + C2); + ssimValue = numerator / denominator; + } + else + { + // General form with individual exponents + const RealType luminance = + (RealType{ 2.0 } * muX * muY + C1) / (muX * muX + muY * muY + C1); + const RealType sigmaX = std::sqrt(std::max(sigmaXX, RealType{})); + const RealType sigmaY = std::sqrt(std::max(sigmaYY, RealType{})); + const RealType contrast = + (RealType{ 2.0 } * sigmaX * sigmaY + C2) / (sigmaXX + sigmaYY + C2); + const RealType structure = (sigmaXY + C3) / (sigmaX * sigmaY + C3); + ssimValue = + std::pow(luminance, m_LuminanceWeight) * std::pow(contrast, m_ContrastWeight) * std::pow(structure, m_StructureWeight); + } + *outputIt = static_cast(ssimValue); + ++outputIt; + localSum += ssimValue; + ++localCount; + } + // Accumulate thread-local results under the mutex + { + const std::lock_guard lock(m_Mutex); + m_AccumulatedSSIM += localSum; + m_PixelCount += localCount; + } + }; + // Tag types to pass policy as template argument through the lambda + struct PolicyTag1 + { + using type = BufferedImageNeighborhoodPixelAccessPolicy; + }; + struct PolicyTag2 + { + using type = ZeroFluxNeumannImageNeighborhoodPixelAccessPolicy; + }; + // Process interior region (fast path, no boundary checks) + const auto & nonBoundaryRegion = calculatorResult.GetNonBoundaryRegion(); + if (nonBoundaryRegion.GetNumberOfPixels() > 0) + { + processSubregion(nonBoundaryRegion, PolicyTag1{}, PolicyTag1{}); + } + // Process boundary faces (with boundary extrapolation) + for (const auto & face : calculatorResult.GetBoundaryFaces()) + { + processSubregion(face, PolicyTag2{}, PolicyTag2{}); + } +} +template +void +StructuralSimilarityImageFilter::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + os << indent << "Radius: " << m_Radius << std::endl; + os << indent << "LuminanceWeight: " << m_LuminanceWeight << std::endl; + os << indent << "ContrastWeight: " << m_ContrastWeight << std::endl; + os << indent << "StructureWeight: " << m_StructureWeight << std::endl; + os << indent << "DynamicRange: " << m_DynamicRange << std::endl; + os << indent << "K1: " << m_K1 << std::endl; + os << indent << "K2: " << m_K2 << std::endl; + os << indent << "SSIM: " << m_SSIM << std::endl; +} +} // end namespace itk +#endif diff --git a/Modules/Filtering/ImageCompare/test/CMakeLists.txt b/Modules/Filtering/ImageCompare/test/CMakeLists.txt index a7caaba4f4f..0122af87a2f 100644 --- a/Modules/Filtering/ImageCompare/test/CMakeLists.txt +++ b/Modules/Filtering/ImageCompare/test/CMakeLists.txt @@ -13,6 +13,13 @@ set( createtestdriver(ITKImageCompare "${ITKImageCompare-Test_LIBRARIES}" "${ITKImageCompareTests}") +set( + ITKImageCompareGTests + itkStructuralSimilarityImageFilterGTest.cxx +) +creategoogletestdriver(ITKImageCompare "${ITKImageCompare-Test_LIBRARIES}" "${ITKImageCompareGTests}") + + itk_add_test( NAME itkAbsoluteValueDifferenceImageFilterTest COMMAND diff --git a/Modules/Filtering/ImageCompare/test/itkStructuralSimilarityImageFilterGTest.cxx b/Modules/Filtering/ImageCompare/test/itkStructuralSimilarityImageFilterGTest.cxx new file mode 100644 index 00000000000..6974c8c17cc --- /dev/null +++ b/Modules/Filtering/ImageCompare/test/itkStructuralSimilarityImageFilterGTest.cxx @@ -0,0 +1,197 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#include "itkGTest.h" +#include "itkImage.h" +#include "itkStructuralSimilarityImageFilter.h" +namespace +{ +using PixelType = float; +constexpr unsigned int Dimension = 2; +using ImageType = itk::Image; +using FilterType = itk::StructuralSimilarityImageFilter; +ImageType::Pointer +CreateConstantImage(PixelType value, unsigned int size = 32) +{ + auto image = ImageType::New(); + ImageType::SizeType imageSize; + imageSize.Fill(size); + const ImageType::RegionType region(imageSize); + image->SetRegions(region); + image->Allocate(); + image->FillBuffer(value); + return image; +} +ImageType::Pointer +CreateGradientImage(unsigned int size = 32) +{ + auto image = ImageType::New(); + ImageType::SizeType imageSize; + imageSize.Fill(size); + const ImageType::RegionType region(imageSize); + image->SetRegions(region); + image->Allocate(); + itk::ImageRegionIteratorWithIndex it(image, region); + for (it.GoToBegin(); !it.IsAtEnd(); ++it) + { + const auto idx = it.GetIndex(); + it.Set(static_cast(idx[0] + idx[1])); + } + return image; +} +} // namespace +TEST(StructuralSimilarityImageFilter, BasicObjectMethods) +{ + auto filter = FilterType::New(); + ITK_GTEST_EXERCISE_BASIC_OBJECT_METHODS(filter, StructuralSimilarityImageFilter, ImageToImageFilter); +} +TEST(StructuralSimilarityImageFilter, IdenticalImagesGiveSSIMOfOne) +{ + auto image = CreateConstantImage(100.0f); + auto filter = FilterType::New(); + filter->SetInput1(image); + filter->SetInput2(image); + filter->Update(); + EXPECT_NEAR(filter->GetSSIM(), 1.0, 1e-6); +} +TEST(StructuralSimilarityImageFilter, IdenticalGradientImagesGiveSSIMOfOne) +{ + auto image = CreateGradientImage(); + auto filter = FilterType::New(); + filter->SetInput1(image); + filter->SetInput2(image); + filter->Update(); + EXPECT_NEAR(filter->GetSSIM(), 1.0, 1e-6); +} +TEST(StructuralSimilarityImageFilter, DifferentConstantImagesHaveLowSSIM) +{ + auto image1 = CreateConstantImage(0.0f); + auto image2 = CreateConstantImage(255.0f); + auto filter = FilterType::New(); + filter->SetInput1(image1); + filter->SetInput2(image2); + filter->Update(); + // Two maximally different constant images should have a very low SSIM + EXPECT_LT(filter->GetSSIM(), 0.01); +} +TEST(StructuralSimilarityImageFilter, OutputMapHasCorrectSize) +{ + constexpr unsigned int imageSize = 16; + auto image1 = CreateConstantImage(100.0f, imageSize); + auto image2 = CreateConstantImage(100.0f, imageSize); + auto filter = FilterType::New(); + filter->SetInput1(image1); + filter->SetInput2(image2); + filter->Update(); + auto output = filter->GetOutput(); + EXPECT_EQ(output->GetLargestPossibleRegion().GetSize()[0], imageSize); + EXPECT_EQ(output->GetLargestPossibleRegion().GetSize()[1], imageSize); +} +TEST(StructuralSimilarityImageFilter, SetGetParameters) +{ + auto filter = FilterType::New(); + filter->SetLuminanceWeight(2.0); + EXPECT_DOUBLE_EQ(filter->GetLuminanceWeight(), 2.0); + filter->SetContrastWeight(3.0); + EXPECT_DOUBLE_EQ(filter->GetContrastWeight(), 3.0); + filter->SetStructureWeight(0.5); + EXPECT_DOUBLE_EQ(filter->GetStructureWeight(), 0.5); + filter->SetDynamicRange(1.0); + EXPECT_DOUBLE_EQ(filter->GetDynamicRange(), 1.0); + filter->SetK1(0.05); + EXPECT_DOUBLE_EQ(filter->GetK1(), 0.05); + filter->SetK2(0.1); + EXPECT_DOUBLE_EQ(filter->GetK2(), 0.1); + FilterType::RadiusType radius; + radius.Fill(3); + filter->SetRadius(radius); + EXPECT_EQ(filter->GetRadius()[0], 3u); + EXPECT_EQ(filter->GetRadius()[1], 3u); + filter->SetRadius(5u); + EXPECT_EQ(filter->GetRadius()[0], 5u); + EXPECT_EQ(filter->GetRadius()[1], 5u); +} +TEST(StructuralSimilarityImageFilter, ScalarRadiusSetGet) +{ + auto filter = FilterType::New(); + filter->SetRadius(4u); + for (unsigned int d = 0; d < Dimension; ++d) + { + EXPECT_EQ(filter->GetRadius()[d], 4u); + } +} +TEST(StructuralSimilarityImageFilter, SlightDifferenceGivesHighSSIM) +{ + auto image1 = CreateGradientImage(32); + auto image2 = CreateGradientImage(32); + // Perturb one pixel slightly + ImageType::IndexType idx; + idx.Fill(16); + image2->SetPixel(idx, image2->GetPixel(idx) + 1.0f); + auto filter = FilterType::New(); + filter->SetInput1(image1); + filter->SetInput2(image2); + filter->Update(); + // A single-pixel perturbation in a gradient image should still yield high SSIM + EXPECT_GT(filter->GetSSIM(), 0.99); +} +TEST(StructuralSimilarityImageFilter, LargerRadiusWorks) +{ + auto image = CreateGradientImage(32); + auto filter = FilterType::New(); + filter->SetInput1(image); + filter->SetInput2(image); + filter->SetRadius(3u); + filter->Update(); + EXPECT_NEAR(filter->GetSSIM(), 1.0, 1e-6); +} +TEST(StructuralSimilarityImageFilter, CustomWeightsWork) +{ + auto image1 = CreateGradientImage(32); + auto image2 = CreateGradientImage(32); + // Add noise + ImageType::IndexType idx; + idx[0] = 10; + idx[1] = 10; + image2->SetPixel(idx, image2->GetPixel(idx) + 50.0f); + auto filter = FilterType::New(); + filter->SetInput1(image1); + filter->SetInput2(image2); + filter->SetLuminanceWeight(1.0); + filter->SetContrastWeight(1.0); + filter->SetStructureWeight(1.0); + filter->Update(); + // Should produce a valid SSIM between 0 and 1 + EXPECT_GT(filter->GetSSIM(), 0.0); + EXPECT_LE(filter->GetSSIM(), 1.0); +} +TEST(StructuralSimilarityImageFilter, ThreeDimensionalImage) +{ + using Image3DType = itk::Image; + using Filter3DType = itk::StructuralSimilarityImageFilter; + auto image = Image3DType::New(); + Image3DType::SizeType size; + size.Fill(8); + image->SetRegions(Image3DType::RegionType(size)); + image->Allocate(); + image->FillBuffer(42.0f); + auto filter = Filter3DType::New(); + filter->SetInput1(image); + filter->SetInput2(image); + filter->Update(); + EXPECT_NEAR(filter->GetSSIM(), 1.0, 1e-6); +} diff --git a/Modules/Filtering/ImageCompare/wrapping/itkStructuralSimilarityImageFilter.wrap b/Modules/Filtering/ImageCompare/wrapping/itkStructuralSimilarityImageFilter.wrap new file mode 100644 index 00000000000..53ed5d34a26 --- /dev/null +++ b/Modules/Filtering/ImageCompare/wrapping/itkStructuralSimilarityImageFilter.wrap @@ -0,0 +1,3 @@ +itk_wrap_class("itk::StructuralSimilarityImageFilter" POINTER) +itk_wrap_image_filter("${WRAP_ITK_SCALAR}" 2) +itk_end_wrap_class()