X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=segmentation%2FclitkExtractLungFilter.txx;h=52276f4478bb0cd55a6040efd07a763bd291d120;hb=a24b0a699298efe54b53c53cb215455fecd633fe;hp=3bb57b8ffb0cf290dc18c3fe2892fc84ebe52696;hpb=6e16222234a90c6079a8f4696c92de7349a496bb;p=clitk.git diff --git a/segmentation/clitkExtractLungFilter.txx b/segmentation/clitkExtractLungFilter.txx index 3bb57b8..52276f4 100644 --- a/segmentation/clitkExtractLungFilter.txx +++ b/segmentation/clitkExtractLungFilter.txx @@ -24,6 +24,8 @@ #include "clitkSetBackgroundImageFilter.h" #include "clitkSegmentationUtils.h" #include "clitkAutoCropFilter.h" +#include "clitkCropLikeImageFilter.h" +#include "clitkFillMaskFilter.h" // itk #include "itkBinaryThresholdImageFilter.h" @@ -32,20 +34,22 @@ #include "itkOtsuThresholdImageFilter.h" #include "itkBinaryThinningImageFilter3D.h" #include "itkImageIteratorWithIndex.h" - +#include "itkBinaryMorphologicalOpeningImageFilter.h" +#include "itkBinaryMorphologicalClosingImageFilter.h" //-------------------------------------------------------------------- -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); @@ -83,16 +87,20 @@ ExtractLungFilter(): p3->UseLastKeepOff(); SetLabelizeParameters3(p3); - // Step 5 : find bronchial bifurcations - FindBronchialBifurcationsOn(); + // Step 5 + OpenCloseOff(); + SetOpenCloseRadius(1); + + // Step 6 + FillHolesOn(); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -template +template void -clitk::ExtractLungFilter:: +clitk::ExtractLungFilter:: SetInput(const ImageType * image) { this->SetNthInput(0, const_cast(image)); @@ -101,21 +109,9 @@ SetInput(const ImageType * image) //-------------------------------------------------------------------- -template -void -clitk::ExtractLungFilter:: -SetInputPatientMask(MaskImageType * image, MaskImagePixelType bg ) -{ - this->SetNthInput(1, const_cast(image)); - SetPatientMaskBackgroundValue(bg); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template +template void -clitk::ExtractLungFilter:: +clitk::ExtractLungFilter:: AddSeed(InternalIndexType s) { m_Seeds.push_back(s); @@ -124,10 +120,10 @@ AddSeed(InternalIndexType s) //-------------------------------------------------------------------- -template +template template void -clitk::ExtractLungFilter:: +clitk::ExtractLungFilter:: SetArgsInfo(ArgsInfoType mArgsInfo) { SetVerboseOption_GGO(mArgsInfo); @@ -135,6 +131,10 @@ SetArgsInfo(ArgsInfoType 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); @@ -156,44 +156,57 @@ SetArgsInfo(ArgsInfoType mArgsInfo) SetRadiusForTrachea_GGO(mArgsInfo); SetLabelizeParameters3_GGO(mArgsInfo); + + SetOpenCloseRadius_GGO(mArgsInfo); + SetOpenClose_GGO(mArgsInfo); + + SetFillHoles_GGO(mArgsInfo); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -template +template void -clitk::ExtractLungFilter:: +clitk::ExtractLungFilter:: GenerateOutputInformation() { Superclass::GenerateOutputInformation(); - //this->GetOutput(0)->SetRequestedRegion(this->GetOutput(0)->GetLargestPossibleRegion()); + + // 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"); - // 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) + StartNewStep("Crop input image to 'patient' extends"); + typedef clitk::CropLikeImageFilter CropImageFilter; + typename CropImageFilter::Pointer cropFilter = CropImageFilter::New(); + cropFilter->SetInput(input); + cropFilter->SetCropLikeImage(patient); + cropFilter->Update(); + working_input = cropFilter->GetOutput(); + DD(working_input->GetLargestPossibleRegion()); + StopCurrentStep(working_input); + + //-------------------------------------------------------------------- + //-------------------------------------------------------------------- + StartNewStep("Set background to initial image"); working_input = SetBackground - (input, patient, GetPatientMaskBackgroundValue(), -1000); + (working_input, patient, GetPatientMaskBackgroundValue(), -1000); StopCurrentStep(working_input); //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStepOrStop("Remove Air"); + StartNewStep("Remove Air"); // Check threshold if (m_UseLowerThreshold) { if (m_LowerThreshold > m_UpperThreshold) { - this->SetLastError("ERROR: lower threshold cannot be greater than upper threshold."); - return; + clitkExceptionMacro("lower threshold cannot be greater than upper threshold."); } } // Threshold to get air @@ -230,12 +243,12 @@ GenerateOutputInformation() //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStepOrStop("Search for the trachea"); + StartNewStep("Search for the trachea"); SearchForTrachea(); //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStepOrStop("Extract the lung with Otsu filter"); + StartNewStep("Extract the lung with Otsu filter"); // Automated Otsu thresholding and relabeling typedef itk::OtsuThresholdImageFilter OtsuThresholdImageFilterType; typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New(); @@ -251,7 +264,7 @@ GenerateOutputInformation() //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStepOrStop("Select labels"); + StartNewStep("Select labels"); // Keep right labels working_image = LabelizeAndSelectLabels (working_image, @@ -267,7 +280,7 @@ GenerateOutputInformation() //-------------------------------------------------------------------- //-------------------------------------------------------------------- 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); @@ -311,15 +324,15 @@ GenerateOutputInformation() //-------------------------------------------------------------------- //-------------------------------------------------------------------- - typedef clitk::AutoCropFilter CropFilterType; - typename CropFilterType::Pointer cropFilter = CropFilterType::New(); + typedef clitk::AutoCropFilter AutoCropFilterType; + typename AutoCropFilterType::Pointer autocropFilter = AutoCropFilterType::New(); if (m_Seeds.size() != 0) { // if ==0 ->no trachea found - StartNewStepOrStop("Croping trachea"); - cropFilter->SetInput(trachea_tmp); - cropFilter->Update(); // Needed + StartNewStep("Cropping trachea"); + autocropFilter->SetInput(trachea_tmp); + autocropFilter->Update(); // Needed typedef itk::CastImageFilter CastImageFilterType; typename CastImageFilterType::Pointer caster= CastImageFilterType::New(); - caster->SetInput(cropFilter->GetOutput()); + caster->SetInput(autocropFilter->GetOutput()); caster->Update(); trachea = caster->GetOutput(); StopCurrentStep(trachea); @@ -327,16 +340,63 @@ GenerateOutputInformation() //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStepOrStop("Croping lung"); - typename CropFilterType::Pointer cropFilter2 = CropFilterType::New(); // Needed to reset pipeline - cropFilter2->SetInput(working_image); - cropFilter2->Update(); - working_image = cropFilter2->GetOutput(); + 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); //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStepOrStop("Separate Left/Right lungs"); + // Final OpenClose + if (GetOpenClose()) { + StartNewStep("Open/Close"); + + // 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_image); + 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_image = closeFilter->GetOutput(); + } + + //-------------------------------------------------------------------- + //-------------------------------------------------------------------- + // Fill Lungs + if (GetFillHoles()) { + StartNewStep("Fill Holes"); + typedef clitk::FillMaskFilter FillMaskFilterType; + typename FillMaskFilterType::Pointer fillMaskFilter = FillMaskFilterType::New(); + fillMaskFilter->SetInput(working_image); + fillMaskFilter->AddDirection(2); + //fillMaskFilter->AddDirection(1); + fillMaskFilter->Update(); + working_image = fillMaskFilter->GetOutput(); + StopCurrentStep(working_image); + } + + //-------------------------------------------------------------------- + //-------------------------------------------------------------------- + StartNewStep("Separate Left/Right lungs"); // Initial label working_image = Labelize(working_image, GetBackgroundValue(), @@ -354,13 +414,14 @@ GenerateOutputInformation() // 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; typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New(); decomposeAndReconstructFilter->SetInput(working_image); - decomposeAndReconstructFilter->SetVerbose(false); + decomposeAndReconstructFilter->SetVerbose(true); decomposeAndReconstructFilter->SetRadius(radius); decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2); decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize()); @@ -384,7 +445,7 @@ GenerateOutputInformation() StopCurrentStep (working_image); // Final Cast - StartNewStepOrStop("Final cast"); + StartNewStep("Cast the lung mask"); typedef itk::CastImageFilter CastImageFilterType; typename CastImageFilterType::Pointer caster= CastImageFilterType::New(); caster->SetInput(working_image); @@ -393,153 +454,30 @@ GenerateOutputInformation() // Update output info this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion()); - - // Try to extract bifurcation in the trachea (bronchi) - // STILL EXPERIMENTAL - if (m_Seeds.size() != 0) { // if ==0 ->no trachea found - - 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; - } - if (it.IsAtEnd()) { - this->SetLastError("ERROR: first point in the skeleton not found ! Abort"); - return; - } - 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()); - - // 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); - } - - } - } } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -template +template void -clitk::ExtractLungFilter:: +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; - - // If everything goes well, set the output + // Set the output this->GraftOutput(output); // not SetNthOutput + // Store image filenames into AFDB + GetAFDB()->SetImageFilename("lungs", this->GetOutputLungFilename()); + GetAFDB()->SetImageFilename("trachea", this->GetOutputTracheaFilename()); + WriteAFDB(); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -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; - bool stop = false; - while (!stop) { - nit.SetLocation(index); - // DD((int)nit.GetCenterPixel()); - nit.SetCenterPixel(label); - listOfTrackedPoint.clear(); - for(unsigned int i=0; i 2) { - std::cerr << "too much bifurcation points ... ?" << std::endl; - exit(0); - } - // Else this it the end of the tracking - } - stop = true; - } - } -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template +template bool -clitk::ExtractLungFilter:: +clitk::ExtractLungFilter:: SearchForTracheaSeed(int skip) { if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user) @@ -587,9 +525,9 @@ SearchForTracheaSeed(int skip) //-------------------------------------------------------------------- -template +template void -clitk::ExtractLungFilter:: +clitk::ExtractLungFilter:: TracheaRegionGrowing() { // Explosion controlled region growing @@ -600,7 +538,8 @@ TracheaRegionGrowing() f->SetLower(-2000); f->SetUpper(GetUpperThresholdForTrachea()); f->SetMinimumLowerThreshold(-2000); - f->SetMaximumUpperThreshold(0); + // f->SetMaximumUpperThreshold(0); // MAYBE TO CHANGE ??? + f->SetMaximumUpperThreshold(-800); // MAYBE TO CHANGE ??? f->SetAdaptLowerBorder(false); f->SetAdaptUpperBorder(true); f->SetMinimumSize(5000); @@ -608,32 +547,32 @@ TracheaRegionGrowing() f->SetMultiplier(GetMultiplierForTrachea()); f->SetThresholdStepSize(GetThresholdStepSizeForTrachea()); f->SetMinimumThresholdStepSize(1); + f->VerboseOn(); for(unsigned int i=0; iAddSeed(m_Seeds[i]); // DD(m_Seeds[i]); } f->Update(); + writeImage(f->GetOutput(), "trg.mhd"); + // take first (main) connected component - writeImage(f->GetOutput(), "t1.mhd"); trachea_tmp = Labelize(f->GetOutput(), GetBackgroundValue(), true, GetMinimalComponentSize()); - writeImage(trachea_tmp, "t2.mhd"); trachea_tmp = KeepLabels(trachea_tmp, GetBackgroundValue(), GetForegroundValue(), 1, 1, false); - writeImage(trachea_tmp, "t3.mhd"); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -template +template double -clitk::ExtractLungFilter:: +clitk::ExtractLungFilter:: ComputeTracheaVolume() { typedef itk::ImageRegionConstIterator IteratorType; @@ -652,9 +591,9 @@ ComputeTracheaVolume() //-------------------------------------------------------------------- -template +template void -clitk::ExtractLungFilter:: +clitk::ExtractLungFilter:: SearchForTrachea() { // Search for seed among n slices, skip some slices before starting @@ -671,8 +610,8 @@ SearchForTrachea() stop = SearchForTracheaSeed(skip); if (stop) { TracheaRegionGrowing(); - volume = ComputeTracheaVolume(); // assume mm3 - if ((volume > 10000) && (volume < 55000 )) { // it is ok + 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;