Skip to content

Commit

Permalink
ENH: StatisticsLabelMapFilter use improve integer histogram
Browse files Browse the repository at this point in the history 
For 8 and 16-bit integer set the number of histogram bins to the
number of values. Also add bin margins so that the center of the bins
exactly match the integer values.

Fix bin selected when there are an odd number of elements. Now
correctly choosing the middle bin.

This combination results in exact median computation for 8-bit and
16-bit integers.
  • Loading branch information
@blowekamp @Leengit
blowekamp authored and Leengit committed May 4, 2021
1 parent 3abace1 commit daa2a20
Show file tree
Hide file tree
Showing 3 changed files with 62 additions and 15 deletions.
13 changes: 9 additions & 4 deletions Modules/Filtering/LabelMap/include/itkStatisticsLabelMapFilter.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -147,10 +147,15 @@ class ITK_TEMPLATE_EXPORT StatisticsLabelMapFilter
PrintSelf(std::ostream & os, Indent indent) const override;

private:
FeatureImagePixelType m_Minimum;
FeatureImagePixelType m_Maximum;
unsigned int m_NumberOfBins;
bool m_ComputeHistogram;
FeatureImagePixelType m_Minimum{ NumericTraits<FeatureImagePixelType>::ZeroValue() };
FeatureImagePixelType m_Maximum{ NumericTraits<FeatureImagePixelType>::ZeroValue() };
// Set the number of bins to match the number of values for 8,16-bit integers.
static constexpr bool m_DefaultNumberOfBinsToNumberOfValues =
NumericTraits<typename Self::FeatureImagePixelType>::IsInteger && sizeof(typename Self::FeatureImagePixelType) <= 2;
unsigned int m_NumberOfBins{ m_DefaultNumberOfBinsToNumberOfValues
? 1 << (8 * sizeof(typename Self::FeatureImagePixelType))
: 128 };
bool m_ComputeHistogram{ true };
}; // end of class
} // end namespace itk

Expand Down
22 changes: 15 additions & 7 deletions Modules/Filtering/LabelMap/include/itkStatisticsLabelMapFilter.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,10 +31,6 @@ namespace itk
template <typename TImage, typename TFeatureImage>
StatisticsLabelMapFilter<TImage, TFeatureImage>::StatisticsLabelMapFilter()
{
m_Minimum = NumericTraits<FeatureImagePixelType>::ZeroValue();
m_Maximum = NumericTraits<FeatureImagePixelType>::ZeroValue();
m_NumberOfBins = 128;
m_ComputeHistogram = true;
this->SetNumberOfRequiredInputs(2);
}

