InsightSoftwareConsortium / ITK

Insight Toolkit (ITK) -- Official Repository. ITK builds on a proven, spatially-oriented architecture for processing, segmentation, and registration of scientific images in two, three, or more dimensions.
https://itk.org
Apache License 2.0
1.37k stars 660 forks source link

DOC: Add See Also entries for alternatives to LabelGeometryImageFilter #4624

Closed dzenanz closed 2 months ago

dzenanz commented 2 months ago

This filter is buggy, and sometimes produces wrong results. Follow-up to: #4530.

PR Checklist

dzenanz commented 2 months ago

Today I realized its IntegratedIntensity is buggy too. Should we deprecate this filter?

blowekamp commented 2 months ago

Today I realized its IntegratedIntensity is buggy too. Should we deprecate this filter?

Yes, I have been a proponent of that for nearly 10 years.

dzenanz commented 2 months ago

@blowekamp do you want to propose a PR which does that? Then we can decide when to do the deprecation. Most likely right after release? And should we go future legacy remove route? Or direct deprecation?

blowekamp commented 2 months ago

We could add the depreciation warning now, and move it to Compatibility/Deprecated for ITKv5.

cookpa commented 2 months ago

Hi, Thanks for the heads up on the bugginess of this filter.

I've been working on replacing LabelGeometryImageFilter in ANTs tools. The suggested replacements work well for most functions, but I have not been able to find an equivalent for the "Eccentricity" and the "WeightedCentroid". Are there options for this in any other ITK filters?

dzenanz commented 2 months ago

Perhaps Roundness from ShapeLabelObject, which get produced by LabelImageToShapeLabelMapFilter? I don't know of a replacement for WeightedCentroid.

blowekamp commented 2 months ago

Hi, Thanks for the heads up on the bugginess of this filter.

I've been working on replacing LabelGeometryImageFilter in ANTs tools. The suggested replacements work well for most functions, but I have not been able to find an equivalent for the "Eccentricity" and the "WeightedCentroid". Are there options for this in any other ITK filters?

The other filters to mention is the LabelImageToStatisticsLabelMapFilter and relates BinaryToStatisticLabelMapFilter and StatisticsLabelMapFilter. The take a label image and an intensity image to compute the attributes here: https://itk.org/Doxygen/html/classitk_1_1StatisticsLabelObject.html

cookpa commented 2 months ago

Thanks @dzenanz and @blowekamp. I switched the shape measures over to LabelImageToStatisticsLabelMapFilter, and the the grayscale measures to LabelStatisticsImageFilter.

blowekamp commented 2 months ago

Thanks @dzenanz and @blowekamp. I switched the shape measures over to LabelImageToStatisticsLabelMapFilter, and the the grayscale measures to LabelStatisticsImageFilter.

If you could post specifics of the metrics migration, it could help writing a migration guide.

cookpa commented 2 months ago

Specifically, we were calling these LabelGeometryImageFilter functions to get shape statistics:

GetAxesLength GetBoundingBox GetCentroid GetEccentricity GetElongation GetOrientation

LabelImageToShapeLabelMapFilter provides a LabelMap which in turn has ShapeLabelObject objects for each unique label, which provide

GetBoundingBox GetCentroid GetElongation

Of these, only BoundingBox was the same in my tests. GetCentroid differs and appears to provide the centroid in physical space. GetElongation returns different numbers, I've not verified which is correct.

There's no GetEccentricity, but there is GetRoundness, which appears to be a more general way of measuring something inversely proportional to eccentricity.

AFAIK, there's not a replacement for GetAxesLength or GetOrientation.

For grayscale statistics, LabelGeometryMeasures has

GetIntegratedIntensity GetWeightedCentroid

I replaced GetIntegratedIntensity with GetSum from the LabelStatisticsImageFilter, which appears to be the same. I don't have a replacement for GetWeightedCentroid.

cookpa commented 2 months ago

Update, I was able to use ShapeLabelObject::GetPrincipalMoments() to reproduce axes length as

axesLengths[idx] = 4.0 * std::sqrt(principalMoments[idx]);

giving the same result as LabelGeometryImageFilter for unit voxels (EDIT: LabelGeometryImageFilter does not scale axes length with the spacing, so the same label will show different axes lengths if the voxel spacing changes).

The moments can also be used to define eccentricity and elongation, and in 2D, orientation. From LabelGeometryImageFilter

https://github.com/InsightSoftwareConsortium/ITK/blob/e5a84a08a3092fe9db81623c7691ba557c2a832b/Modules/Nonunit/Review/include/itkLabelGeometryImageFilter.hxx#L392-L400

Eccentricity is defined as

mapIt->second.m_Eccentricity =
      std::sqrt((eigenvalues[ImageDimension - 1] - eigenvalues[0]) / eigenvalues[ImageDimension - 1]);

which simplifies to sqrt( 1 - V_0 / V1). This differs from the [ellipsoid formula](https://en.wikipedia.org/wiki/Eccentricity(mathematics)#Ellipses) sqrt(1 - {V_0}^2 / {V_1}^2).

In 3D there is no single definition but I decided to use the ratio of the largest to smallest eigenvalues, ie eccentricity = sqrt(1 - {V_0}^2 / {V_2}^2).

For elongation, LabelGeometryImageFilter uses

mapIt->second.m_Elongation = axesLength[ImageDimension - 1] / axesLength[0];

Which is (4 sqrt(V_1)) / (4 sqrt(V_0)) = sqrt(V_1 / V_0).

itkShapeLabelMapFilter uses

https://github.com/InsightSoftwareConsortium/ITK/blob/e5a84a08a3092fe9db81623c7691ba557c2a832b/Modules/Filtering/LabelMap/include/itkShapeLabelMapFilter.hxx#L364-L373

which is the same in 2D but in 3D, elongation becomes sqrt(V_2 / V_1) rather than sqrt(V_2 / V_0). This is arguably better for some applications, for example looking at cortical regions we expect V_2 >= V_1 >> V_0, and might be more interested in the elongation of the first and second moments. So I guess the difference is just an implementation choice.