From e0464c9f7bdde185755d60e9dead188291fc3550 Mon Sep 17 00:00:00 2001 From: dsarrut Date: Thu, 22 Jul 2010 05:50:21 +0000 Subject: [PATCH] add draft bronchus decomposition --- segmentation/clitkExtractLungFilter.h | 91 +++--- segmentation/clitkExtractLungFilter.txx | 260 +++++++++++++----- .../clitkExtractLungGenericFilter.txx | 1 + .../clitkExtractLymphStationsFilter.txx | 4 +- 4 files changed, 257 insertions(+), 99 deletions(-) diff --git a/segmentation/clitkExtractLungFilter.h b/segmentation/clitkExtractLungFilter.h index 781a885..c62afe5 100644 --- a/segmentation/clitkExtractLungFilter.h +++ b/segmentation/clitkExtractLungFilter.h @@ -54,15 +54,36 @@ namespace clitk { */ //-------------------------------------------------------------------- - template + + //-------------------------------------------------------------------- +template +class Bifurcation +{ +public: + Bifurcation(IndexType _index, PixelType _l, PixelType _l1, PixelType _l2) { + index = _index; + _l = l; + _l1 = l1; + _l2 = l2; + } + IndexType index; + PixelType l; + PixelType l1; + PixelType l2; +}; + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template class ITK_EXPORT ExtractLungFilter: public clitk::FilterBase, - public itk::ImageToImageFilter + public itk::ImageToImageFilter { public: /** Standard class typedefs. */ - typedef itk::ImageToImageFilter Superclass; + typedef itk::ImageToImageFilter Superclass; typedef ExtractLungFilter Self; typedef itk::SmartPointer Pointer; typedef itk::SmartPointer ConstPointer; @@ -75,13 +96,13 @@ namespace clitk { FILTERBASE_INIT; /** Some convenient typedefs */ - typedef TInputImageType InputImageType; - typedef typename InputImageType::ConstPointer InputImageConstPointer; - typedef typename InputImageType::Pointer InputImagePointer; - typedef typename InputImageType::RegionType InputImageRegionType; - typedef typename InputImageType::PixelType InputImagePixelType; - typedef typename InputImageType::SizeType InputImageSizeType; - typedef typename InputImageType::IndexType InputImageIndexType; + typedef TImageType ImageType; + typedef typename ImageType::ConstPointer InputImageConstPointer; + typedef typename ImageType::Pointer InputImagePointer; + typedef typename ImageType::RegionType InputImageRegionType; + typedef typename ImageType::PixelType InputImagePixelType; + typedef typename ImageType::SizeType InputImageSizeType; + typedef typename ImageType::IndexType InputImageIndexType; typedef TMaskImageType MaskImageType; typedef typename MaskImageType::ConstPointer MaskImageConstPointer; @@ -91,15 +112,15 @@ namespace clitk { typedef typename MaskImageType::SizeType MaskImageSizeType; typedef typename MaskImageType::IndexType MaskImageIndexType; - itkStaticConstMacro(ImageDimension, unsigned int, InputImageType::ImageDimension); + itkStaticConstMacro(ImageDimension, unsigned int, ImageType::ImageDimension); typedef int InternalPixelType; - typedef itk::Image InternalImageType; - typedef typename InternalImageType::Pointer InternalImagePointer; - typedef typename InternalImageType::IndexType InternalIndexType; - typedef LabelizeParameters LabelParamType; + typedef itk::Image InternalImageType; + typedef typename InternalImageType::Pointer InternalImagePointer; + typedef typename InternalImageType::IndexType InternalIndexType; + typedef LabelizeParameters LabelParamType; /** Connect inputs */ - void SetInput(const InputImageType * image); + void SetInput(const ImageType * image); void SetInputPatientMask(MaskImageType * mask, MaskImagePixelType BG); itkSetMacro(PatientMaskBackgroundValue, MaskImagePixelType); itkGetConstMacro(PatientMaskBackgroundValue, MaskImagePixelType); @@ -152,7 +173,7 @@ namespace clitk { void AddSeed(InternalIndexType s); std::vector & GetSeeds() { return m_Seeds; } - GGO_DefineOption_Vector(seed, AddSeed, InternalIndexType, InputImageType::ImageDimension, true); + GGO_DefineOption_Vector(seed, AddSeed, InternalIndexType, ImageType::ImageDimension, true); // Step 3 options ExtractLung itkSetMacro(NumberOfHistogramBins, int); @@ -177,12 +198,25 @@ namespace clitk { // itkGetConstMacro(FinalOpenClose, bool); // itkBooleanMacro(FinalOpenClose); - // virtual void Update(); + // Bronchial bifurcations + itkSetMacro(FindBronchialBifurcations, bool); + itkGetConstMacro(FindBronchialBifurcations, bool); + itkBooleanMacro(FindBronchialBifurcations); protected: ExtractLungFilter(); virtual ~ExtractLungFilter() {} - + + // Main members + InputImageConstPointer input; + MaskImageConstPointer patient; + InputImagePointer working_input; + typename InternalImageType::Pointer working_image; + typename InternalImageType::Pointer trachea_tmp; + MaskImagePointer trachea; + MaskImagePointer output; + unsigned int m_MaxSeedNumber; + // Global options itkSetMacro(BackgroundValue, MaskImagePixelType); itkSetMacro(ForegroundValue, MaskImagePixelType); @@ -214,21 +248,16 @@ namespace clitk { // Step 5 // bool m_FinalOpenClose; + bool m_FindBronchialBifurcations; + virtual void GenerateOutputInformation(); virtual void GenerateData(); - // Steps - void RemoveAir(); - void FindTrachea(); - void ExtractLung(); - void RemoveTrachea(); - void LungSeparation(); - InputImageConstPointer input; - MaskImageConstPointer patient; - InputImagePointer working_input; - typename InternalImageType::Pointer working_image; - typename InternalImageType::Pointer trachea_tmp; - MaskImagePointer trachea; + typedef Bifurcation BifurcationType; + void TrackFromThisIndex(std::vector & listOfBifurcations, + MaskImagePointer skeleton, + MaskImageIndexType index, + MaskImagePixelType label); private: ExtractLungFilter(const Self&); //purposely not implemented diff --git a/segmentation/clitkExtractLungFilter.txx b/segmentation/clitkExtractLungFilter.txx index e4aa902..55e9e6c 100644 --- a/segmentation/clitkExtractLungFilter.txx +++ b/segmentation/clitkExtractLungFilter.txx @@ -30,20 +30,26 @@ #include "itkConnectedComponentImageFilter.h" #include "itkRelabelComponentImageFilter.h" #include "itkOtsuThresholdImageFilter.h" +#include "itkBinaryThinningImageFilter3D.h" +#include "itkImageIteratorWithIndex.h" + //-------------------------------------------------------------------- -template -clitk::ExtractLungFilter:: +template +clitk::ExtractLungFilter:: ExtractLungFilter(): clitk::FilterBase(), - itk::ImageToImageFilter() + itk::ImageToImageFilter() { SetNumberOfSteps(10); + m_MaxSeedNumber = 500; + // Default global options this->SetNumberOfRequiredInputs(2); SetPatientMaskBackgroundValue(0); SetBackgroundValue(0); // Must be zero SetForegroundValue(1); + SetMinimalComponentSize(100); // Step 1 default values SetUpperThreshold(-300); @@ -75,37 +81,40 @@ ExtractLungFilter(): p3->SetLastKeep(2); p3->UseLastKeepOff(); SetLabelizeParameters3(p3); + + // Step 5 : find bronchial bifurcations + FindBronchialBifurcationsOn(); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -template +template void -clitk::ExtractLungFilter:: -SetInput(const TInputImageType * image) +clitk::ExtractLungFilter:: +SetInput(const ImageType * image) { - this->SetNthInput(0, const_cast(image)); + this->SetNthInput(0, const_cast(image)); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -template +template void -clitk::ExtractLungFilter:: -SetInputPatientMask(TMaskImageType * image, MaskImagePixelType bg ) +clitk::ExtractLungFilter:: +SetInputPatientMask(MaskImageType * image, MaskImagePixelType bg ) { - this->SetNthInput(1, const_cast(image)); + this->SetNthInput(1, const_cast(image)); SetPatientMaskBackgroundValue(bg); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -template +template void -clitk::ExtractLungFilter:: +clitk::ExtractLungFilter:: AddSeed(InternalIndexType s) { m_Seeds.push_back(s); @@ -114,10 +123,10 @@ AddSeed(InternalIndexType s) //-------------------------------------------------------------------- -template +template template void -clitk::ExtractLungFilter:: +clitk::ExtractLungFilter:: SetArgsInfo(ArgsInfoType mArgsInfo) { SetVerboseOption_GGO(mArgsInfo); @@ -138,6 +147,8 @@ SetArgsInfo(ArgsInfoType mArgsInfo) SetThresholdStepSizeForTrachea_GGO(mArgsInfo); AddSeed_GGO(mArgsInfo); + SetMinimalComponentSize_GGO(mArgsInfo); + SetNumberOfHistogramBins_GGO(mArgsInfo); SetLabelizeParameters2_GGO(mArgsInfo); @@ -148,39 +159,36 @@ SetArgsInfo(ArgsInfoType mArgsInfo) //-------------------------------------------------------------------- -template +template void -clitk::ExtractLungFilter:: +clitk::ExtractLungFilter:: GenerateOutputInformation() { - input = dynamic_cast(itk::ProcessObject::GetInput(0)); - Superclass::GenerateOutputInformation(); + Superclass::GenerateOutputInformation(); // Needed ?? + this->GetOutput(0)->SetRequestedRegion(this->GetOutput(0)->GetLargestPossibleRegion()); // Get input pointers - input = dynamic_cast(itk::ProcessObject::GetInput(0)); - patient = dynamic_cast(itk::ProcessObject::GetInput(1)); + patient = dynamic_cast(itk::ProcessObject::GetInput(1)); + input = dynamic_cast(itk::ProcessObject::GetInput(0)); // Check image - if (!HaveSameSizeAndSpacing(input, patient)) { + if (!HaveSameSizeAndSpacing(input, patient)) { this->SetLastError("* ERROR * the images (input and patient mask) must have the same size & spacing"); return; } - // Set Number of steps - SetNumberOfSteps(9); - //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStep("Set background to initial image"); - working_input = SetBackground + StartNewStepOrStop("Set background to initial image"); + working_input = SetBackground (input, patient, GetPatientMaskBackgroundValue(), -1000); - StopCurrentStep(working_input); + StopCurrentStep(working_input); //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStep("Remove Air"); + StartNewStepOrStop("Remove Air"); // Threshold to get air - typedef itk::BinaryThresholdImageFilter InputBinarizeFilterType; + typedef itk::BinaryThresholdImageFilter InputBinarizeFilterType; typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New(); binarizeFilter->SetInput(working_input); if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(m_LowerThreshold); @@ -192,8 +200,10 @@ GenerateOutputInformation() // Labelize and keep right labels working_image = Labelize(working_image, GetBackgroundValue(), true, GetMinimalComponentSize()); + working_image = RemoveLabels (working_image, GetBackgroundValue(), GetLabelizeParameters1()->GetLabelsToRemove()); + typename InternalImageType::Pointer air = KeepLabels (working_image, GetBackgroundValue(), @@ -203,19 +213,18 @@ GenerateOutputInformation() GetLabelizeParameters1()->GetUseLastKeep()); // Set Air to BG - working_input = SetBackground + working_input = SetBackground (working_input, air, this->GetForegroundValue(), this->GetBackgroundValue()); // End - StopCurrentStep(working_input); + StopCurrentStep(working_input); //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStep("Find the trachea"); - //DD(m_Seeds.size()); + StartNewStepOrStop("Find the trachea"); if (m_Seeds.size() == 0) { // try to find seed // Search seed (parameters = UpperThresholdForTrachea) - static const unsigned int Dim = InputImageType::ImageDimension; + 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(); @@ -223,24 +232,20 @@ GenerateOutputInformation() sliceRegionSize[Dim-1]=5; sliceRegion.SetSize(sliceRegionSize); sliceRegion.SetIndex(sliceRegionIndex); - //DD(GetUpperThresholdForTrachea()); - //DD(sliceRegion); - typedef itk::ImageRegionConstIterator IteratorType; + 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; } } - //DD(m_Seeds.size()); if (m_Seeds.size() != 0) { // Explosion controlled region growing - typedef clitk::ExplosionControlledThresholdConnectedImageFilter ImageFilterType; + typedef clitk::ExplosionControlledThresholdConnectedImageFilter ImageFilterType; typename ImageFilterType::Pointer f= ImageFilterType::New(); f->SetInput(working_input); f->SetVerbose(false); @@ -256,26 +261,23 @@ GenerateOutputInformation() 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."); - // Set output StopCurrentStep(); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStep("Extract the lung with Otsu filter"); + StartNewStepOrStop("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()); @@ -289,7 +291,7 @@ GenerateOutputInformation() //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStep("Select labels"); + StartNewStepOrStop("Select labels"); // Keep right labels working_image = LabelizeAndSelectLabels (working_image, @@ -305,13 +307,13 @@ GenerateOutputInformation() //-------------------------------------------------------------------- //-------------------------------------------------------------------- if (m_Seeds.size() != 0) { // if ==0 ->no trachea found - StartNewStep("Remove the trachea"); + StartNewStepOrStop("Remove the trachea"); // Set the trachea working_image = SetBackground (working_image, trachea_tmp, 1, -1); // Dilate the trachea - static const unsigned int Dim = InputImageType::ImageDimension; + static const unsigned int Dim = ImageType::ImageDimension; typedef itk::BinaryBallStructuringElement KernelType; KernelType structuringElement; structuringElement.SetRadius(GetRadiusForTrachea()); @@ -346,14 +348,13 @@ GenerateOutputInformation() // Set output StopCurrentStep(working_image); } - //-------------------------------------------------------------------- //-------------------------------------------------------------------- typedef clitk::AutoCropFilter CropFilterType; typename CropFilterType::Pointer cropFilter = CropFilterType::New(); if (m_Seeds.size() != 0) { // if ==0 ->no trachea found - StartNewStep("Croping trachea"); + StartNewStepOrStop("Croping trachea"); cropFilter->SetInput(trachea_tmp); cropFilter->Update(); // Needed typedef itk::CastImageFilter CastImageFilterType; @@ -366,7 +367,7 @@ GenerateOutputInformation() //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStep("Croping lung"); + StartNewStepOrStop("Croping lung"); typename CropFilterType::Pointer cropFilter2 = CropFilterType::New(); // Needed to reset pipeline cropFilter2->SetInput(working_image); cropFilter2->Update(); @@ -375,7 +376,7 @@ GenerateOutputInformation() //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStep("Separate Left/Right lungs"); + StartNewStepOrStop("Separate Left/Right lungs"); // Initial label working_image = Labelize(working_image, GetBackgroundValue(), @@ -391,10 +392,10 @@ GenerateOutputInformation() working_image = statisticsImageFilter->GetOutput(); // Decompose the first label - static const unsigned int Dim = InputImageType::ImageDimension; + static const unsigned int Dim = ImageType::ImageDimension; if (initialNumberOfLabels<2) { // Structuring element radius - typename InputImageType::SizeType radius; + typename ImageType::SizeType radius; for (unsigned int i=0;i DecomposeAndReconstructFilterType; typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New(); @@ -421,24 +422,151 @@ GenerateOutputInformation() 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(); + + // 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 ! Abord"); + 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() { - // Final Cast - typedef itk::CastImageFilter CastImageFilterType; - typename CastImageFilterType::Pointer caster= CastImageFilterType::New(); - caster->SetInput(working_image); - caster->Update(); - // Set output - //this->SetNthOutput(0, caster->GetOutput()); // -> no because redo filter otherwise - this->GraftOutput(caster->GetOutput()); - return; + // 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 + this->GraftOutput(output); // not SetNthOutput +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +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; + } + } } //-------------------------------------------------------------------- diff --git a/segmentation/clitkExtractLungGenericFilter.txx b/segmentation/clitkExtractLungGenericFilter.txx index 77b4e8b..2fa40b1 100644 --- a/segmentation/clitkExtractLungGenericFilter.txx +++ b/segmentation/clitkExtractLungGenericFilter.txx @@ -84,6 +84,7 @@ void clitk::ExtractLungGenericFilter::UpdateWithInputImageType() this->SetFilterBase(filter); // Set global Options + filter->SetStopOnError(this->GetStopOnError()); filter->SetArgsInfo(mArgsInfo); filter->SetInput(input); filter->SetInputPatientMask(patient, mArgsInfo.patientBG_arg); diff --git a/segmentation/clitkExtractLymphStationsFilter.txx b/segmentation/clitkExtractLymphStationsFilter.txx index 5994f8f..8fbf86d 100644 --- a/segmentation/clitkExtractLymphStationsFilter.txx +++ b/segmentation/clitkExtractLymphStationsFilter.txx @@ -227,7 +227,7 @@ GenerateOutputInformation() { sliceRelPosFilter->SetInput(m_working_image); sliceRelPosFilter->SetInputObject(temp); sliceRelPosFilter->SetDirection(2); - sliceRelPosFilter->SetFuzzyThreshold(0.8); + sliceRelPosFilter->SetFuzzyThreshold(0.6); sliceRelPosFilter->SetOrientationType(SliceRelPosFilterType::RelPosFilterType::RightTo); sliceRelPosFilter->Update(); m_working_image = sliceRelPosFilter->GetOutput(); @@ -247,7 +247,7 @@ GenerateOutputInformation() { sliceRelPosFilter->SetInput(m_working_image); sliceRelPosFilter->SetInputObject(temp); sliceRelPosFilter->SetDirection(2); - sliceRelPosFilter->SetFuzzyThreshold(0.8); + sliceRelPosFilter->SetFuzzyThreshold(0.6); sliceRelPosFilter->SetOrientationType(SliceRelPosFilterType::RelPosFilterType::LeftTo); sliceRelPosFilter->Update(); m_working_image = sliceRelPosFilter->GetOutput(); -- 2.47.1