X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=segmentation%2FclitkExtractLungFilter.txx;h=b7b28665baff5006bc2de75cfc158fd8be9821fe;hb=464e9a7e4aa4d79ecfd40c3b4b32748f0d8e5d51;hp=55e9e6cbcefc0efa015a2f6047a886917f4904f5;hpb=e0464c9f7bdde185755d60e9dead188291fc3550;p=clitk.git diff --git a/segmentation/clitkExtractLungFilter.txx b/segmentation/clitkExtractLungFilter.txx index 55e9e6c..b7b2866 100644 --- a/segmentation/clitkExtractLungFilter.txx +++ b/segmentation/clitkExtractLungFilter.txx @@ -3,7 +3,7 @@ Authors belong to: - University of LYON http://www.universite-lyon.fr/ - - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr + - Léon Bérard cancer center http://www.centreleonberard.fr - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr This software is distributed WITHOUT ANY WARRANTY; without even @@ -14,7 +14,7 @@ - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html - ======================================================================-====*/ + ===========================================================================**/ #ifndef CLITKEXTRACTLUNGSFILTER_TXX #define CLITKEXTRACTLUNGSFILTER_TXX @@ -24,6 +24,9 @@ #include "clitkSetBackgroundImageFilter.h" #include "clitkSegmentationUtils.h" #include "clitkAutoCropFilter.h" +#include "clitkCropLikeImageFilter.h" +#include "clitkFillMaskFilter.h" +#include "clitkMemoryUsage.h" // itk #include "itkBinaryThresholdImageFilter.h" @@ -32,24 +35,38 @@ #include "itkOtsuThresholdImageFilter.h" #include "itkBinaryThinningImageFilter3D.h" #include "itkImageIteratorWithIndex.h" - +#include "itkBinaryMorphologicalOpeningImageFilter.h" +#include "itkBinaryMorphologicalClosingImageFilter.h" +#include "itkConstantPadImageFilter.h" +#include +#include "itkStatisticsLabelObject.h" +#include "itkLabelMap.h" +#include "itkLabelImageToShapeLabelMapFilter.h" +#include "itkLabelMapToLabelImageFilter.h" +#include "itkExtractImageFilter.h" +#include "itkOrientImageFilter.h" +#include "itkSpatialOrientationAdapter.h" +#include "itkImageDuplicator.h" +#include //-------------------------------------------------------------------- -template -clitk::ExtractLungFilter:: +template +clitk::ExtractLungFilter:: ExtractLungFilter(): clitk::FilterBase(), + clitk::FilterWithAnatomicalFeatureDatabaseManagement(), itk::ImageToImageFilter() { SetNumberOfSteps(10); m_MaxSeedNumber = 500; // Default global options - this->SetNumberOfRequiredInputs(2); + this->SetNumberOfRequiredInputs(1); SetPatientMaskBackgroundValue(0); SetBackgroundValue(0); // Must be zero SetForegroundValue(1); SetMinimalComponentSize(100); + VerboseRegionGrowingFlagOff(); // Step 1 default values SetUpperThreshold(-300); @@ -62,9 +79,16 @@ ExtractLungFilter(): SetLabelizeParameters1(p1); // Step 2 default values - SetUpperThresholdForTrachea(-900); + SetTracheaSeedAlgorithm(0); + SetUpperThresholdForTrachea(-700); SetMultiplierForTrachea(5); SetThresholdStepSizeForTrachea(64); + SetNumberOfSlicesToSkipBeforeSearchingSeed(0); + TracheaVolumeMustBeCheckedFlagOn(); + SetNumSlices(50); + SetMaxElongation(0.5); + SetSeedPreProcessingThreshold(-400); + // Step 3 default values SetNumberOfHistogramBins(500); @@ -82,16 +106,21 @@ ExtractLungFilter(): p3->UseLastKeepOff(); SetLabelizeParameters3(p3); - // Step 5 : find bronchial bifurcations - FindBronchialBifurcationsOn(); + // Step 5 + OpenCloseFlagOff(); + SetOpenCloseRadius(1); + + // Step 6 + FillHolesFlagOn(); + AutoCropOn(); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -template +template void -clitk::ExtractLungFilter:: +clitk::ExtractLungFilter:: SetInput(const ImageType * image) { this->SetNthInput(0, const_cast(image)); @@ -100,21 +129,9 @@ SetInput(const ImageType * image) //-------------------------------------------------------------------- -template +template void -clitk::ExtractLungFilter:: -SetInputPatientMask(MaskImageType * image, MaskImagePixelType bg ) -{ - this->SetNthInput(1, const_cast(image)); - SetPatientMaskBackgroundValue(bg); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void -clitk::ExtractLungFilter:: +clitk::ExtractLungFilter:: AddSeed(InternalIndexType s) { m_Seeds.push_back(s); @@ -123,178 +140,149 @@ AddSeed(InternalIndexType s) //-------------------------------------------------------------------- -template -template +template void -clitk::ExtractLungFilter:: -SetArgsInfo(ArgsInfoType mArgsInfo) -{ - SetVerboseOption_GGO(mArgsInfo); - SetVerboseStep_GGO(mArgsInfo); - SetWriteStep_GGO(mArgsInfo); - SetVerboseWarningOff_GGO(mArgsInfo); - - SetUpperThreshold_GGO(mArgsInfo); - SetLowerThreshold_GGO(mArgsInfo); - SetLabelizeParameters1_GGO(mArgsInfo); - if (!mArgsInfo.remove1_given) { - GetLabelizeParameters1()->AddLabelToRemove(2); - if (GetVerboseOption()) VerboseOption("remove1", 2); - } - - SetUpperThresholdForTrachea_GGO(mArgsInfo); - SetMultiplierForTrachea_GGO(mArgsInfo); - SetThresholdStepSizeForTrachea_GGO(mArgsInfo); - AddSeed_GGO(mArgsInfo); - - SetMinimalComponentSize_GGO(mArgsInfo); - - SetNumberOfHistogramBins_GGO(mArgsInfo); - SetLabelizeParameters2_GGO(mArgsInfo); - - SetRadiusForTrachea_GGO(mArgsInfo); - SetLabelizeParameters3_GGO(mArgsInfo); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void -clitk::ExtractLungFilter:: +clitk::ExtractLungFilter:: GenerateOutputInformation() { - Superclass::GenerateOutputInformation(); // Needed ?? - this->GetOutput(0)->SetRequestedRegion(this->GetOutput(0)->GetLargestPossibleRegion()); + clitk::PrintMemory(GetVerboseMemoryFlag(), "Initial memory"); // OK + Superclass::GenerateOutputInformation(); + + // Read DB + LoadAFDB(); // Get input pointers - patient = dynamic_cast(itk::ProcessObject::GetInput(1)); input = dynamic_cast(itk::ProcessObject::GetInput(0)); + patient = GetAFDB()->template GetImage ("Patient"); + PrintMemory(GetVerboseMemoryFlag(), "After reading patient"); // OK - // Check image - if (!HaveSameSizeAndSpacing(input, patient)) { - this->SetLastError("* ERROR * the images (input and patient mask) must have the same size & spacing"); - return; - } - //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStepOrStop("Set background to initial image"); + // Crop input like patient image (must have the same spacing) + // It copy the input if the same are the same + StartNewStep("Copy and crop input image to 'patient' extends"); + typedef clitk::CropLikeImageFilter CropImageFilter; + typename CropImageFilter::Pointer cropFilter = CropImageFilter::New(); + // cropFilter->ReleaseDataFlagOn(); // NO=seg fault !! + cropFilter->SetInput(input); + cropFilter->SetCropLikeImage(patient); + cropFilter->Update(); + working_input = cropFilter->GetOutput(); + // cropFilter->Delete(); // NO !!!! if yes, sg fault, Cropfilter is buggy !?? + StopCurrentStep(working_input); + PrintMemory(GetVerboseMemoryFlag(), "After crop"); // OK, slightly more than a copy of input + + //-------------------------------------------------------------------- + //-------------------------------------------------------------------- + StartNewStep("Set background to initial image"); working_input = SetBackground - (input, patient, GetPatientMaskBackgroundValue(), -1000); + (working_input, patient, GetPatientMaskBackgroundValue(), -1000, true); + + // Pad images with air to prevent patient touching the image border + static const unsigned int Dim = ImageType::ImageDimension; + typedef itk::ConstantPadImageFilter PadFilterType; + typename PadFilterType::Pointer padFilter = PadFilterType::New(); + padFilter->SetInput(working_input); + padFilter->SetConstant(-1000); + typename ImageType::SizeType bounds; + for (unsigned i = 0; i < Dim - 1; ++i) + bounds[i] = 1; + bounds[Dim - 1] = 0; + padFilter->SetPadLowerBound(bounds); + padFilter->SetPadUpperBound(bounds); + padFilter->Update(); + working_input = padFilter->GetOutput(); + StopCurrentStep(working_input); + PrintMemory(GetVerboseMemoryFlag(), "After set bg"); // OK, additional mem = 0 + + /* + // We do not need patient mask anymore, release memory + GetAFDB()->template ReleaseImage("Patient"); + DD(patient->GetReferenceCount()); + PrintMemory(GetVerboseMemoryFlag(), "After delete patient"); // OK, additional mem = 0 + patient->Delete(); + PrintMemory(GetVerboseMemoryFlag(), "After delete patient"); // OK, additional mem = 0 + */ //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStepOrStop("Remove Air"); + StartNewStep("Remove Air"); + // Check threshold + if (m_UseLowerThreshold) { + if (m_LowerThreshold > m_UpperThreshold) { + clitkExceptionMacro("lower threshold cannot be greater than upper threshold."); + } + } // Threshold to get air - typedef itk::BinaryThresholdImageFilter InputBinarizeFilterType; - typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New(); + typedef itk::BinaryThresholdImageFilter InputBinarizeFilterType; + typename InputBinarizeFilterType::Pointer binarizeFilter = InputBinarizeFilterType::New(); binarizeFilter->SetInput(working_input); if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(m_LowerThreshold); + // binarizeFilter->CanRunInPlace() is false binarizeFilter->SetUpperThreshold(m_UpperThreshold); - binarizeFilter ->SetInsideValue(this->GetForegroundValue()); - binarizeFilter ->SetOutsideValue(this->GetBackgroundValue()); + binarizeFilter->SetInsideValue(this->GetForegroundValue()); + binarizeFilter->SetOutsideValue(this->GetBackgroundValue()); binarizeFilter->Update(); - working_image = binarizeFilter->GetOutput(); - - // Labelize and keep right labels - working_image = Labelize(working_image, GetBackgroundValue(), true, GetMinimalComponentSize()); + working_mask = binarizeFilter->GetOutput(); + PrintMemory(GetVerboseMemoryFlag(), "After Binarizefilter"); // OK, additional mem is one mask image - working_image = RemoveLabels - (working_image, GetBackgroundValue(), GetLabelizeParameters1()->GetLabelsToRemove()); + // Labelize and keep right labels + working_mask = + Labelize + (working_mask, GetBackgroundValue(), true, GetMinimalComponentSize()); + PrintMemory(GetVerboseMemoryFlag(), "After Labelize"); // BUG ? additional mem around 1 time the input ? + + working_mask = RemoveLabels + (working_mask, GetBackgroundValue(), GetLabelizeParameters1()->GetLabelsToRemove()); + PrintMemory(GetVerboseMemoryFlag(), "After RemoveLabels"); // OK additional mem = 0 - typename InternalImageType::Pointer air = KeepLabels - (working_image, + working_mask = KeepLabels + (working_mask, GetBackgroundValue(), GetForegroundValue(), GetLabelizeParameters1()->GetFirstKeep(), GetLabelizeParameters1()->GetLastKeep(), GetLabelizeParameters1()->GetUseLastKeep()); + PrintMemory(GetVerboseMemoryFlag(), "After KeepLabels to create the 'air'"); // Set Air to BG - working_input = SetBackground - (working_input, air, this->GetForegroundValue(), this->GetBackgroundValue()); + working_input = SetBackground + (working_input, working_mask, this->GetForegroundValue(), this->GetBackgroundValue(), true); + PrintMemory(GetVerboseMemoryFlag(), "After SetBackground"); // End StopCurrentStep(working_input); //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStepOrStop("Find the trachea"); - if (m_Seeds.size() == 0) { // try to find seed - // Search seed (parameters = UpperThresholdForTrachea) - static const unsigned int Dim = ImageType::ImageDimension; - typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion(); - typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize(); - typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex(); - sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-5; - sliceRegionSize[Dim-1]=5; - sliceRegion.SetSize(sliceRegionSize); - sliceRegion.SetIndex(sliceRegionIndex); - typedef itk::ImageRegionConstIterator IteratorType; - IteratorType it(working_input, sliceRegion); - it.GoToBegin(); - while (!it.IsAtEnd()) { - if(it.Get() < GetUpperThresholdForTrachea() ) { - AddSeed(it.GetIndex()); - } - ++it; - } - } - - if (m_Seeds.size() != 0) { - // Explosion controlled region growing - typedef clitk::ExplosionControlledThresholdConnectedImageFilter ImageFilterType; - typename ImageFilterType::Pointer f= ImageFilterType::New(); - f->SetInput(working_input); - f->SetVerbose(false); - f->SetLower(-2000); - f->SetUpper(GetUpperThresholdForTrachea()); - f->SetMinimumLowerThreshold(-2000); - f->SetMaximumUpperThreshold(0); - f->SetAdaptLowerBorder(false); - f->SetAdaptUpperBorder(true); - f->SetMinimumSize(5000); - f->SetReplaceValue(1); - f->SetMultiplier(GetMultiplierForTrachea()); - f->SetThresholdStepSize(GetThresholdStepSizeForTrachea()); - f->SetMinimumThresholdStepSize(1); - for(unsigned int i=0; iAddSeed(m_Seeds[i]); - } - f->Update(); - trachea_tmp = f->GetOutput(); - // Set output - StopCurrentStep(trachea_tmp); - } - else { // Trachea not found - this->SetWarning("* WARNING * No seed found for trachea."); - StopCurrentStep(); - } + StartNewStep("Search for the trachea"); + SearchForTrachea(); + PrintMemory(GetVerboseMemoryFlag(), "After SearchForTrachea"); //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStepOrStop("Extract the lung with Otsu filter"); + StartNewStep("Extract the lung with Otsu filter"); // Automated Otsu thresholding and relabeling - typedef itk::OtsuThresholdImageFilter OtsuThresholdImageFilterType; + typedef itk::OtsuThresholdImageFilter OtsuThresholdImageFilterType; typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New(); otsuFilter->SetInput(working_input); otsuFilter->SetNumberOfHistogramBins(GetNumberOfHistogramBins()); otsuFilter->SetInsideValue(this->GetForegroundValue()); otsuFilter->SetOutsideValue(this->GetBackgroundValue()); otsuFilter->Update(); - working_image = otsuFilter->GetOutput(); + working_mask = otsuFilter->GetOutput(); // Set output - StopCurrentStep(working_image); + StopCurrentStep(working_mask); + PrintMemory(GetVerboseMemoryFlag(), "After Otsufilter"); //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStepOrStop("Select labels"); + StartNewStep("Select labels"); // Keep right labels - working_image = LabelizeAndSelectLabels - (working_image, + working_mask = LabelizeAndSelectLabels + (working_mask, GetBackgroundValue(), GetForegroundValue(), false, @@ -302,43 +290,46 @@ GenerateOutputInformation() GetLabelizeParameters2()); // Set output - StopCurrentStep(working_image); + StopCurrentStep(working_mask); + PrintMemory(GetVerboseMemoryFlag(), "After LabelizeAndSelectLabels"); //-------------------------------------------------------------------- //-------------------------------------------------------------------- if (m_Seeds.size() != 0) { // if ==0 ->no trachea found - StartNewStepOrStop("Remove the trachea"); + StartNewStep("Remove the trachea"); // Set the trachea - working_image = SetBackground - (working_image, trachea_tmp, 1, -1); + working_mask = SetBackground + (working_mask, trachea, 1, -1, true); + PrintMemory(GetVerboseMemoryFlag(), "After SetBackground"); - // Dilate the trachea + // Dilate the trachea static const unsigned int Dim = ImageType::ImageDimension; typedef itk::BinaryBallStructuringElement KernelType; KernelType structuringElement; structuringElement.SetRadius(GetRadiusForTrachea()); structuringElement.CreateStructuringElement(); - typedef clitk::ConditionalBinaryDilateImageFilter ConditionalBinaryDilateImageFilterType; + typedef clitk::ConditionalBinaryDilateImageFilter ConditionalBinaryDilateImageFilterType; typename ConditionalBinaryDilateImageFilterType::Pointer dilateFilter = ConditionalBinaryDilateImageFilterType::New(); dilateFilter->SetBoundaryToForeground(false); dilateFilter->SetKernel(structuringElement); dilateFilter->SetBackgroundValue (1); dilateFilter->SetForegroundValue (-1); - dilateFilter->SetInput (working_image); + dilateFilter->SetInput (working_mask); dilateFilter->Update(); - working_image = dilateFilter->GetOutput(); + working_mask = dilateFilter->GetOutput(); + PrintMemory(GetVerboseMemoryFlag(), "After dilate"); // Set trachea with dilatation - trachea_tmp = SetBackground - (trachea_tmp, working_image, -1, this->GetForegroundValue()); + trachea = SetBackground + (trachea, working_mask, -1, this->GetForegroundValue(), true); // Remove the trachea - working_image = SetBackground - (working_image, working_image, -1, this->GetBackgroundValue()); + working_mask = SetBackground + (working_mask, working_mask, -1, this->GetBackgroundValue(), true); // Label - working_image = LabelizeAndSelectLabels - (working_image, + working_mask = LabelizeAndSelectLabels + (working_mask, GetBackgroundValue(), GetForegroundValue(), false, @@ -346,229 +337,593 @@ GenerateOutputInformation() GetLabelizeParameters3()); // Set output - StopCurrentStep(working_image); + StopCurrentStep(working_mask); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- - typedef clitk::AutoCropFilter CropFilterType; - typename CropFilterType::Pointer cropFilter = CropFilterType::New(); + PrintMemory(GetVerboseMemoryFlag(), "before autocropfilter"); if (m_Seeds.size() != 0) { // if ==0 ->no trachea found - StartNewStepOrStop("Croping trachea"); - cropFilter->SetInput(trachea_tmp); - cropFilter->Update(); // Needed - typedef itk::CastImageFilter CastImageFilterType; - typename CastImageFilterType::Pointer caster= CastImageFilterType::New(); - caster->SetInput(cropFilter->GetOutput()); - caster->Update(); - trachea = caster->GetOutput(); + if (GetAutoCrop()) + trachea = clitk::AutoCrop(trachea, GetBackgroundValue()); + else + { + // Remove Padding region + typedef itk::CropImageFilter CropFilterType; + typename CropFilterType::Pointer cropFilter = CropFilterType::New(); + cropFilter->SetInput(trachea); + cropFilter->SetLowerBoundaryCropSize(bounds); + cropFilter->SetUpperBoundaryCropSize(bounds); + cropFilter->Update(); + trachea = cropFilter->GetOutput(); + } StopCurrentStep(trachea); + PrintMemory(GetVerboseMemoryFlag(), "after delete trachea"); } + PrintMemory(GetVerboseMemoryFlag(), "after delete trachea"); //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStepOrStop("Croping lung"); - typename CropFilterType::Pointer cropFilter2 = CropFilterType::New(); // Needed to reset pipeline - cropFilter2->SetInput(working_image); - cropFilter2->Update(); - working_image = cropFilter2->GetOutput(); - StopCurrentStep(working_image); + StartNewStep("Cropping lung"); + PrintMemory(GetVerboseMemoryFlag(), "Before Autocropfilter"); + if (GetAutoCrop()) + working_mask = clitk::AutoCrop(working_mask, GetBackgroundValue()); + else + { + // Remove Padding region + typedef itk::CropImageFilter CropFilterType; + typename CropFilterType::Pointer cropFilter = CropFilterType::New(); + cropFilter->SetInput(working_mask); + cropFilter->SetLowerBoundaryCropSize(bounds); + cropFilter->SetUpperBoundaryCropSize(bounds); + cropFilter->Update(); + working_mask = cropFilter->GetOutput(); + } + StopCurrentStep(working_mask); //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStepOrStop("Separate Left/Right lungs"); - // Initial label - working_image = Labelize(working_image, - GetBackgroundValue(), - false, - GetMinimalComponentSize()); - - // Count the labels - typedef itk::StatisticsImageFilter StatisticsImageFilterType; - typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New(); - statisticsImageFilter->SetInput(working_image); - statisticsImageFilter->Update(); - unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum(); - working_image = statisticsImageFilter->GetOutput(); - - // Decompose the first label - static const unsigned int Dim = ImageType::ImageDimension; - if (initialNumberOfLabels<2) { - // Structuring element radius - typename ImageType::SizeType radius; - for (unsigned int i=0;i DecomposeAndReconstructFilterType; - typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New(); - decomposeAndReconstructFilter->SetInput(working_image); - decomposeAndReconstructFilter->SetVerbose(false); - decomposeAndReconstructFilter->SetRadius(radius); - decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2); - decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize()); - decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1); - decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue()); - decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue()); - decomposeAndReconstructFilter->SetFullyConnected(true); - decomposeAndReconstructFilter->SetNumberOfNewLabels(1); - decomposeAndReconstructFilter->Update(); - working_image = decomposeAndReconstructFilter->GetOutput(); + // Final OpenClose + if (GetOpenCloseFlag()) { + StartNewStep("Open/Close"); + PrintMemory(GetVerboseMemoryFlag(), "Before OpenClose"); + + // Structuring element + typedef itk::BinaryBallStructuringElement KernelType; + KernelType structuringElement; + structuringElement.SetRadius(GetOpenCloseRadius()); + structuringElement.CreateStructuringElement(); + + // Open + typedef itk::BinaryMorphologicalOpeningImageFilter OpenFilterType; + typename OpenFilterType::Pointer openFilter = OpenFilterType::New(); + openFilter->SetInput(working_mask); + openFilter->SetBackgroundValue(GetBackgroundValue()); + openFilter->SetForegroundValue(GetForegroundValue()); + openFilter->SetKernel(structuringElement); + + // Close + typedef itk::BinaryMorphologicalClosingImageFilter CloseFilterType; + typename CloseFilterType::Pointer closeFilter = CloseFilterType::New(); + closeFilter->SetInput(openFilter->GetOutput()); + closeFilter->SetSafeBorder(true); + closeFilter->SetForegroundValue(GetForegroundValue()); + closeFilter->SetKernel(structuringElement); + closeFilter->Update(); + working_mask = closeFilter->GetOutput(); + } + + //-------------------------------------------------------------------- + //-------------------------------------------------------------------- + // Fill Lungs + if (GetFillHolesFlag()) { + StartNewStep("Fill Holes"); + PrintMemory(GetVerboseMemoryFlag(), "Before Fill Holes"); + typedef clitk::FillMaskFilter FillMaskFilterType; + typename FillMaskFilterType::Pointer fillMaskFilter = FillMaskFilterType::New(); + fillMaskFilter->SetInput(working_mask); + fillMaskFilter->AddDirection(2); + //fillMaskFilter->AddDirection(1); + fillMaskFilter->Update(); + working_mask = fillMaskFilter->GetOutput(); + StopCurrentStep(working_mask); + } + + if (GetSeparateLungsFlag()) { + //-------------------------------------------------------------------- + //-------------------------------------------------------------------- + StartNewStep("Separate Left/Right lungs"); + PrintMemory(GetVerboseMemoryFlag(), "Before Separate"); + // Initial label + working_mask = Labelize(working_mask, + GetBackgroundValue(), + false, + GetMinimalComponentSize()); + + PrintMemory(GetVerboseMemoryFlag(), "After Labelize"); + + // Count the labels + typedef itk::StatisticsImageFilter StatisticsImageFilterType; + typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New(); + statisticsImageFilter->SetInput(working_mask); + statisticsImageFilter->Update(); + unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum(); + working_mask = statisticsImageFilter->GetOutput(); + + PrintMemory(GetVerboseMemoryFlag(), "After count label"); + + // Decompose the first label + if (initialNumberOfLabels<2) { + // Structuring element radius + typename ImageType::SizeType radius; + for (unsigned int i=0;i DecomposeAndReconstructFilterType; + typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New(); + decomposeAndReconstructFilter->SetInput(working_mask); + decomposeAndReconstructFilter->SetVerbose(false); + decomposeAndReconstructFilter->SetRadius(radius); + decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2); + decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize()); + decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1); + decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue()); + decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue()); + decomposeAndReconstructFilter->SetFullyConnected(true); + decomposeAndReconstructFilter->SetNumberOfNewLabels(1); + decomposeAndReconstructFilter->Update(); + working_mask = decomposeAndReconstructFilter->GetOutput(); + } + PrintMemory(GetVerboseMemoryFlag(), "After decomposeAndReconstructFilter"); } - // Retain labels (lungs) - typedef itk::ThresholdImageFilter ThresholdImageFilterType; + // Retain labels ('1' is largset lung, so right. '2' is left) + typedef itk::ThresholdImageFilter ThresholdImageFilterType; typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New(); - thresholdFilter->SetInput(working_image); + thresholdFilter->SetInput(working_mask); thresholdFilter->ThresholdAbove(2); thresholdFilter->SetOutsideValue(this->GetBackgroundValue()); thresholdFilter->Update(); - working_image = thresholdFilter->GetOutput(); - StopCurrentStep (working_image); - - // Final Cast - StartNewStepOrStop("Final cast"); - typedef itk::CastImageFilter CastImageFilterType; - typename CastImageFilterType::Pointer caster= CastImageFilterType::New(); - caster->SetInput(working_image); - caster->Update(); - output = caster->GetOutput(); + working_mask = thresholdFilter->GetOutput(); + StopCurrentStep (working_mask); + PrintMemory(GetVerboseMemoryFlag(), "After Thresholdfilter"); // Update output info - this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion()); - - // Try to extract bifurcation in the trachea (bronchi) - // STILL EXPERIMENTAL - if (GetFindBronchialBifurcations()) { - StartNewStepOrStop("Find bronchial bifurcations"); - // Step 1 : extract skeleton - // Define the thinning filter - typedef itk::BinaryThinningImageFilter3D ThinningFilterType; - typename ThinningFilterType::Pointer thinningFilter = ThinningFilterType::New(); - thinningFilter->SetInput(trachea); - thinningFilter->Update(); - typename MaskImageType::Pointer skeleton = thinningFilter->GetOutput(); - writeImage(skeleton, "skeleton.mhd"); - - // Step 2 : tracking - DD("tracking"); - - // Step 2.1 : find first point for tracking - typedef itk::ImageRegionConstIteratorWithIndex IteratorType; - IteratorType it(skeleton, skeleton->GetLargestPossibleRegion()); - it.GoToReverseBegin(); - while ((!it.IsAtEnd()) && (it.Get() == GetBackgroundValue())) { - --it; + // output = working_mask; + //this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion()); + + // this->GetOutput(0)->SetRegions(working_mask->GetLargestPossibleRegion()); + +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLungFilter:: +GenerateInputRequestedRegion() { + // DD("GenerateInputRequestedRegion (nothing?)"); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLungFilter:: +GenerateData() +{ + // Set the output + // this->GraftOutput(output); // not SetNthOutput + this->GraftOutput(working_mask); // not SetNthOutput + // Store image filenames into AFDB + GetAFDB()->SetImageFilename("Lungs", this->GetOutputLungFilename()); + GetAFDB()->SetImageFilename("Trachea", this->GetOutputTracheaFilename()); + WriteAFDB(); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +bool +clitk::ExtractLungFilter:: +SearchForTracheaSeed(int skip) +{ + if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user) + // Restart the search until a seed is found, skipping more and more slices + bool stop = false; + while (!stop) { + // Search seed (parameters = UpperThresholdForTrachea) + static const unsigned int Dim = ImageType::ImageDimension; + typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion(); + typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize(); + typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex(); + sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-skip-5; + sliceRegionSize[Dim-1]=5; + sliceRegion.SetSize(sliceRegionSize); + sliceRegion.SetIndex(sliceRegionIndex); + typedef itk::ImageRegionConstIterator IteratorType; + IteratorType it(working_input, sliceRegion); + it.GoToBegin(); + while (!it.IsAtEnd()) { + if(it.Get() < GetUpperThresholdForTrachea() ) { + AddSeed(it.GetIndex()); + // DD(it.GetIndex()); + } + ++it; + } + + // if we do not found : restart + stop = (m_Seeds.size() != 0); + if (!stop) { + if (GetVerboseStepFlag()) { + std::cout << "\t No seed found this time. I skip some slices and restart." << std::endl; + } + if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) { + // we want to skip more than a half of the image, it is probably a bug + std::cerr << "Number of slices to skip to find trachea to high = " << skip << std::endl; + stop = true; + } + skip += 5; + } + else { + // DD(m_Seeds[0]); + // DD(m_Seeds.size()); + } } - if (it.IsAtEnd()) { - this->SetLastError("ERROR: first point in the skeleton not found ! Abord"); - return; + } + return (m_Seeds.size() != 0); +} +//-------------------------------------------------------------------- + + +bool is_orientation_superior(itk::SpatialOrientation::ValidCoordinateOrientationFlags orientation) +{ + itk::SpatialOrientation::CoordinateTerms sup = itk::SpatialOrientation::ITK_COORDINATE_Superior; + std::cout << "orientation: " << std::hex << orientation << "; superior: " << std::hex << sup << std::endl; + std::cout << std::dec; + + bool primary = (orientation & 0x0000ff) == sup; + bool secondary = ((orientation & 0x00ff00) >> 8) == sup; + bool tertiary = ((orientation & 0xff0000) >> 16) == sup; + return primary || secondary || tertiary; +} + +//-------------------------------------------------------------------- +template +bool +clitk::ExtractLungFilter:: +SearchForTracheaSeed2(int numberOfSlices) +{ + if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user) + if (GetVerboseFlag()) + std::cout << "SearchForTracheaSeed2(" << numberOfSlices << ", " << GetMaxElongation() << ")" << std::endl; + + typedef unsigned char MaskPixelType; + typedef itk::Image MaskImageType; + typedef itk::Image MaskImageType2D; + typedef itk::BinaryThresholdImageFilter ThresholdFilterType; + typedef itk::BinaryBallStructuringElement KernelType; + typedef itk::BinaryMorphologicalClosingImageFilter ClosingFilterType; + typedef itk::BinaryMorphologicalOpeningImageFilter OpeningFilterType; + typedef itk::ExtractImageFilter SlicerFilterType; + typedef itk::ConnectedComponentImageFilter LabelFilterType; + typedef itk::ShapeLabelObject ShapeLabelType; + typedef itk::LabelMap LabelImageType; + typedef itk::LabelImageToShapeLabelMapFilter ImageToLabelMapFilterType; + typedef itk::LabelMapToLabelImageFilter LabelMapToImageFilterType; + typedef itk::ImageFileWriter FileWriterType; + + // threshold to isolate airawys and lungs + typename ThresholdFilterType::Pointer threshold = ThresholdFilterType::New(); + threshold->SetLowerThreshold(-2000); + threshold->SetUpperThreshold(GetSeedPreProcessingThreshold()); + threshold->SetInput(working_input); + threshold->Update(); + + KernelType kernel_closing, kernel_opening; + + // remove small noise + typename OpeningFilterType::Pointer opening = OpeningFilterType::New(); + kernel_opening.SetRadius(1); + opening->SetKernel(kernel_opening); + opening->SetInput(threshold->GetOutput()); + opening->Update(); + + typename SlicerFilterType::Pointer slicer = SlicerFilterType::New(); + slicer->SetInput(opening->GetOutput()); + + // label result + typename LabelFilterType::Pointer label_filter = LabelFilterType::New(); + label_filter->SetInput(slicer->GetOutput()); + + // extract shape information from labels + typename ImageToLabelMapFilterType::Pointer label_to_map_filter = ImageToLabelMapFilterType::New(); + label_to_map_filter->SetInput(label_filter->GetOutput()); + + typename LabelMapToImageFilterType::Pointer map_to_label_filter = LabelMapToImageFilterType::New(); + typename FileWriterType::Pointer writer = FileWriterType::New(); + + typename ImageType::IndexType index; + typename ImageType::RegionType region = working_input->GetLargestPossibleRegion(); + typename ImageType::SizeType size = region.GetSize(); + typename ImageType::SpacingType spacing = working_input->GetSpacing(); + typename ImageType::PointType origin = working_input->GetOrigin(); + + int nslices = min(numberOfSlices, size[2]); + int start = 0, increment = 1; + itk::SpatialOrientationAdapter orientation; + typename ImageType::DirectionType dir = working_input->GetDirection(); + if (!is_orientation_superior(orientation.FromDirectionCosines(dir))) { + start = size[2]-1; + increment = -1; } - DD(skeleton->GetLargestPossibleRegion().GetIndex()); - typename MaskImageType::IndexType index = it.GetIndex(); - DD(index); - // Step 2.2 : initialize neighborhooditerator - typedef itk::NeighborhoodIterator NeighborhoodIteratorType; - typename NeighborhoodIteratorType::SizeType radius; - radius.Fill(1); - NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion()); - DD(nit.GetSize()); - DD(nit.Size()); + typename MaskImageType::PointType image_centre; + image_centre[0] = size[0]/2; + image_centre[1] = size[1]/2; + image_centre[2] = 0; + + typedef InternalIndexType SeedType; + SeedType trachea_centre, shape_centre, max_e_centre, prev_e_centre; + typedef std::list PointListType; + typedef std::list SequenceListType; + PointListType* current_sequence = NULL; + SequenceListType sequence_list; + + prev_e_centre.Fill(0); + std::ostringstream file_name; + index[0] = index[1] = 0; + size[0] = size[1] = 512; + size[2] = 0; + while (nslices--) { + index[2] = start; + start += increment; + + region.SetIndex(index); + region.SetSize(size); + slicer->SetExtractionRegion(region); + slicer->Update(); + label_filter->SetInput(slicer->GetOutput()); + label_filter->Update(); + + label_to_map_filter->SetInput(label_filter->GetOutput()); + label_to_map_filter->Update(); + typename LabelImageType::Pointer label_map = label_to_map_filter->GetOutput(); + + if (GetWriteStepFlag()) { + map_to_label_filter->SetInput(label_map); + writer->SetInput(map_to_label_filter->GetOutput()); + file_name.str(""); + file_name << "labels_"; + file_name.width(3); + file_name.fill('0'); + file_name << index[2] << ".mhd"; + writer->SetFileName(file_name.str().c_str()); + writer->Update(); + } + + typename LabelImageType::LabelObjectContainerType shapes_map = label_map->GetLabelObjectContainer(); + typename LabelImageType::LabelObjectContainerType::const_iterator s; + typename ShapeLabelType::Pointer shape, max_e_shape; + double max_elongation = GetMaxElongation(); + double max_size = size[0]*size[1]/128; + double max_e = 0; + int nshapes = 0; + for (s = shapes_map.begin(); s != shapes_map.end(); s++) { + shape = s->second; + if (shape->GetSize() > 150 && shape->GetSize() <= max_size) { + double e = 1 - 1/shape->GetBinaryElongation(); + //double area = 1 - r->Area() ; + if (e < max_elongation) { + nshapes++; + shape_centre[0] = (shape->GetCentroid()[0] - origin[0])/spacing[0]; + shape_centre[1] = (shape->GetCentroid()[1] - origin[1])/spacing[1]; + shape_centre[2] = index[2]; + //double d = 1 - (shape_centre - image_centre).Magnitude()/max_dist; + double dx = shape_centre[0] - image_centre[0]; + double d = 1 - dx*2/size[0]; + e = e + d; + if (e > max_e) + { + max_e = e; + max_e_shape = shape; + max_e_centre = shape_centre; + } + } + } + } + + if (nshapes > 0) + { + itk::Point p1, p2; + p1[0] = max_e_centre[0]; + p1[1] = max_e_centre[1]; + p1[2] = max_e_centre[2]; + + p2[0] = prev_e_centre[0]; + p2[1] = prev_e_centre[1]; + p2[2] = prev_e_centre[2]; + + double mag = (p2 - p1).GetNorm(); + if (GetVerboseFlag()) { + cout.precision(3); + cout << index[2] << ": "; + cout << "region(" << max_e_centre[0] << ", " << max_e_centre[1] << ", " << max_e_centre[2] << "); "; + cout << "prev_region(" << prev_e_centre[0] << ", " << prev_e_centre[1] << ", " << prev_e_centre[2] << "); "; + cout << "mag(" << mag << "); " << endl; + } + + if (mag > 5) + { + PointListType point_list; + point_list.push_back(max_e_centre); + sequence_list.push_back(point_list); + current_sequence = &sequence_list.back(); + } + else if (current_sequence) + current_sequence->push_back(max_e_centre); + + prev_e_centre= max_e_centre; + } + } - // Find first label number (must be different from BG and FG) - typename MaskImageType::PixelType label = GetForegroundValue()+1; - while ((label == GetBackgroundValue()) || (label == GetForegroundValue())) { label++; } - DD(label); - - // Track from the first point - std::vector listOfBifurcations; - TrackFromThisIndex(listOfBifurcations, skeleton, index, label); - DD("end track"); - DD(listOfBifurcations.size()); - writeImage(skeleton, "skeleton2.mhd"); - - for(unsigned int i=0; iTransformIndexToPhysicalPoint(listOfBifurcations[i].index, p); - DD(p); + size_t longest = 0; + for (typename SequenceListType::iterator s = sequence_list.begin(); s != sequence_list.end(); s++) + { + if (s->size() > longest) + { + longest = s->size(); + trachea_centre = s->front(); + } } - + + if (GetVerboseFlag()) + std::cout << "seed at: " << trachea_centre << std::endl; + m_Seeds.push_back(trachea_centre); } + + return (m_Seeds.size() != 0); } //-------------------------------------------------------------------- - //-------------------------------------------------------------------- -template +template void -clitk::ExtractLungFilter:: -GenerateData() { - // Do not put some "startnewstep" here, because the object if - // modified and the filter's pipeline it do two times. But it is - // required to quit if MustStop was set before. - if (GetMustStop()) return; +clitk::ExtractLungFilter:: +TracheaRegionGrowing() +{ + // Explosion controlled region growing + PrintMemory(GetVerboseMemoryFlag(), "Before ExplosionControlledThresholdConnectedImageFilter"); + typedef clitk::ExplosionControlledThresholdConnectedImageFilter ImageFilterType; + typename ImageFilterType::Pointer f= ImageFilterType::New(); + f->SetInput(working_input); + f->SetLower(-2000); + f->SetUpper(GetUpperThresholdForTrachea()); + f->SetMinimumLowerThreshold(-2000); + // f->SetMaximumUpperThreshold(0); // MAYBE TO CHANGE ??? + f->SetMaximumUpperThreshold(-700); // MAYBE TO CHANGE ??? + f->SetAdaptLowerBorder(false); + f->SetAdaptUpperBorder(true); + f->SetMinimumSize(5000); + f->SetReplaceValue(1); + f->SetMultiplier(GetMultiplierForTrachea()); + f->SetThresholdStepSize(GetThresholdStepSizeForTrachea()); + f->SetMinimumThresholdStepSize(1); + f->SetVerbose(GetVerboseRegionGrowingFlag()); + for(unsigned int i=0; iAddSeed(m_Seeds[i]); + // DD(m_Seeds[i]); + } + f->Update(); + PrintMemory(GetVerboseMemoryFlag(), "After RG update"); - // If everything goes well, set the output - this->GraftOutput(output); // not SetNthOutput + // take first (main) connected component + trachea = Labelize(f->GetOutput(), + GetBackgroundValue(), + true, + 1000);//GetMinimalComponentSize()); + PrintMemory(GetVerboseMemoryFlag(), "After Labelize"); + trachea = KeepLabels(trachea, + GetBackgroundValue(), + GetForegroundValue(), + 1, 1, false); + PrintMemory(GetVerboseMemoryFlag(), "After KeepLabels"); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -template -void -clitk::ExtractLungFilter:: -TrackFromThisIndex(std::vector & listOfBifurcations, - MaskImagePointer skeleton, - MaskImageIndexType index, - MaskImagePixelType label) { - DD("TrackFromThisIndex"); - DD(index); - DD((int)label); - // Create NeighborhoodIterator - typedef itk::NeighborhoodIterator NeighborhoodIteratorType; - typename NeighborhoodIteratorType::SizeType radius; - radius.Fill(1); - NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion()); - - // Track - std::vector listOfTrackedPoint; +template +double +clitk::ExtractLungFilter:: +ComputeTracheaVolume() +{ + typedef itk::ImageRegionConstIterator IteratorType; + IteratorType iter(trachea, trachea->GetLargestPossibleRegion()); + iter.GoToBegin(); + double volume = 0.0; + while (!iter.IsAtEnd()) { + if (iter.Get() == this->GetForegroundValue()) volume++; + ++iter; + } + + double voxelsize = trachea->GetSpacing()[0]*trachea->GetSpacing()[1]*trachea->GetSpacing()[2]; + return volume*voxelsize; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLungFilter:: +SearchForTrachea() +{ + // Search for seed among n slices, skip some slices before starting + // if not found -> skip more and restart + + // when seed found : perform region growing + // compute trachea volume + // if volume not plausible -> skip more slices and restart + + bool has_seed; bool stop = false; + double volume = 0.0; + int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed(); while (!stop) { - nit.SetLocation(index); - // DD((int)nit.GetCenterPixel()); - nit.SetCenterPixel(label); - listOfTrackedPoint.clear(); - for(unsigned int i=0; i(trachea, "step-trachea-"+toString(skip)+".mhd"); } - } - // DD(listOfTrackedPoint.size()); - if (listOfTrackedPoint.size() == 1) { - index = listOfTrackedPoint[0]; - } - else { - if (listOfTrackedPoint.size() == 2) { - BifurcationType bif(index, label, label+1, label+2); - listOfBifurcations.push_back(bif); - TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[0], label+1); - TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[1], label+2); + if (GetTracheaVolumeMustBeCheckedFlag()) { + if ((volume > 10) && (volume < 65 )) { // depend on image size ... + // Typical volume 22.59 cm 3 (± 7.69 cm 3 ) [ Leader 2004 ] + if (GetVerboseStepFlag()) + { + std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl; + } + } + else { + if (GetVerboseStepFlag()) { + std::cout << "\t The volume of the trachea (" << volume + << " cc) seems not correct. I skip some slices (" << skip << ") and restart to find seeds." + << std::endl; + } + skip += 5; + stop = false; + // empty the list of seed + m_Seeds.clear(); + } + if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) { + // we want to skip more than a half of the image, it is probably a bug + std::cerr << "2 : Number of slices to skip to find trachea too high = " << skip << std::endl; + stop = true; + } } else { - if (listOfTrackedPoint.size() > 2) { - std::cerr << "too much bifurcation points ... ?" << std::endl; - exit(0); - } - // Else this it the end of the tracking + stop = true; } - stop = true; } } + + if (volume != 0.0) { + // Set output + StopCurrentStep(trachea); + } + else { // Trachea not found + this->SetWarning("* WARNING * No seed found for trachea."); + StopCurrentStep(); + } } //-------------------------------------------------------------------- - #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX