Skip to content

Commit

Permalink
Merge pull request #1741 from ANTsX/remove_deprecated_label_geometry
Browse files Browse the repository at this point in the history 
Remove deprecated label geometry filter
  • Loading branch information
@cookpa
cookpa authored May 7, 2024
2 parents 4c68a48 + a89a499 commit 3a788f8
Show file tree
Hide file tree
Showing 7 changed files with 96 additions and 503 deletions.
92 changes: 56 additions & 36 deletions Examples/CreateDTICohort.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
#include "itkImageFileWriter.h"
#include "itkImageRegionConstIterator.h"
#include "itkImageRegionIterator.h"
#include "itkLabelGeometryImageFilter.h"
#include "itkLabelImageToShapeLabelMapFilter.h"
#include "itkMersenneTwisterRandomVariateGenerator.h"
#include "itkNumericSeriesFileNames.h"
#include "itkTimeProbe.h"
Expand Down Expand Up @@ -222,29 +222,27 @@ CreateDTICohort(itk::ants::CommandLineParser * parser)
//
// Get label information for pathology option
//
using LabelGeometryFilterType = itk::LabelGeometryImageFilter<MaskImageType, ImageType>;
typename LabelGeometryFilterType::Pointer labelGeometry = LabelGeometryFilterType::New();
labelGeometry->SetInput(maskImage);
labelGeometry->CalculatePixelIndicesOff();
labelGeometry->CalculateOrientedBoundingBoxOff();
labelGeometry->CalculateOrientedLabelRegionsOff();
labelGeometry->Update();

typename LabelGeometryFilterType::LabelsType labels = labelGeometry->GetLabels();
std::sort(labels.begin(), labels.end());
using LabelMapFilterType = itk::LabelImageToShapeLabelMapFilter<MaskImageType>;
typename LabelMapFilterType::Pointer labelMapper = LabelMapFilterType::New();
labelMapper->SetInput(maskImage);
labelMapper->ComputeOrientedBoundingBoxOff();
labelMapper->ComputePerimeterOff();
labelMapper->SetBackgroundValue(-1); // include 0 in output
labelMapper->Update();
auto labelObjects = labelMapper->GetOutput()->GetLabelObjects();

unsigned int totalMaskVolume = 0;
for (unsigned int n = 0; n < labels.size(); n++)
for (unsigned int n = 0; n < labelObjects.size(); n++)
{
totalMaskVolume += static_cast<unsigned int>(labelGeometry->GetVolume(labels[n]));
totalMaskVolume += static_cast<unsigned int>(labelObjects[n]->GetNumberOfPixels());
}

// Fill in default values per label:
// column 1: percentage change in longitudinal eigenvalue
// column 2: percentage change in average of transverse eigenvalue(s)
// column 3: percentage of affected voxels
itk::Array2D<RealType> pathologyParameters(labels.size(), 3);
for (unsigned int i = 0; i < labels.size(); i++)
itk::Array2D<RealType> pathologyParameters(labelObjects.size(), 3);
for (unsigned int i = 0; i < labelObjects.size(); i++)
{
pathologyParameters(i, 0) = 0.0;
pathologyParameters(i, 1) = 0.0;
Expand Down Expand Up @@ -275,12 +273,14 @@ CreateDTICohort(itk::ants::CommandLineParser * parser)
{
percentageVoxels = parser->Convert<float>(pathologyOption->GetFunction(0)->GetParameter(2));
}
for (unsigned int n = 0; n < labels.size(); n++)
for (unsigned int n = 0; n < labelObjects.size(); n++)
{
RealType percentage = percentageVoxels;
if (percentage > itk::NumericTraits<RealType>::OneValue())
{
percentage /= static_cast<RealType>(labelGeometry->GetVolume(labels[n]));
percentage /= static_cast<RealType>(labelObjects[n]->GetNumberOfPixels());
std::cout << "WARNING: Percentage of affected voxels must be less than or equal to 1." << std::endl;
std::cout << " Dividing by size of region, percentage = " << percentage << std::endl;
}

pathologyParameters(n, 0) = pathologyDeltaEig1;
Expand All @@ -294,9 +294,17 @@ CreateDTICohort(itk::ants::CommandLineParser * parser)
{
auto whichClass = parser->Convert<LabelType>(pathologyOption->GetFunction(n)->GetName());

auto it = std::find(labels.begin(), labels.end(), whichClass);
auto it = labelObjects.begin();

if (it == labels.end())
for (it = labelObjects.begin(); it != labelObjects.end(); ++it)
{
if ((*it)->GetLabel() == whichClass)
{
break;
}
}

if (it == labelObjects.end())
{
continue;
}
Expand All @@ -321,11 +329,13 @@ CreateDTICohort(itk::ants::CommandLineParser * parser)
RealType percentage = percentageVoxels;
if (percentage > itk::NumericTraits<RealType>::OneValue())
{
percentage /= static_cast<RealType>(labelGeometry->GetVolume(*it));
percentage /= static_cast<RealType>((*it)->GetNumberOfPixels());
std::cout << "WARNING: Percentage of affected voxels must be less than or equal to 1." << std::endl;
std::cout << " Dividing by size of region, percentage = " << percentage << std::endl;
}
pathologyParameters(it - labels.begin(), 0) = pathologyDeltaEig1;
pathologyParameters(it - labels.begin(), 1) = pathologyDeltaEig2_Eig3;
pathologyParameters(it - labels.begin(), 2) = percentage;
pathologyParameters(it - labelObjects.begin(), 0) = pathologyDeltaEig1;
pathologyParameters(it - labelObjects.begin(), 1) = pathologyDeltaEig2_Eig3;
pathologyParameters(it - labelObjects.begin(), 2) = percentage;
}
}
}
Expand Down Expand Up @@ -539,7 +549,7 @@ CreateDTICohort(itk::ants::CommandLineParser * parser)
//
itksys::SystemTools::MakeDirectory(outputDirectory.c_str());

itk::Array2D<RealType> meanFAandMD(labels.size(), 5);
itk::Array2D<RealType> meanFAandMD(labelObjects.size(), 5);
meanFAandMD.Fill(0.0);
for (unsigned n = 0; n <= numberOfControls + numberOfExperimentals; n++)
{
Expand Down Expand Up @@ -609,12 +619,22 @@ CreateDTICohort(itk::ants::CommandLineParser * parser)
eigenvalues[2] = eigenvalues[1];
}

auto it = std::find(labels.begin(), labels.end(), label);
if (it == labels.end())
auto it = labelObjects.begin();

for (it = labelObjects.begin(); it != labelObjects.end(); ++it)
{
if ((*it)->GetLabel() == label)
{
break;
}
}

if (it == labelObjects.end())
{
std::cout << "ERROR: unknown label." << std::endl;
}
unsigned int labelIndex = it - labels.begin();

unsigned int labelIndex = it - labelObjects.begin();

typename TensorType::EigenValuesArrayType newEigenvalues;

