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 ======================================================================-====*/
18 #ifndef clitkMotionMaskGenericFilter_txx
19 #define clitkMotionMaskGenericFilter_txx
21 /* =================================================
22 * @file clitkMotionMaskGenericFilter.txx
28 ===================================================*/
34 //-------------------------------------------------------------------
35 // Update with the number of dimensions
36 //-------------------------------------------------------------------
37 template<unsigned int Dimension>
39 MotionMaskGenericFilter::UpdateWithDim(std::string PixelType)
41 if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<"..."<<std::endl;
43 if(PixelType == "short") {
44 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
45 UpdateWithDimAndPixelType<Dimension, signed short>();
47 // else if(PixelType == "unsigned_short"){
48 // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
49 // UpdateWithDimAndPixelType<Dimension, unsigned short>();
52 // else if (PixelType == "unsigned_char"){
53 // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
54 // UpdateWithDimAndPixelType<Dimension, unsigned char>();
57 // else if (PixelType == "char"){
58 // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
59 // UpdateWithDimAndPixelType<Dimension, signed char>();
62 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
63 UpdateWithDimAndPixelType<Dimension, float>();
67 //-------------------------------------------------------------------
69 //-------------------------------------------------------------------
70 template <unsigned int Dimension, class PixelType>
71 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
72 MotionMaskGenericFilter::GetAirImage(typename itk::Image<PixelType, Dimension>::Pointer input,
73 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer lungs)
77 typedef int InternalPixelType;
78 typedef itk::Image<PixelType, Dimension> InputImageType;
79 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
81 //----------------------------------------------------------------------------------------------------
82 // Typedef for Processing
83 //----------------------------------------------------------------------------------------------------
84 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
85 typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
86 typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> InversionFilterType;
87 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdFilterType;
88 typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
89 typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
90 typedef itk::MirrorPadImageFilter<InternalImageType, InternalImageType> MirrorPadImageFilterType;
93 typename InternalImageType::Pointer air;
94 if (m_ArgsInfo.featureAir_given) {
95 typename FeatureReaderType::Pointer featureReader =FeatureReaderType::New();
96 featureReader->SetFileName(m_ArgsInfo.featureAir_arg);
97 if (m_Verbose) std::cout<<"Reading the air feature image..."<<std::endl;
98 featureReader->Update();
99 air=featureReader->GetOutput();
101 if (m_Verbose) std::cout<<"Extracting the air feature image..."<<std::endl;
102 //---------------------------------
103 // Binarize the image
104 //---------------------------------
105 typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
106 binarizeFilter->SetInput(input);
107 binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdAir_arg));
108 binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdAir_arg));
109 if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<<m_ArgsInfo.lowerThresholdAir_arg
110 <<", "<<m_ArgsInfo.upperThresholdAir_arg<<"..."<<std::endl;
111 binarizeFilter->Update();
113 //---------------------------------
114 // Remove lungs (in place)
115 //---------------------------------
116 typedef itk::ImageRegionIterator<InternalImageType> IteratorType;
117 IteratorType itAir(binarizeFilter->GetOutput(), binarizeFilter->GetOutput()->GetLargestPossibleRegion());
118 IteratorType itLungs(lungs, binarizeFilter->GetOutput()->GetLargestPossibleRegion()); //lungs padded, used air region
121 while(!itAir.IsAtEnd()) {
122 if(!itLungs.Get()) // The lungs are set to 0 at that stage
128 //---------------------------------
129 // Label the connected components
130 //---------------------------------
131 typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
132 connectFilter->SetInput(binarizeFilter->GetOutput());
133 connectFilter->SetBackgroundValue(0);
134 connectFilter->SetFullyConnected(false);
135 if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
137 //---------------------------------
138 // Sort the labels according to size
139 //---------------------------------
140 typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
141 relabelFilter->SetInput(connectFilter->GetOutput());
142 if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
144 //---------------------------------
146 //---------------------------------
147 typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
148 thresholdFilter->SetInput(relabelFilter->GetOutput());
149 thresholdFilter->SetUpper(1);
150 if (m_Verbose) std::cout<<"Keeping the first label..."<<std::endl;
151 thresholdFilter->Update();
152 air=thresholdFilter->GetOutput();
155 //---------------------------------
157 //---------------------------------
158 typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
159 inversionFilter->SetInput(air);
160 inversionFilter->SetLowerThreshold(0);
161 inversionFilter->SetUpperThreshold(0);
162 inversionFilter->SetInsideValue(1);
163 inversionFilter->Update();
165 //---------------------------------
167 //---------------------------------
168 typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
169 mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
170 typename InternalImageType::SizeType padSize;
172 padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
173 mirrorPadImageFilter->SetPadLowerBound(padSize);
174 if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
175 mirrorPadImageFilter->Update();
176 air=mirrorPadImageFilter->GetOutput();
177 //writeImage<InternalImageType>(air,"/home/srit/tmp/air.mhd");
179 //---------------------------------
181 //---------------------------------
182 if(m_ArgsInfo.pad_flag) {
183 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
184 IteratorType it(air, air->GetLargestPossibleRegion());
185 typename InternalImageType::IndexType index;
186 while(!it.IsAtEnd()) {
188 for (unsigned int i=0; i<Dimension; i++)
189 if(index[i]==air->GetLargestPossibleRegion().GetIndex()[i]
190 || index[i]==(unsigned int)air->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air->GetLargestPossibleRegion().GetSize()[i]-1)
200 //-------------------------------------------------------------------
202 //-------------------------------------------------------------------
203 template <unsigned int Dimension, class PixelType>
204 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
205 MotionMaskGenericFilter::GetBonesImage(typename itk::Image<PixelType, Dimension>::Pointer input )
209 typedef int InternalPixelType;
210 typedef itk::Image<PixelType, Dimension> InputImageType;
211 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
213 //----------------------------------------------------------------------------------------------------
214 // Typedef for Processing
215 //----------------------------------------------------------------------------------------------------
216 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
217 typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
218 typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> InversionFilterType;
219 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdFilterType;
220 typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
221 typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
222 typedef itk::MirrorPadImageFilter<InternalImageType, InternalImageType> MirrorPadImageFilterType;
225 typename InternalImageType::Pointer bones;
226 if (m_ArgsInfo.featureBones_given) {
227 typename FeatureReaderType::Pointer featureReader =FeatureReaderType::New();
228 featureReader->SetFileName(m_ArgsInfo.featureBones_arg);
229 if (m_Verbose) std::cout<<"Reading the bones feature image..."<<std::endl;
230 featureReader->Update();
231 bones=featureReader->GetOutput();
233 if (m_Verbose) std::cout<<"Extracting the bones feature image..."<<std::endl;
234 //---------------------------------
235 // Binarize the image
236 //---------------------------------
237 typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
238 binarizeFilter->SetInput(input);
239 binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdBones_arg));
240 binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdBones_arg));
241 if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<<m_ArgsInfo.lowerThresholdBones_arg
242 <<", "<<m_ArgsInfo.upperThresholdBones_arg<<"..."<<std::endl;
244 //---------------------------------
245 // Label the connected components
246 //---------------------------------
247 typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
248 connectFilter->SetInput(binarizeFilter->GetOutput());
249 connectFilter->SetBackgroundValue(0);
250 connectFilter->SetFullyConnected(false);
251 if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
253 //---------------------------------
254 // Sort the labels according to size
255 //---------------------------------
256 typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
257 relabelFilter->SetInput(connectFilter->GetOutput());
258 if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
260 //---------------------------------
262 //---------------------------------
263 typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
264 thresholdFilter->SetInput(relabelFilter->GetOutput());
265 thresholdFilter->SetUpper(1);
266 if (m_Verbose) std::cout<<"Keeping the first label..."<<std::endl;
267 thresholdFilter->Update();
268 bones=thresholdFilter->GetOutput();
272 //---------------------------------
274 //---------------------------------
275 typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
276 inversionFilter->SetInput(bones);
277 inversionFilter->SetLowerThreshold(0);
278 inversionFilter->SetUpperThreshold(0);
279 inversionFilter->SetInsideValue(1);
280 inversionFilter->Update();
282 //---------------------------------
284 //---------------------------------
285 typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
286 mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
287 typename InternalImageType::SizeType padSize;
289 padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
290 mirrorPadImageFilter->SetPadLowerBound(padSize);
291 if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
292 mirrorPadImageFilter->Update();
293 bones=mirrorPadImageFilter->GetOutput();
294 // writeImage<InternalImageType>(bones,"/home/jef/tmp/bones.mhd");
302 //-------------------------------------------------------------------
304 //-------------------------------------------------------------------
305 template <unsigned int Dimension, class PixelType>
306 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
307 MotionMaskGenericFilter::GetLungsImage(typename itk::Image<PixelType, Dimension>::Pointer input )
310 typedef int InternalPixelType;
311 typedef itk::Image<PixelType, Dimension> InputImageType;
312 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
314 //----------------------------------------------------------------------------------------------------
315 // Typedef for Processing
316 //----------------------------------------------------------------------------------------------------
317 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
318 typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
319 typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> InversionFilterType;
320 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdFilterType;
321 typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
322 typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
323 typedef itk::MirrorPadImageFilter<InternalImageType, InternalImageType> MirrorPadImageFilterType;
325 typename InternalImageType::Pointer lungs;
326 if (m_ArgsInfo.featureLungs_given) {
327 typename FeatureReaderType::Pointer featureReader =FeatureReaderType::New();
328 featureReader->SetFileName(m_ArgsInfo.featureLungs_arg);
329 if (m_Verbose) std::cout<<"Reading the lungs feature image..."<<std::endl;
330 featureReader->Update();
331 lungs=featureReader->GetOutput();
333 if (m_Verbose) std::cout<<"Extracting the lungs feature image..."<<std::endl;
334 //---------------------------------
335 // Binarize the image
336 //---------------------------------
337 typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
338 binarizeFilter->SetInput(input);
339 binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdLungs_arg));
340 binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdLungs_arg));
341 if (m_Verbose) std::cout<<"Binarizing the image using a threshold "<<m_ArgsInfo.lowerThresholdLungs_arg
342 <<", "<<m_ArgsInfo.upperThresholdLungs_arg<<"..."<<std::endl;
344 //---------------------------------
345 // Label the connected components
346 //---------------------------------
347 typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
348 connectFilter->SetInput(binarizeFilter->GetOutput());
349 connectFilter->SetBackgroundValue(0);
350 connectFilter->SetFullyConnected(true);
351 if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
353 //---------------------------------
354 // Sort the labels according to size
355 //---------------------------------
356 typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
357 relabelFilter->SetInput(connectFilter->GetOutput());
358 if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
359 // writeImage<InternalImageType> (relabelFilter->GetOutput(), "/home/jef/tmp/labels.mhd");
361 //---------------------------------
363 //---------------------------------
364 typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
365 thresholdFilter->SetInput(relabelFilter->GetOutput());
366 thresholdFilter->SetLower(1);
367 thresholdFilter->SetUpper(2);
368 if (m_Verbose) std::cout<<"Keeping the first two labels..."<<std::endl;
369 thresholdFilter->Update();
370 lungs=thresholdFilter->GetOutput();
375 //---------------------------------
377 //---------------------------------
378 typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
379 inversionFilter->SetInput(lungs);
380 inversionFilter->SetLowerThreshold(0);
381 inversionFilter->SetUpperThreshold(0);
382 inversionFilter->SetInsideValue(1);
383 inversionFilter->Update();
385 //---------------------------------
387 //---------------------------------
388 typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
389 mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
390 typename InternalImageType::SizeType padSize;
392 padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
393 mirrorPadImageFilter->SetPadLowerBound(padSize);
394 if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
395 mirrorPadImageFilter->Update();
396 lungs=mirrorPadImageFilter->GetOutput();
397 // writeImage<InternalImageType>(lungs,"/home/jef/tmp/lungs.mhd");
402 //-------------------------------------------------------------------
404 //-------------------------------------------------------------------
405 template <unsigned int Dimension, class PixelType>
406 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
407 MotionMaskGenericFilter::Resample( typename itk::Image<InternalPixelType,Dimension>::Pointer input )
410 typedef int InternalPixelType;
411 typedef itk::Image<PixelType, Dimension> InputImageType;
412 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
414 //---------------------------------
415 // Resample to new spacing
416 //---------------------------------
417 typename InternalImageType::SizeType size_low;
418 typename InternalImageType::SpacingType spacing_low;
419 for (unsigned int i=0; i<Dimension; i++)
420 if (m_ArgsInfo.spacing_given==Dimension)
421 for (unsigned int i=0; i<Dimension; i++) {
422 spacing_low[i]=m_ArgsInfo.spacing_arg[i];
423 size_low[i]=ceil(input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i]/spacing_low[i]);
426 for (unsigned int i=0; i<Dimension; i++) {
427 spacing_low[i]=m_ArgsInfo.spacing_arg[0];
428 size_low[i]=ceil(input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i]/spacing_low[i]);
431 typename InternalImageType::PointType origin;
432 input->TransformIndexToPhysicalPoint(input->GetLargestPossibleRegion().GetIndex(), origin);
433 typedef itk::ResampleImageFilter<InternalImageType, InternalImageType> ResampleImageFilterType;
434 typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
435 typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
436 typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
437 resampler->SetInterpolator(interpolator);
438 resampler->SetOutputOrigin(origin);
439 resampler->SetOutputSpacing(spacing_low);
440 resampler->SetSize(size_low);
441 resampler->SetInput(input);
443 typename InternalImageType::Pointer output =resampler->GetOutput();
448 //-------------------------------------------------------------------
450 //-------------------------------------------------------------------
451 template <unsigned int Dimension, class PixelType>
452 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
453 MotionMaskGenericFilter::InitializeEllips( typename itk::Vector<double,Dimension> center, typename itk::Image<InternalPixelType,Dimension>::Pointer bones_low )
457 typedef int InternalPixelType;
458 typedef itk::Image<PixelType, Dimension> InputImageType;
459 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
460 typedef itk::Image< unsigned char , Dimension> OutputImageType;
461 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
462 typedef clitk::SetBackgroundImageFilter<InternalImageType,InternalImageType, InternalImageType> SetBackgroundFilterType;
463 typedef itk::LabelStatisticsImageFilter<InternalImageType,InternalImageType> LabelStatisticsImageFilterType;
464 typedef itk::CastImageFilter<InternalImageType,OutputImageType> CastImageFilterType;
467 typename InternalImageType::Pointer ellips;
468 if (m_ArgsInfo.ellips_given || m_ArgsInfo.grownEllips_given || m_ArgsInfo.filledRibs_given) {
469 if(m_ArgsInfo.ellips_given) {
470 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
471 typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
472 featureReader->SetFileName(m_ArgsInfo.ellips_arg);
473 featureReader->Update();
474 ellips=featureReader->GetOutput();
478 std::cout<<std::endl;
479 std::cout<<"=========================================="<<std::endl;
480 std::cout<<"|| Initializing ellipsoide ||"<<std::endl;
481 std::cout<<"=========================================="<<std::endl;
484 //---------------------------------
485 // Offset from center
486 //---------------------------------
487 typename itk::Vector<double,Dimension> offset;
488 if(m_ArgsInfo.offset_given== Dimension) {
489 for(unsigned int i=0; i<Dimension; i++)
490 offset[i]=m_ArgsInfo.offset_arg[i];
495 itk::Vector<double,Dimension> centerEllips=center+offset;
497 std::cout<<"Placing the center of the initial ellipsoide at (mm) "<<std::endl;
498 std::cout<<centerEllips[0];
499 for (unsigned int i=1; i<Dimension; i++)
500 std::cout<<" "<<centerEllips[i];
501 std::cout<<std::endl;
504 //---------------------------------
506 //---------------------------------
507 ellips=InternalImageType::New();
508 ellips->SetRegions (bones_low->GetLargestPossibleRegion());
509 ellips->SetOrigin(bones_low->GetOrigin());
510 ellips->SetSpacing(bones_low->GetSpacing());
512 ellips->FillBuffer(0);
515 typename itk::Vector<double, Dimension> axes;
516 if (m_ArgsInfo.axes_given==Dimension)
517 for(unsigned int i=0; i<Dimension; i++)
518 axes[i]=m_ArgsInfo.axes_arg[i];
526 IteratorType itEllips(ellips, ellips->GetLargestPossibleRegion());
527 itEllips.GoToBegin();
528 typename InternalImageType::PointType point;
529 typename InternalImageType::IndexType index;
532 if (m_Verbose) std::cout<<"Drawing the initial ellipsoide..."<<std::endl;
533 while (!itEllips.IsAtEnd()) {
534 index=itEllips.GetIndex();
535 ellips->TransformIndexToPhysicalPoint(index, point);
537 for(unsigned int i=0; i<Dimension; i++)
538 distance+=powf( ( (centerEllips[i]-point[i])/axes[i] ), 2);
547 //---------------------------------
549 //---------------------------------
550 if (m_ArgsInfo.writeEllips_given) {
551 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
552 caster->SetInput(ellips);
553 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeEllips_arg, m_Verbose);
561 //-------------------------------------------------------------------
562 // Update with the number of dimensions and the pixeltype
563 //-------------------------------------------------------------------
564 template <unsigned int Dimension, class PixelType>
566 MotionMaskGenericFilter::UpdateWithDimAndPixelType()
570 typedef int InternalPixelType;
571 typedef itk::Image<PixelType, Dimension> InputImageType;
572 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
573 typedef itk::Image< float, Dimension> LevelSetImageType; //labeling doesn't work with unsigned char?
574 typedef itk::Image<unsigned char, Dimension> OutputImageType;
577 //----------------------------------------------------------------------------------------------------
578 // Typedef for Processing
579 //----------------------------------------------------------------------------------------------------
580 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
581 typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
582 typedef itk::BinaryThresholdImageFilter< LevelSetImageType,InternalImageType> LevelSetBinarizeFilterType;
583 typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> InversionFilterType;
584 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdFilterType;
585 typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
586 typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
587 typedef itk::MirrorPadImageFilter<InternalImageType, InternalImageType> MirrorPadImageFilterType;
588 typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
589 typedef itk::ResampleImageFilter<InternalImageType, InternalImageType> ResampleImageFilterType;
590 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
591 typedef itk::GeodesicActiveContourLevelSetImageFilter<LevelSetImageType, InternalImageType> LevelSetImageFilterType;
592 typedef itk::SignedDanielssonDistanceMapImageFilter<InternalImageType, LevelSetImageType> DistanceMapImageFilterType;
593 typedef clitk::SetBackgroundImageFilter<InternalImageType,InternalImageType, InternalImageType> SetBackgroundFilterType;
594 typedef itk::LabelStatisticsImageFilter<InternalImageType,InternalImageType> LabelStatisticsImageFilterType;
595 typedef itk::CastImageFilter<InternalImageType,OutputImageType> CastImageFilterType;
599 typedef itk::ImageFileReader<InputImageType> InputReaderType;
600 typename InputReaderType::Pointer reader = InputReaderType::New();
601 reader->SetFileName( m_InputFileName);
603 typename InputImageType::Pointer input= reader->GetOutput();
605 // // globals for avi
606 // unsigned int number=0;
610 std::cout<<std::endl;
611 std::cout<<"=========================================="<<std::endl;
612 std::cout<<"|| Making feature images ||"<<std::endl;
613 std::cout<<"=========================================="<<std::endl;
616 //--------------------------------------------------------------------------------
618 //-------------------------------------------------------------------------------
619 typename InternalImageType::Pointer lungs = this->GetLungsImage<Dimension, PixelType>(input);
621 //-------------------------------------------------------------------------------
623 //-------------------------------------------------------------------------------
624 typename InternalImageType::Pointer air = this->GetAirImage<Dimension, PixelType>(input, lungs);
626 //-------------------------------------------------------------------------------
628 //-------------------------------------------------------------------------------
629 typename InternalImageType::Pointer bones = this->GetBonesImage<Dimension, PixelType>(input);
631 //----------------------------------------------------------------------------------------------------
632 // Write feature images
633 //----------------------------------------------------------------------------------------------------
634 if(m_ArgsInfo.writeFeature_given) {
635 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
636 setBackgroundFilter->SetInput(air);
637 setBackgroundFilter->SetInput2(bones);
638 setBackgroundFilter->SetMaskValue(0);
639 setBackgroundFilter->SetOutsideValue(2);
640 setBackgroundFilter->Update();
641 typename InternalImageType::Pointer bones_air =setBackgroundFilter->GetOutput();
643 typename SetBackgroundFilterType::Pointer setBackgroundFilter2 = SetBackgroundFilterType::New();
644 setBackgroundFilter2->SetInput(bones_air);
645 setBackgroundFilter2->SetInput2(lungs);
646 setBackgroundFilter2->SetMaskValue(0);
647 setBackgroundFilter2->SetOutsideValue(3);
648 setBackgroundFilter2->Update();
649 typename InternalImageType::Pointer lungs_bones_air =setBackgroundFilter2->GetOutput();
651 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
652 caster->SetInput(lungs_bones_air);
653 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeFeature_arg, m_Verbose);
656 //----------------------------------------------------------------------------------------------------
657 // Low dimensional versions
658 //----------------------------------------------------------------------------------------------------
659 typename InternalImageType::Pointer bones_low =Resample<Dimension , PixelType>(bones);
660 typename InternalImageType::Pointer lungs_low =Resample<Dimension , PixelType>(lungs);
661 typename InternalImageType::Pointer air_low =Resample<Dimension , PixelType>(air);
663 //---------------------------------
665 //---------------------------------
666 if(m_ArgsInfo.pad_flag) {
667 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
668 IteratorType it(air_low, air_low->GetLargestPossibleRegion());
669 typename InternalImageType::IndexType index;
670 while(!it.IsAtEnd()) {
672 for (unsigned int i=0; i<Dimension; i++)
673 if(index[i]==air_low->GetLargestPossibleRegion().GetIndex()[i]
674 || index[i]==(unsigned int)air_low->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air_low->GetLargestPossibleRegion().GetSize()[i]-1)
680 //---------------------------------
682 //---------------------------------
683 typename itk::Vector<double,Dimension> center;
684 typedef itk::ImageMomentsCalculator<InternalImageType> MomentsCalculatorType;
685 typename MomentsCalculatorType::Pointer momentsCalculator=MomentsCalculatorType::New();
686 momentsCalculator->SetImage(air);
687 momentsCalculator->Compute();
688 center=momentsCalculator->GetCenterOfGravity();
690 std::cout<<"The center of gravity of the patient body is located at (mm) "<<std::endl;
691 std::cout<<center[0];
692 for (unsigned int i=1; i<Dimension; i++)
693 std::cout<<" "<<center[i];
694 std::cout<<std::endl;
697 //----------------------------------------------------------------------------------------------------
698 // Ellipsoide initialization
699 //----------------------------------------------------------------------------------------------------
700 typename InternalImageType::Pointer m_Ellips= InitializeEllips<Dimension, PixelType>(center,bones_low);
702 //----------------------------------------------------------------------------------------------------
704 //----------------------------------------------------------------------------------------------------
705 typename InternalImageType::Pointer grownEllips;
706 if (m_ArgsInfo.grownEllips_given || m_ArgsInfo.filledRibs_given) {
707 if (m_ArgsInfo.grownEllips_given) {
709 typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
710 featureReader->SetFileName(m_ArgsInfo.grownEllips_arg);
711 featureReader->Update();
712 grownEllips=featureReader->GetOutput();
717 std::cout<<std::endl;
718 std::cout<<"=========================================="<<std::endl;
719 std::cout<<"|| Growing ellipsoide ||"<<std::endl;
720 std::cout<<"=========================================="<<std::endl;
723 //---------------------------------
725 //---------------------------------
726 typename InternalImageType::PointType dPoint;
727 if (m_ArgsInfo.detectionPoint_given) {
728 for (unsigned int i=0; i<Dimension; i++)
729 dPoint[i]=m_ArgsInfo.detectionPoint_arg[i];
731 typename InternalImageType::RegionType searchRegion=air->GetLargestPossibleRegion();
732 typename InternalImageType::SizeType searchRegionSize = searchRegion.GetSize();
733 typename InternalImageType::IndexType searchRegionIndex = searchRegion.GetIndex();
734 searchRegionIndex[2]+=searchRegionSize[2]/2;
735 searchRegionSize[2]=1;
736 searchRegion.SetSize(searchRegionSize);
737 searchRegion.SetIndex(searchRegionIndex);
738 IteratorType itAbdomen(air, searchRegion);
739 itAbdomen.GoToBegin();
741 typename InternalImageType::PointType aPoint;
742 typename InternalImageType::IndexType aIndex;
744 if (m_Verbose) std::cout<<"Detecting the abdomen..."<<std::endl;
745 while (!itAbdomen.IsAtEnd()) {
746 if(itAbdomen.Get()) break;
749 aIndex=itAbdomen.GetIndex();
750 air->TransformIndexToPhysicalPoint(aIndex,aPoint);
751 if (m_Verbose) std::cout<<"Detected the abdomen at "<<aPoint<<".."<<std::endl;
754 //---------------------------------
755 // Detect abdomen in additional images?
756 //---------------------------------
757 if (m_ArgsInfo.detectionPairs_given) {
758 for (unsigned int i=0; i<m_ArgsInfo.detectionPairs_given; i++) {
759 typename InternalImageType::Pointer airAdd;
760 //---------------------------------
762 //--------------------------------
763 typedef itk::ImageFileReader<InputImageType> InputReaderType;
764 typename InputReaderType::Pointer reader = InputReaderType::New();
765 reader->SetFileName( m_ArgsInfo.detectionPairs_arg[i]);
767 typename InputImageType::Pointer additional= reader->GetOutput();
768 if (m_Verbose) std::cout<<"Extracting the air from additional image "<< i<<"..."<<std::endl;
770 //---------------------------------
771 // Binarize the image
772 //---------------------------------
773 typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
774 binarizeFilter->SetInput(additional);
775 binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdAir_arg));
776 binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdAir_arg));
777 if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<<m_ArgsInfo.lowerThresholdAir_arg
778 <<", "<<m_ArgsInfo.upperThresholdAir_arg<<"..."<<std::endl;
780 //---------------------------------
781 // Label the connected components
782 //---------------------------------
783 typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
784 connectFilter->SetInput(binarizeFilter->GetOutput());
785 connectFilter->SetBackgroundValue(0);
786 connectFilter->SetFullyConnected(false);
787 if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
789 //---------------------------------
790 // Sort the labels according to size
791 //---------------------------------
792 typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
793 relabelFilter->SetInput(connectFilter->GetOutput());
794 if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
796 //---------------------------------
798 //---------------------------------
799 typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
800 thresholdFilter->SetInput(relabelFilter->GetOutput());
801 thresholdFilter->SetUpper(1);
802 if (m_Verbose) std::cout<<"Keeping the first label..."<<std::endl;
803 thresholdFilter->Update();
804 airAdd=thresholdFilter->GetOutput();
807 //---------------------------------
809 //---------------------------------
810 typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
811 inversionFilter->SetInput(airAdd);
812 inversionFilter->SetLowerThreshold(0);
813 inversionFilter->SetUpperThreshold(0);
814 inversionFilter->SetInsideValue(1);
815 inversionFilter->Update();
817 //---------------------------------
819 //---------------------------------
820 typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
821 mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
822 typename InternalImageType::SizeType padSize;
824 padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
825 mirrorPadImageFilter->SetPadLowerBound(padSize);
826 if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
827 mirrorPadImageFilter->Update();
828 airAdd=mirrorPadImageFilter->GetOutput();
830 //---------------------------------
832 //---------------------------------
833 typename InternalImageType::RegionType searchRegion=airAdd->GetLargestPossibleRegion();
834 typename InternalImageType::SizeType searchRegionSize = searchRegion.GetSize();
835 typename InternalImageType::IndexType searchRegionIndex = searchRegion.GetIndex();
836 searchRegionIndex[2]+=searchRegionSize[2]/2;
837 searchRegionSize[2]=1;
838 searchRegion.SetSize(searchRegionSize);
839 searchRegion.SetIndex(searchRegionIndex);
840 IteratorType itAbdomen(airAdd, searchRegion);
841 itAbdomen.GoToBegin();
843 typename InternalImageType::PointType additionalPoint;
844 typename InternalImageType::IndexType aIndex;
846 if (m_Verbose) std::cout<<"Detecting the abdomen..."<<std::endl;
847 while (!itAbdomen.IsAtEnd()) {
848 if(itAbdomen.Get()) break;
851 aIndex=itAbdomen.GetIndex();
852 airAdd->TransformIndexToPhysicalPoint(aIndex,additionalPoint);
853 if (m_Verbose) std::cout<<"Detected the abdomen in the additional image at "<<additionalPoint<<".."<<std::endl;
855 if(additionalPoint[1]< aPoint[1]) {
856 aPoint=additionalPoint;
857 if (m_Verbose) std::cout<<"Modifying the detected abdomen to "<<aPoint<<".."<<std::endl;
864 // Determine the detection point
868 if(m_ArgsInfo.offsetDetect_given==Dimension)
869 for(unsigned int i=0; i <Dimension; i++)
870 dPoint[i]+=m_ArgsInfo.offsetDetect_arg[i];
875 if (m_Verbose) std::cout<<"Setting the detection point to "<<dPoint<<".."<<std::endl;
878 //---------------------------------
879 // Pad the rib image and ellips image
880 //---------------------------------
881 typename InternalImageType::Pointer padded_ellips;
882 typename InternalImageType::Pointer padded_bones_low;
884 // If detection point not inside the image: pad
885 typename InternalImageType::IndexType dIndex;
886 if (!bones_low->TransformPhysicalPointToIndex(dPoint, dIndex)) {
887 typename InternalImageType::SizeType padSize;
889 padSize[1]=abs(dIndex[1])+1;
890 if (m_Verbose) std::cout<<"Padding the images with padding size "<<padSize<<" to include the detection point..."<<std::endl;
892 typename MirrorPadImageFilterType::Pointer padBonesFilter=MirrorPadImageFilterType::New();
893 padBonesFilter->SetInput(bones_low);
894 padBonesFilter->SetPadLowerBound(padSize);
895 padBonesFilter->Update();
896 padded_bones_low=padBonesFilter->GetOutput();
898 typename MirrorPadImageFilterType::Pointer padEllipsFilter=MirrorPadImageFilterType::New();
899 padEllipsFilter->SetInput(m_Ellips);
900 padEllipsFilter->SetPadLowerBound(padSize);
901 padEllipsFilter->Update();
902 padded_ellips=padEllipsFilter->GetOutput();
906 padded_bones_low=bones_low;
907 padded_ellips=m_Ellips;
911 //---------------------------------
912 // Calculate distance map
913 //---------------------------------
914 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
915 distanceMapImageFilter->SetInput(padded_ellips);
916 distanceMapImageFilter->SetInsideIsPositive(false);
917 if (m_Verbose) std::cout<<"Calculating the distance map..."<<std::endl;
918 distanceMapImageFilter->Update();
920 //---------------------------------
921 // Grow while monitoring detection point
922 //---------------------------------
923 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
924 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
925 levelSetFilter->SetFeatureImage( padded_bones_low );
926 levelSetFilter->SetPropagationScaling( 1 );
927 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve1_arg );
928 levelSetFilter->SetAdvectionScaling( 0 );
929 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS1_arg );
930 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter1_arg );
931 levelSetFilter->SetUseImageSpacing(true);
933 // //---------------------------------
934 // // Script for making movie
935 // //---------------------------------
936 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
939 if (m_Verbose) std::cout<<"Starting level set segmentation..."<<std::endl;
940 unsigned int totalNumberOfIterations=0;
942 levelSetFilter->Update();
943 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
945 if ( levelSetFilter->GetOutput()->GetPixel(dIndex) < 0 ) {
946 if (m_Verbose) std::cout<<"Detection point reached!"<<std::endl;
949 if (m_Verbose) std::cout<<"Detection point not reached after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
950 levelSetFilter->SetInput(levelSetFilter->GetOutput());
951 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose);
954 // std::ostringstream number_str;
955 // number_str << number;
956 // std::string param = number_str.str();
957 // system((script+param).c_str());
961 if ( (totalNumberOfIterations> 10000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS1_arg) ) break;
965 if (m_Verbose) std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
966 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
967 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
970 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
971 thresholder->SetUpperThreshold( 0.0 );
972 thresholder->SetOutsideValue( 0 );
973 thresholder->SetInsideValue( 1 );
974 thresholder->SetInput( levelSetFilter->GetOutput() );
975 if (m_Verbose) std::cout<<"Thresholding the output level set..."<<std::endl;
976 thresholder->Update();
977 grownEllips=thresholder->GetOutput();
980 //---------------------------------
981 // Write the grown ellips
982 //---------------------------------
983 if (m_ArgsInfo.writeGrownEllips_given) {
984 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
985 caster->SetInput(grownEllips);
986 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeGrownEllips_arg, m_Verbose);
990 //----------------------------------------------------------------------------------------------------
992 //----------------------------------------------------------------------------------------------------
993 typename InternalImageType::Pointer filledRibs;
994 if (m_ArgsInfo.filledRibs_given) {
995 typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
996 featureReader->SetFileName(m_ArgsInfo.filledRibs_arg);
997 if (m_Verbose) std::cout<<"Reading filled ribs mask..."<<std::endl;
998 featureReader->Update();
999 filledRibs=featureReader->GetOutput();
1002 std::cout<<std::endl;
1003 std::cout<<"=========================================="<<std::endl;
1004 std::cout<<"|| Filling the ribs image ||"<<std::endl;
1005 std::cout<<"=========================================="<<std::endl;
1007 //---------------------------------
1008 // Make feature image air+bones
1009 //---------------------------------
1010 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
1011 setBackgroundFilter->SetInput(air_low);
1012 setBackgroundFilter->SetInput2(bones_low);
1013 setBackgroundFilter->SetMaskValue(0);
1014 setBackgroundFilter->SetOutsideValue(0);
1015 setBackgroundFilter->Update();
1016 typename InternalImageType::Pointer bones_air_low =setBackgroundFilter->GetOutput();
1018 //---------------------------------
1019 // Resampling previous solution
1020 //---------------------------------
1021 if (m_Verbose) std::cout<<"Resampling previous solution..."<<std::endl;
1022 typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
1023 typename InternalImageType::PointType origin;
1024 bones_low->TransformIndexToPhysicalPoint(bones_low->GetLargestPossibleRegion().GetIndex(), origin);
1025 resampler->SetOutputOrigin(origin);
1026 resampler->SetOutputSpacing(bones_low->GetSpacing());
1027 typename InternalImageType::SizeType size_low= bones_low->GetLargestPossibleRegion().GetSize();
1028 resampler->SetSize(size_low);
1029 typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
1030 resampler->SetInterpolator(interpolator);
1031 resampler->SetInput(grownEllips);
1032 resampler->Update();
1033 typename InternalImageType::Pointer grownEllips =resampler->GetOutput();
1036 //---------------------------------
1037 // Calculate Distance Map
1038 //---------------------------------
1039 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
1040 distanceMapImageFilter->SetInput(grownEllips);
1041 distanceMapImageFilter->SetInsideIsPositive(false);
1042 if (m_Verbose) std::cout<<"Calculating distance map..."<<std::endl;
1043 distanceMapImageFilter->Update();
1045 //---------------------------------
1046 // Grow while monitoring lung volume coverage
1047 //---------------------------------
1048 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
1049 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
1050 levelSetFilter->SetFeatureImage( bones_air_low );
1051 levelSetFilter->SetPropagationScaling( 1 );
1052 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve2_arg );
1053 levelSetFilter->SetAdvectionScaling( 0 );
1054 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS2_arg );
1055 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter2_arg );
1056 levelSetFilter->SetUseImageSpacing(true);
1058 //---------------------------------
1060 //---------------------------------
1061 typename InversionFilterType::Pointer invertor= InversionFilterType::New();
1062 invertor->SetInput(lungs_low);
1063 invertor->SetLowerThreshold(0);
1064 invertor->SetUpperThreshold(0);
1065 invertor->SetInsideValue(1);
1067 typename InternalImageType::Pointer lungs_low_inv=invertor->GetOutput();
1069 //---------------------------------
1070 // Calculate lung volume
1071 //---------------------------------
1072 typename LabelStatisticsImageFilterType::Pointer statisticsFilter= LabelStatisticsImageFilterType::New();
1073 statisticsFilter->SetInput(lungs_low_inv);
1074 statisticsFilter->SetLabelInput(lungs_low);
1075 statisticsFilter->Update();
1076 unsigned int volume=statisticsFilter->GetSum(0);
1078 // Prepare threshold
1079 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1080 thresholder->SetUpperThreshold( 0.0 );
1081 thresholder->SetOutsideValue( 0 );
1082 thresholder->SetInsideValue( 1 );
1086 unsigned int totalNumberOfIterations=0;
1087 unsigned int coverage=0;
1090 // //---------------------------------
1091 // // Script for making movie
1092 // //---------------------------------
1093 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1096 levelSetFilter->Update();
1097 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1099 thresholder->SetInput( levelSetFilter->GetOutput() );
1100 thresholder->Update();
1101 statisticsFilter->SetInput(thresholder->GetOutput());
1102 statisticsFilter->Update();
1103 coverage=statisticsFilter->GetSum(0);
1105 // Compare the volumes
1106 if ( coverage >= m_ArgsInfo.fillingLevel_arg * volume/100 ) {
1107 if (m_Verbose) std::cout<<"Lungs filled for "<< m_ArgsInfo.fillingLevel_arg<<"% !"<<std::endl;
1110 if (m_Verbose) std::cout<<"After "<<totalNumberOfIterations<<" iterations, lungs are covered for "
1111 <<(double)coverage/volume*100.0<<"%..."<<std::endl;
1112 levelSetFilter->SetInput(levelSetFilter->GetOutput());
1113 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose);
1116 // std::ostringstream number_str;
1117 // number_str << number;
1118 // std::string param = number_str.str();
1119 // system((script+param).c_str());
1123 if ( (totalNumberOfIterations> 30000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS2_arg) ) break;
1127 if (m_Verbose) std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1128 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1129 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1132 thresholder->SetInput( levelSetFilter->GetOutput() );
1133 thresholder->Update();
1134 filledRibs=thresholder->GetOutput();
1135 // writeImage<InternalImageType>(filledRibs,"/home/jef/tmp/filled_ribs_before.mhd", m_Verbose);
1139 //---------------------------------
1140 // Write the filled ribs
1141 //---------------------------------
1142 if (m_ArgsInfo.writeFilledRibs_given) {
1143 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1144 caster->SetInput(filledRibs);
1145 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeFilledRibs_arg, m_Verbose);
1149 //----------------------------------------------------------------------------------------------------
1150 // Collapse to the lungs
1151 //----------------------------------------------------------------------------------------------------
1153 std::cout<<std::endl;
1154 std::cout<<"=========================================="<<std::endl;
1155 std::cout<<"|| Collapsing to the lung image ||"<<std::endl;
1156 std::cout<<"=========================================="<<std::endl;
1159 //---------------------------------
1160 // Make feature image air+bones
1161 //---------------------------------
1162 if (m_Verbose) std::cout<<"Making feature images..."<<std::endl;
1163 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
1164 setBackgroundFilter->SetInput(air);
1165 setBackgroundFilter->SetInput2(bones);
1166 setBackgroundFilter->SetMaskValue(0);
1167 setBackgroundFilter->SetOutsideValue(0);
1168 setBackgroundFilter->Update();
1169 typename InternalImageType::Pointer bones_air =setBackgroundFilter->GetOutput();
1171 typename SetBackgroundFilterType::Pointer setBackgroundFilter2 = SetBackgroundFilterType::New();
1172 setBackgroundFilter2->SetInput(bones_air);
1173 setBackgroundFilter2->SetInput2(lungs);
1174 setBackgroundFilter2->SetMaskValue(0);
1175 setBackgroundFilter2->SetOutsideValue(0);
1176 setBackgroundFilter2->Update();
1177 typename InternalImageType::Pointer lungs_bones_air =setBackgroundFilter2->GetOutput();
1179 //---------------------------------
1180 // Prepare previous solution
1181 //---------------------------------
1182 if (m_Verbose) std::cout<<"Resampling previous solution..."<<std::endl;
1183 typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
1184 typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
1185 typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
1186 resampler->SetOutputStartIndex(bones->GetLargestPossibleRegion().GetIndex());
1187 resampler->SetInput(filledRibs);
1188 resampler->SetInterpolator(interpolator);
1189 resampler->SetOutputParametersFromImage(bones);
1190 resampler->Update();
1191 filledRibs =resampler->GetOutput();
1193 if (m_Verbose) std::cout<<"Setting lungs to 1..."<<std::endl;
1194 typename SetBackgroundFilterType::Pointer setBackgroundFilter3 = SetBackgroundFilterType::New();
1195 setBackgroundFilter3->SetInput(filledRibs);
1196 setBackgroundFilter3->SetInput2(lungs);
1197 setBackgroundFilter3->SetMaskValue(0);
1198 setBackgroundFilter3->SetOutsideValue(1);
1199 setBackgroundFilter3->Update();
1200 filledRibs=setBackgroundFilter3->GetOutput();
1202 if (m_Verbose) std::cout<<"Removing overlap with bones..."<<std::endl;
1203 typename SetBackgroundFilterType::Pointer setBackgroundFilter4 = SetBackgroundFilterType::New();
1204 setBackgroundFilter4->SetInput(filledRibs);
1205 setBackgroundFilter4->SetInput2(bones);
1206 setBackgroundFilter4->SetMaskValue(0);
1207 setBackgroundFilter4->SetOutsideValue(0);
1208 setBackgroundFilter4->Update();
1209 filledRibs =setBackgroundFilter4->GetOutput();
1210 // writeImage<InternalImageType>(filledRibs,"/home/jef/tmp/filledRibs_pp.mhd");
1211 //---------------------------------
1212 // Calculate Distance Map
1213 //---------------------------------
1214 // typedef itk::ApproximateSignedDistanceMapImageFilter <InternalImageType, LevelSetImageType> DistanceMapImageFilterType2;
1215 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
1216 distanceMapImageFilter->SetInput(filledRibs);
1217 distanceMapImageFilter->SetInsideIsPositive(false);
1218 // distanceMapImageFilter->SetInsideValue(0);
1219 // distanceMapImageFilter->SetOutsideValue(1);
1220 if (m_Verbose) std::cout<<"Calculating distance map..."<<std::endl;
1221 distanceMapImageFilter->Update();
1223 //---------------------------------
1225 //---------------------------------
1226 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
1227 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
1228 levelSetFilter->SetFeatureImage( lungs_bones_air );
1229 levelSetFilter->SetPropagationScaling(m_ArgsInfo.prop3_arg);
1230 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve3_arg );
1231 levelSetFilter->SetAdvectionScaling( 0 );
1232 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS3_arg );
1233 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter3_arg );
1234 levelSetFilter->SetUseImageSpacing(true);
1237 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1240 if (m_Verbose) std::cout<<"Starting the levelset..."<<std::endl;
1241 unsigned int totalNumberOfIterations=0;
1243 levelSetFilter->Update();
1246 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1247 levelSetFilter->SetInput(levelSetFilter->GetOutput());
1248 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg);
1249 std::cout << "After "<<totalNumberOfIterations<<" iteration the RMS change is " << levelSetFilter->GetRMSChange() <<"..."<< std::endl;
1252 // std::ostringstream number_str;
1253 // number_str << number;
1254 // std::string param = number_str.str();
1255 // system((script+param).c_str());
1258 // stopping condition
1259 if ( (totalNumberOfIterations>= (unsigned int) m_ArgsInfo.maxIter3_arg) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS3_arg) ) break;
1260 levelSetFilter->SetNumberOfIterations( std::min ((unsigned int) m_ArgsInfo.iter3_arg, (unsigned int) m_ArgsInfo.maxIter3_arg-totalNumberOfIterations ) );
1265 std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1266 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1267 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1271 typedef itk::BinaryThresholdImageFilter< LevelSetImageType,InternalImageType> LevelSetBinarizeFilterType;
1272 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1273 thresholder->SetUpperThreshold( 0.0 );
1274 thresholder->SetOutsideValue( 0 );
1275 thresholder->SetInsideValue( 1 );
1276 thresholder->SetInput( levelSetFilter->GetOutput() );
1277 thresholder->Update();
1278 typename InternalImageType::Pointer output = thresholder->GetOutput();
1281 //----------------------------------------------------------------------------------------------------
1282 // Prepare the output
1283 //----------------------------------------------------------------------------------------------------
1285 //---------------------------------
1287 //---------------------------------
1288 if (m_Verbose) std::cout<<"Removing overlap mask with air..."<<std::endl;
1289 typename SetBackgroundFilterType::Pointer setBackgroundFilter5 = SetBackgroundFilterType::New();
1290 setBackgroundFilter5->SetInput(output);
1291 setBackgroundFilter5->SetInput2(air);
1292 setBackgroundFilter5->SetMaskValue(0);
1293 setBackgroundFilter5->SetOutsideValue(0);
1294 setBackgroundFilter5->Update();
1295 output=setBackgroundFilter5->GetOutput();
1297 //---------------------------------
1298 // Open & close to cleanup
1299 //---------------------------------
1300 if ( m_ArgsInfo.openClose_flag) {
1302 //---------------------------------
1303 // Structuring element
1304 //---------------------------------
1305 typedef itk::BinaryBallStructuringElement<InternalPixelType,Dimension > KernelType;
1306 KernelType structuringElement;
1307 structuringElement.SetRadius(1);
1308 structuringElement.CreateStructuringElement();
1310 //---------------------------------
1312 //---------------------------------
1313 typedef itk::BinaryMorphologicalOpeningImageFilter<InternalImageType, InternalImageType , KernelType> OpenFilterType;
1314 typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
1315 openFilter->SetInput(output);
1316 openFilter->SetBackgroundValue(0);
1317 openFilter->SetForegroundValue(1);
1318 openFilter->SetKernel(structuringElement);
1319 if(m_Verbose) std::cout<<"Opening the output image..."<<std::endl;
1321 //---------------------------------
1323 //---------------------------------
1324 typedef itk::BinaryMorphologicalClosingImageFilter<InternalImageType, InternalImageType , KernelType> CloseFilterType;
1325 typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
1326 closeFilter->SetInput(openFilter->GetOutput());
1327 closeFilter->SetSafeBorder(true);
1328 closeFilter->SetForegroundValue(1);
1329 closeFilter->SetKernel(structuringElement);
1330 if(m_Verbose) std::cout<<"Closing the output image..."<<std::endl;
1331 closeFilter->Update();
1332 output=closeFilter->GetOutput();
1335 //writeImage<InternalImageType>(output,"/home/jef/tmp/mm_double.mhd");
1337 // Extract the upper part
1338 typedef itk::CropImageFilter<InternalImageType, InternalImageType> CropImageFilterType;
1339 typename CropImageFilterType::Pointer cropFilter=CropImageFilterType::New();
1340 cropFilter->SetInput(output);
1341 typename InputImageType::SizeType lSize, uSize;
1344 lSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
1345 cropFilter->SetLowerBoundaryCropSize(lSize);
1346 cropFilter->SetUpperBoundaryCropSize(uSize);
1347 cropFilter->Update();
1350 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1351 caster->SetInput(cropFilter->GetOutput());
1352 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.output_arg, m_Verbose);
1358 #endif //#define clitkMotionMaskGenericFilter_txx