BUG: Fix LabelMap and label-statistics correctness bugs from issue 6575

Modules/Filtering/ImageFusion/include/itkLabelMapContourOverlayImageFilter.hxx

@@ -155,7 +155,7 @@ LabelMapContourOverlayImageFilter<TLabelMap, TFeatureImage, TOutputImage>::Befor
   RadiusType srad{};
   for (unsigned int i = 0, j = 0; i < ImageDimension; ++i)
   {
-    if (j != static_cast<unsigned int>(m_SliceDimension) && (j < (ImageDimension - 1)))
+    if (i != static_cast<unsigned int>(m_SliceDimension) && (j < (ImageDimension - 1)))
     {
       srad[j] = m_ContourThickness[i];
       ++j;

Modules/Filtering/LabelMap/include/itkAttributeUniqueLabelMapFilter.hxx

@@ -236,7 +236,7 @@ AttributeUniqueLabelMapFilter<TImage, TAttributeAccessor>::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();

Modules/Filtering/LabelMap/include/itkLabelMap.hxx

@@ -414,6 +414,7 @@ LabelMap<TLabelObject>::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<TLabelObject>::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.");
       }

Modules/Filtering/LabelMap/include/itkLabelMapMaskImageFilter.hxx

@@ -140,10 +140,19 @@ LabelMapMaskImageFilter<TInputImage, TOutputImage>::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<TInputImage, TOutputImage>::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 };
       }

Modules/Filtering/LabelMap/include/itkMergeLabelMapFilter.hxx

@@ -71,12 +71,12 @@ MergeLabelMapFilter<TImage>::MergeWithKeep()
   ImageType * output = this->GetOutput();
 
   using VectorType = std::deque<LabelObjectPointer>;
-  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())
     {

Modules/Filtering/LabelMap/include/itkStatisticsLabelMapFilter.hxx

@@ -233,7 +233,7 @@ StatisticsLabelMapFilter<TImage, TFeatureImage>::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<TImage, TFeatureImage>::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

Modules/Filtering/LabelMap/test/CMakeLists.txt

@@ -1796,6 +1796,7 @@ itk_add_test(
 set(
   ITKLabelMapGTests
   itkShapeLabelMapFilterGTest.cxx
+  itkLabelMapBugfixGTest.cxx
   itkStatisticsLabelMapFilterGTest.cxx
   itkUniqueLabelMapFiltersGTest.cxx
 )

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<unsigned char, 2>;
+  using LabelMapType = itk::LabelMap<LabelObjectType>;
+
+  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<unsigned char>(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<unsigned char, 2>;
+  using LabelMapType = itk::LabelMap<LabelObjectType>;
+  using FeatureImageType = itk::Image<float, 2>;
+
+  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<LabelMapType, FeatureImageType>::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<unsigned char, 2>;
+  using LabelMapType = itk::LabelMap<LabelObjectType>;
+
+  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<LabelMapType>::New();
+  unique->SetInput(labelMap);
+  unique->Update();
+
+  EXPECT_EQ(unique->GetOutput()->GetNumberOfLabelObjects(), 1u);
+}