diff --git a/Examples/CreateDTICohort.cxx b/Examples/CreateDTICohort.cxx index 522f0c732..dbda673e0 100644 --- a/Examples/CreateDTICohort.cxx +++ b/Examples/CreateDTICohort.cxx @@ -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" @@ -222,29 +222,27 @@ CreateDTICohort(itk::ants::CommandLineParser * parser) // // Get label information for pathology option // - using LabelGeometryFilterType = itk::LabelGeometryImageFilter; - 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; + 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(labelGeometry->GetVolume(labels[n])); + totalMaskVolume += static_cast(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 pathologyParameters(labels.size(), 3); - for (unsigned int i = 0; i < labels.size(); i++) + itk::Array2D pathologyParameters(labelObjects.size(), 3); + for (unsigned int i = 0; i < labelObjects.size(); i++) { pathologyParameters(i, 0) = 0.0; pathologyParameters(i, 1) = 0.0; @@ -275,12 +273,14 @@ CreateDTICohort(itk::ants::CommandLineParser * parser) { percentageVoxels = parser->Convert(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::OneValue()) { - percentage /= static_cast(labelGeometry->GetVolume(labels[n])); + percentage /= static_cast(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; @@ -294,9 +294,17 @@ CreateDTICohort(itk::ants::CommandLineParser * parser) { auto whichClass = parser->Convert(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; } @@ -321,11 +329,13 @@ CreateDTICohort(itk::ants::CommandLineParser * parser) RealType percentage = percentageVoxels; if (percentage > itk::NumericTraits::OneValue()) { - percentage /= static_cast(labelGeometry->GetVolume(*it)); + percentage /= static_cast((*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; } } } @@ -539,7 +549,7 @@ CreateDTICohort(itk::ants::CommandLineParser * parser) // itksys::SystemTools::MakeDirectory(outputDirectory.c_str()); - itk::Array2D meanFAandMD(labels.size(), 5); + itk::Array2D meanFAandMD(labelObjects.size(), 5); meanFAandMD.Fill(0.0); for (unsigned n = 0; n <= numberOfControls + numberOfExperimentals; n++) { @@ -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; @@ -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) @@ -904,13 +924,13 @@ 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. "); @@ -918,8 +938,8 @@ InitializeCommandLineOptions(itk::ants::CommandLineParser * parser) OptionType::Pointer option = OptionType::New(); option->SetLongName("pathology"); option->SetUsageOption(0, - "label[,,]"); + "label[,,]"); option->SetShortName('p'); option->SetDescription(description); parser->AddOption(option); diff --git a/Examples/GetConnectedComponentsFeatureImages.cxx b/Examples/GetConnectedComponentsFeatureImages.cxx index 71badd16e..648ff76c2 100644 --- a/Examples/GetConnectedComponentsFeatureImages.cxx +++ b/Examples/GetConnectedComponentsFeatureImages.cxx @@ -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" @@ -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(spacing[d]); - } using RelabelerType = itk::RelabelComponentImageFilter; typename RelabelerType::Pointer relabeler = RelabelerType::New(); @@ -78,18 +72,13 @@ GetConnectedComponentsFeatureImages(int itkNotUsed(argc), char * argv[]) relabeler2->SetInput(filter->GetOutput()); relabeler2->Update(); - using GeometryFilterType = itk::LabelGeometryImageFilter; + using GeometryFilterType = itk::LabelImageToShapeLabelMapFilter; 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; - typename AreaFilterType::Pointer area = AreaFilterType::New(); - area->SetImage(relabeler2->GetOutput()); - area->Compute(); + geometry->Update(); itk::ImageRegionIteratorWithIndex It(relabeler->GetOutput(), relabeler->GetOutput()->GetRequestedRegion()); @@ -108,12 +97,32 @@ GetConnectedComponentsFeatureImages(int itkNotUsed(argc), char * argv[]) // [2] = eccentricity // [3] = elongation - float volume = prefactor * static_cast(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(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(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; + } } } } diff --git a/Examples/ImageMath.cxx b/Examples/ImageMath.cxx index 231637f0f..b3ad3a74f 100644 --- a/Examples/ImageMath.cxx +++ b/Examples/ImageMath.cxx @@ -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" diff --git a/Examples/ImageMath_Templates.hxx b/Examples/ImageMath_Templates.hxx index 6a6b83181..453902cd6 100644 --- a/Examples/ImageMath_Templates.hxx +++ b/Examples/ImageMath_Templates.hxx @@ -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" @@ -10779,93 +10777,6 @@ LabelThickness(int argc, char * argv[]) return EXIT_SUCCESS; } -template -int -LabelThickness2(int argc, char * argv[]) -{ - typedef unsigned int LabelType; - typedef itk::Image LabelImageType; - typedef float RealType; - typedef itk::Image 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(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(spacing[i]); - } - - typedef itk::LabelGeometryImageFilter GeometryFilterType; - typename GeometryFilterType::Pointer geometryFilter = GeometryFilterType::New(); - geometryFilter->SetInput(labelImage); - geometryFilter->CalculatePixelIndicesOff(); - geometryFilter->CalculateOrientedBoundingBoxOff(); - geometryFilter->CalculateOrientedLabelRegionsOff(); - geometryFilter->Update(); - - typedef itk::LabelPerimeterEstimationCalculator 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 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(thicknessPriorImage, outname.c_str()); - return EXIT_SUCCESS; -} - // now words can be accessed like this WordList[n]; where 'n' is the index template int @@ -14936,11 +14847,6 @@ ImageMathHelperAll(int argc, char ** argv) LabelThickness(argc, argv); return EXIT_SUCCESS; } - if (operation == "LabelThickness2") - { - LabelThickness2(argc, argv); - return EXIT_SUCCESS; - } if (operation == "DiceAndMinDistSum") { DiceAndMinDistSum(argc, argv); diff --git a/ImageSegmentation/antsAtroposSegmentationImageFilter.hxx b/ImageSegmentation/antsAtroposSegmentationImageFilter.hxx index fa5e07a14..ff8480bce 100644 --- a/ImageSegmentation/antsAtroposSegmentationImageFilter.hxx +++ b/ImageSegmentation/antsAtroposSegmentationImageFilter.hxx @@ -33,7 +33,7 @@ #include "itkImageRegionIteratorWithIndex.h" #include "itkIterationReporter.h" #include "itkKdTreeBasedKmeansEstimator.h" -#include "itkLabelGeometryImageFilter.h" +#include "itkLabelImageToShapeLabelMapFilter.h" #include "itkLabelStatisticsImageFilter.h" #include "itkMinimumDecisionRule.h" #include "itkMultiplyImageFilter.h" @@ -1306,14 +1306,18 @@ AtroposSegmentationImageFilter::Updat NumericTraits::OneValue() / static_cast(totalNumberOfClasses); } } + auto labelMapFilter = itk::LabelImageToShapeLabelMapFilter::New(); + labelMapFilter->SetInput(maxLabels); + labelMapFilter->SetComputeOrientedBoundingBox(false); + labelMapFilter->SetComputePerimeter(false); + labelMapFilter->SetComputeFeretDiameter(false); + labelMapFilter->SetInput(maxLabels); + labelMapFilter->Update(); - typedef LabelGeometryImageFilter GeometryType; - typename GeometryType::Pointer geom = GeometryType::New(); - geom->SetInput(maxLabels); - geom->Update(); for (unsigned int n = 0; n < totalNumberOfClasses; n++) { - this->m_LabelVolumes[n] = geom->GetVolume(n + 1); + auto labelObject = labelMapFilter->GetOutput()->GetNthLabelObject(n); + this->m_LabelVolumes[n] = labelObject->GetNumberOfPixels(); } this->SetNthOutput(0, maxLabels); diff --git a/Utilities/itkLabelPerimeterEstimationCalculator.h b/Utilities/itkLabelPerimeterEstimationCalculator.h deleted file mode 100644 index d0a52e63a..000000000 --- a/Utilities/itkLabelPerimeterEstimationCalculator.h +++ /dev/null @@ -1,144 +0,0 @@ -/*========================================================================= - - Program: Insight Segmentation & Registration Toolkit - Module: itkLabelPerimeterEstimationCalculator.h - Date: $Date$ - Version: $Revision$ - - Copyright (c) Insight Software Consortium. All rights reserved. - See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. - - This software is distributed WITHOUT ANY WARRANTY; without even - the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR - PURPOSE. See the above copyright notices for more information. - -=========================================================================*/ -#ifndef __itkLabelPerimeterEstimationCalculator_h -#define __itkLabelPerimeterEstimationCalculator_h - -#include "itkObject.h" - -namespace itk -{ -/** \class LabelPerimeterEstimationCalculator - * \brief Estimates the perimeter of a label object. - * - * The LabelPerimeterEstimationCalculator takes a label object and calculates - * an estimated perimeter based on pixels' and their neighbors. - * - * This implementation was taken from the Insight Journal paper: - * http://hdl.handle.net/1926/584 or - * http://www.insight-journal.org/browse/publication/176 - * - * \author Gaetan Lehmann. Biologie du Developpement et de la Reproduction, INRA de Jouy-en-Josas, France. - * - * \sa - */ -template -class LabelPerimeterEstimationCalculator final : public Object -{ -public: - /** Standard class typedefs. */ - typedef LabelPerimeterEstimationCalculator Self; - typedef Object Superclass; - typedef SmartPointer Pointer; - typedef SmartPointer ConstPointer; - - /** Some convenient typedefs. */ - typedef TInputImage InputImageType; - typedef typename InputImageType::Pointer InputImagePointer; - typedef typename InputImageType::ConstPointer InputImageConstPointer; - typedef typename InputImageType::RegionType InputImageRegionType; - typedef typename InputImageType::PixelType InputImagePixelType; - - typedef typename InputImageType::RegionType RegionType; - typedef typename InputImageType::SizeType SizeType; - typedef typename InputImageType::IndexType IndexType; - - typedef typename std::map PerimetersType; - - /** ImageDimension constants */ - static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; - - /** Standard New method. */ - itkNewMacro(Self); - - /** Runtime information support. */ - itkTypeMacro(LabelPerimeterEstimationCalculator, Object); - - /** - * Set/Get whether the connected components are defined strictly by - * face connectivity or by face+edge+vertex connectivity. Default is - * FullyConnectedOff. For objects that are 1 pixel wide, use - * FullyConnectedOn. - */ - itkSetMacro(FullyConnected, bool); - itkGetConstReferenceMacro(FullyConnected, bool); - itkBooleanMacro(FullyConnected); - - void - SetImage(const InputImageType * img) - { - m_Image = img; - } - - const InputImageType * - GetImage() const - { - return m_Image; - } - - void - Compute(); - - const PerimetersType & - GetPerimeters() const - { - return m_Perimeters; - } - - const double & - GetPerimeter(const InputImagePixelType & label) const - { - if (m_Perimeters.find(label) != m_Perimeters.end()) - { - return m_Perimeters.find(label)->second; - } - itkExceptionMacro(<< "Unknown label: " - << static_cast::PrintType>(label)); - } - - bool - HasLabel(const InputImagePixelType & label) const - { - if (m_Perimeters.find(label) != m_Perimeters.end()) - { - return true; - } - return false; - } - -protected: - LabelPerimeterEstimationCalculator(); - ~LabelPerimeterEstimationCalculator() override = default; - void - PrintSelf(std::ostream & os, Indent indent) const override; - -private: - LabelPerimeterEstimationCalculator(const Self &) = delete; - void - operator=(const Self &) = delete; - - bool m_FullyConnected; - - const InputImageType * m_Image; - - PerimetersType m_Perimeters; -}; // end of class -} // end namespace itk - -#ifndef ITK_MANUAL_INSTANTIATION -# include "itkLabelPerimeterEstimationCalculator.hxx" -#endif - -#endif diff --git a/Utilities/itkLabelPerimeterEstimationCalculator.hxx b/Utilities/itkLabelPerimeterEstimationCalculator.hxx deleted file mode 100644 index 58798fbfc..000000000 --- a/Utilities/itkLabelPerimeterEstimationCalculator.hxx +++ /dev/null @@ -1,201 +0,0 @@ -/*========================================================================= - - Program: Insight Segmentation & Registration Toolkit - Module: itkLabelPerimeterEstimationCalculator.txx - Date: $Date$ - Version: $Revision$ - - Copyright (c) Insight Software Consortium. All rights reserved. - See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. - - This software is distributed WITHOUT ANY WARRANTY; without even - the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR - PURPOSE. See the above copyright notices for more information. - -=========================================================================*/ -#ifndef __itkLabelPerimeterEstimationCalculator_hxx -#define __itkLabelPerimeterEstimationCalculator_hxx - -#include "itkProgressReporter.h" -#include "itkImageRegionIterator.h" -#include "itkImageRegionConstIteratorWithIndex.h" -#include "itkConstShapedNeighborhoodIterator.h" -#include "itkConstantBoundaryCondition.h" -#include "itkSize.h" -#include "itkConnectedComponentAlgorithm.h" -#include - -namespace itk -{ -template -LabelPerimeterEstimationCalculator::LabelPerimeterEstimationCalculator() -{ - m_FullyConnected = false; -} - -template -void -LabelPerimeterEstimationCalculator::Compute() -{ - m_Perimeters.clear(); - - // ProgressReporter progress( this, 0, - // this->GetImage()->GetRequestedRegion().GetNumberOfPixels() ); - - // reduce the region to avoid reading outside - RegionType region = this->GetImage()->GetRequestedRegion(); - SizeType size = region.GetSize(); - for (unsigned int i = 0; i < ImageDimension; i++) - { - size[i]--; - } - region.SetSize(size); - - // the radius which will be used for all the shaped iterators - Size radius; - radius.Fill(1); - - // set up the iterator - typedef ConstShapedNeighborhoodIterator IteratorType; - typename IteratorType::ConstIterator nIt; - IteratorType iIt(radius, this->GetImage(), region); - // we want to search the neighbors with offset >= 0 - // 2D -> 4 neighbors - // 3D -> 8 neighbors - typename IteratorType::OffsetType offset; - unsigned int centerIndex = iIt.GetCenterNeighborhoodIndex(); - // store the offsets to reuse them to evaluate the contributions of the - // configurations - typename std::vector indexes; - IndexType idx0; - idx0.Fill(0); - for (unsigned int d = centerIndex; d < 2 * centerIndex + 1; d++) - { - offset = iIt.GetOffset(d); - bool deactivate = false; - for (unsigned int j = 0; j < ImageDimension && !deactivate; j++) - { - if (offset[j] < 0) - { - deactivate = true; - } - } - if (deactivate) - { - iIt.DeactivateOffset(offset); - } - else - { - iIt.ActivateOffset(offset); - indexes.push_back(idx0 + offset); - } - } - - // to store the configurations count for all the labels - typedef typename std::map MapType; - typedef typename std::map LabelMapType; - LabelMapType confCount; - - for (iIt.GoToBegin(); !iIt.IsAtEnd(); ++iIt) - { - // 2 pass - find the labels in the neighborhood - // - count the configurations for all the labels - - typedef typename std::set LabelSetType; - LabelSetType labelSet; - for (nIt = iIt.Begin(); nIt != iIt.End(); nIt++) - { - labelSet.insert(nIt.Get()); - } - - for (typename LabelSetType::const_iterator it = labelSet.begin(); it != labelSet.end(); it++) - { - unsigned long conf = 0; - int i = 0; - - for (nIt = iIt.Begin(); nIt != iIt.End(); nIt++, i++) - { - if (nIt.Get() == *it) - { - conf += 1 << i; - } - } - - confCount[*it][conf]++; - } - - // progress.CompletedPixel(); - } - - // compute the participation to the perimeter for all the configurations - double physicalSize = 1.0; - for (unsigned int i = 0; i < ImageDimension; i++) - { - physicalSize *= this->GetImage()->GetSpacing()[i]; - } - typedef typename std::map ContributionMapType; - ContributionMapType contributions; - const unsigned int numberOfNeighbors = static_cast(std::pow(2.0, static_cast(ImageDimension))); - const unsigned int numberOfConfigurations = - static_cast(std::pow(2.0, static_cast(numberOfNeighbors))); - // create an image to store the neighbors - typedef typename itk::Image ImageType; - typename ImageType::Pointer neighborsImage = ImageType::New(); - // typename ImageType::SizeType size; - size.Fill(2); - neighborsImage->SetRegions(size); - neighborsImage->Allocate(); - for (unsigned int i = 0; i < numberOfConfigurations; i++) - { - neighborsImage->FillBuffer(false); - for (unsigned int j = 0; j < numberOfNeighbors; j++) - { - if (i & 1 << j) - { - neighborsImage->SetPixel(indexes[j], true); - } - } - // the image is created - we can now compute the contributions of the pixels - // for that configuration - contributions[i] = 0; - for (unsigned int j = 0; j < numberOfNeighbors; j++) - { - IndexType currentIdx = indexes[j]; - if (neighborsImage->GetPixel(currentIdx)) - { - for (unsigned int k = 0; k < ImageDimension; k++) - { - IndexType idx = currentIdx; - idx[k] = std::abs(idx[k] - 1); - if (!neighborsImage->GetPixel(idx)) - { - contributions[i] += physicalSize / this->GetImage()->GetSpacing()[k] / 2.0; - } - } - } - } - contributions[i] /= ImageDimension; - } - - // and use those contributions to found the perimeter - m_Perimeters.clear(); - for (typename LabelMapType::const_iterator it = confCount.begin(); it != confCount.end(); it++) - { - m_Perimeters[it->first] = 0; - for (typename MapType::const_iterator it2 = it->second.begin(); it2 != it->second.end(); it2++) - { - m_Perimeters[it->first] += contributions[it2->first] * it2->second; - } - } -} - -template -void -LabelPerimeterEstimationCalculator::PrintSelf(std::ostream & os, Indent indent) const -{ - Superclass::PrintSelf(os, indent); - - os << indent << "FullyConnected: " << m_FullyConnected << std::endl; -} -} // end namespace itk -#endif