Deprecation of LabelGeometryImageFilter (used in LabelGeometryMeasures)

[cookpa]

InsightSoftwareConsortium/ITK#4630

Thanks @gdevenyi for the heads up.

I don't use LabelGeometryMeasures much because it doesn't do what most users want, which is estimate cortical surface area.

I was actually looking at this issue recently. If I understand this paper correctly,

https://onlinelibrary.wiley.com/doi/abs/10.1002/1361-6374%28199303%291%3A1%3C6%3A%3AAID-BIO3%3E3.0.CO%3B2-3

it's possible to approximate surface area from a voxel image, if the local geometry is known. What we need though is a filter that computes this only for voxels connected to background.

[stnava]

hmm - this would be a problem for me. @cookpa thickness = Volume / area .... so we compute area on the fly when we need to ... A is the surface area of one side of the sheet.

[cookpa]

Looks like this is a suggested replacement

https://itk.org/Doxygen/html/classitk_1_1LabelImageToShapeLabelMapFilter.html

[ntustison]

yes, thx @gdevenyi . @cookpa (or anyone else) do you want to work on replacement or would you like me to?

[cookpa]

I'll work on it, and do some comparisons on the results

[gdevenyi]

I'm happy to test as I use LabelGeometryMeasures for some work.

[cookpa]

As well as LabelGeometryMeasures, this filter is used in

ImageSegmentation/antsAtroposSegmentationImageFilter.hxx
Examples/ImageMath_Templates.hxx
Examples/GetConnectedComponentsFeatureImages.cxx
Examples/CreateDTICohort.cxx

[cookpa]

LabelGeometryImageFilter is used in the undocumented "LabelThickness2" function in ImageMath. It appears to compute average thickness based on the ratio of the volume to perimeter, then replaces the input labels with the average thickness for each label.

@ntustison @stnava do we still need this?

ANTs/Examples/ImageMath_Templates.hxx

Lines 10782 to 10867 in 4c68a48

	template <unsigned int ImageDimension>
	int
	LabelThickness2(int argc, char * argv[])
	{
	typedef unsigned int LabelType;
	typedef itk::Image<LabelType, ImageDimension> LabelImageType;
	typedef float RealType;
	typedef itk::Image<RealType, ImageDimension> RealImageType;
	if (argc < 5)
	{
	// std::cout << " Not enough inputs " << std::endl;
	return EXIT_FAILURE;
	}
	int argct = 2;
	const std::string outname = std::string(argv[argct]);
	argct += 2;
	std::string fn = std::string(argv[argct++]);
	typename LabelImageType::Pointer labelImage = nullptr;
	ReadImage<LabelImageType>(labelImage, fn.c_str());
	typename LabelImageType::SpacingType spacing = labelImage->GetSpacing();
	float volumeElement = 1.0;
	for (unsigned int i = 0; i < spacing.Size(); i++)
	{
	volumeElement *= static_cast<float>(spacing[i]);
	}
	typedef itk::LabelGeometryImageFilter<LabelImageType, RealImageType> GeometryFilterType;
	typename GeometryFilterType::Pointer geometryFilter = GeometryFilterType::New();
	geometryFilter->SetInput(labelImage);
	geometryFilter->CalculatePixelIndicesOff();
	geometryFilter->CalculateOrientedBoundingBoxOff();
	geometryFilter->CalculateOrientedLabelRegionsOff();
	geometryFilter->Update();
	typedef itk::LabelPerimeterEstimationCalculator<LabelImageType> AreaFilterType;
	typename AreaFilterType::Pointer areaFilter = AreaFilterType::New();
	areaFilter->SetImage(labelImage);
	areaFilter->SetFullyConnected(true);
	areaFilter->Compute();
	typename GeometryFilterType::LabelsType allLabels = geometryFilter->GetLabels();
	typename GeometryFilterType::LabelsType::iterator allLabelsIt;
	std::sort(allLabels.begin(), allLabels.end());
	for (allLabelsIt = allLabels.begin(); allLabelsIt != allLabels.end(); allLabelsIt++)
	{
	if (*allLabelsIt == 0)
	{
	continue;
	}
	// RealType volume = geometryFilter->GetVolume( *allLabelsIt ) * volumeElement;
	// RealType perimeter = areaFilter->GetPerimeter( *allLabelsIt );
	// RealType thicknessPrior = volume / perimeter;
	// std::cout << "Label " << *allLabelsIt << ": ";
	// std::cout << "volume = " << volume << ", ";
	// std::cout << "area = " << perimeter << ", ";
	// std::cout << "thickness = " << thicknessPrior << std::endl;
	}
	typename RealImageType::Pointer thicknessPriorImage = RealImageType::New();
	thicknessPriorImage->CopyInformation(labelImage);
	thicknessPriorImage->SetRegions(labelImage->GetLargestPossibleRegion());
	thicknessPriorImage->AllocateInitialized();
	itk::ImageRegionIteratorWithIndex<LabelImageType> It(labelImage, labelImage->GetLargestPossibleRegion());
	for (It.GoToBegin(); !It.IsAtEnd(); ++It)
	{
	LabelType label = It.Get();
	if (label == 0)
	{
	continue;
	}
	RealType volume = geometryFilter->GetVolume(label) * volumeElement;
	RealType perimeter = areaFilter->GetPerimeter(label);
	RealType thicknessPrior = volume / perimeter;
	thicknessPriorImage->SetPixel(It.GetIndex(), thicknessPrior);
	}
	ANTs::WriteImage<RealImageType>(thicknessPriorImage, outname.c_str());
	return EXIT_SUCCESS;
	}

[cookpa]

It seems that the pseudo-thickness image from ImageMath LabelThickness2 should be provided by GetConnectedComponentsFeatureImages - only there's a bug in the latter and it outputs area / volume, not volume / area.

[ntustison]

Thanks @cookpa . I'm pretty sure that the only use that had was around the time the 2014 cortical thickness paper was written. As far as I'm concerned, we can get rid of it.

[cookpa]

I talk about the changes a bit more here

InsightSoftwareConsortium/ITK#4624 (comment)

[cookpa]

OK I think we are good to go now, all uses of LabelGeometryImageFilter have been replaced.

[ntustison]

Awesome. Thanks @cookpa.