1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://www.centreleonberard.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
9 This software is distributed WITHOUT ANY WARRANTY; without even
10 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 PURPOSE. See the copyright notices for more information.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ===========================================================================**/
19 #ifndef CLITKEXTRACTLUNGSFILTER_TXX
20 #define CLITKEXTRACTLUNGSFILTER_TXX
23 #include "clitkImageCommon.h"
24 #include "clitkSetBackgroundImageFilter.h"
25 #include "clitkSegmentationUtils.h"
26 #include "clitkAutoCropFilter.h"
27 #include "clitkCropLikeImageFilter.h"
28 #include "clitkFillMaskFilter.h"
29 #include "clitkMemoryUsage.h"
32 #include "itkBinaryThresholdImageFilter.h"
33 #include "itkConnectedComponentImageFilter.h"
34 #include "itkRelabelComponentImageFilter.h"
35 #include "itkOtsuThresholdImageFilter.h"
36 #include "itkBinaryThinningImageFilter3D.h"
37 #include "itkImageIteratorWithIndex.h"
38 #include "itkBinaryMorphologicalOpeningImageFilter.h"
39 #include "itkBinaryMorphologicalClosingImageFilter.h"
41 //--------------------------------------------------------------------
42 template <class ImageType>
43 clitk::ExtractLungFilter<ImageType>::
46 clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
47 itk::ImageToImageFilter<ImageType, MaskImageType>()
50 m_MaxSeedNumber = 500;
52 // Default global options
53 this->SetNumberOfRequiredInputs(1);
54 SetPatientMaskBackgroundValue(0);
55 SetBackgroundValue(0); // Must be zero
56 SetForegroundValue(1);
57 SetMinimalComponentSize(100);
58 VerboseRegionGrowingFlagOff();
60 // Step 1 default values
61 SetUpperThreshold(-300);
62 SetLowerThreshold(-1000);
63 UseLowerThresholdOff();
64 LabelParamType * p1 = new LabelParamType;
67 p1->AddLabelToRemove(2);
68 SetLabelizeParameters1(p1);
70 // Step 2 default values
71 SetUpperThresholdForTrachea(-900);
72 SetMultiplierForTrachea(5);
73 SetThresholdStepSizeForTrachea(64);
74 SetNumberOfSlicesToSkipBeforeSearchingSeed(0);
75 TracheaVolumeMustBeCheckedFlagOn();
77 // Step 3 default values
78 SetNumberOfHistogramBins(500);
79 LabelParamType * p2 = new LabelParamType;
82 // p->AddLabelToRemove(2); // No label to remove by default
83 SetLabelizeParameters2(p2);
85 // Step 4 default values
86 SetRadiusForTrachea(1);
87 LabelParamType * p3 = new LabelParamType;
91 SetLabelizeParameters3(p3);
95 SetOpenCloseRadius(1);
100 //--------------------------------------------------------------------
103 //--------------------------------------------------------------------
104 template <class ImageType>
106 clitk::ExtractLungFilter<ImageType>::
107 SetInput(const ImageType * image)
109 this->SetNthInput(0, const_cast<ImageType *>(image));
111 //--------------------------------------------------------------------
114 //--------------------------------------------------------------------
115 template <class ImageType>
117 clitk::ExtractLungFilter<ImageType>::
118 AddSeed(InternalIndexType s)
120 m_Seeds.push_back(s);
122 //--------------------------------------------------------------------
125 //--------------------------------------------------------------------
126 template <class ImageType>
128 clitk::ExtractLungFilter<ImageType>::
129 GenerateOutputInformation()
131 clitk::PrintMemory(GetVerboseMemoryFlag(), "Initial memory"); // OK
132 Superclass::GenerateOutputInformation();
137 // Get input pointers
138 input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
139 patient = GetAFDB()->template GetImage <MaskImageType>("Patient");
140 PrintMemory(GetVerboseMemoryFlag(), "After reading patient"); // OK
142 //--------------------------------------------------------------------
143 //--------------------------------------------------------------------
144 // Crop input like patient image (must have the same spacing)
145 // It copy the input if the same are the same
146 StartNewStep("Copy and crop input image to 'patient' extends");
147 typedef clitk::CropLikeImageFilter<ImageType> CropImageFilter;
148 typename CropImageFilter::Pointer cropFilter = CropImageFilter::New();
149 // cropFilter->ReleaseDataFlagOn(); // NO=seg fault !!
150 cropFilter->SetInput(input);
151 cropFilter->SetCropLikeImage(patient);
152 cropFilter->Update();
153 working_input = cropFilter->GetOutput();
154 // cropFilter->Delete(); // NO !!!! if yes, sg fault, Cropfilter is buggy !??
155 StopCurrentStep<ImageType>(working_input);
156 PrintMemory(GetVerboseMemoryFlag(), "After crop"); // OK, slightly more than a copy of input
158 //--------------------------------------------------------------------
159 //--------------------------------------------------------------------
160 StartNewStep("Set background to initial image");
161 working_input = SetBackground<ImageType, MaskImageType>
162 (working_input, patient, GetPatientMaskBackgroundValue(), -1000, true);
163 StopCurrentStep<ImageType>(working_input);
164 PrintMemory(GetVerboseMemoryFlag(), "After set bg"); // OK, additional mem = 0
167 // We do not need patient mask anymore, release memory
168 GetAFDB()->template ReleaseImage<MaskImageType>("Patient");
169 DD(patient->GetReferenceCount());
170 PrintMemory(GetVerboseMemoryFlag(), "After delete patient"); // OK, additional mem = 0
172 PrintMemory(GetVerboseMemoryFlag(), "After delete patient"); // OK, additional mem = 0
175 //--------------------------------------------------------------------
176 //--------------------------------------------------------------------
177 StartNewStep("Remove Air");
179 if (m_UseLowerThreshold) {
180 if (m_LowerThreshold > m_UpperThreshold) {
181 clitkExceptionMacro("lower threshold cannot be greater than upper threshold.");
184 // Threshold to get air
185 typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> InputBinarizeFilterType;
186 typename InputBinarizeFilterType::Pointer binarizeFilter = InputBinarizeFilterType::New();
187 binarizeFilter->SetInput(working_input);
188 if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(m_LowerThreshold);
189 // binarizeFilter->CanRunInPlace() is false
190 binarizeFilter->SetUpperThreshold(m_UpperThreshold);
191 binarizeFilter->SetInsideValue(this->GetForegroundValue());
192 binarizeFilter->SetOutsideValue(this->GetBackgroundValue());
193 binarizeFilter->Update();
194 working_mask = binarizeFilter->GetOutput();
195 PrintMemory(GetVerboseMemoryFlag(), "After Binarizefilter"); // OK, additional mem is one mask image
197 // Labelize and keep right labels
199 Labelize<MaskImageType>
200 (working_mask, GetBackgroundValue(), true, GetMinimalComponentSize());
201 PrintMemory(GetVerboseMemoryFlag(), "After Labelize"); // BUG ? additional mem around 1 time the input ?
203 working_mask = RemoveLabels<MaskImageType>
204 (working_mask, GetBackgroundValue(), GetLabelizeParameters1()->GetLabelsToRemove());
205 PrintMemory(GetVerboseMemoryFlag(), "After RemoveLabels"); // OK additional mem = 0
207 working_mask = KeepLabels<MaskImageType>
209 GetBackgroundValue(),
210 GetForegroundValue(),
211 GetLabelizeParameters1()->GetFirstKeep(),
212 GetLabelizeParameters1()->GetLastKeep(),
213 GetLabelizeParameters1()->GetUseLastKeep());
214 PrintMemory(GetVerboseMemoryFlag(), "After KeepLabels to create the 'air'");
217 working_input = SetBackground<ImageType, MaskImageType>
218 (working_input, working_mask, this->GetForegroundValue(), this->GetBackgroundValue(), true);
219 PrintMemory(GetVerboseMemoryFlag(), "After SetBackground");
222 StopCurrentStep<ImageType>(working_input);
224 //--------------------------------------------------------------------
225 //--------------------------------------------------------------------
226 StartNewStep("Search for the trachea");
228 PrintMemory(GetVerboseMemoryFlag(), "After SearchForTrachea");
230 //--------------------------------------------------------------------
231 //--------------------------------------------------------------------
232 StartNewStep("Extract the lung with Otsu filter");
233 // Automated Otsu thresholding and relabeling
234 typedef itk::OtsuThresholdImageFilter<ImageType,MaskImageType> OtsuThresholdImageFilterType;
235 typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New();
236 otsuFilter->SetInput(working_input);
237 otsuFilter->SetNumberOfHistogramBins(GetNumberOfHistogramBins());
238 otsuFilter->SetInsideValue(this->GetForegroundValue());
239 otsuFilter->SetOutsideValue(this->GetBackgroundValue());
240 otsuFilter->Update();
241 working_mask = otsuFilter->GetOutput();
244 StopCurrentStep<MaskImageType>(working_mask);
245 PrintMemory(GetVerboseMemoryFlag(), "After Otsufilter");
247 //--------------------------------------------------------------------
248 //--------------------------------------------------------------------
249 StartNewStep("Select labels");
251 working_mask = LabelizeAndSelectLabels<MaskImageType>
253 GetBackgroundValue(),
254 GetForegroundValue(),
256 GetMinimalComponentSize(),
257 GetLabelizeParameters2());
260 StopCurrentStep<MaskImageType>(working_mask);
261 PrintMemory(GetVerboseMemoryFlag(), "After LabelizeAndSelectLabels");
263 //--------------------------------------------------------------------
264 //--------------------------------------------------------------------
265 if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
266 StartNewStep("Remove the trachea");
268 working_mask = SetBackground<MaskImageType, MaskImageType>
269 (working_mask, trachea, 1, -1, true);
270 PrintMemory(GetVerboseMemoryFlag(), "After SetBackground");
272 // Dilate the trachea
273 static const unsigned int Dim = ImageType::ImageDimension;
274 typedef itk::BinaryBallStructuringElement<InternalPixelType, Dim> KernelType;
275 KernelType structuringElement;
276 structuringElement.SetRadius(GetRadiusForTrachea());
277 structuringElement.CreateStructuringElement();
278 typedef clitk::ConditionalBinaryDilateImageFilter<MaskImageType, MaskImageType, KernelType> ConditionalBinaryDilateImageFilterType;
279 typename ConditionalBinaryDilateImageFilterType::Pointer dilateFilter = ConditionalBinaryDilateImageFilterType::New();
280 dilateFilter->SetBoundaryToForeground(false);
281 dilateFilter->SetKernel(structuringElement);
282 dilateFilter->SetBackgroundValue (1);
283 dilateFilter->SetForegroundValue (-1);
284 dilateFilter->SetInput (working_mask);
285 dilateFilter->Update();
286 working_mask = dilateFilter->GetOutput();
287 PrintMemory(GetVerboseMemoryFlag(), "After dilate");
289 // Set trachea with dilatation
290 trachea = SetBackground<MaskImageType, MaskImageType>
291 (trachea, working_mask, -1, this->GetForegroundValue(), true);
293 // Remove the trachea
294 working_mask = SetBackground<MaskImageType, MaskImageType>
295 (working_mask, working_mask, -1, this->GetBackgroundValue(), true);
298 working_mask = LabelizeAndSelectLabels<MaskImageType>
300 GetBackgroundValue(),
301 GetForegroundValue(),
303 GetMinimalComponentSize(),
304 GetLabelizeParameters3());
307 StopCurrentStep<MaskImageType>(working_mask);
310 //--------------------------------------------------------------------
311 //--------------------------------------------------------------------
312 PrintMemory(GetVerboseMemoryFlag(), "before autocropfilter");
313 if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
314 trachea = clitk::AutoCrop<MaskImageType>(trachea, GetBackgroundValue());
315 StopCurrentStep<MaskImageType>(trachea);
316 PrintMemory(GetVerboseMemoryFlag(), "after delete trachea");
318 PrintMemory(GetVerboseMemoryFlag(), "after delete trachea");
320 //--------------------------------------------------------------------
321 //--------------------------------------------------------------------
322 StartNewStep("Cropping lung");
323 PrintMemory(GetVerboseMemoryFlag(), "Before Autocropfilter");
324 working_mask = clitk::AutoCrop<MaskImageType>(working_mask, GetBackgroundValue());
325 StopCurrentStep<MaskImageType>(working_mask);
327 //--------------------------------------------------------------------
328 //--------------------------------------------------------------------
330 if (GetOpenCloseFlag()) {
331 StartNewStep("Open/Close");
332 PrintMemory(GetVerboseMemoryFlag(), "Before OpenClose");
334 // Structuring element
335 typedef itk::BinaryBallStructuringElement<InternalPixelType, ImageDimension> KernelType;
336 KernelType structuringElement;
337 structuringElement.SetRadius(GetOpenCloseRadius());
338 structuringElement.CreateStructuringElement();
341 typedef itk::BinaryMorphologicalOpeningImageFilter<MaskImageType, InternalImageType, KernelType> OpenFilterType;
342 typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
343 openFilter->SetInput(working_mask);
344 openFilter->SetBackgroundValue(GetBackgroundValue());
345 openFilter->SetForegroundValue(GetForegroundValue());
346 openFilter->SetKernel(structuringElement);
349 typedef itk::BinaryMorphologicalClosingImageFilter<MaskImageType, MaskImageType, KernelType> CloseFilterType;
350 typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
351 closeFilter->SetInput(openFilter->GetOutput());
352 closeFilter->SetSafeBorder(true);
353 closeFilter->SetForegroundValue(GetForegroundValue());
354 closeFilter->SetKernel(structuringElement);
355 closeFilter->Update();
356 working_mask = closeFilter->GetOutput();
359 //--------------------------------------------------------------------
360 //--------------------------------------------------------------------
362 if (GetFillHolesFlag()) {
363 StartNewStep("Fill Holes");
364 PrintMemory(GetVerboseMemoryFlag(), "Before Fill Holes");
365 typedef clitk::FillMaskFilter<MaskImageType> FillMaskFilterType;
366 typename FillMaskFilterType::Pointer fillMaskFilter = FillMaskFilterType::New();
367 fillMaskFilter->SetInput(working_mask);
368 fillMaskFilter->AddDirection(2);
369 //fillMaskFilter->AddDirection(1);
370 fillMaskFilter->Update();
371 working_mask = fillMaskFilter->GetOutput();
372 StopCurrentStep<MaskImageType>(working_mask);
375 //--------------------------------------------------------------------
376 //--------------------------------------------------------------------
377 StartNewStep("Separate Left/Right lungs");
378 PrintMemory(GetVerboseMemoryFlag(), "Before Separate");
380 working_mask = Labelize<MaskImageType>(working_mask,
381 GetBackgroundValue(),
383 GetMinimalComponentSize());
385 PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
388 typedef itk::StatisticsImageFilter<MaskImageType> StatisticsImageFilterType;
389 typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
390 statisticsImageFilter->SetInput(working_mask);
391 statisticsImageFilter->Update();
392 unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
393 working_mask = statisticsImageFilter->GetOutput();
395 PrintMemory(GetVerboseMemoryFlag(), "After count label");
397 // Decompose the first label
398 static const unsigned int Dim = ImageType::ImageDimension;
399 if (initialNumberOfLabels<2) {
400 // Structuring element radius
401 typename ImageType::SizeType radius;
402 for (unsigned int i=0;i<Dim;i++) radius[i]=1;
403 typedef clitk::DecomposeAndReconstructImageFilter<MaskImageType,MaskImageType> DecomposeAndReconstructFilterType;
404 typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
405 decomposeAndReconstructFilter->SetInput(working_mask);
406 decomposeAndReconstructFilter->SetVerbose(false);
407 decomposeAndReconstructFilter->SetRadius(radius);
408 decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
409 decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
410 decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
411 decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
412 decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
413 decomposeAndReconstructFilter->SetFullyConnected(true);
414 decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
415 decomposeAndReconstructFilter->Update();
416 working_mask = decomposeAndReconstructFilter->GetOutput();
418 PrintMemory(GetVerboseMemoryFlag(), "After decomposeAndReconstructFilter");
420 // Retain labels ('1' is largset lung, so right. '2' is left)
421 typedef itk::ThresholdImageFilter<MaskImageType> ThresholdImageFilterType;
422 typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
423 thresholdFilter->SetInput(working_mask);
424 thresholdFilter->ThresholdAbove(2);
425 thresholdFilter->SetOutsideValue(this->GetBackgroundValue());
426 thresholdFilter->Update();
427 working_mask = thresholdFilter->GetOutput();
428 StopCurrentStep<MaskImageType> (working_mask);
429 PrintMemory(GetVerboseMemoryFlag(), "After Thresholdfilter");
431 // Update output info
432 // output = working_mask;
433 //this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
435 // this->GetOutput(0)->SetRegions(working_mask->GetLargestPossibleRegion());
438 //--------------------------------------------------------------------
441 //--------------------------------------------------------------------
442 template <class TImageType>
444 clitk::ExtractLungFilter<TImageType>::
445 GenerateInputRequestedRegion() {
446 // DD("GenerateInputRequestedRegion (nothing?)");
448 //--------------------------------------------------------------------
451 //--------------------------------------------------------------------
452 template <class ImageType>
454 clitk::ExtractLungFilter<ImageType>::
458 // this->GraftOutput(output); // not SetNthOutput
459 this->GraftOutput(working_mask); // not SetNthOutput
460 // Store image filenames into AFDB
461 GetAFDB()->SetImageFilename("Lungs", this->GetOutputLungFilename());
462 GetAFDB()->SetImageFilename("Trachea", this->GetOutputTracheaFilename());
465 //--------------------------------------------------------------------
468 //--------------------------------------------------------------------
469 template <class ImageType>
471 clitk::ExtractLungFilter<ImageType>::
472 SearchForTracheaSeed(int skip)
474 if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)
475 // Restart the search until a seed is found, skipping more and more slices
478 // Search seed (parameters = UpperThresholdForTrachea)
479 static const unsigned int Dim = ImageType::ImageDimension;
480 typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
481 typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
482 typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
483 sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-skip-5;
484 sliceRegionSize[Dim-1]=5;
485 sliceRegion.SetSize(sliceRegionSize);
486 sliceRegion.SetIndex(sliceRegionIndex);
487 typedef itk::ImageRegionConstIterator<ImageType> IteratorType;
488 IteratorType it(working_input, sliceRegion);
490 while (!it.IsAtEnd()) {
491 if(it.Get() < GetUpperThresholdForTrachea() ) {
492 AddSeed(it.GetIndex());
493 // DD(it.GetIndex());
498 // if we do not found : restart
499 stop = (m_Seeds.size() != 0);
501 if (GetVerboseStepFlag()) {
502 std::cout << "\t No seed found this time. I skip some slices and restart." << std::endl;
504 if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
505 // we want to skip more than a half of the image, it is probably a bug
506 std::cerr << "Number of slices to skip to find trachea to high = " << skip << std::endl;
513 // DD(m_Seeds.size());
517 return (m_Seeds.size() != 0);
519 //--------------------------------------------------------------------
522 //--------------------------------------------------------------------
523 template <class ImageType>
525 clitk::ExtractLungFilter<ImageType>::
526 TracheaRegionGrowing()
528 // Explosion controlled region growing
529 PrintMemory(GetVerboseMemoryFlag(), "Before ExplosionControlledThresholdConnectedImageFilter");
530 typedef clitk::ExplosionControlledThresholdConnectedImageFilter<ImageType, MaskImageType> ImageFilterType;
531 typename ImageFilterType::Pointer f= ImageFilterType::New();
532 f->SetInput(working_input);
534 f->SetUpper(GetUpperThresholdForTrachea());
535 f->SetMinimumLowerThreshold(-2000);
536 // f->SetMaximumUpperThreshold(0); // MAYBE TO CHANGE ???
537 f->SetMaximumUpperThreshold(-800); // MAYBE TO CHANGE ???
538 f->SetAdaptLowerBorder(false);
539 f->SetAdaptUpperBorder(true);
540 f->SetMinimumSize(5000);
541 f->SetReplaceValue(1);
542 f->SetMultiplier(GetMultiplierForTrachea());
543 f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
544 f->SetMinimumThresholdStepSize(1);
545 f->SetVerbose(GetVerboseRegionGrowingFlag());
546 for(unsigned int i=0; i<m_Seeds.size();i++) {
547 f->AddSeed(m_Seeds[i]);
551 PrintMemory(GetVerboseMemoryFlag(), "After RG update");
553 // take first (main) connected component
554 trachea = Labelize<MaskImageType>(f->GetOutput(),
555 GetBackgroundValue(),
557 1000);//GetMinimalComponentSize());
558 PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
559 trachea = KeepLabels<MaskImageType>(trachea,
560 GetBackgroundValue(),
561 GetForegroundValue(),
563 PrintMemory(GetVerboseMemoryFlag(), "After KeepLabels");
565 //--------------------------------------------------------------------
568 //--------------------------------------------------------------------
569 template <class ImageType>
571 clitk::ExtractLungFilter<ImageType>::
572 ComputeTracheaVolume()
574 typedef itk::ImageRegionConstIterator<InternalImageType> IteratorType;
575 IteratorType iter(trachea, trachea->GetLargestPossibleRegion());
578 while (!iter.IsAtEnd()) {
579 if (iter.Get() == this->GetForegroundValue()) volume++;
583 double voxelsize = trachea->GetSpacing()[0]*trachea->GetSpacing()[1]*trachea->GetSpacing()[2];
584 return volume*voxelsize;
586 //--------------------------------------------------------------------
589 //--------------------------------------------------------------------
590 template <class ImageType>
592 clitk::ExtractLungFilter<ImageType>::
595 // Search for seed among n slices, skip some slices before starting
596 // if not found -> skip more and restart
598 // when seed found : perform region growing
599 // compute trachea volume
600 // if volume not plausible -> skip more slices and restart
604 int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed();
606 stop = SearchForTracheaSeed(skip);
608 TracheaRegionGrowing();
609 volume = ComputeTracheaVolume()/1000; // assume mm3, so divide by 1000 to get cc
610 if (GetWriteStepFlag()) {
611 writeImage<MaskImageType>(trachea, "step-trachea-"+toString(skip)+".mhd");
613 if (GetTracheaVolumeMustBeCheckedFlag()) {
614 if ((volume > 10) && (volume < 65 )) { // depend on image size ...
615 // Typical volume 22.59 cm 3 (± 7.69 cm 3 ) [ Leader 2004 ]
616 if (GetVerboseStepFlag()) {
617 std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl;
622 if (GetVerboseStepFlag()) {
623 std::cout << "\t The volume of the trachea (" << volume
624 << " cc) seems not correct. I skip some slices (" << skip << ") and restart to find seeds."
629 // empty the list of seed
641 StopCurrentStep<MaskImageType>(trachea);
643 else { // Trachea not found
644 this->SetWarning("* WARNING * No seed found for trachea.");
648 //--------------------------------------------------------------------
650 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX