X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=segmentation%2FclitkExtractLungFilter.txx;h=0954c2a3de46a11249c21237b9644e9a4a123404;hb=5038b96891adaefcde5b5520bb875b9378a0f81a;hp=52276f4478bb0cd55a6040efd07a763bd291d120;hpb=5e2af376544fce0c6dc46bb3c3227d35b501c1f1;p=clitk.git diff --git a/segmentation/clitkExtractLungFilter.txx b/segmentation/clitkExtractLungFilter.txx index 52276f4..0954c2a 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 @@ -26,6 +26,7 @@ #include "clitkAutoCropFilter.h" #include "clitkCropLikeImageFilter.h" #include "clitkFillMaskFilter.h" +#include "clitkMemoryUsage.h" // itk #include "itkBinaryThresholdImageFilter.h" @@ -36,6 +37,7 @@ #include "itkImageIteratorWithIndex.h" #include "itkBinaryMorphologicalOpeningImageFilter.h" #include "itkBinaryMorphologicalClosingImageFilter.h" +#include "itkConstantPadImageFilter.h" //-------------------------------------------------------------------- template @@ -54,6 +56,7 @@ ExtractLungFilter(): SetBackgroundValue(0); // Must be zero SetForegroundValue(1); SetMinimalComponentSize(100); + VerboseRegionGrowingFlagOff(); // Step 1 default values SetUpperThreshold(-300); @@ -70,6 +73,7 @@ ExtractLungFilter(): SetMultiplierForTrachea(5); SetThresholdStepSizeForTrachea(64); SetNumberOfSlicesToSkipBeforeSearchingSeed(0); + TracheaVolumeMustBeCheckedFlagOn(); // Step 3 default values SetNumberOfHistogramBins(500); @@ -88,11 +92,12 @@ ExtractLungFilter(): SetLabelizeParameters3(p3); // Step 5 - OpenCloseOff(); + OpenCloseFlagOff(); SetOpenCloseRadius(1); // Step 6 - FillHolesOn(); + FillHolesFlagOn(); + AutoCropOn(); } //-------------------------------------------------------------------- @@ -119,58 +124,13 @@ AddSeed(InternalIndexType s) //-------------------------------------------------------------------- -//-------------------------------------------------------------------- -template -template -void -clitk::ExtractLungFilter:: -SetArgsInfo(ArgsInfoType mArgsInfo) -{ - SetVerboseOption_GGO(mArgsInfo); - SetVerboseStep_GGO(mArgsInfo); - SetWriteStep_GGO(mArgsInfo); - SetVerboseWarningOff_GGO(mArgsInfo); - - SetAFDBFilename_GGO(mArgsInfo); - SetOutputLungFilename_GGO(mArgsInfo); - SetOutputTracheaFilename_GGO(mArgsInfo); - - SetUpperThreshold_GGO(mArgsInfo); - SetLowerThreshold_GGO(mArgsInfo); - SetNumberOfSlicesToSkipBeforeSearchingSeed_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); - - SetOpenCloseRadius_GGO(mArgsInfo); - SetOpenClose_GGO(mArgsInfo); - - SetFillHoles_GGO(mArgsInfo); -} -//-------------------------------------------------------------------- - - //-------------------------------------------------------------------- template void clitk::ExtractLungFilter:: GenerateOutputInformation() { + clitk::PrintMemory(GetVerboseMemoryFlag(), "Initial memory"); // OK Superclass::GenerateOutputInformation(); // Read DB @@ -178,27 +138,57 @@ GenerateOutputInformation() // Get input pointers input = dynamic_cast(itk::ProcessObject::GetInput(0)); - patient = GetAFDB()->template GetImage ("patient"); + patient = GetAFDB()->template GetImage ("Patient"); + PrintMemory(GetVerboseMemoryFlag(), "After reading patient"); // OK //-------------------------------------------------------------------- //-------------------------------------------------------------------- // Crop input like patient image (must have the same spacing) - StartNewStep("Crop input image to 'patient' extends"); + // 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(); - DD(working_input->GetLargestPossibleRegion()); + // 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 - (working_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 + */ //-------------------------------------------------------------------- //-------------------------------------------------------------------- @@ -210,33 +200,41 @@ GenerateOutputInformation() } } // 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); @@ -245,29 +243,31 @@ GenerateOutputInformation() //-------------------------------------------------------------------- StartNewStep("Search for the trachea"); SearchForTrachea(); + PrintMemory(GetVerboseMemoryFlag(), "After SearchForTrachea"); //-------------------------------------------------------------------- //-------------------------------------------------------------------- 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"); //-------------------------------------------------------------------- //-------------------------------------------------------------------- StartNewStep("Select labels"); // Keep right labels - working_image = LabelizeAndSelectLabels - (working_image, + working_mask = LabelizeAndSelectLabels + (working_mask, GetBackgroundValue(), GetForegroundValue(), false, @@ -275,15 +275,17 @@ GenerateOutputInformation() GetLabelizeParameters2()); // Set output - StopCurrentStep(working_image); + StopCurrentStep(working_mask); + PrintMemory(GetVerboseMemoryFlag(), "After LabelizeAndSelectLabels"); //-------------------------------------------------------------------- //-------------------------------------------------------------------- if (m_Seeds.size() != 0) { // if ==0 ->no trachea found 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 static const unsigned int Dim = ImageType::ImageDimension; @@ -291,27 +293,28 @@ GenerateOutputInformation() 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, @@ -319,41 +322,57 @@ GenerateOutputInformation() GetLabelizeParameters3()); // Set output - StopCurrentStep(working_image); + StopCurrentStep(working_mask); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- - typedef clitk::AutoCropFilter AutoCropFilterType; - typename AutoCropFilterType::Pointer autocropFilter = AutoCropFilterType::New(); + PrintMemory(GetVerboseMemoryFlag(), "before autocropfilter"); if (m_Seeds.size() != 0) { // if ==0 ->no trachea found - StartNewStep("Cropping trachea"); - autocropFilter->SetInput(trachea_tmp); - autocropFilter->Update(); // Needed - typedef itk::CastImageFilter CastImageFilterType; - typename CastImageFilterType::Pointer caster= CastImageFilterType::New(); - caster->SetInput(autocropFilter->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"); //-------------------------------------------------------------------- //-------------------------------------------------------------------- StartNewStep("Cropping lung"); - typename AutoCropFilterType::Pointer autocropFilter2 = AutoCropFilterType::New(); // Needed to reset pipeline - autocropFilter2->SetInput(working_image); - autocropFilter2->Update(); - working_image = autocropFilter2->GetOutput(); - DD(working_image->GetLargestPossibleRegion()); - StopCurrentStep(working_image); + 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); //-------------------------------------------------------------------- //-------------------------------------------------------------------- // Final OpenClose - if (GetOpenClose()) { + if (GetOpenCloseFlag()) { StartNewStep("Open/Close"); - + PrintMemory(GetVerboseMemoryFlag(), "Before OpenClose"); + // Structuring element typedef itk::BinaryBallStructuringElement KernelType; KernelType structuringElement; @@ -361,67 +380,71 @@ GenerateOutputInformation() structuringElement.CreateStructuringElement(); // Open - typedef itk::BinaryMorphologicalOpeningImageFilter OpenFilterType; + typedef itk::BinaryMorphologicalOpeningImageFilter OpenFilterType; typename OpenFilterType::Pointer openFilter = OpenFilterType::New(); - openFilter->SetInput(working_image); + openFilter->SetInput(working_mask); openFilter->SetBackgroundValue(GetBackgroundValue()); openFilter->SetForegroundValue(GetForegroundValue()); openFilter->SetKernel(structuringElement); // Close - typedef itk::BinaryMorphologicalClosingImageFilter CloseFilterType; + 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_image = closeFilter->GetOutput(); + working_mask = closeFilter->GetOutput(); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- // Fill Lungs - if (GetFillHoles()) { + if (GetFillHolesFlag()) { StartNewStep("Fill Holes"); - typedef clitk::FillMaskFilter FillMaskFilterType; + PrintMemory(GetVerboseMemoryFlag(), "Before Fill Holes"); + typedef clitk::FillMaskFilter FillMaskFilterType; typename FillMaskFilterType::Pointer fillMaskFilter = FillMaskFilterType::New(); - fillMaskFilter->SetInput(working_image); + fillMaskFilter->SetInput(working_mask); fillMaskFilter->AddDirection(2); //fillMaskFilter->AddDirection(1); fillMaskFilter->Update(); - working_image = fillMaskFilter->GetOutput(); - StopCurrentStep(working_image); + working_mask = fillMaskFilter->GetOutput(); + StopCurrentStep(working_mask); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- StartNewStep("Separate Left/Right lungs"); + PrintMemory(GetVerboseMemoryFlag(), "Before Separate"); // Initial label - working_image = Labelize(working_image, - GetBackgroundValue(), - false, - GetMinimalComponentSize()); + working_mask = Labelize(working_mask, + GetBackgroundValue(), + false, + GetMinimalComponentSize()); + + PrintMemory(GetVerboseMemoryFlag(), "After Labelize"); // Count the labels - typedef itk::StatisticsImageFilter StatisticsImageFilterType; + typedef itk::StatisticsImageFilter StatisticsImageFilterType; typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New(); - statisticsImageFilter->SetInput(working_image); + statisticsImageFilter->SetInput(working_mask); statisticsImageFilter->Update(); unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum(); - working_image = statisticsImageFilter->GetOutput(); + working_mask = statisticsImageFilter->GetOutput(); + + PrintMemory(GetVerboseMemoryFlag(), "After count label"); // Decompose the first label - static const unsigned int Dim = ImageType::ImageDimension; if (initialNumberOfLabels<2) { - DD(initialNumberOfLabels); // Structuring element radius typename ImageType::SizeType radius; for (unsigned int i=0;i DecomposeAndReconstructFilterType; + typedef clitk::DecomposeAndReconstructImageFilter DecomposeAndReconstructFilterType; typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New(); - decomposeAndReconstructFilter->SetInput(working_image); - decomposeAndReconstructFilter->SetVerbose(true); + decomposeAndReconstructFilter->SetInput(working_mask); + decomposeAndReconstructFilter->SetVerbose(false); decomposeAndReconstructFilter->SetRadius(radius); decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2); decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize()); @@ -431,29 +454,37 @@ GenerateOutputInformation() decomposeAndReconstructFilter->SetFullyConnected(true); decomposeAndReconstructFilter->SetNumberOfNewLabels(1); decomposeAndReconstructFilter->Update(); - working_image = decomposeAndReconstructFilter->GetOutput(); + working_mask = decomposeAndReconstructFilter->GetOutput(); } + PrintMemory(GetVerboseMemoryFlag(), "After decomposeAndReconstructFilter"); // Retain labels ('1' is largset lung, so right. '2' is left) - typedef itk::ThresholdImageFilter ThresholdImageFilterType; + 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 - StartNewStep("Cast the lung mask"); - 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()); + // output = working_mask; + //this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion()); + + // this->GetOutput(0)->SetRegions(working_mask->GetLargestPossibleRegion()); + +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLungFilter:: +GenerateInputRequestedRegion() { + // DD("GenerateInputRequestedRegion (nothing?)"); } //-------------------------------------------------------------------- @@ -465,10 +496,11 @@ clitk::ExtractLungFilter:: GenerateData() { // Set the output - this->GraftOutput(output); // not SetNthOutput + // 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()); + GetAFDB()->SetImageFilename("Lungs", this->GetOutputLungFilename()); + GetAFDB()->SetImageFilename("Trachea", this->GetOutputTracheaFilename()); WriteAFDB(); } //-------------------------------------------------------------------- @@ -507,7 +539,7 @@ SearchForTracheaSeed(int skip) // if we do not found : restart stop = (m_Seeds.size() != 0); if (!stop) { - if (GetVerboseStep()) { + 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]) { @@ -517,6 +549,10 @@ SearchForTracheaSeed(int skip) } skip += 5; } + else { + // DD(m_Seeds[0]); + // DD(m_Seeds.size()); + } } } return (m_Seeds.size() != 0); @@ -531,10 +567,10 @@ clitk::ExtractLungFilter:: TracheaRegionGrowing() { // Explosion controlled region growing - typedef clitk::ExplosionControlledThresholdConnectedImageFilter ImageFilterType; + PrintMemory(GetVerboseMemoryFlag(), "Before ExplosionControlledThresholdConnectedImageFilter"); + 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); @@ -547,24 +583,25 @@ TracheaRegionGrowing() f->SetMultiplier(GetMultiplierForTrachea()); f->SetThresholdStepSize(GetThresholdStepSizeForTrachea()); f->SetMinimumThresholdStepSize(1); - f->VerboseOn(); + f->SetVerbose(GetVerboseRegionGrowingFlag()); for(unsigned int i=0; iAddSeed(m_Seeds[i]); // DD(m_Seeds[i]); } f->Update(); - - writeImage(f->GetOutput(), "trg.mhd"); - + PrintMemory(GetVerboseMemoryFlag(), "After RG update"); + // take first (main) connected component - trachea_tmp = Labelize(f->GetOutput(), - GetBackgroundValue(), - true, - GetMinimalComponentSize()); - trachea_tmp = KeepLabels(trachea_tmp, + trachea = Labelize(f->GetOutput(), + GetBackgroundValue(), + true, + 1000);//GetMinimalComponentSize()); + PrintMemory(GetVerboseMemoryFlag(), "After Labelize"); + trachea = KeepLabels(trachea, GetBackgroundValue(), GetForegroundValue(), 1, 1, false); + PrintMemory(GetVerboseMemoryFlag(), "After KeepLabels"); } //-------------------------------------------------------------------- @@ -576,7 +613,7 @@ clitk::ExtractLungFilter:: ComputeTracheaVolume() { typedef itk::ImageRegionConstIterator IteratorType; - IteratorType iter(trachea_tmp, trachea_tmp->GetLargestPossibleRegion()); + IteratorType iter(trachea, trachea->GetLargestPossibleRegion()); iter.GoToBegin(); double volume = 0.0; while (!iter.IsAtEnd()) { @@ -584,7 +621,7 @@ ComputeTracheaVolume() ++iter; } - double voxelsize = trachea_tmp->GetSpacing()[0]*trachea_tmp->GetSpacing()[1]*trachea_tmp->GetSpacing()[2]; + double voxelsize = trachea->GetSpacing()[0]*trachea->GetSpacing()[1]*trachea->GetSpacing()[2]; return volume*voxelsize; } //-------------------------------------------------------------------- @@ -592,7 +629,7 @@ ComputeTracheaVolume() //-------------------------------------------------------------------- template -void +void clitk::ExtractLungFilter:: SearchForTrachea() { @@ -607,34 +644,46 @@ SearchForTrachea() double volume = 0.0; int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed(); while (!stop) { - stop = SearchForTracheaSeed(skip); - if (stop) { + stop = true; + if (SearchForTracheaSeed(skip)) { TracheaRegionGrowing(); volume = ComputeTracheaVolume()/1000; // assume mm3, so divide by 1000 to get cc - if ((volume > 10) && (volume < 55 )) { // it is ok - // Typical volume 22.59 cm 3 (± 7.69 cm 3 ) [ Leader 2004 ] - if (GetVerboseStep()) { - std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl; - } - stop = true; + if (GetWriteStepFlag()) { + writeImage(trachea, "step-trachea-"+toString(skip)+".mhd"); + } + 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 (GetVerboseStep()) { - 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(); + stop = true; } } } if (volume != 0.0) { // Set output - StopCurrentStep(trachea_tmp); + StopCurrentStep(trachea); } else { // Trachea not found this->SetWarning("* WARNING * No seed found for trachea.");