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);
95 // Step 3 default values
96 SetNumberOfHistogramBins(500);
97 LabelParamType * p2 = new LabelParamType;
100 // p->AddLabelToRemove(2); // No label to remove by default
101 SetLabelizeParameters2(p2);
103 // Step 4 default values
104 SetRadiusForTrachea(1);
105 LabelParamType * p3 = new LabelParamType;
108 p3->UseLastKeepOff();
109 SetLabelizeParameters3(p3);
113 SetOpenCloseRadius(1);
119 //--------------------------------------------------------------------
122 //--------------------------------------------------------------------
123 template <class ImageType>
125 clitk::ExtractLungFilter<ImageType>::
126 SetInput(const ImageType * image)
128 this->SetNthInput(0, const_cast<ImageType *>(image));
130 //--------------------------------------------------------------------
133 //--------------------------------------------------------------------
134 template <class ImageType>
136 clitk::ExtractLungFilter<ImageType>::
137 //AddSeed(InternalIndexType s)
138 AddSeed(InputImagePointType s)
140 m_SeedsInMM.push_back(s);
142 //--------------------------------------------------------------------
145 //--------------------------------------------------------------------
146 template <class ImageType>
148 clitk::ExtractLungFilter<ImageType>::
149 AddSeedInPixels(InternalIndexType s)
151 m_Seeds.push_back(s);
153 //--------------------------------------------------------------------
156 //--------------------------------------------------------------------
157 template <class ImageType>
159 clitk::ExtractLungFilter<ImageType>::
160 GenerateOutputInformation()
162 clitk::PrintMemory(GetVerboseMemoryFlag(), "Initial memory"); // OK
163 Superclass::GenerateOutputInformation();
168 // Get input pointers
169 input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
170 patient = GetAFDB()->template GetImage <MaskImageType>("Patient");
171 PrintMemory(GetVerboseMemoryFlag(), "After reading patient"); // OK
173 //--------------------------------------------------------------------
174 //--------------------------------------------------------------------
175 // Crop input like patient image (must have the same spacing)
176 // It copy the input if the same are the same
177 StartNewStep("Copy and crop input image to 'patient' extends");
178 typedef clitk::CropLikeImageFilter<ImageType> CropImageFilter;
179 typename CropImageFilter::Pointer cropFilter = CropImageFilter::New();
180 // cropFilter->ReleaseDataFlagOn(); // NO=seg fault !!
181 cropFilter->SetInput(input);
182 cropFilter->SetCropLikeImage(patient);
183 cropFilter->Update();
184 working_input = cropFilter->GetOutput();
185 // cropFilter->Delete(); // NO !!!! if yes, sg fault, Cropfilter is buggy !??
186 StopCurrentStep<ImageType>(working_input);
187 PrintMemory(GetVerboseMemoryFlag(), "After crop"); // OK, slightly more than a copy of input
189 //--------------------------------------------------------------------
190 //--------------------------------------------------------------------
191 StartNewStep("Set background to initial image");
192 working_input = SetBackground<ImageType, MaskImageType>
193 (working_input, patient, GetPatientMaskBackgroundValue(), -1000, true);
195 // Pad images with air to prevent patient touching the image border
196 static const unsigned int Dim = ImageType::ImageDimension;
197 typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadFilterType;
198 typename PadFilterType::Pointer padFilter = PadFilterType::New();
199 padFilter->SetInput(working_input);
200 padFilter->SetConstant(-1000);
201 typename ImageType::SizeType bounds;
202 for (unsigned i = 0; i < Dim - 1; ++i)
205 padFilter->SetPadLowerBound(bounds);
206 padFilter->SetPadUpperBound(bounds);
208 working_input = padFilter->GetOutput();
210 StopCurrentStep<ImageType>(working_input);
211 PrintMemory(GetVerboseMemoryFlag(), "After set bg"); // OK, additional mem = 0
214 // We do not need patient mask anymore, release memory
215 GetAFDB()->template ReleaseImage<MaskImageType>("Patient");
216 DD(patient->GetReferenceCount());
217 PrintMemory(GetVerboseMemoryFlag(), "After delete patient"); // OK, additional mem = 0
219 PrintMemory(GetVerboseMemoryFlag(), "After delete patient"); // OK, additional mem = 0
222 //--------------------------------------------------------------------
223 //--------------------------------------------------------------------
224 StartNewStep("Remove Air");
226 if (m_UseLowerThreshold) {
227 if (m_LowerThreshold > m_UpperThreshold) {
228 clitkExceptionMacro("lower threshold cannot be greater than upper threshold.");
231 // Threshold to get air
232 typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> InputBinarizeFilterType;
233 typename InputBinarizeFilterType::Pointer binarizeFilter = InputBinarizeFilterType::New();
234 binarizeFilter->SetInput(working_input);
235 if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(m_LowerThreshold);
236 // binarizeFilter->CanRunInPlace() is false
237 binarizeFilter->SetUpperThreshold(m_UpperThreshold);
238 binarizeFilter->SetInsideValue(this->GetForegroundValue());
239 binarizeFilter->SetOutsideValue(this->GetBackgroundValue());
240 binarizeFilter->Update();
241 working_mask = binarizeFilter->GetOutput();
242 PrintMemory(GetVerboseMemoryFlag(), "After Binarizefilter"); // OK, additional mem is one mask image
244 // Labelize and keep right labels
246 Labelize<MaskImageType>
247 (working_mask, GetBackgroundValue(), true, GetMinimalComponentSize());
248 PrintMemory(GetVerboseMemoryFlag(), "After Labelize"); // BUG ? additional mem around 1 time the input ?
250 working_mask = RemoveLabels<MaskImageType>
251 (working_mask, GetBackgroundValue(), GetLabelizeParameters1()->GetLabelsToRemove());
252 PrintMemory(GetVerboseMemoryFlag(), "After RemoveLabels"); // OK additional mem = 0
254 working_mask = KeepLabels<MaskImageType>
256 GetBackgroundValue(),
257 GetForegroundValue(),
258 GetLabelizeParameters1()->GetFirstKeep(),
259 GetLabelizeParameters1()->GetLastKeep(),
260 GetLabelizeParameters1()->GetUseLastKeep());
261 PrintMemory(GetVerboseMemoryFlag(), "After KeepLabels to create the 'air'");
264 working_input = SetBackground<ImageType, MaskImageType>
265 (working_input, working_mask, this->GetForegroundValue(), this->GetBackgroundValue(), true);
266 PrintMemory(GetVerboseMemoryFlag(), "After SetBackground");
269 StopCurrentStep<ImageType>(working_input);
271 //--------------------------------------------------------------------
272 //--------------------------------------------------------------------
273 StartNewStep("Search for the trachea");
275 PrintMemory(GetVerboseMemoryFlag(), "After SearchForTrachea");
276 if (m_Seeds.empty()) {
277 clitkExceptionMacro("No seeds for trachea... Aborting.");
280 //--------------------------------------------------------------------
281 //--------------------------------------------------------------------
282 StartNewStep("Extract the lung with Otsu filter");
283 // Automated Otsu thresholding and relabeling
284 typedef itk::OtsuThresholdImageFilter<ImageType,MaskImageType> OtsuThresholdImageFilterType;
285 typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New();
286 otsuFilter->SetInput(working_input);
287 otsuFilter->SetNumberOfHistogramBins(GetNumberOfHistogramBins());
288 otsuFilter->SetInsideValue(this->GetForegroundValue());
289 otsuFilter->SetOutsideValue(this->GetBackgroundValue());
290 otsuFilter->Update();
291 working_mask = otsuFilter->GetOutput();
294 StopCurrentStep<MaskImageType>(working_mask);
295 PrintMemory(GetVerboseMemoryFlag(), "After Otsufilter");
297 //--------------------------------------------------------------------
298 //--------------------------------------------------------------------
299 StartNewStep("Select labels");
301 working_mask = LabelizeAndSelectLabels<MaskImageType>
303 GetBackgroundValue(),
304 GetForegroundValue(),
306 GetMinimalComponentSize(),
307 GetLabelizeParameters2());
310 StopCurrentStep<MaskImageType>(working_mask);
311 PrintMemory(GetVerboseMemoryFlag(), "After LabelizeAndSelectLabels");
313 //--------------------------------------------------------------------
314 //--------------------------------------------------------------------
315 if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
316 StartNewStep("Remove the trachea");
318 working_mask = SetBackground<MaskImageType, MaskImageType>
319 (working_mask, trachea, 1, -1, true);
320 PrintMemory(GetVerboseMemoryFlag(), "After SetBackground");
322 // Dilate the trachea
323 static const unsigned int Dim = ImageType::ImageDimension;
324 typedef itk::BinaryBallStructuringElement<InternalPixelType, Dim> KernelType;
325 KernelType structuringElement;
326 structuringElement.SetRadius(GetRadiusForTrachea());
327 structuringElement.CreateStructuringElement();
328 typedef clitk::ConditionalBinaryDilateImageFilter<MaskImageType, MaskImageType, KernelType> ConditionalBinaryDilateImageFilterType;
329 typename ConditionalBinaryDilateImageFilterType::Pointer dilateFilter = ConditionalBinaryDilateImageFilterType::New();
330 dilateFilter->SetBoundaryToForeground(false);
331 dilateFilter->SetKernel(structuringElement);
332 dilateFilter->SetBackgroundValue (1);
333 dilateFilter->SetForegroundValue (-1);
334 dilateFilter->SetInput (working_mask);
335 dilateFilter->Update();
336 working_mask = dilateFilter->GetOutput();
337 PrintMemory(GetVerboseMemoryFlag(), "After dilate");
339 // Set trachea with dilatation
340 trachea = SetBackground<MaskImageType, MaskImageType>
341 (trachea, working_mask, -1, this->GetForegroundValue(), true);
343 // Remove the trachea
344 working_mask = SetBackground<MaskImageType, MaskImageType>
345 (working_mask, working_mask, -1, this->GetBackgroundValue(), true);
348 working_mask = LabelizeAndSelectLabels<MaskImageType>
350 GetBackgroundValue(),
351 GetForegroundValue(),
353 GetMinimalComponentSize(),
354 GetLabelizeParameters3());
357 StopCurrentStep<MaskImageType>(working_mask);
360 //--------------------------------------------------------------------
361 //--------------------------------------------------------------------
362 PrintMemory(GetVerboseMemoryFlag(), "before autocropfilter");
363 if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
365 trachea = clitk::AutoCrop<MaskImageType>(trachea, GetBackgroundValue());
368 // Remove Padding region
369 typedef itk::CropImageFilter<MaskImageType, MaskImageType> CropFilterType;
370 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
371 cropFilter->SetInput(trachea);
372 cropFilter->SetLowerBoundaryCropSize(bounds);
373 cropFilter->SetUpperBoundaryCropSize(bounds);
374 cropFilter->Update();
375 trachea = cropFilter->GetOutput();
377 StopCurrentStep<MaskImageType>(trachea);
378 PrintMemory(GetVerboseMemoryFlag(), "after delete trachea");
380 PrintMemory(GetVerboseMemoryFlag(), "after delete trachea");
382 //--------------------------------------------------------------------
383 //--------------------------------------------------------------------
384 StartNewStep("Cropping lung");
385 PrintMemory(GetVerboseMemoryFlag(), "Before Autocropfilter");
387 working_mask = clitk::AutoCrop<MaskImageType>(working_mask, GetBackgroundValue());
390 // Remove Padding region
391 typedef itk::CropImageFilter<MaskImageType, MaskImageType> CropFilterType;
392 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
393 cropFilter->SetInput(working_mask);
394 cropFilter->SetLowerBoundaryCropSize(bounds);
395 cropFilter->SetUpperBoundaryCropSize(bounds);
396 cropFilter->Update();
397 working_mask = cropFilter->GetOutput();
399 StopCurrentStep<MaskImageType>(working_mask);
401 //--------------------------------------------------------------------
402 //--------------------------------------------------------------------
404 if (GetOpenCloseFlag()) {
405 StartNewStep("Open/Close");
406 PrintMemory(GetVerboseMemoryFlag(), "Before OpenClose");
408 // Structuring element
409 typedef itk::BinaryBallStructuringElement<InternalPixelType, ImageDimension> KernelType;
410 KernelType structuringElement;
411 structuringElement.SetRadius(GetOpenCloseRadius());
412 structuringElement.CreateStructuringElement();
415 typedef itk::BinaryMorphologicalOpeningImageFilter<MaskImageType, InternalImageType, KernelType> OpenFilterType;
416 typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
417 openFilter->SetInput(working_mask);
418 openFilter->SetBackgroundValue(GetBackgroundValue());
419 openFilter->SetForegroundValue(GetForegroundValue());
420 openFilter->SetKernel(structuringElement);
423 typedef itk::BinaryMorphologicalClosingImageFilter<MaskImageType, MaskImageType, KernelType> CloseFilterType;
424 typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
425 closeFilter->SetInput(openFilter->GetOutput());
426 closeFilter->SetSafeBorder(true);
427 closeFilter->SetForegroundValue(GetForegroundValue());
428 closeFilter->SetKernel(structuringElement);
429 closeFilter->Update();
430 working_mask = closeFilter->GetOutput();
433 //--------------------------------------------------------------------
434 //--------------------------------------------------------------------
436 if (GetFillHolesFlag()) {
437 StartNewStep("Fill Holes");
438 PrintMemory(GetVerboseMemoryFlag(), "Before Fill Holes");
439 typedef clitk::FillMaskFilter<MaskImageType> FillMaskFilterType;
440 typename FillMaskFilterType::Pointer fillMaskFilter = FillMaskFilterType::New();
441 fillMaskFilter->SetInput(working_mask);
442 fillMaskFilter->AddDirection(2);
443 //fillMaskFilter->AddDirection(1);
444 fillMaskFilter->Update();
445 working_mask = fillMaskFilter->GetOutput();
446 StopCurrentStep<MaskImageType>(working_mask);
449 if (GetSeparateLungsFlag()) {
450 //--------------------------------------------------------------------
451 //--------------------------------------------------------------------
452 StartNewStep("Separate Left/Right lungs");
453 PrintMemory(GetVerboseMemoryFlag(), "Before Separate");
455 working_mask = Labelize<MaskImageType>(working_mask,
456 GetBackgroundValue(),
458 GetMinimalComponentSize());
460 PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
463 typedef itk::StatisticsImageFilter<MaskImageType> StatisticsImageFilterType;
464 typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
465 statisticsImageFilter->SetInput(working_mask);
466 statisticsImageFilter->Update();
467 unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
468 working_mask = statisticsImageFilter->GetOutput();
469 PrintMemory(GetVerboseMemoryFlag(), "After count label");
471 // If already 2 labels, but a too big differences, remove the
472 // smalest one (sometimes appends with the stomach
473 if (initialNumberOfLabels >1) {
474 if (GetRemoveSmallLabelBeforeSeparationFlag()) {
475 typedef itk::RelabelComponentImageFilter<MaskImageType, MaskImageType> RelabelFilterType;
476 typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
477 relabelFilter->SetInput(working_mask);
478 relabelFilter->SetMinimumObjectSize(10);
479 relabelFilter->Update();
480 const std::vector<float> & a = relabelFilter->GetSizeOfObjectsInPhysicalUnits();
481 std::vector<MaskImagePixelType> remove_label;
482 for(unsigned int i=1; i<a.size(); i++) {
483 if (a[i] < 0.5*a[0]) { // more than 0.5 difference
484 remove_label.push_back(i+1); // label zero is BG
488 clitk::RemoveLabels<MaskImageType>(working_mask, GetBackgroundValue(), remove_label);
489 statisticsImageFilter->SetInput(working_mask);
490 statisticsImageFilter->Update();
491 initialNumberOfLabels = statisticsImageFilter->GetMaximum();
495 // Decompose the first label
496 if (initialNumberOfLabels<2) {
497 // Structuring element radius
498 typename ImageType::SizeType radius;
499 for (unsigned int i=0;i<Dim;i++) radius[i]=1;
500 typedef clitk::DecomposeAndReconstructImageFilter<MaskImageType,MaskImageType> DecomposeAndReconstructFilterType;
501 typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
502 decomposeAndReconstructFilter->SetInput(working_mask);
503 decomposeAndReconstructFilter->SetVerbose(false);
504 decomposeAndReconstructFilter->SetRadius(radius);
505 decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
506 decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
507 decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
508 decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
509 decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
510 decomposeAndReconstructFilter->SetFullyConnected(true);
511 decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
512 decomposeAndReconstructFilter->Update();
513 working_mask = decomposeAndReconstructFilter->GetOutput();
515 PrintMemory(GetVerboseMemoryFlag(), "After decomposeAndReconstructFilter");
518 // Retain labels ('1' is largset lung, so right. '2' is left)
519 typedef itk::ThresholdImageFilter<MaskImageType> ThresholdImageFilterType;
520 typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
521 thresholdFilter->SetInput(working_mask);
522 thresholdFilter->ThresholdAbove(2);
523 thresholdFilter->SetOutsideValue(this->GetBackgroundValue());
524 thresholdFilter->Update();
525 working_mask = thresholdFilter->GetOutput();
526 StopCurrentStep<MaskImageType> (working_mask);
527 PrintMemory(GetVerboseMemoryFlag(), "After Thresholdfilter");
529 // Update output info
530 // output = working_mask;
531 //this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
533 // this->GetOutput(0)->SetRegions(working_mask->GetLargestPossibleRegion());
536 //--------------------------------------------------------------------
539 //--------------------------------------------------------------------
540 template <class TImageType>
542 clitk::ExtractLungFilter<TImageType>::
543 GenerateInputRequestedRegion() {
544 // DD("GenerateInputRequestedRegion (nothing?)");
546 //--------------------------------------------------------------------
549 //--------------------------------------------------------------------
550 template <class ImageType>
552 clitk::ExtractLungFilter<ImageType>::
556 // this->GraftOutput(output); // not SetNthOutput
557 this->GraftOutput(working_mask); // not SetNthOutput
558 // Store image filenames into AFDB
559 GetAFDB()->SetImageFilename("Lungs", this->GetOutputLungFilename());
560 GetAFDB()->SetImageFilename("Trachea", this->GetOutputTracheaFilename());
563 //--------------------------------------------------------------------
566 //--------------------------------------------------------------------
567 template <class ImageType>
569 clitk::ExtractLungFilter<ImageType>::
570 SearchForTracheaSeed(int skip)
572 if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)
573 // Restart the search until a seed is found, skipping more and more slices
576 // Search seed (parameters = UpperThresholdForTrachea)
577 static const unsigned int Dim = ImageType::ImageDimension;
578 typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
579 typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
580 typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
581 sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-skip-5;
582 sliceRegionSize[Dim-1]=5;
583 sliceRegion.SetSize(sliceRegionSize);
584 sliceRegion.SetIndex(sliceRegionIndex);
585 typedef itk::ImageRegionConstIterator<ImageType> IteratorType;
586 IteratorType it(working_input, sliceRegion);
588 while (!it.IsAtEnd()) {
589 if(it.Get() < GetUpperThresholdForTrachea() ) {
590 AddSeedInPixels(it.GetIndex());
591 // DD(it.GetIndex());
596 // if we do not found : restart
597 stop = (m_Seeds.size() != 0);
599 if (GetVerboseStepFlag()) {
600 std::cout << "\t No seed found this time. I skip some slices and restart." << std::endl;
602 if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
603 // we want to skip more than a half of the image, it is probably a bug
604 std::cerr << "Number of slices to skip to find trachea to high = " << skip << std::endl;
611 // DD(m_Seeds.size());
615 return (m_Seeds.size() != 0);
617 //--------------------------------------------------------------------
620 bool is_orientation_superior(itk::SpatialOrientation::ValidCoordinateOrientationFlags orientation)
622 itk::SpatialOrientation::CoordinateTerms sup = itk::SpatialOrientation::ITK_COORDINATE_Superior;
623 bool primary = (orientation & 0x0000ff) == sup;
624 bool secondary = ((orientation & 0x00ff00) >> 8) == sup;
625 bool tertiary = ((orientation & 0xff0000) >> 16) == sup;
626 return primary || secondary || tertiary;
629 //--------------------------------------------------------------------
630 template <class ImageType>
632 clitk::ExtractLungFilter<ImageType>::
633 SearchForTracheaSeed2(int numberOfSlices)
635 if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)
636 if (GetVerboseRegionGrowingFlag())
637 std::cout << "SearchForTracheaSeed2(" << numberOfSlices << ", " << GetMaxElongation() << ")" << std::endl;
639 typedef unsigned char MaskPixelType;
640 typedef itk::Image<MaskPixelType, ImageType::ImageDimension> MaskImageType;
641 typedef itk::Image<typename MaskImageType::PixelType, 2> MaskImageType2D;
642 typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> ThresholdFilterType;
643 typedef itk::BinaryBallStructuringElement<MaskPixelType, MaskImageType::ImageDimension> KernelType;
644 typedef itk::BinaryMorphologicalClosingImageFilter<MaskImageType, MaskImageType, KernelType> ClosingFilterType;
645 typedef itk::BinaryMorphologicalOpeningImageFilter<MaskImageType, MaskImageType, KernelType> OpeningFilterType;
646 typedef itk::ExtractImageFilter<MaskImageType, MaskImageType2D> SlicerFilterType;
647 typedef itk::ConnectedComponentImageFilter<MaskImageType2D, MaskImageType2D> LabelFilterType;
648 typedef itk::ShapeLabelObject<MaskPixelType, MaskImageType2D::ImageDimension> ShapeLabelType;
649 typedef itk::LabelMap<ShapeLabelType> LabelImageType;
650 typedef itk::LabelImageToShapeLabelMapFilter<MaskImageType2D, LabelImageType> ImageToLabelMapFilterType;
651 typedef itk::LabelMapToLabelImageFilter<LabelImageType, MaskImageType2D> LabelMapToImageFilterType;
652 typedef itk::ImageFileWriter<MaskImageType2D> FileWriterType;
654 // threshold to isolate airawys and lungs
655 typename ThresholdFilterType::Pointer threshold = ThresholdFilterType::New();
656 threshold->SetLowerThreshold(-2000);
657 threshold->SetUpperThreshold(GetSeedPreProcessingThreshold());
658 threshold->SetInput(working_input);
661 KernelType kernel_closing, kernel_opening;
663 // remove small noise
664 typename OpeningFilterType::Pointer opening = OpeningFilterType::New();
665 kernel_opening.SetRadius(1);
666 opening->SetKernel(kernel_opening);
667 opening->SetInput(threshold->GetOutput());
670 typename SlicerFilterType::Pointer slicer = SlicerFilterType::New();
671 slicer->SetDirectionCollapseToIdentity();
672 slicer->SetInput(opening->GetOutput());
675 typename LabelFilterType::Pointer label_filter = LabelFilterType::New();
676 label_filter->SetInput(slicer->GetOutput());
678 // extract shape information from labels
679 typename ImageToLabelMapFilterType::Pointer label_to_map_filter = ImageToLabelMapFilterType::New();
680 label_to_map_filter->SetInput(label_filter->GetOutput());
682 typename LabelMapToImageFilterType::Pointer map_to_label_filter = LabelMapToImageFilterType::New();
683 typename FileWriterType::Pointer writer = FileWriterType::New();
685 typename ImageType::IndexType index;
686 typename ImageType::RegionType region = working_input->GetLargestPossibleRegion();
687 typename ImageType::SizeType size = region.GetSize();
688 typename ImageType::SpacingType spacing = working_input->GetSpacing();
689 typename ImageType::PointType origin = working_input->GetOrigin();
691 int nslices = min(numberOfSlices, size[2]);
692 int start = 0, increment = 1;
693 itk::SpatialOrientationAdapter orientation;
694 typename ImageType::DirectionType dir = working_input->GetDirection();
695 if (!is_orientation_superior(orientation.FromDirectionCosines(dir))) {
700 typename MaskImageType::PointType image_centre;
701 image_centre[0] = size[0]/2;
702 image_centre[1] = size[1]/2;
705 typedef InternalIndexType SeedType;
706 SeedType trachea_centre, shape_centre, max_e_centre, prev_e_centre;
707 typedef std::list<SeedType> PointListType;
708 typedef std::list<PointListType> SequenceListType;
709 PointListType* current_sequence = NULL;
710 SequenceListType sequence_list;
712 prev_e_centre.Fill(0);
713 std::ostringstream file_name;
714 index[0] = index[1] = 0;
715 size[0] = size[1] = 512;
721 region.SetIndex(index);
722 region.SetSize(size);
723 slicer->SetExtractionRegion(region);
725 label_filter->SetInput(slicer->GetOutput());
726 label_filter->Update();
728 label_to_map_filter->SetInput(label_filter->GetOutput());
729 label_to_map_filter->Update();
730 typename LabelImageType::Pointer label_map = label_to_map_filter->GetOutput();
732 if (GetWriteStepFlag()) {
733 map_to_label_filter->SetInput(label_map);
734 writer->SetInput(map_to_label_filter->GetOutput());
736 file_name << "labels_";
739 file_name << index[2] << ".mhd";
740 writer->SetFileName(file_name.str().c_str());
744 typename ShapeLabelType::Pointer shape, max_e_shape;
745 double max_elongation = GetMaxElongation();
746 double max_size = size[0]*size[1]/128;
749 unsigned int nlables = label_map->GetNumberOfLabelObjects();
750 for (unsigned int j = 0; j < nlables; j++) {
751 shape = label_map->GetNthLabelObject(j);
752 if (shape->Size() > 150 && shape->Size() <= max_size) {
753 double e = 1 - 1/shape->GetElongation();
754 //double area = 1 - r->Area() ;
755 if (e < max_elongation) {
757 shape_centre[0] = (shape->GetCentroid()[0] - origin[0])/spacing[0];
758 shape_centre[1] = (shape->GetCentroid()[1] - origin[1])/spacing[1];
759 shape_centre[2] = index[2];
760 //double d = 1 - (shape_centre - image_centre).Magnitude()/max_dist;
761 double dx = shape_centre[0] - image_centre[0];
762 double d = 1 - dx*2/size[0];
768 max_e_centre = shape_centre;
776 itk::Point<typename SeedType::IndexValueType, ImageType::ImageDimension> p1, p2;
777 p1[0] = max_e_centre[0];
778 p1[1] = max_e_centre[1];
779 p1[2] = max_e_centre[2];
781 p2[0] = prev_e_centre[0];
782 p2[1] = prev_e_centre[1];
783 p2[2] = prev_e_centre[2];
785 double mag = (p2 - p1).GetNorm();
786 if (GetVerboseRegionGrowingFlag()) {
788 cout << index[2] << ": ";
789 cout << "region(" << max_e_centre[0] << ", " << max_e_centre[1] << ", " << max_e_centre[2] << "); ";
790 cout << "prev_region(" << prev_e_centre[0] << ", " << prev_e_centre[1] << ", " << prev_e_centre[2] << "); ";
791 cout << "mag(" << mag << "); " << endl;
796 PointListType point_list;
797 point_list.push_back(max_e_centre);
798 sequence_list.push_back(point_list);
799 current_sequence = &sequence_list.back();
801 else if (current_sequence)
802 current_sequence->push_back(max_e_centre);
804 prev_e_centre= max_e_centre;
807 if (GetVerboseRegionGrowingFlag()) {
808 cout << "No shapes found at slice " << index[2] << std::endl;
814 for (typename SequenceListType::iterator s = sequence_list.begin(); s != sequence_list.end(); s++)
816 if (s->size() > longest)
819 trachea_centre = s->front();
824 if (GetVerboseRegionGrowingFlag())
825 std::cout << "seed at: " << trachea_centre << std::endl;
826 m_Seeds.push_back(trachea_centre);
830 return (m_Seeds.size() != 0);
832 //--------------------------------------------------------------------
834 //--------------------------------------------------------------------
835 template <class ImageType>
837 clitk::ExtractLungFilter<ImageType>::
838 TracheaRegionGrowing()
840 // Explosion controlled region growing
841 PrintMemory(GetVerboseMemoryFlag(), "Before ExplosionControlledThresholdConnectedImageFilter");
842 typedef clitk::ExplosionControlledThresholdConnectedImageFilter<ImageType, MaskImageType> ImageFilterType;
843 typename ImageFilterType::Pointer f= ImageFilterType::New();
844 f->SetInput(working_input);
846 f->SetUpper(GetUpperThresholdForTrachea());
847 f->SetMinimumLowerThreshold(-2000);
848 // f->SetMaximumUpperThreshold(0); // MAYBE TO CHANGE ???
849 f->SetMaximumUpperThreshold(-300); // MAYBE TO CHANGE ???
850 f->SetAdaptLowerBorder(false);
851 f->SetAdaptUpperBorder(true);
852 f->SetMinimumSize(5000);
853 f->SetReplaceValue(1);
854 f->SetMultiplier(GetMultiplierForTrachea());
855 f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
856 f->SetMinimumThresholdStepSize(1);
857 f->SetVerbose(GetVerboseRegionGrowingFlag());
858 for(unsigned int i=0; i<m_Seeds.size();i++) {
859 f->AddSeed(m_Seeds[i]);
863 PrintMemory(GetVerboseMemoryFlag(), "After RG update");
865 // take first (main) connected component
866 trachea = Labelize<MaskImageType>(f->GetOutput(),
867 GetBackgroundValue(),
869 1000);//GetMinimalComponentSize());
870 PrintMemory(GetVerboseMemoryFlag(), "After Labelize");
871 trachea = KeepLabels<MaskImageType>(trachea,
872 GetBackgroundValue(),
873 GetForegroundValue(),
875 PrintMemory(GetVerboseMemoryFlag(), "After KeepLabels");
877 //--------------------------------------------------------------------
880 //--------------------------------------------------------------------
881 template <class ImageType>
883 clitk::ExtractLungFilter<ImageType>::
884 ComputeTracheaVolume()
886 typedef itk::ImageRegionConstIterator<InternalImageType> IteratorType;
887 IteratorType iter(trachea, trachea->GetLargestPossibleRegion());
890 while (!iter.IsAtEnd()) {
891 if (iter.Get() == this->GetForegroundValue()) volume++;
895 double voxelsize = trachea->GetSpacing()[0]*trachea->GetSpacing()[1]*trachea->GetSpacing()[2];
896 return volume*voxelsize;
898 //--------------------------------------------------------------------
901 //--------------------------------------------------------------------
902 template <class ImageType>
904 clitk::ExtractLungFilter<ImageType>::
907 // Search for seed among n slices, skip some slices before starting
908 // if not found -> skip more and restart
910 // when seed found : perform region growing
911 // compute trachea volume
912 // if volume not plausible -> skip more slices and restart
914 // If initial seed, convert from mm to pixels
915 if (m_SeedsInMM.size() > 0) {
916 for(unsigned int i=0; i<m_SeedsInMM.size(); i++) {
917 InputImageIndexType index;
918 working_input->TransformPhysicalPointToIndex(m_SeedsInMM[i], index);
919 m_Seeds.push_back(index);
926 int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed();
929 if (GetTracheaSeedAlgorithm() == 0)
930 has_seed = SearchForTracheaSeed(skip);
932 has_seed = SearchForTracheaSeed2(GetNumSlices());
935 TracheaRegionGrowing();
936 volume = ComputeTracheaVolume()/1000; // assume mm3, so divide by 1000 to get cc
937 if (GetWriteStepFlag()) {
938 writeImage<MaskImageType>(trachea, "step-trachea-"+toString(skip)+".mhd");
940 if (GetTracheaVolumeMustBeCheckedFlag()) {
941 if ((volume > 10) && (volume < 65 )) { // depend on image size ...
942 // Typical volume 22.59 cm 3 (± 7.69 cm 3 ) [ Leader 2004 ]
943 if (GetVerboseStepFlag())
945 std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl;
949 if (GetTracheaSeedAlgorithm() == 0) {
950 if (GetVerboseStepFlag()) {
951 std::cout << "\t The volume of the trachea (" << volume
952 << " cc) seems not correct. I skip some slices (" << skip << ") and restart to find seeds."
957 // empty the list of seed
960 if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
961 // we want to skip more than a half of the image, it is probably a bug
962 std::cerr << "2 : Number of slices to skip to find trachea too high = " << skip << std::endl;
974 StopCurrentStep<MaskImageType>(trachea);
976 else { // Trachea not found
977 this->SetWarning("* WARNING * No seed found for trachea.");
981 //--------------------------------------------------------------------
983 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX