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"
41 #include <itkBinaryBallStructuringElement.h>
42 #include "itkStatisticsLabelObject.h"
43 #include "itkLabelMap.h"
44 #include "itkLabelImageToShapeLabelMapFilter.h"
45 #include "itkLabelMapToLabelImageFilter.h"
46 #include "itkExtractImageFilter.h"
47 #include "itkOrientImageFilter.h"
48 #include "itkSpatialOrientationAdapter.h"
49 #include "itkImageDuplicator.h"
52 //--------------------------------------------------------------------
53 template <class ImageType>
54 clitk::ExtractLungFilter<ImageType>::
57 clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
58 itk::ImageToImageFilter<ImageType, MaskImageType>()
61 m_MaxSeedNumber = 500;
63 // Default global options
64 this->SetNumberOfRequiredInputs(1);
65 SetPatientMaskBackgroundValue(0);
66 SetBackgroundValue(0); // Must be zero
67 SetForegroundValue(1);
68 SetMinimalComponentSize(100);
69 VerboseRegionGrowingFlagOff();
71 // Step 1 default values
72 SetUpperThreshold(-300);
73 SetLowerThreshold(-1000);
74 UseLowerThresholdOff();
75 LabelParamType * p1 = new LabelParamType;
78 p1->AddLabelToRemove(2);
79 SetLabelizeParameters1(p1);
81 // Step 2 default values
82 SetTracheaSeedAlgorithm(0);
83 SetUpperThresholdForTrachea(-700);
84 SetMultiplierForTrachea(5);
85 SetThresholdStepSizeForTrachea(64);
86 SetNumberOfSlicesToSkipBeforeSearchingSeed(0);
87 TracheaVolumeMustBeCheckedFlagOn();
89 SetMaxElongation(0.5);
90 SetSeedPreProcessingThreshold(-400);
93 // Step 3 default values
94 SetNumberOfHistogramBins(500);
95 LabelParamType * p2 = new LabelParamType;
98 // p->AddLabelToRemove(2); // No label to remove by default
99 SetLabelizeParameters2(p2);
101 // Step 4 default values
102 SetRadiusForTrachea(1);
103 LabelParamType * p3 = new LabelParamType;
106 p3->UseLastKeepOff();
107 SetLabelizeParameters3(p3);
111 SetOpenCloseRadius(1);
117 //--------------------------------------------------------------------
120 //--------------------------------------------------------------------
121 template <class ImageType>
123 clitk::ExtractLungFilter<ImageType>::
124 SetInput(const ImageType * image)
126 this->SetNthInput(0, const_cast<ImageType *>(image));
128 //--------------------------------------------------------------------
131 //--------------------------------------------------------------------
132 template <class ImageType>
134 clitk::ExtractLungFilter<ImageType>::
135 AddSeed(InternalIndexType s)
137 m_Seeds.push_back(s);
139 //--------------------------------------------------------------------
142 //--------------------------------------------------------------------
143 template <class ImageType>
145 clitk::ExtractLungFilter<ImageType>::
146 GenerateOutputInformation()
148 clitk::PrintMemory(GetVerboseMemoryFlag(), "Initial memory"); // OK
149 Superclass::GenerateOutputInformation();
154 // Get input pointers
155 input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
156 patient = GetAFDB()->template GetImage <MaskImageType>("Patient");
157 PrintMemory(GetVerboseMemoryFlag(), "After reading patient"); // OK
159 //--------------------------------------------------------------------
160 //--------------------------------------------------------------------
161 // Crop input like patient image (must have the same spacing)
162 // It copy the input if the same are the same
163 StartNewStep("Copy and crop input image to 'patient' extends");
164 typedef clitk::CropLikeImageFilter<ImageType> CropImageFilter;
165 typename CropImageFilter::Pointer cropFilter = CropImageFilter::New();
166 // cropFilter->ReleaseDataFlagOn(); // NO=seg fault !!
167 cropFilter->SetInput(input);
168 cropFilter->SetCropLikeImage(patient);
169 cropFilter->Update();
170 working_input = cropFilter->GetOutput();
171 // cropFilter->Delete(); // NO !!!! if yes, sg fault, Cropfilter is buggy !??
172 StopCurrentStep<ImageType>(working_input);
173 PrintMemory(GetVerboseMemoryFlag(), "After crop"); // OK, slightly more than a copy of input
175 //--------------------------------------------------------------------
176 //--------------------------------------------------------------------
177 StartNewStep("Set background to initial image");
178 working_input = SetBackground<ImageType, MaskImageType>
179 (working_input, patient, GetPatientMaskBackgroundValue(), -1000, true);
181 // Pad images with air to prevent patient touching the image border
182 static const unsigned int Dim = ImageType::ImageDimension;
183 typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadFilterType;
184 typename PadFilterType::Pointer padFilter = PadFilterType::New();
185 padFilter->SetInput(working_input);
186 padFilter->SetConstant(-1000);
187 typename ImageType::SizeType bounds;
188 for (unsigned i = 0; i < Dim - 1; ++i)
191 padFilter->SetPadLowerBound(bounds);
192 padFilter->SetPadUpperBound(bounds);
194 working_input = padFilter->GetOutput();
196 StopCurrentStep<ImageType>(working_input);
197 PrintMemory(GetVerboseMemoryFlag(), "After set bg"); // OK, additional mem = 0
200 // We do not need patient mask anymore, release memory
201 GetAFDB()->template ReleaseImage<MaskImageType>("Patient");
202 DD(patient->GetReferenceCount());
203 PrintMemory(GetVerboseMemoryFlag(), "After delete patient"); // OK, additional mem = 0
205 PrintMemory(GetVerboseMemoryFlag(), "After delete patient"); // OK, additional mem = 0
208 //--------------------------------------------------------------------
209 //--------------------------------------------------------------------
210 StartNewStep("Remove Air");
212 if (m_UseLowerThreshold) {
213 if (m_LowerThreshold > m_UpperThreshold) {
214 clitkExceptionMacro("lower threshold cannot be greater than upper threshold.");
217 // Threshold to get air
218 typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> InputBinarizeFilterType;
219 typename InputBinarizeFilterType::Pointer binarizeFilter = InputBinarizeFilterType::New();
220 binarizeFilter->SetInput(working_input);
221 if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(m_LowerThreshold);
222 // binarizeFilter->CanRunInPlace() is false
223 binarizeFilter->SetUpperThreshold(m_UpperThreshold);
224 binarizeFilter->SetInsideValue(this->GetForegroundValue());
225 binarizeFilter->SetOutsideValue(this->GetBackgroundValue());
226 binarizeFilter->Update();
227 working_mask = binarizeFilter->GetOutput();
228 PrintMemory(GetVerboseMemoryFlag(), "After Binarizefilter"); // OK, additional mem is one mask image
230 // Labelize and keep right labels
232 Labelize<MaskImageType>
233 (working_mask, GetBackgroundValue(), true, GetMinimalComponentSize());
234 PrintMemory(GetVerboseMemoryFlag(), "After Labelize"); // BUG ? additional mem around 1 time the input ?
236 working_mask = RemoveLabels<MaskImageType>
237 (working_mask, GetBackgroundValue(), GetLabelizeParameters1()->GetLabelsToRemove());
238 PrintMemory(GetVerboseMemoryFlag(), "After RemoveLabels"); // OK additional mem = 0
240 working_mask = KeepLabels<MaskImageType>
242 GetBackgroundValue(),
243 GetForegroundValue(),
244 GetLabelizeParameters1()->GetFirstKeep(),
245 GetLabelizeParameters1()->GetLastKeep(),
246 GetLabelizeParameters1()->GetUseLastKeep());
247 PrintMemory(GetVerboseMemoryFlag(), "After KeepLabels to create the 'air'");
250 working_input = SetBackground<ImageType, MaskImageType>
251 (working_input, working_mask, this->GetForegroundValue(), this->GetBackgroundValue(), true);
252 PrintMemory(GetVerboseMemoryFlag(), "After SetBackground");
255 StopCurrentStep<ImageType>(working_input);
257 //--------------------------------------------------------------------
258 //--------------------------------------------------------------------
259 StartNewStep("Search for the trachea");
261 PrintMemory(GetVerboseMemoryFlag(), "After SearchForTrachea");
262 if (m_Seeds.empty()) {
263 clitkExceptionMacro("No seeds for trachea... Aborting.");
266 //--------------------------------------------------------------------
267 //--------------------------------------------------------------------
268 StartNewStep("Extract the lung with Otsu filter");
269 // Automated Otsu thresholding and relabeling
270 typedef itk::OtsuThresholdImageFilter<ImageType,MaskImageType> OtsuThresholdImageFilterType;
271 typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New();
272 otsuFilter->SetInput(working_input);
273 otsuFilter->SetNumberOfHistogramBins(GetNumberOfHistogramBins());
274 otsuFilter->SetInsideValue(this->GetForegroundValue());
275 otsuFilter->SetOutsideValue(this->GetBackgroundValue());
276 otsuFilter->Update();
277 working_mask = otsuFilter->GetOutput();
280 StopCurrentStep<MaskImageType>(working_mask);
281 PrintMemory(GetVerboseMemoryFlag(), "After Otsufilter");
283 //--------------------------------------------------------------------
284 //--------------------------------------------------------------------
285 StartNewStep("Select labels");
287 working_mask = LabelizeAndSelectLabels<MaskImageType>
289 GetBackgroundValue(),
290 GetForegroundValue(),
292 GetMinimalComponentSize(),
293 GetLabelizeParameters2());
296 StopCurrentStep<MaskImageType>(working_mask);
297 PrintMemory(GetVerboseMemoryFlag(), "After LabelizeAndSelectLabels");
299 //--------------------------------------------------------------------
300 //--------------------------------------------------------------------
301 if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
302 StartNewStep("Remove the trachea");
304 working_mask = SetBackground<MaskImageType, MaskImageType>
305 (working_mask, trachea, 1, -1, true);
306 PrintMemory(GetVerboseMemoryFlag(), "After SetBackground");
308 // Dilate the trachea
309 static const unsigned int Dim = ImageType::ImageDimension;
310 typedef itk::BinaryBallStructuringElement<InternalPixelType, Dim> KernelType;
311 KernelType structuringElement;
312 structuringElement.SetRadius(GetRadiusForTrachea());
313 structuringElement.CreateStructuringElement();
314 typedef clitk::ConditionalBinaryDilateImageFilter<MaskImageType, MaskImageType, KernelType> ConditionalBinaryDilateImageFilterType;
315 typename ConditionalBinaryDilateImageFilterType::Pointer dilateFilter = ConditionalBinaryDilateImageFilterType::New();
316 dilateFilter->SetBoundaryToForeground(false);
317 dilateFilter->SetKernel(structuringElement);
318 dilateFilter->SetBackgroundValue (1);
319 dilateFilter->SetForegroundValue (-1);
320 dilateFilter->SetInput (working_mask);
321 dilateFilter->Update();
322 working_mask = dilateFilter->GetOutput();
323 PrintMemory(GetVerboseMemoryFlag(), "After dilate");
325 // Set trachea with dilatation
326 trachea = SetBackground<MaskImageType, MaskImageType>
327 (trachea, working_mask, -1, this->GetForegroundValue(), true);
329 // Remove the trachea
330 working_mask = SetBackground<MaskImageType, MaskImageType>
331 (working_mask, working_mask, -1, this->GetBackgroundValue(), true);
334 working_mask = LabelizeAndSelectLabels<MaskImageType>
336 GetBackgroundValue(),
337 GetForegroundValue(),
339 GetMinimalComponentSize(),
340 GetLabelizeParameters3());
343 StopCurrentStep<MaskImageType>(working_mask);
346 //--------------------------------------------------------------------
347 //--------------------------------------------------------------------
348 PrintMemory(GetVerboseMemoryFlag(), "before autocropfilter");
349 if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
351 trachea = clitk::AutoCrop<MaskImageType>(trachea, GetBackgroundValue());
354 // Remove Padding region
355 typedef itk::CropImageFilter<MaskImageType, MaskImageType> CropFilterType;
356 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
357 cropFilter->SetInput(trachea);
358 cropFilter->SetLowerBoundaryCropSize(bounds);
359 cropFilter->SetUpperBoundaryCropSize(bounds);
360 cropFilter->Update();
361 trachea = cropFilter->GetOutput();
363 StopCurrentStep<MaskImageType>(trachea);
364 PrintMemory(GetVerboseMemoryFlag(), "after delete trachea");
366 PrintMemory(GetVerboseMemoryFlag(), "after delete trachea");
368 //--------------------------------------------------------------------
369 //--------------------------------------------------------------------
370 StartNewStep("Cropping lung");
371 PrintMemory(GetVerboseMemoryFlag(), "Before Autocropfilter");
373 working_mask = clitk::AutoCrop<MaskImageType>(working_mask, GetBackgroundValue());
376 // Remove Padding region
377 typedef itk::CropImageFilter<MaskImageType, MaskImageType> CropFilterType;
378 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
379 cropFilter->SetInput(working_mask);
380 cropFilter->SetLowerBoundaryCropSize(bounds);
381 cropFilter->SetUpperBoundaryCropSize(bounds);
382 cropFilter->Update();
383 working_mask = cropFilter->GetOutput();
385 StopCurrentStep<MaskImageType>(working_mask);
387 //--------------------------------------------------------------------
388 //--------------------------------------------------------------------
390 if (GetOpenCloseFlag()) {
391 StartNewStep("Open/Close");
392 PrintMemory(GetVerboseMemoryFlag(), "Before OpenClose");
394 // Structuring element
395 typedef itk::BinaryBallStructuringElement<InternalPixelType, ImageDimension> KernelType;
396 KernelType structuringElement;
397 structuringElement.SetRadius(GetOpenCloseRadius());
398 structuringElement.CreateStructuringElement();
401 typedef itk::BinaryMorphologicalOpeningImageFilter<MaskImageType, InternalImageType, KernelType> OpenFilterType;
402 typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
403 openFilter->SetInput(working_mask);
404 openFilter->SetBackgroundValue(GetBackgroundValue());
405 openFilter->SetForegroundValue(GetForegroundValue());
406 openFilter->SetKernel(structuringElement);
409 typedef itk::BinaryMorphologicalClosingImageFilter<MaskImageType, MaskImageType, KernelType> CloseFilterType;
410 typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
411 closeFilter->SetInput(openFilter->GetOutput());
412 closeFilter->SetSafeBorder(true);
413 closeFilter->SetForegroundValue(GetForegroundValue());
414 closeFilter->SetKernel(structuringElement);
415 closeFilter->Update();
416 working_mask = closeFilter->GetOutput();
419 //--------------------------------------------------------------------
420 //--------------------------------------------------------------------
422 if (GetFillHolesFlag()) {
423 StartNewStep("Fill Holes");
424 PrintMemory(GetVerboseMemoryFlag(), "Before Fill Holes");
425 typedef clitk::FillMaskFilter<MaskImageType> FillMaskFilterType;
426 typename FillMaskFilterType::Pointer fillMaskFilter = FillMaskFilterType::New();
427 fillMaskFilter->SetInput(working_mask);
428 fillMaskFilter->AddDirection(2);
429 //fillMaskFilter->AddDirection(1);
430 fillMaskFilter->Update();
431 working_mask = fillMaskFilter->GetOutput();
432 StopCurrentStep<MaskImageType>(working_mask);
435 if (GetSeparateLungsFlag()) {
436 //--------------------------------------------------------------------
437 //--------------------------------------------------------------------
438 StartNewStep("Separate Left/Right lungs");
439 PrintMemory(GetVerboseMemoryFlag(), "Before Separate");
441 working_mask = Labelize<MaskImageType>(working_mask,
442 GetBackgroundValue(),
444 GetMinimalComponentSize());
446 PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
449 typedef itk::StatisticsImageFilter<MaskImageType> StatisticsImageFilterType;
450 typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
451 statisticsImageFilter->SetInput(working_mask);
452 statisticsImageFilter->Update();
453 unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
454 working_mask = statisticsImageFilter->GetOutput();
456 PrintMemory(GetVerboseMemoryFlag(), "After count label");
458 // Decompose the first label
459 if (initialNumberOfLabels<2) {
460 // Structuring element radius
461 typename ImageType::SizeType radius;
462 for (unsigned int i=0;i<Dim;i++) radius[i]=1;
463 typedef clitk::DecomposeAndReconstructImageFilter<MaskImageType,MaskImageType> DecomposeAndReconstructFilterType;
464 typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
465 decomposeAndReconstructFilter->SetInput(working_mask);
466 decomposeAndReconstructFilter->SetVerbose(false);
467 decomposeAndReconstructFilter->SetRadius(radius);
468 decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
469 decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
470 decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
471 decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
472 decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
473 decomposeAndReconstructFilter->SetFullyConnected(true);
474 decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
475 decomposeAndReconstructFilter->Update();
476 working_mask = decomposeAndReconstructFilter->GetOutput();
478 PrintMemory(GetVerboseMemoryFlag(), "After decomposeAndReconstructFilter");
481 // Retain labels ('1' is largset lung, so right. '2' is left)
482 typedef itk::ThresholdImageFilter<MaskImageType> ThresholdImageFilterType;
483 typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
484 thresholdFilter->SetInput(working_mask);
485 thresholdFilter->ThresholdAbove(2);
486 thresholdFilter->SetOutsideValue(this->GetBackgroundValue());
487 thresholdFilter->Update();
488 working_mask = thresholdFilter->GetOutput();
489 StopCurrentStep<MaskImageType> (working_mask);
490 PrintMemory(GetVerboseMemoryFlag(), "After Thresholdfilter");
492 // Update output info
493 // output = working_mask;
494 //this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
496 // this->GetOutput(0)->SetRegions(working_mask->GetLargestPossibleRegion());
499 //--------------------------------------------------------------------
502 //--------------------------------------------------------------------
503 template <class TImageType>
505 clitk::ExtractLungFilter<TImageType>::
506 GenerateInputRequestedRegion() {
507 // DD("GenerateInputRequestedRegion (nothing?)");
509 //--------------------------------------------------------------------
512 //--------------------------------------------------------------------
513 template <class ImageType>
515 clitk::ExtractLungFilter<ImageType>::
519 // this->GraftOutput(output); // not SetNthOutput
520 this->GraftOutput(working_mask); // not SetNthOutput
521 // Store image filenames into AFDB
522 GetAFDB()->SetImageFilename("Lungs", this->GetOutputLungFilename());
523 GetAFDB()->SetImageFilename("Trachea", this->GetOutputTracheaFilename());
526 //--------------------------------------------------------------------
529 //--------------------------------------------------------------------
530 template <class ImageType>
532 clitk::ExtractLungFilter<ImageType>::
533 SearchForTracheaSeed(int skip)
535 if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)
536 // Restart the search until a seed is found, skipping more and more slices
539 // Search seed (parameters = UpperThresholdForTrachea)
540 static const unsigned int Dim = ImageType::ImageDimension;
541 typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
542 typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
543 typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
544 sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-skip-5;
545 sliceRegionSize[Dim-1]=5;
546 sliceRegion.SetSize(sliceRegionSize);
547 sliceRegion.SetIndex(sliceRegionIndex);
548 typedef itk::ImageRegionConstIterator<ImageType> IteratorType;
549 IteratorType it(working_input, sliceRegion);
551 while (!it.IsAtEnd()) {
552 if(it.Get() < GetUpperThresholdForTrachea() ) {
553 AddSeed(it.GetIndex());
554 // DD(it.GetIndex());
559 // if we do not found : restart
560 stop = (m_Seeds.size() != 0);
562 if (GetVerboseStepFlag()) {
563 std::cout << "\t No seed found this time. I skip some slices and restart." << std::endl;
565 if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
566 // we want to skip more than a half of the image, it is probably a bug
567 std::cerr << "Number of slices to skip to find trachea to high = " << skip << std::endl;
574 // DD(m_Seeds.size());
578 return (m_Seeds.size() != 0);
580 //--------------------------------------------------------------------
583 bool is_orientation_superior(itk::SpatialOrientation::ValidCoordinateOrientationFlags orientation)
585 itk::SpatialOrientation::CoordinateTerms sup = itk::SpatialOrientation::ITK_COORDINATE_Superior;
586 bool primary = (orientation & 0x0000ff) == sup;
587 bool secondary = ((orientation & 0x00ff00) >> 8) == sup;
588 bool tertiary = ((orientation & 0xff0000) >> 16) == sup;
589 return primary || secondary || tertiary;
592 //--------------------------------------------------------------------
593 template <class ImageType>
595 clitk::ExtractLungFilter<ImageType>::
596 SearchForTracheaSeed2(int numberOfSlices)
598 if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)
599 if (GetVerboseRegionGrowingFlag())
600 std::cout << "SearchForTracheaSeed2(" << numberOfSlices << ", " << GetMaxElongation() << ")" << std::endl;
602 typedef unsigned char MaskPixelType;
603 typedef itk::Image<MaskPixelType, ImageType::ImageDimension> MaskImageType;
604 typedef itk::Image<typename MaskImageType::PixelType, 2> MaskImageType2D;
605 typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> ThresholdFilterType;
606 typedef itk::BinaryBallStructuringElement<MaskPixelType, MaskImageType::ImageDimension> KernelType;
607 typedef itk::BinaryMorphologicalClosingImageFilter<MaskImageType, MaskImageType, KernelType> ClosingFilterType;
608 typedef itk::BinaryMorphologicalOpeningImageFilter<MaskImageType, MaskImageType, KernelType> OpeningFilterType;
609 typedef itk::ExtractImageFilter<MaskImageType, MaskImageType2D> SlicerFilterType;
610 typedef itk::ConnectedComponentImageFilter<MaskImageType2D, MaskImageType2D> LabelFilterType;
611 typedef itk::ShapeLabelObject<MaskPixelType, MaskImageType2D::ImageDimension> ShapeLabelType;
612 typedef itk::LabelMap<ShapeLabelType> LabelImageType;
613 typedef itk::LabelImageToShapeLabelMapFilter<MaskImageType2D, LabelImageType> ImageToLabelMapFilterType;
614 typedef itk::LabelMapToLabelImageFilter<LabelImageType, MaskImageType2D> LabelMapToImageFilterType;
615 typedef itk::ImageFileWriter<MaskImageType2D> FileWriterType;
617 // threshold to isolate airawys and lungs
618 typename ThresholdFilterType::Pointer threshold = ThresholdFilterType::New();
619 threshold->SetLowerThreshold(-2000);
620 threshold->SetUpperThreshold(GetSeedPreProcessingThreshold());
621 threshold->SetInput(working_input);
624 KernelType kernel_closing, kernel_opening;
626 // remove small noise
627 typename OpeningFilterType::Pointer opening = OpeningFilterType::New();
628 kernel_opening.SetRadius(1);
629 opening->SetKernel(kernel_opening);
630 opening->SetInput(threshold->GetOutput());
633 typename SlicerFilterType::Pointer slicer = SlicerFilterType::New();
634 slicer->SetInput(opening->GetOutput());
637 typename LabelFilterType::Pointer label_filter = LabelFilterType::New();
638 label_filter->SetInput(slicer->GetOutput());
640 // extract shape information from labels
641 typename ImageToLabelMapFilterType::Pointer label_to_map_filter = ImageToLabelMapFilterType::New();
642 label_to_map_filter->SetInput(label_filter->GetOutput());
644 typename LabelMapToImageFilterType::Pointer map_to_label_filter = LabelMapToImageFilterType::New();
645 typename FileWriterType::Pointer writer = FileWriterType::New();
647 typename ImageType::IndexType index;
648 typename ImageType::RegionType region = working_input->GetLargestPossibleRegion();
649 typename ImageType::SizeType size = region.GetSize();
650 typename ImageType::SpacingType spacing = working_input->GetSpacing();
651 typename ImageType::PointType origin = working_input->GetOrigin();
653 int nslices = min(numberOfSlices, size[2]);
654 int start = 0, increment = 1;
655 itk::SpatialOrientationAdapter orientation;
656 typename ImageType::DirectionType dir = working_input->GetDirection();
657 if (!is_orientation_superior(orientation.FromDirectionCosines(dir))) {
662 typename MaskImageType::PointType image_centre;
663 image_centre[0] = size[0]/2;
664 image_centre[1] = size[1]/2;
667 typedef InternalIndexType SeedType;
668 SeedType trachea_centre, shape_centre, max_e_centre, prev_e_centre;
669 typedef std::list<SeedType> PointListType;
670 typedef std::list<PointListType> SequenceListType;
671 PointListType* current_sequence = NULL;
672 SequenceListType sequence_list;
674 prev_e_centre.Fill(0);
675 std::ostringstream file_name;
676 index[0] = index[1] = 0;
677 size[0] = size[1] = 512;
683 region.SetIndex(index);
684 region.SetSize(size);
685 slicer->SetExtractionRegion(region);
687 label_filter->SetInput(slicer->GetOutput());
688 label_filter->Update();
690 label_to_map_filter->SetInput(label_filter->GetOutput());
691 label_to_map_filter->Update();
692 typename LabelImageType::Pointer label_map = label_to_map_filter->GetOutput();
694 if (GetWriteStepFlag()) {
695 map_to_label_filter->SetInput(label_map);
696 writer->SetInput(map_to_label_filter->GetOutput());
698 file_name << "labels_";
701 file_name << index[2] << ".mhd";
702 writer->SetFileName(file_name.str().c_str());
706 typename ShapeLabelType::Pointer shape, max_e_shape;
707 double max_elongation = GetMaxElongation();
708 double max_size = size[0]*size[1]/128;
711 unsigned int nlables = label_map->GetNumberOfLabelObjects();
712 for (unsigned int j = 0; j < nlables; j++) {
713 shape = label_map->GetNthLabelObject(j);
714 if (shape->Size() > 150 && shape->Size() <= max_size) {
715 #if ITK_VERSION_MAJOR < 4
716 double e = 1 - 1/shape->GetBinaryElongation();
718 double e = 1 - 1/shape->GetElongation();
720 //double area = 1 - r->Area() ;
721 if (e < max_elongation) {
723 shape_centre[0] = (shape->GetCentroid()[0] - origin[0])/spacing[0];
724 shape_centre[1] = (shape->GetCentroid()[1] - origin[1])/spacing[1];
725 shape_centre[2] = index[2];
726 //double d = 1 - (shape_centre - image_centre).Magnitude()/max_dist;
727 double dx = shape_centre[0] - image_centre[0];
728 double d = 1 - dx*2/size[0];
734 max_e_centre = shape_centre;
742 itk::Point<typename SeedType::IndexValueType, ImageType::ImageDimension> p1, p2;
743 p1[0] = max_e_centre[0];
744 p1[1] = max_e_centre[1];
745 p1[2] = max_e_centre[2];
747 p2[0] = prev_e_centre[0];
748 p2[1] = prev_e_centre[1];
749 p2[2] = prev_e_centre[2];
751 double mag = (p2 - p1).GetNorm();
752 if (GetVerboseRegionGrowingFlag()) {
754 cout << index[2] << ": ";
755 cout << "region(" << max_e_centre[0] << ", " << max_e_centre[1] << ", " << max_e_centre[2] << "); ";
756 cout << "prev_region(" << prev_e_centre[0] << ", " << prev_e_centre[1] << ", " << prev_e_centre[2] << "); ";
757 cout << "mag(" << mag << "); " << endl;
762 PointListType point_list;
763 point_list.push_back(max_e_centre);
764 sequence_list.push_back(point_list);
765 current_sequence = &sequence_list.back();
767 else if (current_sequence)
768 current_sequence->push_back(max_e_centre);
770 prev_e_centre= max_e_centre;
773 if (GetVerboseRegionGrowingFlag()) {
774 cout << "No shapes found at slice " << index[2] << std::endl;
780 for (typename SequenceListType::iterator s = sequence_list.begin(); s != sequence_list.end(); s++)
782 if (s->size() > longest)
785 trachea_centre = s->front();
790 if (GetVerboseRegionGrowingFlag())
791 std::cout << "seed at: " << trachea_centre << std::endl;
792 m_Seeds.push_back(trachea_centre);
796 return (m_Seeds.size() != 0);
798 //--------------------------------------------------------------------
800 //--------------------------------------------------------------------
801 template <class ImageType>
803 clitk::ExtractLungFilter<ImageType>::
804 TracheaRegionGrowing()
806 // Explosion controlled region growing
807 PrintMemory(GetVerboseMemoryFlag(), "Before ExplosionControlledThresholdConnectedImageFilter");
808 typedef clitk::ExplosionControlledThresholdConnectedImageFilter<ImageType, MaskImageType> ImageFilterType;
809 typename ImageFilterType::Pointer f= ImageFilterType::New();
810 f->SetInput(working_input);
812 f->SetUpper(GetUpperThresholdForTrachea());
813 f->SetMinimumLowerThreshold(-2000);
814 // f->SetMaximumUpperThreshold(0); // MAYBE TO CHANGE ???
815 f->SetMaximumUpperThreshold(-700); // MAYBE TO CHANGE ???
816 f->SetAdaptLowerBorder(false);
817 f->SetAdaptUpperBorder(true);
818 f->SetMinimumSize(5000);
819 f->SetReplaceValue(1);
820 f->SetMultiplier(GetMultiplierForTrachea());
821 f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
822 f->SetMinimumThresholdStepSize(1);
823 f->SetVerbose(GetVerboseRegionGrowingFlag());
824 for(unsigned int i=0; i<m_Seeds.size();i++) {
825 f->AddSeed(m_Seeds[i]);
829 PrintMemory(GetVerboseMemoryFlag(), "After RG update");
831 // take first (main) connected component
832 trachea = Labelize<MaskImageType>(f->GetOutput(),
833 GetBackgroundValue(),
835 1000);//GetMinimalComponentSize());
836 PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
837 trachea = KeepLabels<MaskImageType>(trachea,
838 GetBackgroundValue(),
839 GetForegroundValue(),
841 PrintMemory(GetVerboseMemoryFlag(), "After KeepLabels");
843 //--------------------------------------------------------------------
846 //--------------------------------------------------------------------
847 template <class ImageType>
849 clitk::ExtractLungFilter<ImageType>::
850 ComputeTracheaVolume()
852 typedef itk::ImageRegionConstIterator<InternalImageType> IteratorType;
853 IteratorType iter(trachea, trachea->GetLargestPossibleRegion());
856 while (!iter.IsAtEnd()) {
857 if (iter.Get() == this->GetForegroundValue()) volume++;
861 double voxelsize = trachea->GetSpacing()[0]*trachea->GetSpacing()[1]*trachea->GetSpacing()[2];
862 return volume*voxelsize;
864 //--------------------------------------------------------------------
867 //--------------------------------------------------------------------
868 template <class ImageType>
870 clitk::ExtractLungFilter<ImageType>::
873 // Search for seed among n slices, skip some slices before starting
874 // if not found -> skip more and restart
876 // when seed found : perform region growing
877 // compute trachea volume
878 // if volume not plausible -> skip more slices and restart
883 int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed();
886 if (GetTracheaSeedAlgorithm() == 0)
887 has_seed = SearchForTracheaSeed(skip);
889 has_seed = SearchForTracheaSeed2(GetNumSlices());
892 TracheaRegionGrowing();
893 volume = ComputeTracheaVolume()/1000; // assume mm3, so divide by 1000 to get cc
894 if (GetWriteStepFlag()) {
895 writeImage<MaskImageType>(trachea, "step-trachea-"+toString(skip)+".mhd");
897 if (GetTracheaVolumeMustBeCheckedFlag()) {
898 if ((volume > 10) && (volume < 65 )) { // depend on image size ...
899 // Typical volume 22.59 cm 3 (± 7.69 cm 3 ) [ Leader 2004 ]
900 if (GetVerboseStepFlag())
902 std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl;
906 if (GetTracheaSeedAlgorithm() == 0) {
907 if (GetVerboseStepFlag()) {
908 std::cout << "\t The volume of the trachea (" << volume
909 << " cc) seems not correct. I skip some slices (" << skip << ") and restart to find seeds."
914 // empty the list of seed
917 if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
918 // we want to skip more than a half of the image, it is probably a bug
919 std::cerr << "2 : Number of slices to skip to find trachea too high = " << skip << std::endl;
931 StopCurrentStep<MaskImageType>(trachea);
933 else { // Trachea not found
934 this->SetWarning("* WARNING * No seed found for trachea.");
938 //--------------------------------------------------------------------
940 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX