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"
40 #include "itkConstantPadImageFilter.h"
42 //--------------------------------------------------------------------
43 template <class ImageType>
44 clitk::ExtractLungFilter<ImageType>::
47 clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
48 itk::ImageToImageFilter<ImageType, MaskImageType>()
51 m_MaxSeedNumber = 500;
53 // Default global options
54 this->SetNumberOfRequiredInputs(1);
55 SetPatientMaskBackgroundValue(0);
56 SetBackgroundValue(0); // Must be zero
57 SetForegroundValue(1);
58 SetMinimalComponentSize(100);
59 VerboseRegionGrowingFlagOff();
61 // Step 1 default values
62 SetUpperThreshold(-300);
63 SetLowerThreshold(-1000);
64 UseLowerThresholdOff();
65 LabelParamType * p1 = new LabelParamType;
68 p1->AddLabelToRemove(2);
69 SetLabelizeParameters1(p1);
71 // Step 2 default values
72 SetUpperThresholdForTrachea(-900);
73 SetMultiplierForTrachea(5);
74 SetThresholdStepSizeForTrachea(64);
75 SetNumberOfSlicesToSkipBeforeSearchingSeed(0);
76 TracheaVolumeMustBeCheckedFlagOn();
78 // Step 3 default values
79 SetNumberOfHistogramBins(500);
80 LabelParamType * p2 = new LabelParamType;
83 // p->AddLabelToRemove(2); // No label to remove by default
84 SetLabelizeParameters2(p2);
86 // Step 4 default values
87 SetRadiusForTrachea(1);
88 LabelParamType * p3 = new LabelParamType;
92 SetLabelizeParameters3(p3);
96 SetOpenCloseRadius(1);
102 //--------------------------------------------------------------------
105 //--------------------------------------------------------------------
106 template <class ImageType>
108 clitk::ExtractLungFilter<ImageType>::
109 SetInput(const ImageType * image)
111 this->SetNthInput(0, const_cast<ImageType *>(image));
113 //--------------------------------------------------------------------
116 //--------------------------------------------------------------------
117 template <class ImageType>
119 clitk::ExtractLungFilter<ImageType>::
120 AddSeed(InternalIndexType s)
122 m_Seeds.push_back(s);
124 //--------------------------------------------------------------------
127 //--------------------------------------------------------------------
128 template <class ImageType>
130 clitk::ExtractLungFilter<ImageType>::
131 GenerateOutputInformation()
133 clitk::PrintMemory(GetVerboseMemoryFlag(), "Initial memory"); // OK
134 Superclass::GenerateOutputInformation();
139 // Get input pointers
140 input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
141 patient = GetAFDB()->template GetImage <MaskImageType>("Patient");
142 PrintMemory(GetVerboseMemoryFlag(), "After reading patient"); // OK
144 //--------------------------------------------------------------------
145 //--------------------------------------------------------------------
146 // Crop input like patient image (must have the same spacing)
147 // It copy the input if the same are the same
148 StartNewStep("Copy and crop input image to 'patient' extends");
149 typedef clitk::CropLikeImageFilter<ImageType> CropImageFilter;
150 typename CropImageFilter::Pointer cropFilter = CropImageFilter::New();
151 // cropFilter->ReleaseDataFlagOn(); // NO=seg fault !!
152 cropFilter->SetInput(input);
153 cropFilter->SetCropLikeImage(patient);
154 cropFilter->Update();
155 working_input = cropFilter->GetOutput();
156 // cropFilter->Delete(); // NO !!!! if yes, sg fault, Cropfilter is buggy !??
157 StopCurrentStep<ImageType>(working_input);
158 PrintMemory(GetVerboseMemoryFlag(), "After crop"); // OK, slightly more than a copy of input
160 //--------------------------------------------------------------------
161 //--------------------------------------------------------------------
162 StartNewStep("Set background to initial image");
163 working_input = SetBackground<ImageType, MaskImageType>
164 (working_input, patient, GetPatientMaskBackgroundValue(), -1000, true);
166 // Pad images with air to prevent patient touching the image border
167 static const unsigned int Dim = ImageType::ImageDimension;
168 typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadFilterType;
169 typename PadFilterType::Pointer padFilter = PadFilterType::New();
170 padFilter->SetInput(working_input);
171 padFilter->SetConstant(-1000);
172 typename ImageType::SizeType bounds;
173 for (unsigned i = 0; i < Dim - 1; ++i)
176 padFilter->SetPadLowerBound(bounds);
177 padFilter->SetPadUpperBound(bounds);
179 working_input = padFilter->GetOutput();
181 StopCurrentStep<ImageType>(working_input);
182 PrintMemory(GetVerboseMemoryFlag(), "After set bg"); // OK, additional mem = 0
185 // We do not need patient mask anymore, release memory
186 GetAFDB()->template ReleaseImage<MaskImageType>("Patient");
187 DD(patient->GetReferenceCount());
188 PrintMemory(GetVerboseMemoryFlag(), "After delete patient"); // OK, additional mem = 0
190 PrintMemory(GetVerboseMemoryFlag(), "After delete patient"); // OK, additional mem = 0
193 //--------------------------------------------------------------------
194 //--------------------------------------------------------------------
195 StartNewStep("Remove Air");
197 if (m_UseLowerThreshold) {
198 if (m_LowerThreshold > m_UpperThreshold) {
199 clitkExceptionMacro("lower threshold cannot be greater than upper threshold.");
202 // Threshold to get air
203 typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> InputBinarizeFilterType;
204 typename InputBinarizeFilterType::Pointer binarizeFilter = InputBinarizeFilterType::New();
205 binarizeFilter->SetInput(working_input);
206 if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(m_LowerThreshold);
207 // binarizeFilter->CanRunInPlace() is false
208 binarizeFilter->SetUpperThreshold(m_UpperThreshold);
209 binarizeFilter->SetInsideValue(this->GetForegroundValue());
210 binarizeFilter->SetOutsideValue(this->GetBackgroundValue());
211 binarizeFilter->Update();
212 working_mask = binarizeFilter->GetOutput();
213 PrintMemory(GetVerboseMemoryFlag(), "After Binarizefilter"); // OK, additional mem is one mask image
215 // Labelize and keep right labels
217 Labelize<MaskImageType>
218 (working_mask, GetBackgroundValue(), true, GetMinimalComponentSize());
219 PrintMemory(GetVerboseMemoryFlag(), "After Labelize"); // BUG ? additional mem around 1 time the input ?
221 working_mask = RemoveLabels<MaskImageType>
222 (working_mask, GetBackgroundValue(), GetLabelizeParameters1()->GetLabelsToRemove());
223 PrintMemory(GetVerboseMemoryFlag(), "After RemoveLabels"); // OK additional mem = 0
225 working_mask = KeepLabels<MaskImageType>
227 GetBackgroundValue(),
228 GetForegroundValue(),
229 GetLabelizeParameters1()->GetFirstKeep(),
230 GetLabelizeParameters1()->GetLastKeep(),
231 GetLabelizeParameters1()->GetUseLastKeep());
232 PrintMemory(GetVerboseMemoryFlag(), "After KeepLabels to create the 'air'");
235 working_input = SetBackground<ImageType, MaskImageType>
236 (working_input, working_mask, this->GetForegroundValue(), this->GetBackgroundValue(), true);
237 PrintMemory(GetVerboseMemoryFlag(), "After SetBackground");
240 StopCurrentStep<ImageType>(working_input);
242 //--------------------------------------------------------------------
243 //--------------------------------------------------------------------
244 StartNewStep("Search for the trachea");
246 PrintMemory(GetVerboseMemoryFlag(), "After SearchForTrachea");
248 //--------------------------------------------------------------------
249 //--------------------------------------------------------------------
250 StartNewStep("Extract the lung with Otsu filter");
251 // Automated Otsu thresholding and relabeling
252 typedef itk::OtsuThresholdImageFilter<ImageType,MaskImageType> OtsuThresholdImageFilterType;
253 typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New();
254 otsuFilter->SetInput(working_input);
255 otsuFilter->SetNumberOfHistogramBins(GetNumberOfHistogramBins());
256 otsuFilter->SetInsideValue(this->GetForegroundValue());
257 otsuFilter->SetOutsideValue(this->GetBackgroundValue());
258 otsuFilter->Update();
259 working_mask = otsuFilter->GetOutput();
262 StopCurrentStep<MaskImageType>(working_mask);
263 PrintMemory(GetVerboseMemoryFlag(), "After Otsufilter");
265 //--------------------------------------------------------------------
266 //--------------------------------------------------------------------
267 StartNewStep("Select labels");
269 working_mask = LabelizeAndSelectLabels<MaskImageType>
271 GetBackgroundValue(),
272 GetForegroundValue(),
274 GetMinimalComponentSize(),
275 GetLabelizeParameters2());
278 StopCurrentStep<MaskImageType>(working_mask);
279 PrintMemory(GetVerboseMemoryFlag(), "After LabelizeAndSelectLabels");
281 //--------------------------------------------------------------------
282 //--------------------------------------------------------------------
283 if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
284 StartNewStep("Remove the trachea");
286 working_mask = SetBackground<MaskImageType, MaskImageType>
287 (working_mask, trachea, 1, -1, true);
288 PrintMemory(GetVerboseMemoryFlag(), "After SetBackground");
290 // Dilate the trachea
291 static const unsigned int Dim = ImageType::ImageDimension;
292 typedef itk::BinaryBallStructuringElement<InternalPixelType, Dim> KernelType;
293 KernelType structuringElement;
294 structuringElement.SetRadius(GetRadiusForTrachea());
295 structuringElement.CreateStructuringElement();
296 typedef clitk::ConditionalBinaryDilateImageFilter<MaskImageType, MaskImageType, KernelType> ConditionalBinaryDilateImageFilterType;
297 typename ConditionalBinaryDilateImageFilterType::Pointer dilateFilter = ConditionalBinaryDilateImageFilterType::New();
298 dilateFilter->SetBoundaryToForeground(false);
299 dilateFilter->SetKernel(structuringElement);
300 dilateFilter->SetBackgroundValue (1);
301 dilateFilter->SetForegroundValue (-1);
302 dilateFilter->SetInput (working_mask);
303 dilateFilter->Update();
304 working_mask = dilateFilter->GetOutput();
305 PrintMemory(GetVerboseMemoryFlag(), "After dilate");
307 // Set trachea with dilatation
308 trachea = SetBackground<MaskImageType, MaskImageType>
309 (trachea, working_mask, -1, this->GetForegroundValue(), true);
311 // Remove the trachea
312 working_mask = SetBackground<MaskImageType, MaskImageType>
313 (working_mask, working_mask, -1, this->GetBackgroundValue(), true);
316 working_mask = LabelizeAndSelectLabels<MaskImageType>
318 GetBackgroundValue(),
319 GetForegroundValue(),
321 GetMinimalComponentSize(),
322 GetLabelizeParameters3());
325 StopCurrentStep<MaskImageType>(working_mask);
328 //--------------------------------------------------------------------
329 //--------------------------------------------------------------------
330 PrintMemory(GetVerboseMemoryFlag(), "before autocropfilter");
331 if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
333 trachea = clitk::AutoCrop<MaskImageType>(trachea, GetBackgroundValue());
336 // Remove Padding region
337 typedef itk::CropImageFilter<MaskImageType, MaskImageType> CropFilterType;
338 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
339 cropFilter->SetInput(trachea);
340 cropFilter->SetLowerBoundaryCropSize(bounds);
341 cropFilter->SetUpperBoundaryCropSize(bounds);
342 cropFilter->Update();
343 trachea = cropFilter->GetOutput();
345 StopCurrentStep<MaskImageType>(trachea);
346 PrintMemory(GetVerboseMemoryFlag(), "after delete trachea");
348 PrintMemory(GetVerboseMemoryFlag(), "after delete trachea");
350 //--------------------------------------------------------------------
351 //--------------------------------------------------------------------
352 StartNewStep("Cropping lung");
353 PrintMemory(GetVerboseMemoryFlag(), "Before Autocropfilter");
355 working_mask = clitk::AutoCrop<MaskImageType>(working_mask, GetBackgroundValue());
358 // Remove Padding region
359 typedef itk::CropImageFilter<MaskImageType, MaskImageType> CropFilterType;
360 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
361 cropFilter->SetInput(working_mask);
362 cropFilter->SetLowerBoundaryCropSize(bounds);
363 cropFilter->SetUpperBoundaryCropSize(bounds);
364 cropFilter->Update();
365 working_mask = cropFilter->GetOutput();
367 StopCurrentStep<MaskImageType>(working_mask);
369 //--------------------------------------------------------------------
370 //--------------------------------------------------------------------
372 if (GetOpenCloseFlag()) {
373 StartNewStep("Open/Close");
374 PrintMemory(GetVerboseMemoryFlag(), "Before OpenClose");
376 // Structuring element
377 typedef itk::BinaryBallStructuringElement<InternalPixelType, ImageDimension> KernelType;
378 KernelType structuringElement;
379 structuringElement.SetRadius(GetOpenCloseRadius());
380 structuringElement.CreateStructuringElement();
383 typedef itk::BinaryMorphologicalOpeningImageFilter<MaskImageType, InternalImageType, KernelType> OpenFilterType;
384 typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
385 openFilter->SetInput(working_mask);
386 openFilter->SetBackgroundValue(GetBackgroundValue());
387 openFilter->SetForegroundValue(GetForegroundValue());
388 openFilter->SetKernel(structuringElement);
391 typedef itk::BinaryMorphologicalClosingImageFilter<MaskImageType, MaskImageType, KernelType> CloseFilterType;
392 typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
393 closeFilter->SetInput(openFilter->GetOutput());
394 closeFilter->SetSafeBorder(true);
395 closeFilter->SetForegroundValue(GetForegroundValue());
396 closeFilter->SetKernel(structuringElement);
397 closeFilter->Update();
398 working_mask = closeFilter->GetOutput();
401 //--------------------------------------------------------------------
402 //--------------------------------------------------------------------
404 if (GetFillHolesFlag()) {
405 StartNewStep("Fill Holes");
406 PrintMemory(GetVerboseMemoryFlag(), "Before Fill Holes");
407 typedef clitk::FillMaskFilter<MaskImageType> FillMaskFilterType;
408 typename FillMaskFilterType::Pointer fillMaskFilter = FillMaskFilterType::New();
409 fillMaskFilter->SetInput(working_mask);
410 fillMaskFilter->AddDirection(2);
411 //fillMaskFilter->AddDirection(1);
412 fillMaskFilter->Update();
413 working_mask = fillMaskFilter->GetOutput();
414 StopCurrentStep<MaskImageType>(working_mask);
417 if (GetSeparateLungsFlag()) {
418 //--------------------------------------------------------------------
419 //--------------------------------------------------------------------
420 StartNewStep("Separate Left/Right lungs");
421 PrintMemory(GetVerboseMemoryFlag(), "Before Separate");
423 working_mask = Labelize<MaskImageType>(working_mask,
424 GetBackgroundValue(),
426 GetMinimalComponentSize());
428 PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
431 typedef itk::StatisticsImageFilter<MaskImageType> StatisticsImageFilterType;
432 typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
433 statisticsImageFilter->SetInput(working_mask);
434 statisticsImageFilter->Update();
435 unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
436 working_mask = statisticsImageFilter->GetOutput();
438 PrintMemory(GetVerboseMemoryFlag(), "After count label");
440 // Decompose the first label
441 if (initialNumberOfLabels<2) {
442 // Structuring element radius
443 typename ImageType::SizeType radius;
444 for (unsigned int i=0;i<Dim;i++) radius[i]=1;
445 typedef clitk::DecomposeAndReconstructImageFilter<MaskImageType,MaskImageType> DecomposeAndReconstructFilterType;
446 typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
447 decomposeAndReconstructFilter->SetInput(working_mask);
448 decomposeAndReconstructFilter->SetVerbose(false);
449 decomposeAndReconstructFilter->SetRadius(radius);
450 decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
451 decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
452 decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
453 decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
454 decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
455 decomposeAndReconstructFilter->SetFullyConnected(true);
456 decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
457 decomposeAndReconstructFilter->Update();
458 working_mask = decomposeAndReconstructFilter->GetOutput();
460 PrintMemory(GetVerboseMemoryFlag(), "After decomposeAndReconstructFilter");
463 // Retain labels ('1' is largset lung, so right. '2' is left)
464 typedef itk::ThresholdImageFilter<MaskImageType> ThresholdImageFilterType;
465 typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
466 thresholdFilter->SetInput(working_mask);
467 thresholdFilter->ThresholdAbove(2);
468 thresholdFilter->SetOutsideValue(this->GetBackgroundValue());
469 thresholdFilter->Update();
470 working_mask = thresholdFilter->GetOutput();
471 StopCurrentStep<MaskImageType> (working_mask);
472 PrintMemory(GetVerboseMemoryFlag(), "After Thresholdfilter");
474 // Update output info
475 // output = working_mask;
476 //this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
478 // this->GetOutput(0)->SetRegions(working_mask->GetLargestPossibleRegion());
481 //--------------------------------------------------------------------
484 //--------------------------------------------------------------------
485 template <class TImageType>
487 clitk::ExtractLungFilter<TImageType>::
488 GenerateInputRequestedRegion() {
489 // DD("GenerateInputRequestedRegion (nothing?)");
491 //--------------------------------------------------------------------
494 //--------------------------------------------------------------------
495 template <class ImageType>
497 clitk::ExtractLungFilter<ImageType>::
501 // this->GraftOutput(output); // not SetNthOutput
502 this->GraftOutput(working_mask); // not SetNthOutput
503 // Store image filenames into AFDB
504 GetAFDB()->SetImageFilename("Lungs", this->GetOutputLungFilename());
505 GetAFDB()->SetImageFilename("Trachea", this->GetOutputTracheaFilename());
508 //--------------------------------------------------------------------
511 //--------------------------------------------------------------------
512 template <class ImageType>
514 clitk::ExtractLungFilter<ImageType>::
515 SearchForTracheaSeed(int skip)
517 if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)
518 // Restart the search until a seed is found, skipping more and more slices
521 // Search seed (parameters = UpperThresholdForTrachea)
522 static const unsigned int Dim = ImageType::ImageDimension;
523 typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
524 typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
525 typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
526 sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-skip-5;
527 sliceRegionSize[Dim-1]=5;
528 sliceRegion.SetSize(sliceRegionSize);
529 sliceRegion.SetIndex(sliceRegionIndex);
530 typedef itk::ImageRegionConstIterator<ImageType> IteratorType;
531 IteratorType it(working_input, sliceRegion);
533 while (!it.IsAtEnd()) {
534 if(it.Get() < GetUpperThresholdForTrachea() ) {
535 AddSeed(it.GetIndex());
536 // DD(it.GetIndex());
541 // if we do not found : restart
542 stop = (m_Seeds.size() != 0);
544 if (GetVerboseStepFlag()) {
545 std::cout << "\t No seed found this time. I skip some slices and restart." << std::endl;
547 if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
548 // we want to skip more than a half of the image, it is probably a bug
549 std::cerr << "Number of slices to skip to find trachea to high = " << skip << std::endl;
556 // DD(m_Seeds.size());
560 return (m_Seeds.size() != 0);
562 //--------------------------------------------------------------------
565 //--------------------------------------------------------------------
566 template <class ImageType>
568 clitk::ExtractLungFilter<ImageType>::
569 TracheaRegionGrowing()
571 // Explosion controlled region growing
572 PrintMemory(GetVerboseMemoryFlag(), "Before ExplosionControlledThresholdConnectedImageFilter");
573 typedef clitk::ExplosionControlledThresholdConnectedImageFilter<ImageType, MaskImageType> ImageFilterType;
574 typename ImageFilterType::Pointer f= ImageFilterType::New();
575 f->SetInput(working_input);
577 f->SetUpper(GetUpperThresholdForTrachea());
578 f->SetMinimumLowerThreshold(-2000);
579 // f->SetMaximumUpperThreshold(0); // MAYBE TO CHANGE ???
580 f->SetMaximumUpperThreshold(-800); // MAYBE TO CHANGE ???
581 f->SetAdaptLowerBorder(false);
582 f->SetAdaptUpperBorder(true);
583 f->SetMinimumSize(5000);
584 f->SetReplaceValue(1);
585 f->SetMultiplier(GetMultiplierForTrachea());
586 f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
587 f->SetMinimumThresholdStepSize(1);
588 f->SetVerbose(GetVerboseRegionGrowingFlag());
589 for(unsigned int i=0; i<m_Seeds.size();i++) {
590 f->AddSeed(m_Seeds[i]);
594 PrintMemory(GetVerboseMemoryFlag(), "After RG update");
596 // take first (main) connected component
597 trachea = Labelize<MaskImageType>(f->GetOutput(),
598 GetBackgroundValue(),
600 1000);//GetMinimalComponentSize());
601 PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
602 trachea = KeepLabels<MaskImageType>(trachea,
603 GetBackgroundValue(),
604 GetForegroundValue(),
606 PrintMemory(GetVerboseMemoryFlag(), "After KeepLabels");
608 //--------------------------------------------------------------------
611 //--------------------------------------------------------------------
612 template <class ImageType>
614 clitk::ExtractLungFilter<ImageType>::
615 ComputeTracheaVolume()
617 typedef itk::ImageRegionConstIterator<InternalImageType> IteratorType;
618 IteratorType iter(trachea, trachea->GetLargestPossibleRegion());
621 while (!iter.IsAtEnd()) {
622 if (iter.Get() == this->GetForegroundValue()) volume++;
626 double voxelsize = trachea->GetSpacing()[0]*trachea->GetSpacing()[1]*trachea->GetSpacing()[2];
627 return volume*voxelsize;
629 //--------------------------------------------------------------------
632 //--------------------------------------------------------------------
633 template <class ImageType>
635 clitk::ExtractLungFilter<ImageType>::
638 // Search for seed among n slices, skip some slices before starting
639 // if not found -> skip more and restart
641 // when seed found : perform region growing
642 // compute trachea volume
643 // if volume not plausible -> skip more slices and restart
647 int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed();
650 if (SearchForTracheaSeed(skip)) {
651 TracheaRegionGrowing();
652 volume = ComputeTracheaVolume()/1000; // assume mm3, so divide by 1000 to get cc
653 if (GetWriteStepFlag()) {
654 writeImage<MaskImageType>(trachea, "step-trachea-"+toString(skip)+".mhd");
656 if (GetTracheaVolumeMustBeCheckedFlag()) {
657 if ((volume > 10) && (volume < 65 )) { // depend on image size ...
658 // Typical volume 22.59 cm 3 (± 7.69 cm 3 ) [ Leader 2004 ]
659 if (GetVerboseStepFlag()) {
660 std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl;
664 if (GetVerboseStepFlag()) {
665 std::cout << "\t The volume of the trachea (" << volume
666 << " cc) seems not correct. I skip some slices (" << skip << ") and restart to find seeds."
671 // empty the list of seed
674 if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
675 // we want to skip more than a half of the image, it is probably a bug
676 std::cerr << "2 : Number of slices to skip to find trachea too high = " << skip << std::endl;
688 StopCurrentStep<MaskImageType>(trachea);
690 else { // Trachea not found
691 this->SetWarning("* WARNING * No seed found for trachea.");
695 //--------------------------------------------------------------------
697 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX