BUG: Fix StatisticsLabelMap median for even number of pixels

When there are an even number of elements, average the bins measurements vectors of the middle two.
InsightSoftwareConsortium · May 4, 2021 · 3abace1 · 3abace1
1 parent 55d0fdf
commit 3abace1
Show file tree

Hide file tree

Showing 2 changed files with 52 additions and 3 deletions.
diff --git a/Modules/Filtering/LabelMap/include/itkStatisticsLabelMapFilter.hxx b/Modules/Filtering/LabelMap/include/itkStatisticsLabelMapFilter.hxx
@@ -177,15 +177,28 @@ StatisticsLabelMapFilter<TImage, TFeatureImage>::ThreadedProcessLabelObject(Labe
   }
 
   // the median
-  double median = 0;
-  double count = 0; // will not be fully set, so do not use later !
-  for (SizeValueType i = 0; i < histogram->Size(); i++)
+  double median = 0.0;
+  double count = 0.0; // will not be fully set, so do not use later !
+  for (SizeValueType i = 0; i < histogram->Size(); ++i)
   {
     count += histogram->GetFrequency(i);
 
     if (count >= (totalFreq / 2))
     {
       median = histogram->GetMeasurementVector(i)[0];
+      // If there are an even number of elements average with the next bin with elements
+      if (labelObject->Size() % 2 == 0 && count == totalFreq / 2)
+      {
+        while (++i < histogram->Size())
+        {
+          if (histogram->GetFrequency(i) > 0)
+          {
+            median += histogram->GetMeasurementVector(i)[0];
+            median *= 0.5;
+            break;
+          }
+        }
+      }
       break;
     }
   }

diff --git a/Modules/Filtering/LabelMap/test/itkStatisticsLabelMapFilterGTest.cxx b/Modules/Filtering/LabelMap/test/itkStatisticsLabelMapFilterGTest.cxx
@@ -249,3 +249,39 @@ TEST_F(StatisticsLabelMapFixture, 2D_rand_with_outliers)
     labelObject->Print(std::cout);
   }
 }
+
+
+TEST_F(StatisticsLabelMapFixture, 2D_even)
+{
+  using Utils = FixtureUtilities<2, unsigned char>;
+
+  auto image = Utils::CreateImage();
+  auto labelImage = Utils ::CreateLabelImage();
+
+  // Set label with two elements far apart, the median should be average
+  image->SetPixel({ 0, 0 }, 10);
+  image->SetPixel({ 0, 1 }, 100);
+  image->SetPixel({ 0, 2 }, 1);
+  image->SetPixel({ 0, 3 }, 200);
+
+  Utils::LabelPixelType label = 1;
+  labelImage->SetPixel({ 0, 0 }, label);
+  labelImage->SetPixel({ 0, 1 }, label);
+  labelImage->SetPixel({ 0, 2 }, label);
+  labelImage->SetPixel({ 0, 3 }, label);
+
+
+  Utils::LabelObjectType::ConstPointer labelObject = Utils::ComputeLabelObject(labelImage, image, label, 1 << 8);
+
+  EXPECT_NEAR(1.0, labelObject->GetMinimum(), 1e-12);
+  EXPECT_NEAR(200.0, labelObject->GetMaximum(), 1e-12);
+  EXPECT_NEAR(Utils::ComputeExactMedian(labelObject, image), labelObject->GetMedian(), 0.5);
+  EXPECT_NEAR(311.0, labelObject->GetSum(), 1e-12);
+  EXPECT_NEAR(8640.25, labelObject->GetVariance(), 1e-10);
+  EXPECT_NEAR(92.95294, labelObject->GetStandardDeviation(), 1e-5);
+
+  if (::testing::Test::HasFailure())
+  {
+    labelObject->Print(std::cout);
+  }
+}