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(); // Needed ??
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");
190 // Threshold to get air
191 typedef itk::BinaryThresholdImageFilter<ImageType, InternalImageType> InputBinarizeFilterType;
192 typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
193 binarizeFilter->SetInput(working_input);
194 if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(m_LowerThreshold);
195 binarizeFilter->SetUpperThreshold(m_UpperThreshold);
196 binarizeFilter ->SetInsideValue(this->GetForegroundValue());
197 binarizeFilter ->SetOutsideValue(this->GetBackgroundValue());
198 binarizeFilter->Update();
199 working_image = binarizeFilter->GetOutput();
201 // Labelize and keep right labels
202 working_image = Labelize<InternalImageType>(working_image, GetBackgroundValue(), true, GetMinimalComponentSize());
204 working_image = RemoveLabels<InternalImageType>
205 (working_image, GetBackgroundValue(), GetLabelizeParameters1()->GetLabelsToRemove());
207 typename InternalImageType::Pointer air = KeepLabels<InternalImageType>
209 GetBackgroundValue(),
210 GetForegroundValue(),
211 GetLabelizeParameters1()->GetFirstKeep(),
212 GetLabelizeParameters1()->GetLastKeep(),
213 GetLabelizeParameters1()->GetUseLastKeep());
216 working_input = SetBackground<ImageType, InternalImageType>
217 (working_input, air, this->GetForegroundValue(), this->GetBackgroundValue());
220 StopCurrentStep<ImageType>(working_input);
222 //--------------------------------------------------------------------
223 //--------------------------------------------------------------------
224 StartNewStepOrStop("Find the trachea");
225 if (m_Seeds.size() == 0) { // try to find seed
226 // Search seed (parameters = UpperThresholdForTrachea)
227 static const unsigned int Dim = ImageType::ImageDimension;
228 typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
229 typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
230 typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
231 sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-5;
232 sliceRegionSize[Dim-1]=5;
233 sliceRegion.SetSize(sliceRegionSize);
234 sliceRegion.SetIndex(sliceRegionIndex);
235 typedef itk::ImageRegionConstIterator<ImageType> IteratorType;
236 IteratorType it(working_input, sliceRegion);
238 while (!it.IsAtEnd()) {
239 if(it.Get() < GetUpperThresholdForTrachea() ) {
240 AddSeed(it.GetIndex());
246 if (m_Seeds.size() != 0) {
247 // Explosion controlled region growing
248 typedef clitk::ExplosionControlledThresholdConnectedImageFilter<ImageType, InternalImageType> ImageFilterType;
249 typename ImageFilterType::Pointer f= ImageFilterType::New();
250 f->SetInput(working_input);
251 f->SetVerbose(false);
253 f->SetUpper(GetUpperThresholdForTrachea());
254 f->SetMinimumLowerThreshold(-2000);
255 f->SetMaximumUpperThreshold(0);
256 f->SetAdaptLowerBorder(false);
257 f->SetAdaptUpperBorder(true);
258 f->SetMinimumSize(5000);
259 f->SetReplaceValue(1);
260 f->SetMultiplier(GetMultiplierForTrachea());
261 f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
262 f->SetMinimumThresholdStepSize(1);
263 for(unsigned int i=0; i<m_Seeds.size();i++) {
264 f->AddSeed(m_Seeds[i]);
267 trachea_tmp = f->GetOutput();
269 StopCurrentStep<InternalImageType>(trachea_tmp);
271 else { // Trachea not found
272 this->SetWarning("* WARNING * No seed found for trachea.");
276 //--------------------------------------------------------------------
277 //--------------------------------------------------------------------
278 StartNewStepOrStop("Extract the lung with Otsu filter");
279 // Automated Otsu thresholding and relabeling
280 typedef itk::OtsuThresholdImageFilter<ImageType,InternalImageType> OtsuThresholdImageFilterType;
281 typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New();
282 otsuFilter->SetInput(working_input);
283 otsuFilter->SetNumberOfHistogramBins(GetNumberOfHistogramBins());
284 otsuFilter->SetInsideValue(this->GetForegroundValue());
285 otsuFilter->SetOutsideValue(this->GetBackgroundValue());
286 otsuFilter->Update();
287 working_image = otsuFilter->GetOutput();
290 StopCurrentStep<InternalImageType>(working_image);
292 //--------------------------------------------------------------------
293 //--------------------------------------------------------------------
294 StartNewStepOrStop("Select labels");
296 working_image = LabelizeAndSelectLabels<InternalImageType>
298 GetBackgroundValue(),
299 GetForegroundValue(),
301 GetMinimalComponentSize(),
302 GetLabelizeParameters2());
305 StopCurrentStep<InternalImageType>(working_image);
307 //--------------------------------------------------------------------
308 //--------------------------------------------------------------------
309 if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
310 StartNewStepOrStop("Remove the trachea");
312 working_image = SetBackground<InternalImageType, InternalImageType>
313 (working_image, trachea_tmp, 1, -1);
315 // Dilate the trachea
316 static const unsigned int Dim = ImageType::ImageDimension;
317 typedef itk::BinaryBallStructuringElement<InternalPixelType, Dim> KernelType;
318 KernelType structuringElement;
319 structuringElement.SetRadius(GetRadiusForTrachea());
320 structuringElement.CreateStructuringElement();
321 typedef clitk::ConditionalBinaryDilateImageFilter<InternalImageType, InternalImageType, KernelType> ConditionalBinaryDilateImageFilterType;
322 typename ConditionalBinaryDilateImageFilterType::Pointer dilateFilter = ConditionalBinaryDilateImageFilterType::New();
323 dilateFilter->SetBoundaryToForeground(false);
324 dilateFilter->SetKernel(structuringElement);
325 dilateFilter->SetBackgroundValue (1);
326 dilateFilter->SetForegroundValue (-1);
327 dilateFilter->SetInput (working_image);
328 dilateFilter->Update();
329 working_image = dilateFilter->GetOutput();
331 // Set trachea with dilatation
332 trachea_tmp = SetBackground<InternalImageType, InternalImageType>
333 (trachea_tmp, working_image, -1, this->GetForegroundValue());
335 // Remove the trachea
336 working_image = SetBackground<InternalImageType, InternalImageType>
337 (working_image, working_image, -1, this->GetBackgroundValue());
340 working_image = LabelizeAndSelectLabels<InternalImageType>
342 GetBackgroundValue(),
343 GetForegroundValue(),
345 GetMinimalComponentSize(),
346 GetLabelizeParameters3());
349 StopCurrentStep<InternalImageType>(working_image);
352 //--------------------------------------------------------------------
353 //--------------------------------------------------------------------
354 typedef clitk::AutoCropFilter<InternalImageType> CropFilterType;
355 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
356 if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
357 StartNewStepOrStop("Croping trachea");
358 cropFilter->SetInput(trachea_tmp);
359 cropFilter->Update(); // Needed
360 typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
361 typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
362 caster->SetInput(cropFilter->GetOutput());
364 trachea = caster->GetOutput();
365 StopCurrentStep<MaskImageType>(trachea);
368 //--------------------------------------------------------------------
369 //--------------------------------------------------------------------
370 StartNewStepOrStop("Croping lung");
371 typename CropFilterType::Pointer cropFilter2 = CropFilterType::New(); // Needed to reset pipeline
372 cropFilter2->SetInput(working_image);
373 cropFilter2->Update();
374 working_image = cropFilter2->GetOutput();
375 StopCurrentStep<InternalImageType>(working_image);
377 //--------------------------------------------------------------------
378 //--------------------------------------------------------------------
379 StartNewStepOrStop("Separate Left/Right lungs");
381 working_image = Labelize<InternalImageType>(working_image,
382 GetBackgroundValue(),
384 GetMinimalComponentSize());
387 typedef itk::StatisticsImageFilter<InternalImageType> StatisticsImageFilterType;
388 typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
389 statisticsImageFilter->SetInput(working_image);
390 statisticsImageFilter->Update();
391 unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
392 working_image = statisticsImageFilter->GetOutput();
394 // Decompose the first label
395 static const unsigned int Dim = ImageType::ImageDimension;
396 if (initialNumberOfLabels<2) {
397 // Structuring element radius
398 typename ImageType::SizeType radius;
399 for (unsigned int i=0;i<Dim;i++) radius[i]=1;
400 typedef clitk::DecomposeAndReconstructImageFilter<InternalImageType,InternalImageType> DecomposeAndReconstructFilterType;
401 typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
402 decomposeAndReconstructFilter->SetInput(working_image);
403 decomposeAndReconstructFilter->SetVerbose(false);
404 decomposeAndReconstructFilter->SetRadius(radius);
405 decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
406 decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
407 decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
408 decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
409 decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
410 decomposeAndReconstructFilter->SetFullyConnected(true);
411 decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
412 decomposeAndReconstructFilter->Update();
413 working_image = decomposeAndReconstructFilter->GetOutput();
416 // Retain labels (lungs)
417 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
418 typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
419 thresholdFilter->SetInput(working_image);
420 thresholdFilter->ThresholdAbove(2);
421 thresholdFilter->SetOutsideValue(this->GetBackgroundValue());
422 thresholdFilter->Update();
423 working_image = thresholdFilter->GetOutput();
424 StopCurrentStep<InternalImageType> (working_image);
427 StartNewStepOrStop("Final cast");
428 typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
429 typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
430 caster->SetInput(working_image);
432 output = caster->GetOutput();
434 // Update output info
435 this->GetOutput(0)->SetRegions(output->GetLargestPossibleRegion());
437 // Try to extract bifurcation in the trachea (bronchi)
438 // STILL EXPERIMENTAL
439 if (GetFindBronchialBifurcations()) {
440 StartNewStepOrStop("Find bronchial bifurcations");
441 // Step 1 : extract skeleton
442 // Define the thinning filter
443 typedef itk::BinaryThinningImageFilter3D<MaskImageType, MaskImageType> ThinningFilterType;
444 typename ThinningFilterType::Pointer thinningFilter = ThinningFilterType::New();
445 thinningFilter->SetInput(trachea);
446 thinningFilter->Update();
447 typename MaskImageType::Pointer skeleton = thinningFilter->GetOutput();
448 writeImage<MaskImageType>(skeleton, "skeleton.mhd");
453 // Step 2.1 : find first point for tracking
454 typedef itk::ImageRegionConstIteratorWithIndex<MaskImageType> IteratorType;
455 IteratorType it(skeleton, skeleton->GetLargestPossibleRegion());
456 it.GoToReverseBegin();
457 while ((!it.IsAtEnd()) && (it.Get() == GetBackgroundValue())) {
461 this->SetLastError("ERROR: first point in the skeleton not found ! Abord");
464 DD(skeleton->GetLargestPossibleRegion().GetIndex());
465 typename MaskImageType::IndexType index = it.GetIndex();
468 // Step 2.2 : initialize neighborhooditerator
469 typedef itk::NeighborhoodIterator<MaskImageType> NeighborhoodIteratorType;
470 typename NeighborhoodIteratorType::SizeType radius;
472 NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
476 // Find first label number (must be different from BG and FG)
477 typename MaskImageType::PixelType label = GetForegroundValue()+1;
478 while ((label == GetBackgroundValue()) || (label == GetForegroundValue())) { label++; }
481 // Track from the first point
482 std::vector<BifurcationType> listOfBifurcations;
483 TrackFromThisIndex(listOfBifurcations, skeleton, index, label);
485 DD(listOfBifurcations.size());
486 writeImage<MaskImageType>(skeleton, "skeleton2.mhd");
488 for(unsigned int i=0; i<listOfBifurcations.size(); i++) {
489 typename MaskImageType::PointType p;
490 skeleton->TransformIndexToPhysicalPoint(listOfBifurcations[i].index, p);
496 //--------------------------------------------------------------------
499 //--------------------------------------------------------------------
500 template <class ImageType, class MaskImageType>
502 clitk::ExtractLungFilter<ImageType, MaskImageType>::
504 // Do not put some "startnewstep" here, because the object if
505 // modified and the filter's pipeline it do two times. But it is
506 // required to quit if MustStop was set before.
507 if (GetMustStop()) return;
509 // If everything goes well, set the output
510 this->GraftOutput(output); // not SetNthOutput
512 //--------------------------------------------------------------------
515 //--------------------------------------------------------------------
516 template <class ImageType, class MaskImageType>
518 clitk::ExtractLungFilter<ImageType, MaskImageType>::
519 TrackFromThisIndex(std::vector<BifurcationType> & listOfBifurcations,
520 MaskImagePointer skeleton,
521 MaskImageIndexType index,
522 MaskImagePixelType label) {
523 DD("TrackFromThisIndex");
526 // Create NeighborhoodIterator
527 typedef itk::NeighborhoodIterator<MaskImageType> NeighborhoodIteratorType;
528 typename NeighborhoodIteratorType::SizeType radius;
530 NeighborhoodIteratorType nit(radius, skeleton, skeleton->GetLargestPossibleRegion());
533 std::vector<typename NeighborhoodIteratorType::IndexType> listOfTrackedPoint;
536 nit.SetLocation(index);
537 // DD((int)nit.GetCenterPixel());
538 nit.SetCenterPixel(label);
539 listOfTrackedPoint.clear();
540 for(unsigned int i=0; i<nit.Size(); i++) {
541 if (i != nit.GetCenterNeighborhoodIndex ()) { // Do not observe the current point
542 // DD(nit.GetIndex(i));
543 if (nit.GetPixel(i) == GetForegroundValue()) { // if this is foreground, we continue the tracking
544 // DD(nit.GetIndex(i));
545 listOfTrackedPoint.push_back(nit.GetIndex(i));
549 // DD(listOfTrackedPoint.size());
550 if (listOfTrackedPoint.size() == 1) {
551 index = listOfTrackedPoint[0];
554 if (listOfTrackedPoint.size() == 2) {
555 BifurcationType bif(index, label, label+1, label+2);
556 listOfBifurcations.push_back(bif);
557 TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[0], label+1);
558 TrackFromThisIndex(listOfBifurcations, skeleton, listOfTrackedPoint[1], label+2);
561 if (listOfTrackedPoint.size() > 2) {
562 std::cerr << "too much bifurcation points ... ?" << std::endl;
565 // Else this it the end of the tracking
571 //--------------------------------------------------------------------
574 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX