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"
50 #include "itkRelabelComponentImageFilter.h"
54 //--------------------------------------------------------------------
55 template <class ImageType>
56 clitk::ExtractLungFilter<ImageType>::
59 clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
60 itk::ImageToImageFilter<ImageType, MaskImageType>()
63 m_MaxSeedNumber = 500;
65 // Default global options
66 this->SetNumberOfRequiredInputs(1);
67 SetPatientMaskBackgroundValue(0);
68 SetBackgroundValue(0); // Must be zero
69 SetForegroundValue(1);
70 SetMinimalComponentSize(100);
71 VerboseRegionGrowingFlagOff();
72 RemoveSmallLabelBeforeSeparationFlagOn();
74 // Step 1 default values
75 SetUpperThreshold(-300);
76 SetLowerThreshold(-1000);
77 UseLowerThresholdOff();
78 LabelParamType * p1 = new LabelParamType;
81 p1->AddLabelToRemove(2);
82 SetLabelizeParameters1(p1);
84 // Step 2 default values
85 SetTracheaSeedAlgorithm(0);
86 SetUpperThresholdForTrachea(-700);
87 SetMultiplierForTrachea(5);
88 SetThresholdStepSizeForTrachea(64);
89 SetNumberOfSlicesToSkipBeforeSearchingSeed(0);
90 TracheaVolumeMustBeCheckedFlagOn();
92 SetMaxElongation(0.5);
93 SetSeedPreProcessingThreshold(-400);
96 // Step 3 default values
97 SetNumberOfHistogramBins(500);
98 LabelParamType * p2 = new LabelParamType;
100 p2->UseLastKeepOff();
101 // p->AddLabelToRemove(2); // No label to remove by default
102 SetLabelizeParameters2(p2);
104 // Step 4 default values
105 SetRadiusForTrachea(1);
106 LabelParamType * p3 = new LabelParamType;
109 p3->UseLastKeepOff();
110 SetLabelizeParameters3(p3);
114 SetOpenCloseRadius(1);
120 //--------------------------------------------------------------------
123 //--------------------------------------------------------------------
124 template <class ImageType>
126 clitk::ExtractLungFilter<ImageType>::
127 SetInput(const ImageType * image)
129 this->SetNthInput(0, const_cast<ImageType *>(image));
131 //--------------------------------------------------------------------
134 //--------------------------------------------------------------------
135 template <class ImageType>
137 clitk::ExtractLungFilter<ImageType>::
138 AddSeed(InternalIndexType s)
140 m_Seeds.push_back(s);
142 //--------------------------------------------------------------------
145 //--------------------------------------------------------------------
146 template <class ImageType>
148 clitk::ExtractLungFilter<ImageType>::
149 GenerateOutputInformation()
151 clitk::PrintMemory(GetVerboseMemoryFlag(), "Initial memory"); // OK
152 Superclass::GenerateOutputInformation();
157 // Get input pointers
158 input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
159 patient = GetAFDB()->template GetImage <MaskImageType>("Patient");
160 PrintMemory(GetVerboseMemoryFlag(), "After reading patient"); // OK
162 //--------------------------------------------------------------------
163 //--------------------------------------------------------------------
164 // Crop input like patient image (must have the same spacing)
165 // It copy the input if the same are the same
166 StartNewStep("Copy and crop input image to 'patient' extends");
167 typedef clitk::CropLikeImageFilter<ImageType> CropImageFilter;
168 typename CropImageFilter::Pointer cropFilter = CropImageFilter::New();
169 // cropFilter->ReleaseDataFlagOn(); // NO=seg fault !!
170 cropFilter->SetInput(input);
171 cropFilter->SetCropLikeImage(patient);
172 cropFilter->Update();
173 working_input = cropFilter->GetOutput();
174 // cropFilter->Delete(); // NO !!!! if yes, sg fault, Cropfilter is buggy !??
175 StopCurrentStep<ImageType>(working_input);
176 PrintMemory(GetVerboseMemoryFlag(), "After crop"); // OK, slightly more than a copy of input
178 //--------------------------------------------------------------------
179 //--------------------------------------------------------------------
180 StartNewStep("Set background to initial image");
181 working_input = SetBackground<ImageType, MaskImageType>
182 (working_input, patient, GetPatientMaskBackgroundValue(), -1000, true);
184 // Pad images with air to prevent patient touching the image border
185 static const unsigned int Dim = ImageType::ImageDimension;
186 typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadFilterType;
187 typename PadFilterType::Pointer padFilter = PadFilterType::New();
188 padFilter->SetInput(working_input);
189 padFilter->SetConstant(-1000);
190 typename ImageType::SizeType bounds;
191 for (unsigned i = 0; i < Dim - 1; ++i)
194 padFilter->SetPadLowerBound(bounds);
195 padFilter->SetPadUpperBound(bounds);
197 working_input = padFilter->GetOutput();
199 StopCurrentStep<ImageType>(working_input);
200 PrintMemory(GetVerboseMemoryFlag(), "After set bg"); // OK, additional mem = 0
203 // We do not need patient mask anymore, release memory
204 GetAFDB()->template ReleaseImage<MaskImageType>("Patient");
205 DD(patient->GetReferenceCount());
206 PrintMemory(GetVerboseMemoryFlag(), "After delete patient"); // OK, additional mem = 0
208 PrintMemory(GetVerboseMemoryFlag(), "After delete patient"); // OK, additional mem = 0
211 //--------------------------------------------------------------------
212 //--------------------------------------------------------------------
213 StartNewStep("Remove Air");
215 if (m_UseLowerThreshold) {
216 if (m_LowerThreshold > m_UpperThreshold) {
217 clitkExceptionMacro("lower threshold cannot be greater than upper threshold.");
220 // Threshold to get air
221 typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> InputBinarizeFilterType;
222 typename InputBinarizeFilterType::Pointer binarizeFilter = InputBinarizeFilterType::New();
223 binarizeFilter->SetInput(working_input);
224 if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(m_LowerThreshold);
225 // binarizeFilter->CanRunInPlace() is false
226 binarizeFilter->SetUpperThreshold(m_UpperThreshold);
227 binarizeFilter->SetInsideValue(this->GetForegroundValue());
228 binarizeFilter->SetOutsideValue(this->GetBackgroundValue());
229 binarizeFilter->Update();
230 working_mask = binarizeFilter->GetOutput();
231 PrintMemory(GetVerboseMemoryFlag(), "After Binarizefilter"); // OK, additional mem is one mask image
233 // Labelize and keep right labels
235 Labelize<MaskImageType>
236 (working_mask, GetBackgroundValue(), true, GetMinimalComponentSize());
237 PrintMemory(GetVerboseMemoryFlag(), "After Labelize"); // BUG ? additional mem around 1 time the input ?
239 working_mask = RemoveLabels<MaskImageType>
240 (working_mask, GetBackgroundValue(), GetLabelizeParameters1()->GetLabelsToRemove());
241 PrintMemory(GetVerboseMemoryFlag(), "After RemoveLabels"); // OK additional mem = 0
243 working_mask = KeepLabels<MaskImageType>
245 GetBackgroundValue(),
246 GetForegroundValue(),
247 GetLabelizeParameters1()->GetFirstKeep(),
248 GetLabelizeParameters1()->GetLastKeep(),
249 GetLabelizeParameters1()->GetUseLastKeep());
250 PrintMemory(GetVerboseMemoryFlag(), "After KeepLabels to create the 'air'");
253 working_input = SetBackground<ImageType, MaskImageType>
254 (working_input, working_mask, this->GetForegroundValue(), this->GetBackgroundValue(), true);
255 PrintMemory(GetVerboseMemoryFlag(), "After SetBackground");
258 StopCurrentStep<ImageType>(working_input);
260 //--------------------------------------------------------------------
261 //--------------------------------------------------------------------
262 StartNewStep("Search for the trachea");
264 PrintMemory(GetVerboseMemoryFlag(), "After SearchForTrachea");
265 if (m_Seeds.empty()) {
266 clitkExceptionMacro("No seeds for trachea... Aborting.");
269 //--------------------------------------------------------------------
270 //--------------------------------------------------------------------
271 StartNewStep("Extract the lung with Otsu filter");
272 // Automated Otsu thresholding and relabeling
273 typedef itk::OtsuThresholdImageFilter<ImageType,MaskImageType> OtsuThresholdImageFilterType;
274 typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New();
275 otsuFilter->SetInput(working_input);
276 otsuFilter->SetNumberOfHistogramBins(GetNumberOfHistogramBins());
277 otsuFilter->SetInsideValue(this->GetForegroundValue());
278 otsuFilter->SetOutsideValue(this->GetBackgroundValue());
279 otsuFilter->Update();
280 working_mask = otsuFilter->GetOutput();
283 StopCurrentStep<MaskImageType>(working_mask);
284 PrintMemory(GetVerboseMemoryFlag(), "After Otsufilter");
286 //--------------------------------------------------------------------
287 //--------------------------------------------------------------------
288 StartNewStep("Select labels");
290 working_mask = LabelizeAndSelectLabels<MaskImageType>
292 GetBackgroundValue(),
293 GetForegroundValue(),
295 GetMinimalComponentSize(),
296 GetLabelizeParameters2());
299 StopCurrentStep<MaskImageType>(working_mask);
300 PrintMemory(GetVerboseMemoryFlag(), "After LabelizeAndSelectLabels");
302 //--------------------------------------------------------------------
303 //--------------------------------------------------------------------
304 if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
305 StartNewStep("Remove the trachea");
307 working_mask = SetBackground<MaskImageType, MaskImageType>
308 (working_mask, trachea, 1, -1, true);
309 PrintMemory(GetVerboseMemoryFlag(), "After SetBackground");
311 // Dilate the trachea
312 static const unsigned int Dim = ImageType::ImageDimension;
313 typedef itk::BinaryBallStructuringElement<InternalPixelType, Dim> KernelType;
314 KernelType structuringElement;
315 structuringElement.SetRadius(GetRadiusForTrachea());
316 structuringElement.CreateStructuringElement();
317 typedef clitk::ConditionalBinaryDilateImageFilter<MaskImageType, MaskImageType, KernelType> ConditionalBinaryDilateImageFilterType;
318 typename ConditionalBinaryDilateImageFilterType::Pointer dilateFilter = ConditionalBinaryDilateImageFilterType::New();
319 dilateFilter->SetBoundaryToForeground(false);
320 dilateFilter->SetKernel(structuringElement);
321 dilateFilter->SetBackgroundValue (1);
322 dilateFilter->SetForegroundValue (-1);
323 dilateFilter->SetInput (working_mask);
324 dilateFilter->Update();
325 working_mask = dilateFilter->GetOutput();
326 PrintMemory(GetVerboseMemoryFlag(), "After dilate");
328 // Set trachea with dilatation
329 trachea = SetBackground<MaskImageType, MaskImageType>
330 (trachea, working_mask, -1, this->GetForegroundValue(), true);
332 // Remove the trachea
333 working_mask = SetBackground<MaskImageType, MaskImageType>
334 (working_mask, working_mask, -1, this->GetBackgroundValue(), true);
337 working_mask = LabelizeAndSelectLabels<MaskImageType>
339 GetBackgroundValue(),
340 GetForegroundValue(),
342 GetMinimalComponentSize(),
343 GetLabelizeParameters3());
346 StopCurrentStep<MaskImageType>(working_mask);
349 //--------------------------------------------------------------------
350 //--------------------------------------------------------------------
351 PrintMemory(GetVerboseMemoryFlag(), "before autocropfilter");
352 if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
354 trachea = clitk::AutoCrop<MaskImageType>(trachea, GetBackgroundValue());
357 // Remove Padding region
358 typedef itk::CropImageFilter<MaskImageType, MaskImageType> CropFilterType;
359 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
360 cropFilter->SetInput(trachea);
361 cropFilter->SetLowerBoundaryCropSize(bounds);
362 cropFilter->SetUpperBoundaryCropSize(bounds);
363 cropFilter->Update();
364 trachea = cropFilter->GetOutput();
366 StopCurrentStep<MaskImageType>(trachea);
367 PrintMemory(GetVerboseMemoryFlag(), "after delete trachea");
369 PrintMemory(GetVerboseMemoryFlag(), "after delete trachea");
371 //--------------------------------------------------------------------
372 //--------------------------------------------------------------------
373 StartNewStep("Cropping lung");
374 PrintMemory(GetVerboseMemoryFlag(), "Before Autocropfilter");
376 working_mask = clitk::AutoCrop<MaskImageType>(working_mask, GetBackgroundValue());
379 // Remove Padding region
380 typedef itk::CropImageFilter<MaskImageType, MaskImageType> CropFilterType;
381 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
382 cropFilter->SetInput(working_mask);
383 cropFilter->SetLowerBoundaryCropSize(bounds);
384 cropFilter->SetUpperBoundaryCropSize(bounds);
385 cropFilter->Update();
386 working_mask = cropFilter->GetOutput();
388 StopCurrentStep<MaskImageType>(working_mask);
390 //--------------------------------------------------------------------
391 //--------------------------------------------------------------------
393 if (GetOpenCloseFlag()) {
394 StartNewStep("Open/Close");
395 PrintMemory(GetVerboseMemoryFlag(), "Before OpenClose");
397 // Structuring element
398 typedef itk::BinaryBallStructuringElement<InternalPixelType, ImageDimension> KernelType;
399 KernelType structuringElement;
400 structuringElement.SetRadius(GetOpenCloseRadius());
401 structuringElement.CreateStructuringElement();
404 typedef itk::BinaryMorphologicalOpeningImageFilter<MaskImageType, InternalImageType, KernelType> OpenFilterType;
405 typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
406 openFilter->SetInput(working_mask);
407 openFilter->SetBackgroundValue(GetBackgroundValue());
408 openFilter->SetForegroundValue(GetForegroundValue());
409 openFilter->SetKernel(structuringElement);
412 typedef itk::BinaryMorphologicalClosingImageFilter<MaskImageType, MaskImageType, KernelType> CloseFilterType;
413 typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
414 closeFilter->SetInput(openFilter->GetOutput());
415 closeFilter->SetSafeBorder(true);
416 closeFilter->SetForegroundValue(GetForegroundValue());
417 closeFilter->SetKernel(structuringElement);
418 closeFilter->Update();
419 working_mask = closeFilter->GetOutput();
422 //--------------------------------------------------------------------
423 //--------------------------------------------------------------------
425 if (GetFillHolesFlag()) {
426 StartNewStep("Fill Holes");
427 PrintMemory(GetVerboseMemoryFlag(), "Before Fill Holes");
428 typedef clitk::FillMaskFilter<MaskImageType> FillMaskFilterType;
429 typename FillMaskFilterType::Pointer fillMaskFilter = FillMaskFilterType::New();
430 fillMaskFilter->SetInput(working_mask);
431 fillMaskFilter->AddDirection(2);
432 //fillMaskFilter->AddDirection(1);
433 fillMaskFilter->Update();
434 working_mask = fillMaskFilter->GetOutput();
435 StopCurrentStep<MaskImageType>(working_mask);
438 if (GetSeparateLungsFlag()) {
439 //--------------------------------------------------------------------
440 //--------------------------------------------------------------------
441 StartNewStep("Separate Left/Right lungs");
442 PrintMemory(GetVerboseMemoryFlag(), "Before Separate");
444 working_mask = Labelize<MaskImageType>(working_mask,
445 GetBackgroundValue(),
447 GetMinimalComponentSize());
449 PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
452 typedef itk::StatisticsImageFilter<MaskImageType> StatisticsImageFilterType;
453 typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
454 statisticsImageFilter->SetInput(working_mask);
455 statisticsImageFilter->Update();
456 unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
457 working_mask = statisticsImageFilter->GetOutput();
458 PrintMemory(GetVerboseMemoryFlag(), "After count label");
460 // If already 2 labels, but a too big differences, remove the
461 // smalest one (sometimes appends with the stomach
462 if (initialNumberOfLabels >1) {
463 if (GetRemoveSmallLabelBeforeSeparationFlag()) {
464 DD(GetRemoveSmallLabelBeforeSeparationFlag());
465 typedef itk::RelabelComponentImageFilter<MaskImageType, MaskImageType> RelabelFilterType;
466 typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
467 relabelFilter->SetInput(working_mask);
468 relabelFilter->SetMinimumObjectSize(10);
469 relabelFilter->Update();
470 const std::vector<float> & a = relabelFilter->GetSizeOfObjectsInPhysicalUnits();
471 std::vector<MaskImagePixelType> remove_label;
472 for(unsigned int i=1; i<a.size(); i++) {
473 if (a[i] < 0.5*a[0]) { // more than 0.5 difference
474 remove_label.push_back(i+1); // label zero is BG
478 clitk::RemoveLabels<MaskImageType>(working_mask, GetBackgroundValue(), remove_label);
479 statisticsImageFilter->SetInput(working_mask);
480 statisticsImageFilter->Update();
481 initialNumberOfLabels = statisticsImageFilter->GetMaximum();
485 // Decompose the first label
486 if (initialNumberOfLabels<2) {
487 // Structuring element radius
488 typename ImageType::SizeType radius;
489 for (unsigned int i=0;i<Dim;i++) radius[i]=1;
490 typedef clitk::DecomposeAndReconstructImageFilter<MaskImageType,MaskImageType> DecomposeAndReconstructFilterType;
491 typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
492 decomposeAndReconstructFilter->SetInput(working_mask);
493 decomposeAndReconstructFilter->SetVerbose(false);
494 decomposeAndReconstructFilter->SetRadius(radius);
495 decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
496 decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
497 decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
498 decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
499 decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
500 decomposeAndReconstructFilter->SetFullyConnected(true);
501 decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
502 decomposeAndReconstructFilter->Update();
503 working_mask = decomposeAndReconstructFilter->GetOutput();
505 PrintMemory(GetVerboseMemoryFlag(), "After decomposeAndReconstructFilter");
508 // Retain labels ('1' is largset lung, so right. '2' is left)
509 typedef itk::ThresholdImageFilter<MaskImageType> ThresholdImageFilterType;
510 typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
511 thresholdFilter->SetInput(working_mask);
512 thresholdFilter->ThresholdAbove(2);
513 thresholdFilter->SetOutsideValue(this->GetBackgroundValue());
514 thresholdFilter->Update();
515 working_mask = thresholdFilter->GetOutput();
516 StopCurrentStep<MaskImageType> (working_mask);
517 PrintMemory(GetVerboseMemoryFlag(), "After Thresholdfilter");
519 // Update output info
520 // output = working_mask;
521 //this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
523 // this->GetOutput(0)->SetRegions(working_mask->GetLargestPossibleRegion());
526 //--------------------------------------------------------------------
529 //--------------------------------------------------------------------
530 template <class TImageType>
532 clitk::ExtractLungFilter<TImageType>::
533 GenerateInputRequestedRegion() {
534 // DD("GenerateInputRequestedRegion (nothing?)");
536 //--------------------------------------------------------------------
539 //--------------------------------------------------------------------
540 template <class ImageType>
542 clitk::ExtractLungFilter<ImageType>::
546 // this->GraftOutput(output); // not SetNthOutput
547 this->GraftOutput(working_mask); // not SetNthOutput
548 // Store image filenames into AFDB
549 GetAFDB()->SetImageFilename("Lungs", this->GetOutputLungFilename());
550 GetAFDB()->SetImageFilename("Trachea", this->GetOutputTracheaFilename());
553 //--------------------------------------------------------------------
556 //--------------------------------------------------------------------
557 template <class ImageType>
559 clitk::ExtractLungFilter<ImageType>::
560 SearchForTracheaSeed(int skip)
562 if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)
563 // Restart the search until a seed is found, skipping more and more slices
566 // Search seed (parameters = UpperThresholdForTrachea)
567 static const unsigned int Dim = ImageType::ImageDimension;
568 typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
569 typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
570 typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
571 sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-skip-5;
572 sliceRegionSize[Dim-1]=5;
573 sliceRegion.SetSize(sliceRegionSize);
574 sliceRegion.SetIndex(sliceRegionIndex);
575 typedef itk::ImageRegionConstIterator<ImageType> IteratorType;
576 IteratorType it(working_input, sliceRegion);
578 while (!it.IsAtEnd()) {
579 if(it.Get() < GetUpperThresholdForTrachea() ) {
580 AddSeed(it.GetIndex());
581 // DD(it.GetIndex());
586 // if we do not found : restart
587 stop = (m_Seeds.size() != 0);
589 if (GetVerboseStepFlag()) {
590 std::cout << "\t No seed found this time. I skip some slices and restart." << std::endl;
592 if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
593 // we want to skip more than a half of the image, it is probably a bug
594 std::cerr << "Number of slices to skip to find trachea to high = " << skip << std::endl;
601 // DD(m_Seeds.size());
605 return (m_Seeds.size() != 0);
607 //--------------------------------------------------------------------
610 bool is_orientation_superior(itk::SpatialOrientation::ValidCoordinateOrientationFlags orientation)
612 itk::SpatialOrientation::CoordinateTerms sup = itk::SpatialOrientation::ITK_COORDINATE_Superior;
613 bool primary = (orientation & 0x0000ff) == sup;
614 bool secondary = ((orientation & 0x00ff00) >> 8) == sup;
615 bool tertiary = ((orientation & 0xff0000) >> 16) == sup;
616 return primary || secondary || tertiary;
619 //--------------------------------------------------------------------
620 template <class ImageType>
622 clitk::ExtractLungFilter<ImageType>::
623 SearchForTracheaSeed2(int numberOfSlices)
625 if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)
626 if (GetVerboseRegionGrowingFlag())
627 std::cout << "SearchForTracheaSeed2(" << numberOfSlices << ", " << GetMaxElongation() << ")" << std::endl;
629 typedef unsigned char MaskPixelType;
630 typedef itk::Image<MaskPixelType, ImageType::ImageDimension> MaskImageType;
631 typedef itk::Image<typename MaskImageType::PixelType, 2> MaskImageType2D;
632 typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> ThresholdFilterType;
633 typedef itk::BinaryBallStructuringElement<MaskPixelType, MaskImageType::ImageDimension> KernelType;
634 typedef itk::BinaryMorphologicalClosingImageFilter<MaskImageType, MaskImageType, KernelType> ClosingFilterType;
635 typedef itk::BinaryMorphologicalOpeningImageFilter<MaskImageType, MaskImageType, KernelType> OpeningFilterType;
636 typedef itk::ExtractImageFilter<MaskImageType, MaskImageType2D> SlicerFilterType;
637 typedef itk::ConnectedComponentImageFilter<MaskImageType2D, MaskImageType2D> LabelFilterType;
638 typedef itk::ShapeLabelObject<MaskPixelType, MaskImageType2D::ImageDimension> ShapeLabelType;
639 typedef itk::LabelMap<ShapeLabelType> LabelImageType;
640 typedef itk::LabelImageToShapeLabelMapFilter<MaskImageType2D, LabelImageType> ImageToLabelMapFilterType;
641 typedef itk::LabelMapToLabelImageFilter<LabelImageType, MaskImageType2D> LabelMapToImageFilterType;
642 typedef itk::ImageFileWriter<MaskImageType2D> FileWriterType;
644 // threshold to isolate airawys and lungs
645 typename ThresholdFilterType::Pointer threshold = ThresholdFilterType::New();
646 threshold->SetLowerThreshold(-2000);
647 threshold->SetUpperThreshold(GetSeedPreProcessingThreshold());
648 threshold->SetInput(working_input);
651 KernelType kernel_closing, kernel_opening;
653 // remove small noise
654 typename OpeningFilterType::Pointer opening = OpeningFilterType::New();
655 kernel_opening.SetRadius(1);
656 opening->SetKernel(kernel_opening);
657 opening->SetInput(threshold->GetOutput());
660 typename SlicerFilterType::Pointer slicer = SlicerFilterType::New();
661 slicer->SetInput(opening->GetOutput());
664 typename LabelFilterType::Pointer label_filter = LabelFilterType::New();
665 label_filter->SetInput(slicer->GetOutput());
667 // extract shape information from labels
668 typename ImageToLabelMapFilterType::Pointer label_to_map_filter = ImageToLabelMapFilterType::New();
669 label_to_map_filter->SetInput(label_filter->GetOutput());
671 typename LabelMapToImageFilterType::Pointer map_to_label_filter = LabelMapToImageFilterType::New();
672 typename FileWriterType::Pointer writer = FileWriterType::New();
674 typename ImageType::IndexType index;
675 typename ImageType::RegionType region = working_input->GetLargestPossibleRegion();
676 typename ImageType::SizeType size = region.GetSize();
677 typename ImageType::SpacingType spacing = working_input->GetSpacing();
678 typename ImageType::PointType origin = working_input->GetOrigin();
680 int nslices = min(numberOfSlices, size[2]);
681 int start = 0, increment = 1;
682 itk::SpatialOrientationAdapter orientation;
683 typename ImageType::DirectionType dir = working_input->GetDirection();
684 if (!is_orientation_superior(orientation.FromDirectionCosines(dir))) {
689 typename MaskImageType::PointType image_centre;
690 image_centre[0] = size[0]/2;
691 image_centre[1] = size[1]/2;
694 typedef InternalIndexType SeedType;
695 SeedType trachea_centre, shape_centre, max_e_centre, prev_e_centre;
696 typedef std::list<SeedType> PointListType;
697 typedef std::list<PointListType> SequenceListType;
698 PointListType* current_sequence = NULL;
699 SequenceListType sequence_list;
701 prev_e_centre.Fill(0);
702 std::ostringstream file_name;
703 index[0] = index[1] = 0;
704 size[0] = size[1] = 512;
710 region.SetIndex(index);
711 region.SetSize(size);
712 slicer->SetExtractionRegion(region);
714 label_filter->SetInput(slicer->GetOutput());
715 label_filter->Update();
717 label_to_map_filter->SetInput(label_filter->GetOutput());
718 label_to_map_filter->Update();
719 typename LabelImageType::Pointer label_map = label_to_map_filter->GetOutput();
721 if (GetWriteStepFlag()) {
722 map_to_label_filter->SetInput(label_map);
723 writer->SetInput(map_to_label_filter->GetOutput());
725 file_name << "labels_";
728 file_name << index[2] << ".mhd";
729 writer->SetFileName(file_name.str().c_str());
733 typename ShapeLabelType::Pointer shape, max_e_shape;
734 double max_elongation = GetMaxElongation();
735 double max_size = size[0]*size[1]/128;
738 unsigned int nlables = label_map->GetNumberOfLabelObjects();
739 for (unsigned int j = 0; j < nlables; j++) {
740 shape = label_map->GetNthLabelObject(j);
741 if (shape->Size() > 150 && shape->Size() <= max_size) {
742 #if ITK_VERSION_MAJOR < 4
743 double e = 1 - 1/shape->GetBinaryElongation();
745 double e = 1 - 1/shape->GetElongation();
747 //double area = 1 - r->Area() ;
748 if (e < max_elongation) {
750 shape_centre[0] = (shape->GetCentroid()[0] - origin[0])/spacing[0];
751 shape_centre[1] = (shape->GetCentroid()[1] - origin[1])/spacing[1];
752 shape_centre[2] = index[2];
753 //double d = 1 - (shape_centre - image_centre).Magnitude()/max_dist;
754 double dx = shape_centre[0] - image_centre[0];
755 double d = 1 - dx*2/size[0];
761 max_e_centre = shape_centre;
769 itk::Point<typename SeedType::IndexValueType, ImageType::ImageDimension> p1, p2;
770 p1[0] = max_e_centre[0];
771 p1[1] = max_e_centre[1];
772 p1[2] = max_e_centre[2];
774 p2[0] = prev_e_centre[0];
775 p2[1] = prev_e_centre[1];
776 p2[2] = prev_e_centre[2];
778 double mag = (p2 - p1).GetNorm();
779 if (GetVerboseRegionGrowingFlag()) {
781 cout << index[2] << ": ";
782 cout << "region(" << max_e_centre[0] << ", " << max_e_centre[1] << ", " << max_e_centre[2] << "); ";
783 cout << "prev_region(" << prev_e_centre[0] << ", " << prev_e_centre[1] << ", " << prev_e_centre[2] << "); ";
784 cout << "mag(" << mag << "); " << endl;
789 PointListType point_list;
790 point_list.push_back(max_e_centre);
791 sequence_list.push_back(point_list);
792 current_sequence = &sequence_list.back();
794 else if (current_sequence)
795 current_sequence->push_back(max_e_centre);
797 prev_e_centre= max_e_centre;
800 if (GetVerboseRegionGrowingFlag()) {
801 cout << "No shapes found at slice " << index[2] << std::endl;
807 for (typename SequenceListType::iterator s = sequence_list.begin(); s != sequence_list.end(); s++)
809 if (s->size() > longest)
812 trachea_centre = s->front();
817 if (GetVerboseRegionGrowingFlag())
818 std::cout << "seed at: " << trachea_centre << std::endl;
819 m_Seeds.push_back(trachea_centre);
823 return (m_Seeds.size() != 0);
825 //--------------------------------------------------------------------
827 //--------------------------------------------------------------------
828 template <class ImageType>
830 clitk::ExtractLungFilter<ImageType>::
831 TracheaRegionGrowing()
833 // Explosion controlled region growing
834 PrintMemory(GetVerboseMemoryFlag(), "Before ExplosionControlledThresholdConnectedImageFilter");
835 typedef clitk::ExplosionControlledThresholdConnectedImageFilter<ImageType, MaskImageType> ImageFilterType;
836 typename ImageFilterType::Pointer f= ImageFilterType::New();
837 f->SetInput(working_input);
839 f->SetUpper(GetUpperThresholdForTrachea());
840 f->SetMinimumLowerThreshold(-2000);
841 // f->SetMaximumUpperThreshold(0); // MAYBE TO CHANGE ???
842 f->SetMaximumUpperThreshold(-300); // MAYBE TO CHANGE ???
843 f->SetAdaptLowerBorder(false);
844 f->SetAdaptUpperBorder(true);
845 f->SetMinimumSize(5000);
846 f->SetReplaceValue(1);
847 f->SetMultiplier(GetMultiplierForTrachea());
848 f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
849 f->SetMinimumThresholdStepSize(1);
850 f->SetVerbose(GetVerboseRegionGrowingFlag());
851 for(unsigned int i=0; i<m_Seeds.size();i++) {
852 f->AddSeed(m_Seeds[i]);
856 PrintMemory(GetVerboseMemoryFlag(), "After RG update");
858 // take first (main) connected component
859 trachea = Labelize<MaskImageType>(f->GetOutput(),
860 GetBackgroundValue(),
862 1000);//GetMinimalComponentSize());
863 PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
864 trachea = KeepLabels<MaskImageType>(trachea,
865 GetBackgroundValue(),
866 GetForegroundValue(),
868 PrintMemory(GetVerboseMemoryFlag(), "After KeepLabels");
870 //--------------------------------------------------------------------
873 //--------------------------------------------------------------------
874 template <class ImageType>
876 clitk::ExtractLungFilter<ImageType>::
877 ComputeTracheaVolume()
879 typedef itk::ImageRegionConstIterator<InternalImageType> IteratorType;
880 IteratorType iter(trachea, trachea->GetLargestPossibleRegion());
883 while (!iter.IsAtEnd()) {
884 if (iter.Get() == this->GetForegroundValue()) volume++;
888 double voxelsize = trachea->GetSpacing()[0]*trachea->GetSpacing()[1]*trachea->GetSpacing()[2];
889 return volume*voxelsize;
891 //--------------------------------------------------------------------
894 //--------------------------------------------------------------------
895 template <class ImageType>
897 clitk::ExtractLungFilter<ImageType>::
900 // Search for seed among n slices, skip some slices before starting
901 // if not found -> skip more and restart
903 // when seed found : perform region growing
904 // compute trachea volume
905 // if volume not plausible -> skip more slices and restart
910 int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed();
913 if (GetTracheaSeedAlgorithm() == 0)
914 has_seed = SearchForTracheaSeed(skip);
916 has_seed = SearchForTracheaSeed2(GetNumSlices());
919 TracheaRegionGrowing();
920 volume = ComputeTracheaVolume()/1000; // assume mm3, so divide by 1000 to get cc
921 if (GetWriteStepFlag()) {
922 writeImage<MaskImageType>(trachea, "step-trachea-"+toString(skip)+".mhd");
924 if (GetTracheaVolumeMustBeCheckedFlag()) {
925 if ((volume > 10) && (volume < 65 )) { // depend on image size ...
926 // Typical volume 22.59 cm 3 (± 7.69 cm 3 ) [ Leader 2004 ]
927 if (GetVerboseStepFlag())
929 std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl;
933 if (GetTracheaSeedAlgorithm() == 0) {
934 if (GetVerboseStepFlag()) {
935 std::cout << "\t The volume of the trachea (" << volume
936 << " cc) seems not correct. I skip some slices (" << skip << ") and restart to find seeds."
941 // empty the list of seed
944 if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
945 // we want to skip more than a half of the image, it is probably a bug
946 std::cerr << "2 : Number of slices to skip to find trachea too high = " << skip << std::endl;
958 StopCurrentStep<MaskImageType>(trachea);
960 else { // Trachea not found
961 this->SetWarning("* WARNING * No seed found for trachea.");
965 //--------------------------------------------------------------------
967 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX