X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=segmentation%2FclitkExtractLungFilter.txx;h=3de027cc3b8210d0dbbac0e5d2a24d1a4aa02cfa;hb=cc7d70464c9ece01281a987973607eebba30d8dc;hp=4baee45296f5555d49b1748195b19d3757106e13;hpb=c5669e2c629b8248e22d8892c6576ada7a14cc25;p=clitk.git diff --git a/segmentation/clitkExtractLungFilter.txx b/segmentation/clitkExtractLungFilter.txx index 4baee45..3de027c 100644 --- a/segmentation/clitkExtractLungFilter.txx +++ b/segmentation/clitkExtractLungFilter.txx @@ -32,7 +32,8 @@ #include "itkOtsuThresholdImageFilter.h" #include "itkBinaryThinningImageFilter3D.h" #include "itkImageIteratorWithIndex.h" - +#include "itkBinaryMorphologicalOpeningImageFilter.h" +#include "itkBinaryMorphologicalClosingImageFilter.h" //-------------------------------------------------------------------- template @@ -65,6 +66,7 @@ ExtractLungFilter(): SetUpperThresholdForTrachea(-900); SetMultiplierForTrachea(5); SetThresholdStepSizeForTrachea(64); + SetNumberOfSlicesToSkipBeforeSearchingSeed(0); // Step 3 default values SetNumberOfHistogramBins(500); @@ -82,8 +84,9 @@ ExtractLungFilter(): p3->UseLastKeepOff(); SetLabelizeParameters3(p3); - // Step 5 : find bronchial bifurcations - FindBronchialBifurcationsOn(); + // Step 5 + FinalOpenCloseOff(); + SetFinalOpenCloseRadius(1); } //-------------------------------------------------------------------- @@ -136,6 +139,7 @@ SetArgsInfo(ArgsInfoType mArgsInfo) SetUpperThreshold_GGO(mArgsInfo); SetLowerThreshold_GGO(mArgsInfo); + SetNumberOfSlicesToSkipBeforeSearchingSeed_GGO(mArgsInfo); SetLabelizeParameters1_GGO(mArgsInfo); if (!mArgsInfo.remove1_given) { GetLabelizeParameters1()->AddLabelToRemove(2); @@ -154,6 +158,9 @@ SetArgsInfo(ArgsInfoType mArgsInfo) SetRadiusForTrachea_GGO(mArgsInfo); SetLabelizeParameters3_GGO(mArgsInfo); + + SetFinalOpenCloseRadius_GGO(mArgsInfo); + SetFinalOpenClose_GGO(mArgsInfo); } //-------------------------------------------------------------------- @@ -165,7 +172,6 @@ clitk::ExtractLungFilter:: GenerateOutputInformation() { Superclass::GenerateOutputInformation(); - //this->GetOutput(0)->SetRequestedRegion(this->GetOutput(0)->GetLargestPossibleRegion()); // Get input pointers patient = dynamic_cast(itk::ProcessObject::GetInput(1)); @@ -173,25 +179,23 @@ GenerateOutputInformation() // Check image if (!HaveSameSizeAndSpacing(input, patient)) { - this->SetLastError("* ERROR * the images (input and patient mask) must have the same size & spacing"); - return; + clitkExceptionMacro("the 'input' and 'patient' masks must have the same size & spacing."); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStepOrStop("Set background to initial image"); + StartNewStep("Set background to initial image"); working_input = SetBackground (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 @@ -228,61 +232,12 @@ GenerateOutputInformation() //-------------------------------------------------------------------- //-------------------------------------------------------------------- - 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(); //-------------------------------------------------------------------- //-------------------------------------------------------------------- - 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(); @@ -298,7 +253,7 @@ GenerateOutputInformation() //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStepOrStop("Select labels"); + StartNewStep("Select labels"); // Keep right labels working_image = LabelizeAndSelectLabels (working_image, @@ -314,12 +269,12 @@ 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); - // Dilate the trachea + // Dilate the trachea static const unsigned int Dim = ImageType::ImageDimension; typedef itk::BinaryBallStructuringElement KernelType; KernelType structuringElement; @@ -361,7 +316,7 @@ GenerateOutputInformation() typedef clitk::AutoCropFilter CropFilterType; typename CropFilterType::Pointer cropFilter = CropFilterType::New(); if (m_Seeds.size() != 0) { // if ==0 ->no trachea found - StartNewStepOrStop("Croping trachea"); + StartNewStep("Croping trachea"); cropFilter->SetInput(trachea_tmp); cropFilter->Update(); // Needed typedef itk::CastImageFilter CastImageFilterType; @@ -374,7 +329,7 @@ GenerateOutputInformation() //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStepOrStop("Croping lung"); + StartNewStep("Croping lung"); typename CropFilterType::Pointer cropFilter2 = CropFilterType::New(); // Needed to reset pipeline cropFilter2->SetInput(working_image); cropFilter2->Update(); @@ -383,7 +338,38 @@ GenerateOutputInformation() //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStepOrStop("Separate Left/Right lungs"); + // Final OpenClose + if (GetFinalOpenClose()) { + StartNewStep("Open/Close"); + + // Structuring element + typedef itk::BinaryBallStructuringElement KernelType; + KernelType structuringElement; + structuringElement.SetRadius(GetFinalOpenCloseRadius()); + 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(); + } + + //-------------------------------------------------------------------- + //-------------------------------------------------------------------- + StartNewStep("Separate Left/Right lungs"); // Initial label working_image = Labelize(working_image, GetBackgroundValue(), @@ -420,7 +406,7 @@ GenerateOutputInformation() working_image = decomposeAndReconstructFilter->GetOutput(); } - // Retain labels (lungs) + // 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); @@ -431,7 +417,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); @@ -441,64 +427,7 @@ GenerateOutputInformation() // 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; - } - 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); - } - } } //-------------------------------------------------------------------- @@ -507,13 +436,9 @@ GenerateOutputInformation() 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; - - // If everything goes well, set the output +GenerateData() +{ + // Set the output this->GraftOutput(output); // not SetNthOutput } //-------------------------------------------------------------------- @@ -521,61 +446,167 @@ GenerateData() { //-------------------------------------------------------------------- template -void +bool 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()); +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; + } - // 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 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; } } - // 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); + } + return (m_Seeds.size() != 0); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLungFilter:: +TracheaRegionGrowing() +{ + // 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]); + // DD(m_Seeds[i]); + } + f->Update(); + + // take first (main) connected component + trachea_tmp = Labelize(f->GetOutput(), + GetBackgroundValue(), + true, + GetMinimalComponentSize()); + trachea_tmp = KeepLabels(trachea_tmp, + GetBackgroundValue(), + GetForegroundValue(), + 1, 1, false); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +double +clitk::ExtractLungFilter:: +ComputeTracheaVolume() +{ + typedef itk::ImageRegionConstIterator IteratorType; + IteratorType iter(trachea_tmp, trachea_tmp->GetLargestPossibleRegion()); + iter.GoToBegin(); + double volume = 0.0; + while (!iter.IsAtEnd()) { + if (iter.Get() == this->GetForegroundValue()) volume++; + ++iter; + } + + double voxelsize = trachea_tmp->GetSpacing()[0]*trachea_tmp->GetSpacing()[1]*trachea_tmp->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 stop = false; + double volume = 0.0; + int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed(); + while (!stop) { + stop = SearchForTracheaSeed(skip); + if (stop) { + 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; } else { - if (listOfTrackedPoint.size() > 2) { - std::cerr << "too much bifurcation points ... ?" << std::endl; - exit(0); - } - // Else this it the end of the tracking + 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); + } + else { // Trachea not found + this->SetWarning("* WARNING * No seed found for trachea."); + StopCurrentStep(); + } } //-------------------------------------------------------------------- - #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX