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://oncora1.lyon.fnclcc.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"
29 #include "itkBinaryThresholdImageFilter.h"
30 #include "itkConnectedComponentImageFilter.h"
31 #include "itkRelabelComponentImageFilter.h"
32 #include "itkOtsuThresholdImageFilter.h"
33 #include "itkBinaryThinningImageFilter3D.h"
34 #include "itkImageIteratorWithIndex.h"
37 //--------------------------------------------------------------------
38 template <class ImageType, class MaskImageType>
39 clitk::ExtractLungFilter<ImageType, MaskImageType>::
42 itk::ImageToImageFilter<ImageType, MaskImageType>()
45 m_MaxSeedNumber = 500;
47 // Default global options
48 this->SetNumberOfRequiredInputs(2);
49 SetPatientMaskBackgroundValue(0);
50 SetBackgroundValue(0); // Must be zero
51 SetForegroundValue(1);
52 SetMinimalComponentSize(100);
54 // Step 1 default values
55 SetUpperThreshold(-300);
56 SetLowerThreshold(-1000);
57 UseLowerThresholdOff();
58 LabelParamType * p1 = new LabelParamType;
61 p1->AddLabelToRemove(2);
62 SetLabelizeParameters1(p1);
64 // Step 2 default values
65 SetUpperThresholdForTrachea(-900);
66 SetMultiplierForTrachea(5);
67 SetThresholdStepSizeForTrachea(64);
69 // Step 3 default values
70 SetNumberOfHistogramBins(500);
71 LabelParamType * p2 = new LabelParamType;
74 // p->AddLabelToRemove(2); // No label to remove by default
75 SetLabelizeParameters2(p2);
77 // Step 4 default values
78 SetRadiusForTrachea(1);
79 LabelParamType * p3 = new LabelParamType;
83 SetLabelizeParameters3(p3);
85 // Step 5 : find bronchial bifurcations
86 FindBronchialBifurcationsOn();
88 //--------------------------------------------------------------------
91 //--------------------------------------------------------------------
92 template <class ImageType, class MaskImageType>
94 clitk::ExtractLungFilter<ImageType, MaskImageType>::
95 SetInput(const ImageType * image)
97 this->SetNthInput(0, const_cast<ImageType *>(image));
99 //--------------------------------------------------------------------
102 //--------------------------------------------------------------------
103 template <class ImageType, class MaskImageType>
105 clitk::ExtractLungFilter<ImageType, MaskImageType>::
106 SetInputPatientMask(MaskImageType * image, MaskImagePixelType bg )
108 this->SetNthInput(1, const_cast<MaskImageType *>(image));
109 SetPatientMaskBackgroundValue(bg);
111 //--------------------------------------------------------------------
114 //--------------------------------------------------------------------
115 template <class ImageType, class MaskImageType>
117 clitk::ExtractLungFilter<ImageType, MaskImageType>::
118 AddSeed(InternalIndexType s)
120 m_Seeds.push_back(s);
122 //--------------------------------------------------------------------
125 //--------------------------------------------------------------------
126 template <class ImageType, class MaskImageType>
127 template<class ArgsInfoType>
129 clitk::ExtractLungFilter<ImageType, MaskImageType>::
130 SetArgsInfo(ArgsInfoType mArgsInfo)
132 SetVerboseOption_GGO(mArgsInfo);
133 SetVerboseStep_GGO(mArgsInfo);
134 SetWriteStep_GGO(mArgsInfo);
135 SetVerboseWarningOff_GGO(mArgsInfo);
137 SetUpperThreshold_GGO(mArgsInfo);
138 SetLowerThreshold_GGO(mArgsInfo);
139 SetLabelizeParameters1_GGO(mArgsInfo);
140 if (!mArgsInfo.remove1_given) {
141 GetLabelizeParameters1()->AddLabelToRemove(2);
142 if (GetVerboseOption()) VerboseOption("remove1", 2);
145 SetUpperThresholdForTrachea_GGO(mArgsInfo);
146 SetMultiplierForTrachea_GGO(mArgsInfo);
147 SetThresholdStepSizeForTrachea_GGO(mArgsInfo);
148 AddSeed_GGO(mArgsInfo);
150 SetMinimalComponentSize_GGO(mArgsInfo);
152 SetNumberOfHistogramBins_GGO(mArgsInfo);
153 SetLabelizeParameters2_GGO(mArgsInfo);
155 SetRadiusForTrachea_GGO(mArgsInfo);
156 SetLabelizeParameters3_GGO(mArgsInfo);
158 //--------------------------------------------------------------------
161 //--------------------------------------------------------------------
162 template <class ImageType, class MaskImageType>
164 clitk::ExtractLungFilter<ImageType, MaskImageType>::
165 GenerateOutputInformation()
167 Superclass::GenerateOutputInformation();
168 //this->GetOutput(0)->SetRequestedRegion(this->GetOutput(0)->GetLargestPossibleRegion());
170 // Get input pointers
171 patient = dynamic_cast<const MaskImageType*>(itk::ProcessObject::GetInput(1));
172 input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
175 if (!HaveSameSizeAndSpacing<ImageType, MaskImageType>(input, patient)) {
176 this->SetLastError("* ERROR * the images (input and patient mask) must have the same size & spacing");
180 //--------------------------------------------------------------------
181 //--------------------------------------------------------------------
182 StartNewStepOrStop("Set background to initial image");
183 working_input = SetBackground<ImageType, MaskImageType>
184 (input, patient, GetPatientMaskBackgroundValue(), -1000);
185 StopCurrentStep<ImageType>(working_input);
187 //--------------------------------------------------------------------
188 //--------------------------------------------------------------------
189 StartNewStepOrStop("Remove Air");
191 if (m_UseLowerThreshold) {
192 if (m_LowerThreshold > m_UpperThreshold) {
193 this->SetLastError("ERROR: lower threshold cannot be greater than upper threshold.");
197 // Threshold to get air
198 typedef itk::BinaryThresholdImageFilter<ImageType, InternalImageType> InputBinarizeFilterType;
199 typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
200 binarizeFilter->SetInput(working_input);
201 if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(m_LowerThreshold);
202 binarizeFilter->SetUpperThreshold(m_UpperThreshold);
203 binarizeFilter ->SetInsideValue(this->GetForegroundValue());
204 binarizeFilter ->SetOutsideValue(this->GetBackgroundValue());
205 binarizeFilter->Update();
206 working_image = binarizeFilter->GetOutput();
208 // Labelize and keep right labels
209 working_image = Labelize<InternalImageType>(working_image, GetBackgroundValue(), true, GetMinimalComponentSize());
211 working_image = RemoveLabels<InternalImageType>
212 (working_image, GetBackgroundValue(), GetLabelizeParameters1()->GetLabelsToRemove());
214 typename InternalImageType::Pointer air = KeepLabels<InternalImageType>
216 GetBackgroundValue(),
217 GetForegroundValue(),
218 GetLabelizeParameters1()->GetFirstKeep(),
219 GetLabelizeParameters1()->GetLastKeep(),
220 GetLabelizeParameters1()->GetUseLastKeep());
223 working_input = SetBackground<ImageType, InternalImageType>
224 (working_input, air, this->GetForegroundValue(), this->GetBackgroundValue());
227 StopCurrentStep<ImageType>(working_input);
229 //--------------------------------------------------------------------
230 //--------------------------------------------------------------------
231 StartNewStepOrStop("Find the trachea");
232 if (m_Seeds.size() == 0) { // try to find seed
233 // Search seed (parameters = UpperThresholdForTrachea)
234 static const unsigned int Dim = ImageType::ImageDimension;
235 typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
236 typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
237 typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
238 sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-5;
239 sliceRegionSize[Dim-1]=5;
240 sliceRegion.SetSize(sliceRegionSize);
241 sliceRegion.SetIndex(sliceRegionIndex);
242 typedef itk::ImageRegionConstIterator<ImageType> IteratorType;
243 IteratorType it(working_input, sliceRegion);
245 while (!it.IsAtEnd()) {
246 if(it.Get() < GetUpperThresholdForTrachea() ) {
247 AddSeed(it.GetIndex());
253 if (m_Seeds.size() != 0) {
254 // Explosion controlled region growing
255 typedef clitk::ExplosionControlledThresholdConnectedImageFilter<ImageType, InternalImageType> ImageFilterType;
256 typename ImageFilterType::Pointer f= ImageFilterType::New();
257 f->SetInput(working_input);
258 f->SetVerbose(false);
260 f->SetUpper(GetUpperThresholdForTrachea());
261 f->SetMinimumLowerThreshold(-2000);
262 f->SetMaximumUpperThreshold(0);
263 f->SetAdaptLowerBorder(false);
264 f->SetAdaptUpperBorder(true);
265 f->SetMinimumSize(5000);
266 f->SetReplaceValue(1);
267 f->SetMultiplier(GetMultiplierForTrachea());
268 f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
269 f->SetMinimumThresholdStepSize(1);
270 for(unsigned int i=0; i<m_Seeds.size();i++) {
271 f->AddSeed(m_Seeds[i]);
274 trachea_tmp = f->GetOutput();
276 StopCurrentStep<InternalImageType>(trachea_tmp);
278 else { // Trachea not found
279 this->SetWarning("* WARNING * No seed found for trachea.");
283 //--------------------------------------------------------------------
284 //--------------------------------------------------------------------
285 StartNewStepOrStop("Extract the lung with Otsu filter");
286 // Automated Otsu thresholding and relabeling
287 typedef itk::OtsuThresholdImageFilter<ImageType,InternalImageType> OtsuThresholdImageFilterType;
288 typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New();
289 otsuFilter->SetInput(working_input);
290 otsuFilter->SetNumberOfHistogramBins(GetNumberOfHistogramBins());
291 otsuFilter->SetInsideValue(this->GetForegroundValue());
292 otsuFilter->SetOutsideValue(this->GetBackgroundValue());
293 otsuFilter->Update();
294 working_image = otsuFilter->GetOutput();
297 StopCurrentStep<InternalImageType>(working_image);
299 //--------------------------------------------------------------------
300 //--------------------------------------------------------------------
301 StartNewStepOrStop("Select labels");
303 working_image = LabelizeAndSelectLabels<InternalImageType>
305 GetBackgroundValue(),
306 GetForegroundValue(),
308 GetMinimalComponentSize(),
309 GetLabelizeParameters2());
312 StopCurrentStep<InternalImageType>(working_image);
314 //--------------------------------------------------------------------
315 //--------------------------------------------------------------------
316 if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
317 StartNewStepOrStop("Remove the trachea");
319 working_image = SetBackground<InternalImageType, InternalImageType>
320 (working_image, trachea_tmp, 1, -1);
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<InternalImageType, InternalImageType, 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_image);
335 dilateFilter->Update();
336 working_image = dilateFilter->GetOutput();
338 // Set trachea with dilatation
339 trachea_tmp = SetBackground<InternalImageType, InternalImageType>
340 (trachea_tmp, working_image, -1, this->GetForegroundValue());
342 // Remove the trachea
343 working_image = SetBackground<InternalImageType, InternalImageType>
344 (working_image, working_image, -1, this->GetBackgroundValue());
347 working_image = LabelizeAndSelectLabels<InternalImageType>
349 GetBackgroundValue(),
350 GetForegroundValue(),
352 GetMinimalComponentSize(),
353 GetLabelizeParameters3());
356 StopCurrentStep<InternalImageType>(working_image);
359 //--------------------------------------------------------------------
360 //--------------------------------------------------------------------
361 typedef clitk::AutoCropFilter<InternalImageType> CropFilterType;
362 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
363 if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
364 StartNewStepOrStop("Croping trachea");
365 cropFilter->SetInput(trachea_tmp);
366 cropFilter->Update(); // Needed
367 typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
368 typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
369 caster->SetInput(cropFilter->GetOutput());
371 trachea = caster->GetOutput();
372 StopCurrentStep<MaskImageType>(trachea);
375 //--------------------------------------------------------------------
376 //--------------------------------------------------------------------
377 StartNewStepOrStop("Croping lung");
378 typename CropFilterType::Pointer cropFilter2 = CropFilterType::New(); // Needed to reset pipeline
379 cropFilter2->SetInput(working_image);
380 cropFilter2->Update();
381 working_image = cropFilter2->GetOutput();
382 StopCurrentStep<InternalImageType>(working_image);
384 //--------------------------------------------------------------------
385 //--------------------------------------------------------------------
386 StartNewStepOrStop("Separate Left/Right lungs");
388 working_image = Labelize<InternalImageType>(working_image,
389 GetBackgroundValue(),
391 GetMinimalComponentSize());
394 typedef itk::StatisticsImageFilter<InternalImageType> StatisticsImageFilterType;
395 typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
396 statisticsImageFilter->SetInput(working_image);
397 statisticsImageFilter->Update();
398 unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
399 working_image = statisticsImageFilter->GetOutput();
401 // Decompose the first label
402 static const unsigned int Dim = ImageType::ImageDimension;
403 if (initialNumberOfLabels<2) {
404 // Structuring element radius
405 typename ImageType::SizeType radius;
406 for (unsigned int i=0;i<Dim;i++) radius[i]=1;
407 typedef clitk::DecomposeAndReconstructImageFilter<InternalImageType,InternalImageType> DecomposeAndReconstructFilterType;
408 typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
409 decomposeAndReconstructFilter->SetInput(working_image);
410 decomposeAndReconstructFilter->SetVerbose(false);
411 decomposeAndReconstructFilter->SetRadius(radius);
412 decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
413 decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
414 decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
415 decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
416 decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
417 decomposeAndReconstructFilter->SetFullyConnected(true);
418 decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
419 decomposeAndReconstructFilter->Update();
420 working_image = decomposeAndReconstructFilter->GetOutput();
423 // Retain labels (lungs)
424 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
425 typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
426 thresholdFilter->SetInput(working_image);
427 thresholdFilter->ThresholdAbove(2);
428 thresholdFilter->SetOutsideValue(this->GetBackgroundValue());
429 thresholdFilter->Update();
430 working_image = thresholdFilter->GetOutput();
431 StopCurrentStep<InternalImageType> (working_image);
434 StartNewStepOrStop("Final cast");
435 typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
436 typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
437 caster->SetInput(working_image);
439 output = caster->GetOutput();
441 // Update output info
442 this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
444 // Try to extract bifurcation in the trachea (bronchi)
445 // STILL EXPERIMENTAL
446 if (GetFindBronchialBifurcations()) {
447 StartNewStepOrStop("Find bronchial bifurcations");
448 // Step 1 : extract skeleton
449 // Define the thinning filter
450 typedef itk::BinaryThinningImageFilter3D<MaskImageType, MaskImageType> ThinningFilterType;
451 typename ThinningFilterType::Pointer thinningFilter = ThinningFilterType::New();
452 thinningFilter->SetInput(trachea);
453 thinningFilter->Update();
454 typename MaskImageType::Pointer skeleton = thinningFilter->GetOutput();
455 writeImage<MaskImageType>(skeleton, "skeleton.mhd");
460 // Step 2.1 : find first point for tracking
461 typedef itk::ImageRegionConstIteratorWithIndex<MaskImageType> IteratorType;
462 IteratorType it(skeleton, skeleton->GetLargestPossibleRegion());
463 it.GoToReverseBegin();
464 while ((!it.IsAtEnd()) && (it.Get() == GetBackgroundValue())) {
468 this->SetLastError("ERROR: first point in the skeleton not found ! Abort");
471 DD(skeleton->GetLargestPossibleRegion().GetIndex());
472 typename MaskImageType::IndexType index = it.GetIndex();
475 // Step 2.2 : initialize neighborhooditerator
476 typedef itk::NeighborhoodIterator<MaskImageType> NeighborhoodIteratorType;
477 typename NeighborhoodIteratorType::SizeType radius;
479 NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
483 // Find first label number (must be different from BG and FG)
484 typename MaskImageType::PixelType label = GetForegroundValue()+1;
485 while ((label == GetBackgroundValue()) || (label == GetForegroundValue())) { label++; }
488 // Track from the first point
489 std::vector<BifurcationType> listOfBifurcations;
490 TrackFromThisIndex(listOfBifurcations, skeleton, index, label);
492 DD(listOfBifurcations.size());
493 writeImage<MaskImageType>(skeleton, "skeleton2.mhd");
495 for(unsigned int i=0; i<listOfBifurcations.size(); i++) {
496 typename MaskImageType::PointType p;
497 skeleton->TransformIndexToPhysicalPoint(listOfBifurcations[i].index, p);
503 //--------------------------------------------------------------------
506 //--------------------------------------------------------------------
507 template <class ImageType, class MaskImageType>
509 clitk::ExtractLungFilter<ImageType, MaskImageType>::
511 // Do not put some "startnewstep" here, because the object if
512 // modified and the filter's pipeline it do two times. But it is
513 // required to quit if MustStop was set before.
514 if (GetMustStop()) return;
516 // If everything goes well, set the output
517 this->GraftOutput(output); // not SetNthOutput
519 //--------------------------------------------------------------------
522 //--------------------------------------------------------------------
523 template <class ImageType, class MaskImageType>
525 clitk::ExtractLungFilter<ImageType, MaskImageType>::
526 TrackFromThisIndex(std::vector<BifurcationType> & listOfBifurcations,
527 MaskImagePointer skeleton,
528 MaskImageIndexType index,
529 MaskImagePixelType label) {
530 DD("TrackFromThisIndex");
533 // Create NeighborhoodIterator
534 typedef itk::NeighborhoodIterator<MaskImageType> NeighborhoodIteratorType;
535 typename NeighborhoodIteratorType::SizeType radius;
537 NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
540 std::vector<typename NeighborhoodIteratorType::IndexType> listOfTrackedPoint;
543 nit.SetLocation(index);
544 // DD((int)nit.GetCenterPixel());
545 nit.SetCenterPixel(label);
546 listOfTrackedPoint.clear();
547 for(unsigned int i=0; i<nit.Size(); i++) {
548 if (i != nit.GetCenterNeighborhoodIndex ()) { // Do not observe the current point
549 // DD(nit.GetIndex(i));
550 if (nit.GetPixel(i) == GetForegroundValue()) { // if this is foreground, we continue the tracking
551 // DD(nit.GetIndex(i));
552 listOfTrackedPoint.push_back(nit.GetIndex(i));
556 // DD(listOfTrackedPoint.size());
557 if (listOfTrackedPoint.size() == 1) {
558 index = listOfTrackedPoint[0];
561 if (listOfTrackedPoint.size() == 2) {
562 BifurcationType bif(index, label, label+1, label+2);
563 listOfBifurcations.push_back(bif);
564 TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[0], label+1);
565 TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[1], label+2);
568 if (listOfTrackedPoint.size() > 2) {
569 std::cerr << "too much bifurcation points ... ?" << std::endl;
572 // Else this it the end of the tracking
578 //--------------------------------------------------------------------
581 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX