diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/include/itkStructureTensorImageFilter.hxx b/Modules/Filtering/AnisotropicDiffusionLBR/include/itkStructureTensorImageFilter.hxx index b6d1cfaa77d..bf3e6389668 100644 --- a/Modules/Filtering/AnisotropicDiffusionLBR/include/itkStructureTensorImageFilter.hxx +++ b/Modules/Filtering/AnisotropicDiffusionLBR/include/itkStructureTensorImageFilter.hxx @@ -150,7 +150,8 @@ StructureTensorImageFilter::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()); diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/CMakeLists.txt b/Modules/Filtering/AnisotropicDiffusionLBR/test/CMakeLists.txt index 6be49446c8d..c47d032e291 100644 --- a/Modules/Filtering/AnisotropicDiffusionLBR/test/CMakeLists.txt +++ b/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}") diff --git a/Modules/Filtering/AnisotropicDiffusionLBR/test/itkStructureTensorImageFilterGTest.cxx b/Modules/Filtering/AnisotropicDiffusionLBR/test/itkStructureTensorImageFilterGTest.cxx new file mode 100644 index 00000000000..dac2ac7f429 --- /dev/null +++ b/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; + auto image = ImageType::New(); + image->SetRegions(ImageType::RegionType{ ImageType::IndexType{}, ImageType::SizeType::Filled(16) }); + image->Allocate(); + image->FillBuffer(0.0F); + + using FilterType = itk::StructureTensorImageFilter; + auto filter = FilterType::New(); + filter->SetInput(image); + filter->SetRescaleForUnitMaximumTrace(true); + filter->Update(); + + using TensorImageType = FilterType::TensorImageType; + itk::ImageRegionConstIterator 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); + } + } +} diff --git a/Modules/Filtering/BiasCorrection/include/itkN4BiasFieldCorrectionImageFilter.hxx b/Modules/Filtering/BiasCorrection/include/itkN4BiasFieldCorrectionImageFilter.hxx index aa4fd7d6ae5..70729bb01fa 100644 --- a/Modules/Filtering/BiasCorrection/include/itkN4BiasFieldCorrectionImageFilter.hxx +++ b/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 // For clamp and copy. ITK_GCC_PRAGMA_POP namespace itk @@ -285,12 +286,21 @@ N4BiasFieldCorrectionImageFilter::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(this->m_NumberOfHistogramBins - 1); // Create the intensity profile (within the masked region, if applicable) diff --git a/Modules/Filtering/BiasCorrection/test/CMakeLists.txt b/Modules/Filtering/BiasCorrection/test/CMakeLists.txt index 96d7d3a8a45..7b525cece3e 100644 --- a/Modules/Filtering/BiasCorrection/test/CMakeLists.txt +++ b/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}") diff --git a/Modules/Filtering/BiasCorrection/test/itkN4BiasFieldCorrectionImageFilterGTest.cxx b/Modules/Filtering/BiasCorrection/test/itkN4BiasFieldCorrectionImageFilterGTest.cxx new file mode 100644 index 00000000000..77df8e5f868 --- /dev/null +++ b/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; + auto image = ImageType::New(); + image->SetRegions(ImageType::RegionType{ ImageType::IndexType{}, ImageType::SizeType::Filled(32) }); + image->Allocate(); + image->FillBuffer(1.0F); + + using FilterType = itk::N4BiasFieldCorrectionImageFilter; + auto filter = FilterType::New(); + filter->SetInput(image); + filter->SetMaximumNumberOfIterations(FilterType::VariableSizeArrayType(1, 2)); + filter->Update(); + + itk::ImageRegionConstIterator 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"; + } +} diff --git a/Modules/Filtering/Deconvolution/include/itkWienerDeconvolutionImageFilter.h b/Modules/Filtering/Deconvolution/include/itkWienerDeconvolutionImageFilter.h index 0c992790222..9494d46183b 100644 --- a/Modules/Filtering/Deconvolution/include/itkWienerDeconvolutionImageFilter.h +++ b/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; diff --git a/Modules/Filtering/Deconvolution/test/CMakeLists.txt b/Modules/Filtering/Deconvolution/test/CMakeLists.txt index bb8e3ed6319..25aa381788a 100644 --- a/Modules/Filtering/Deconvolution/test/CMakeLists.txt +++ b/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}") diff --git a/Modules/Filtering/Deconvolution/test/itkWienerDeconvolutionImageFilterGTest.cxx b/Modules/Filtering/Deconvolution/test/itkWienerDeconvolutionImageFilterGTest.cxx new file mode 100644 index 00000000000..0337734cc6c --- /dev/null +++ b/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 + +// 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; + itk::Functor::WienerDeconvolutionFunctor 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)); +} diff --git a/Modules/Filtering/Denoising/include/itkPatchBasedDenoisingImageFilter.hxx b/Modules/Filtering/Denoising/include/itkPatchBasedDenoisingImageFilter.hxx index a929aa42d16..bb2b1ef4b67 100644 --- a/Modules/Filtering/Denoising/include/itkPatchBasedDenoisingImageFilter.hxx +++ b/Modules/Filtering/Denoising/include/itkPatchBasedDenoisingImageFilter.hxx @@ -30,6 +30,7 @@ #include "itkSpatialNeighborSubsampler.h" #include "itkMacro.h" #include "itkMath.h" +#include // For clamp and copy. namespace itk { @@ -2075,6 +2076,14 @@ PatchBasedDenoisingImageFilter::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(NumericTraits::NonpositiveMin()); + const auto highest = static_cast(NumericTraits::max()); + this->SetComponent(result, pc, std::clamp(this->GetComponent(result, pc), lowest, highest)); + } updateIt.Set(static_cast(result)); ++updateIt; diff --git a/Modules/Filtering/ImageCompare/include/itkSTAPLEImageFilter.hxx b/Modules/Filtering/ImageCompare/include/itkSTAPLEImageFilter.hxx index 844df742b33..90647ea6a65 100644 --- a/Modules/Filtering/ImageCompare/include/itkSTAPLEImageFilter.hxx +++ b/Modules/Filtering/ImageCompare/include/itkSTAPLEImageFilter.hxx @@ -73,6 +73,8 @@ STAPLEImageFilter::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::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 diff --git a/Modules/Filtering/ImageCompare/test/CMakeLists.txt b/Modules/Filtering/ImageCompare/test/CMakeLists.txt index 87260430bd1..a7dc12111f3 100644 --- a/Modules/Filtering/ImageCompare/test/CMakeLists.txt +++ b/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 diff --git a/Modules/Filtering/ImageCompare/test/itkSTAPLEImageFilterGTest.cxx b/Modules/Filtering/ImageCompare/test/itkSTAPLEImageFilterGTest.cxx new file mode 100644 index 00000000000..5d4e4fd7c58 --- /dev/null +++ b/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; + using OutputImageType = itk::Image; + + 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; + 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 it(filter->GetOutput(), filter->GetOutput()->GetBufferedRegion()); + for (it.GoToBegin(); !it.IsAtEnd(); ++it) + { + ASSERT_TRUE(std::isfinite(it.Get())); + } +} diff --git a/Modules/Filtering/ImageStatistics/include/itkAdaptiveEqualizationHistogram.h b/Modules/Filtering/ImageStatistics/include/itkAdaptiveEqualizationHistogram.h index fd84726eb19..7c587a1a6a4 100644 --- a/Modules/Filtering/ImageStatistics/include/itkAdaptiveEqualizationHistogram.h +++ b/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(m_Maximum) - m_Minimum; + if (iscale == 0.0) + { + // Constant image: equalization is an identity mapping. + return static_cast(pixel); + } double sum = 0.0; auto itMap = m_Map.begin(); diff --git a/Modules/Filtering/ImageStatistics/include/itkLabelOverlapMeasuresImageFilter.hxx b/Modules/Filtering/ImageStatistics/include/itkLabelOverlapMeasuresImageFilter.hxx index bbcac5cbdb1..a565f4d2dc7 100644 --- a/Modules/Filtering/ImageStatistics/include/itkLabelOverlapMeasuresImageFilter.hxx +++ b/Modules/Filtering/ImageStatistics/include/itkLabelOverlapMeasuresImageFilter.hxx @@ -385,17 +385,16 @@ LabelOverlapMeasuresImageFilter::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::max(); } else { - auto nComplementIntersection = nVox - mapIt->second.m_Union; // TN - - value = static_cast(mapIt->second.m_SourceComplement) / - static_cast(mapIt->second.m_SourceComplement + nComplementIntersection); + value = static_cast(mapIt->second.m_SourceComplement) / static_cast(denominator); } return value; diff --git a/Modules/Filtering/ImageStatistics/test/CMakeLists.txt b/Modules/Filtering/ImageStatistics/test/CMakeLists.txt index 41c49a34924..5eb1da61d20 100644 --- a/Modules/Filtering/ImageStatistics/test/CMakeLists.txt +++ b/Modules/Filtering/ImageStatistics/test/CMakeLists.txt @@ -467,6 +467,7 @@ itk_add_test( set( ITKImageStatisticsGTests + itkAdaptiveHistogramEqualizationImageFilterGTest.cxx itkLabelOverlapMeasuresImageFilterGTest.cxx itkMinimumMaximumImageFilterGTest.cxx ) diff --git a/Modules/Filtering/ImageStatistics/test/itkAdaptiveHistogramEqualizationImageFilterGTest.cxx b/Modules/Filtering/ImageStatistics/test/itkAdaptiveHistogramEqualizationImageFilterGTest.cxx new file mode 100644 index 00000000000..20b55421576 --- /dev/null +++ b/Modules/Filtering/ImageStatistics/test/itkAdaptiveHistogramEqualizationImageFilterGTest.cxx @@ -0,0 +1,44 @@ +/*========================================================================= + * + * 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 "itkAdaptiveHistogramEqualizationImageFilter.h" +#include "itkImageRegionConstIterator.h" +#include "itkGTest.h" + +// A constant image has zero intensity range; equalization must be an identity mapping +// instead of dividing by zero and producing NaN everywhere (issue #6575, B9). +TEST(AdaptiveHistogramEqualizationImageFilter, ConstantImageIsIdentity) +{ + using ImageType = itk::Image; + auto image = ImageType::New(); + image->SetRegions(ImageType::RegionType{ ImageType::IndexType{}, ImageType::SizeType::Filled(16) }); + image->Allocate(); + image->FillBuffer(7.0F); + + using FilterType = itk::AdaptiveHistogramEqualizationImageFilter; + auto filter = FilterType::New(); + filter->SetInput(image); + filter->SetRadius(2); + filter->Update(); + + itk::ImageRegionConstIterator it(filter->GetOutput(), filter->GetOutput()->GetBufferedRegion()); + for (it.GoToBegin(); !it.IsAtEnd(); ++it) + { + ASSERT_EQ(it.Get(), 7.0F); + } +} diff --git a/Modules/Filtering/ImageStatistics/test/itkLabelOverlapMeasuresImageFilterGTest.cxx b/Modules/Filtering/ImageStatistics/test/itkLabelOverlapMeasuresImageFilterGTest.cxx index eddd1960e08..e2638d7f571 100644 --- a/Modules/Filtering/ImageStatistics/test/itkLabelOverlapMeasuresImageFilterGTest.cxx +++ b/Modules/Filtering/ImageStatistics/test/itkLabelOverlapMeasuresImageFilterGTest.cxx @@ -204,3 +204,19 @@ TEST_F(LabelOverlapMeasuresImageFilterFixture, test3) EXPECT_NEAR(filter->GetFalsePositiveError(2), 4.76837612950e-07, 1e-17); EXPECT_NEAR(filter->GetFalseDiscoveryRate(2), 1.0 / 3.0, 0.0); } + +// A label covering the whole image has FP = TN = 0; the zero denominator must map to +// the degenerate-value convention, not NaN (issue #6575, B12). +TEST_F(LabelOverlapMeasuresImageFilterFixture, WholeImageLabelFalsePositiveError) +{ + auto source = Utils::CreateImage(1); + auto target = Utils::CreateImage(1); + + auto filter = Utils::FilterType::New(); + filter->SetSourceImage(source); + filter->SetTargetImage(target); + filter->Update(); + + using RealType = Utils::FilterType::RealType; + EXPECT_EQ(filter->GetFalsePositiveError(1), itk::NumericTraits::max()); +} diff --git a/Modules/Filtering/Smoothing/include/itkBoxUtilities.h b/Modules/Filtering/Smoothing/include/itkBoxUtilities.h index e7476b0679e..151c49d3a54 100644 --- a/Modules/Filtering/Smoothing/include/itkBoxUtilities.h +++ b/Modules/Filtering/Smoothing/include/itkBoxUtilities.h @@ -473,7 +473,9 @@ BoxSigmaCalculatorFunction(const TInputImage * accImage, ++(cornerItVec[k]); } - oIt.Set(static_cast(std::sqrt((squareSum - sum * sum / pixelscount) / (pixelscount - 1)))); + oIt.Set(pixelscount > 1 + ? static_cast(std::sqrt((squareSum - sum * sum / pixelscount) / (pixelscount - 1))) + : OutputPixelType{}); } } else @@ -538,8 +540,9 @@ BoxSigmaCalculatorFunction(const TInputImage * accImage, } } - oIt.Set( - static_cast(std::sqrt((squareSum - sum * sum / edgepixelscount) / (edgepixelscount - 1)))); + oIt.Set(edgepixelscount > 1 ? static_cast( + std::sqrt((squareSum - sum * sum / edgepixelscount) / (edgepixelscount - 1))) + : OutputPixelType{}); } } } diff --git a/Modules/Filtering/Smoothing/test/CMakeLists.txt b/Modules/Filtering/Smoothing/test/CMakeLists.txt index d3864da347e..19f3cd69ba6 100644 --- a/Modules/Filtering/Smoothing/test/CMakeLists.txt +++ b/Modules/Filtering/Smoothing/test/CMakeLists.txt @@ -396,6 +396,7 @@ itk_add_test( set( ITKSmoothingGTests + itkBoxSigmaImageFilterGTest.cxx itkMeanImageFilterGTest.cxx itkMedianImageFilterGTest.cxx ) diff --git a/Modules/Filtering/Smoothing/test/itkBoxSigmaImageFilterGTest.cxx b/Modules/Filtering/Smoothing/test/itkBoxSigmaImageFilterGTest.cxx new file mode 100644 index 00000000000..7746e833d55 --- /dev/null +++ b/Modules/Filtering/Smoothing/test/itkBoxSigmaImageFilterGTest.cxx @@ -0,0 +1,47 @@ +/*========================================================================= + * + * 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 "itkBoxSigmaImageFilter.h" +#include "itkImageRegionConstIteratorWithIndex.h" +#include "itkGTest.h" + +// A radius-0 window holds a single sample, so the (unbiased) standard deviation must be +// zero, not sqrt(0/0) = NaN (issue #6575, B20). +TEST(BoxSigmaImageFilter, RadiusZeroYieldsZeroNotNaN) +{ + using ImageType = itk::Image; + auto image = ImageType::New(); + image->SetRegions(ImageType::RegionType{ ImageType::IndexType{}, ImageType::SizeType::Filled(8) }); + image->Allocate(); + for (itk::ImageRegionIteratorWithIndex it(image, image->GetBufferedRegion()); !it.IsAtEnd(); ++it) + { + it.Set(static_cast(it.GetIndex()[0] + 10 * it.GetIndex()[1])); + } + + using FilterType = itk::BoxSigmaImageFilter; + auto filter = FilterType::New(); + filter->SetInput(image); + filter->SetRadius(0); + filter->Update(); + + itk::ImageRegionConstIteratorWithIndex it(filter->GetOutput(), filter->GetOutput()->GetBufferedRegion()); + for (it.GoToBegin(); !it.IsAtEnd(); ++it) + { + EXPECT_EQ(it.Get(), 0.0F) << " at index " << it.GetIndex(); + } +} diff --git a/Modules/Numerics/Statistics/include/itkGaussianRandomSpatialNeighborSubsampler.hxx b/Modules/Numerics/Statistics/include/itkGaussianRandomSpatialNeighborSubsampler.hxx index 987d519963b..44a47daeb8d 100644 --- a/Modules/Numerics/Statistics/include/itkGaussianRandomSpatialNeighborSubsampler.hxx +++ b/Modules/Numerics/Statistics/include/itkGaussianRandomSpatialNeighborSubsampler.hxx @@ -48,14 +48,14 @@ GaussianRandomSpatialNeighborSubsampler::GetIntegerVariate(Ran itkExceptionMacro("upperBound (" << upperBound << ") not >= to lowerBound(" << lowerBound << ')'); } - RandomIntType randInt = 0; - + // Reject out-of-bound variates before the integer conversion: converting a + // negative value to an unsigned type is undefined behavior. + RealType floorVar = 0; do { - const RealType randVar = this->m_RandomNumberGenerator->GetNormalVariate(mean, m_Variance); - randInt = static_cast(std::floor(randVar)); - } while ((randInt < lowerBound) || (randInt > upperBound)); - return randInt; + floorVar = std::floor(this->m_RandomNumberGenerator->GetNormalVariate(mean, m_Variance)); + } while ((floorVar < static_cast(lowerBound)) || (floorVar > static_cast(upperBound))); + return static_cast(floorVar); } template diff --git a/Modules/Segmentation/SuperPixel/include/itkSLICImageFilter.hxx b/Modules/Segmentation/SuperPixel/include/itkSLICImageFilter.hxx index 443c02af2bd..576217ca9e4 100644 --- a/Modules/Segmentation/SuperPixel/include/itkSLICImageFilter.hxx +++ b/Modules/Segmentation/SuperPixel/include/itkSLICImageFilter.hxx @@ -642,10 +642,17 @@ SLICImageFilter::GenerateData() for (size_t i = 0; i * numberOfClusterComponents < m_Clusters.size(); ++i) { - ClusterType cluster(numberOfClusterComponents, &m_Clusters[i * numberOfClusterComponents]); - cluster /= clusterCount[i]; - + ClusterType cluster(numberOfClusterComponents, &m_Clusters[i * numberOfClusterComponents]); const ClusterType oldCluster(numberOfClusterComponents, &m_OldClusters[i * numberOfClusterComponents]); + if (clusterCount[i] > 0) + { + cluster /= clusterCount[i]; + } + else + { + // Empty cluster: keep the previous centroid to avoid a NaN center. + cluster.copy_in(oldCluster.data_block()); + } l1Residual += Distance(cluster, oldCluster); }