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)
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"
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);
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;
InputImageSizeType m_Radius2;
int m_SampleRate2;
- bool m_AutoCrop;
-
virtual void GenerateOutputInformation();
virtual void GenerateData();
void RemoveTrachea();
void BonesSeparation();
InputImageConstPointer input;
+ InputImagePointer filtered_input;
OutputImageConstPointer patient;
InputImagePointer working_input;
typename InternalImageType::Pointer working_image;
#include "itkConnectedComponentImageFilter.h"
#include "itkRelabelComponentImageFilter.h"
#include "itkNeighborhoodConnectedImageFilter.h"
+#include "itkCurvatureAnisotropicDiffusionImageFilter.h"
//--------------------------------------------------------------------
template <class TInputImageType, class TOutputImageType>
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);
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);
typedef itk::CastImageFilter<InternalImageType,OutputImageType> CastImageFilterType;
typedef itk::ImageFileWriter<OutputImageType> WriterType;
+ //---------------------------------
+ // Smoothing [Optional]
+ //---------------------------------
+ if (GetInitialSmoothing()) {
+ StartNewStep("Initial Smoothing");
+ typedef itk::CurvatureAnisotropicDiffusionImageFilter<InputImageType, InputImageType> 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<InputImageType>(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());
binarizeFilter2->Update();
firstLabelImage = binarizeFilter2->GetOutput();
+ StopCurrentStep<InternalImageType>(firstLabelImage);
//--------------------------------------------------------------------
//--------------------------------------------------------------------
neighborhoodConnectedImageFilter->SetUpper(GetUpperThreshold2());
neighborhoodConnectedImageFilter->SetReplaceValue(this->GetForegroundValue());
neighborhoodConnectedImageFilter->SetRadius(GetRadius2());
- neighborhoodConnectedImageFilter->SetInput(input);
+ neighborhoodConnectedImageFilter->SetInput(filtered_input);
// Seeds from label image
typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
//--------------------------------------------------------------------
//--------------------------------------------------------------------
- StartNewStep("Combine de images");
+ StartNewStep("Combine the images");
typedef clitk::SetBackgroundImageFilter<InternalImageType, InternalImageType, InternalImageType>
SetBackgroundImageFilterType;
typename SetBackgroundImageFilterType::Pointer setBackgroundFilter=SetBackgroundImageFilterType::New();
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"
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);
InputImagePixelType m_ThresholdStepSizeForTrachea;
double m_MultiplierForTrachea;
std::vector<InternalIndexType> m_Seeds;
+ int m_NumberOfSlicesToSkipBeforeSearchingSeed;
// Step 3
int m_NumberOfHistogramBins;
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
SetUpperThresholdForTrachea(-900);
SetMultiplierForTrachea(5);
SetThresholdStepSizeForTrachea(64);
+ SetNumberOfSlicesToSkipBeforeSearchingSeed(0);
// Step 3 default values
SetNumberOfHistogramBins(500);
SetUpperThreshold_GGO(mArgsInfo);
SetLowerThreshold_GGO(mArgsInfo);
+ SetNumberOfSlicesToSkipBeforeSearchingSeed_GGO(mArgsInfo);
SetLabelizeParameters1_GGO(mArgsInfo);
if (!mArgsInfo.remove1_given) {
GetLabelizeParameters1()->AddLabelToRemove(2);
//--------------------------------------------------------------------
//--------------------------------------------------------------------
- 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<ImageType> 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<ImageType, InternalImageType> 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; i<m_Seeds.size();i++) {
- f->AddSeed(m_Seeds[i]);
- }
- f->Update();
- trachea_tmp = f->GetOutput();
- // Set output
- StopCurrentStep<InternalImageType>(trachea_tmp);
- }
- else { // Trachea not found
- this->SetWarning("* WARNING * No seed found for trachea.");
- StopCurrentStep();
- }
+ StartNewStepOrStop("Search for the trachea");
+ SearchForTrachea();
//--------------------------------------------------------------------
//--------------------------------------------------------------------
working_image = SetBackground<InternalImageType, InternalImageType>
(working_image, trachea_tmp, 1, -1);
- // Dilate the trachea
+ // Dilate the trachea
static const unsigned int Dim = ImageType::ImageDimension;
typedef itk::BinaryBallStructuringElement<InternalPixelType, Dim> KernelType;
KernelType structuringElement;
working_image = decomposeAndReconstructFilter->GetOutput();
}
- // Retain labels (lungs)
+ // Retain labels ('1' is largset lung, so right. '2' is left)
typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
thresholdFilter->SetInput(working_image);
// 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<MaskImageType, MaskImageType> ThinningFilterType;
- typename ThinningFilterType::Pointer thinningFilter = ThinningFilterType::New();
- thinningFilter->SetInput(trachea);
- thinningFilter->Update();
- typename MaskImageType::Pointer skeleton = thinningFilter->GetOutput();
- writeImage<MaskImageType>(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<MaskImageType, MaskImageType> ThinningFilterType;
+ typename ThinningFilterType::Pointer thinningFilter = ThinningFilterType::New();
+ thinningFilter->SetInput(trachea);
+ thinningFilter->Update();
+ typename MaskImageType::Pointer skeleton = thinningFilter->GetOutput();
+ writeImage<MaskImageType>(skeleton, "skeleton.mhd");
+
+ // Step 2 : tracking
+ DD("tracking");
- // Step 2.1 : find first point for tracking
- typedef itk::ImageRegionConstIteratorWithIndex<MaskImageType> 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<MaskImageType> 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<MaskImageType> 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<MaskImageType> 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<BifurcationType> listOfBifurcations;
- TrackFromThisIndex(listOfBifurcations, skeleton, index, label);
- DD("end track");
- DD(listOfBifurcations.size());
- writeImage<MaskImageType>(skeleton, "skeleton2.mhd");
-
- for(unsigned int i=0; i<listOfBifurcations.size(); i++) {
- typename MaskImageType::PointType p;
- skeleton->TransformIndexToPhysicalPoint(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<BifurcationType> listOfBifurcations;
+ TrackFromThisIndex(listOfBifurcations, skeleton, index, label);
+ DD("end track");
+ DD(listOfBifurcations.size());
+ writeImage<MaskImageType>(skeleton, "skeleton2.mhd");
+
+ for(unsigned int i=0; i<listOfBifurcations.size(); i++) {
+ typename MaskImageType::PointType p;
+ skeleton->TransformIndexToPhysicalPoint(listOfBifurcations[i].index, p);
+ DD(p);
+ }
+ }
}
}
//--------------------------------------------------------------------
template <class ImageType, class MaskImageType>
void
clitk::ExtractLungFilter<ImageType, MaskImageType>::
-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.
TrackFromThisIndex(std::vector<BifurcationType> & listOfBifurcations,
MaskImagePointer skeleton,
MaskImageIndexType index,
- MaskImagePixelType label) {
+ MaskImagePixelType label)
+{
DD("TrackFromThisIndex");
DD(index);
DD((int)label);
}
//--------------------------------------------------------------------
+
+//--------------------------------------------------------------------
+template <class ImageType, class MaskImageType>
+bool
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+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<ImageType> 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 <class ImageType, class MaskImageType>
+void
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+TracheaRegionGrowing()
+{
+ // Explosion controlled region growing
+ typedef clitk::ExplosionControlledThresholdConnectedImageFilter<ImageType, InternalImageType> 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; i<m_Seeds.size();i++) {
+ f->AddSeed(m_Seeds[i]);
+ // DD(m_Seeds[i]);
+ }
+ f->Update();
+
+ // take first (main) connected component
+ writeImage<InternalImageType>(f->GetOutput(), "t1.mhd");
+ trachea_tmp = Labelize<InternalImageType>(f->GetOutput(),
+ GetBackgroundValue(),
+ true,
+ GetMinimalComponentSize());
+ writeImage<InternalImageType>(trachea_tmp, "t2.mhd");
+ trachea_tmp = KeepLabels<InternalImageType>(trachea_tmp,
+ GetBackgroundValue(),
+ GetForegroundValue(),
+ 1, 1, false);
+ writeImage<InternalImageType>(trachea_tmp, "t3.mhd");
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType, class MaskImageType>
+double
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+ComputeTracheaVolume()
+{
+ typedef itk::ImageRegionConstIterator<InternalImageType> 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 <class ImageType, class MaskImageType>
+void
+clitk::ExtractLungFilter<ImageType, MaskImageType>::
+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<InternalImageType>(trachea_tmp);
+ }
+ else { // Trachea not found
+ this->SetWarning("* WARNING * No seed found for trachea.");
+ StopCurrentStep();
+ }
+}
+//--------------------------------------------------------------------
+
#endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX
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();
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();
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
- Patient label image
- Lungs label image
- Bones label image
+ - Trachea label image
*/
//--------------------------------------------------------------------
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);
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);
ImagePixelType m_BackgroundValuePatient;
ImagePixelType m_BackgroundValueLung;
ImagePixelType m_BackgroundValueBones;
+ ImagePixelType m_BackgroundValueTrachea;
ImagePixelType m_ForegroundValueLeftLung;
ImagePixelType m_ForegroundValueRightLung;
#include "clitkExtractMediastinumFilter.h"
#include "clitkAddRelativePositionConstraintToLabelImageFilter.h"
#include "clitkSegmentationUtils.h"
+#include "clitkExtractAirwayTreeInfoFilter.h"
// itk
#include <deque>
#include "RelativePositionPropImageFilter.h"
//--------------------------------------------------------------------
-template <class TImageType>
-clitk::ExtractMediastinumFilter<TImageType>::
+template <class ImageType>
+clitk::ExtractMediastinumFilter<ImageType>::
ExtractMediastinumFilter():
clitk::FilterBase(),
- itk::ImageToImageFilter<TImageType, TImageType>()
+ itk::ImageToImageFilter<ImageType, ImageType>()
{
- this->SetNumberOfRequiredInputs(3);
+ this->SetNumberOfRequiredInputs(4);
SetBackgroundValuePatient(0);
SetBackgroundValueLung(0);
//--------------------------------------------------------------------
-template <class TImageType>
+template <class ImageType>
void
-clitk::ExtractMediastinumFilter<TImageType>::
-SetInputPatientLabelImage(const TImageType * image, ImagePixelType bg) {
- this->SetNthInput(0, const_cast<TImageType *>(image));
+clitk::ExtractMediastinumFilter<ImageType>::
+SetInputPatientLabelImage(const ImageType * image, ImagePixelType bg)
+{
+ this->SetNthInput(0, const_cast<ImageType *>(image));
m_BackgroundValuePatient = bg;
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
-template <class TImageType>
+template <class ImageType>
void
-clitk::ExtractMediastinumFilter<TImageType>::
-SetInputLungLabelImage(const TImageType * image, ImagePixelType bg,
- ImagePixelType fgLeft, ImagePixelType fgRight) {
- this->SetNthInput(1, const_cast<TImageType *>(image));
+clitk::ExtractMediastinumFilter<ImageType>::
+SetInputLungLabelImage(const ImageType * image, ImagePixelType bg,
+ ImagePixelType fgLeft, ImagePixelType fgRight)
+{
+ this->SetNthInput(1, const_cast<ImageType *>(image));
m_BackgroundValueLung = bg;
m_ForegroundValueLeftLung = fgLeft;
m_ForegroundValueRightLung = fgRight;
//--------------------------------------------------------------------
-template <class TImageType>
+template <class ImageType>
void
-clitk::ExtractMediastinumFilter<TImageType>::
-SetInputBonesLabelImage(const TImageType * image, ImagePixelType bg) {
- this->SetNthInput(2, const_cast<TImageType *>(image));
+clitk::ExtractMediastinumFilter<ImageType>::
+SetInputBonesLabelImage(const ImageType * image, ImagePixelType bg)
+{
+ this->SetNthInput(2, const_cast<ImageType *>(image));
m_BackgroundValueBones = bg;
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
-template <class TImageType>
+template <class ImageType>
+void
+clitk::ExtractMediastinumFilter<ImageType>::
+SetInputTracheaLabelImage(const ImageType * image, ImagePixelType bg)
+{
+ this->SetNthInput(3, const_cast<ImageType *>(image));
+ m_BackgroundValueTrachea = bg;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
template<class ArgsInfoType>
void
-clitk::ExtractMediastinumFilter<TImageType>::
+clitk::ExtractMediastinumFilter<ImageType>::
SetArgsInfo(ArgsInfoType mArgsInfo)
{
SetVerboseOption_GGO(mArgsInfo);
SetBackgroundValuePatient_GGO(mArgsInfo);
SetBackgroundValueLung_GGO(mArgsInfo);
- SetBackgroundValueBones_GGO(mArgsInfo);
+ SetBackgroundValueTrachea_GGO(mArgsInfo);
SetForegroundValueLeftLung_GGO(mArgsInfo);
SetForegroundValueRightLung_GGO(mArgsInfo);
//--------------------------------------------------------------------
-template <class TImageType>
+template <class ImageType>
void
-clitk::ExtractMediastinumFilter<TImageType>::
+clitk::ExtractMediastinumFilter<ImageType>::
GenerateOutputInformation() {
- ImagePointer input = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(0));
+ ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
ImagePointer outputImage = this->GetOutput(0);
outputImage->SetRegions(outputImage->GetLargestPossibleRegion());
}
//--------------------------------------------------------------------
-template <class TImageType>
+template <class ImageType>
void
-clitk::ExtractMediastinumFilter<TImageType>::
-GenerateInputRequestedRegion() {
+clitk::ExtractMediastinumFilter<ImageType>::
+GenerateInputRequestedRegion()
+{
// Call default
Superclass::GenerateInputRequestedRegion();
// Get input pointers
- ImagePointer patient = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(0));
- ImagePointer lung = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(1));
- ImagePointer bones = dynamic_cast<TImageType*>(itk::ProcessObject::GetInput(2));
+ ImagePointer patient = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
+ ImagePointer lung = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
+ ImagePointer bones = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(2));
+ ImagePointer trachea = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(3));
patient->SetRequestedRegion(patient->GetLargestPossibleRegion());
lung->SetRequestedRegion(lung->GetLargestPossibleRegion());
bones->SetRequestedRegion(bones->GetLargestPossibleRegion());
+ trachea->SetRequestedRegion(trachea->GetLargestPossibleRegion());
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
-template <class TImageType>
+template <class ImageType>
void
-clitk::ExtractMediastinumFilter<TImageType>::
-GenerateData() {
+clitk::ExtractMediastinumFilter<ImageType>::
+GenerateData()
+{
// Get input pointers
- ImageConstPointer patient = dynamic_cast<const TImageType*>(itk::ProcessObject::GetInput(0));
- ImageConstPointer lung = dynamic_cast<const TImageType*>(itk::ProcessObject::GetInput(1));
- ImageConstPointer bones = dynamic_cast<const TImageType*>(itk::ProcessObject::GetInput(2));
+ ImageConstPointer patient = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
+ ImageConstPointer lung = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(1));
+ ImageConstPointer bones = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(2));
+ ImageConstPointer trachea = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(3));
// Get output pointer
ImagePointer output;
boolFilter->Update();
output = boolFilter->GetOutput();
+ // Auto crop to gain support area
output = clitk::AutoCrop<ImageType>(output, GetBackgroundValue());
- ////autoCropFilter->GetOutput(); typedef clitk::AutoCropFilter<ImageType> CropFilterType;
- //typename CropFilterType::Pointer cropFilter = CropFilterType::New();
- //cropFilter->SetInput(output);
- //cropFilter->Update();
- //output = cropFilter->GetOutput();
-
- this->template StopCurrentStep<TImageType>(output);
+ this->template StopCurrentStep<ImageType>(output);
// Step 2: LR limits from lung (need separate lung ?)
StartNewStep("Left/Right limits with lungs");
- typedef clitk::AddRelativePositionConstraintToLabelImageFilter<TImageType> 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<ImageType, ImageType>(lung, lung, 2, 0);
+ ImagePointer left_lung = clitk::SetBackground<ImageType, ImageType>(lung, lung, 1, 0);
+ writeImage<ImageType>(right_lung, "right.mhd");
+ writeImage<ImageType>(left_lung, "left.mhd");
+ */
+
+ typedef clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType> 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<TImageType>(output);
+ DD(output->GetLargestPossibleRegion());
+ this->template StopCurrentStep<ImageType>(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<ImageType> 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<ImageType, ImageType> CropFilterType;
+ // typedef itk::ExtractImageFilter<ImageType, ImageType> 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<ImageType>(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<ImageType>(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<ImageType>(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<TImageType>(output);
-
+ this->template StopCurrentStep<ImageType>(output);
// Get CCL
- output = clitk::Labelize<TImageType>(output, GetBackgroundValue(), true, 100);
- // output = RemoveLabels<TImageType>(output, BG, param->GetLabelsToRemove());
- output = clitk::KeepLabels<TImageType>(output, GetBackgroundValue(),
- GetForegroundValue(), 1, 1, 0);
+ output = clitk::Labelize<ImageType>(output, GetBackgroundValue(), true, 100);
+ // output = RemoveLabels<ImageType>(output, BG, param->GetLabelsToRemove());
+ output = clitk::KeepLabels<ImageType>(output, GetBackgroundValue(),
+ GetForegroundValue(), 1, 1, 0);
output = clitk::AutoCrop<ImageType>(output, GetBackgroundValue());
// cropFilter = CropFilterType::New();
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);
}
//--------------------------------------------------------------------
{
// Reading input
typename ImageType::Pointer patient = this->template GetInput<ImageType>(0);
- typename ImageType::Pointer lung = this->template GetInput<ImageType>(1);
- typename ImageType::Pointer bones = this->template GetInput<ImageType>(2);
+ typename ImageType::Pointer lung = this->template GetInput<ImageType>(1);
+ typename ImageType::Pointer bones = this->template GetInput<ImageType>(2);
+ typename ImageType::Pointer trachea = this->template GetInput<ImageType>(3);
// Create filter
typedef clitk::ExtractMediastinumFilter<ImageType> FilterType;
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 !
git://git.creatis.insa-lyon.fr / clitk.git / commitdiff