1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://www.centreleonberard.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
9 This software is distributed WITHOUT ANY WARRANTY; without even
10 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 PURPOSE. See the copyright notices for more information.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ===========================================================================**/
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 //---------------------------------
103 //---------------------------------
104 if(m_ArgsInfo.pad_flag) {
105 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
106 IteratorType it(air, air->GetLargestPossibleRegion());
107 typename InternalImageType::IndexType index;
108 while(!it.IsAtEnd()) {
111 //Pad borders of all dimensions but 2 (cranio-caudal)
112 for (unsigned int i=0; i<Dimension; i++){
115 if(index[i]==air->GetLargestPossibleRegion().GetIndex()[i]
116 || index[i]==(unsigned int)air->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air->GetLargestPossibleRegion().GetSize()[i]-1)
123 if (m_Verbose) std::cout<<"Extracting the air feature image..."<<std::endl;
124 //---------------------------------
125 // Binarize the image
126 //---------------------------------
127 typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
128 binarizeFilter->SetInput(input);
129 binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdAir_arg));
130 binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdAir_arg));
131 if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<<m_ArgsInfo.lowerThresholdAir_arg
132 <<", "<<m_ArgsInfo.upperThresholdAir_arg<<"..."<<std::endl;
133 binarizeFilter->Update();
134 air = binarizeFilter->GetOutput();
136 //---------------------------------
138 //---------------------------------
139 if(m_ArgsInfo.pad_flag) {
140 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
141 IteratorType it(air, air->GetLargestPossibleRegion());
142 typename InternalImageType::IndexType index;
143 while(!it.IsAtEnd()) {
146 //Pad borders of all dimensions but 2 (cranio-caudal)
147 for (unsigned int i=0; i<Dimension; i++){
150 if(index[i]==air->GetLargestPossibleRegion().GetIndex()[i]
151 || index[i]==(unsigned int)air->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air->GetLargestPossibleRegion().GetSize()[i]-1)
152 it.Set(binarizeFilter->GetInsideValue());
158 //---------------------------------
159 // Remove lungs (in place)
160 //---------------------------------
161 typedef itk::ImageRegionIterator<InternalImageType> IteratorType;
162 IteratorType itAir(air, binarizeFilter->GetOutput()->GetLargestPossibleRegion());
163 IteratorType itLungs(lungs, binarizeFilter->GetOutput()->GetLargestPossibleRegion()); //lungs padded, used air region
166 while(!itAir.IsAtEnd()) {
167 if(!itLungs.Get()) // The lungs are set to 0 at that stage
173 //---------------------------------
174 // Label the connected components
175 //---------------------------------
176 typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
177 connectFilter->SetInput(binarizeFilter->GetOutput());
178 connectFilter->SetBackgroundValue(0);
179 connectFilter->SetFullyConnected(false);
180 if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
182 //---------------------------------
183 // Sort the labels according to size
184 //---------------------------------
185 typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
186 relabelFilter->SetInput(connectFilter->GetOutput());
187 if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
189 //---------------------------------
191 //---------------------------------
192 typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
193 thresholdFilter->SetInput(relabelFilter->GetOutput());
194 thresholdFilter->SetUpper(1);
195 if (m_Verbose) std::cout<<"Keeping the first label..."<<std::endl;
196 thresholdFilter->Update();
197 air=thresholdFilter->GetOutput();
200 //---------------------------------
202 //---------------------------------
203 typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
204 inversionFilter->SetInput(air);
205 inversionFilter->SetLowerThreshold(0);
206 inversionFilter->SetUpperThreshold(0);
207 inversionFilter->SetInsideValue(1);
208 inversionFilter->Update();
210 //---------------------------------
212 //---------------------------------
213 typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
214 mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
215 typename InternalImageType::SizeType padSize;
217 padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
218 mirrorPadImageFilter->SetPadLowerBound(padSize);
219 if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
220 mirrorPadImageFilter->Update();
221 air=mirrorPadImageFilter->GetOutput();
222 //writeImage<InternalImageType>(air,"/home/srit/tmp/air.mhd");
228 //-------------------------------------------------------------------
230 //-------------------------------------------------------------------
231 template <unsigned int Dimension, class PixelType>
232 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
233 MotionMaskGenericFilter::GetBonesImage(typename itk::Image<PixelType, Dimension>::Pointer input )
237 typedef int InternalPixelType;
238 typedef itk::Image<PixelType, Dimension> InputImageType;
239 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
241 //----------------------------------------------------------------------------------------------------
242 // Typedef for Processing
243 //----------------------------------------------------------------------------------------------------
244 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
245 typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
246 typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> InversionFilterType;
247 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdFilterType;
248 typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
249 typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
250 typedef itk::MirrorPadImageFilter<InternalImageType, InternalImageType> MirrorPadImageFilterType;
253 typename InternalImageType::Pointer bones;
254 if (m_ArgsInfo.featureBones_given) {
255 typename FeatureReaderType::Pointer featureReader =FeatureReaderType::New();
256 featureReader->SetFileName(m_ArgsInfo.featureBones_arg);
257 if (m_Verbose) std::cout<<"Reading the bones feature image..."<<std::endl;
258 featureReader->Update();
259 bones=featureReader->GetOutput();
261 if (m_Verbose) std::cout<<"Extracting the bones feature image..."<<std::endl;
262 //---------------------------------
263 // Binarize the image
264 //---------------------------------
265 typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
266 binarizeFilter->SetInput(input);
267 binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdBones_arg));
268 binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdBones_arg));
269 if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<<m_ArgsInfo.lowerThresholdBones_arg
270 <<", "<<m_ArgsInfo.upperThresholdBones_arg<<"..."<<std::endl;
272 //---------------------------------
273 // Label the connected components
274 //---------------------------------
275 typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
276 connectFilter->SetInput(binarizeFilter->GetOutput());
277 connectFilter->SetBackgroundValue(0);
278 connectFilter->SetFullyConnected(false);
279 if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
281 //---------------------------------
282 // Sort the labels according to size
283 //---------------------------------
284 typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
285 relabelFilter->SetInput(connectFilter->GetOutput());
286 if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
288 //---------------------------------
290 //---------------------------------
291 typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
292 thresholdFilter->SetInput(relabelFilter->GetOutput());
293 thresholdFilter->SetUpper(1);
294 if (m_Verbose) std::cout<<"Keeping the first label..."<<std::endl;
295 thresholdFilter->Update();
296 bones=thresholdFilter->GetOutput();
300 //---------------------------------
302 //---------------------------------
303 typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
304 inversionFilter->SetInput(bones);
305 inversionFilter->SetLowerThreshold(0);
306 inversionFilter->SetUpperThreshold(0);
307 inversionFilter->SetInsideValue(1);
308 inversionFilter->Update();
310 //---------------------------------
312 //---------------------------------
313 typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
314 mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
315 typename InternalImageType::SizeType padSize;
317 padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
318 mirrorPadImageFilter->SetPadLowerBound(padSize);
319 if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
320 mirrorPadImageFilter->Update();
321 bones=mirrorPadImageFilter->GetOutput();
322 // writeImage<InternalImageType>(bones,"/home/jef/tmp/bones.mhd");
330 //-------------------------------------------------------------------
332 //-------------------------------------------------------------------
333 template <unsigned int Dimension, class PixelType>
334 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
335 MotionMaskGenericFilter::GetLungsImage(typename itk::Image<PixelType, Dimension>::Pointer input )
338 typedef int InternalPixelType;
339 typedef itk::Image<PixelType, Dimension> InputImageType;
340 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
342 //----------------------------------------------------------------------------------------------------
343 // Typedef for Processing
344 //----------------------------------------------------------------------------------------------------
345 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
346 typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
347 typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> InversionFilterType;
348 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdFilterType;
349 typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
350 typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
351 typedef itk::MirrorPadImageFilter<InternalImageType, InternalImageType> MirrorPadImageFilterType;
353 typename InternalImageType::Pointer lungs;
354 if (m_ArgsInfo.featureLungs_given) {
355 typename FeatureReaderType::Pointer featureReader =FeatureReaderType::New();
356 featureReader->SetFileName(m_ArgsInfo.featureLungs_arg);
357 if (m_Verbose) std::cout<<"Reading the lungs feature image..."<<std::endl;
358 featureReader->Update();
359 lungs=featureReader->GetOutput();
361 if (m_Verbose) std::cout<<"Extracting the lungs feature image..."<<std::endl;
362 //---------------------------------
363 // Binarize the image
364 //---------------------------------
365 typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
366 binarizeFilter->SetInput(input);
367 binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdLungs_arg));
368 binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdLungs_arg));
369 if (m_Verbose) std::cout<<"Binarizing the image using a threshold "<<m_ArgsInfo.lowerThresholdLungs_arg
370 <<", "<<m_ArgsInfo.upperThresholdLungs_arg<<"..."<<std::endl;
372 //---------------------------------
373 // Label the connected components
374 //---------------------------------
375 typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
376 connectFilter->SetInput(binarizeFilter->GetOutput());
377 connectFilter->SetBackgroundValue(0);
378 connectFilter->SetFullyConnected(true);
379 if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
381 //---------------------------------
382 // Sort the labels according to size
383 //---------------------------------
384 typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
385 relabelFilter->SetInput(connectFilter->GetOutput());
386 if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
387 // writeImage<InternalImageType> (relabelFilter->GetOutput(), "/home/vdelmon/tmp/labels.mhd");
389 //---------------------------------
391 //---------------------------------
392 typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
393 thresholdFilter->SetInput(relabelFilter->GetOutput());
394 thresholdFilter->SetLower(1);
395 thresholdFilter->SetUpper(2);
396 if (m_Verbose) std::cout<<"Keeping the first two labels..."<<std::endl;
397 thresholdFilter->Update();
398 lungs=thresholdFilter->GetOutput();
403 //---------------------------------
405 //---------------------------------
406 typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
407 inversionFilter->SetInput(lungs);
408 inversionFilter->SetLowerThreshold(0);
409 inversionFilter->SetUpperThreshold(0);
410 inversionFilter->SetInsideValue(1);
411 inversionFilter->Update();
413 //---------------------------------
415 //---------------------------------
416 typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
417 mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
418 typename InternalImageType::SizeType padSize;
420 padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
421 mirrorPadImageFilter->SetPadLowerBound(padSize);
422 if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
423 mirrorPadImageFilter->Update();
424 lungs=mirrorPadImageFilter->GetOutput();
425 // writeImage<InternalImageType>(lungs,"/home/jef/tmp/lungs.mhd");
430 //-------------------------------------------------------------------
432 //-------------------------------------------------------------------
433 template <unsigned int Dimension, class PixelType>
434 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
435 MotionMaskGenericFilter::Resample( typename itk::Image<InternalPixelType,Dimension>::Pointer input )
438 typedef int InternalPixelType;
439 typedef itk::Image<PixelType, Dimension> InputImageType;
440 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
442 //---------------------------------
443 // Resample to new spacing
444 //---------------------------------
445 typename InternalImageType::SizeType size_low;
446 typename InternalImageType::SpacingType spacing_low;
447 for (unsigned int i=0; i<Dimension; i++)
448 if (m_ArgsInfo.spacing_given==Dimension)
449 for (unsigned int i=0; i<Dimension; i++) {
450 spacing_low[i]=m_ArgsInfo.spacing_arg[i];
451 size_low[i]=ceil(input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i]/spacing_low[i]);
454 for (unsigned int i=0; i<Dimension; i++) {
455 spacing_low[i]=m_ArgsInfo.spacing_arg[0];
456 size_low[i]=ceil(input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i]/spacing_low[i]);
459 typename InternalImageType::PointType origin;
460 input->TransformIndexToPhysicalPoint(input->GetLargestPossibleRegion().GetIndex(), origin);
461 typedef itk::ResampleImageFilter<InternalImageType, InternalImageType> ResampleImageFilterType;
462 typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
463 typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
464 typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
465 resampler->SetInterpolator(interpolator);
466 resampler->SetOutputOrigin(origin);
467 resampler->SetOutputSpacing(spacing_low);
468 resampler->SetSize(size_low);
469 resampler->SetInput(input);
471 typename InternalImageType::Pointer output =resampler->GetOutput();
476 //-------------------------------------------------------------------
478 //-------------------------------------------------------------------
479 template <unsigned int Dimension, class PixelType>
480 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
481 MotionMaskGenericFilter::InitializeEllips( typename itk::Vector<double,Dimension> center, typename itk::Image<InternalPixelType,Dimension>::Pointer bones_low )
485 typedef int InternalPixelType;
486 typedef itk::Image<PixelType, Dimension> InputImageType;
487 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
488 typedef itk::Image< unsigned char , Dimension> OutputImageType;
489 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
490 typedef clitk::SetBackgroundImageFilter<InternalImageType,InternalImageType, InternalImageType> SetBackgroundFilterType;
491 typedef itk::LabelStatisticsImageFilter<InternalImageType,InternalImageType> LabelStatisticsImageFilterType;
492 typedef itk::CastImageFilter<InternalImageType,OutputImageType> CastImageFilterType;
495 typename InternalImageType::Pointer ellips;
496 if (m_ArgsInfo.ellips_given || m_ArgsInfo.grownEllips_given || m_ArgsInfo.filledRibs_given) {
497 if(m_ArgsInfo.ellips_given) {
498 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
499 typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
500 featureReader->SetFileName(m_ArgsInfo.ellips_arg);
501 featureReader->Update();
502 ellips=featureReader->GetOutput();
506 std::cout<<std::endl;
507 std::cout<<"=========================================="<<std::endl;
508 std::cout<<"|| Initializing ellipsoide ||"<<std::endl;
509 std::cout<<"=========================================="<<std::endl;
512 //---------------------------------
513 // Offset from center
514 //---------------------------------
515 typename itk::Vector<double,Dimension> offset;
516 if(m_ArgsInfo.offset_given== Dimension) {
517 for(unsigned int i=0; i<Dimension; i++)
518 offset[i]=m_ArgsInfo.offset_arg[i];
523 itk::Vector<double,Dimension> centerEllips=center+offset;
525 std::cout<<"Placing the center of the initial ellipsoide at (mm) "<<std::endl;
526 std::cout<<centerEllips[0];
527 for (unsigned int i=1; i<Dimension; i++)
528 std::cout<<" "<<centerEllips[i];
529 std::cout<<std::endl;
532 //---------------------------------
534 //---------------------------------
535 ellips=InternalImageType::New();
536 ellips->SetRegions (bones_low->GetLargestPossibleRegion());
537 ellips->SetOrigin(bones_low->GetOrigin());
538 ellips->SetSpacing(bones_low->GetSpacing());
540 ellips->FillBuffer(0);
543 typename itk::Vector<double, Dimension> axes;
544 if (m_ArgsInfo.axes_given==Dimension)
545 for(unsigned int i=0; i<Dimension; i++)
546 axes[i]=m_ArgsInfo.axes_arg[i];
554 IteratorType itEllips(ellips, ellips->GetLargestPossibleRegion());
555 itEllips.GoToBegin();
556 typename InternalImageType::PointType point;
557 typename InternalImageType::IndexType index;
560 if (m_Verbose) std::cout<<"Drawing the initial ellipsoide..."<<std::endl;
561 while (!itEllips.IsAtEnd()) {
562 index=itEllips.GetIndex();
563 ellips->TransformIndexToPhysicalPoint(index, point);
565 for(unsigned int i=0; i<Dimension; i++)
566 distance+=powf( ( (centerEllips[i]-point[i])/axes[i] ), 2);
575 //---------------------------------
577 //---------------------------------
578 if (m_ArgsInfo.writeEllips_given) {
579 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
580 caster->SetInput(ellips);
581 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeEllips_arg, m_Verbose);
589 //-------------------------------------------------------------------
590 // Update with the number of dimensions and the pixeltype
591 //-------------------------------------------------------------------
592 template <unsigned int Dimension, class PixelType>
594 MotionMaskGenericFilter::UpdateWithDimAndPixelType()
598 typedef int InternalPixelType;
599 typedef itk::Image<PixelType, Dimension> InputImageType;
600 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
601 typedef itk::Image< float, Dimension> LevelSetImageType; //labeling doesn't work with unsigned char?
602 typedef itk::Image<unsigned char, Dimension> OutputImageType;
605 //----------------------------------------------------------------------------------------------------
606 // Typedef for Processing
607 //----------------------------------------------------------------------------------------------------
608 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
609 typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
610 typedef itk::BinaryThresholdImageFilter< LevelSetImageType,InternalImageType> LevelSetBinarizeFilterType;
611 typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> InversionFilterType;
612 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdFilterType;
613 typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
614 typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
615 typedef itk::MirrorPadImageFilter<InternalImageType, InternalImageType> MirrorPadImageFilterType;
616 typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
617 typedef itk::ResampleImageFilter<InternalImageType, InternalImageType> ResampleImageFilterType;
618 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
619 typedef itk::GeodesicActiveContourLevelSetImageFilter<LevelSetImageType, InternalImageType> LevelSetImageFilterType;
620 typedef itk::SignedDanielssonDistanceMapImageFilter<InternalImageType, LevelSetImageType> DistanceMapImageFilterType;
621 typedef clitk::SetBackgroundImageFilter<InternalImageType,InternalImageType, InternalImageType> SetBackgroundFilterType;
622 typedef itk::LabelStatisticsImageFilter<InternalImageType,InternalImageType> LabelStatisticsImageFilterType;
623 typedef itk::CastImageFilter<InternalImageType,OutputImageType> CastImageFilterType;
627 typedef itk::ImageFileReader<InputImageType> InputReaderType;
628 typename InputReaderType::Pointer reader = InputReaderType::New();
629 reader->SetFileName( m_InputFileName);
631 typename InputImageType::Pointer input= reader->GetOutput();
633 // // globals for avi
634 // unsigned int number=0;
638 std::cout<<std::endl;
639 std::cout<<"=========================================="<<std::endl;
640 std::cout<<"|| Making feature images ||"<<std::endl;
641 std::cout<<"=========================================="<<std::endl;
644 //--------------------------------------------------------------------------------
646 //-------------------------------------------------------------------------------
647 typename InternalImageType::Pointer lungs = this->GetLungsImage<Dimension, PixelType>(input);
649 //-------------------------------------------------------------------------------
651 //-------------------------------------------------------------------------------
652 typename InternalImageType::Pointer air = this->GetAirImage<Dimension, PixelType>(input, lungs);
654 //-------------------------------------------------------------------------------
656 //-------------------------------------------------------------------------------
657 typename InternalImageType::Pointer bones = this->GetBonesImage<Dimension, PixelType>(input);
659 //----------------------------------------------------------------------------------------------------
660 // Write feature images
661 //----------------------------------------------------------------------------------------------------
662 if(m_ArgsInfo.writeFeature_given) {
663 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
664 setBackgroundFilter->SetInput(air);
665 setBackgroundFilter->SetInput2(bones);
666 setBackgroundFilter->SetMaskValue(0);
667 setBackgroundFilter->SetOutsideValue(2);
668 setBackgroundFilter->Update();
669 typename InternalImageType::Pointer bones_air =setBackgroundFilter->GetOutput();
671 typename SetBackgroundFilterType::Pointer setBackgroundFilter2 = SetBackgroundFilterType::New();
672 setBackgroundFilter2->SetInput(bones_air);
673 setBackgroundFilter2->SetInput2(lungs);
674 setBackgroundFilter2->SetMaskValue(0);
675 setBackgroundFilter2->SetOutsideValue(3);
676 setBackgroundFilter2->Update();
677 typename InternalImageType::Pointer lungs_bones_air =setBackgroundFilter2->GetOutput();
679 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
680 caster->SetInput(lungs_bones_air);
681 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeFeature_arg, m_Verbose);
684 //----------------------------------------------------------------------------------------------------
685 // Low dimensional versions
686 //----------------------------------------------------------------------------------------------------
687 typename InternalImageType::Pointer bones_low =Resample<Dimension , PixelType>(bones);
688 typename InternalImageType::Pointer lungs_low =Resample<Dimension , PixelType>(lungs);
689 typename InternalImageType::Pointer air_low =Resample<Dimension , PixelType>(air);
691 //---------------------------------
693 //---------------------------------
694 if(m_ArgsInfo.pad_flag) {
695 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
696 IteratorType it(air_low, air_low->GetLargestPossibleRegion());
697 typename InternalImageType::IndexType index;
698 while(!it.IsAtEnd()) {
700 for (unsigned int i=0; i<Dimension; i++)
701 if(index[i]==air_low->GetLargestPossibleRegion().GetIndex()[i]
702 || index[i]==(unsigned int)air_low->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air_low->GetLargestPossibleRegion().GetSize()[i]-1)
708 //---------------------------------
710 //---------------------------------
711 typename itk::Vector<double,Dimension> center;
712 typedef itk::ImageMomentsCalculator<InternalImageType> MomentsCalculatorType;
713 typename MomentsCalculatorType::Pointer momentsCalculator=MomentsCalculatorType::New();
714 momentsCalculator->SetImage(air);
715 momentsCalculator->Compute();
716 center=momentsCalculator->GetCenterOfGravity();
718 std::cout<<"The center of gravity of the patient body is located at (mm) "<<std::endl;
719 std::cout<<center[0];
720 for (unsigned int i=1; i<Dimension; i++)
721 std::cout<<" "<<center[i];
722 std::cout<<std::endl;
725 //----------------------------------------------------------------------------------------------------
726 // Ellipsoide initialization
727 //----------------------------------------------------------------------------------------------------
728 typename InternalImageType::Pointer m_Ellips= InitializeEllips<Dimension, PixelType>(center,bones_low);
730 //----------------------------------------------------------------------------------------------------
732 //----------------------------------------------------------------------------------------------------
733 typename InternalImageType::Pointer grownEllips;
734 if (m_ArgsInfo.grownEllips_given || m_ArgsInfo.filledRibs_given) {
735 if (m_ArgsInfo.grownEllips_given) {
737 typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
738 featureReader->SetFileName(m_ArgsInfo.grownEllips_arg);
739 featureReader->Update();
740 grownEllips=featureReader->GetOutput();
745 std::cout<<std::endl;
746 std::cout<<"=========================================="<<std::endl;
747 std::cout<<"|| Growing ellipsoide ||"<<std::endl;
748 std::cout<<"=========================================="<<std::endl;
751 //---------------------------------
753 //---------------------------------
754 typename InternalImageType::PointType dPoint;
755 if (m_ArgsInfo.detectionPoint_given) {
756 for (unsigned int i=0; i<Dimension; i++)
757 dPoint[i]=m_ArgsInfo.detectionPoint_arg[i];
759 typename InternalImageType::RegionType searchRegion=air->GetLargestPossibleRegion();
760 typename InternalImageType::SizeType searchRegionSize = searchRegion.GetSize();
761 typename InternalImageType::IndexType searchRegionIndex = searchRegion.GetIndex();
762 searchRegionIndex[2]+=searchRegionSize[2]/2;
763 searchRegionSize[2]=1;
764 searchRegion.SetSize(searchRegionSize);
765 searchRegion.SetIndex(searchRegionIndex);
766 IteratorType itAbdomen(air, searchRegion);
767 itAbdomen.GoToBegin();
769 typename InternalImageType::PointType aPoint;
770 typename InternalImageType::IndexType aIndex;
772 if (m_Verbose) std::cout<<"Detecting the abdomen..."<<std::endl;
773 while (!itAbdomen.IsAtEnd()) {
774 if(itAbdomen.Get()) break;
777 aIndex=itAbdomen.GetIndex();
778 air->TransformIndexToPhysicalPoint(aIndex,aPoint);
779 if (m_Verbose) std::cout<<"Detected the abdomen at "<<aPoint<<".."<<std::endl;
782 //---------------------------------
783 // Detect abdomen in additional images?
784 //---------------------------------
785 if (m_ArgsInfo.detectionPairs_given) {
786 for (unsigned int i=0; i<m_ArgsInfo.detectionPairs_given; i++) {
787 typename InternalImageType::Pointer airAdd;
788 //---------------------------------
790 //--------------------------------
791 typedef itk::ImageFileReader<InputImageType> InputReaderType;
792 typename InputReaderType::Pointer reader = InputReaderType::New();
793 reader->SetFileName( m_ArgsInfo.detectionPairs_arg[i]);
795 typename InputImageType::Pointer additional= reader->GetOutput();
796 if (m_Verbose) std::cout<<"Extracting the air from additional image "<< i<<"..."<<std::endl;
798 //---------------------------------
799 // Binarize the image
800 //---------------------------------
801 typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
802 binarizeFilter->SetInput(additional);
803 binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdAir_arg));
804 binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdAir_arg));
805 if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<<m_ArgsInfo.lowerThresholdAir_arg
806 <<", "<<m_ArgsInfo.upperThresholdAir_arg<<"..."<<std::endl;
808 //---------------------------------
809 // Label the connected components
810 //---------------------------------
811 typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
812 connectFilter->SetInput(binarizeFilter->GetOutput());
813 connectFilter->SetBackgroundValue(0);
814 connectFilter->SetFullyConnected(false);
815 if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
817 //---------------------------------
818 // Sort the labels according to size
819 //---------------------------------
820 typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
821 relabelFilter->SetInput(connectFilter->GetOutput());
822 if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
824 //---------------------------------
826 //---------------------------------
827 typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
828 thresholdFilter->SetInput(relabelFilter->GetOutput());
829 thresholdFilter->SetUpper(1);
830 if (m_Verbose) std::cout<<"Keeping the first label..."<<std::endl;
831 thresholdFilter->Update();
832 airAdd=thresholdFilter->GetOutput();
835 //---------------------------------
837 //---------------------------------
838 typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
839 inversionFilter->SetInput(airAdd);
840 inversionFilter->SetLowerThreshold(0);
841 inversionFilter->SetUpperThreshold(0);
842 inversionFilter->SetInsideValue(1);
843 inversionFilter->Update();
845 //---------------------------------
847 //---------------------------------
848 typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
849 mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
850 typename InternalImageType::SizeType padSize;
852 padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
853 mirrorPadImageFilter->SetPadLowerBound(padSize);
854 if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
855 mirrorPadImageFilter->Update();
856 airAdd=mirrorPadImageFilter->GetOutput();
858 //---------------------------------
860 //---------------------------------
861 typename InternalImageType::RegionType searchRegion=airAdd->GetLargestPossibleRegion();
862 typename InternalImageType::SizeType searchRegionSize = searchRegion.GetSize();
863 typename InternalImageType::IndexType searchRegionIndex = searchRegion.GetIndex();
864 searchRegionIndex[2]+=searchRegionSize[2]/2;
865 searchRegionSize[2]=1;
866 searchRegion.SetSize(searchRegionSize);
867 searchRegion.SetIndex(searchRegionIndex);
868 IteratorType itAbdomen(airAdd, searchRegion);
869 itAbdomen.GoToBegin();
871 typename InternalImageType::PointType additionalPoint;
872 typename InternalImageType::IndexType aIndex;
874 if (m_Verbose) std::cout<<"Detecting the abdomen..."<<std::endl;
875 while (!itAbdomen.IsAtEnd()) {
876 if(itAbdomen.Get()) break;
879 aIndex=itAbdomen.GetIndex();
880 airAdd->TransformIndexToPhysicalPoint(aIndex,additionalPoint);
881 if (m_Verbose) std::cout<<"Detected the abdomen in the additional image at "<<additionalPoint<<".."<<std::endl;
883 if(additionalPoint[1]< aPoint[1]) {
884 aPoint=additionalPoint;
885 if (m_Verbose) std::cout<<"Modifying the detected abdomen to "<<aPoint<<".."<<std::endl;
892 // Determine the detection point
896 if(m_ArgsInfo.offsetDetect_given==Dimension)
897 for(unsigned int i=0; i <Dimension; i++)
898 dPoint[i]+=m_ArgsInfo.offsetDetect_arg[i];
903 if (m_Verbose) std::cout<<"Setting the detection point to "<<dPoint<<".."<<std::endl;
906 //---------------------------------
907 // Pad the rib image and ellips image
908 //---------------------------------
909 typename InternalImageType::Pointer padded_ellips;
910 typename InternalImageType::Pointer padded_bones_low;
912 // If detection point not inside the image: pad
913 typename InternalImageType::IndexType dIndex;
914 if (!bones_low->TransformPhysicalPointToIndex(dPoint, dIndex)) {
915 typename InternalImageType::SizeType padSize;
917 padSize[1]=abs(dIndex[1])+1;
918 if (m_Verbose) std::cout<<"Padding the images with padding size "<<padSize<<" to include the detection point..."<<std::endl;
920 typename MirrorPadImageFilterType::Pointer padBonesFilter=MirrorPadImageFilterType::New();
921 padBonesFilter->SetInput(bones_low);
922 padBonesFilter->SetPadLowerBound(padSize);
923 padBonesFilter->Update();
924 padded_bones_low=padBonesFilter->GetOutput();
926 typename MirrorPadImageFilterType::Pointer padEllipsFilter=MirrorPadImageFilterType::New();
927 padEllipsFilter->SetInput(m_Ellips);
928 padEllipsFilter->SetPadLowerBound(padSize);
929 padEllipsFilter->Update();
930 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;
970 levelSetFilter->Update();
971 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
973 if ( levelSetFilter->GetOutput()->GetPixel(dIndex) < 0 ) {
974 if (m_Verbose) std::cout<<"Detection point reached!"<<std::endl;
977 if (m_Verbose) std::cout<<"Detection point not reached after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
978 levelSetFilter->SetInput(levelSetFilter->GetOutput());
979 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose);
982 // std::ostringstream number_str;
983 // number_str << number;
984 // std::string param = number_str.str();
985 // system((script+param).c_str());
989 if ( (totalNumberOfIterations> 10000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS1_arg) ) break;
993 if (m_Verbose) std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
994 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
995 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
998 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
999 thresholder->SetUpperThreshold( 0.0 );
1000 thresholder->SetOutsideValue( 0 );
1001 thresholder->SetInsideValue( 1 );
1002 thresholder->SetInput( levelSetFilter->GetOutput() );
1003 if (m_Verbose) std::cout<<"Thresholding the output level set..."<<std::endl;
1004 thresholder->Update();
1005 grownEllips=thresholder->GetOutput();
1008 //---------------------------------
1009 // Write the grown ellips
1010 //---------------------------------
1011 if (m_ArgsInfo.writeGrownEllips_given) {
1012 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1013 caster->SetInput(grownEllips);
1014 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeGrownEllips_arg, m_Verbose);
1018 //----------------------------------------------------------------------------------------------------
1020 //----------------------------------------------------------------------------------------------------
1021 typename InternalImageType::Pointer filledRibs;
1022 if (m_ArgsInfo.filledRibs_given) {
1023 typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
1024 featureReader->SetFileName(m_ArgsInfo.filledRibs_arg);
1025 if (m_Verbose) std::cout<<"Reading filled ribs mask..."<<std::endl;
1026 featureReader->Update();
1027 filledRibs=featureReader->GetOutput();
1030 std::cout<<std::endl;
1031 std::cout<<"=========================================="<<std::endl;
1032 std::cout<<"|| Filling the ribs image ||"<<std::endl;
1033 std::cout<<"=========================================="<<std::endl;
1035 //---------------------------------
1036 // Make feature image air+bones
1037 //---------------------------------
1038 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
1039 setBackgroundFilter->SetInput(air_low);
1040 setBackgroundFilter->SetInput2(bones_low);
1041 setBackgroundFilter->SetMaskValue(0);
1042 setBackgroundFilter->SetOutsideValue(0);
1043 setBackgroundFilter->Update();
1044 typename InternalImageType::Pointer bones_air_low =setBackgroundFilter->GetOutput();
1046 //---------------------------------
1047 // Resampling previous solution
1048 //---------------------------------
1049 if (m_Verbose) std::cout<<"Resampling previous solution..."<<std::endl;
1050 typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
1051 typename InternalImageType::PointType origin;
1052 bones_low->TransformIndexToPhysicalPoint(bones_low->GetLargestPossibleRegion().GetIndex(), origin);
1053 resampler->SetOutputOrigin(origin);
1054 resampler->SetOutputSpacing(bones_low->GetSpacing());
1055 typename InternalImageType::SizeType size_low= bones_low->GetLargestPossibleRegion().GetSize();
1056 resampler->SetSize(size_low);
1057 typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
1058 resampler->SetInterpolator(interpolator);
1059 resampler->SetInput(grownEllips);
1060 resampler->Update();
1061 typename InternalImageType::Pointer grownEllips =resampler->GetOutput();
1064 //---------------------------------
1065 // Calculate Distance Map
1066 //---------------------------------
1067 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
1068 distanceMapImageFilter->SetInput(grownEllips);
1069 distanceMapImageFilter->SetInsideIsPositive(false);
1070 if (m_Verbose) std::cout<<"Calculating distance map..."<<std::endl;
1071 distanceMapImageFilter->Update();
1073 //---------------------------------
1074 // Grow while monitoring lung volume coverage
1075 //---------------------------------
1076 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
1077 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
1078 levelSetFilter->SetFeatureImage( bones_air_low );
1079 levelSetFilter->SetPropagationScaling( 1 );
1080 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve2_arg );
1081 levelSetFilter->SetAdvectionScaling( 0 );
1082 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS2_arg );
1083 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter2_arg );
1084 levelSetFilter->SetUseImageSpacing(true);
1086 //---------------------------------
1088 //---------------------------------
1089 typename InversionFilterType::Pointer invertor= InversionFilterType::New();
1090 invertor->SetInput(lungs_low);
1091 invertor->SetLowerThreshold(0);
1092 invertor->SetUpperThreshold(0);
1093 invertor->SetInsideValue(1);
1095 typename InternalImageType::Pointer lungs_low_inv=invertor->GetOutput();
1097 //---------------------------------
1098 // Calculate lung volume
1099 //---------------------------------
1100 typename LabelStatisticsImageFilterType::Pointer statisticsFilter= LabelStatisticsImageFilterType::New();
1101 statisticsFilter->SetInput(lungs_low_inv);
1102 statisticsFilter->SetLabelInput(lungs_low);
1103 statisticsFilter->Update();
1104 unsigned int volume=statisticsFilter->GetSum(0);
1106 // Prepare threshold
1107 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1108 thresholder->SetUpperThreshold( 0.0 );
1109 thresholder->SetOutsideValue( 0 );
1110 thresholder->SetInsideValue( 1 );
1114 unsigned int totalNumberOfIterations=0;
1115 unsigned int coverage=0;
1118 // //---------------------------------
1119 // // Script for making movie
1120 // //---------------------------------
1121 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1124 levelSetFilter->Update();
1125 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1127 thresholder->SetInput( levelSetFilter->GetOutput() );
1128 thresholder->Update();
1129 statisticsFilter->SetInput(thresholder->GetOutput());
1130 statisticsFilter->Update();
1131 coverage=statisticsFilter->GetSum(0);
1133 // Compare the volumes
1134 if ( coverage >= m_ArgsInfo.fillingLevel_arg * volume/100 ) {
1135 if (m_Verbose) std::cout<<"Lungs filled for "<< m_ArgsInfo.fillingLevel_arg<<"% !"<<std::endl;
1138 if (m_Verbose) std::cout<<"After "<<totalNumberOfIterations<<" iterations, lungs are covered for "
1139 <<(double)coverage/volume*100.0<<"%..."<<std::endl;
1140 levelSetFilter->SetInput(levelSetFilter->GetOutput());
1141 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose);
1144 // std::ostringstream number_str;
1145 // number_str << number;
1146 // std::string param = number_str.str();
1147 // system((script+param).c_str());
1151 if ( (totalNumberOfIterations> 30000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS2_arg) ) break;
1155 if (m_Verbose) std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1156 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1157 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1160 thresholder->SetInput( levelSetFilter->GetOutput() );
1161 thresholder->Update();
1162 filledRibs=thresholder->GetOutput();
1163 // writeImage<InternalImageType>(filledRibs,"/home/jef/tmp/filled_ribs_before.mhd", m_Verbose);
1167 //---------------------------------
1168 // Write the filled ribs
1169 //---------------------------------
1170 if (m_ArgsInfo.writeFilledRibs_given) {
1171 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1172 caster->SetInput(filledRibs);
1173 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeFilledRibs_arg, m_Verbose);
1177 //----------------------------------------------------------------------------------------------------
1178 // Collapse to the lungs
1179 //----------------------------------------------------------------------------------------------------
1181 std::cout<<std::endl;
1182 std::cout<<"=========================================="<<std::endl;
1183 std::cout<<"|| Collapsing to the lung image ||"<<std::endl;
1184 std::cout<<"=========================================="<<std::endl;
1187 //---------------------------------
1188 // Make feature image air+bones
1189 //---------------------------------
1190 if (m_Verbose) std::cout<<"Making feature images..."<<std::endl;
1191 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
1192 setBackgroundFilter->SetInput(air);
1193 setBackgroundFilter->SetInput2(bones);
1194 setBackgroundFilter->SetMaskValue(0);
1195 setBackgroundFilter->SetOutsideValue(0);
1196 setBackgroundFilter->Update();
1197 typename InternalImageType::Pointer bones_air =setBackgroundFilter->GetOutput();
1199 typename SetBackgroundFilterType::Pointer setBackgroundFilter2 = SetBackgroundFilterType::New();
1200 setBackgroundFilter2->SetInput(bones_air);
1201 setBackgroundFilter2->SetInput2(lungs);
1202 setBackgroundFilter2->SetMaskValue(0);
1203 setBackgroundFilter2->SetOutsideValue(0);
1204 setBackgroundFilter2->Update();
1205 typename InternalImageType::Pointer lungs_bones_air =setBackgroundFilter2->GetOutput();
1207 //---------------------------------
1208 // Prepare previous solution
1209 //---------------------------------
1210 if (m_Verbose) std::cout<<"Resampling previous solution..."<<std::endl;
1211 typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
1212 typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
1213 typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
1214 resampler->SetOutputStartIndex(bones->GetLargestPossibleRegion().GetIndex());
1215 resampler->SetInput(filledRibs);
1216 resampler->SetInterpolator(interpolator);
1217 resampler->SetOutputParametersFromImage(bones);
1218 resampler->Update();
1219 filledRibs =resampler->GetOutput();
1221 if (m_Verbose) std::cout<<"Setting lungs to 1..."<<std::endl;
1222 typename SetBackgroundFilterType::Pointer setBackgroundFilter3 = SetBackgroundFilterType::New();
1223 setBackgroundFilter3->SetInput(filledRibs);
1224 setBackgroundFilter3->SetInput2(lungs);
1225 setBackgroundFilter3->SetMaskValue(0);
1226 setBackgroundFilter3->SetOutsideValue(1);
1227 setBackgroundFilter3->Update();
1228 filledRibs=setBackgroundFilter3->GetOutput();
1230 if (m_Verbose) std::cout<<"Removing overlap with bones..."<<std::endl;
1231 typename SetBackgroundFilterType::Pointer setBackgroundFilter4 = SetBackgroundFilterType::New();
1232 setBackgroundFilter4->SetInput(filledRibs);
1233 setBackgroundFilter4->SetInput2(bones);
1234 setBackgroundFilter4->SetMaskValue(0);
1235 setBackgroundFilter4->SetOutsideValue(0);
1236 setBackgroundFilter4->Update();
1237 filledRibs =setBackgroundFilter4->GetOutput();
1238 // writeImage<InternalImageType>(filledRibs,"/home/jef/tmp/filledRibs_pp.mhd");
1239 //---------------------------------
1240 // Calculate Distance Map
1241 //---------------------------------
1242 // typedef itk::ApproximateSignedDistanceMapImageFilter <InternalImageType, LevelSetImageType> DistanceMapImageFilterType2;
1243 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
1244 distanceMapImageFilter->SetInput(filledRibs);
1245 distanceMapImageFilter->SetInsideIsPositive(false);
1246 // distanceMapImageFilter->SetInsideValue(0);
1247 // distanceMapImageFilter->SetOutsideValue(1);
1248 if (m_Verbose) std::cout<<"Calculating distance map..."<<std::endl;
1249 distanceMapImageFilter->Update();
1251 //---------------------------------
1253 //---------------------------------
1254 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
1255 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
1256 levelSetFilter->SetFeatureImage( lungs_bones_air );
1257 levelSetFilter->SetPropagationScaling(m_ArgsInfo.prop3_arg);
1258 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve3_arg );
1259 levelSetFilter->SetAdvectionScaling( 0 );
1260 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS3_arg );
1261 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter3_arg );
1262 levelSetFilter->SetUseImageSpacing(true);
1265 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1268 if (m_Verbose) std::cout<<"Starting the levelset..."<<std::endl;
1269 unsigned int totalNumberOfIterations=0;
1271 levelSetFilter->Update();
1274 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1275 levelSetFilter->SetInput(levelSetFilter->GetOutput());
1276 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg);
1277 std::cout << "After "<<totalNumberOfIterations<<" iteration the RMS change is " << levelSetFilter->GetRMSChange() <<"..."<< std::endl;
1280 // std::ostringstream number_str;
1281 // number_str << number;
1282 // std::string param = number_str.str();
1283 // system((script+param).c_str());
1286 // stopping condition
1287 if ( (totalNumberOfIterations>= (unsigned int) m_ArgsInfo.maxIter3_arg) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS3_arg) ) break;
1288 levelSetFilter->SetNumberOfIterations( std::min ((unsigned int) m_ArgsInfo.iter3_arg, (unsigned int) m_ArgsInfo.maxIter3_arg-totalNumberOfIterations ) );
1293 std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1294 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1295 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1299 typedef itk::BinaryThresholdImageFilter< LevelSetImageType,InternalImageType> LevelSetBinarizeFilterType;
1300 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1301 thresholder->SetUpperThreshold( 0.0 );
1302 thresholder->SetOutsideValue( 0 );
1303 thresholder->SetInsideValue( 1 );
1304 thresholder->SetInput( levelSetFilter->GetOutput() );
1305 thresholder->Update();
1306 typename InternalImageType::Pointer output = thresholder->GetOutput();
1309 //----------------------------------------------------------------------------------------------------
1310 // Prepare the output
1311 //----------------------------------------------------------------------------------------------------
1313 //---------------------------------
1315 //---------------------------------
1316 if (m_Verbose) std::cout<<"Removing overlap mask with air..."<<std::endl;
1317 typename SetBackgroundFilterType::Pointer setBackgroundFilter5 = SetBackgroundFilterType::New();
1318 setBackgroundFilter5->SetInput(output);
1319 setBackgroundFilter5->SetInput2(air);
1320 setBackgroundFilter5->SetMaskValue(0);
1321 setBackgroundFilter5->SetOutsideValue(0);
1322 setBackgroundFilter5->Update();
1323 output=setBackgroundFilter5->GetOutput();
1325 //---------------------------------
1326 // Open & close to cleanup
1327 //---------------------------------
1328 if ( m_ArgsInfo.openClose_flag) {
1330 //---------------------------------
1331 // Structuring element
1332 //---------------------------------
1333 typedef itk::BinaryBallStructuringElement<InternalPixelType,Dimension > KernelType;
1334 KernelType structuringElement;
1335 structuringElement.SetRadius(1);
1336 structuringElement.CreateStructuringElement();
1338 //---------------------------------
1340 //---------------------------------
1341 typedef itk::BinaryMorphologicalOpeningImageFilter<InternalImageType, InternalImageType , KernelType> OpenFilterType;
1342 typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
1343 openFilter->SetInput(output);
1344 openFilter->SetBackgroundValue(0);
1345 openFilter->SetForegroundValue(1);
1346 openFilter->SetKernel(structuringElement);
1347 if(m_Verbose) std::cout<<"Opening the output image..."<<std::endl;
1349 //---------------------------------
1351 //---------------------------------
1352 typedef itk::BinaryMorphologicalClosingImageFilter<InternalImageType, InternalImageType , KernelType> CloseFilterType;
1353 typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
1354 closeFilter->SetInput(openFilter->GetOutput());
1355 closeFilter->SetSafeBorder(true);
1356 closeFilter->SetForegroundValue(1);
1357 closeFilter->SetKernel(structuringElement);
1358 if(m_Verbose) std::cout<<"Closing the output image..."<<std::endl;
1359 closeFilter->Update();
1360 output=closeFilter->GetOutput();
1363 //writeImage<InternalImageType>(output,"/home/jef/tmp/mm_double.mhd");
1365 // Extract the upper part
1366 typedef itk::CropImageFilter<InternalImageType, InternalImageType> CropImageFilterType;
1367 typename CropImageFilterType::Pointer cropFilter=CropImageFilterType::New();
1368 cropFilter->SetInput(output);
1369 typename InputImageType::SizeType lSize, uSize;
1372 lSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
1373 cropFilter->SetLowerBoundaryCropSize(lSize);
1374 cropFilter->SetUpperBoundaryCropSize(uSize);
1375 cropFilter->Update();
1378 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1379 caster->SetInput(cropFilter->GetOutput());
1380 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.output_arg, m_Verbose);
1386 #endif //#define clitkMotionMaskGenericFilter_txx