diff --git a/Modules/Filtering/ImageFusion/include/itkLabelMapContourOverlayImageFilter.hxx b/Modules/Filtering/ImageFusion/include/itkLabelMapContourOverlayImageFilter.hxx index a7c8403a5d5..663ec0ad17b 100644 --- a/Modules/Filtering/ImageFusion/include/itkLabelMapContourOverlayImageFilter.hxx +++ b/Modules/Filtering/ImageFusion/include/itkLabelMapContourOverlayImageFilter.hxx @@ -155,7 +155,7 @@ LabelMapContourOverlayImageFilter::Befor RadiusType srad{}; for (unsigned int i = 0, j = 0; i < ImageDimension; ++i) { - if (j != static_cast(m_SliceDimension) && (j < (ImageDimension - 1))) + if (i != static_cast(m_SliceDimension) && (j < (ImageDimension - 1))) { srad[j] = m_ContourThickness[i]; ++j; diff --git a/Modules/Filtering/LabelMap/include/itkAttributeUniqueLabelMapFilter.hxx b/Modules/Filtering/LabelMap/include/itkAttributeUniqueLabelMapFilter.hxx index 23403e62db7..9f0ae5b5b6d 100644 --- a/Modules/Filtering/LabelMap/include/itkAttributeUniqueLabelMapFilter.hxx +++ b/Modules/Filtering/LabelMap/include/itkAttributeUniqueLabelMapFilter.hxx @@ -236,7 +236,7 @@ AttributeUniqueLabelMapFilter::GenerateData() // remove objects without lines typename ImageType::Iterator it(this->GetLabelMap()); - while (it.IsAtEnd()) + while (!it.IsAtEnd()) { const typename LabelObjectType::LabelType label = it.GetLabel(); LabelObjectType * labelObject = it.GetLabelObject(); diff --git a/Modules/Filtering/LabelMap/include/itkLabelMap.hxx b/Modules/Filtering/LabelMap/include/itkLabelMap.hxx index ccd460da2f4..6a5beb72ac4 100644 --- a/Modules/Filtering/LabelMap/include/itkLabelMap.hxx +++ b/Modules/Filtering/LabelMap/include/itkLabelMap.hxx @@ -414,6 +414,7 @@ LabelMap::PushLabelObject(LabelObjectType * labelObject) // search for an unused label LabelType label = firstLabel; LabelObjectContainerConstIterator it; + bool foundFreeLabel = false; for (it = m_LabelObjectContainer.begin(); it != m_LabelObjectContainer.end(); ++it, ++label) { assert((it->second.IsNotNull())); @@ -424,10 +425,11 @@ LabelMap::PushLabelObject(LabelObjectType * labelObject) if (label != it->first) { labelObject->SetLabel(label); + foundFreeLabel = true; break; } } - if (label == lastLabel) + if (!foundFreeLabel) { itkExceptionStringMacro("Can't push the label object: the label map is full."); } diff --git a/Modules/Filtering/LabelMap/include/itkLabelMapMaskImageFilter.hxx b/Modules/Filtering/LabelMap/include/itkLabelMapMaskImageFilter.hxx index 9263c7b30b9..152d98f6034 100644 --- a/Modules/Filtering/LabelMap/include/itkLabelMapMaskImageFilter.hxx +++ b/Modules/Filtering/LabelMap/include/itkLabelMapMaskImageFilter.hxx @@ -140,10 +140,19 @@ LabelMapMaskImageFilter::GenerateOutputInformation() } // Final computation - SizeType regionSize; - for (unsigned int i = 0; i < ImageDimension; ++i) + SizeType regionSize{}; + if (maxs[0] >= mins[0]) { - regionSize[i] = maxs[i] - mins[i] + 1; + for (unsigned int i = 0; i < ImageDimension; ++i) + { + regionSize[i] = maxs[i] - mins[i] + 1; + } + } + else + { + // No matching object lines: use an empty region rather than overflowing + // the size with the uninitialized min/max sentinels. + mins.Fill(0); } cropRegion = { mins, regionSize }; } @@ -194,10 +203,19 @@ LabelMapMaskImageFilter::GenerateOutputInformation() ++lit; } // Final computation - SizeType regionSize; - for (unsigned int i = 0; i < ImageDimension; ++i) + SizeType regionSize{}; + if (maxs[0] >= mins[0]) + { + for (unsigned int i = 0; i < ImageDimension; ++i) + { + regionSize[i] = maxs[i] - mins[i] + 1; + } + } + else { - regionSize[i] = maxs[i] - mins[i] + 1; + // No matching object lines: use an empty region rather than overflowing + // the size with the uninitialized min/max sentinels. + mins.Fill(0); } cropRegion = { mins, regionSize }; } diff --git a/Modules/Filtering/LabelMap/include/itkMergeLabelMapFilter.hxx b/Modules/Filtering/LabelMap/include/itkMergeLabelMapFilter.hxx index 15ef7f74b9e..22e9df3ee42 100644 --- a/Modules/Filtering/LabelMap/include/itkMergeLabelMapFilter.hxx +++ b/Modules/Filtering/LabelMap/include/itkMergeLabelMapFilter.hxx @@ -71,12 +71,12 @@ MergeLabelMapFilter::MergeWithKeep() ImageType * output = this->GetOutput(); using VectorType = std::deque; - VectorType labelObjects; ProgressReporter progress(this, 0, 1); for (unsigned int i = 1; i < this->GetNumberOfIndexedInputs(); ++i) { + VectorType labelObjects; typename ImageType::ConstIterator it2(this->GetInput(i)); while (!it2.IsAtEnd()) { diff --git a/Modules/Filtering/LabelMap/include/itkStatisticsLabelMapFilter.hxx b/Modules/Filtering/LabelMap/include/itkStatisticsLabelMapFilter.hxx index 0e99d00ee2f..1e563fc752c 100644 --- a/Modules/Filtering/LabelMap/include/itkStatisticsLabelMapFilter.hxx +++ b/Modules/Filtering/LabelMap/include/itkStatisticsLabelMapFilter.hxx @@ -233,7 +233,7 @@ StatisticsLabelMapFilter::ThreadedProcessLabelObject(Labe for (unsigned int i = 0; i < ImageDimension; ++i) { // principalMoments[i] = 4 * std::sqrt( pm(i,i) ); - principalMoments[i] = pm(i); + principalMoments[i] = std::max(pm(i), 0.0); } principalAxes = eigen.V.transpose(); @@ -258,12 +258,18 @@ StatisticsLabelMapFilter::ThreadedProcessLabelObject(Labe elongation = 1; flatness = 1; } - else if (Math::NotAlmostEquals(principalMoments[0], typename VectorType::ValueType{})) + else { - // elongation = principalMoments[ImageDimension-1] / - // principalMoments[0]; - elongation = std::sqrt(principalMoments[ImageDimension - 1] / principalMoments[ImageDimension - 2]); - flatness = std::sqrt(principalMoments[1] / principalMoments[0]); + if (Math::NotAlmostEquals(principalMoments[0], typename VectorType::ValueType{})) + { + const double flatnessRatio = principalMoments[1] / principalMoments[0]; + flatness = (flatnessRatio > 0.0) ? std::sqrt(flatnessRatio) : 0.0; + } + if (Math::NotAlmostEquals(principalMoments[ImageDimension - 2], typename VectorType::ValueType{})) + { + const double elongationRatio = principalMoments[ImageDimension - 1] / principalMoments[ImageDimension - 2]; + elongation = (elongationRatio > 0.0) ? std::sqrt(elongationRatio) : 0.0; + } } } else diff --git a/Modules/Filtering/LabelMap/test/CMakeLists.txt b/Modules/Filtering/LabelMap/test/CMakeLists.txt index 55e224e5816..eda0f129c6a 100644 --- a/Modules/Filtering/LabelMap/test/CMakeLists.txt +++ b/Modules/Filtering/LabelMap/test/CMakeLists.txt @@ -1796,6 +1796,7 @@ itk_add_test( set( ITKLabelMapGTests itkShapeLabelMapFilterGTest.cxx + itkLabelMapBugfixGTest.cxx itkStatisticsLabelMapFilterGTest.cxx itkUniqueLabelMapFiltersGTest.cxx ) diff --git a/Modules/Filtering/LabelMap/test/itkLabelMapBugfixGTest.cxx b/Modules/Filtering/LabelMap/test/itkLabelMapBugfixGTest.cxx new file mode 100644 index 00000000000..5761b10c185 --- /dev/null +++ b/Modules/Filtering/LabelMap/test/itkLabelMapBugfixGTest.cxx @@ -0,0 +1,110 @@ +/*========================================================================= + * + * 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 "itkLabelMap.h" +#include "itkLabelObject.h" +#include "itkStatisticsLabelObject.h" +#include "itkStatisticsLabelMapFilter.h" +#include "itkLabelUniqueLabelMapFilter.h" +#include "itkGTest.h" + +// A completely full label range must make PushLabelObject throw, not silently +// overwrite an existing entry (issue #6575, B44). +TEST(LabelMap, PushLabelObjectThrowsWhenFull) +{ + using LabelObjectType = itk::LabelObject; + using LabelMapType = itk::LabelMap; + + auto labelMap = LabelMapType::New(); + labelMap->SetRegions(LabelMapType::RegionType{ LabelMapType::IndexType{}, LabelMapType::SizeType::Filled(4) }); + + // Background is 0; occupy every remaining label 1..255. + for (unsigned int label = 1; label <= 255; ++label) + { + labelMap->SetLine(LabelMapType::IndexType{}, 1, static_cast(label)); + } + ASSERT_EQ(labelMap->GetNumberOfLabelObjects(), 255u); + + auto extra = LabelObjectType::New(); + extra->AddLine(itk::MakeIndex(0, 1), 1); + EXPECT_THROW(labelMap->PushLabelObject(extra), itk::ExceptionObject); +} + +// Weighted elongation/flatness must stay finite when a signed feature image makes +// the weighted covariance indefinite (issue #6575, B45). +TEST(StatisticsLabelMapFilter, WeightedShapeFinelyDefinedForSignedFeature) +{ + using LabelObjectType = itk::StatisticsLabelObject; + using LabelMapType = itk::LabelMap; + using FeatureImageType = itk::Image; + + auto labelMap = LabelMapType::New(); + labelMap->SetRegions(LabelMapType::RegionType{ LabelMapType::IndexType{}, LabelMapType::SizeType::Filled(9) }); + // A cross: horizontal arm (y == 4, x in 2..6) and vertical arm (x == 4, y in 2,3,5,6). + labelMap->SetLine(itk::MakeIndex(2, 4), 5, 1); + labelMap->SetLine(itk::MakeIndex(4, 2), 1, 1); + labelMap->SetLine(itk::MakeIndex(4, 3), 1, 1); + labelMap->SetLine(itk::MakeIndex(4, 5), 1, 1); + labelMap->SetLine(itk::MakeIndex(4, 6), 1, 1); + + // Positive weights along the x arm, negative along the y arm: the weighted + // second moment is positive in xx and negative in yy -> indefinite covariance. + auto feature = FeatureImageType::New(); + feature->SetRegions(labelMap->GetLargestPossibleRegion()); + feature->Allocate(); + feature->FillBuffer(0.0F); + for (int x = 2; x <= 6; ++x) + { + feature->SetPixel(itk::MakeIndex(x, 4), 2.0F); + } + for (int y : { 2, 3, 5, 6 }) + { + feature->SetPixel(itk::MakeIndex(4, y), -1.0F); + } + + auto filter = itk::StatisticsLabelMapFilter::New(); + filter->SetInput(labelMap); + filter->SetFeatureImage(feature); + filter->Update(); + + const auto * object = filter->GetOutput()->GetLabelObject(1); + EXPECT_TRUE(std::isfinite(object->GetWeightedElongation())) << " weighted elongation"; + EXPECT_TRUE(std::isfinite(object->GetWeightedFlatness())) << " weighted flatness"; +} + +// A label object whose every pixel is claimed by a higher-priority object becomes +// empty and must be removed; the removal pass used an inverted loop condition and +// never ran (issue #6575, B43). +TEST(LabelUniqueLabelMapFilter, EmptiedObjectIsRemoved) +{ + using LabelObjectType = itk::LabelObject; + using LabelMapType = itk::LabelMap; + + auto labelMap = LabelMapType::New(); + labelMap->SetRegions(LabelMapType::RegionType{ LabelMapType::IndexType{}, LabelMapType::SizeType::Filled(8) }); + // Two objects covering exactly the same pixels: after uniqueness one is emptied. + labelMap->SetLine(itk::MakeIndex(0, 0), 4, 1); + labelMap->SetLine(itk::MakeIndex(0, 0), 4, 2); + ASSERT_EQ(labelMap->GetNumberOfLabelObjects(), 2u); + + auto unique = itk::LabelUniqueLabelMapFilter::New(); + unique->SetInput(labelMap); + unique->Update(); + + EXPECT_EQ(unique->GetOutput()->GetNumberOfLabelObjects(), 1u); +}