Expand Down Expand Up @@ -726,12 +746,12 @@ CreateDTICohort(itk::ants::CommandLineParser * parser)
if (n == 0)
{
std::cout << " " << std::left << std::setw(7) << "Region" << std::left << std::setw(15) << "FA (original)"
<< std::left << std::setw(15) << "FA (path+isv)" << std::left << std::setw(15) << "FA (% change)"
<< std::left << std::setw(15) << "FA (path+isv)" << std::left << std::setw(15) << "FA (prop. change)"
<< std::left << std::setw(15) << "MD (original)" << std::left << std::setw(15) << "MD (path+isv)"
<< std::left << std::setw(15) << "MD (% change)" << std::endl;
for (unsigned int l = 1; l < labels.size(); l++)
<< std::left << std::setw(15) << "MD (prop. change)" << std::endl;
for (unsigned int l = 1; l < labelObjects.size(); l++)
{
std::cout << " " << std::left << std::setw(7) << labels[l] << std::left << std::setw(15)
std::cout << " " << std::left << std::setw(7) << labelObjects[l]->GetLabel() << std::left << std::setw(15)
<< meanFAandMD(l, 0) / meanFAandMD(l, 4) << std::left << std::setw(15)
<< meanFAandMD(l, 2) / meanFAandMD(l, 4) << std::left << std::setw(15)
<< (meanFAandMD(l, 2) - meanFAandMD(l, 0)) / meanFAandMD(l, 0) << std::left << std::setw(15)
Expand Down Expand Up @@ -904,22 +924,22 @@ InitializeCommandLineOptions(itk::ants::CommandLineParser * parser)
std::string("Pathology is simulated by changing the eigenvalues. Typically ") +
std::string("this involves a decrease in the largest eigenvalue and an ") +
std::string("increase in the average of the remaining eigenvalues. ") +
std::string("Change is specified as a percentage of the current eigenvalues. ") +
std::string("Change is specified as a proportion of the current eigenvalues. ") +
std::string("However, care is taken ") +
std::string("to ensure that diffusion direction does not change. ") +
std::string("Additionally, one can specify the number of voxels affected ") +
std::string("in each region or one can specify the percentage of voxels ") +
std::string("in each region or one can specify the proportion of voxels ") +
std::string("affected. Default is to change all voxels. Note that the ") +
std::string("percentages must be specified in the range [0,1]. For ") +
std::string("proportions must be specified in the range [0,1]. For ") +
std::string("dimension=3 where the average transverse diffusion eigenvalues ") +
std::string("are altered, this change is propagated to the distinct eigenvalues ") +
std::string("by forcing the ratio to be the same before the change. ");

OptionType::Pointer option = OptionType::New();
option->SetLongName("pathology");
option->SetUsageOption(0,
"label[<percentageChangeEig1=-0.05>,<percentageChangeAvgEig2andEig3=0.05>,<numberOfVoxels="
"all or percentageOfvoxels>]");
"label[<propChangeEig1=-0.05>,<propChangeAvgEig2andEig3=0.05>,<numberOfVoxels="
"all or propOfVoxels>]");
option->SetShortName('p');
option->SetDescription(description);
parser->AddOption(option);
Expand Down
51 changes: 30 additions & 21 deletions Examples/GetConnectedComponentsFeatureImages.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,7 @@
#include "itkBinaryThresholdImageFilter.h"
#include "itkConnectedComponentImageFilter.h"
#include "itkRelabelComponentImageFilter.h"
#include "itkLabelGeometryImageFilter.h"
#include "itkLabelPerimeterEstimationCalculator.h"
#include "itkLabelImageToShapeLabelMapFilter.h"
#include "itkLabelStatisticsImageFilter.h"
#include "itkMultiplyImageFilter.h"
#include "itkSignedMaurerDistanceMapImageFilter.h"
Expand Down Expand Up @@ -47,11 +46,6 @@ GetConnectedComponentsFeatureImages(int itkNotUsed(argc), char * argv[])
}

typename ImageType::SpacingType spacing = inputImage->GetSpacing();
float prefactor = 1.0;
for (unsigned int d = 0; d < ImageDimension; d++)
{
prefactor *= static_cast<float>(spacing[d]);
}

using RelabelerType = itk::RelabelComponentImageFilter<ImageType, ImageType>;
typename RelabelerType::Pointer relabeler = RelabelerType::New();
Expand All @@ -78,18 +72,13 @@ GetConnectedComponentsFeatureImages(int itkNotUsed(argc), char * argv[])
relabeler2->SetInput(filter->GetOutput());
relabeler2->Update();

using GeometryFilterType = itk::LabelGeometryImageFilter<ImageType, RealImageType>;
using GeometryFilterType = itk::LabelImageToShapeLabelMapFilter<ImageType>;
typename GeometryFilterType::Pointer geometry = GeometryFilterType::New();
geometry->SetInput(relabeler2->GetOutput());
geometry->CalculatePixelIndicesOff();
geometry->CalculateOrientedBoundingBoxOff();
geometry->CalculateOrientedLabelRegionsOff();
geometry->Update();
geometry->ComputeOrientedBoundingBoxOff();
geometry->ComputePerimeterOn();

