ENH: Add StructuralSimilarityImageFilter for SSIM image quality

Implements the Structural Similarity Index Measure as an N-dimensional,
multi-threaded ITK filter in the ImageCompare module.  Closes #6030,
supersedes #6031.

== Algorithm and reference sources ==

The implementation follows the canonical SSIM formulation of Wang, Bovik,
Sheikh, and Simoncelli, "Image Quality Assessment: From Error Visibility
to Structural Similarity," IEEE Trans. Image Processing 13(4), 2004.
The reference materials consulted while building this filter were:

  - Wang et al. 2004 paper (full text):
    https://www.cns.nyu.edu/pub/eero/wang03-reprint.pdf
  - Wang et al. SSIM MATLAB reference (utlive/ssim/ssim_index.m):
    https://github.com/utlive/ssim/blob/main/ssim_index.m
  - scikit-image v0.25 structural_similarity implementation:
    https://github.com/scikit-image/scikit-image/blob/v0.25.0/skimage/metrics/_structural_similarity.py
  - Wang, Simoncelli, Bovik, "Multi-Scale Structural Similarity for Image
    Quality Assessment," Asilomar 2003 (MS-SSIM, future extension):
    https://www.cns.nyu.edu/pub/eero/wang03b.pdf
  - Wikipedia SSIM article (formulas, defaults):
    https://en.wikipedia.org/wiki/Structural_similarity_index_measure

For two images x and y the filter computes local statistics by convolving
with a discrete Gaussian (sigma=1.5, 11x11 default, matching Wang et al.
and the default of skimage.metrics.structural_similarity):

  mu_x  = G_sigma * x
  mu_y  = G_sigma * y
  var_x = G_sigma * (x*x) - mu_x^2
  var_y = G_sigma * (y*y) - mu_y^2
  cov   = G_sigma * (x*y) - mu_x * mu_y

The three SSIM components are
  l(x,y) = (2*mu_x*mu_y + C1) / (mu_x^2 + mu_y^2 + C1)
  c(x,y) = (2*sigma_x*sigma_y + C2) / (var_x + var_y + C2)
  s(x,y) = (cov + C3) / (sigma_x*sigma_y + C3)
with C1 = (K1*L)^2, C2 = (K2*L)^2, C3 = C2/2, K1=0.01, K2=0.03, L the
dynamic range.  The combined SSIM is l^alpha * c^beta * s^gamma.  When
alpha=beta=gamma=1 (the default), the filter takes the simplified product
fast path
  SSIM = (2*mu_x*mu_y + C1)*(2*cov + C2) /
         ((mu_x^2 + mu_y^2 + C1)*(var_x + var_y + C2))
which is what Wang et al.'s reference distributes and what skimage uses
by default.

== Filter architecture ==

The filter is structured as a composite ImageToImageFilter:

  1. Internal sub-pipeline of five DiscreteGaussianImageFilter passes
     (mu_x, mu_y, mu_xx, mu_yy, mu_xy) reusing ITK's well-optimized,
     multi-threaded smoothing.
  2. A parallelized per-pixel combination via
     MultiThreaderBase::ParallelizeImageRegion that reads from the five
     smoothed buffers and writes the SSIM map.
  3. The mean SSIM is accumulated only over the interior region (cropped
     by half the Gaussian kernel width), matching scikit-image and the
     MATLAB reference, since pixels within the kernel half-width use
     boundary-extended values inside the convolution and are less
     reliable.

Inputs:  two images of the same template type and identical region.
Outputs:
  - a per-pixel SSIM map (TOutputImage, default Image<float, D>)
  - GetMeanSSIM() returns the scalar after Update()

