X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=segmentation%2FclitkExtractLungFilter.txx;h=92b2d4f75cf2e8daa0c6633e64c23c7812886ffd;hb=090511076dfb931319d0ec16fa5a08f2467362b6;hp=f291739d09cab2f346ca6f2aae46ab72493c0fca;hpb=c8d53d40ac8f5f83843193c51ac1787a8f38ee90;p=clitk.git diff --git a/segmentation/clitkExtractLungFilter.txx b/segmentation/clitkExtractLungFilter.txx index f291739..92b2d4f 100644 --- a/segmentation/clitkExtractLungFilter.txx +++ b/segmentation/clitkExtractLungFilter.txx @@ -24,25 +24,33 @@ #include "clitkSetBackgroundImageFilter.h" #include "clitkSegmentationUtils.h" #include "clitkAutoCropFilter.h" +#include "clitkExtractSliceFilter.h" // itk #include "itkBinaryThresholdImageFilter.h" #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() + clitk::FilterWithAnatomicalFeatureDatabaseManagement(), + 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); @@ -58,6 +66,7 @@ ExtractLungFilter(): SetUpperThresholdForTrachea(-900); SetMultiplierForTrachea(5); SetThresholdStepSizeForTrachea(64); + SetNumberOfSlicesToSkipBeforeSearchingSeed(0); // Step 3 default values SetNumberOfHistogramBins(500); @@ -74,37 +83,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); @@ -113,10 +125,10 @@ AddSeed(InternalIndexType s) //-------------------------------------------------------------------- -template +template template void -clitk::ExtractLungFilter:: +clitk::ExtractLungFilter:: SetArgsInfo(ArgsInfoType mArgsInfo) { SetVerboseOption_GGO(mArgsInfo); @@ -126,6 +138,7 @@ SetArgsInfo(ArgsInfoType mArgsInfo) SetUpperThreshold_GGO(mArgsInfo); SetLowerThreshold_GGO(mArgsInfo); + SetNumberOfSlicesToSkipBeforeSearchingSeed_GGO(mArgsInfo); SetLabelizeParameters1_GGO(mArgsInfo); if (!mArgsInfo.remove1_given) { GetLabelizeParameters1()->AddLabelToRemove(2); @@ -137,47 +150,55 @@ SetArgsInfo(ArgsInfoType mArgsInfo) SetThresholdStepSizeForTrachea_GGO(mArgsInfo); AddSeed_GGO(mArgsInfo); + SetMinimalComponentSize_GGO(mArgsInfo); + SetNumberOfHistogramBins_GGO(mArgsInfo); SetLabelizeParameters2_GGO(mArgsInfo); SetRadiusForTrachea_GGO(mArgsInfo); SetLabelizeParameters3_GGO(mArgsInfo); + + SetAFDBFilename_GGO(mArgsInfo); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -template +template void -clitk::ExtractLungFilter:: +clitk::ExtractLungFilter:: GenerateOutputInformation() { - input = dynamic_cast(itk::ProcessObject::GetInput(0)); Superclass::GenerateOutputInformation(); -// MaskImagePointer output = this->GetOutput(0); + //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)) { - this->SetLastError("* ERROR * the images (input and patient mask) must have the same size & spacing"); - return; + if (!HaveSameSizeAndSpacing(input, patient)) { + clitkExceptionMacro("the 'input' and 'patient' masks must have the same size & spacing."); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- StartNewStep("Set background to initial image"); - working_input = SetBackground + working_input = SetBackground (input, patient, GetPatientMaskBackgroundValue(), -1000); - StopCurrentStep(working_input); + StopCurrentStep(working_input); //-------------------------------------------------------------------- //-------------------------------------------------------------------- StartNewStep("Remove Air"); + // Check threshold + if (m_UseLowerThreshold) { + if (m_LowerThreshold > m_UpperThreshold) { + clitkExceptionMacro("lower threshold cannot be greater than upper threshold."); + } + } // 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); @@ -189,8 +210,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(), @@ -200,79 +223,22 @@ 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()); - if (m_Seeds.size() == 0) { // try to find seed - // Search seed (parameters = UpperThresholdForTrachea) - static const unsigned int Dim = InputImageType::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); - //DD(GetUpperThresholdForTrachea()); - //DD(sliceRegion); - 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; - 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."); - // Set output - StopCurrentStep(); - } + StartNewStep("Search for the trachea"); + SearchForTrachea(); //-------------------------------------------------------------------- //-------------------------------------------------------------------- StartNewStep("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()); @@ -307,8 +273,8 @@ GenerateOutputInformation() working_image = SetBackground (working_image, trachea_tmp, 1, -1); - // Dilate the trachea - static const unsigned int Dim = InputImageType::ImageDimension; + // Dilate the trachea + static const unsigned int Dim = ImageType::ImageDimension; typedef itk::BinaryBallStructuringElement KernelType; KernelType structuringElement; structuringElement.SetRadius(GetRadiusForTrachea()); @@ -343,7 +309,6 @@ GenerateOutputInformation() // Set output StopCurrentStep(working_image); } - //-------------------------------------------------------------------- //-------------------------------------------------------------------- @@ -388,10 +353,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(); @@ -409,7 +374,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); @@ -418,26 +383,376 @@ GenerateOutputInformation() thresholdFilter->Update(); working_image = thresholdFilter->GetOutput(); StopCurrentStep (working_image); + + // Final Cast + StartNewStep("Cast the lung mask"); + 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) + if (m_Seeds.size() != 0) { // if ==0 ->no trachea found + + if (GetFindBronchialBifurcations()) { + StartNewStep("Find bronchial bifurcations"); + // Step 1 : extract skeleton + typedef itk::BinaryThinningImageFilter3D ThinningFilterType; + typename ThinningFilterType::Pointer thinningFilter = ThinningFilterType::New(); + thinningFilter->SetInput(trachea); + thinningFilter->Update(); + typename MaskImageType::Pointer skeleton = thinningFilter->GetOutput(); + + // 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()) { + clitkExceptionMacro("first point in the trachea skeleton not found."); + } + typename MaskImageType::IndexType index = it.GetIndex(); + + // Step 2.2 : initialize neighborhooditerator + typedef itk::NeighborhoodIterator NeighborhoodIteratorType; + typename NeighborhoodIteratorType::SizeType radius; + radius.Fill(1); + NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion()); + + // Find first label number (must be different from BG and FG) + typename MaskImageType::PixelType label = GetForegroundValue()+1; + while ((label == GetBackgroundValue()) || (label == GetForegroundValue())) { label++; } + + // Track from the first point + std::vector listOfBifurcations; + m_SkeletonTree.set_head(index); + TrackFromThisIndex(listOfBifurcations, skeleton, index, label, m_SkeletonTree.begin()); + DD("end track"); + DD(listOfBifurcations.size()); + DD(m_SkeletonTree.size()); + + for(unsigned int i=0; iTransformIndexToPhysicalPoint(listOfBifurcations[i].index, + listOfBifurcations[i].point); + } + + // Search for the first slice that separate the bronchus (carena) + typedef clitk::ExtractSliceFilter ExtractSliceFilterType; + typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New(); + extractSliceFilter->SetInput(trachea); + extractSliceFilter->SetDirection(2); + extractSliceFilter->Update(); + typedef typename ExtractSliceFilterType::SliceType SliceType; + std::vector mInputSlices; + extractSliceFilter->GetOutputSlices(mInputSlices); + + bool stop = false; + DD(listOfBifurcations[0].index); + DD(listOfBifurcations[1].index); + int slice_index = listOfBifurcations[0].index[2]; // first slice from carena in skeleton + int i=0; + TreeIterator firstIter = m_SkeletonTree.child(listOfBifurcations[1].treeIter, 0); + TreeIterator secondIter = m_SkeletonTree.child(listOfBifurcations[1].treeIter, 1); + typename SliceType::IndexType in1; + typename SliceType::IndexType in2; + while (!stop) { + DD(slice_index); + + // Labelize the current slice + typename SliceType::Pointer temp = Labelize(mInputSlices[slice_index], + GetBackgroundValue(), + true, + GetMinimalComponentSize()); + // Check the value of the two skeleton points; + in1[0] = (*firstIter)[0]; + in1[1] = (*firstIter)[1]; + typename SliceType::PixelType v1 = temp->GetPixel(in1); + DD(in1); + DD((int)v1); + in2[0] = (*secondIter)[0]; + in2[1] = (*secondIter)[1]; + typename SliceType::PixelType v2 = temp->GetPixel(in2); + DD(in2); + DD((int)v2); + + // TODO IF NOT FOUND ???? + + if (v1 != v2) { + stop = true; + } + else { + i++; + --slice_index; + ++firstIter; + ++secondIter; + } + } + MaskImageIndexType carena_index; + carena_index[0] = lrint(in2[0] + in1[0])/2.0; + carena_index[1] = lrint(in2[1] + in1[1])/2.0; + carena_index[2] = slice_index; + MaskImagePointType carena_position; + DD(carena_index); + skeleton->TransformIndexToPhysicalPoint(carena_index, + carena_position); + DD(carena_position); + + // Set and save Carina position + StartNewStep("Save carina position"); + // Try to load the DB + try { + LoadAFDB(); + } catch (clitk::ExceptionObject e) { + // Do nothing if not found, it will be used anyway to write the result + } + GetAFDB()->SetPoint3D("Carena", carena_position); + } + } } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -template +template void -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; +clitk::ExtractLungFilter:: +GenerateData() +{ + // 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, + TreeIterator currentNode) +{ + // 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); + nit.SetCenterPixel(label); + listOfTrackedPoint.clear(); + for(unsigned int i=0; iAdd(listOfTrackedPoint[0], index); // the parent is 'index' + // m_SkeletonTree->Add(listOfTrackedPoint[1], index); // the parent is 'index' + BifurcationType bif(index, label, label+1, label+2); + bif.treeIter = currentNode; + listOfBifurcations.push_back(bif); + TreeIterator firstNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[0]); + TreeIterator secondNode = m_SkeletonTree.append_child(currentNode, listOfTrackedPoint[1]); + TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[0], label+1, firstNode); + TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[1], label+2, secondNode); + } + else { + if (listOfTrackedPoint.size() > 2) { + clitkExceptionMacro("error while tracking trachea bifurcation. Too much bifurcation points ... ?"); + } + // Else this it the end of the tracking + } + stop = true; + } + } +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +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 + trachea_tmp = Labelize(f->GetOutput(), + GetBackgroundValue(), + true, + GetMinimalComponentSize()); + trachea_tmp = KeepLabels(trachea_tmp, + GetBackgroundValue(), + GetForegroundValue(), + 1, 1, false); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +double +clitk::ExtractLungFilter:: +ComputeTracheaVolume() +{ + typedef itk::ImageRegionConstIterator IteratorType; + IteratorType iter(trachea_tmp, trachea_tmp->GetLargestPossibleRegion()); + iter.GoToBegin(); + double volume = 0.0; + while (!iter.IsAtEnd()) { + if (iter.Get() == this->GetForegroundValue()) volume++; + ++iter; + } + + double voxelsize = trachea_tmp->GetSpacing()[0]*trachea_tmp->GetSpacing()[1]*trachea_tmp->GetSpacing()[2]; + return volume*voxelsize; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLungFilter:: +SearchForTrachea() +{ + // Search for seed among n slices, skip some slices before starting + // if not found -> skip more and restart + + // when seed found : perform region growing + // compute trachea volume + // if volume not plausible -> skip more slices and restart + + bool stop = false; + double volume = 0.0; + int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed(); + while (!stop) { + stop = SearchForTracheaSeed(skip); + if (stop) { + TracheaRegionGrowing(); + volume = ComputeTracheaVolume()/1000; // assume mm3, so divide by 1000 to get cc + if ((volume > 10) && (volume < 55 )) { // it is ok + // Typical volume 22.59 cm 3 (± 7.69 cm 3 ) [ Leader 2004 ] + if (GetVerboseStep()) { + std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl; + } + stop = true; + } + else { + if (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