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 CLITKEXTRACTPATIENTFILTER_TXX
20 #define CLITKEXTRACTPATIENTFILTER_TXX
23 #include "clitkImageCommon.h"
24 #include "clitkSetBackgroundImageFilter.h"
25 #include "clitkDecomposeAndReconstructImageFilter.h"
26 #include "clitkAutoCropFilter.h"
27 #include "clitkMemoryUsage.h"
28 #include "clitkSegmentationUtils.h"
31 #include "itkBinaryThresholdImageFilter.h"
32 #include "itkConnectedComponentImageFilter.h"
33 #include "itkRelabelComponentImageFilter.h"
34 #include "itkBinaryMorphologicalClosingImageFilter.h"
35 #include "itkBinaryMorphologicalOpeningImageFilter.h"
36 #include "itkBinaryBallStructuringElement.h"
37 #include "itkCastImageFilter.h"
38 #include "itkConstantPadImageFilter.h"
40 //--------------------------------------------------------------------
41 template <class TInputImageType>
42 clitk::ExtractPatientFilter<TInputImageType>::
43 ExtractPatientFilter():
45 clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
46 itk::ImageToImageFilter<TInputImageType, MaskImageType>()
48 this->SetNumberOfRequiredInputs(1);
49 SetBackgroundValue(0); // Must be zero
50 SetForegroundValue(1);
51 SetPrimaryOpeningRadius(0);
53 // Step 1: Threshold + CC + sort (Find low density areas)
54 SetUpperThreshold(-300);
55 SetLowerThreshold(-1000);
56 UseLowerThresholdOff();
58 // Step 2: DecomposeAndReconstructImageFilter (optional)
59 DecomposeAndReconstructDuringFirstStepOff();
60 InternalImageSizeType r;
63 SetMaximumNumberOfLabels1(2);
64 SetNumberOfNewLabels1(1);
66 // Step 3: Remove the air (largest area).
68 // Step 4: 2nd DecomposeAndReconstructImageFilter
69 DecomposeAndReconstructDuringSecondStepOff();
71 SetMaximumNumberOfLabels2(2);
72 SetNumberOfNewLabels2(1);
74 // Step 5: Only keep label corresponding (Keep patient's labels)
78 // Step 4: OpenClose (option)
82 //--------------------------------------------------------------------
85 //--------------------------------------------------------------------
86 template <class TInputImageType>
88 clitk::ExtractPatientFilter<TInputImageType>::
89 SetInput(const TInputImageType * image)
91 this->SetNthInput(0, const_cast<TInputImageType *>(image));
93 //--------------------------------------------------------------------
96 //--------------------------------------------------------------------
97 template <class TInputImageType>
99 clitk::ExtractPatientFilter<TInputImageType>::
100 GenerateOutputInformation() {
102 clitk::PrintMemory(GetVerboseMemoryFlag(), "Initial memory"); // OK
104 Superclass::GenerateOutputInformation();
105 input = dynamic_cast<const TInputImageType*>(itk::ProcessObject::GetInput(0));
107 // MaskImagePointer outputImage = this->GetOutput(0);
108 // outputImage->SetRegions(input->GetLargestPossibleRegion());
110 // Get input pointers
111 static const unsigned int Dim = InputImageType::ImageDimension;
112 //input = dynamic_cast<const TInputImageType*>(itk::ProcessObject::GetInput(0));
114 //--------------------------------------------------------------------
115 //--------------------------------------------------------------------
117 StartNewStep("Find low densities areas");
119 // Pad images with air to prevent patient touching the image border
120 // But the lungs can touch the outside air. In such case, use primaryOpeningRadius with 1 kernel
121 typedef itk::ConstantPadImageFilter<InputImageType, InputImageType> PadFilterType;
122 typename PadFilterType::Pointer padFilter = PadFilterType::New();
123 padFilter->SetInput(input);
124 padFilter->SetConstant(GetUpperThreshold() - 1);
125 typename InputImageType::SizeType bounds;
126 for (unsigned i = 0; i < Dim - 1; ++i)
129 padFilter->SetPadLowerBound(bounds);
130 padFilter->SetPadUpperBound(bounds);
133 typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> BinarizeFilterType;
134 typename BinarizeFilterType::Pointer binarizeFilter=BinarizeFilterType::New();
135 binarizeFilter->SetInput(padFilter->GetOutput());
136 if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(GetLowerThreshold());
137 binarizeFilter->SetUpperThreshold(GetUpperThreshold());
138 binarizeFilter ->SetInsideValue(this->GetForegroundValue());
139 binarizeFilter ->SetOutsideValue(this->GetBackgroundValue());
140 padFilter->GetOutput()->ReleaseData();
141 working_image = binarizeFilter->GetOutput();
143 typedef itk::BinaryBallStructuringElement<InternalPixelType,Dim> KernelType;
144 unsigned int radius = this->GetPrimaryOpeningRadius();
147 if (this->GetVerboseOptionFlag()) std::cout << ("Opening after threshold; R = ") << radius << std::endl;
149 kernel.SetRadius(radius);
151 typedef itk::BinaryMorphologicalOpeningImageFilter<InternalImageType, InternalImageType , KernelType> OpenFilterType2;
152 typename OpenFilterType2::Pointer openFilter2 = OpenFilterType2::New();
153 openFilter2->SetInput(working_image);
154 openFilter2->SetBackgroundValue(0);
155 openFilter2->SetForegroundValue(1);
156 openFilter2->SetKernel(kernel);
157 openFilter2->Update();
158 working_image->ReleaseData();
159 working_image = openFilter2->GetOutput();
162 if (this->GetVerboseOptionFlag()) std::cout << ("Labelling") << std::endl;
163 // Connected component labeling
164 typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
165 typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
166 connectFilter->SetInput(working_image);
167 connectFilter->SetBackgroundValue(this->GetBackgroundValue());
168 connectFilter->SetFullyConnected(false);
169 connectFilter->Update();
170 working_image->ReleaseData();
171 working_image = connectFilter->GetOutput();
173 if (this->GetVerboseOptionFlag()) std::cout << ("RelabelComponentImageFilter") << std::endl;
174 // Sort labels according to size
175 typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
176 typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
177 relabelFilter->InPlaceOn();
178 relabelFilter->SetInput(connectFilter->GetOutput());
179 relabelFilter->Update();
180 working_image->ReleaseData();
181 working_image = relabelFilter->GetOutput();
184 StopCurrentStep<InternalImageType>(working_image);
186 //--------------------------------------------------------------------
187 //--------------------------------------------------------------------
189 if (GetDecomposeAndReconstructDuringFirstStep()) {
190 StartNewStep("First Decompose & Reconstruct step");
191 typedef clitk::DecomposeAndReconstructImageFilter<InternalImageType,InternalImageType> FilterType;
192 typename FilterType::Pointer f = FilterType::New();
193 f->SetInput(working_image);
194 // f->SetVerbose(m_Verbose);
195 f->SetRadius(GetRadius1());
196 f->SetMaximumNumberOfLabels(GetMaximumNumberOfLabels1());
197 f->SetBackgroundValue(this->GetBackgroundValue());
198 f->SetForegroundValue(this->GetForegroundValue());
199 f->SetFullyConnected(true);
200 f->SetNumberOfNewLabels(GetNumberOfNewLabels1());
202 working_image->ReleaseData();
203 working_image = f->GetOutput();
204 StopCurrentStep<InternalImageType>(working_image);
207 //--------------------------------------------------------------------
208 //--------------------------------------------------------------------
209 if (this->GetVerboseOptionFlag()) std::cout << ("Remove the air (largest area)") << std::endl;
210 StartNewStep("Remove the air (largest area)");
211 typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> iBinarizeFilterType;
212 typename iBinarizeFilterType::Pointer binarizeFilter2 = iBinarizeFilterType::New();
213 binarizeFilter2->SetInput(working_image);
214 binarizeFilter2->SetLowerThreshold(GetFirstKeep());
215 binarizeFilter2->SetUpperThreshold(GetLastKeep());
216 binarizeFilter2 ->SetInsideValue(0);
217 binarizeFilter2 ->SetOutsideValue(1);
218 binarizeFilter2 ->Update();
219 working_image->ReleaseData();
220 working_image = binarizeFilter2->GetOutput();
222 typename ConnectFilterType::Pointer connectFilter2 = ConnectFilterType::New();
223 connectFilter2->SetInput(working_image);
224 connectFilter2->SetBackgroundValue(this->GetBackgroundValue());
225 connectFilter2->SetFullyConnected(false);
226 connectFilter2->Update();
227 working_image->ReleaseData();
228 working_image = connectFilter2->GetOutput();
230 typename RelabelFilterType::Pointer relabelFilter2 = RelabelFilterType::New();
231 relabelFilter2->SetInput(working_image);
232 relabelFilter2->Update();
233 working_image->ReleaseData();
234 working_image = relabelFilter2->GetOutput();
237 working_image = KeepLabels<InternalImageType>
238 (working_image, GetBackgroundValue(), GetForegroundValue(), 1, 1, true);
239 StopCurrentStep<InternalImageType>(working_image);
241 //--------------------------------------------------------------------
242 //--------------------------------------------------------------------
244 if (GetDecomposeAndReconstructDuringSecondStep()) {
245 StartNewStep("Second Decompose & Reconstruct step");
246 typedef clitk::DecomposeAndReconstructImageFilter<InternalImageType,InternalImageType> FilterType;
247 typename FilterType::Pointer f = FilterType::New();
248 f->SetInput(working_image);
249 // f->SetVerbose(m_Verbose);
250 f->SetRadius(GetRadius2());
251 f->SetMaximumNumberOfLabels(GetMaximumNumberOfLabels2());
252 f->SetBackgroundValue(this->GetBackgroundValue());
253 f->SetForegroundValue(this->GetForegroundValue());
254 f->SetFullyConnected(true);
255 f->SetNumberOfNewLabels(GetNumberOfNewLabels2());
257 working_image->ReleaseData();
258 working_image = f->GetOutput();
259 StopCurrentStep<InternalImageType>(working_image);
262 //--------------------------------------------------------------------
263 //--------------------------------------------------------------------
265 if (GetFinalOpenClose()) {
266 StartNewStep("Final OpenClose");
268 KernelType structuringElement;
269 structuringElement.SetRadius(1);
270 structuringElement.CreateStructuringElement();
271 typedef itk::BinaryMorphologicalOpeningImageFilter<InternalImageType, InternalImageType , KernelType> OpenFilterType;
272 typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
273 openFilter->SetInput(working_image);
274 openFilter->SetBackgroundValue(this->GetBackgroundValue());
275 openFilter->SetForegroundValue(this->GetForegroundValue());
276 openFilter->SetKernel(structuringElement);
278 typedef itk::BinaryMorphologicalClosingImageFilter<InternalImageType, InternalImageType , KernelType> CloseFilterType;
279 typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
280 closeFilter->SetInput(openFilter->GetOutput());
281 closeFilter->SetSafeBorder(true);
282 closeFilter->SetForegroundValue(this->GetForegroundValue());
283 // closeFilter->SetBackgroundValue(SetBackgroundValue());
284 closeFilter->SetKernel(structuringElement);
285 closeFilter->Update();
286 working_image->ReleaseData();
287 working_image = closeFilter->GetOutput();
288 StopCurrentStep<InternalImageType>(working_image);
291 //--------------------------------------------------------------------
292 //--------------------------------------------------------------------
294 typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
295 typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
296 caster->SetInput(working_image);
298 working_image->ReleaseData();
299 output = caster->GetOutput();
301 //--------------------------------------------------------------------
302 //--------------------------------------------------------------------
305 StartNewStep("AutoCrop");
306 typedef clitk::AutoCropFilter<MaskImageType> CropFilterType;
307 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
308 cropFilter->SetInput(output);
309 cropFilter->SetBackgroundValue(GetBackgroundValue());
310 cropFilter->Update();
311 output->ReleaseData();
312 output = cropFilter->GetOutput();
313 StopCurrentStep<MaskImageType>(output);
317 // Remove Padding region
318 typedef itk::CropImageFilter<MaskImageType, MaskImageType> CropFilterType;
319 typename CropFilterType::Pointer cropFilter = CropFilterType::New();
320 cropFilter->SetInput(output);
321 cropFilter->SetLowerBoundaryCropSize(bounds);
322 cropFilter->SetUpperBoundaryCropSize(bounds);
323 cropFilter->Update();
324 output->ReleaseData();
325 output = cropFilter->GetOutput();
328 //--------------------------------------------------------------------
331 //--------------------------------------------------------------------
332 template <class TInputImageType>
334 clitk::ExtractPatientFilter<TInputImageType>::
337 this->GraftOutput(output);
338 // Store image filename into AFDB
339 GetAFDB()->SetImageFilename("Patient", this->GetOutputPatientFilename());
342 //--------------------------------------------------------------------
345 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX