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");
263 //--------------------------------------------------------------------
264 //--------------------------------------------------------------------
265 StartNewStep("Extract the lung with Otsu filter");
266 // Automated Otsu thresholding and relabeling
267 typedef itk::OtsuThresholdImageFilter<ImageType,MaskImageType> OtsuThresholdImageFilterType;
268 typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New();
269 otsuFilter->SetInput(working_input);
270 otsuFilter->SetNumberOfHistogramBins(GetNumberOfHistogramBins());
271 otsuFilter->SetInsideValue(this->GetForegroundValue());
272 otsuFilter->SetOutsideValue(this->GetBackgroundValue());
273 otsuFilter->Update();
274 working_mask = otsuFilter->GetOutput();
277 StopCurrentStep<MaskImageType>(working_mask);
278 PrintMemory(GetVerboseMemoryFlag(), "After Otsufilter");
280 //--------------------------------------------------------------------
281 //--------------------------------------------------------------------
282 StartNewStep("Select labels");
284 working_mask = LabelizeAndSelectLabels<MaskImageType>
286 GetBackgroundValue(),
287 GetForegroundValue(),
289 GetMinimalComponentSize(),
290 GetLabelizeParameters2());
293 StopCurrentStep<MaskImageType>(working_mask);
294 PrintMemory(GetVerboseMemoryFlag(), "After LabelizeAndSelectLabels");
296 //--------------------------------------------------------------------
297 //--------------------------------------------------------------------
298 if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
299 StartNewStep("Remove the trachea");
301 working_mask = SetBackground<MaskImageType, MaskImageType>
302 (working_mask, trachea, 1, -1, true);
303 PrintMemory(GetVerboseMemoryFlag(), "After SetBackground");
305 // Dilate the trachea
306 static const unsigned int Dim = ImageType::ImageDimension;
307 typedef itk::BinaryBallStructuringElement<InternalPixelType, Dim> KernelType;
308 KernelType structuringElement;
309 structuringElement.SetRadius(GetRadiusForTrachea());
310 structuringElement.CreateStructuringElement();
311 typedef clitk::ConditionalBinaryDilateImageFilter<MaskImageType, MaskImageType, KernelType> ConditionalBinaryDilateImageFilterType;
312 typename ConditionalBinaryDilateImageFilterType::Pointer dilateFilter = ConditionalBinaryDilateImageFilterType::New();
313 dilateFilter->SetBoundaryToForeground(false);
314 dilateFilter->SetKernel(structuringElement);
315 dilateFilter->SetBackgroundValue (1);
316 dilateFilter->SetForegroundValue (-1);
317 dilateFilter->SetInput (working_mask);
318 dilateFilter->Update();
319 working_mask = dilateFilter->GetOutput();
320 PrintMemory(GetVerboseMemoryFlag(), "After dilate");
322 // Set trachea with dilatation
323 trachea = SetBackground<MaskImageType, MaskImageType>
324 (trachea, working_mask, -1, this->GetForegroundValue(), true);
326 // Remove the trachea
327 working_mask = SetBackground<MaskImageType, MaskImageType>
328 (working_mask, working_mask, -1, this->GetBackgroundValue(), true);
331 working_mask = LabelizeAndSelectLabels<MaskImageType>
333 GetBackgroundValue(),
334 GetForegroundValue(),
336 GetMinimalComponentSize(),
337 GetLabelizeParameters3());
340 StopCurrentStep<MaskImageType>(working_mask);
343 //--------------------------------------------------------------------
344 //--------------------------------------------------------------------
345 PrintMemory(GetVerboseMemoryFlag(), "before autocropfilter");
346 if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
348 trachea = clitk::AutoCrop<MaskImageType>(trachea, GetBackgroundValue());
351 // Remove Padding region
352 typedef itk::CropImageFilter<MaskImageType, MaskImageType> CropFilterType;
353 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
354 cropFilter->SetInput(trachea);
355 cropFilter->SetLowerBoundaryCropSize(bounds);
356 cropFilter->SetUpperBoundaryCropSize(bounds);
357 cropFilter->Update();
358 trachea = cropFilter->GetOutput();
360 StopCurrentStep<MaskImageType>(trachea);
361 PrintMemory(GetVerboseMemoryFlag(), "after delete trachea");
363 PrintMemory(GetVerboseMemoryFlag(), "after delete trachea");
365 //--------------------------------------------------------------------
366 //--------------------------------------------------------------------
367 StartNewStep("Cropping lung");
368 PrintMemory(GetVerboseMemoryFlag(), "Before Autocropfilter");
370 working_mask = clitk::AutoCrop<MaskImageType>(working_mask, GetBackgroundValue());
373 // Remove Padding region
374 typedef itk::CropImageFilter<MaskImageType, MaskImageType> CropFilterType;
375 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
376 cropFilter->SetInput(working_mask);
377 cropFilter->SetLowerBoundaryCropSize(bounds);
378 cropFilter->SetUpperBoundaryCropSize(bounds);
379 cropFilter->Update();
380 working_mask = cropFilter->GetOutput();
382 StopCurrentStep<MaskImageType>(working_mask);
384 //--------------------------------------------------------------------
385 //--------------------------------------------------------------------
387 if (GetOpenCloseFlag()) {
388 StartNewStep("Open/Close");
389 PrintMemory(GetVerboseMemoryFlag(), "Before OpenClose");
391 // Structuring element
392 typedef itk::BinaryBallStructuringElement<InternalPixelType, ImageDimension> KernelType;
393 KernelType structuringElement;
394 structuringElement.SetRadius(GetOpenCloseRadius());
395 structuringElement.CreateStructuringElement();
398 typedef itk::BinaryMorphologicalOpeningImageFilter<MaskImageType, InternalImageType, KernelType> OpenFilterType;
399 typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
400 openFilter->SetInput(working_mask);
401 openFilter->SetBackgroundValue(GetBackgroundValue());
402 openFilter->SetForegroundValue(GetForegroundValue());
403 openFilter->SetKernel(structuringElement);
406 typedef itk::BinaryMorphologicalClosingImageFilter<MaskImageType, MaskImageType, KernelType> CloseFilterType;
407 typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
408 closeFilter->SetInput(openFilter->GetOutput());
409 closeFilter->SetSafeBorder(true);
410 closeFilter->SetForegroundValue(GetForegroundValue());
411 closeFilter->SetKernel(structuringElement);
412 closeFilter->Update();
413 working_mask = closeFilter->GetOutput();
416 //--------------------------------------------------------------------
417 //--------------------------------------------------------------------
419 if (GetFillHolesFlag()) {
420 StartNewStep("Fill Holes");
421 PrintMemory(GetVerboseMemoryFlag(), "Before Fill Holes");
422 typedef clitk::FillMaskFilter<MaskImageType> FillMaskFilterType;
423 typename FillMaskFilterType::Pointer fillMaskFilter = FillMaskFilterType::New();
424 fillMaskFilter->SetInput(working_mask);
425 fillMaskFilter->AddDirection(2);
426 //fillMaskFilter->AddDirection(1);
427 fillMaskFilter->Update();
428 working_mask = fillMaskFilter->GetOutput();
429 StopCurrentStep<MaskImageType>(working_mask);
432 if (GetSeparateLungsFlag()) {
433 //--------------------------------------------------------------------
434 //--------------------------------------------------------------------
435 StartNewStep("Separate Left/Right lungs");
436 PrintMemory(GetVerboseMemoryFlag(), "Before Separate");
438 working_mask = Labelize<MaskImageType>(working_mask,
439 GetBackgroundValue(),
441 GetMinimalComponentSize());
443 PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
446 typedef itk::StatisticsImageFilter<MaskImageType> StatisticsImageFilterType;
447 typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
448 statisticsImageFilter->SetInput(working_mask);
449 statisticsImageFilter->Update();
450 unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
451 working_mask = statisticsImageFilter->GetOutput();
453 PrintMemory(GetVerboseMemoryFlag(), "After count label");
455 // Decompose the first label
456 if (initialNumberOfLabels<2) {
457 // Structuring element radius
458 typename ImageType::SizeType radius;
459 for (unsigned int i=0;i<Dim;i++) radius[i]=1;
460 typedef clitk::DecomposeAndReconstructImageFilter<MaskImageType,MaskImageType> DecomposeAndReconstructFilterType;
461 typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
462 decomposeAndReconstructFilter->SetInput(working_mask);
463 decomposeAndReconstructFilter->SetVerbose(false);
464 decomposeAndReconstructFilter->SetRadius(radius);
465 decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
466 decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
467 decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
468 decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
469 decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
470 decomposeAndReconstructFilter->SetFullyConnected(true);
471 decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
472 decomposeAndReconstructFilter->Update();
473 working_mask = decomposeAndReconstructFilter->GetOutput();
475 PrintMemory(GetVerboseMemoryFlag(), "After decomposeAndReconstructFilter");
478 // Retain labels ('1' is largset lung, so right. '2' is left)
479 typedef itk::ThresholdImageFilter<MaskImageType> ThresholdImageFilterType;
480 typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
481 thresholdFilter->SetInput(working_mask);
482 thresholdFilter->ThresholdAbove(2);
483 thresholdFilter->SetOutsideValue(this->GetBackgroundValue());
484 thresholdFilter->Update();
485 working_mask = thresholdFilter->GetOutput();
486 StopCurrentStep<MaskImageType> (working_mask);
487 PrintMemory(GetVerboseMemoryFlag(), "After Thresholdfilter");
489 // Update output info
490 // output = working_mask;
491 //this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
493 // this->GetOutput(0)->SetRegions(working_mask->GetLargestPossibleRegion());
496 //--------------------------------------------------------------------
499 //--------------------------------------------------------------------
500 template <class TImageType>
502 clitk::ExtractLungFilter<TImageType>::
503 GenerateInputRequestedRegion() {
504 // DD("GenerateInputRequestedRegion (nothing?)");
506 //--------------------------------------------------------------------
509 //--------------------------------------------------------------------
510 template <class ImageType>
512 clitk::ExtractLungFilter<ImageType>::
516 // this->GraftOutput(output); // not SetNthOutput
517 this->GraftOutput(working_mask); // not SetNthOutput
518 // Store image filenames into AFDB
519 GetAFDB()->SetImageFilename("Lungs", this->GetOutputLungFilename());
520 GetAFDB()->SetImageFilename("Trachea", this->GetOutputTracheaFilename());
523 //--------------------------------------------------------------------
526 //--------------------------------------------------------------------
527 template <class ImageType>
529 clitk::ExtractLungFilter<ImageType>::
530 SearchForTracheaSeed(int skip)
532 if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)
533 // Restart the search until a seed is found, skipping more and more slices
536 // Search seed (parameters = UpperThresholdForTrachea)
537 static const unsigned int Dim = ImageType::ImageDimension;
538 typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
539 typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
540 typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
541 sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-skip-5;
542 sliceRegionSize[Dim-1]=5;
543 sliceRegion.SetSize(sliceRegionSize);
544 sliceRegion.SetIndex(sliceRegionIndex);
545 typedef itk::ImageRegionConstIterator<ImageType> IteratorType;
546 IteratorType it(working_input, sliceRegion);
548 while (!it.IsAtEnd()) {
549 if(it.Get() < GetUpperThresholdForTrachea() ) {
550 AddSeed(it.GetIndex());
551 // DD(it.GetIndex());
556 // if we do not found : restart
557 stop = (m_Seeds.size() != 0);
559 if (GetVerboseStepFlag()) {
560 std::cout << "\t No seed found this time. I skip some slices and restart." << std::endl;
562 if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
563 // we want to skip more than a half of the image, it is probably a bug
564 std::cerr << "Number of slices to skip to find trachea to high = " << skip << std::endl;
571 // DD(m_Seeds.size());
575 return (m_Seeds.size() != 0);
577 //--------------------------------------------------------------------
580 bool is_orientation_superior(itk::SpatialOrientation::ValidCoordinateOrientationFlags orientation)
582 itk::SpatialOrientation::CoordinateTerms sup = itk::SpatialOrientation::ITK_COORDINATE_Superior;
583 std::cout << "orientation: " << std::hex << orientation << "; superior: " << std::hex << sup << std::endl;
584 std::cout << std::dec;
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 (GetVerboseFlag())
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 LabelImageType::LabelObjectContainerType shapes_map = label_map->GetLabelObjectContainer();
707 typename LabelImageType::LabelObjectContainerType::const_iterator s;
708 typename ShapeLabelType::Pointer shape, max_e_shape;
709 double max_elongation = GetMaxElongation();
710 double max_size = size[0]*size[1]/128;
713 for (s = shapes_map.begin(); s != shapes_map.end(); s++) {
715 if (shape->GetSize() > 150 && shape->GetSize() <= max_size) {
716 double e = 1 - 1/shape->GetBinaryElongation();
717 //double area = 1 - r->Area() ;
718 if (e < max_elongation) {
720 shape_centre[0] = (shape->GetCentroid()[0] - origin[0])/spacing[0];
721 shape_centre[1] = (shape->GetCentroid()[1] - origin[1])/spacing[1];
722 shape_centre[2] = index[2];
723 //double d = 1 - (shape_centre - image_centre).Magnitude()/max_dist;
724 double dx = shape_centre[0] - image_centre[0];
725 double d = 1 - dx*2/size[0];
731 max_e_centre = shape_centre;
739 itk::Point<typename SeedType::IndexValueType, ImageType::ImageDimension> p1, p2;
740 p1[0] = max_e_centre[0];
741 p1[1] = max_e_centre[1];
742 p1[2] = max_e_centre[2];
744 p2[0] = prev_e_centre[0];
745 p2[1] = prev_e_centre[1];
746 p2[2] = prev_e_centre[2];
748 double mag = (p2 - p1).GetNorm();
749 if (GetVerboseFlag()) {
751 cout << index[2] << ": ";
752 cout << "region(" << max_e_centre[0] << ", " << max_e_centre[1] << ", " << max_e_centre[2] << "); ";
753 cout << "prev_region(" << prev_e_centre[0] << ", " << prev_e_centre[1] << ", " << prev_e_centre[2] << "); ";
754 cout << "mag(" << mag << "); " << endl;
759 PointListType point_list;
760 point_list.push_back(max_e_centre);
761 sequence_list.push_back(point_list);
762 current_sequence = &sequence_list.back();
764 else if (current_sequence)
765 current_sequence->push_back(max_e_centre);
767 prev_e_centre= max_e_centre;
772 for (typename SequenceListType::iterator s = sequence_list.begin(); s != sequence_list.end(); s++)
774 if (s->size() > longest)
777 trachea_centre = s->front();
781 if (GetVerboseFlag())
782 std::cout << "seed at: " << trachea_centre << std::endl;
783 m_Seeds.push_back(trachea_centre);
786 return (m_Seeds.size() != 0);
788 //--------------------------------------------------------------------
790 //--------------------------------------------------------------------
791 template <class ImageType>
793 clitk::ExtractLungFilter<ImageType>::
794 TracheaRegionGrowing()
796 // Explosion controlled region growing
797 PrintMemory(GetVerboseMemoryFlag(), "Before ExplosionControlledThresholdConnectedImageFilter");
798 typedef clitk::ExplosionControlledThresholdConnectedImageFilter<ImageType, MaskImageType> ImageFilterType;
799 typename ImageFilterType::Pointer f= ImageFilterType::New();
800 f->SetInput(working_input);
802 f->SetUpper(GetUpperThresholdForTrachea());
803 f->SetMinimumLowerThreshold(-2000);
804 // f->SetMaximumUpperThreshold(0); // MAYBE TO CHANGE ???
805 f->SetMaximumUpperThreshold(-700); // MAYBE TO CHANGE ???
806 f->SetAdaptLowerBorder(false);
807 f->SetAdaptUpperBorder(true);
808 f->SetMinimumSize(5000);
809 f->SetReplaceValue(1);
810 f->SetMultiplier(GetMultiplierForTrachea());
811 f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
812 f->SetMinimumThresholdStepSize(1);
813 f->SetVerbose(GetVerboseRegionGrowingFlag());
814 for(unsigned int i=0; i<m_Seeds.size();i++) {
815 f->AddSeed(m_Seeds[i]);
819 PrintMemory(GetVerboseMemoryFlag(), "After RG update");
821 // take first (main) connected component
822 trachea = Labelize<MaskImageType>(f->GetOutput(),
823 GetBackgroundValue(),
825 1000);//GetMinimalComponentSize());
826 PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
827 trachea = KeepLabels<MaskImageType>(trachea,
828 GetBackgroundValue(),
829 GetForegroundValue(),
831 PrintMemory(GetVerboseMemoryFlag(), "After KeepLabels");
833 //--------------------------------------------------------------------
836 //--------------------------------------------------------------------
837 template <class ImageType>
839 clitk::ExtractLungFilter<ImageType>::
840 ComputeTracheaVolume()
842 typedef itk::ImageRegionConstIterator<InternalImageType> IteratorType;
843 IteratorType iter(trachea, trachea->GetLargestPossibleRegion());
846 while (!iter.IsAtEnd()) {
847 if (iter.Get() == this->GetForegroundValue()) volume++;
851 double voxelsize = trachea->GetSpacing()[0]*trachea->GetSpacing()[1]*trachea->GetSpacing()[2];
852 return volume*voxelsize;
854 //--------------------------------------------------------------------
857 //--------------------------------------------------------------------
858 template <class ImageType>
860 clitk::ExtractLungFilter<ImageType>::
863 // Search for seed among n slices, skip some slices before starting
864 // if not found -> skip more and restart
866 // when seed found : perform region growing
867 // compute trachea volume
868 // if volume not plausible -> skip more slices and restart
873 int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed();
876 if (GetTracheaSeedAlgorithm() == 0)
877 has_seed = SearchForTracheaSeed(skip);
879 has_seed = SearchForTracheaSeed2(GetNumSlices());
882 TracheaRegionGrowing();
883 volume = ComputeTracheaVolume()/1000; // assume mm3, so divide by 1000 to get cc
884 if (GetWriteStepFlag()) {
885 writeImage<MaskImageType>(trachea, "step-trachea-"+toString(skip)+".mhd");
887 if (GetTracheaVolumeMustBeCheckedFlag()) {
888 if ((volume > 10) && (volume < 65 )) { // depend on image size ...
889 // Typical volume 22.59 cm 3 (± 7.69 cm 3 ) [ Leader 2004 ]
890 if (GetVerboseStepFlag())
892 std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl;
896 if (GetVerboseStepFlag()) {
897 std::cout << "\t The volume of the trachea (" << volume
898 << " cc) seems not correct. I skip some slices (" << skip << ") and restart to find seeds."
903 // empty the list of seed
906 if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
907 // we want to skip more than a half of the image, it is probably a bug
908 std::cerr << "2 : Number of slices to skip to find trachea too high = " << skip << std::endl;
920 StopCurrentStep<MaskImageType>(trachea);
922 else { // Trachea not found
923 this->SetWarning("* WARNING * No seed found for trachea.");
927 //--------------------------------------------------------------------
929 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX