From 6e16222234a90c6079a8f4696c92de7349a496bb Mon Sep 17 00:00:00 2001 From: dsarrut Date: Wed, 8 Sep 2010 12:38:12 +0000 Subject: [PATCH] add smooth option to extract bones --- segmentation/CMakeLists.txt | 4 + segmentation/clitkExtractBones.ggo | 9 + segmentation/clitkExtractBonesFilter.h | 34 +- segmentation/clitkExtractBonesFilter.txx | 41 ++- segmentation/clitkExtractLung.ggo | 1 + segmentation/clitkExtractLungFilter.h | 11 + segmentation/clitkExtractLungFilter.txx | 338 ++++++++++++------ .../clitkExtractLymphStationsFilter.txx | 4 +- segmentation/clitkExtractMediastinum.ggo | 2 + segmentation/clitkExtractMediastinumFilter.h | 8 + .../clitkExtractMediastinumFilter.txx | 189 +++++++--- .../clitkExtractMediastinumGenericFilter.txx | 7 +- 12 files changed, 479 insertions(+), 169 deletions(-) diff --git a/segmentation/CMakeLists.txt b/segmentation/CMakeLists.txt index ceb2857..df23691 100644 --- a/segmentation/CMakeLists.txt +++ b/segmentation/CMakeLists.txt @@ -28,6 +28,10 @@ IF(CLITK_BUILD_SEGMENTATION) ADD_EXECUTABLE(clitkExtractLung clitkExtractLung.cxx ${clitkExtractLung_GGO_C}) TARGET_LINK_LIBRARIES(clitkExtractLung clitkCommon ITKIO) + WRAP_GGO(clitkExtractAirwayTreeInfo_GGO_C clitkExtractAirwayTreeInfo.ggo) + ADD_EXECUTABLE(clitkExtractAirwayTreeInfo clitkExtractAirwayTreeInfo.cxx ${clitkExtractAirwayTreeInfo_GGO_C}) + TARGET_LINK_LIBRARIES(clitkExtractAirwayTreeInfo clitkCommon ITKIO) + WRAP_GGO(clitkExtractBones_GGO_C clitkExtractBones.ggo) ADD_EXECUTABLE(clitkExtractBones clitkExtractBones.cxx ${clitkExtractBones_GGO_C}) TARGET_LINK_LIBRARIES(clitkExtractBones clitkCommon ITKIO) diff --git a/segmentation/clitkExtractBones.ggo b/segmentation/clitkExtractBones.ggo index 581392b..4938aa5 100644 --- a/segmentation/clitkExtractBones.ggo +++ b/segmentation/clitkExtractBones.ggo @@ -17,6 +17,15 @@ option "input" i "Input image filename" string yes option "output" o "Output image filename" string yes option "like" l "Resample like this image" string no +section "Smoothing (curvature anistropic diffusion)" + +option "smooth" - "smooth input image" flag off +option "spacing" - "Use image spacing" flag off +option "cond" - "Conductance parameter" double no default="3.0" +option "time" - "Time step (0.125 for 2D and 0.0625 for 3D)" double no default="0.0625" +option "iter" - "Iterations" double no default="5" + + section "Initial connected component labelling (CCL)" option "lower1" - "Lower threshold for CCL" double no default="100" diff --git a/segmentation/clitkExtractBonesFilter.h b/segmentation/clitkExtractBonesFilter.h index 0122d03..70daa48 100644 --- a/segmentation/clitkExtractBonesFilter.h +++ b/segmentation/clitkExtractBonesFilter.h @@ -95,6 +95,29 @@ namespace clitk { itkGetConstMacro(MinimalComponentSize, int); GGO_DefineOption(minSize, SetMinimalComponentSize, int); + // Step 0 + itkBooleanMacro(InitialSmoothing); + itkSetMacro(InitialSmoothing, bool); + itkGetMacro(InitialSmoothing, bool); + GGO_DefineOption_Flag(smooth, SetInitialSmoothing); + + itkSetMacro(SmoothingConductanceParameter, double); + itkGetConstMacro(SmoothingConductanceParameter, double); + GGO_DefineOption(cond, SetSmoothingConductanceParameter, double); + + itkSetMacro(SmoothingNumberOfIterations, int); + itkGetConstMacro(SmoothingNumberOfIterations, int); + GGO_DefineOption(iter, SetSmoothingNumberOfIterations, int); + + itkSetMacro(SmoothingTimeStep, double); + itkGetConstMacro(SmoothingTimeStep, double); + GGO_DefineOption(time, SetSmoothingTimeStep, double); + + itkSetMacro(SmoothingUseImageSpacing, bool); + itkGetConstMacro(SmoothingUseImageSpacing, bool); + itkBooleanMacro(SmoothingUseImageSpacing); + GGO_DefineOption_Flag(spacing, SetSmoothingUseImageSpacing); + // Step 1 itkSetMacro(UpperThreshold1, InputImagePixelType); itkGetMacro(UpperThreshold1, InputImagePixelType); @@ -141,6 +164,14 @@ namespace clitk { itkSetMacro(ForegroundValue, OutputImagePixelType); OutputImagePixelType m_BackgroundValue; OutputImagePixelType m_ForegroundValue; + bool m_AutoCrop; + + // Step 0 : Initial Filtering + bool m_InitialSmoothing; + double m_SmoothingConductanceParameter; + int m_SmoothingNumberOfIterations; + double m_SmoothingTimeStep; + bool m_SmoothingUseImageSpacing; // Step 1 InputImagePixelType m_UpperThreshold1; @@ -154,8 +185,6 @@ namespace clitk { InputImageSizeType m_Radius2; int m_SampleRate2; - bool m_AutoCrop; - virtual void GenerateOutputInformation(); virtual void GenerateData(); @@ -166,6 +195,7 @@ namespace clitk { void RemoveTrachea(); void BonesSeparation(); InputImageConstPointer input; + InputImagePointer filtered_input; OutputImageConstPointer patient; InputImagePointer working_input; typename InternalImageType::Pointer working_image; diff --git a/segmentation/clitkExtractBonesFilter.txx b/segmentation/clitkExtractBonesFilter.txx index 00f72f2..ef3f944 100644 --- a/segmentation/clitkExtractBonesFilter.txx +++ b/segmentation/clitkExtractBonesFilter.txx @@ -30,6 +30,7 @@ #include "itkConnectedComponentImageFilter.h" #include "itkRelabelComponentImageFilter.h" #include "itkNeighborhoodConnectedImageFilter.h" +#include "itkCurvatureAnisotropicDiffusionImageFilter.h" //-------------------------------------------------------------------- template @@ -42,6 +43,12 @@ ExtractBonesFilter(): this->SetNumberOfRequiredInputs(1); SetBackgroundValue(0); // Must be zero SetForegroundValue(1); + SetInitialSmoothing(false); + + SetSmoothingConductanceParameter(3.0); + SetSmoothingNumberOfIterations(5); + SetSmoothingTimeStep(0.0625); + SetSmoothingUseImageSpacing(false); SetMinimalComponentSize(100); SetUpperThreshold1(1500); @@ -81,6 +88,12 @@ SetArgsInfo(ArgsInfoType mArgsInfo) SetVerboseStep_GGO(mArgsInfo); SetWriteStep_GGO(mArgsInfo); SetVerboseWarningOff_GGO(mArgsInfo); + + SetInitialSmoothing_GGO(mArgsInfo); + SetSmoothingConductanceParameter_GGO(mArgsInfo); + SetSmoothingNumberOfIterations_GGO(mArgsInfo); + SetSmoothingTimeStep_GGO(mArgsInfo); + SetSmoothingUseImageSpacing_GGO(mArgsInfo); SetMinimalComponentSize_GGO(mArgsInfo); SetUpperThreshold1_GGO(mArgsInfo); @@ -117,17 +130,36 @@ GenerateOutputInformation() { typedef itk::CastImageFilter CastImageFilterType; typedef itk::ImageFileWriter WriterType; + //--------------------------------- + // Smoothing [Optional] + //--------------------------------- + if (GetInitialSmoothing()) { + StartNewStep("Initial Smoothing"); + typedef itk::CurvatureAnisotropicDiffusionImageFilter FilterType; + typename FilterType::Pointer df = FilterType::New(); + df->SetConductanceParameter(GetSmoothingConductanceParameter()); + df->SetNumberOfIterations(GetSmoothingNumberOfIterations()); + df->SetTimeStep(GetSmoothingTimeStep()); + df->SetUseImageSpacing(GetSmoothingUseImageSpacing()); + df->SetInput(input); + df->Update(); + filtered_input = df->GetOutput(); + StopCurrentStep(filtered_input); + } + else { + filtered_input = input; + } + //-------------------------------------------------------------------- //-------------------------------------------------------------------- StartNewStep("Initial Labeling"); - typename InternalImageType::Pointer firstLabelImage; //--------------------------------- // Binarize the image //--------------------------------- typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New(); - binarizeFilter->SetInput(input); + binarizeFilter->SetInput(filtered_input); binarizeFilter->SetLowerThreshold(GetLowerThreshold1()); binarizeFilter->SetUpperThreshold(GetUpperThreshold1()); binarizeFilter->SetInsideValue(this->GetForegroundValue()); @@ -160,6 +192,7 @@ GenerateOutputInformation() { binarizeFilter2->Update(); firstLabelImage = binarizeFilter2->GetOutput(); + StopCurrentStep(firstLabelImage); //-------------------------------------------------------------------- //-------------------------------------------------------------------- @@ -179,7 +212,7 @@ GenerateOutputInformation() { neighborhoodConnectedImageFilter->SetUpper(GetUpperThreshold2()); neighborhoodConnectedImageFilter->SetReplaceValue(this->GetForegroundValue()); neighborhoodConnectedImageFilter->SetRadius(GetRadius2()); - neighborhoodConnectedImageFilter->SetInput(input); + neighborhoodConnectedImageFilter->SetInput(filtered_input); // Seeds from label image typedef itk::ImageRegionIteratorWithIndex IteratorType; @@ -209,7 +242,7 @@ GenerateOutputInformation() { //-------------------------------------------------------------------- //-------------------------------------------------------------------- - StartNewStep("Combine de images"); + StartNewStep("Combine the images"); typedef clitk::SetBackgroundImageFilter SetBackgroundImageFilterType; typename SetBackgroundImageFilterType::Pointer setBackgroundFilter=SetBackgroundImageFilterType::New(); diff --git a/segmentation/clitkExtractLung.ggo b/segmentation/clitkExtractLung.ggo index 883a013..08f066a 100644 --- a/segmentation/clitkExtractLung.ggo +++ b/segmentation/clitkExtractLung.ggo @@ -30,6 +30,7 @@ option "lastKeep1" - "Last label to keep" int no section "Step 2 : find trachea" +option "skipslices" - "Number of slices to skip before searching seed" int no default="0" option "upperThresholdForTrachea" - "Initial upper threshold for trachea" double no default="-900" option "multiplierForTrachea" - "Multiplier for the region growing" double no default="5" option "thresholdStepSizeForTrachea" - "Threshold step size" int no default="64" diff --git a/segmentation/clitkExtractLungFilter.h b/segmentation/clitkExtractLungFilter.h index c62afe5..6a316d2 100644 --- a/segmentation/clitkExtractLungFilter.h +++ b/segmentation/clitkExtractLungFilter.h @@ -147,6 +147,10 @@ public: itkGetConstMacro(UpperThreshold, InputImagePixelType); GGO_DefineOption(upper, SetUpperThreshold, InputImagePixelType); + itkSetMacro(NumberOfSlicesToSkipBeforeSearchingSeed, int); + itkGetConstMacro(NumberOfSlicesToSkipBeforeSearchingSeed, int); + GGO_DefineOption(skipslices, SetNumberOfSlicesToSkipBeforeSearchingSeed, int); + itkSetMacro(LowerThreshold, InputImagePixelType); itkGetConstMacro(LowerThreshold, InputImagePixelType); itkSetMacro(UseLowerThreshold, bool); @@ -236,6 +240,7 @@ public: InputImagePixelType m_ThresholdStepSizeForTrachea; double m_MultiplierForTrachea; std::vector m_Seeds; + int m_NumberOfSlicesToSkipBeforeSearchingSeed; // Step 3 int m_NumberOfHistogramBins; @@ -259,6 +264,12 @@ public: MaskImageIndexType index, MaskImagePixelType label); + + bool SearchForTracheaSeed(int skip); + void SearchForTrachea(); + void TracheaRegionGrowing(); + double ComputeTracheaVolume(); + private: ExtractLungFilter(const Self&); //purposely not implemented void operator=(const Self&); //purposely not implemented diff --git a/segmentation/clitkExtractLungFilter.txx b/segmentation/clitkExtractLungFilter.txx index 4baee45..3bb57b8 100644 --- a/segmentation/clitkExtractLungFilter.txx +++ b/segmentation/clitkExtractLungFilter.txx @@ -65,6 +65,7 @@ ExtractLungFilter(): SetUpperThresholdForTrachea(-900); SetMultiplierForTrachea(5); SetThresholdStepSizeForTrachea(64); + SetNumberOfSlicesToSkipBeforeSearchingSeed(0); // Step 3 default values SetNumberOfHistogramBins(500); @@ -136,6 +137,7 @@ SetArgsInfo(ArgsInfoType mArgsInfo) SetUpperThreshold_GGO(mArgsInfo); SetLowerThreshold_GGO(mArgsInfo); + SetNumberOfSlicesToSkipBeforeSearchingSeed_GGO(mArgsInfo); SetLabelizeParameters1_GGO(mArgsInfo); if (!mArgsInfo.remove1_given) { GetLabelizeParameters1()->AddLabelToRemove(2); @@ -228,57 +230,8 @@ 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(); - } + StartNewStepOrStop("Search for the trachea"); + SearchForTrachea(); //-------------------------------------------------------------------- //-------------------------------------------------------------------- @@ -319,7 +272,7 @@ GenerateOutputInformation() 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; @@ -420,7 +373,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); @@ -443,61 +396,64 @@ GenerateOutputInformation() // 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"); + 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.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()); + // 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); - } + // 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,7 +463,8 @@ GenerateOutputInformation() template void clitk::ExtractLungFilter:: -GenerateData() { +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. @@ -526,7 +483,8 @@ clitk::ExtractLungFilter:: TrackFromThisIndex(std::vector & listOfBifurcations, MaskImagePointer skeleton, MaskImageIndexType index, - MaskImagePixelType label) { + MaskImagePixelType label) +{ DD("TrackFromThisIndex"); DD(index); DD((int)label); @@ -577,5 +535,173 @@ TrackFromThisIndex(std::vector & listOfBifurcations, } //-------------------------------------------------------------------- + +//-------------------------------------------------------------------- +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 (GetVerboseStep()) { + 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; + } + } + } + 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 + 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 +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(); // assume mm3 + if ((volume > 10000) && (volume < 55000 )) { // 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 (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(); + } + } + } + + 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 diff --git a/segmentation/clitkExtractLymphStationsFilter.txx b/segmentation/clitkExtractLymphStationsFilter.txx index 8fbf86d..aab78d8 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.6); + sliceRelPosFilter->SetFuzzyThreshold(0.5); 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.6); + sliceRelPosFilter->SetFuzzyThreshold(0.5); sliceRelPosFilter->SetOrientationType(SliceRelPosFilterType::RelPosFilterType::LeftTo); sliceRelPosFilter->Update(); m_working_image = sliceRelPosFilter->GetOutput(); diff --git a/segmentation/clitkExtractMediastinum.ggo b/segmentation/clitkExtractMediastinum.ggo index 490c673..2f57cae 100644 --- a/segmentation/clitkExtractMediastinum.ggo +++ b/segmentation/clitkExtractMediastinum.ggo @@ -21,6 +21,8 @@ option "lung" l "Input lung mask filename" string yes option "lungBG" - "Lung Background" int default="0" no option "lungLeft" - "Lung Left value" int default="1" no option "lungRight" - "Lung Right value" int default="2" no +option "trachea" t "Input trachea mask filename" string yes +option "tracheaBG" - "Trachea Background" int default="0" no option "output" o "Output lungs mask filename" string yes diff --git a/segmentation/clitkExtractMediastinumFilter.h b/segmentation/clitkExtractMediastinumFilter.h index 6c252d4..f02833d 100644 --- a/segmentation/clitkExtractMediastinumFilter.h +++ b/segmentation/clitkExtractMediastinumFilter.h @@ -30,6 +30,7 @@ namespace clitk { - Patient label image - Lungs label image - Bones label image + - Trachea label image */ //-------------------------------------------------------------------- @@ -61,12 +62,14 @@ namespace clitk { typedef typename ImageType::PixelType ImagePixelType; typedef typename ImageType::SizeType ImageSizeType; typedef typename ImageType::IndexType ImageIndexType; + typedef typename ImageType::PointType ImagePointType; /** Connect inputs */ void SetInputPatientLabelImage(const TImageType * image, ImagePixelType bg=0); void SetInputLungLabelImage(const TImageType * image, ImagePixelType bg=0, ImagePixelType fgLeftLung=1, ImagePixelType fgRightLung=2); void SetInputBonesLabelImage(const TImageType * image, ImagePixelType bg=0); + void SetInputTracheaLabelImage(const TImageType * image, ImagePixelType bg=0); /** ImageDimension constants */ itkStaticConstMacro(ImageDimension, unsigned int, TImageType::ImageDimension); @@ -99,6 +102,10 @@ namespace clitk { itkGetConstMacro(ForegroundValueRightLung, ImagePixelType); GGO_DefineOption(lungRight, SetForegroundValueRightLung, ImagePixelType); + itkSetMacro(BackgroundValueTrachea, ImagePixelType); + itkGetConstMacro(BackgroundValueTrachea, ImagePixelType); + GGO_DefineOption(lungBG, SetBackgroundValueTrachea, ImagePixelType); + itkSetMacro(IntermediateSpacing, double); itkGetConstMacro(IntermediateSpacing, double); GGO_DefineOption(spacing, SetIntermediateSpacing, double); @@ -125,6 +132,7 @@ namespace clitk { ImagePixelType m_BackgroundValuePatient; ImagePixelType m_BackgroundValueLung; ImagePixelType m_BackgroundValueBones; + ImagePixelType m_BackgroundValueTrachea; ImagePixelType m_ForegroundValueLeftLung; ImagePixelType m_ForegroundValueRightLung; diff --git a/segmentation/clitkExtractMediastinumFilter.txx b/segmentation/clitkExtractMediastinumFilter.txx index 854a804..fde335c 100644 --- a/segmentation/clitkExtractMediastinumFilter.txx +++ b/segmentation/clitkExtractMediastinumFilter.txx @@ -24,6 +24,7 @@ #include "clitkExtractMediastinumFilter.h" #include "clitkAddRelativePositionConstraintToLabelImageFilter.h" #include "clitkSegmentationUtils.h" +#include "clitkExtractAirwayTreeInfoFilter.h" // itk #include @@ -36,13 +37,13 @@ #include "RelativePositionPropImageFilter.h" //-------------------------------------------------------------------- -template -clitk::ExtractMediastinumFilter:: +template +clitk::ExtractMediastinumFilter:: ExtractMediastinumFilter(): clitk::FilterBase(), - itk::ImageToImageFilter() + itk::ImageToImageFilter() { - this->SetNumberOfRequiredInputs(3); + this->SetNumberOfRequiredInputs(4); SetBackgroundValuePatient(0); SetBackgroundValueLung(0); @@ -60,23 +61,25 @@ ExtractMediastinumFilter(): //-------------------------------------------------------------------- -template +template void -clitk::ExtractMediastinumFilter:: -SetInputPatientLabelImage(const TImageType * image, ImagePixelType bg) { - this->SetNthInput(0, const_cast(image)); +clitk::ExtractMediastinumFilter:: +SetInputPatientLabelImage(const ImageType * image, ImagePixelType bg) +{ + this->SetNthInput(0, const_cast(image)); m_BackgroundValuePatient = bg; } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -template +template void -clitk::ExtractMediastinumFilter:: -SetInputLungLabelImage(const TImageType * image, ImagePixelType bg, - ImagePixelType fgLeft, ImagePixelType fgRight) { - this->SetNthInput(1, const_cast(image)); +clitk::ExtractMediastinumFilter:: +SetInputLungLabelImage(const ImageType * image, ImagePixelType bg, + ImagePixelType fgLeft, ImagePixelType fgRight) +{ + this->SetNthInput(1, const_cast(image)); m_BackgroundValueLung = bg; m_ForegroundValueLeftLung = fgLeft; m_ForegroundValueRightLung = fgRight; @@ -85,21 +88,34 @@ SetInputLungLabelImage(const TImageType * image, ImagePixelType bg, //-------------------------------------------------------------------- -template +template void -clitk::ExtractMediastinumFilter:: -SetInputBonesLabelImage(const TImageType * image, ImagePixelType bg) { - this->SetNthInput(2, const_cast(image)); +clitk::ExtractMediastinumFilter:: +SetInputBonesLabelImage(const ImageType * image, ImagePixelType bg) +{ + this->SetNthInput(2, const_cast(image)); m_BackgroundValueBones = bg; } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -template +template +void +clitk::ExtractMediastinumFilter:: +SetInputTracheaLabelImage(const ImageType * image, ImagePixelType bg) +{ + this->SetNthInput(3, const_cast(image)); + m_BackgroundValueTrachea = bg; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template template void -clitk::ExtractMediastinumFilter:: +clitk::ExtractMediastinumFilter:: SetArgsInfo(ArgsInfoType mArgsInfo) { SetVerboseOption_GGO(mArgsInfo); @@ -109,7 +125,7 @@ SetArgsInfo(ArgsInfoType mArgsInfo) SetBackgroundValuePatient_GGO(mArgsInfo); SetBackgroundValueLung_GGO(mArgsInfo); - SetBackgroundValueBones_GGO(mArgsInfo); + SetBackgroundValueTrachea_GGO(mArgsInfo); SetForegroundValueLeftLung_GGO(mArgsInfo); SetForegroundValueRightLung_GGO(mArgsInfo); @@ -122,11 +138,11 @@ SetArgsInfo(ArgsInfoType mArgsInfo) //-------------------------------------------------------------------- -template +template void -clitk::ExtractMediastinumFilter:: +clitk::ExtractMediastinumFilter:: GenerateOutputInformation() { - ImagePointer input = dynamic_cast(itk::ProcessObject::GetInput(0)); + ImagePointer input = dynamic_cast(itk::ProcessObject::GetInput(0)); ImagePointer outputImage = this->GetOutput(0); outputImage->SetRegions(outputImage->GetLargestPossibleRegion()); } @@ -134,33 +150,38 @@ GenerateOutputInformation() { //-------------------------------------------------------------------- -template +template void -clitk::ExtractMediastinumFilter:: -GenerateInputRequestedRegion() { +clitk::ExtractMediastinumFilter:: +GenerateInputRequestedRegion() +{ // Call default Superclass::GenerateInputRequestedRegion(); // Get input pointers - ImagePointer patient = dynamic_cast(itk::ProcessObject::GetInput(0)); - ImagePointer lung = dynamic_cast(itk::ProcessObject::GetInput(1)); - ImagePointer bones = dynamic_cast(itk::ProcessObject::GetInput(2)); + ImagePointer patient = dynamic_cast(itk::ProcessObject::GetInput(0)); + ImagePointer lung = dynamic_cast(itk::ProcessObject::GetInput(1)); + ImagePointer bones = dynamic_cast(itk::ProcessObject::GetInput(2)); + ImagePointer trachea = dynamic_cast(itk::ProcessObject::GetInput(3)); patient->SetRequestedRegion(patient->GetLargestPossibleRegion()); lung->SetRequestedRegion(lung->GetLargestPossibleRegion()); bones->SetRequestedRegion(bones->GetLargestPossibleRegion()); + trachea->SetRequestedRegion(trachea->GetLargestPossibleRegion()); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -template +template void -clitk::ExtractMediastinumFilter:: -GenerateData() { +clitk::ExtractMediastinumFilter:: +GenerateData() +{ // Get input pointers - ImageConstPointer patient = dynamic_cast(itk::ProcessObject::GetInput(0)); - ImageConstPointer lung = dynamic_cast(itk::ProcessObject::GetInput(1)); - ImageConstPointer bones = dynamic_cast(itk::ProcessObject::GetInput(2)); + ImageConstPointer patient = dynamic_cast(itk::ProcessObject::GetInput(0)); + ImageConstPointer lung = dynamic_cast(itk::ProcessObject::GetInput(1)); + ImageConstPointer bones = dynamic_cast(itk::ProcessObject::GetInput(2)); + ImageConstPointer trachea = dynamic_cast(itk::ProcessObject::GetInput(3)); // Get output pointer ImagePointer output; @@ -180,73 +201,135 @@ GenerateData() { boolFilter->Update(); output = boolFilter->GetOutput(); + // Auto crop to gain support area output = clitk::AutoCrop(output, GetBackgroundValue()); - ////autoCropFilter->GetOutput(); typedef clitk::AutoCropFilter CropFilterType; - //typename CropFilterType::Pointer cropFilter = CropFilterType::New(); - //cropFilter->SetInput(output); - //cropFilter->Update(); - //output = cropFilter->GetOutput(); - - this->template StopCurrentStep(output); + this->template StopCurrentStep(output); // Step 2: LR limits from lung (need separate lung ?) StartNewStep("Left/Right limits with lungs"); - typedef clitk::AddRelativePositionConstraintToLabelImageFilter RelPosFilterType; + + /* // WE DO NOT NEED THE FOLLOWING + // Get separate lung images to get only the right and left lung + // (label must be '1' because right is greater than left). + ImagePointer right_lung = clitk::SetBackground(lung, lung, 2, 0); + ImagePointer left_lung = clitk::SetBackground(lung, lung, 1, 0); + writeImage(right_lung, "right.mhd"); + writeImage(left_lung, "left.mhd"); + */ + + typedef clitk::AddRelativePositionConstraintToLabelImageFilter RelPosFilterType; typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New(); relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId()); relPosFilter->VerboseStepOff(); relPosFilter->WriteStepOff(); relPosFilter->SetInput(output); + DD(output->GetLargestPossibleRegion().GetIndex()); + // relPosFilter->SetInputObject(left_lung); relPosFilter->SetInputObject(lung); - relPosFilter->SetOrientationType(RelPosFilterType::LeftTo); + relPosFilter->SetOrientationType(RelPosFilterType::LeftTo); // warning left lung is at right ;) relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing()); relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1()); relPosFilter->Update(); output = relPosFilter->GetOutput(); + DD(output->GetLargestPossibleRegion()); relPosFilter->SetInput(output); relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId()); relPosFilter->VerboseStepOff(); relPosFilter->WriteStepOff(); + // relPosFilter->SetInputObject(right_lung); relPosFilter->SetInputObject(lung); relPosFilter->SetOrientationType(RelPosFilterType::RightTo); relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing()); relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold1()); relPosFilter->Update(); output = relPosFilter->GetOutput(); - this->template StopCurrentStep(output); + DD(output->GetLargestPossibleRegion()); + this->template StopCurrentStep(output); // Step 3: AP limits from bones StartNewStep("Ant/Post limits with bones"); + ImageConstPointer bones_ant; + ImageConstPointer bones_post; + + // Find ant-post coordinate of trachea (assume the first point is a + // good estimation of the ant-post position of the trachea) + typedef clitk::ExtractAirwayTreeInfoFilter AirwayFilter; + typename AirwayFilter::Pointer airwayfilter = AirwayFilter::New(); + airwayfilter->SetInput(trachea); + airwayfilter->Update(); + DD(airwayfilter->GetFirstTracheaPoint()); + ImagePointType point_trachea = airwayfilter->GetFirstTracheaPoint(); + ImageIndexType index_trachea; + bones->TransformPhysicalPointToIndex(point_trachea, index_trachea); + DD(index_trachea); + + // Split bone image first into two parts (ant and post) + typedef itk::RegionOfInterestImageFilter CropFilterType; + // typedef itk::ExtractImageFilter CropFilterType; + typename CropFilterType::Pointer cropFilter = CropFilterType::New(); + ImageRegionType region = bones->GetLargestPossibleRegion(); + ImageSizeType size = region.GetSize(); + DD(size); + size[1] = index_trachea[1]; //size[1]/2.0; + DD(size); + region.SetSize(size); + cropFilter->SetInput(bones); + // cropFilter->SetExtractionRegion(region); + cropFilter->SetRegionOfInterest(region); + cropFilter->ReleaseDataFlagOff(); + cropFilter->Update(); + bones_ant = cropFilter->GetOutput(); + writeImage(bones_ant, "b_ant.mhd"); + + // cropFilter->ResetPipeline();// = CropFilterType::New(); + cropFilter = CropFilterType::New(); + ImageIndexType index = region.GetIndex(); + size[1] = bones->GetLargestPossibleRegion().GetSize()[1] - size[1]; + index[1] = bones->GetLargestPossibleRegion().GetIndex()[1] + size[1]-1; + DD(size); + region.SetIndex(index); + region.SetSize(size); + cropFilter->SetInput(bones); + // cropFilter->SetExtractionRegion(region); + cropFilter->SetRegionOfInterest(region); + cropFilter->ReleaseDataFlagOff(); + cropFilter->Update(); + bones_post = cropFilter->GetOutput(); + writeImage(bones_post, "b_post.mhd"); + + // Go ! relPosFilter->SetCurrentStepNumber(0); relPosFilter->ResetPipeline();// = RelPosFilterType::New(); relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId()); relPosFilter->VerboseStepOff(); relPosFilter->WriteStepOff(); relPosFilter->SetInput(output); - relPosFilter->SetInputObject(bones); + relPosFilter->SetInputObject(bones_post); relPosFilter->SetOrientationType(RelPosFilterType::AntTo); relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing()); relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold2()); relPosFilter->Update(); + output = relPosFilter->GetOutput(); + writeImage(output, "post.mhd"); relPosFilter->SetInput(relPosFilter->GetOutput()); relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId()); relPosFilter->VerboseStepOff(); relPosFilter->WriteStepOff(); - relPosFilter->SetInputObject(bones); + relPosFilter->SetInput(output); + relPosFilter->SetInputObject(bones_ant); relPosFilter->SetOrientationType(RelPosFilterType::PostTo); relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing()); relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold2()); relPosFilter->Update(); output = relPosFilter->GetOutput(); - this->template StopCurrentStep(output); - + this->template StopCurrentStep(output); // Get CCL - output = clitk::Labelize(output, GetBackgroundValue(), true, 100); - // output = RemoveLabels(output, BG, param->GetLabelsToRemove()); - output = clitk::KeepLabels(output, GetBackgroundValue(), - GetForegroundValue(), 1, 1, 0); + output = clitk::Labelize(output, GetBackgroundValue(), true, 100); + // output = RemoveLabels(output, BG, param->GetLabelsToRemove()); + output = clitk::KeepLabels(output, GetBackgroundValue(), + GetForegroundValue(), 1, 1, 0); output = clitk::AutoCrop(output, GetBackgroundValue()); // cropFilter = CropFilterType::New(); diff --git a/segmentation/clitkExtractMediastinumGenericFilter.txx b/segmentation/clitkExtractMediastinumGenericFilter.txx index 5a2cead..8632222 100644 --- a/segmentation/clitkExtractMediastinumGenericFilter.txx +++ b/segmentation/clitkExtractMediastinumGenericFilter.txx @@ -56,6 +56,7 @@ void clitk::ExtractMediastinumGenericFilter::SetArgsInfo(const Arg if (mArgsInfo.patient_given) AddInputFilename(mArgsInfo.patient_arg); if (mArgsInfo.lung_given) AddInputFilename(mArgsInfo.lung_arg); if (mArgsInfo.bones_given) AddInputFilename(mArgsInfo.bones_arg); + if (mArgsInfo.trachea_given) AddInputFilename(mArgsInfo.trachea_arg); if (mArgsInfo.output_given) AddOutputFilename(mArgsInfo.output_arg); } //-------------------------------------------------------------------- @@ -70,8 +71,9 @@ void clitk::ExtractMediastinumGenericFilter::UpdateWithInputImageT { // Reading input typename ImageType::Pointer patient = this->template GetInput(0); - typename ImageType::Pointer lung = this->template GetInput(1); - typename ImageType::Pointer bones = this->template GetInput(2); + typename ImageType::Pointer lung = this->template GetInput(1); + typename ImageType::Pointer bones = this->template GetInput(2); + typename ImageType::Pointer trachea = this->template GetInput(3); // Create filter typedef clitk::ExtractMediastinumFilter FilterType; @@ -81,6 +83,7 @@ void clitk::ExtractMediastinumGenericFilter::UpdateWithInputImageT filter->SetInputPatientLabelImage(patient, mArgsInfo.patientBG_arg); filter->SetInputLungLabelImage(lung, mArgsInfo.lungBG_arg, mArgsInfo.lungRight_arg, mArgsInfo.lungLeft_arg); filter->SetInputBonesLabelImage(bones, mArgsInfo.bonesBG_arg); + filter->SetInputTracheaLabelImage(trachea, mArgsInfo.tracheaBG_arg); filter->SetArgsInfo(mArgsInfo); // Go ! -- 2.45.1