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 "clitkSegmentationFunctions.h"
26 #include "clitkAutoCropFilter.h"
29 #include "itkBinaryThresholdImageFilter.h"
30 #include "itkConnectedComponentImageFilter.h"
31 #include "itkRelabelComponentImageFilter.h"
32 #include "itkOtsuThresholdImageFilter.h"
34 //--------------------------------------------------------------------
35 template <class TInputImageType, class TMaskImageType>
36 clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
39 itk::ImageToImageFilter<TInputImageType, TMaskImageType>()
41 // Default global options
42 this->SetNumberOfRequiredInputs(2);
43 SetPatientMaskBackgroundValue(0);
44 SetBackgroundValue(0); // Must be zero
45 SetForegroundValue(1);
47 // Step 1 default values
48 SetUpperThreshold(-300);
49 SetLowerThreshold(-1000);
50 UseLowerThresholdOff();
51 LabelParamType * p1 = new LabelParamType;
54 p1->AddLabelToRemove(2);
55 SetLabelizeParameters1(p1);
57 // Step 2 default values
58 SetUpperThresholdForTrachea(-900);
59 SetMultiplierForTrachea(5);
60 SetThresholdStepSizeForTrachea(64);
62 // Step 3 default values
63 SetNumberOfHistogramBins(500);
64 LabelParamType * p2 = new LabelParamType;
67 // p->AddLabelToRemove(2); // No label to remove by default
68 SetLabelizeParameters2(p2);
70 // Step 4 default values
71 SetRadiusForTrachea(1);
72 LabelParamType * p3 = new LabelParamType;
76 SetLabelizeParameters3(p3);
78 //--------------------------------------------------------------------
81 //--------------------------------------------------------------------
82 template <class TInputImageType, class TMaskImageType>
84 clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
85 SetInput(const TInputImageType * image)
87 this->SetNthInput(0, const_cast<TInputImageType *>(image));
89 //--------------------------------------------------------------------
92 //--------------------------------------------------------------------
93 template <class TInputImageType, class TMaskImageType>
95 clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
96 SetInputPatientMask(TMaskImageType * image, MaskImagePixelType bg )
98 this->SetNthInput(1, const_cast<TMaskImageType *>(image));
99 SetPatientMaskBackgroundValue(bg);
101 //--------------------------------------------------------------------
104 //--------------------------------------------------------------------
105 template <class TInputImageType, class TMaskImageType>
107 clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
108 AddSeed(InternalIndexType s)
110 m_Seeds.push_back(s);
112 //--------------------------------------------------------------------
115 //--------------------------------------------------------------------
116 template <class TInputImageType, class TMaskImageType>
117 template<class ArgsInfoType>
119 clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
120 SetArgsInfo(ArgsInfoType mArgsInfo)
122 SetVerboseOption_GGO(mArgsInfo);
123 SetVerboseStep_GGO(mArgsInfo);
124 SetWriteStep_GGO(mArgsInfo);
125 SetVerboseWarningOff_GGO(mArgsInfo);
127 SetUpperThreshold_GGO(mArgsInfo);
128 SetLowerThreshold_GGO(mArgsInfo);
129 SetLabelizeParameters1_GGO(mArgsInfo);
130 if (!mArgsInfo.remove1_given) {
131 GetLabelizeParameters1()->AddLabelToRemove(2);
132 if (GetVerboseOption()) VerboseOption("remove1", 2);
135 SetUpperThresholdForTrachea_GGO(mArgsInfo);
136 SetMultiplierForTrachea_GGO(mArgsInfo);
137 SetThresholdStepSizeForTrachea_GGO(mArgsInfo);
138 AddSeed_GGO(mArgsInfo);
140 SetNumberOfHistogramBins_GGO(mArgsInfo);
141 SetLabelizeParameters2_GGO(mArgsInfo);
143 SetRadiusForTrachea_GGO(mArgsInfo);
144 SetLabelizeParameters3_GGO(mArgsInfo);
146 //--------------------------------------------------------------------
149 //--------------------------------------------------------------------
150 template <class TInputImageType, class TMaskImageType>
152 clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
153 GenerateOutputInformation()
155 input = dynamic_cast<const TInputImageType*>(itk::ProcessObject::GetInput(0));
156 Superclass::GenerateOutputInformation();
157 // MaskImagePointer output = this->GetOutput(0);
159 // Get input pointers
160 input = dynamic_cast<const TInputImageType*>(itk::ProcessObject::GetInput(0));
161 patient = dynamic_cast<const TMaskImageType*>(itk::ProcessObject::GetInput(1));
164 if (!HasSameSizeAndSpacing<TInputImageType, TMaskImageType>(input, patient)) {
165 this->SetLastError("* ERROR * the images (input and patient mask) must have the same size & spacing");
169 //--------------------------------------------------------------------
170 //--------------------------------------------------------------------
171 StartNewStep("Set background to initial image");
172 working_input = SetBackground<TInputImageType, TMaskImageType>
173 (input, patient, GetPatientMaskBackgroundValue(), -1000);
174 StopCurrentStep<InputImageType>(working_input);
176 //--------------------------------------------------------------------
177 //--------------------------------------------------------------------
178 StartNewStep("Remove Air");
179 // Threshold to get air
180 typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
181 typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
182 binarizeFilter->SetInput(working_input);
183 if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(m_LowerThreshold);
184 binarizeFilter->SetUpperThreshold(m_UpperThreshold);
185 binarizeFilter ->SetInsideValue(this->GetForegroundValue());
186 binarizeFilter ->SetOutsideValue(this->GetBackgroundValue());
187 binarizeFilter->Update();
188 working_image = binarizeFilter->GetOutput();
190 // Labelize and keep right labels
191 working_image = Labelize<InternalImageType>(working_image, GetBackgroundValue(), true, GetMinimalComponentSize());
192 working_image = RemoveLabels<InternalImageType>
193 (working_image, GetBackgroundValue(), GetLabelizeParameters1()->GetLabelsToRemove());
194 typename InternalImageType::Pointer air = KeepLabels<InternalImageType>
196 GetBackgroundValue(),
197 GetForegroundValue(),
198 GetLabelizeParameters1()->GetFirstKeep(),
199 GetLabelizeParameters1()->GetLastKeep(),
200 GetLabelizeParameters1()->GetUseLastKeep());
203 working_input = SetBackground<TInputImageType, InternalImageType>
204 (working_input, air, this->GetForegroundValue(), this->GetBackgroundValue());
207 StopCurrentStep<InputImageType>(working_input);
209 //--------------------------------------------------------------------
210 //--------------------------------------------------------------------
211 StartNewStep("Find the trachea");
212 if (m_Seeds.size() == 0) { // try to find seed
213 // Search seed (parameters = UpperThresholdForTrachea)
214 static const unsigned int Dim = InputImageType::ImageDimension;
215 typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
216 typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
217 typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
218 sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-5;
219 sliceRegionSize[Dim-1]=5;
220 sliceRegion.SetSize(sliceRegionSize);
221 sliceRegion.SetIndex(sliceRegionIndex);
223 typedef itk::ImageRegionConstIterator<InputImageType> IteratorType;
224 IteratorType it(working_input, sliceRegion);
226 while (!it.IsAtEnd()) {
227 if(it.Get() < GetUpperThresholdForTrachea() ) {
228 AddSeed(it.GetIndex());
234 if (m_Seeds.size() != 0) {
235 // Explosion controlled region growing
236 typedef clitk::ExplosionControlledThresholdConnectedImageFilter<InputImageType, InternalImageType> ImageFilterType;
237 typename ImageFilterType::Pointer f= ImageFilterType::New();
238 f->SetInput(working_input);
239 f->SetVerbose(false);
241 f->SetUpper(GetUpperThresholdForTrachea());
242 f->SetMinimumLowerThreshold(-2000);
243 f->SetMaximumUpperThreshold(0);
244 f->SetAdaptLowerBorder(false);
245 f->SetAdaptUpperBorder(true);
246 f->SetMinimumSize(5000);
247 f->SetReplaceValue(1);
248 f->SetMultiplier(GetMultiplierForTrachea());
249 f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
250 f->SetMinimumThresholdStepSize(1);
251 for(unsigned int i=0; i<m_Seeds.size();i++) {
252 // std::cout<<"Adding seed " <<m_Seeds[i]<<"..."<<std::endl;
253 f->AddSeed(m_Seeds[i]);
256 trachea_tmp = f->GetOutput();
258 StopCurrentStep<InternalImageType>(trachea_tmp);
261 else { // Trachea not found
262 this->SetWarning("* WARNING * No seed found for trachea.");
267 //--------------------------------------------------------------------
268 //--------------------------------------------------------------------
269 StartNewStep("Extract the lung with Otsu filter");
270 // Automated Otsu thresholding and relabeling
271 typedef itk::OtsuThresholdImageFilter<InputImageType,InternalImageType> OtsuThresholdImageFilterType;
272 typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New();
273 otsuFilter->SetInput(working_input);
274 otsuFilter->SetNumberOfHistogramBins(GetNumberOfHistogramBins());
275 otsuFilter->SetInsideValue(this->GetForegroundValue());
276 otsuFilter->SetOutsideValue(this->GetBackgroundValue());
277 otsuFilter->Update();
278 working_image = otsuFilter->GetOutput();
281 StopCurrentStep<InternalImageType>(working_image);
283 //--------------------------------------------------------------------
284 //--------------------------------------------------------------------
285 StartNewStep("Select labels");
287 working_image = LabelizeAndSelectLabels<InternalImageType>
289 GetBackgroundValue(),
290 GetForegroundValue(),
292 GetMinimalComponentSize(),
293 GetLabelizeParameters2());
296 StopCurrentStep<InternalImageType>(working_image);
298 //--------------------------------------------------------------------
299 //--------------------------------------------------------------------
300 if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
301 StartNewStep("Remove the trachea");
303 working_image = SetBackground<InternalImageType, InternalImageType>
304 (working_image, trachea_tmp, 1, -1);
306 // Dilate the trachea
307 static const unsigned int Dim = InputImageType::ImageDimension;
308 typedef itk::BinaryBallStructuringElement<InternalPixelType, Dim> KernelType;
309 KernelType structuringElement;
310 structuringElement.SetRadius(GetRadiusForTrachea());
311 structuringElement.CreateStructuringElement();
312 typedef clitk::ConditionalBinaryDilateImageFilter<InternalImageType, InternalImageType, KernelType> ConditionalBinaryDilateImageFilterType;
313 typename ConditionalBinaryDilateImageFilterType::Pointer dilateFilter = ConditionalBinaryDilateImageFilterType::New();
314 dilateFilter->SetBoundaryToForeground(false);
315 dilateFilter->SetKernel(structuringElement);
316 dilateFilter->SetBackgroundValue (1);
317 dilateFilter->SetForegroundValue (-1);
318 dilateFilter->SetInput (working_image);
319 dilateFilter->Update();
320 working_image = dilateFilter->GetOutput();
322 // Set trachea with dilatation
323 trachea_tmp = SetBackground<InternalImageType, InternalImageType>
324 (trachea_tmp, working_image, -1, this->GetForegroundValue());
326 // Remove the trachea
327 working_image = SetBackground<InternalImageType, InternalImageType>
328 (working_image, working_image, -1, this->GetBackgroundValue());
331 working_image = LabelizeAndSelectLabels<InternalImageType>
333 GetBackgroundValue(),
334 GetForegroundValue(),
336 GetMinimalComponentSize(),
337 GetLabelizeParameters3());
340 StopCurrentStep<InternalImageType>(working_image);
344 //--------------------------------------------------------------------
345 //--------------------------------------------------------------------
346 typedef clitk::AutoCropFilter<InternalImageType> CropFilterType;
347 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
348 if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
349 StartNewStep("Croping trachea");
350 cropFilter->SetInput(trachea_tmp);
351 cropFilter->Update(); // Needed
352 typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
353 typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
354 caster->SetInput(cropFilter->GetOutput());
356 trachea = caster->GetOutput();
357 StopCurrentStep<MaskImageType>(trachea);
360 //--------------------------------------------------------------------
361 //--------------------------------------------------------------------
362 StartNewStep("Croping lung");
363 cropFilter = CropFilterType::New(); // Needed to reset pipeline
364 cropFilter->SetInput(working_image);
365 cropFilter->Update();
366 working_image = cropFilter->GetOutput();
367 StopCurrentStep<InternalImageType>(working_image);
369 //--------------------------------------------------------------------
370 //--------------------------------------------------------------------
371 StartNewStep("Separate Left/Right lungs");
373 working_image = Labelize<InternalImageType>(working_image,
374 GetBackgroundValue(),
376 GetMinimalComponentSize());
379 typedef itk::StatisticsImageFilter<InternalImageType> StatisticsImageFilterType;
380 typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
381 statisticsImageFilter->SetInput(working_image);
382 statisticsImageFilter->Update();
383 unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
384 working_image = statisticsImageFilter->GetOutput();
386 // Decompose the first label
387 static const unsigned int Dim = InputImageType::ImageDimension;
388 if (initialNumberOfLabels<2) {
389 // Structuring element radius
390 typename InputImageType::SizeType radius;
391 for (unsigned int i=0;i<Dim;i++) radius[i]=1;
392 typedef clitk::DecomposeAndReconstructImageFilter<InternalImageType,InternalImageType> DecomposeAndReconstructFilterType;
393 typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
394 decomposeAndReconstructFilter->SetInput(working_image);
395 decomposeAndReconstructFilter->SetVerbose(false);
396 decomposeAndReconstructFilter->SetRadius(radius);
397 decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
398 decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
399 decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
400 decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
401 decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
402 decomposeAndReconstructFilter->SetFullyConnected(true);
403 decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
404 decomposeAndReconstructFilter->Update();
405 working_image = decomposeAndReconstructFilter->GetOutput();
408 // Retain labels (lungs)
409 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
410 typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
411 thresholdFilter->SetInput(working_image);
412 thresholdFilter->ThresholdAbove(2);
413 thresholdFilter->SetOutsideValue(this->GetBackgroundValue());
414 thresholdFilter->Update();
415 working_image = thresholdFilter->GetOutput();
416 StopCurrentStep<InternalImageType> (working_image);
418 //--------------------------------------------------------------------
421 //--------------------------------------------------------------------
422 template <class TInputImageType, class TMaskImageType>
424 clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
427 typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
428 typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
429 caster->SetInput(working_image);
432 //this->SetNthOutput(0, caster->GetOutput()); // -> no because redo filter otherwise
433 this->GraftOutput(caster->GetOutput());
436 //--------------------------------------------------------------------
439 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX