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 )
76 typedef int InternalPixelType;
77 typedef itk::Image<PixelType, Dimension> InputImageType;
78 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
80 //----------------------------------------------------------------------------------------------------
81 // Typedef for Processing
82 //----------------------------------------------------------------------------------------------------
83 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
84 typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
85 typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> InversionFilterType;
86 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdFilterType;
87 typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
88 typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
89 typedef itk::MirrorPadImageFilter<InternalImageType, InternalImageType> MirrorPadImageFilterType;
92 typename InternalImageType::Pointer air;
93 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();
103 if (m_Verbose) std::cout<<"Extracting the air feature image..."<<std::endl;
104 //---------------------------------
105 // Binarize the image
106 //---------------------------------
107 typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
108 binarizeFilter->SetInput(input);
109 binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdAir_arg));
110 binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdAir_arg));
111 if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<<m_ArgsInfo.lowerThresholdAir_arg
112 <<", "<<m_ArgsInfo.upperThresholdAir_arg<<"..."<<std::endl;
114 //---------------------------------
115 // Label the connected components
116 //---------------------------------
117 typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
118 connectFilter->SetInput(binarizeFilter->GetOutput());
119 connectFilter->SetBackgroundValue(0);
120 connectFilter->SetFullyConnected(false);
121 if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
123 //---------------------------------
124 // Sort the labels according to size
125 //---------------------------------
126 typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
127 relabelFilter->SetInput(connectFilter->GetOutput());
128 if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
130 //---------------------------------
132 //---------------------------------
133 typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
134 thresholdFilter->SetInput(relabelFilter->GetOutput());
135 thresholdFilter->SetUpper(1);
136 if (m_Verbose) std::cout<<"Keeping the first label..."<<std::endl;
137 thresholdFilter->Update();
138 air=thresholdFilter->GetOutput();
141 //---------------------------------
143 //---------------------------------
144 typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
145 inversionFilter->SetInput(air);
146 inversionFilter->SetLowerThreshold(0);
147 inversionFilter->SetUpperThreshold(0);
148 inversionFilter->SetInsideValue(1);
149 inversionFilter->Update();
151 //---------------------------------
153 //---------------------------------
154 typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
155 mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
156 typename InternalImageType::SizeType padSize;
158 padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
159 mirrorPadImageFilter->SetPadLowerBound(padSize);
160 if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
161 mirrorPadImageFilter->Update();
162 air=mirrorPadImageFilter->GetOutput();
163 // writeImage<InternalImageType>(air,"/home/jef/tmp/air.mhd");
165 //---------------------------------
167 //---------------------------------
168 if(m_ArgsInfo.pad_flag)
170 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
171 IteratorType it(air, air->GetLargestPossibleRegion());
172 typename InternalImageType::IndexType index;
176 for (unsigned int i=0;i<Dimension;i++)
177 if(index[i]==air->GetLargestPossibleRegion().GetIndex()[i]
178 || index[i]==(unsigned int)air->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air->GetLargestPossibleRegion().GetSize()[i]-1)
188 //-------------------------------------------------------------------
190 //-------------------------------------------------------------------
191 template <unsigned int Dimension, class PixelType>
192 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
193 MotionMaskGenericFilter::GetBonesImage(typename itk::Image<PixelType, Dimension>::Pointer input )
197 typedef int InternalPixelType;
198 typedef itk::Image<PixelType, Dimension> InputImageType;
199 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
201 //----------------------------------------------------------------------------------------------------
202 // Typedef for Processing
203 //----------------------------------------------------------------------------------------------------
204 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
205 typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
206 typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> InversionFilterType;
207 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdFilterType;
208 typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
209 typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
210 typedef itk::MirrorPadImageFilter<InternalImageType, InternalImageType> MirrorPadImageFilterType;
213 typename InternalImageType::Pointer bones;
214 if (m_ArgsInfo.featureBones_given)
216 typename FeatureReaderType::Pointer featureReader =FeatureReaderType::New();
217 featureReader->SetFileName(m_ArgsInfo.featureBones_arg);
218 if (m_Verbose) std::cout<<"Reading the bones feature image..."<<std::endl;
219 featureReader->Update();
220 bones=featureReader->GetOutput();
224 if (m_Verbose) std::cout<<"Extracting the bones feature image..."<<std::endl;
225 //---------------------------------
226 // Binarize the image
227 //---------------------------------
228 typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
229 binarizeFilter->SetInput(input);
230 binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdBones_arg));
231 binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdBones_arg));
232 if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<<m_ArgsInfo.lowerThresholdBones_arg
233 <<", "<<m_ArgsInfo.upperThresholdBones_arg<<"..."<<std::endl;
235 //---------------------------------
236 // Label the connected components
237 //---------------------------------
238 typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
239 connectFilter->SetInput(binarizeFilter->GetOutput());
240 connectFilter->SetBackgroundValue(0);
241 connectFilter->SetFullyConnected(false);
242 if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
244 //---------------------------------
245 // Sort the labels according to size
246 //---------------------------------
247 typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
248 relabelFilter->SetInput(connectFilter->GetOutput());
249 if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
251 //---------------------------------
253 //---------------------------------
254 typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
255 thresholdFilter->SetInput(relabelFilter->GetOutput());
256 thresholdFilter->SetUpper(1);
257 if (m_Verbose) std::cout<<"Keeping the first label..."<<std::endl;
258 thresholdFilter->Update();
259 bones=thresholdFilter->GetOutput();
263 //---------------------------------
265 //---------------------------------
266 typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
267 inversionFilter->SetInput(bones);
268 inversionFilter->SetLowerThreshold(0);
269 inversionFilter->SetUpperThreshold(0);
270 inversionFilter->SetInsideValue(1);
271 inversionFilter->Update();
273 //---------------------------------
275 //---------------------------------
276 typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
277 mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
278 typename InternalImageType::SizeType padSize;
280 padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
281 mirrorPadImageFilter->SetPadLowerBound(padSize);
282 if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
283 mirrorPadImageFilter->Update();
284 bones=mirrorPadImageFilter->GetOutput();
285 // writeImage<InternalImageType>(bones,"/home/jef/tmp/bones.mhd");
293 //-------------------------------------------------------------------
295 //-------------------------------------------------------------------
296 template <unsigned int Dimension, class PixelType>
297 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
298 MotionMaskGenericFilter::GetLungsImage(typename itk::Image<PixelType, Dimension>::Pointer input )
301 typedef int InternalPixelType;
302 typedef itk::Image<PixelType, Dimension> InputImageType;
303 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
305 //----------------------------------------------------------------------------------------------------
306 // Typedef for Processing
307 //----------------------------------------------------------------------------------------------------
308 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
309 typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
310 typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> InversionFilterType;
311 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdFilterType;
312 typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
313 typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
314 typedef itk::MirrorPadImageFilter<InternalImageType, InternalImageType> MirrorPadImageFilterType;
316 typename InternalImageType::Pointer lungs;
317 if (m_ArgsInfo.featureLungs_given)
319 typename FeatureReaderType::Pointer featureReader =FeatureReaderType::New();
320 featureReader->SetFileName(m_ArgsInfo.featureLungs_arg);
321 if (m_Verbose) std::cout<<"Reading the lungs feature image..."<<std::endl;
322 featureReader->Update();
323 lungs=featureReader->GetOutput();
327 if (m_Verbose) std::cout<<"Extracting the lungs feature image..."<<std::endl;
328 //---------------------------------
329 // Binarize the image
330 //---------------------------------
331 typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
332 binarizeFilter->SetInput(input);
333 binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdLungs_arg));
334 binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdLungs_arg));
335 if (m_Verbose) std::cout<<"Binarizing the image using a threshold "<<m_ArgsInfo.lowerThresholdLungs_arg
336 <<", "<<m_ArgsInfo.upperThresholdLungs_arg<<"..."<<std::endl;
338 //---------------------------------
339 // Label the connected components
340 //---------------------------------
341 typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
342 connectFilter->SetInput(binarizeFilter->GetOutput());
343 connectFilter->SetBackgroundValue(0);
344 connectFilter->SetFullyConnected(true);
345 if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
347 //---------------------------------
348 // Sort the labels according to size
349 //---------------------------------
350 typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
351 relabelFilter->SetInput(connectFilter->GetOutput());
352 if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
353 // writeImage<InternalImageType> (relabelFilter->GetOutput(), "/home/jef/tmp/labels.mhd");
355 //---------------------------------
357 //---------------------------------
358 typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
359 thresholdFilter->SetInput(relabelFilter->GetOutput());
360 thresholdFilter->SetLower(1);
361 thresholdFilter->SetUpper(2);
362 if (m_Verbose) std::cout<<"Keeping the first two labels..."<<std::endl;
363 thresholdFilter->Update();
364 lungs=thresholdFilter->GetOutput();
369 //---------------------------------
371 //---------------------------------
372 typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
373 inversionFilter->SetInput(lungs);
374 inversionFilter->SetLowerThreshold(0);
375 inversionFilter->SetUpperThreshold(0);
376 inversionFilter->SetInsideValue(1);
377 inversionFilter->Update();
379 //---------------------------------
381 //---------------------------------
382 typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
383 mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
384 typename InternalImageType::SizeType padSize;
386 padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
387 mirrorPadImageFilter->SetPadLowerBound(padSize);
388 if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
389 mirrorPadImageFilter->Update();
390 lungs=mirrorPadImageFilter->GetOutput();
391 // writeImage<InternalImageType>(lungs,"/home/jef/tmp/lungs.mhd");
396 //-------------------------------------------------------------------
398 //-------------------------------------------------------------------
399 template <unsigned int Dimension, class PixelType>
400 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
401 MotionMaskGenericFilter::Resample( typename itk::Image<InternalPixelType,Dimension>::Pointer input )
404 typedef int InternalPixelType;
405 typedef itk::Image<PixelType, Dimension> InputImageType;
406 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
408 //---------------------------------
409 // Resample to new spacing
410 //---------------------------------
411 typename InternalImageType::SizeType size_low;
412 typename InternalImageType::SpacingType spacing_low;
413 for (unsigned int i=0; i<Dimension;i++)
414 if (m_ArgsInfo.spacing_given==Dimension)
415 for (unsigned int i=0; i<Dimension;i++)
417 spacing_low[i]=m_ArgsInfo.spacing_arg[i];
418 size_low[i]=ceil(input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i]/spacing_low[i]);
421 for (unsigned int i=0; i<Dimension;i++)
423 spacing_low[i]=m_ArgsInfo.spacing_arg[0];
424 size_low[i]=ceil(input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i]/spacing_low[i]);
427 typename InternalImageType::PointType origin;
428 input->TransformIndexToPhysicalPoint(input->GetLargestPossibleRegion().GetIndex(), origin);
429 typedef itk::ResampleImageFilter<InternalImageType, InternalImageType> ResampleImageFilterType;
430 typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
431 typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
432 typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
433 resampler->SetInterpolator(interpolator);
434 resampler->SetOutputOrigin(origin);
435 resampler->SetOutputSpacing(spacing_low);
436 resampler->SetSize(size_low);
437 resampler->SetInput(input);
439 typename InternalImageType::Pointer output =resampler->GetOutput();
444 //-------------------------------------------------------------------
446 //-------------------------------------------------------------------
447 template <unsigned int Dimension, class PixelType>
448 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
449 MotionMaskGenericFilter::InitializeEllips( typename itk::Vector<double,Dimension> center, typename itk::Image<InternalPixelType,Dimension>::Pointer bones_low )
453 typedef int InternalPixelType;
454 typedef itk::Image<PixelType, Dimension> InputImageType;
455 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
456 typedef itk::Image< unsigned char , Dimension> OutputImageType;
457 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
458 typedef clitk::SetBackgroundImageFilter<InternalImageType,InternalImageType, InternalImageType> SetBackgroundFilterType;
459 typedef itk::LabelStatisticsImageFilter<InternalImageType,InternalImageType> LabelStatisticsImageFilterType;
460 typedef itk::CastImageFilter<InternalImageType,OutputImageType> CastImageFilterType;
463 typename InternalImageType::Pointer ellips;
464 if (m_ArgsInfo.ellips_given || m_ArgsInfo.grownEllips_given || m_ArgsInfo.filledRibs_given)
466 if(m_ArgsInfo.ellips_given)
468 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
469 typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
470 featureReader->SetFileName(m_ArgsInfo.ellips_arg);
471 featureReader->Update();
472 ellips=featureReader->GetOutput();
479 std::cout<<std::endl;
480 std::cout<<"=========================================="<<std::endl;
481 std::cout<<"|| Initializing ellipsoide ||"<<std::endl;
482 std::cout<<"=========================================="<<std::endl;
485 //---------------------------------
486 // Offset from center
487 //---------------------------------
488 typename itk::Vector<double,Dimension> offset;
489 if(m_ArgsInfo.offset_given== Dimension)
491 for(unsigned int i=0; i<Dimension; i++)
492 offset[i]=m_ArgsInfo.offset_arg[i];
499 itk::Vector<double,Dimension> centerEllips=center+offset;
502 std::cout<<"Placing the center of the initial ellipsoide at (mm) "<<std::endl;
503 std::cout<<centerEllips[0];
504 for (unsigned int i=1; i<Dimension;i++)
505 std::cout<<" "<<centerEllips[i];
506 std::cout<<std::endl;
509 //---------------------------------
511 //---------------------------------
512 ellips=InternalImageType::New();
513 ellips->SetRegions (bones_low->GetLargestPossibleRegion());
514 ellips->SetOrigin(bones_low->GetOrigin());
515 ellips->SetSpacing(bones_low->GetSpacing());
517 ellips->FillBuffer(0);
520 typename itk::Vector<double, Dimension> axes;
521 if (m_ArgsInfo.axes_given==Dimension)
522 for(unsigned int i=0; i<Dimension; i++)
523 axes[i]=m_ArgsInfo.axes_arg[i];
532 IteratorType itEllips(ellips, ellips->GetLargestPossibleRegion());
533 itEllips.GoToBegin();
534 typename InternalImageType::PointType point;
535 typename InternalImageType::IndexType index;
538 if (m_Verbose) std::cout<<"Drawing the initial ellipsoide..."<<std::endl;
539 while (!itEllips.IsAtEnd())
541 index=itEllips.GetIndex();
542 ellips->TransformIndexToPhysicalPoint(index, point);
544 for(unsigned int i=0; i<Dimension; i++)
545 distance+=powf( ( (centerEllips[i]-point[i])/axes[i] ), 2);
554 //---------------------------------
556 //---------------------------------
557 if (m_ArgsInfo.writeEllips_given)
559 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
560 caster->SetInput(ellips);
561 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeEllips_arg, m_Verbose);
569 //-------------------------------------------------------------------
570 // Update with the number of dimensions and the pixeltype
571 //-------------------------------------------------------------------
572 template <unsigned int Dimension, class PixelType>
574 MotionMaskGenericFilter::UpdateWithDimAndPixelType()
578 typedef int InternalPixelType;
579 typedef itk::Image<PixelType, Dimension> InputImageType;
580 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
581 typedef itk::Image< float, Dimension> LevelSetImageType; //labeling doesn't work with unsigned char?
582 typedef itk::Image<unsigned char, Dimension> OutputImageType;
585 //----------------------------------------------------------------------------------------------------
586 // Typedef for Processing
587 //----------------------------------------------------------------------------------------------------
588 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
589 typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
590 typedef itk::BinaryThresholdImageFilter< LevelSetImageType,InternalImageType> LevelSetBinarizeFilterType;
591 typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> InversionFilterType;
592 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdFilterType;
593 typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
594 typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
595 typedef itk::MirrorPadImageFilter<InternalImageType, InternalImageType> MirrorPadImageFilterType;
596 typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
597 typedef itk::ResampleImageFilter<InternalImageType, InternalImageType> ResampleImageFilterType;
598 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
599 typedef itk::GeodesicActiveContourLevelSetImageFilter<LevelSetImageType, InternalImageType> LevelSetImageFilterType;
600 typedef itk::SignedDanielssonDistanceMapImageFilter<InternalImageType, LevelSetImageType> DistanceMapImageFilterType;
601 typedef clitk::SetBackgroundImageFilter<InternalImageType,InternalImageType, InternalImageType> SetBackgroundFilterType;
602 typedef itk::LabelStatisticsImageFilter<InternalImageType,InternalImageType> LabelStatisticsImageFilterType;
603 typedef itk::CastImageFilter<InternalImageType,OutputImageType> CastImageFilterType;
607 typedef itk::ImageFileReader<InputImageType> InputReaderType;
608 typename InputReaderType::Pointer reader = InputReaderType::New();
609 reader->SetFileName( m_InputFileName);
611 typename InputImageType::Pointer input= reader->GetOutput();
613 // // globals for avi
614 // unsigned int number=0;
619 std::cout<<std::endl;
620 std::cout<<"=========================================="<<std::endl;
621 std::cout<<"|| Making feaure images ||"<<std::endl;
622 std::cout<<"=========================================="<<std::endl;
625 //-------------------------------------------------------------------------------
627 //-------------------------------------------------------------------------------
628 typename InternalImageType::Pointer air = this->GetAirImage<Dimension, PixelType>(input);
630 //-------------------------------------------------------------------------------
632 //-------------------------------------------------------------------------------
633 typename InternalImageType::Pointer bones = this->GetBonesImage<Dimension, PixelType>(input);
635 //--------------------------------------------------------------------------------
637 //-------------------------------------------------------------------------------
638 typename InternalImageType::Pointer lungs = this->GetLungsImage<Dimension, PixelType>(input);
640 //----------------------------------------------------------------------------------------------------
641 // Write feature images
642 //----------------------------------------------------------------------------------------------------
643 if(m_ArgsInfo.writeFeature_given)
645 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
646 setBackgroundFilter->SetInput(air);
647 setBackgroundFilter->SetInput2(bones);
648 setBackgroundFilter->SetMaskValue(0);
649 setBackgroundFilter->SetOutsideValue(2);
650 setBackgroundFilter->Update();
651 typename InternalImageType::Pointer bones_air =setBackgroundFilter->GetOutput();
653 typename SetBackgroundFilterType::Pointer setBackgroundFilter2 = SetBackgroundFilterType::New();
654 setBackgroundFilter2->SetInput(bones_air);
655 setBackgroundFilter2->SetInput2(lungs);
656 setBackgroundFilter2->SetMaskValue(0);
657 setBackgroundFilter2->SetOutsideValue(3);
658 setBackgroundFilter2->Update();
659 typename InternalImageType::Pointer lungs_bones_air =setBackgroundFilter2->GetOutput();
661 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
662 caster->SetInput(lungs_bones_air);
663 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeFeature_arg, m_Verbose);
666 //----------------------------------------------------------------------------------------------------
667 // Low dimensional versions
668 //----------------------------------------------------------------------------------------------------
669 typename InternalImageType::Pointer bones_low =Resample<Dimension , PixelType>(bones);
670 typename InternalImageType::Pointer lungs_low =Resample<Dimension , PixelType>(lungs);
671 typename InternalImageType::Pointer air_low =Resample<Dimension , PixelType>(air);
673 //---------------------------------
675 //---------------------------------
676 if(m_ArgsInfo.pad_flag)
678 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
679 IteratorType it(air_low, air_low->GetLargestPossibleRegion());
680 typename InternalImageType::IndexType index;
684 for (unsigned int i=0;i<Dimension;i++)
685 if(index[i]==air_low->GetLargestPossibleRegion().GetIndex()[i]
686 || index[i]==(unsigned int)air_low->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air_low->GetLargestPossibleRegion().GetSize()[i]-1)
692 //---------------------------------
694 //---------------------------------
695 typename itk::Vector<double,Dimension> center;
696 typedef itk::ImageMomentsCalculator<InternalImageType> MomentsCalculatorType;
697 typename MomentsCalculatorType::Pointer momentsCalculator=MomentsCalculatorType::New();
698 momentsCalculator->SetImage(air);
699 momentsCalculator->Compute();
700 center=momentsCalculator->GetCenterOfGravity();
703 std::cout<<"The center of gravity of the patient body is located at (mm) "<<std::endl;
704 std::cout<<center[0];
705 for (unsigned int i=1; i<Dimension;i++)
706 std::cout<<" "<<center[i];
707 std::cout<<std::endl;
710 //----------------------------------------------------------------------------------------------------
711 // Ellipsoide initialization
712 //----------------------------------------------------------------------------------------------------
713 typename InternalImageType::Pointer m_Ellips= InitializeEllips<Dimension, PixelType>(center,bones_low);
715 //----------------------------------------------------------------------------------------------------
717 //----------------------------------------------------------------------------------------------------
718 typename InternalImageType::Pointer grownEllips;
719 if (m_ArgsInfo.grownEllips_given || m_ArgsInfo.filledRibs_given)
721 if (m_ArgsInfo.grownEllips_given)
724 typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
725 featureReader->SetFileName(m_ArgsInfo.grownEllips_arg);
726 featureReader->Update();
727 grownEllips=featureReader->GetOutput();
735 std::cout<<std::endl;
736 std::cout<<"=========================================="<<std::endl;
737 std::cout<<"|| Growing ellipsoide ||"<<std::endl;
738 std::cout<<"=========================================="<<std::endl;
741 //---------------------------------
743 //---------------------------------
744 typename InternalImageType::PointType dPoint;
745 if (m_ArgsInfo.detectionPoint_given)
747 for (unsigned int i=0;i<Dimension;i++)
748 dPoint[i]=m_ArgsInfo.detectionPoint_arg[i];
752 typename InternalImageType::RegionType searchRegion=air->GetLargestPossibleRegion();
753 typename InternalImageType::SizeType searchRegionSize = searchRegion.GetSize();
754 typename InternalImageType::IndexType searchRegionIndex = searchRegion.GetIndex();
755 searchRegionIndex[2]+=searchRegionSize[2]/2;
756 searchRegionSize[2]=1;
757 searchRegion.SetSize(searchRegionSize);
758 searchRegion.SetIndex(searchRegionIndex);
759 IteratorType itAbdomen(air, searchRegion);
760 itAbdomen.GoToBegin();
762 typename InternalImageType::PointType aPoint;
763 typename InternalImageType::IndexType aIndex;
765 if (m_Verbose) std::cout<<"Detecting the abdomen..."<<std::endl;
766 while (!itAbdomen.IsAtEnd())
768 if(itAbdomen.Get()) break;
771 aIndex=itAbdomen.GetIndex();
772 air->TransformIndexToPhysicalPoint(aIndex,aPoint);
773 if (m_Verbose) std::cout<<"Detected the abdomen at "<<aPoint<<".."<<std::endl;
776 //---------------------------------
777 // Detect abdomen in additional images?
778 //---------------------------------
779 if (m_ArgsInfo.detectionPairs_given)
781 for (unsigned int i=0; i<m_ArgsInfo.detectionPairs_given; i++)
783 typename InternalImageType::Pointer airAdd;
784 //---------------------------------
786 //--------------------------------
787 typedef itk::ImageFileReader<InputImageType> InputReaderType;
788 typename InputReaderType::Pointer reader = InputReaderType::New();
789 reader->SetFileName( m_ArgsInfo.detectionPairs_arg[i]);
791 typename InputImageType::Pointer additional= reader->GetOutput();
792 if (m_Verbose) std::cout<<"Extracting the air from additional image "<< i<<"..."<<std::endl;
794 //---------------------------------
795 // Binarize the image
796 //---------------------------------
797 typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
798 binarizeFilter->SetInput(additional);
799 binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdAir_arg));
800 binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdAir_arg));
801 if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<<m_ArgsInfo.lowerThresholdAir_arg
802 <<", "<<m_ArgsInfo.upperThresholdAir_arg<<"..."<<std::endl;
804 //---------------------------------
805 // Label the connected components
806 //---------------------------------
807 typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
808 connectFilter->SetInput(binarizeFilter->GetOutput());
809 connectFilter->SetBackgroundValue(0);
810 connectFilter->SetFullyConnected(false);
811 if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
813 //---------------------------------
814 // Sort the labels according to size
815 //---------------------------------
816 typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
817 relabelFilter->SetInput(connectFilter->GetOutput());
818 if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
820 //---------------------------------
822 //---------------------------------
823 typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
824 thresholdFilter->SetInput(relabelFilter->GetOutput());
825 thresholdFilter->SetUpper(1);
826 if (m_Verbose) std::cout<<"Keeping the first label..."<<std::endl;
827 thresholdFilter->Update();
828 airAdd=thresholdFilter->GetOutput();
831 //---------------------------------
833 //---------------------------------
834 typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
835 inversionFilter->SetInput(airAdd);
836 inversionFilter->SetLowerThreshold(0);
837 inversionFilter->SetUpperThreshold(0);
838 inversionFilter->SetInsideValue(1);
839 inversionFilter->Update();
841 //---------------------------------
843 //---------------------------------
844 typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
845 mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
846 typename InternalImageType::SizeType padSize;
848 padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
849 mirrorPadImageFilter->SetPadLowerBound(padSize);
850 if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
851 mirrorPadImageFilter->Update();
852 airAdd=mirrorPadImageFilter->GetOutput();
854 //---------------------------------
856 //---------------------------------
857 typename InternalImageType::RegionType searchRegion=airAdd->GetLargestPossibleRegion();
858 typename InternalImageType::SizeType searchRegionSize = searchRegion.GetSize();
859 typename InternalImageType::IndexType searchRegionIndex = searchRegion.GetIndex();
860 searchRegionIndex[2]+=searchRegionSize[2]/2;
861 searchRegionSize[2]=1;
862 searchRegion.SetSize(searchRegionSize);
863 searchRegion.SetIndex(searchRegionIndex);
864 IteratorType itAbdomen(airAdd, searchRegion);
865 itAbdomen.GoToBegin();
867 typename InternalImageType::PointType additionalPoint;
868 typename InternalImageType::IndexType aIndex;
870 if (m_Verbose) std::cout<<"Detecting the abdomen..."<<std::endl;
871 while (!itAbdomen.IsAtEnd())
873 if(itAbdomen.Get()) break;
876 aIndex=itAbdomen.GetIndex();
877 airAdd->TransformIndexToPhysicalPoint(aIndex,additionalPoint);
878 if (m_Verbose) std::cout<<"Detected the abdomen in the additional image at "<<additionalPoint<<".."<<std::endl;
880 if(additionalPoint[1]< aPoint[1])
882 aPoint=additionalPoint;
883 if (m_Verbose) std::cout<<"Modifying the detected abdomen to "<<aPoint<<".."<<std::endl;
890 // Determine the detection point
894 if(m_ArgsInfo.offsetDetect_given==Dimension)
895 for(unsigned int i=0; i <Dimension; i++)
896 dPoint[i]+=m_ArgsInfo.offsetDetect_arg[i];
901 if (m_Verbose) std::cout<<"Setting the detection point to "<<dPoint<<".."<<std::endl;
904 //---------------------------------
905 // Pad the rib image and ellips image
906 //---------------------------------
907 typename InternalImageType::Pointer padded_ellips;
908 typename InternalImageType::Pointer padded_bones_low;
910 // If detection point not inside the image: pad
911 typename InternalImageType::IndexType dIndex;
912 if (!bones_low->TransformPhysicalPointToIndex(dPoint, dIndex))
914 typename InternalImageType::SizeType padSize;
916 padSize[1]=abs(dIndex[1])+1;
917 if (m_Verbose) std::cout<<"Padding the images with padding size "<<padSize<<" to include the detection point..."<<std::endl;
919 typename MirrorPadImageFilterType::Pointer padBonesFilter=MirrorPadImageFilterType::New();
920 padBonesFilter->SetInput(bones_low);
921 padBonesFilter->SetPadLowerBound(padSize);
922 padBonesFilter->Update();
923 padded_bones_low=padBonesFilter->GetOutput();
925 typename MirrorPadImageFilterType::Pointer padEllipsFilter=MirrorPadImageFilterType::New();
926 padEllipsFilter->SetInput(m_Ellips);
927 padEllipsFilter->SetPadLowerBound(padSize);
928 padEllipsFilter->Update();
929 padded_ellips=padEllipsFilter->GetOutput();
934 padded_bones_low=bones_low;
935 padded_ellips=m_Ellips;
939 //---------------------------------
940 // Calculate distance map
941 //---------------------------------
942 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
943 distanceMapImageFilter->SetInput(padded_ellips);
944 distanceMapImageFilter->SetInsideIsPositive(false);
945 if (m_Verbose) std::cout<<"Calculating the distance map..."<<std::endl;
946 distanceMapImageFilter->Update();
948 //---------------------------------
949 // Grow while monitoring detection point
950 //---------------------------------
951 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
952 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
953 levelSetFilter->SetFeatureImage( padded_bones_low );
954 levelSetFilter->SetPropagationScaling( 1 );
955 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve1_arg );
956 levelSetFilter->SetAdvectionScaling( 0 );
957 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS1_arg );
958 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter1_arg );
959 levelSetFilter->SetUseImageSpacing(true);
961 // //---------------------------------
962 // // Script for making movie
963 // //---------------------------------
964 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
967 if (m_Verbose) std::cout<<"Starting level set segmentation..."<<std::endl;
968 unsigned int totalNumberOfIterations=0;
971 levelSetFilter->Update();
972 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
974 if ( levelSetFilter->GetOutput()->GetPixel(dIndex) < 0 )
976 if (m_Verbose) std::cout<<"Detection point reached!"<<std::endl;
981 if (m_Verbose) std::cout<<"Detection point not reached after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
982 levelSetFilter->SetInput(levelSetFilter->GetOutput());
983 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose);
986 // std::ostringstream number_str;
987 // number_str << number;
988 // std::string param = number_str.str();
989 // system((script+param).c_str());
993 if ( (totalNumberOfIterations> 10000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS1_arg) ) break;
997 if (m_Verbose) std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
998 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
999 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1002 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1003 thresholder->SetUpperThreshold( 0.0 );
1004 thresholder->SetOutsideValue( 0 );
1005 thresholder->SetInsideValue( 1 );
1006 thresholder->SetInput( levelSetFilter->GetOutput() );
1007 if (m_Verbose) std::cout<<"Thresholding the output level set..."<<std::endl;
1008 thresholder->Update();
1009 grownEllips=thresholder->GetOutput();
1012 //---------------------------------
1013 // Write the grown ellips
1014 //---------------------------------
1015 if (m_ArgsInfo.writeGrownEllips_given)
1017 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1018 caster->SetInput(grownEllips);
1019 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeGrownEllips_arg, m_Verbose);
1023 //----------------------------------------------------------------------------------------------------
1025 //----------------------------------------------------------------------------------------------------
1026 typename InternalImageType::Pointer filledRibs;
1027 if (m_ArgsInfo.filledRibs_given)
1029 typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
1030 featureReader->SetFileName(m_ArgsInfo.filledRibs_arg);
1031 if (m_Verbose) std::cout<<"Reading filled ribs mask..."<<std::endl;
1032 featureReader->Update();
1033 filledRibs=featureReader->GetOutput();
1039 std::cout<<std::endl;
1040 std::cout<<"=========================================="<<std::endl;
1041 std::cout<<"|| Filling the ribs image ||"<<std::endl;
1042 std::cout<<"=========================================="<<std::endl;
1044 //---------------------------------
1045 // Make feature image air+bones
1046 //---------------------------------
1047 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
1048 setBackgroundFilter->SetInput(air_low);
1049 setBackgroundFilter->SetInput2(bones_low);
1050 setBackgroundFilter->SetMaskValue(0);
1051 setBackgroundFilter->SetOutsideValue(0);
1052 setBackgroundFilter->Update();
1053 typename InternalImageType::Pointer bones_air_low =setBackgroundFilter->GetOutput();
1055 //---------------------------------
1056 // Resampling previous solution
1057 //---------------------------------
1058 if (m_Verbose) std::cout<<"Resampling previous solution..."<<std::endl;
1059 typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
1060 typename InternalImageType::PointType origin;
1061 bones_low->TransformIndexToPhysicalPoint(bones_low->GetLargestPossibleRegion().GetIndex(), origin);
1062 resampler->SetOutputOrigin(origin);
1063 resampler->SetOutputSpacing(bones_low->GetSpacing());
1064 typename InternalImageType::SizeType size_low= bones_low->GetLargestPossibleRegion().GetSize();
1065 resampler->SetSize(size_low);
1066 typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
1067 resampler->SetInterpolator(interpolator);
1068 resampler->SetInput(grownEllips);
1069 resampler->Update();
1070 typename InternalImageType::Pointer grownEllips =resampler->GetOutput();
1073 //---------------------------------
1074 // Calculate Distance Map
1075 //---------------------------------
1076 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
1077 distanceMapImageFilter->SetInput(grownEllips);
1078 distanceMapImageFilter->SetInsideIsPositive(false);
1079 if (m_Verbose) std::cout<<"Calculating distance map..."<<std::endl;
1080 distanceMapImageFilter->Update();
1082 //---------------------------------
1083 // Grow while monitoring lung volume coverage
1084 //---------------------------------
1085 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
1086 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
1087 levelSetFilter->SetFeatureImage( bones_air_low );
1088 levelSetFilter->SetPropagationScaling( 1 );
1089 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve2_arg );
1090 levelSetFilter->SetAdvectionScaling( 0 );
1091 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS2_arg );
1092 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter2_arg );
1093 levelSetFilter->SetUseImageSpacing(true);
1095 //---------------------------------
1097 //---------------------------------
1098 typename InversionFilterType::Pointer invertor= InversionFilterType::New();
1099 invertor->SetInput(lungs_low);
1100 invertor->SetLowerThreshold(0);
1101 invertor->SetUpperThreshold(0);
1102 invertor->SetInsideValue(1);
1104 typename InternalImageType::Pointer lungs_low_inv=invertor->GetOutput();
1106 //---------------------------------
1107 // Calculate lung volume
1108 //---------------------------------
1109 typename LabelStatisticsImageFilterType::Pointer statisticsFilter= LabelStatisticsImageFilterType::New();
1110 statisticsFilter->SetInput(lungs_low_inv);
1111 statisticsFilter->SetLabelInput(lungs_low);
1112 statisticsFilter->Update();
1113 unsigned int volume=statisticsFilter->GetSum(0);
1115 // Prepare threshold
1116 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1117 thresholder->SetUpperThreshold( 0.0 );
1118 thresholder->SetOutsideValue( 0 );
1119 thresholder->SetInsideValue( 1 );
1123 unsigned int totalNumberOfIterations=0;
1124 unsigned int coverage=0;
1127 // //---------------------------------
1128 // // Script for making movie
1129 // //---------------------------------
1130 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1134 levelSetFilter->Update();
1135 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1137 thresholder->SetInput( levelSetFilter->GetOutput() );
1138 thresholder->Update();
1139 statisticsFilter->SetInput(thresholder->GetOutput());
1140 statisticsFilter->Update();
1141 coverage=statisticsFilter->GetSum(0);
1143 // Compare the volumes
1144 if ( coverage >= m_ArgsInfo.fillingLevel_arg * volume/100 )
1146 if (m_Verbose) std::cout<<"Lungs filled for "<< m_ArgsInfo.fillingLevel_arg<<"% !"<<std::endl;
1151 if (m_Verbose) std::cout<<"After "<<totalNumberOfIterations<<" iterations, lungs are covered for "
1152 <<(double)coverage/volume*100.0<<"%..."<<std::endl;
1153 levelSetFilter->SetInput(levelSetFilter->GetOutput());
1154 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose);
1157 // std::ostringstream number_str;
1158 // number_str << number;
1159 // std::string param = number_str.str();
1160 // system((script+param).c_str());
1164 if ( (totalNumberOfIterations> 30000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS2_arg) ) break;
1168 if (m_Verbose) std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1169 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1170 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1173 thresholder->SetInput( levelSetFilter->GetOutput() );
1174 thresholder->Update();
1175 filledRibs=thresholder->GetOutput();
1176 // writeImage<InternalImageType>(filledRibs,"/home/jef/tmp/filled_ribs_before.mhd", m_Verbose);
1180 //---------------------------------
1181 // Write the filled ribs
1182 //---------------------------------
1183 if (m_ArgsInfo.writeFilledRibs_given)
1185 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1186 caster->SetInput(filledRibs);
1187 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeFilledRibs_arg, m_Verbose);
1191 //----------------------------------------------------------------------------------------------------
1192 // Collapse to the lungs
1193 //----------------------------------------------------------------------------------------------------
1196 std::cout<<std::endl;
1197 std::cout<<"=========================================="<<std::endl;
1198 std::cout<<"|| Collapsing to the lung image ||"<<std::endl;
1199 std::cout<<"=========================================="<<std::endl;
1202 //---------------------------------
1203 // Make feature image air+bones
1204 //---------------------------------
1205 if (m_Verbose) std::cout<<"Making feature images..."<<std::endl;
1206 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
1207 setBackgroundFilter->SetInput(air);
1208 setBackgroundFilter->SetInput2(bones);
1209 setBackgroundFilter->SetMaskValue(0);
1210 setBackgroundFilter->SetOutsideValue(0);
1211 setBackgroundFilter->Update();
1212 typename InternalImageType::Pointer bones_air =setBackgroundFilter->GetOutput();
1214 typename SetBackgroundFilterType::Pointer setBackgroundFilter2 = SetBackgroundFilterType::New();
1215 setBackgroundFilter2->SetInput(bones_air);
1216 setBackgroundFilter2->SetInput2(lungs);
1217 setBackgroundFilter2->SetMaskValue(0);
1218 setBackgroundFilter2->SetOutsideValue(0);
1219 setBackgroundFilter2->Update();
1220 typename InternalImageType::Pointer lungs_bones_air =setBackgroundFilter2->GetOutput();
1222 //---------------------------------
1223 // Prepare previous solution
1224 //---------------------------------
1225 if (m_Verbose) std::cout<<"Resampling previous solution..."<<std::endl;
1226 typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
1227 typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
1228 typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
1229 resampler->SetOutputStartIndex(bones->GetLargestPossibleRegion().GetIndex());
1230 resampler->SetInput(filledRibs);
1231 resampler->SetInterpolator(interpolator);
1232 resampler->SetOutputParametersFromImage(bones);
1233 resampler->Update();
1234 filledRibs =resampler->GetOutput();
1236 if (m_Verbose) std::cout<<"Setting lungs to 1..."<<std::endl;
1237 typename SetBackgroundFilterType::Pointer setBackgroundFilter3 = SetBackgroundFilterType::New();
1238 setBackgroundFilter3->SetInput(filledRibs);
1239 setBackgroundFilter3->SetInput2(lungs);
1240 setBackgroundFilter3->SetMaskValue(0);
1241 setBackgroundFilter3->SetOutsideValue(1);
1242 setBackgroundFilter3->Update();
1243 filledRibs=setBackgroundFilter3->GetOutput();
1245 if (m_Verbose) std::cout<<"Removing overlap with bones..."<<std::endl;
1246 typename SetBackgroundFilterType::Pointer setBackgroundFilter4 = SetBackgroundFilterType::New();
1247 setBackgroundFilter4->SetInput(filledRibs);
1248 setBackgroundFilter4->SetInput2(bones);
1249 setBackgroundFilter4->SetMaskValue(0);
1250 setBackgroundFilter4->SetOutsideValue(0);
1251 setBackgroundFilter4->Update();
1252 filledRibs =setBackgroundFilter4->GetOutput();
1253 // writeImage<InternalImageType>(filledRibs,"/home/jef/tmp/filledRibs_pp.mhd");
1254 //---------------------------------
1255 // Calculate Distance Map
1256 //---------------------------------
1257 // typedef itk::ApproximateSignedDistanceMapImageFilter <InternalImageType, LevelSetImageType> DistanceMapImageFilterType2;
1258 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
1259 distanceMapImageFilter->SetInput(filledRibs);
1260 distanceMapImageFilter->SetInsideIsPositive(false);
1261 // distanceMapImageFilter->SetInsideValue(0);
1262 // distanceMapImageFilter->SetOutsideValue(1);
1263 if (m_Verbose) std::cout<<"Calculating distance map..."<<std::endl;
1264 distanceMapImageFilter->Update();
1266 //---------------------------------
1268 //---------------------------------
1269 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
1270 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
1271 levelSetFilter->SetFeatureImage( lungs_bones_air );
1272 levelSetFilter->SetPropagationScaling(m_ArgsInfo.prop3_arg);
1273 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve3_arg );
1274 levelSetFilter->SetAdvectionScaling( 0 );
1275 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS3_arg );
1276 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter3_arg );
1277 levelSetFilter->SetUseImageSpacing(true);
1280 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1283 if (m_Verbose) std::cout<<"Starting the levelset..."<<std::endl;
1284 unsigned int totalNumberOfIterations=0;
1287 levelSetFilter->Update();
1290 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1291 levelSetFilter->SetInput(levelSetFilter->GetOutput());
1292 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg);
1293 std::cout << "After "<<totalNumberOfIterations<<" iteration the RMS change is " << levelSetFilter->GetRMSChange() <<"..."<< std::endl;
1296 // std::ostringstream number_str;
1297 // number_str << number;
1298 // std::string param = number_str.str();
1299 // system((script+param).c_str());
1302 // stopping condition
1303 if ( (totalNumberOfIterations>= (unsigned int) m_ArgsInfo.maxIter3_arg) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS3_arg) ) break;
1304 levelSetFilter->SetNumberOfIterations( std::min ((unsigned int) m_ArgsInfo.iter3_arg, (unsigned int) m_ArgsInfo.maxIter3_arg-totalNumberOfIterations ) );
1310 std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1311 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1312 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1316 typedef itk::BinaryThresholdImageFilter< LevelSetImageType,InternalImageType> LevelSetBinarizeFilterType;
1317 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1318 thresholder->SetUpperThreshold( 0.0 );
1319 thresholder->SetOutsideValue( 0 );
1320 thresholder->SetInsideValue( 1 );
1321 thresholder->SetInput( levelSetFilter->GetOutput() );
1322 thresholder->Update();
1323 typename InternalImageType::Pointer output = thresholder->GetOutput();
1326 //----------------------------------------------------------------------------------------------------
1327 // Prepare the output
1328 //----------------------------------------------------------------------------------------------------
1330 //---------------------------------
1332 //---------------------------------
1333 if (m_Verbose) std::cout<<"Removing overlap mask with air..."<<std::endl;
1334 typename SetBackgroundFilterType::Pointer setBackgroundFilter5 = SetBackgroundFilterType::New();
1335 setBackgroundFilter5->SetInput(output);
1336 setBackgroundFilter5->SetInput2(air);
1337 setBackgroundFilter5->SetMaskValue(0);
1338 setBackgroundFilter5->SetOutsideValue(0);
1339 setBackgroundFilter5->Update();
1340 output=setBackgroundFilter5->GetOutput();
1342 //---------------------------------
1343 // Open & close to cleanup
1344 //---------------------------------
1345 if ( m_ArgsInfo.openClose_flag)
1348 //---------------------------------
1349 // Structuring element
1350 //---------------------------------
1351 typedef itk::BinaryBallStructuringElement<InternalPixelType,Dimension > KernelType;
1352 KernelType structuringElement;
1353 structuringElement.SetRadius(1);
1354 structuringElement.CreateStructuringElement();
1356 //---------------------------------
1358 //---------------------------------
1359 typedef itk::BinaryMorphologicalOpeningImageFilter<InternalImageType, InternalImageType , KernelType> OpenFilterType;
1360 typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
1361 openFilter->SetInput(output);
1362 openFilter->SetBackgroundValue(0);
1363 openFilter->SetForegroundValue(1);
1364 openFilter->SetKernel(structuringElement);
1365 if(m_Verbose) std::cout<<"Opening the output image..."<<std::endl;
1367 //---------------------------------
1369 //---------------------------------
1370 typedef itk::BinaryMorphologicalClosingImageFilter<InternalImageType, InternalImageType , KernelType> CloseFilterType;
1371 typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
1372 closeFilter->SetInput(openFilter->GetOutput());
1373 closeFilter->SetSafeBorder(true);
1374 closeFilter->SetForegroundValue(1);
1375 closeFilter->SetKernel(structuringElement);
1376 if(m_Verbose) std::cout<<"Closing the output image..."<<std::endl;
1377 closeFilter->Update();
1378 output=closeFilter->GetOutput();
1381 writeImage<InternalImageType>(output,"/home/jef/tmp/mm_double.mhd");
1383 // Extract the upper part
1384 typedef itk::CropImageFilter<InternalImageType, InternalImageType> CropImageFilterType;
1385 typename CropImageFilterType::Pointer cropFilter=CropImageFilterType::New();
1386 cropFilter->SetInput(output);
1387 typename InputImageType::SizeType lSize, uSize;
1390 lSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
1391 cropFilter->SetLowerBoundaryCropSize(lSize);
1392 cropFilter->SetUpperBoundaryCropSize(uSize);
1393 cropFilter->Update();
1396 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1397 caster->SetInput(cropFilter->GetOutput());
1398 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.output_arg, m_Verbose);
1404 #endif //#define clitkMotionMaskGenericFilter_txx