BUG: Guard division-by-zero NaN and UB-cast sites from issue 6575

Modules/Filtering/AnisotropicDiffusionLBR/include/itkStructureTensorImageFilter.hxx

@@ -150,7 +150,8 @@ StructureTensorImageFilter<TImage, TTensorImage>::GenerateData()
 
   traceFilter->Update();
   maximumCalculator->ComputeMaximum();
-  m_PostRescaling = 1. / maximumCalculator->GetMaximum();
+  const auto maximum = maximumCalculator->GetMaximum();
+  m_PostRescaling = (maximum > 0.) ? 1. / maximum : 1.;
   scaleFilter->GetFunctor().scaling = m_PostRescaling;
   scaleFilter->Update();
   this->GraftOutput(scaleFilter->GetOutput());

Modules/Filtering/AnisotropicDiffusionLBR/test/CMakeLists.txt

@@ -7,6 +7,7 @@ createtestdriver(AnisotropicDiffusionLBR  "${AnisotropicDiffusionLBR-Test_LIBRAR
 set(
   AnisotropicDiffusionLBRGTests
   itkLinearAnisotropicDiffusionLBRImageFilterGTest.cxx
+  itkStructureTensorImageFilterGTest.cxx
 )
 
 creategoogletestdriver(AnisotropicDiffusionLBR "${AnisotropicDiffusionLBR-Test_LIBRARIES}" "${AnisotropicDiffusionLBRGTests}")

Modules/Filtering/AnisotropicDiffusionLBR/test/itkStructureTensorImageFilterGTest.cxx

@@ -0,0 +1,51 @@
+/*=========================================================================
+ *
+ *  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 "itkStructureTensorImageFilter.h"
+#include "itkImageRegionConstIterator.h"
+#include "itkGTest.h"
+
+// A constant image has zero gradients, so the tensor-trace maximum is zero; the
+// adimensionizing rescale must not divide by it and poison the tensors with NaN
+// (issue #6575, B23).
+TEST(StructureTensorImageFilter, ConstantImageYieldsFiniteZeroTensors)
+{
+  using ImageType = itk::Image<float, 2>;
+  auto image = ImageType::New();
+  image->SetRegions(ImageType::RegionType{ ImageType::IndexType{}, ImageType::SizeType::Filled(16) });
+  image->Allocate();
+  image->FillBuffer(0.0F);
+
+  using FilterType = itk::StructureTensorImageFilter<ImageType>;
+  auto filter = FilterType::New();
+  filter->SetInput(image);
+  filter->SetRescaleForUnitMaximumTrace(true);
+  filter->Update();
+
+  using TensorImageType = FilterType::TensorImageType;
+  itk::ImageRegionConstIterator<TensorImageType> it(filter->GetOutput(), filter->GetOutput()->GetBufferedRegion());
+  for (it.GoToBegin(); !it.IsAtEnd(); ++it)
+  {
+    const auto & tensor = it.Get();
+    for (unsigned int c = 0; c < tensor.Size(); ++c)
+    {
+      ASSERT_TRUE(std::isfinite(tensor[c]));
+      ASSERT_NEAR(tensor[c], 0.0F, 1e-20);
+    }
+  }
+}

Modules/Filtering/BiasCorrection/include/itkN4BiasFieldCorrectionImageFilter.hxx

@@ -38,6 +38,7 @@ ITK_GCC_SUPPRESS_Wfloat_equal
 #include "vnl/vnl_complex_traits.h"
 #include "complex"
 #include "itkPrintHelper.h"
+#include <algorithm> // For clamp and copy.
 ITK_GCC_PRAGMA_POP
 
 namespace itk
@@ -285,12 +286,21 @@ N4BiasFieldCorrectionImageFilter<TInputImage, TMaskImage, TOutputImage>::Sharpen
       {
         binMaximum = pixel;
       }
-      else if (pixel < binMinimum)
+      if (pixel < binMinimum)
       {
         binMinimum = pixel;
       }
     }
   }
+  if (binMaximum <= binMinimum)
+  {
+    // Degenerate intensity range: sharpening is an identity mapping.
+    sharpenedImage->FillBuffer(0);
+    const ImageBufferRange identityImageBufferRange{ *sharpenedImage };
+    std::copy(
+      unsharpenedImageBufferRange.cbegin(), unsharpenedImageBufferRange.cend(), identityImageBufferRange.begin());
+    return;
+  }
   const RealType histogramSlope = (binMaximum - binMinimum) / static_cast<RealType>(this->m_NumberOfHistogramBins - 1);
 
   // Create the intensity profile (within the masked region, if applicable)

Modules/Filtering/BiasCorrection/test/CMakeLists.txt

@@ -7,7 +7,11 @@ set(
 
 createtestdriver(ITKBiasCorrection "${ITKBiasCorrection-Test_LIBRARIES}" "${ITKBiasCorrectionTests}")
 
-set(ITKBiasCorrectionGTests itkCompositeValleyFunctionGTest.cxx)
+set(
+  ITKBiasCorrectionGTests
+  itkCompositeValleyFunctionGTest.cxx
+  itkN4BiasFieldCorrectionImageFilterGTest.cxx
+)
 
 creategoogletestdriver(ITKBiasCorrection "${ITKBiasCorrection-Test_LIBRARIES}" "${ITKBiasCorrectionGTests}")
 

Modules/Filtering/BiasCorrection/test/itkN4BiasFieldCorrectionImageFilterGTest.cxx

@@ -0,0 +1,45 @@
+/*=========================================================================
+ *
+ *  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 "itkN4BiasFieldCorrectionImageFilter.h"
+#include "itkImageRegionConstIterator.h"
+#include "itkGTest.h"
+
+// A constant image has a degenerate histogram range; sharpening must degrade to an
+// identity mapping instead of producing NaN bin indices (issue #6575, B2).
+TEST(N4BiasFieldCorrectionImageFilter, ConstantImageProducesFiniteOutput)
+{
+  using ImageType = itk::Image<float, 2>;
+  auto image = ImageType::New();
+  image->SetRegions(ImageType::RegionType{ ImageType::IndexType{}, ImageType::SizeType::Filled(32) });
+  image->Allocate();
+  image->FillBuffer(1.0F);
+
+  using FilterType = itk::N4BiasFieldCorrectionImageFilter<ImageType>;
+  auto filter = FilterType::New();
+  filter->SetInput(image);
+  filter->SetMaximumNumberOfIterations(FilterType::VariableSizeArrayType(1, 2));
+  filter->Update();
+
+  itk::ImageRegionConstIterator<ImageType> it(filter->GetOutput(), filter->GetOutput()->GetBufferedRegion());
+  for (it.GoToBegin(); !it.IsAtEnd(); ++it)
+  {
+    ASSERT_TRUE(std::isfinite(it.Get()));
+    ASSERT_NEAR(it.Get(), 1.0F, 1e-2) << " constant input must be (nearly) unchanged";
+  }
+}

Modules/Filtering/Deconvolution/include/itkWienerDeconvolutionImageFilter.h

@@ -170,11 +170,14 @@ class ITK_TEMPLATE_EXPORT WienerDeconvolutionFunctor
     // minus the power spectral density of the noise.
     TPixel Pf = std::norm(I);
 
-    TPixel denominator = std::norm(H) + (Pn / (Pf - Pn));
     TPixel value{};
-    if (itk::Math::Absolute(denominator) >= m_KernelZeroMagnitudeThreshold)
+    if (Pf != Pn)
     {
-      value = I * (std::conj(H) / denominator);
+      const TPixel denominator = std::norm(H) + (Pn / (Pf - Pn));
+      if (itk::Math::Absolute(denominator) >= m_KernelZeroMagnitudeThreshold)
+      {
+        value = I * (std::conj(H) / denominator);
+      }
     }
 
     return value;

Modules/Filtering/Deconvolution/test/CMakeLists.txt

@@ -15,6 +15,7 @@ createtestdriver(ITKDeconvolution "${ITKDeconvolution-Test_LIBRARIES}" "${ITKDec
 set(
   ITKDeconvolutionGTests
   itkProjectedIterativeDeconvolutionImageFilterGTest.cxx
+  itkWienerDeconvolutionImageFilterGTest.cxx
 )
 
 creategoogletestdriver(ITKDeconvolution "${ITKDeconvolution-Test_LIBRARIES}" "${ITKDeconvolutionGTests}")

Modules/Filtering/Deconvolution/test/itkWienerDeconvolutionImageFilterGTest.cxx

@@ -0,0 +1,41 @@
+/*=========================================================================
+ *
+ *  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 "itkWienerDeconvolutionImageFilter.h"
+#include "itkGTest.h"
+#include <complex>
+
+// When the input power at a frequency equals the noise power spectral density, the
+// regularization term divides by zero; the functor must return zero instead of NaN
+// (issue #6575, B4).
+TEST(WienerDeconvolutionFunctor, ZeroSignalToNoiseYieldsZero)
+{
+  using ComplexType = std::complex<double>;
+  itk::Functor::WienerDeconvolutionFunctor<ComplexType> functor;
+  functor.SetNoisePowerSpectralDensityConstant(4.0);
+  functor.SetKernelZeroMagnitudeThreshold(1.0e-4);
+
+  // |I|^2 == Pn triggers the degenerate branch.
+  const ComplexType value = functor(ComplexType(2.0, 0.0), ComplexType(1.0, 0.0));
+  EXPECT_EQ(value, ComplexType(0.0, 0.0));
+
+  // A non-degenerate frequency still deconvolves finitely.
+  const ComplexType regular = functor(ComplexType(3.0, 0.0), ComplexType(1.0, 0.0));
+  EXPECT_TRUE(std::isfinite(regular.real()) && std::isfinite(regular.imag()));
+  EXPECT_NE(regular, ComplexType(0.0, 0.0));
+}

Modules/Filtering/Denoising/include/itkPatchBasedDenoisingImageFilter.hxx

@@ -30,6 +30,7 @@
 #include "itkSpatialNeighborSubsampler.h"
 #include "itkMacro.h"
 #include "itkMath.h"
+#include <algorithm> // For clamp and copy.
 
 namespace itk
 {
@@ -2075,6 +2076,14 @@ PatchBasedDenoisingImageFilter<TInputImage, TOutputImage>::ThreadedComputeImageU
       // this->ApplyUpdate() will then copy the results to outputIt since we
       // have to wait until all threads have finished using outputIt to
       // compute the update
+      for (unsigned int pc = 0; pc < m_NumPixelComponents; ++pc)
+      {
+        // Clamp to the representable range: the float-to-integer conversion of an
+        // out-of-range update is undefined behavior.
+        const auto lowest = static_cast<RealValueType>(NumericTraits<PixelValueType>::NonpositiveMin());
+        const auto highest = static_cast<RealValueType>(NumericTraits<PixelValueType>::max());
+        this->SetComponent(result, pc, std::clamp(this->GetComponent(result, pc), lowest, highest));
+      }
       updateIt.Set(static_cast<PixelType>(result));
 
       ++updateIt;

Modules/Filtering/ImageCompare/include/itkSTAPLEImageFilter.hxx

@@ -73,6 +73,8 @@ STAPLEImageFilter<TInputImage, TOutputImage>::GenerateData()
   {
     last_p[i] = -10.0;
     last_q[i] = -10.0;
+    p[i] = 1.0;
+    q[i] = 1.0;
   }
 
   // Come up with an initial Wi which is simply the average of
@@ -158,8 +160,14 @@ STAPLEImageFilter<TInputImage, TOutputImage>::GenerateData()
         out.NextLine();
       }
 
-      p[i] = p_num / p_denom;
-      q[i] = q_num / q_denom;
+      if (p_denom != 0.0)
+      {
+        p[i] = p_num / p_denom;
+      }
+      if (q_denom != 0.0)
+      {
+        q[i] = q_num / q_denom;
+      }
     }
 
     // Now recreate W using the new p's and q's

Modules/Filtering/ImageCompare/test/CMakeLists.txt

@@ -13,6 +13,9 @@ set(
 
 createtestdriver(ITKImageCompare "${ITKImageCompare-Test_LIBRARIES}" "${ITKImageCompareTests}")
 
+set(ITKImageCompareGTests itkSTAPLEImageFilterGTest.cxx)
+creategoogletestdriver(ITKImageCompare "${ITKImageCompare-Test_LIBRARIES}" "${ITKImageCompareGTests}")
+
 itk_add_test(
   NAME itkAbsoluteValueDifferenceImageFilterTest
   COMMAND

Modules/Filtering/ImageCompare/test/itkSTAPLEImageFilterGTest.cxx

@@ -0,0 +1,56 @@
+/*=========================================================================
+ *
+ *  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 "itkSTAPLEImageFilter.h"
+#include "itkImageRegionConstIterator.h"
+#include "itkGTest.h"
+
+// All-background segmentations zero the weight image, so the sensitivity denominator is
+// zero; the estimates must stay finite instead of becoming NaN (issue #6575, B11).
+TEST(STAPLEImageFilter, AllBackgroundInputsStayFinite)
+{
+  using InputImageType = itk::Image<unsigned char, 2>;
+  using OutputImageType = itk::Image<double, 2>;
+
+  const auto makeImage = [] {
+    auto image = InputImageType::New();
+    image->SetRegions(InputImageType::RegionType{ InputImageType::IndexType{}, InputImageType::SizeType::Filled(16) });
+    image->Allocate();
+    image->FillBuffer(0);
+    return image;
+  };
+
+  using FilterType = itk::STAPLEImageFilter<InputImageType, OutputImageType>;
+  auto filter = FilterType::New();
+  filter->SetInput(0, makeImage());
+  filter->SetInput(1, makeImage());
+  filter->SetForegroundValue(1);
+  filter->SetMaximumIterations(5);
+  filter->Update();
+
+  for (unsigned int i = 0; i < 2; ++i)
+  {
+    EXPECT_TRUE(std::isfinite(filter->GetSensitivity(i))) << " input " << i;
+    EXPECT_TRUE(std::isfinite(filter->GetSpecificity(i))) << " input " << i;
+  }
+  itk::ImageRegionConstIterator<OutputImageType> it(filter->GetOutput(), filter->GetOutput()->GetBufferedRegion());
+  for (it.GoToBegin(); !it.IsAtEnd(); ++it)
+  {
+    ASSERT_TRUE(std::isfinite(it.Get()));
+  }
+}

Modules/Filtering/ImageStatistics/include/itkAdaptiveEqualizationHistogram.h

@@ -74,6 +74,11 @@ class AdaptiveEqualizationHistogram
     // AdaptiveHistogramEqualization compute kernel components with
     // float, but use double for accumulate and temporaries.
     const double iscale = static_cast<double>(m_Maximum) - m_Minimum;
+    if (iscale == 0.0)
+    {
+      // Constant image: equalization is an identity mapping.
+      return static_cast<TOutputPixel>(pixel);
+    }
 
     double         sum = 0.0;
     auto           itMap = m_Map.begin();

Modules/Filtering/ImageStatistics/include/itkLabelOverlapMeasuresImageFilter.hxx

@@ -385,17 +385,16 @@ LabelOverlapMeasuresImageFilter<TLabelImage>::GetFalsePositiveError(LabelType la
     return 0.0;
   }
 
-  RealType value;
-  if (Math::ExactlyEquals(mapIt->second.m_Source, 0.0))
+  RealType   value;
+  const auto nComplementIntersection = nVox - mapIt->second.m_Union; // TN
+  const auto denominator = mapIt->second.m_SourceComplement + nComplementIntersection;
+  if (denominator == 0)
   {
     value = NumericTraits<RealType>::max();
   }
   else
   {
-    auto nComplementIntersection = nVox - mapIt->second.m_Union; // TN
-
-    value = static_cast<RealType>(mapIt->second.m_SourceComplement) /
-            static_cast<RealType>(mapIt->second.m_SourceComplement + nComplementIntersection);
+    value = static_cast<RealType>(mapIt->second.m_SourceComplement) / static_cast<RealType>(denominator);
   }
 
   return value;

Modules/Filtering/ImageStatistics/test/CMakeLists.txt

@@ -467,6 +467,7 @@ itk_add_test(
 
 set(
   ITKImageStatisticsGTests
+  itkAdaptiveHistogramEqualizationImageFilterGTest.cxx
   itkLabelOverlapMeasuresImageFilterGTest.cxx
   itkMinimumMaximumImageFilterGTest.cxx
 )