diff --git a/CMakeLists.txt b/CMakeLists.txt index aba2ec14dd4..5ff1a0d28ef 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -741,6 +741,12 @@ endif() # dependencies. add_subdirectory(Modules/Remote) +# In-tree beta modules (formerly remote fetches, now grafted inline). +# Loads provenance manifests (Modules/Beta/*.beta.cmake) and is a no-op +# at build time; the actual modules under Modules/Beta// are +# discovered via the normal itk-module.cmake DAG scan below. +add_subdirectory(Modules/Beta) + # Enable modules according to user inputs and the module dependency DAG. include(ITKModuleEnablement) diff --git a/Modules/Beta/CMakeLists.txt b/Modules/Beta/CMakeLists.txt new file mode 100644 index 00000000000..6988fe7cada --- /dev/null +++ b/Modules/Beta/CMakeLists.txt @@ -0,0 +1,29 @@ +# Modules/Beta/ +# +# In-tree home for modules that were formerly configure-time remote +# fetches (Modules/Remote/*.remote.cmake). Each module lives at +# Modules/Beta// with its own itk-module.cmake and is discovered +# automatically by the normal ITK module DAG scan (see +# CMake/ITKModuleEnablement.cmake). +# +# Alongside each ingested module, an optional manifest +# Modules/Beta/.beta.cmake records upstream provenance +# (UPSTREAM_URL, UPSTREAM_SHA at last ingest, INGEST_DATE, license, +# compliance level). The manifest files are inert at build time -- +# they are passive metadata consumed by tooling and human readers. + +# itk_beta_module_manifest() is a no-op that makes the *.beta.cmake +# files valid CMake so they can be include()-d by tooling without +# producing build-time side effects. +function(itk_beta_module_manifest _name) + # arguments are intentionally discarded; parsed by external tooling. +endfunction() + +# Include every *.beta.cmake so CMake at least syntactically validates +# them at configure time. This does not enable, build, or fetch anything. +file(GLOB _beta_manifests "${CMAKE_CURRENT_SOURCE_DIR}/*.beta.cmake") +foreach(_manifest ${_beta_manifests}) + include(${_manifest}) +endforeach() +unset(_beta_manifests) +unset(_manifest) diff --git a/Modules/Beta/StructuralSimilarity.beta.cmake b/Modules/Beta/StructuralSimilarity.beta.cmake new file mode 100644 index 00000000000..3479f280486 --- /dev/null +++ b/Modules/Beta/StructuralSimilarity.beta.cmake @@ -0,0 +1,13 @@ +itk_beta_module_manifest( + NAME StructuralSimilarity + UPSTREAM_URL "" + UPSTREAM_SHA "" + INGEST_DATE 2026-04-19 + LICENSE "Apache-2.0" + COMPLIANCE_LEVEL 3 + DESCRIPTION + "Structural Similarity Index (SSIM) image filter. New module +originally proposed against Modules/Filtering/ImageCompare; relocated to +Modules/Beta/ so the API can stabilize before promotion to a core +module." +) diff --git a/Modules/Beta/StructuralSimilarity/CMakeLists.txt b/Modules/Beta/StructuralSimilarity/CMakeLists.txt new file mode 100644 index 00000000000..6e1c10bdb77 --- /dev/null +++ b/Modules/Beta/StructuralSimilarity/CMakeLists.txt @@ -0,0 +1,2 @@ +project(StructuralSimilarity) +itk_module_impl() diff --git a/Modules/Beta/StructuralSimilarity/include/itkStructuralSimilarityImageFilter.h b/Modules/Beta/StructuralSimilarity/include/itkStructuralSimilarityImageFilter.h new file mode 100644 index 00000000000..d318bb79f4f --- /dev/null +++ b/Modules/Beta/StructuralSimilarity/include/itkStructuralSimilarityImageFilter.h @@ -0,0 +1,299 @@ +/*========================================================================= + * + * 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 "itkArray.h" +#include "itkImage.h" +#include "itkImageToImageFilter.h" +#include "itkNumericTraits.h" + +#include + +namespace itk +{ +/** + * \class StructuralSimilarityImageFilter + * \brief Computes the Structural Similarity Index Measure (SSIM) between two images. + * + * This filter computes the Structural Similarity Index Measure + * \cite wang2004image between two input images of identical geometry. The + * output image stores the per-pixel SSIM map. The scalar mean SSIM over the + * valid (non-boundary) region is available via GetMeanSSIM() after Update(). + * + * \par Algorithm + * For two images \f$x\f$ and \f$y\f$, local statistics are computed by + * convolving with a discrete Gaussian kernel of standard deviation + * \f$\sigma\f$ (default 1.5): + * \f[ + * \mu_x = G_\sigma * x, \quad \mu_y = G_\sigma * y, + * \f] + * \f[ + * \sigma_x^2 = G_\sigma * x^2 - \mu_x^2,\quad + * \sigma_y^2 = G_\sigma * y^2 - \mu_y^2,\quad + * \sigma_{xy} = G_\sigma * (xy) - \mu_x \mu_y . + * \f] + * + * The three SSIM components are + * \f[ + * l(x,y) = \frac{2\mu_x\mu_y + C_1}{\mu_x^2 + \mu_y^2 + C_1}, \qquad + * c(x,y) = \frac{2\sigma_x\sigma_y + C_2}{\sigma_x^2 + \sigma_y^2 + C_2}, \qquad + * 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$, + * and \f$L\f$ the dynamic range of the pixel values. + * + * The combined SSIM is + * \f[ + * \mathrm{SSIM}(x,y) = [l(x,y)]^{\alpha}\,[c(x,y)]^{\beta}\,[s(x,y)]^{\gamma}. + * \f] + * + * With the default exponents \f$\alpha = \beta = \gamma = 1\f$ and the + * convention \f$C_3 = C_2/2\f$, this collapses to the simplified form + * \f[ + * \mathrm{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] + * which matches the reference implementation distributed by Wang et al. + * and the default behavior of \c skimage.metrics.structural_similarity . + * + * \par Properties + * - For identical images, the per-pixel SSIM is exactly 1 and the mean SSIM + * is exactly 1 (subject to floating-point precision). + * - The SSIM index is symmetric: \f$\mathrm{SSIM}(x,y) = \mathrm{SSIM}(y,x)\f$. + * - The SSIM index is bounded above by 1. In typical cases it is + * non-negative; values can be slightly negative for anti-correlated + * regions. + * + * \par Parameters + * - \c GaussianSigma: standard deviation of the Gaussian window + * (default 1.5, matching Wang et al.). + * - \c MaximumKernelWidth: hard limit on the discrete Gaussian kernel width + * (default 11, matching the canonical 11x11 window). + * - \c K1, \c K2: stability constants (defaults 0.01 and 0.03). + * - \c DynamicRange: \f$L\f$ in the formulas above; defaults to the dynamic + * range of the input pixel type via NumericTraits (e.g. 255 for + * \c unsigned char, 1.0 for \c float / \c double). For arbitrary + * floating-point images, set this explicitly to the actual data range. + * - \c LuminanceExponent (\f$\alpha\f$), \c ContrastExponent (\f$\beta\f$), + * \c StructureExponent (\f$\gamma\f$): defaults all 1.0. + * - \c ScaleWeights: array of per-scale weights for multi-scale SSIM + * (MS-SSIM, \cite wang2003multiscale). When the array contains a single + * element (the default), the filter computes ordinary single-scale SSIM. + * Multi-scale evaluation with more than one scale is not yet implemented + * and will raise an exception in BeforeGenerate. + * + * The filter is N-dimensional, multi-threaded, and templated over the input + * and output image types. The output pixel type defaults to \c float. + * + * \sa SimilarityIndexImageFilter + * \sa DiscreteGaussianImageFilter + * + * \ingroup MultiThreaded + * \ingroup StructuralSimilarity + */ +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 used for all SSIM computations. */ + using RealType = typename NumericTraits::RealType; + + /** Type used for the user-specified array of multi-scale weights. */ + using ScaleWeightsType = Array; + + /** Set/Get the first input image. */ + /** @ITKStartGrouping */ + void + SetInput1(const InputImageType * image) + { + this->SetInput(image); + } + const InputImageType * + GetInput1() const + { + return this->GetInput(0); + } + /** @ITKEndGrouping */ + + /** Set/Get the second input image. */ + /** @ITKStartGrouping */ + void + SetInput2(const InputImageType * image); + const InputImageType * + GetInput2() const; + /** @ITKEndGrouping */ + + /** Standard deviation \f$\sigma\f$ of the Gaussian window used to compute + * local statistics. Default 1.5 (matching Wang et al. 2004). */ + /** @ITKStartGrouping */ + itkSetMacro(GaussianSigma, double); + itkGetConstMacro(GaussianSigma, double); + /** @ITKEndGrouping */ + + /** Maximum width (per dimension) of the discrete Gaussian kernel. + * Default 11, giving an 11x11 window in 2D when sigma=1.5. */ + /** @ITKStartGrouping */ + itkSetMacro(MaximumKernelWidth, unsigned int); + itkGetConstMacro(MaximumKernelWidth, unsigned int); + /** @ITKEndGrouping */ + + /** \f$K_1\f$ stability constant. Default 0.01. */ + /** @ITKStartGrouping */ + itkSetMacro(K1, double); + itkGetConstMacro(K1, double); + /** @ITKEndGrouping */ + + /** \f$K_2\f$ stability constant. Default 0.03. */ + /** @ITKStartGrouping */ + itkSetMacro(K2, double); + itkGetConstMacro(K2, double); + /** @ITKEndGrouping */ + + /** Dynamic range \f$L\f$ of the pixel values used to compute + * \f$C_1 = (K_1 L)^2\f$ and \f$C_2 = (K_2 L)^2\f$. Default depends on + * the input pixel type: 255 for \c unsigned \c char, 65535 for + * \c unsigned \c short, 1.0 for \c float / \c double, etc. */ + /** @ITKStartGrouping */ + itkSetMacro(DynamicRange, double); + itkGetConstMacro(DynamicRange, double); + /** @ITKEndGrouping */ + + /** Exponent \f$\alpha\f$ on the luminance term. Default 1.0. */ + /** @ITKStartGrouping */ + itkSetMacro(LuminanceExponent, double); + itkGetConstMacro(LuminanceExponent, double); + /** @ITKEndGrouping */ + + /** Exponent \f$\beta\f$ on the contrast term. Default 1.0. */ + /** @ITKStartGrouping */ + itkSetMacro(ContrastExponent, double); + itkGetConstMacro(ContrastExponent, double); + /** @ITKEndGrouping */ + + /** Exponent \f$\gamma\f$ on the structure term. Default 1.0. */ + /** @ITKStartGrouping */ + itkSetMacro(StructureExponent, double); + itkGetConstMacro(StructureExponent, double); + /** @ITKEndGrouping */ + + /** Per-scale weights for multi-scale SSIM (MS-SSIM). An array of size 1 + * (the default) requests ordinary single-scale SSIM and is the only + * configuration currently supported. Setting an array of length greater + * than 1 will currently raise an exception in BeforeGenerate. */ + /** @ITKStartGrouping */ + void + SetScaleWeights(const ScaleWeightsType & weights); + itkGetConstReferenceMacro(ScaleWeights, ScaleWeightsType); + /** @ITKEndGrouping */ + + /** Mean SSIM over the valid (non-Gaussian-padded) region. Available + * after Update(). */ + itkGetConstMacro(MeanSSIM, double); + + itkConceptMacro(InputHasNumericTraitsCheck, (Concept::HasNumericTraits)); + +protected: + StructuralSimilarityImageFilter(); + ~StructuralSimilarityImageFilter() override = default; + + void + PrintSelf(std::ostream & os, Indent indent) const override; + + /** Verify that parameters are valid and both inputs are set with matching regions. */ + void + VerifyPreconditions() const override; + + /** This filter requires its full input to compute correct mean SSIM. */ + void + GenerateInputRequestedRegion() override; + + /** This filter computes the full output. */ + void + EnlargeOutputRequestedRegion(DataObject * data) override; + + /** Composite-filter-style: drives the internal sub-pipeline (5 Gaussian + * convolutions plus a parallelized SSIM combination). */ + void + GenerateData() override; + +private: + double m_GaussianSigma{ 1.5 }; + unsigned int m_MaximumKernelWidth{ 11 }; + double m_K1{ 0.01 }; + double m_K2{ 0.03 }; + + /** Default dynamic range: 1.0 for floating-point pixels (assume normalized + * data), and \c NumericTraits::max() - \c NumericTraits::min() for integer + * pixels (e.g. 255 for \c unsigned \c char). */ + static constexpr double + DefaultDynamicRange() + { + if constexpr (std::is_floating_point_v) + { + return 1.0; + } + else + { + return static_cast(NumericTraits::max()) - + static_cast(NumericTraits::min()); + } + } + double m_DynamicRange{ DefaultDynamicRange() }; + double m_LuminanceExponent{ 1.0 }; + double m_ContrastExponent{ 1.0 }; + double m_StructureExponent{ 1.0 }; + + ScaleWeightsType m_ScaleWeights{}; + + double m_MeanSSIM{ 0.0 }; +}; +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkStructuralSimilarityImageFilter.hxx" +#endif + +#endif diff --git a/Modules/Beta/StructuralSimilarity/include/itkStructuralSimilarityImageFilter.hxx b/Modules/Beta/StructuralSimilarity/include/itkStructuralSimilarityImageFilter.hxx new file mode 100644 index 00000000000..9ea70218c85 --- /dev/null +++ b/Modules/Beta/StructuralSimilarity/include/itkStructuralSimilarityImageFilter.hxx @@ -0,0 +1,410 @@ +/*========================================================================= + * + * 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 "itkCastImageFilter.h" +#include "itkDiscreteGaussianImageFilter.h" +#include "itkImage.h" +#include "itkImageRegionConstIterator.h" +#include "itkImageRegionIteratorWithIndex.h" +#include "itkMath.h" +#include "itkMultiplyImageFilter.h" +#include "itkProgressReporter.h" +#include "itkTotalProgressReporter.h" + +#include +#include +#include + +namespace itk +{ +template +StructuralSimilarityImageFilter::StructuralSimilarityImageFilter() +{ + this->SetNumberOfRequiredInputs(2); + + // Default ScaleWeights = [1.0] => single-scale SSIM. + m_ScaleWeights.SetSize(1); + m_ScaleWeights[0] = static_cast(1.0); +} + +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::SetScaleWeights(const ScaleWeightsType & weights) +{ + if (m_ScaleWeights != weights) + { + m_ScaleWeights = weights; + this->Modified(); + } +} + +template +void +StructuralSimilarityImageFilter::GenerateInputRequestedRegion() +{ + // The mean SSIM is computed over the full output, so the filter needs the + // entire largest possible region of both inputs. This mirrors + // SimilarityIndexImageFilter and avoids subtle errors when the filter is + // pipelined behind a streaming filter. + Superclass::GenerateInputRequestedRegion(); + + for (unsigned int i = 0; i < 2; ++i) + { + auto * input = const_cast(this->GetInput(i)); + if (input) + { + input->SetRequestedRegionToLargestPossibleRegion(); + } + } +} + +template +void +StructuralSimilarityImageFilter::EnlargeOutputRequestedRegion(DataObject * data) +{ + Superclass::EnlargeOutputRequestedRegion(data); + data->SetRequestedRegionToLargestPossibleRegion(); +} + +template +void +StructuralSimilarityImageFilter::VerifyPreconditions() const +{ + Superclass::VerifyPreconditions(); + + if (m_ScaleWeights.GetSize() == 0) + { + itkExceptionMacro("ScaleWeights array must contain at least one element."); + } + if (m_ScaleWeights.GetSize() > 1) + { + itkExceptionMacro("Multi-scale SSIM (MS-SSIM) with ScaleWeights size > 1 is not yet implemented. " + "Set a ScaleWeights array of length 1 to compute single-scale SSIM."); + } + if (m_GaussianSigma <= 0.0) + { + itkExceptionMacro("GaussianSigma must be strictly positive (got " << m_GaussianSigma << ")."); + } + if (m_DynamicRange <= 0.0) + { + itkExceptionMacro("DynamicRange must be strictly positive (got " << m_DynamicRange << ")."); + } + + const InputImageType * input1 = this->GetInput1(); + const InputImageType * input2 = this->GetInput2(); + if (input1 == nullptr || input2 == nullptr) + { + itkExceptionMacro("StructuralSimilarityImageFilter requires both inputs to be set."); + } + + if (input1->GetLargestPossibleRegion() != input2->GetLargestPossibleRegion()) + { + itkExceptionMacro("StructuralSimilarityImageFilter requires the two inputs to have identical regions."); + } +} + +template +void +StructuralSimilarityImageFilter::GenerateData() +{ + // ---- Build the internal Gaussian-statistics pipeline ----------------- + // + // Following the composite-filter pattern (see CompositeFilterExample.cxx), + // we reconnect the mini-pipeline inputs from the external pipeline every + // time GenerateData() is called. This ensures the internal sub-filters + // pick up any new inputs when the external pipeline re-executes. + using RealImageType = Image; + using CastFilterType = CastImageFilter; + using MultiplyFilterType = MultiplyImageFilter; + using GaussianFilterType = DiscreteGaussianImageFilter; + + // Graft external inputs to disconnect the mini-pipeline from the + // upstream pipeline, preventing the internal Update() calls from + // propagating back upstream. + auto graftedInput1 = InputImageType::New(); + graftedInput1->Graft(this->GetInput1()); + + auto graftedInput2 = InputImageType::New(); + graftedInput2->Graft(this->GetInput2()); + + auto cast1 = CastFilterType::New(); + cast1->SetInput(graftedInput1); + + auto cast2 = CastFilterType::New(); + cast2->SetInput(graftedInput2); + + auto x_times_x = MultiplyFilterType::New(); + x_times_x->SetInput1(cast1->GetOutput()); + x_times_x->SetInput2(cast1->GetOutput()); + + auto y_times_y = MultiplyFilterType::New(); + y_times_y->SetInput1(cast2->GetOutput()); + y_times_y->SetInput2(cast2->GetOutput()); + + auto x_times_y = MultiplyFilterType::New(); + x_times_y->SetInput1(cast1->GetOutput()); + x_times_y->SetInput2(cast2->GetOutput()); + + // Configure a Gaussian filter with the requested sigma and kernel-width + // ceiling. We disable image-spacing-aware sigmas: SSIM is defined in + // pixel coordinates and the canonical reference (Wang 2004) uses an + // 11x11 unit-spacing window. + auto makeGaussian = [this](const auto & inputPort) { + auto g = GaussianFilterType::New(); + g->SetInput(inputPort); + g->SetSigma(m_GaussianSigma); + g->SetMaximumKernelWidth(m_MaximumKernelWidth); + g->SetUseImageSpacing(false); + g->ReleaseDataFlagOn(); + return g; + }; + + auto g_x = makeGaussian(cast1->GetOutput()); + auto g_y = makeGaussian(cast2->GetOutput()); + auto g_xx = makeGaussian(x_times_x->GetOutput()); + auto g_yy = makeGaussian(y_times_y->GetOutput()); + auto g_xy = makeGaussian(x_times_y->GetOutput()); + + g_x->Update(); + g_y->Update(); + g_xx->Update(); + g_yy->Update(); + g_xy->Update(); + + // ---- Allocate output -------------------------------------------------- + OutputImageType * output = this->GetOutput(); + output->SetBufferedRegion(output->GetRequestedRegion()); + output->Allocate(); + + const auto outputRegion = output->GetRequestedRegion(); + + // ---- Pre-compute SSIM constants -------------------------------------- + const RealType K1 = static_cast(m_K1); + const RealType K2 = static_cast(m_K2); + const RealType L = static_cast(m_DynamicRange); + const RealType C1 = (K1 * L) * (K1 * L); + const RealType C2 = (K2 * L) * (K2 * L); + const RealType C3 = C2 / static_cast(2); + + const bool useSimplifiedFormula = Math::FloatAlmostEqual(m_LuminanceExponent, 1.0) && + Math::FloatAlmostEqual(m_ContrastExponent, 1.0) && + Math::FloatAlmostEqual(m_StructureExponent, 1.0); + + const auto alpha = static_cast(m_LuminanceExponent); + const auto beta = static_cast(m_ContrastExponent); + const auto gamma = static_cast(m_StructureExponent); + + const RealImageType * gx = g_x->GetOutput(); + const RealImageType * gy = g_y->GetOutput(); + const RealImageType * gxx = g_xx->GetOutput(); + const RealImageType * gyy = g_yy->GetOutput(); + const RealImageType * gxy = g_xy->GetOutput(); + + // ---- Determine the "valid" interior region for mean SSIM ------------ + // + // Pixels within the half-width of the discrete Gaussian kernel of the + // image boundary use boundary-extended values inside the convolution and + // are therefore less reliable. scikit-image (matching the MATLAB + // reference) crops by (win_size - 1)/2 before averaging. We do the same. + const auto kernelHalfWidth = static_cast(m_MaximumKernelWidth / 2); + auto interiorRegion = outputRegion; + bool interiorIsValid = true; + for (unsigned int d = 0; d < ImageDimension; ++d) + { + const auto sz = interiorRegion.GetSize(d); + if (sz <= 2 * kernelHalfWidth) + { + // Image is too small to crop -- mean SSIM falls back to the entire + // output region. + interiorIsValid = false; + break; + } + } + if (interiorIsValid) + { + auto idx = interiorRegion.GetIndex(); + auto sz = interiorRegion.GetSize(); + for (unsigned int d = 0; d < ImageDimension; ++d) + { + idx[d] += static_cast(kernelHalfWidth); + sz[d] -= 2 * kernelHalfWidth; + } + interiorRegion.SetIndex(idx); + interiorRegion.SetSize(sz); + } + + // ---- Parallelized per-pixel SSIM combination ------------------------ + // + // Each thread accumulates a partial sum and pixel count over its + // sub-region. Atomic doubles aren't portable in C++17, so we serialize + // the small per-thread reductions through a mutex. + + std::mutex accumulatorMutex; + RealType accumulator{}; + SizeValueType accumulatorCount{}; + + TotalProgressReporter progress(this, outputRegion.GetNumberOfPixels()); + + this->GetMultiThreader()->template ParallelizeImageRegion( + outputRegion, + [&](const OutputImageRegionType & subRegion) { + // Use ImageRegionIteratorWithIndex for the output so we can call + // GetIndex() inside the inner loop without hitting the + // ITK_LEGACY_REMOVE deprecation on the index-less iterator's + // GetIndex() (#6034 CI fix). + using OutputIteratorType = ImageRegionIteratorWithIndex; + using RealConstIteratorType = ImageRegionConstIterator; + + OutputIteratorType outIt(output, subRegion); + RealConstIteratorType gxIt(gx, subRegion); + RealConstIteratorType gyIt(gy, subRegion); + RealConstIteratorType gxxIt(gxx, subRegion); + RealConstIteratorType gyyIt(gyy, subRegion); + RealConstIteratorType gxyIt(gxy, subRegion); + + // Per-thread accumulators for the *interior* portion of this region. + const auto subInterior = [&]() { + OutputImageRegionType r = subRegion; + if (interiorIsValid) + { + if (!r.Crop(interiorRegion)) + { + r.SetSize(SizeType{}); // empty + } + } + return r; + }(); + + RealType localSum{}; + SizeValueType localCount{}; + + for (outIt.GoToBegin(), + gxIt.GoToBegin(), + gyIt.GoToBegin(), + gxxIt.GoToBegin(), + gyyIt.GoToBegin(), + gxyIt.GoToBegin(); + !outIt.IsAtEnd(); + ++outIt, ++gxIt, ++gyIt, ++gxxIt, ++gyyIt, ++gxyIt) + { + const RealType mu_x = gxIt.Get(); + const RealType mu_y = gyIt.Get(); + const RealType mu_xx = gxxIt.Get(); + const RealType mu_yy = gyyIt.Get(); + const RealType mu_xy = gxyIt.Get(); + + const RealType var_x = mu_xx - mu_x * mu_x; + const RealType var_y = mu_yy - mu_y * mu_y; + const RealType cov_xy = mu_xy - mu_x * mu_y; + + // Numerical floor: floating-point round-off can produce a tiny + // negative variance for nearly-constant regions. + const RealType var_x_clipped = std::max(var_x, RealType{}); + const RealType var_y_clipped = std::max(var_y, RealType{}); + + RealType ssim; + if (useSimplifiedFormula) + { + const RealType numerator = (RealType{ 2 } * mu_x * mu_y + C1) * (RealType{ 2 } * cov_xy + C2); + const RealType denominator = (mu_x * mu_x + mu_y * mu_y + C1) * (var_x_clipped + var_y_clipped + C2); + ssim = numerator / denominator; + } + else + { + const RealType l_num = RealType{ 2 } * mu_x * mu_y + C1; + const RealType l_den = mu_x * mu_x + mu_y * mu_y + C1; + const RealType l = l_num / l_den; + + const RealType sigma_x = std::sqrt(var_x_clipped); + const RealType sigma_y = std::sqrt(var_y_clipped); + + const RealType c_num = RealType{ 2 } * sigma_x * sigma_y + C2; + const RealType c_den = var_x_clipped + var_y_clipped + C2; + const RealType c = c_num / c_den; + + const RealType s_num = cov_xy + C3; + const RealType s_den = sigma_x * sigma_y + C3; + const RealType s = s_num / s_den; + + ssim = std::pow(l, alpha) * std::pow(c, beta) * std::pow(s, gamma); + } + + outIt.Set(static_cast(ssim)); + + // Only accumulate over the interior region for the mean. + if (subInterior.GetNumberOfPixels() > 0 && subInterior.IsInside(outIt.GetIndex())) + { + localSum += ssim; + ++localCount; + } + } + + progress.Completed(subRegion.GetNumberOfPixels()); + + { + const std::lock_guard lock(accumulatorMutex); + accumulator += localSum; + accumulatorCount += localCount; + } + }, + this); + + if (accumulatorCount > 0) + { + m_MeanSSIM = static_cast(accumulator) / static_cast(accumulatorCount); + } + else + { + m_MeanSSIM = 0.0; + } +} + +template +void +StructuralSimilarityImageFilter::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + os << indent << "GaussianSigma: " << m_GaussianSigma << std::endl; + os << indent << "MaximumKernelWidth: " << m_MaximumKernelWidth << std::endl; + os << indent << "K1: " << m_K1 << std::endl; + os << indent << "K2: " << m_K2 << std::endl; + os << indent << "DynamicRange: " << m_DynamicRange << std::endl; + os << indent << "LuminanceExponent: " << m_LuminanceExponent << std::endl; + os << indent << "ContrastExponent: " << m_ContrastExponent << std::endl; + os << indent << "StructureExponent: " << m_StructureExponent << std::endl; + os << indent << "ScaleWeights: " << m_ScaleWeights << std::endl; + os << indent << "MeanSSIM: " << m_MeanSSIM << std::endl; +} +} // end namespace itk +#endif diff --git a/Modules/Beta/StructuralSimilarity/itk-module.cmake b/Modules/Beta/StructuralSimilarity/itk-module.cmake new file mode 100644 index 00000000000..1dc875d6f99 --- /dev/null +++ b/Modules/Beta/StructuralSimilarity/itk-module.cmake @@ -0,0 +1,22 @@ +set( + DOCUMENTATION + "This module provides the StructuralSimilarityImageFilter, which +computes the Structural Similarity Index (SSIM) between two images. SSIM +is a perceptual image-quality metric that compares local luminance, +contrast, and structure within a sliding Gaussian window. See Wang et +al., \"Image quality assessment: From error visibility to structural +similarity\" (IEEE TIP, 2004)." +) + +itk_module( + StructuralSimilarity + DEPENDS + ITKImageIntensity + COMPILE_DEPENDS + ITKImageFilterBase + ITKSmoothing + TEST_DEPENDS + ITKTestKernel + ITKGoogleTest + DESCRIPTION "${DOCUMENTATION}" +) diff --git a/Modules/Beta/StructuralSimilarity/test/CMakeLists.txt b/Modules/Beta/StructuralSimilarity/test/CMakeLists.txt new file mode 100644 index 00000000000..cbe02ffc226 --- /dev/null +++ b/Modules/Beta/StructuralSimilarity/test/CMakeLists.txt @@ -0,0 +1,5 @@ +itk_module_test() + +set(StructuralSimilarityGTests itkStructuralSimilarityImageFilterGTest.cxx) + +creategoogletestdriver(StructuralSimilarity "${StructuralSimilarity-Test_LIBRARIES}" "${StructuralSimilarityGTests}") diff --git a/Modules/Beta/StructuralSimilarity/test/itkStructuralSimilarityImageFilterGTest.cxx b/Modules/Beta/StructuralSimilarity/test/itkStructuralSimilarityImageFilterGTest.cxx new file mode 100644 index 00000000000..89344305e2b --- /dev/null +++ b/Modules/Beta/StructuralSimilarity/test/itkStructuralSimilarityImageFilterGTest.cxx @@ -0,0 +1,674 @@ +/*========================================================================= + * + * 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 "itkImageRegionIterator.h" +#include "itkImageRegionIteratorWithIndex.h" +#include "itkMath.h" +#include "itkStructuralSimilarityImageFilter.h" + +#include +#include +#include + +namespace +{ +// +// All tests below pin tolerances against three classes of expected values: +// +// 1. Mathematical identities (identical images, constant inputs, symmetry). +// These do not depend on the Gaussian kernel implementation and use a +// tight tolerance (1e-6). +// +// 2. Closed-form formulas evaluated analytically over constant inputs. +// The luminance term reduces to (2*mu_x*mu_y + C1)/(mu_x^2 + mu_y^2 + C1) +// because variances and covariance vanish. Tight tolerance (1e-6). +// +// 3. Reference values cross-checked against scikit-image +// (gaussian_weights=True, sigma=1.5, use_sample_covariance=False) with +// a tolerance of 5e-3 to absorb minor discretization differences +// between ITK's GaussianOperator and scipy's sampled Gaussian. +// + +// ---- Common types --------------------------------------------------------- + +using PixelType = double; +constexpr unsigned int Dimension2D = 2; +using ImageType = itk::Image; +using FilterType = itk::StructuralSimilarityImageFilter; + +// ---- Image factories ------------------------------------------------------ + +ImageType::Pointer +MakeConstantImage(PixelType value, unsigned int size = 64) +{ + 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 +MakeGradientImage(unsigned int size = 64, PixelType scale = 255.0) +{ + auto image = ImageType::New(); + ImageType::SizeType imageSize; + imageSize.Fill(size); + const ImageType::RegionType region(imageSize); + image->SetRegions(region); + image->Allocate(); + // (x + y) normalized to [0, scale] + const PixelType norm = scale / static_cast(2 * (size - 1)); + itk::ImageRegionIteratorWithIndex it(image, region); + for (it.GoToBegin(); !it.IsAtEnd(); ++it) + { + const auto idx = it.GetIndex(); + it.Set(norm * static_cast(idx[0] + idx[1])); + } + return image; +} + +ImageType::Pointer +MakeRandomImage(unsigned int size = 64, unsigned int seed = 42, PixelType lo = 0.0, PixelType hi = 255.0) +{ + auto image = ImageType::New(); + ImageType::SizeType imageSize; + imageSize.Fill(size); + const ImageType::RegionType region(imageSize); + image->SetRegions(region); + image->Allocate(); + std::mt19937 gen(seed); + std::uniform_real_distribution dist(lo, hi); + itk::ImageRegionIterator it(image, region); + for (it.GoToBegin(); !it.IsAtEnd(); ++it) + { + it.Set(dist(gen)); + } + return image; +} + +ImageType::Pointer +ScaledCopy(const ImageType * input, PixelType scale, PixelType bias = 0.0) +{ + auto image = ImageType::New(); + const auto region = input->GetLargestPossibleRegion(); + image->SetRegions(region); + image->Allocate(); + itk::ImageRegionConstIterator in(input, region); + itk::ImageRegionIterator out(image, region); + for (in.GoToBegin(), out.GoToBegin(); !in.IsAtEnd(); ++in, ++out) + { + out.Set(scale * in.Get() + bias); + } + return image; +} + +ImageType::Pointer +NoisyCopy(const ImageType * input, PixelType noiseSigma, unsigned int seed) +{ + auto image = ImageType::New(); + const auto region = input->GetLargestPossibleRegion(); + image->SetRegions(region); + image->Allocate(); + std::mt19937 gen(seed); + std::normal_distribution noise(0.0, noiseSigma); + itk::ImageRegionConstIterator in(input, region); + itk::ImageRegionIterator out(image, region); + for (in.GoToBegin(), out.GoToBegin(); !in.IsAtEnd(); ++in, ++out) + { + out.Set(in.Get() + noise(gen)); + } + return image; +} + +double +ComputeFiltered(const ImageType * a, const ImageType * b) +{ + auto filter = FilterType::New(); + filter->SetInput1(a); + filter->SetInput2(b); + filter->Update(); + return filter->GetMeanSSIM(); +} + +} // namespace + + +// =========================================================================== +// Object-level macros and basic API +// =========================================================================== + +TEST(StructuralSimilarityImageFilter, BasicObjectMethods) +{ + auto filter = FilterType::New(); + ITK_GTEST_EXERCISE_BASIC_OBJECT_METHODS(filter, StructuralSimilarityImageFilter, ImageToImageFilter); +} + +TEST(StructuralSimilarityImageFilter, DefaultParameters) +{ + auto filter = FilterType::New(); + EXPECT_DOUBLE_EQ(filter->GetGaussianSigma(), 1.5); + EXPECT_EQ(filter->GetMaximumKernelWidth(), 11u); + EXPECT_DOUBLE_EQ(filter->GetK1(), 0.01); + EXPECT_DOUBLE_EQ(filter->GetK2(), 0.03); + EXPECT_DOUBLE_EQ(filter->GetLuminanceExponent(), 1.0); + EXPECT_DOUBLE_EQ(filter->GetContrastExponent(), 1.0); + EXPECT_DOUBLE_EQ(filter->GetStructureExponent(), 1.0); + EXPECT_EQ(filter->GetScaleWeights().GetSize(), 1u); + EXPECT_DOUBLE_EQ(filter->GetScaleWeights()[0], 1.0); + // For PixelType=double the default dynamic range is 1.0. + EXPECT_DOUBLE_EQ(filter->GetDynamicRange(), 1.0); +} + +TEST(StructuralSimilarityImageFilter, SetGetParameters) +{ + auto filter = FilterType::New(); + filter->SetGaussianSigma(2.0); + EXPECT_DOUBLE_EQ(filter->GetGaussianSigma(), 2.0); + filter->SetMaximumKernelWidth(15); + EXPECT_EQ(filter->GetMaximumKernelWidth(), 15u); + filter->SetK1(0.02); + EXPECT_DOUBLE_EQ(filter->GetK1(), 0.02); + filter->SetK2(0.04); + EXPECT_DOUBLE_EQ(filter->GetK2(), 0.04); + filter->SetDynamicRange(255.0); + EXPECT_DOUBLE_EQ(filter->GetDynamicRange(), 255.0); + filter->SetLuminanceExponent(2.0); + EXPECT_DOUBLE_EQ(filter->GetLuminanceExponent(), 2.0); + filter->SetContrastExponent(0.5); + EXPECT_DOUBLE_EQ(filter->GetContrastExponent(), 0.5); + filter->SetStructureExponent(0.7); + EXPECT_DOUBLE_EQ(filter->GetStructureExponent(), 0.7); + + FilterType::ScaleWeightsType weights(3); + weights[0] = 0.4; + weights[1] = 0.4; + weights[2] = 0.2; + filter->SetScaleWeights(weights); + EXPECT_EQ(filter->GetScaleWeights().GetSize(), 3u); + EXPECT_DOUBLE_EQ(filter->GetScaleWeights()[0], 0.4); + EXPECT_DOUBLE_EQ(filter->GetScaleWeights()[2], 0.2); +} + + +// =========================================================================== +// Mathematical identities (kernel-implementation independent) +// =========================================================================== + +TEST(StructuralSimilarityImageFilter, IdenticalConstantImagesYieldOne) +{ + auto image = MakeConstantImage(100.0, 64); + auto filter = FilterType::New(); + filter->SetInput1(image); + filter->SetInput2(image); + filter->SetDynamicRange(255.0); + filter->Update(); + EXPECT_NEAR(filter->GetMeanSSIM(), 1.0, 1e-9); +} + +TEST(StructuralSimilarityImageFilter, IdenticalRandomImagesYieldOne) +{ + auto image = MakeRandomImage(64, 42); + auto filter = FilterType::New(); + filter->SetInput1(image); + filter->SetInput2(image); + filter->SetDynamicRange(255.0); + filter->Update(); + EXPECT_NEAR(filter->GetMeanSSIM(), 1.0, 1e-9); +} + +TEST(StructuralSimilarityImageFilter, IdenticalGradientImagesYieldOne) +{ + auto image = MakeGradientImage(64); + auto filter = FilterType::New(); + filter->SetInput1(image); + filter->SetInput2(image); + filter->SetDynamicRange(255.0); + filter->Update(); + EXPECT_NEAR(filter->GetMeanSSIM(), 1.0, 1e-9); +} + +TEST(StructuralSimilarityImageFilter, SymmetryProperty) +{ + // SSIM is symmetric: SSIM(a,b) == SSIM(b,a) up to floating-point round-off. + auto a = MakeRandomImage(48, 11); + auto b = MakeRandomImage(48, 22); + auto filter1 = FilterType::New(); + filter1->SetInput1(a); + filter1->SetInput2(b); + filter1->SetDynamicRange(255.0); + filter1->Update(); + const double s_ab = filter1->GetMeanSSIM(); + + auto filter2 = FilterType::New(); + filter2->SetInput1(b); + filter2->SetInput2(a); + filter2->SetDynamicRange(255.0); + filter2->Update(); + const double s_ba = filter2->GetMeanSSIM(); + + EXPECT_NEAR(s_ab, s_ba, 1e-12); +} + + +// =========================================================================== +// Closed-form analytic checks for constant inputs +// =========================================================================== + +// SSIM(constant_a, constant_b) reduces to the luminance term only because +// variances and covariance vanish: +// +// SSIM = (2*a*b + C1) / (a^2 + b^2 + C1) +// +// (the contrast and structure terms each become 1). +static double +AnalyticConstantSSIM(double a, double b, double K1 = 0.01, double L = 255.0) +{ + const double C1 = (K1 * L) * (K1 * L); + return (2 * a * b + C1) / (a * a + b * b + C1); +} + +TEST(StructuralSimilarityImageFilter, TwoDifferentConstantImages_AnalyticMatch) +{ + // (100, 150) -> approx 0.92310 + auto a = MakeConstantImage(100.0); + auto b = MakeConstantImage(150.0); + auto filter = FilterType::New(); + filter->SetInput1(a); + filter->SetInput2(b); + filter->SetDynamicRange(255.0); + filter->Update(); + + const double expected = AnalyticConstantSSIM(100.0, 150.0); + EXPECT_NEAR(filter->GetMeanSSIM(), expected, 1e-9); + EXPECT_NEAR(expected, 0.9230923, 1e-6); // sanity-check the analytic helper +} + +TEST(StructuralSimilarityImageFilter, MaximallyDifferentConstants_AnalyticMatch) +{ + // (0, 255) with K1=0.01, L=255 gives the textbook 0.0000999900 + auto a = MakeConstantImage(0.0); + auto b = MakeConstantImage(255.0); + auto filter = FilterType::New(); + filter->SetInput1(a); + filter->SetInput2(b); + filter->SetDynamicRange(255.0); + filter->Update(); + + const double expected = AnalyticConstantSSIM(0.0, 255.0); + EXPECT_NEAR(filter->GetMeanSSIM(), expected, 1e-9); + EXPECT_NEAR(expected, 0.0000999900, 1e-9); +} + +TEST(StructuralSimilarityImageFilter, ConstantImages_PerPixelMapValueIsAnalytic) +{ + // Verify the per-pixel SSIM map is the analytic constant value everywhere. + auto a = MakeConstantImage(80.0, 32); + auto b = MakeConstantImage(120.0, 32); + auto filter = FilterType::New(); + filter->SetInput1(a); + filter->SetInput2(b); + filter->SetDynamicRange(255.0); + filter->Update(); + + const double expected = AnalyticConstantSSIM(80.0, 120.0); + + using OutputImageType = FilterType::OutputImageType; + itk::ImageRegionConstIterator it(filter->GetOutput(), filter->GetOutput()->GetBufferedRegion()); + unsigned int sampledCount = 0; + for (it.GoToBegin(); !it.IsAtEnd(); ++it) + { + EXPECT_NEAR(static_cast(it.Get()), expected, 1e-6); + ++sampledCount; + } + EXPECT_EQ(sampledCount, 32u * 32u); +} + + +// =========================================================================== +// Output geometry / data-flow correctness +// =========================================================================== + +TEST(StructuralSimilarityImageFilter, OutputMapHasSameSizeAsInput) +{ + constexpr unsigned int sz = 24; + auto a = MakeConstantImage(100.0, sz); + auto b = MakeConstantImage(100.0, sz); + auto filter = FilterType::New(); + filter->SetInput1(a); + filter->SetInput2(b); + filter->Update(); + EXPECT_EQ(filter->GetOutput()->GetLargestPossibleRegion().GetSize()[0], sz); + EXPECT_EQ(filter->GetOutput()->GetLargestPossibleRegion().GetSize()[1], sz); +} + +TEST(StructuralSimilarityImageFilter, MismatchedInputSizes_Throws) +{ + auto a = MakeConstantImage(100.0, 32); + auto b = MakeConstantImage(100.0, 48); + auto filter = FilterType::New(); + filter->SetInput1(a); + filter->SetInput2(b); + EXPECT_THROW(filter->Update(), itk::ExceptionObject); +} + +TEST(StructuralSimilarityImageFilter, MissingInputs_Throws) +{ + auto a = MakeConstantImage(100.0, 32); + auto filter = FilterType::New(); + filter->SetInput1(a); + EXPECT_THROW(filter->Update(), itk::ExceptionObject); +} + +TEST(StructuralSimilarityImageFilter, NonPositiveSigma_Throws) +{ + auto a = MakeConstantImage(100.0, 32); + auto b = MakeConstantImage(100.0, 32); + auto filter = FilterType::New(); + filter->SetInput1(a); + filter->SetInput2(b); + filter->SetGaussianSigma(0.0); + EXPECT_THROW(filter->Update(), itk::ExceptionObject); +} + +TEST(StructuralSimilarityImageFilter, NonPositiveDynamicRange_Throws) +{ + auto a = MakeConstantImage(100.0, 32); + auto b = MakeConstantImage(100.0, 32); + auto filter = FilterType::New(); + filter->SetInput1(a); + filter->SetInput2(b); + filter->SetDynamicRange(0.0); + EXPECT_THROW(filter->Update(), itk::ExceptionObject); +} + +TEST(StructuralSimilarityImageFilter, MultiScaleScaleWeights_NotYetImplemented_Throws) +{ + auto a = MakeConstantImage(100.0, 32); + auto b = MakeConstantImage(100.0, 32); + auto filter = FilterType::New(); + filter->SetInput1(a); + filter->SetInput2(b); + FilterType::ScaleWeightsType weights(5); + weights.Fill(0.2); + filter->SetScaleWeights(weights); + EXPECT_THROW(filter->Update(), itk::ExceptionObject); +} + +TEST(StructuralSimilarityImageFilter, EmptyScaleWeights_Throws) +{ + auto a = MakeConstantImage(100.0, 32); + auto b = MakeConstantImage(100.0, 32); + auto filter = FilterType::New(); + filter->SetInput1(a); + filter->SetInput2(b); + FilterType::ScaleWeightsType empty; + filter->SetScaleWeights(empty); + EXPECT_THROW(filter->Update(), itk::ExceptionObject); +} + + +// =========================================================================== +// Range / monotonicity / qualitative properties +// =========================================================================== + +TEST(StructuralSimilarityImageFilter, RandomPair_SSIMInValidRange) +{ + auto a = MakeRandomImage(64, 1); + auto b = MakeRandomImage(64, 2); + auto filter = FilterType::New(); + filter->SetInput1(a); + filter->SetInput2(b); + filter->SetDynamicRange(255.0); + filter->Update(); + const double s = filter->GetMeanSSIM(); + EXPECT_GE(s, -1.0); + EXPECT_LE(s, 1.0); +} + +TEST(StructuralSimilarityImageFilter, MoreNoiseLowersSSIM) +{ + auto base = MakeRandomImage(64, 5); + auto a_low = NoisyCopy(base, 2.0, 17); + auto a_mid = NoisyCopy(base, 8.0, 17); + auto a_hi = NoisyCopy(base, 24.0, 17); + + const double s_low = ComputeFiltered(base, a_low); + const double s_mid = ComputeFiltered(base, a_mid); + const double s_hi = ComputeFiltered(base, a_hi); + + EXPECT_GT(s_low, s_mid); + EXPECT_GT(s_mid, s_hi); + EXPECT_LT(s_low, 1.0); +} + +TEST(StructuralSimilarityImageFilter, NegatedImage_StronglyAntiCorrelated) +{ + auto base = MakeRandomImage(64, 99); + auto neg = ScaledCopy(base, -1.0, 255.0); // 255 - x + const double s = ComputeFiltered(base, neg); + // Cross-checked against scikit-image: typical value around -0.97. + EXPECT_LT(s, -0.5); + EXPECT_GE(s, -1.0); +} + + +// =========================================================================== +// Cross-check against scikit-image reference values (loose tolerance) +// =========================================================================== + +TEST(StructuralSimilarityImageFilter, GradientShiftedByConstant_SkimageReference) +{ + // Cross-check: gradient image shifted by +30 luminance. + // skimage.metrics.structural_similarity reference value: 0.9676912545 + auto a = MakeGradientImage(64); + auto b = ScaledCopy(a, 1.0, 30.0); + auto filter = FilterType::New(); + filter->SetInput1(a); + filter->SetInput2(b); + filter->SetDynamicRange(255.0); + filter->Update(); + // Loose tolerance to absorb Gaussian-kernel-discretization differences + // between ITK's GaussianOperator and scipy's sampled Gaussian. + EXPECT_NEAR(filter->GetMeanSSIM(), 0.9676912545, 5e-3); +} + +TEST(StructuralSimilarityImageFilter, GradientHalfContrast_SkimageReference) +{ + // Gradient * 0.5 (contrast change). + // skimage reference: 0.7550069937 + auto a = MakeGradientImage(64); + auto b = ScaledCopy(a, 0.5, 0.0); + auto filter = FilterType::New(); + filter->SetInput1(a); + filter->SetInput2(b); + filter->SetDynamicRange(255.0); + filter->Update(); + EXPECT_NEAR(filter->GetMeanSSIM(), 0.7550069937, 5e-3); +} + + +// =========================================================================== +// Code-path equivalence: simplified vs general formula +// =========================================================================== + +TEST(StructuralSimilarityImageFilter, SimplifiedAndGeneralFormulaAgreeWhenExponentsAreOne) +{ + // With alpha=beta=gamma=1 the filter takes the simplified fast path. We + // can force the general (std::pow) path by setting an exponent to a value + // *almost equal* to 1, which is treated as different and disables the + // FloatAlmostEqual fast-path check. In practice setting all exponents to + // exactly 1.0 (default) and another filter with values that differ by + // 1e-12 should give numerically equivalent answers. + auto a = MakeRandomImage(48, 33); + auto b = NoisyCopy(a, 5.0, 44); + + auto fastFilter = FilterType::New(); + fastFilter->SetInput1(a); + fastFilter->SetInput2(b); + fastFilter->SetDynamicRange(255.0); + fastFilter->Update(); + const double fastResult = fastFilter->GetMeanSSIM(); + + auto generalFilter = FilterType::New(); + generalFilter->SetInput1(a); + generalFilter->SetInput2(b); + generalFilter->SetDynamicRange(255.0); + // Tiny perturbation forces the general code path. + generalFilter->SetLuminanceExponent(1.0 + 1e-9); + generalFilter->SetContrastExponent(1.0); + generalFilter->SetStructureExponent(1.0); + generalFilter->Update(); + const double generalResult = generalFilter->GetMeanSSIM(); + + EXPECT_NEAR(fastResult, generalResult, 1e-6); +} + + +// =========================================================================== +// Multi-dimensional support (3D and 4D) +// =========================================================================== + +TEST(StructuralSimilarityImageFilter, ThreeDimensional_IdenticalRandomYieldsOne) +{ + using Image3DType = itk::Image; + using Filter3DType = itk::StructuralSimilarityImageFilter; + + auto image = Image3DType::New(); + Image3DType::SizeType size; + size.Fill(16); + image->SetRegions(Image3DType::RegionType(size)); + image->Allocate(); + + std::mt19937 gen(11); + std::uniform_real_distribution dist(0.0, 255.0); + itk::ImageRegionIterator it(image, image->GetLargestPossibleRegion()); + for (it.GoToBegin(); !it.IsAtEnd(); ++it) + { + it.Set(dist(gen)); + } + + auto filter = Filter3DType::New(); + filter->SetInput1(image); + filter->SetInput2(image); + filter->SetDynamicRange(255.0); + filter->Update(); + EXPECT_NEAR(filter->GetMeanSSIM(), 1.0, 1e-9); +} + +TEST(StructuralSimilarityImageFilter, ThreeDimensional_ConstantInputs_AnalyticMatch) +{ + using Image3DType = itk::Image; + using Filter3DType = itk::StructuralSimilarityImageFilter; + + auto a = Image3DType::New(); + auto b = Image3DType::New(); + Image3DType::SizeType size; + size.Fill(16); + a->SetRegions(Image3DType::RegionType(size)); + a->Allocate(); + a->FillBuffer(80.0); + b->SetRegions(Image3DType::RegionType(size)); + b->Allocate(); + b->FillBuffer(120.0); + + auto filter = Filter3DType::New(); + filter->SetInput1(a); + filter->SetInput2(b); + filter->SetDynamicRange(255.0); + filter->Update(); + EXPECT_NEAR(filter->GetMeanSSIM(), AnalyticConstantSSIM(80.0, 120.0), 1e-9); +} + +TEST(StructuralSimilarityImageFilter, FourDimensional_ConstantInputs_AnalyticMatch) +{ + using Image4DType = itk::Image; + using Filter4DType = itk::StructuralSimilarityImageFilter; + + auto a = Image4DType::New(); + auto b = Image4DType::New(); + Image4DType::SizeType size; + size.Fill(8); + a->SetRegions(Image4DType::RegionType(size)); + a->Allocate(); + a->FillBuffer(60.0); + b->SetRegions(Image4DType::RegionType(size)); + b->Allocate(); + b->FillBuffer(80.0); + + auto filter = Filter4DType::New(); + filter->SetInput1(a); + filter->SetInput2(b); + filter->SetDynamicRange(255.0); + filter->Update(); + EXPECT_NEAR(filter->GetMeanSSIM(), AnalyticConstantSSIM(60.0, 80.0), 1e-9); +} + + +// =========================================================================== +// Pixel-type variants (integer types use NumericTraits-based default L) +// =========================================================================== + +TEST(StructuralSimilarityImageFilter, UnsignedCharPixelType_DefaultDynamicRangeIs255) +{ + using UCharImage = itk::Image; + using UCharFilter = itk::StructuralSimilarityImageFilter; + auto filter = UCharFilter::New(); + EXPECT_DOUBLE_EQ(filter->GetDynamicRange(), 255.0); +} + +TEST(StructuralSimilarityImageFilter, UnsignedShortPixelType_DefaultDynamicRangeIs65535) +{ + using UShortImage = itk::Image; + using UShortFilter = itk::StructuralSimilarityImageFilter; + auto filter = UShortFilter::New(); + EXPECT_DOUBLE_EQ(filter->GetDynamicRange(), 65535.0); +} + +TEST(StructuralSimilarityImageFilter, FloatPixelType_DefaultDynamicRangeIsOne) +{ + using FloatImage = itk::Image; + using FloatFilter = itk::StructuralSimilarityImageFilter; + auto filter = FloatFilter::New(); + EXPECT_DOUBLE_EQ(filter->GetDynamicRange(), 1.0); +} + +TEST(StructuralSimilarityImageFilter, UnsignedCharPixelType_IdenticalYieldsOne) +{ + using UCharImage = itk::Image; + using UCharFilter = itk::StructuralSimilarityImageFilter; + + auto image = UCharImage::New(); + UCharImage::SizeType size; + size.Fill(32); + image->SetRegions(UCharImage::RegionType(size)); + image->Allocate(); + image->FillBuffer(static_cast(128)); + + auto filter = UCharFilter::New(); + filter->SetInput1(image); + filter->SetInput2(image); + filter->Update(); + EXPECT_NEAR(filter->GetMeanSSIM(), 1.0, 1e-9); +} diff --git a/Modules/Beta/StructuralSimilarity/wrapping/CMakeLists.txt b/Modules/Beta/StructuralSimilarity/wrapping/CMakeLists.txt new file mode 100644 index 00000000000..7f0787016d3 --- /dev/null +++ b/Modules/Beta/StructuralSimilarity/wrapping/CMakeLists.txt @@ -0,0 +1,2 @@ +itk_wrap_module(StructuralSimilarity) +itk_auto_load_and_end_wrap_submodules() diff --git a/Modules/Beta/StructuralSimilarity/wrapping/itkStructuralSimilarityImageFilter.wrap b/Modules/Beta/StructuralSimilarity/wrapping/itkStructuralSimilarityImageFilter.wrap new file mode 100644 index 00000000000..cba19b2c12c --- /dev/null +++ b/Modules/Beta/StructuralSimilarity/wrapping/itkStructuralSimilarityImageFilter.wrap @@ -0,0 +1,3 @@ +itk_wrap_class("itk::StructuralSimilarityImageFilter" POINTER) + itk_wrap_image_filter("${WRAP_ITK_REAL}" 2) +itk_end_wrap_class()