Configurable runtime parameters (covering all of issue #6030):
  - GaussianSigma                      (default 1.5)
  - MaximumKernelWidth                 (default 11)
  - K1, K2 stability constants         (defaults 0.01, 0.03)
  - DynamicRange L                     (NumericTraits-derived: 1.0 for
                                         float/double, 255 for uchar,
                                         65535 for ushort, ...)
  - LuminanceExponent  alpha           (default 1.0)
  - ContrastExponent   beta            (default 1.0)
  - StructureExponent  gamma           (default 1.0)
  - ScaleWeights array                 (default {1.0} -> single-scale).
                                       Multi-element arrays reserve API
                                       space for a future MS-SSIM
                                       extension and currently throw a
                                       not-yet-implemented exception in
                                       BeforeGenerate.

== Test strategy ==

The filter ships with 30 GoogleTest assertions in
itkStructuralSimilarityImageFilterGTest.cxx, organized into four classes
of expected values that decouple correctness checks from sensitivity to
the discrete-Gaussian implementation:

  Class 1 -- mathematical identities (kernel-independent, tolerance 1e-9)
    SSIM(x, x) = 1 exactly for any image (constant, random, gradient).
    SSIM is symmetric: SSIM(a,b) == SSIM(b,a).

  Class 2 -- closed-form analytic checks for constant inputs
    Constant inputs make all variances and the covariance vanish, so
    SSIM(constant_a, constant_b) = (2*a*b + C1) / (a^2 + b^2 + C1).
    Tested at (100, 150) -> 0.9230923 and the textbook (0, 255) ->
    0.0000999900.  Verified pixel-wise across the output map (every
    map element matches the closed form).

  Class 3 -- input validation (exception tests)
    Mismatched input sizes, missing inputs, non-positive sigma,
    non-positive dynamic range, empty ScaleWeights, multi-element
    ScaleWeights (MS-SSIM not yet implemented).

  Class 4 -- qualitative properties
    Result range bounded in [-1, 1].
    Monotonic decay of mean SSIM as additive Gaussian noise grows
    (sigma = 2 -> 8 -> 24).
    Strong anti-correlation for negated images (SSIM(x, 255-x) < -0.5).

  Class 5 -- cross-checks against scikit-image (loose tolerance 5e-3)
    Two reference values were computed offline against
    skimage.metrics.structural_similarity with gaussian_weights=True,
    sigma=1.5, use_sample_covariance=False, data_range=255, win_size=11
    (i.e. the canonical Wang configuration):
      gradient + 30 luminance shift -> 0.9676912545
      gradient * 0.5  contrast      -> 0.7550069937
    The 5e-3 tolerance absorbs minor discretization differences between
    ITK's GaussianOperator and scipy's sampled Gaussian (the two
    libraries do not produce bit-identical Gaussian kernels).

  Class 6 -- code-path equivalence
    The simplified-product fast path (alpha=beta=gamma=1) and the
    general l^alpha*c^beta*s^gamma path are exercised on the same
    inputs and required to agree to 1e-6.

  Class 7 -- multi-dimensional and pixel-type coverage
    3D and 4D variants of the identity and constant-input tests.
    Default DynamicRange is correct for unsigned char (255), unsigned
    short (65535), and float (1.0).

The reference values for Class 5 were generated with the following
script (commit history; not shipped):

  import numpy as np
  from skimage.metrics import structural_similarity as ssim
  def ref(x, y, L=255.0):
      return ssim(x, y, gaussian_weights=True, sigma=1.5,
                  use_sample_covariance=False, data_range=L,
                  win_size=11, K1=0.01, K2=0.03)

== Results ==

Local build with GCC 13.3 / Ninja / Release on Ubuntu 24.04:

  $ cmake --build build-ssim -j48 --target ITKImageCompareGTestDriver
  [4/4] Linking CXX executable bin/ITKImageCompareGTestDriver

  $ ./bin/ITKImageCompareGTestDriver \\
        --gtest_filter='StructuralSimilarityImageFilter.*'
  [==========] 30 tests from 1 test suite ran. (60 ms total)
  [  PASSED  ] 30 tests.

  $ ctest -R 'StructuralSimilarity' --output-on-failure
  100% tests passed, 0 tests failed out of 30
  Total Test time (real) = 0.30 sec

ITKImageCompareHeaderTest1 (the auto-generated header self-test) and the
existing ITKImageCompareTestDriver also link cleanly with the new module
dependencies.

pre-commit (gersemi, clang-format, kw-pre-commit, etc.) reports all
checks Passed on every touched file.

== Module dependency changes ==

ITKImageCompare/itk-module.cmake gains:
  - ITKSmoothing as a COMPILE_DEPENDS  (for DiscreteGaussianImageFilter)
  - ITKGoogleTest as a TEST_DEPENDS    (for the new GTest driver)

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: hjmjohnson

Modules/Filtering/ImageCompare/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 <type_traits>
+
+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 ITKImageCompare
+ */
+template <typename TInputImage, typename TOutputImage = Image<float, TInputImage::ImageDimension>>
+class ITK_TEMPLATE_EXPORT StructuralSimilarityImageFilter : public ImageToImageFilter<TInputImage, TOutputImage>
+{
+public:
+  ITK_DISALLOW_COPY_AND_MOVE(StructuralSimilarityImageFilter);
+
+  /** Standard class type aliases. */
+  using Self = StructuralSimilarityImageFilter;
+  using Superclass = ImageToImageFilter<TInputImage, TOutputImage>;
+  using Pointer = SmartPointer<Self>;
+  using ConstPointer = SmartPointer<const Self>;
+
+  /** 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<InputPixelType>::RealType;
+
+  /** Type used for the user-specified array of multi-scale weights. */
+  using ScaleWeightsType = Array<RealType>;
+
+  /** 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<InputPixelType>));
+
+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<InputPixelType>)
+    {
+      return 1.0;
+    }
+    else
+    {
+      return static_cast<double>(NumericTraits<InputPixelType>::max()) -
+             static_cast<double>(NumericTraits<InputPixelType>::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