Expand Down Expand Up @@ -73,10 +69,22 @@ StatisticsLabelMapFilter<TImage, TFeatureImage>::ThreadedProcessLabelObject(Labe
histogramSize.Fill(m_NumberOfBins);

typename HistogramType::MeasurementVectorType featureImageMin(1);
featureImageMin.Fill(m_Minimum);

typename HistogramType::MeasurementVectorType featureImageMax(1);
featureImageMax.Fill(m_Maximum);


if (NumericTraits<typename Self::FeatureImagePixelType>::IsInteger &&
m_NumberOfBins == 1 << (8 * sizeof(typename Self::FeatureImagePixelType)))
{
// Add padding so the center of bins are integers
featureImageMin.Fill(NumericTraits<typename Self::FeatureImagePixelType>::min() - 0.5);
featureImageMax.Fill(NumericTraits<typename Self::FeatureImagePixelType>::max() + 0.5);
}
else
{
featureImageMin.Fill(m_Minimum);
featureImageMax.Fill(m_Maximum);
}

typename HistogramType::Pointer histogram = HistogramType::New();
histogram->SetMeasurementVectorSize(1);
Expand Down Expand Up @@ -183,7 +191,7 @@ StatisticsLabelMapFilter<TImage, TFeatureImage>::ThreadedProcessLabelObject(Labe
{
count += histogram->GetFrequency(i);

if (count >= (totalFreq / 2))
if (count >= ((totalFreq + 1) / 2))
{
median = histogram->GetMeasurementVector(i)[0];
// If there are an even number of elements average with the next bin with elements
Expand Down
42 changes: 38 additions & 4 deletions Modules/Filtering/LabelMap/test/itkStatisticsLabelMapFilterGTest.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -175,7 +175,7 @@ TEST_F(StatisticsLabelMapFixture, 2D_zero)
EXPECT_NEAR(Utils::ComputeExactMedian(labelObject, image), labelObject->GetMedian(), 1e-12);
EXPECT_NEAR(0.0, labelObject->GetSum(), 1e-12);
EXPECT_NEAR(0.0, labelObject->GetVariance(), 1e-12);
EXPECT_NEAR(0.0, labelObject->GetStandardDeviation(), 1e-12);
EXPECT_NEAR(0.0, labelObject->GetStandardDeviation(), 0.5);

if (::testing::Test::HasFailure())
{
Expand Down Expand Up @@ -207,7 +207,7 @@ TEST_F(StatisticsLabelMapFixture, 2D_ones_with_outliers)

EXPECT_NEAR(value, labelObject->GetMinimum(), 1e-12);
EXPECT_NEAR(value, labelObject->GetMaximum(), 1e-12);
EXPECT_NEAR(Utils::ComputeExactMedian(labelObject, image), labelObject->GetMedian(), 0.5);
EXPECT_NEAR(Utils::ComputeExactMedian(labelObject, image), labelObject->GetMedian(), 1e-12);
EXPECT_NEAR(25 * 25 - 2, labelObject->GetSum(), 1e-12);
EXPECT_NEAR(0.0, labelObject->GetVariance(), 1e-12);
EXPECT_NEAR(0.0, labelObject->GetStandardDeviation(), 1e-12);
Expand Down Expand Up @@ -242,7 +242,7 @@ TEST_F(StatisticsLabelMapFixture, 2D_rand_with_outliers)

EXPECT_NEAR(0.0, labelObject->GetMinimum(), 1e-12);
EXPECT_NEAR(500.0, labelObject->GetMaximum(), 1e-12);
EXPECT_NEAR(Utils::ComputeExactMedian(labelObject, image), labelObject->GetMedian(), 0.5);
EXPECT_NEAR(Utils::ComputeExactMedian(labelObject, image), labelObject->GetMedian(), 1e-12);

if (::testing::Test::HasFailure())
{
Expand Down Expand Up @@ -275,7 +275,7 @@ TEST_F(StatisticsLabelMapFixture, 2D_even)

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(Utils::ComputeExactMedian(labelObject, image), labelObject->GetMedian(), 1e-12);
EXPECT_NEAR(311.0, labelObject->GetSum(), 1e-12);
EXPECT_NEAR(8640.25, labelObject->GetVariance(), 1e-10);
EXPECT_NEAR(92.95294, labelObject->GetStandardDeviation(), 1e-5);
Expand All @@ -285,3 +285,37 @@ TEST_F(StatisticsLabelMapFixture, 2D_even)
labelObject->Print(std::cout);
}
}


TEST_F(StatisticsLabelMapFixture, 2D_three)
{
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 }, 1);
image->SetPixel({ 0, 1 }, 3);
image->SetPixel({ 0, 2 }, 10);

Utils::LabelPixelType label = 1;
labelImage->SetPixel({ 0, 0 }, label);
labelImage->SetPixel({ 0, 1 }, label);
labelImage->SetPixel({ 0, 2 }, label);


Utils::LabelObjectType::ConstPointer labelObject = Utils::ComputeLabelObject(labelImage, image, label, 1 << 8);

EXPECT_NEAR(1.0, labelObject->GetMinimum(), 1e-12);
EXPECT_NEAR(10.0, labelObject->GetMaximum(), 1e-12);
EXPECT_NEAR(Utils::ComputeExactMedian(labelObject, image), labelObject->GetMedian(), 1e-12);
EXPECT_NEAR(14.0, labelObject->GetSum(), 1e-12);
EXPECT_NEAR(22.33333, labelObject->GetVariance(), 1e-5);
EXPECT_NEAR(4.725815, labelObject->GetStandardDeviation(), 1e-5);

if (::testing::Test::HasFailure())
{
labelObject->Print(std::cout);
}
}

0 comments on commit daa2a20

Please sign in to comment.