using AreaFilterType = itk::LabelPerimeterEstimationCalculator<ImageType>;
typename AreaFilterType::Pointer area = AreaFilterType::New();
area->SetImage(relabeler2->GetOutput());
area->Compute();
geometry->Update();

itk::ImageRegionIteratorWithIndex<ImageType> It(relabeler->GetOutput(),
relabeler->GetOutput()->GetRequestedRegion());
Expand All @@ -108,12 +97,32 @@ GetConnectedComponentsFeatureImages(int itkNotUsed(argc), char * argv[])
// [2] = eccentricity
// [3] = elongation

float volume = prefactor * static_cast<float>(geometry->GetVolume(label));
try
{
auto labelObject = geometry->GetOutput()->GetLabelObject(label);
float volume = labelObject->GetPhysicalSize();

outputImages[0]->SetPixel(index, volume);
outputImages[1]->SetPixel(index, volume / static_cast<float>(labelObject->GetPerimeter()));

auto principalMoments = labelObject->GetPrincipalMoments();

float lambda1 = principalMoments[0];
float lambdaN = principalMoments[ImageDimension - 1];
float eccentricity = 0.0;

if (!itk::Math::FloatAlmostEqual(lambda1, 0.0f))
{
eccentricity = std::sqrt(1.0 - (lambda1 * lambda1) / (lambdaN * lambdaN));
}

outputImages[0]->SetPixel(index, volume);
outputImages[1]->SetPixel(index, static_cast<float>(area->GetPerimeter(label)) / volume);
outputImages[2]->SetPixel(index, geometry->GetEccentricity(label));
outputImages[3]->SetPixel(index, geometry->GetElongation(label));
outputImages[2]->SetPixel(index, eccentricity);
outputImages[3]->SetPixel(index, labelObject->GetElongation());
}
catch (itk::ExceptionObject & e)
{
std::cerr << "Could not find label " << label << std::endl;
}
}
}
}
Expand Down
1 change: 0 additions & 1 deletion Examples/ImageMath.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,6 @@
#include "itkImageRegionIteratorWithIndex.h"
#include "itkLabelGeometryImageFilter.h"
#include "itkLabelOverlapMeasuresImageFilter.h"
#include "itkLabelPerimeterEstimationCalculator.h"
#include "itkKdTree.h"
#include "itkKdTreeBasedKmeansEstimator.h"
#include "itkLabelContourImageFilter.h"
Expand Down
94 changes: 0 additions & 94 deletions Examples/ImageMath_Templates.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,9 +66,7 @@
#include "itkImageRandomConstIteratorWithIndex.h"
#include "itkImageRegionIterator.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "itkLabelGeometryImageFilter.h"
#include "itkLabelOverlapMeasuresImageFilter.h"
#include "itkLabelPerimeterEstimationCalculator.h"
#include "itkKdTree.h"
#include "itkKdTreeBasedKmeansEstimator.h"
#include "itkLabelContourImageFilter.h"
Expand Down Expand Up @@ -10779,93 +10777,6 @@ LabelThickness(int argc, char * argv[])
return EXIT_SUCCESS;
}

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;
}

// now words can be accessed like this WordList[n]; where 'n' is the index
template <unsigned int ImageDimension>
int
Expand Down Expand Up @@ -14936,11 +14847,6 @@ ImageMathHelperAll(int argc, char ** argv)
LabelThickness<DIM>(argc, argv);
return EXIT_SUCCESS;
}
if (operation == "LabelThickness2")
{
LabelThickness2<DIM>(argc, argv);
return EXIT_SUCCESS;
}
if (operation == "DiceAndMinDistSum")
{
DiceAndMinDistSum<DIM>(argc, argv);
Expand Down
Loading

0 comments on commit 3a788f8

Please sign in to comment.