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;
380 connectFilter->Update();
381 if (m_Verbose) std::cout<<"found "<< connectFilter->GetObjectCount() << std::endl;
383 //---------------------------------
384 // Sort the labels according to size
385 //---------------------------------
386 typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
387 relabelFilter->SetInput(connectFilter->GetOutput());
388 if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
389 // writeImage<InternalImageType> (relabelFilter->GetOutput(), "/home/vdelmon/tmp/labels.mhd");
391 //---------------------------------
393 //---------------------------------
394 typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
395 thresholdFilter->SetInput(relabelFilter->GetOutput());
396 thresholdFilter->SetLower(1);
397 thresholdFilter->SetUpper(2);
398 if (m_Verbose) std::cout<<"Keeping the first two labels..."<<std::endl;
399 thresholdFilter->Update();
400 lungs=thresholdFilter->GetOutput();
405 //---------------------------------
407 //---------------------------------
408 typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
409 inversionFilter->SetInput(lungs);
410 inversionFilter->SetLowerThreshold(0);
411 inversionFilter->SetUpperThreshold(0);
412 inversionFilter->SetInsideValue(1);
413 inversionFilter->Update();
415 //---------------------------------
417 //---------------------------------
418 typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
419 mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
420 typename InternalImageType::SizeType padSize;
422 padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
423 mirrorPadImageFilter->SetPadLowerBound(padSize);
424 if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
425 mirrorPadImageFilter->Update();
426 lungs=mirrorPadImageFilter->GetOutput();
427 // writeImage<InternalImageType>(lungs,"/home/jef/tmp/lungs.mhd");
432 //-------------------------------------------------------------------
434 //-------------------------------------------------------------------
435 template <unsigned int Dimension, class PixelType>
436 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
437 MotionMaskGenericFilter::Resample( typename itk::Image<InternalPixelType,Dimension>::Pointer input )
440 typedef int InternalPixelType;
441 typedef itk::Image<PixelType, Dimension> InputImageType;
442 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
444 //---------------------------------
445 // Resample to new spacing
446 //---------------------------------
447 typename InternalImageType::SizeType size_low;
448 typename InternalImageType::SpacingType spacing_low;
449 for (unsigned int i=0; i<Dimension; i++)
450 if (m_ArgsInfo.spacing_given==Dimension)
451 for (unsigned int i=0; i<Dimension; i++) {
452 spacing_low[i]=m_ArgsInfo.spacing_arg[i];
453 size_low[i]=ceil(input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i]/spacing_low[i]);
456 for (unsigned int i=0; i<Dimension; i++) {
457 spacing_low[i]=m_ArgsInfo.spacing_arg[0];
458 size_low[i]=ceil(input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i]/spacing_low[i]);
461 typename InternalImageType::PointType origin;
462 input->TransformIndexToPhysicalPoint(input->GetLargestPossibleRegion().GetIndex(), origin);
463 typedef itk::ResampleImageFilter<InternalImageType, InternalImageType> ResampleImageFilterType;
464 typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
465 typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
466 typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
467 resampler->SetInterpolator(interpolator);
468 resampler->SetOutputOrigin(origin);
469 resampler->SetOutputSpacing(spacing_low);
470 resampler->SetSize(size_low);
471 resampler->SetInput(input);
473 typename InternalImageType::Pointer output =resampler->GetOutput();
478 //-------------------------------------------------------------------
480 //-------------------------------------------------------------------
481 template <unsigned int Dimension, class PixelType>
482 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
483 MotionMaskGenericFilter::InitializeEllips( typename itk::Vector<double,Dimension> center, typename itk::Image<InternalPixelType,Dimension>::Pointer bones_low )
487 typedef int InternalPixelType;
488 typedef itk::Image<PixelType, Dimension> InputImageType;
489 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
490 typedef itk::Image< unsigned char , Dimension> OutputImageType;
491 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
492 typedef clitk::SetBackgroundImageFilter<InternalImageType,InternalImageType, InternalImageType> SetBackgroundFilterType;
493 typedef itk::LabelStatisticsImageFilter<InternalImageType,InternalImageType> LabelStatisticsImageFilterType;
494 typedef itk::CastImageFilter<InternalImageType,OutputImageType> CastImageFilterType;
497 typename InternalImageType::Pointer ellips;
498 if (m_ArgsInfo.ellips_given || m_ArgsInfo.grownEllips_given || m_ArgsInfo.filledRibs_given) {
499 if(m_ArgsInfo.ellips_given) {
500 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
501 typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
502 featureReader->SetFileName(m_ArgsInfo.ellips_arg);
503 featureReader->Update();
504 ellips=featureReader->GetOutput();
508 std::cout<<std::endl;
509 std::cout<<"=========================================="<<std::endl;
510 std::cout<<"|| Initializing ellipsoide ||"<<std::endl;
511 std::cout<<"=========================================="<<std::endl;
514 //---------------------------------
515 // Offset from center
516 //---------------------------------
517 typename itk::Vector<double,Dimension> offset;
518 if(m_ArgsInfo.offset_given== Dimension) {
519 for(unsigned int i=0; i<Dimension; i++)
520 offset[i]=m_ArgsInfo.offset_arg[i];
525 itk::Vector<double,Dimension> centerEllips=center+offset;
527 std::cout<<"Placing the center of the initial ellipsoide at (mm) "<<std::endl;
528 std::cout<<centerEllips[0];
529 for (unsigned int i=1; i<Dimension; i++)
530 std::cout<<" "<<centerEllips[i];
531 std::cout<<std::endl;
534 //---------------------------------
536 //---------------------------------
537 ellips=InternalImageType::New();
538 ellips->SetRegions (bones_low->GetLargestPossibleRegion());
539 ellips->SetOrigin(bones_low->GetOrigin());
540 ellips->SetSpacing(bones_low->GetSpacing());
542 ellips->FillBuffer(0);
545 typename itk::Vector<double, Dimension> axes;
546 if (m_ArgsInfo.axes_given==Dimension)
547 for(unsigned int i=0; i<Dimension; i++)
548 axes[i]=m_ArgsInfo.axes_arg[i];
556 IteratorType itEllips(ellips, ellips->GetLargestPossibleRegion());
557 itEllips.GoToBegin();
558 typename InternalImageType::PointType point;
559 typename InternalImageType::IndexType index;
562 if (m_Verbose) std::cout<<"Drawing the initial ellipsoide..."<<std::endl;
563 while (!itEllips.IsAtEnd()) {
564 index=itEllips.GetIndex();
565 ellips->TransformIndexToPhysicalPoint(index, point);
567 for(unsigned int i=0; i<Dimension; i++)
568 distance+=powf( ( (centerEllips[i]-point[i])/axes[i] ), 2);
577 //---------------------------------
579 //---------------------------------
580 if (m_ArgsInfo.writeEllips_given) {
581 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
582 caster->SetInput(ellips);
583 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeEllips_arg, m_Verbose);
591 //-------------------------------------------------------------------
592 // Update with the number of dimensions and the pixeltype
593 //-------------------------------------------------------------------
594 template <unsigned int Dimension, class PixelType>
596 MotionMaskGenericFilter::UpdateWithDimAndPixelType()
600 typedef int InternalPixelType;
601 typedef itk::Image<PixelType, Dimension> InputImageType;
602 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
603 typedef itk::Image< float, Dimension> LevelSetImageType; //labeling doesn't work with unsigned char?
604 typedef itk::Image<unsigned char, Dimension> OutputImageType;
607 //----------------------------------------------------------------------------------------------------
608 // Typedef for Processing
609 //----------------------------------------------------------------------------------------------------
610 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
611 typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
612 typedef itk::BinaryThresholdImageFilter< LevelSetImageType,InternalImageType> LevelSetBinarizeFilterType;
613 typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> InversionFilterType;
614 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdFilterType;
615 typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
616 typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
617 typedef itk::MirrorPadImageFilter<InternalImageType, InternalImageType> MirrorPadImageFilterType;
618 typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
619 typedef itk::ResampleImageFilter<InternalImageType, InternalImageType> ResampleImageFilterType;
620 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
621 typedef itk::GeodesicActiveContourLevelSetImageFilter<LevelSetImageType, InternalImageType> LevelSetImageFilterType;
622 typedef itk::SignedDanielssonDistanceMapImageFilter<InternalImageType, LevelSetImageType> DistanceMapImageFilterType;
623 typedef clitk::SetBackgroundImageFilter<InternalImageType,InternalImageType, InternalImageType> SetBackgroundFilterType;
624 typedef itk::LabelStatisticsImageFilter<InternalImageType,InternalImageType> LabelStatisticsImageFilterType;
625 typedef itk::CastImageFilter<InternalImageType,OutputImageType> CastImageFilterType;
629 typedef itk::ImageFileReader<InputImageType> InputReaderType;
630 typename InputReaderType::Pointer reader = InputReaderType::New();
631 reader->SetFileName( m_InputFileName);
633 typename InputImageType::Pointer input= reader->GetOutput();
635 // // globals for avi
636 // unsigned int number=0;
640 std::cout<<std::endl;
641 std::cout<<"=========================================="<<std::endl;
642 std::cout<<"|| Making feature images ||"<<std::endl;
643 std::cout<<"=========================================="<<std::endl;
646 //--------------------------------------------------------------------------------
648 //-------------------------------------------------------------------------------
649 typename InternalImageType::Pointer lungs = this->GetLungsImage<Dimension, PixelType>(input);
651 //-------------------------------------------------------------------------------
653 //-------------------------------------------------------------------------------
654 typename InternalImageType::Pointer air = this->GetAirImage<Dimension, PixelType>(input, lungs);
656 //-------------------------------------------------------------------------------
658 //-------------------------------------------------------------------------------
659 typename InternalImageType::Pointer bones = this->GetBonesImage<Dimension, PixelType>(input);
661 //----------------------------------------------------------------------------------------------------
662 // Write feature images
663 //----------------------------------------------------------------------------------------------------
664 if(m_ArgsInfo.writeFeature_given) {
665 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
666 setBackgroundFilter->SetInput(air);
667 setBackgroundFilter->SetInput2(bones);
668 setBackgroundFilter->SetMaskValue(0);
669 setBackgroundFilter->SetOutsideValue(2);
670 setBackgroundFilter->Update();
671 typename InternalImageType::Pointer bones_air =setBackgroundFilter->GetOutput();
673 typename SetBackgroundFilterType::Pointer setBackgroundFilter2 = SetBackgroundFilterType::New();
674 setBackgroundFilter2->SetInput(bones_air);
675 setBackgroundFilter2->SetInput2(lungs);
676 setBackgroundFilter2->SetMaskValue(0);
677 setBackgroundFilter2->SetOutsideValue(3);
678 setBackgroundFilter2->Update();
679 typename InternalImageType::Pointer lungs_bones_air =setBackgroundFilter2->GetOutput();
681 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
682 caster->SetInput(lungs_bones_air);
683 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeFeature_arg, m_Verbose);
686 //----------------------------------------------------------------------------------------------------
687 // Low dimensional versions
688 //----------------------------------------------------------------------------------------------------
689 typename InternalImageType::Pointer bones_low =Resample<Dimension , PixelType>(bones);
690 typename InternalImageType::Pointer lungs_low =Resample<Dimension , PixelType>(lungs);
691 typename InternalImageType::Pointer air_low =Resample<Dimension , PixelType>(air);
693 //---------------------------------
695 //---------------------------------
696 if(m_ArgsInfo.pad_flag) {
697 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
698 IteratorType it(air_low, air_low->GetLargestPossibleRegion());
699 typename InternalImageType::IndexType index;
700 while(!it.IsAtEnd()) {
702 for (unsigned int i=0; i<Dimension; i++)
703 if(index[i]==air_low->GetLargestPossibleRegion().GetIndex()[i]
704 || index[i]==(unsigned int)air_low->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air_low->GetLargestPossibleRegion().GetSize()[i]-1)
710 //---------------------------------
712 //---------------------------------
713 typename itk::Vector<double,Dimension> center;
714 typedef itk::ImageMomentsCalculator<InternalImageType> MomentsCalculatorType;
715 typename MomentsCalculatorType::Pointer momentsCalculator=MomentsCalculatorType::New();
716 momentsCalculator->SetImage(air);
717 momentsCalculator->Compute();
718 center=momentsCalculator->GetCenterOfGravity();
720 std::cout<<"The center of gravity of the patient body is located at (mm) "<<std::endl;
721 std::cout<<center[0];
722 for (unsigned int i=1; i<Dimension; i++)
723 std::cout<<" "<<center[i];
724 std::cout<<std::endl;
727 //----------------------------------------------------------------------------------------------------
728 // Ellipsoide initialization
729 //----------------------------------------------------------------------------------------------------
730 typename InternalImageType::Pointer m_Ellips= InitializeEllips<Dimension, PixelType>(center,bones_low);
732 //----------------------------------------------------------------------------------------------------
734 //----------------------------------------------------------------------------------------------------
735 typename InternalImageType::Pointer grownEllips;
736 if (m_ArgsInfo.grownEllips_given || m_ArgsInfo.filledRibs_given) {
737 if (m_ArgsInfo.grownEllips_given) {
739 typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
740 featureReader->SetFileName(m_ArgsInfo.grownEllips_arg);
741 featureReader->Update();
742 grownEllips=featureReader->GetOutput();
747 std::cout<<std::endl;
748 std::cout<<"=========================================="<<std::endl;
749 std::cout<<"|| Growing ellipsoide ||"<<std::endl;
750 std::cout<<"=========================================="<<std::endl;
753 //---------------------------------
755 //---------------------------------
756 typename InternalImageType::PointType dPoint;
757 if (m_ArgsInfo.detectionPoint_given) {
758 for (unsigned int i=0; i<Dimension; i++)
759 dPoint[i]=m_ArgsInfo.detectionPoint_arg[i];
761 typename InternalImageType::RegionType searchRegion=air->GetLargestPossibleRegion();
762 typename InternalImageType::SizeType searchRegionSize = searchRegion.GetSize();
763 typename InternalImageType::IndexType searchRegionIndex = searchRegion.GetIndex();
764 searchRegionIndex[2]+=searchRegionSize[2]/2;
765 searchRegionSize[2]=1;
766 searchRegion.SetSize(searchRegionSize);
767 searchRegion.SetIndex(searchRegionIndex);
768 IteratorType itAbdomen(air, searchRegion);
769 itAbdomen.GoToBegin();
771 typename InternalImageType::PointType aPoint;
772 typename InternalImageType::IndexType aIndex;
774 if (m_Verbose) std::cout<<"Detecting the abdomen..."<<std::endl;
775 while (!itAbdomen.IsAtEnd()) {
776 if(itAbdomen.Get()) break;
779 aIndex=itAbdomen.GetIndex();
780 air->TransformIndexToPhysicalPoint(aIndex,aPoint);
781 if (m_Verbose) std::cout<<"Detected the abdomen at "<<aPoint<<".."<<std::endl;
784 //---------------------------------
785 // Detect abdomen in additional images?
786 //---------------------------------
787 if (m_ArgsInfo.detectionPairs_given) {
788 for (unsigned int i=0; i<m_ArgsInfo.detectionPairs_given; i++) {
789 typename InternalImageType::Pointer airAdd;
790 //---------------------------------
792 //--------------------------------
793 typedef itk::ImageFileReader<InputImageType> InputReaderType;
794 typename InputReaderType::Pointer reader = InputReaderType::New();
795 reader->SetFileName( m_ArgsInfo.detectionPairs_arg[i]);
797 typename InputImageType::Pointer additional= reader->GetOutput();
798 if (m_Verbose) std::cout<<"Extracting the air from additional image "<< i<<"..."<<std::endl;
800 //---------------------------------
801 // Binarize the image
802 //---------------------------------
803 typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
804 binarizeFilter->SetInput(additional);
805 binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdAir_arg));
806 binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdAir_arg));
807 if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<<m_ArgsInfo.lowerThresholdAir_arg
808 <<", "<<m_ArgsInfo.upperThresholdAir_arg<<"..."<<std::endl;
810 //---------------------------------
811 // Label the connected components
812 //---------------------------------
813 typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
814 connectFilter->SetInput(binarizeFilter->GetOutput());
815 connectFilter->SetBackgroundValue(0);
816 connectFilter->SetFullyConnected(false);
817 if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
819 //---------------------------------
820 // Sort the labels according to size
821 //---------------------------------
822 typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
823 relabelFilter->SetInput(connectFilter->GetOutput());
824 if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
826 //---------------------------------
828 //---------------------------------
829 typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
830 thresholdFilter->SetInput(relabelFilter->GetOutput());
831 thresholdFilter->SetUpper(1);
832 if (m_Verbose) std::cout<<"Keeping the first label..."<<std::endl;
833 thresholdFilter->Update();
834 airAdd=thresholdFilter->GetOutput();
837 //---------------------------------
839 //---------------------------------
840 typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
841 inversionFilter->SetInput(airAdd);
842 inversionFilter->SetLowerThreshold(0);
843 inversionFilter->SetUpperThreshold(0);
844 inversionFilter->SetInsideValue(1);
845 inversionFilter->Update();
847 //---------------------------------
849 //---------------------------------
850 typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
851 mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
852 typename InternalImageType::SizeType padSize;
854 padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
855 mirrorPadImageFilter->SetPadLowerBound(padSize);
856 if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
857 mirrorPadImageFilter->Update();
858 airAdd=mirrorPadImageFilter->GetOutput();
860 //---------------------------------
862 //---------------------------------
863 typename InternalImageType::RegionType searchRegion=airAdd->GetLargestPossibleRegion();
864 typename InternalImageType::SizeType searchRegionSize = searchRegion.GetSize();
865 typename InternalImageType::IndexType searchRegionIndex = searchRegion.GetIndex();
866 searchRegionIndex[2]+=searchRegionSize[2]/2;
867 searchRegionSize[2]=1;
868 searchRegion.SetSize(searchRegionSize);
869 searchRegion.SetIndex(searchRegionIndex);
870 IteratorType itAbdomen(airAdd, searchRegion);
871 itAbdomen.GoToBegin();
873 typename InternalImageType::PointType additionalPoint;
874 typename InternalImageType::IndexType aIndex;
876 if (m_Verbose) std::cout<<"Detecting the abdomen..."<<std::endl;
877 while (!itAbdomen.IsAtEnd()) {
878 if(itAbdomen.Get()) break;
881 aIndex=itAbdomen.GetIndex();
882 airAdd->TransformIndexToPhysicalPoint(aIndex,additionalPoint);
883 if (m_Verbose) std::cout<<"Detected the abdomen in the additional image at "<<additionalPoint<<".."<<std::endl;
885 if(additionalPoint[1]< aPoint[1]) {
886 aPoint=additionalPoint;
887 if (m_Verbose) std::cout<<"Modifying the detected abdomen to "<<aPoint<<".."<<std::endl;
894 // Determine the detection point
898 if(m_ArgsInfo.offsetDetect_given==Dimension)
899 for(unsigned int i=0; i <Dimension; i++)
900 dPoint[i]+=m_ArgsInfo.offsetDetect_arg[i];
905 if (m_Verbose) std::cout<<"Setting the detection point to "<<dPoint<<".."<<std::endl;
908 //---------------------------------
909 // Pad the rib image and ellips image
910 //---------------------------------
911 typename InternalImageType::Pointer padded_ellips;
912 typename InternalImageType::Pointer padded_bones_low;
914 // If detection point not inside the image: pad
915 typename InternalImageType::IndexType dIndex;
916 if (!bones_low->TransformPhysicalPointToIndex(dPoint, dIndex)) {
917 typename InternalImageType::SizeType padSize;
919 padSize[1]=abs(dIndex[1])+1;
920 if (m_Verbose) std::cout<<"Padding the images with padding size "<<padSize<<" to include the detection point..."<<std::endl;
922 typename MirrorPadImageFilterType::Pointer padBonesFilter=MirrorPadImageFilterType::New();
923 padBonesFilter->SetInput(bones_low);
924 padBonesFilter->SetPadLowerBound(padSize);
925 padBonesFilter->Update();
926 padded_bones_low=padBonesFilter->GetOutput();
928 typename MirrorPadImageFilterType::Pointer padEllipsFilter=MirrorPadImageFilterType::New();
929 padEllipsFilter->SetInput(m_Ellips);
930 padEllipsFilter->SetPadLowerBound(padSize);
931 padEllipsFilter->Update();
932 padded_ellips=padEllipsFilter->GetOutput();
936 padded_bones_low=bones_low;
937 padded_ellips=m_Ellips;
941 //---------------------------------
942 // Calculate distance map
943 //---------------------------------
944 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
945 distanceMapImageFilter->SetInput(padded_ellips);
946 distanceMapImageFilter->SetInsideIsPositive(false);
947 if (m_Verbose) std::cout<<"Calculating the distance map..."<<std::endl;
948 distanceMapImageFilter->Update();
950 //---------------------------------
951 // Grow while monitoring detection point
952 //---------------------------------
953 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
954 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
955 levelSetFilter->SetFeatureImage( padded_bones_low );
956 levelSetFilter->SetPropagationScaling( 1 );
957 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve1_arg );
958 levelSetFilter->SetAdvectionScaling( 0 );
959 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS1_arg );
960 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter1_arg );
961 levelSetFilter->SetUseImageSpacing(true);
963 // //---------------------------------
964 // // Script for making movie
965 // //---------------------------------
966 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
969 if (m_Verbose) std::cout<<"Starting level set segmentation..."<<std::endl;
970 unsigned int totalNumberOfIterations=0;
972 levelSetFilter->Update();
973 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
975 if ( levelSetFilter->GetOutput()->GetPixel(dIndex) < 0 ) {
976 if (m_Verbose) std::cout<<"Detection point reached!"<<std::endl;
979 if (m_Verbose) std::cout<<"Detection point not reached after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
980 levelSetFilter->SetInput(levelSetFilter->GetOutput());
981 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose);
984 // std::ostringstream number_str;
985 // number_str << number;
986 // std::string param = number_str.str();
987 // system((script+param).c_str());
991 if ( (totalNumberOfIterations> 10000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS1_arg) ) break;
995 if (m_Verbose) std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
996 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
997 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1000 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1001 thresholder->SetUpperThreshold( 0.0 );
1002 thresholder->SetOutsideValue( 0 );
1003 thresholder->SetInsideValue( 1 );
1004 thresholder->SetInput( levelSetFilter->GetOutput() );
1005 if (m_Verbose) std::cout<<"Thresholding the output level set..."<<std::endl;
1006 thresholder->Update();
1007 grownEllips=thresholder->GetOutput();
1010 //---------------------------------
1011 // Write the grown ellips
1012 //---------------------------------
1013 if (m_ArgsInfo.writeGrownEllips_given) {
1014 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1015 caster->SetInput(grownEllips);
1016 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeGrownEllips_arg, m_Verbose);
1020 //----------------------------------------------------------------------------------------------------
1022 //----------------------------------------------------------------------------------------------------
1023 typename InternalImageType::Pointer filledRibs;
1024 if (m_ArgsInfo.filledRibs_given) {
1025 typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
1026 featureReader->SetFileName(m_ArgsInfo.filledRibs_arg);
1027 if (m_Verbose) std::cout<<"Reading filled ribs mask..."<<std::endl;
1028 featureReader->Update();
1029 filledRibs=featureReader->GetOutput();
1032 std::cout<<std::endl;
1033 std::cout<<"=========================================="<<std::endl;
1034 std::cout<<"|| Filling the ribs image ||"<<std::endl;
1035 std::cout<<"=========================================="<<std::endl;
1037 //---------------------------------
1038 // Make feature image air+bones
1039 //---------------------------------
1040 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
1041 setBackgroundFilter->SetInput(air_low);
1042 setBackgroundFilter->SetInput2(bones_low);
1043 setBackgroundFilter->SetMaskValue(0);
1044 setBackgroundFilter->SetOutsideValue(0);
1045 setBackgroundFilter->Update();
1046 typename InternalImageType::Pointer bones_air_low =setBackgroundFilter->GetOutput();
1048 //---------------------------------
1049 // Resampling previous solution
1050 //---------------------------------
1051 if (m_Verbose) std::cout<<"Resampling previous solution..."<<std::endl;
1052 typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
1053 typename InternalImageType::PointType origin;
1054 bones_low->TransformIndexToPhysicalPoint(bones_low->GetLargestPossibleRegion().GetIndex(), origin);
1055 resampler->SetOutputOrigin(origin);
1056 resampler->SetOutputSpacing(bones_low->GetSpacing());
1057 typename InternalImageType::SizeType size_low= bones_low->GetLargestPossibleRegion().GetSize();
1058 resampler->SetSize(size_low);
1059 typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
1060 resampler->SetInterpolator(interpolator);
1061 resampler->SetInput(grownEllips);
1062 resampler->Update();
1063 typename InternalImageType::Pointer grownEllips =resampler->GetOutput();
1066 //---------------------------------
1067 // Calculate Distance Map
1068 //---------------------------------
1069 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
1070 distanceMapImageFilter->SetInput(grownEllips);
1071 distanceMapImageFilter->SetInsideIsPositive(false);
1072 if (m_Verbose) std::cout<<"Calculating distance map..."<<std::endl;
1073 distanceMapImageFilter->Update();
1075 //---------------------------------
1076 // Grow while monitoring lung volume coverage
1077 //---------------------------------
1078 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
1079 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
1080 levelSetFilter->SetFeatureImage( bones_air_low );
1081 levelSetFilter->SetPropagationScaling( 1 );
1082 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve2_arg );
1083 levelSetFilter->SetAdvectionScaling( 0 );
1084 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS2_arg );
1085 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter2_arg );
1086 levelSetFilter->SetUseImageSpacing(true);
1088 //---------------------------------
1090 //---------------------------------
1091 typename InversionFilterType::Pointer invertor= InversionFilterType::New();
1092 invertor->SetInput(lungs_low);
1093 invertor->SetLowerThreshold(0);
1094 invertor->SetUpperThreshold(0);
1095 invertor->SetInsideValue(1);
1097 typename InternalImageType::Pointer lungs_low_inv=invertor->GetOutput();
1099 //---------------------------------
1100 // Calculate lung volume
1101 //---------------------------------
1102 typename LabelStatisticsImageFilterType::Pointer statisticsFilter= LabelStatisticsImageFilterType::New();
1103 statisticsFilter->SetInput(lungs_low_inv);
1104 statisticsFilter->SetLabelInput(lungs_low);
1105 statisticsFilter->Update();
1106 unsigned int volume=statisticsFilter->GetSum(0);
1108 // Prepare threshold
1109 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1110 thresholder->SetUpperThreshold( 0.0 );
1111 thresholder->SetOutsideValue( 0 );
1112 thresholder->SetInsideValue( 1 );
1116 unsigned int totalNumberOfIterations=0;
1117 unsigned int coverage=0;
1120 // //---------------------------------
1121 // // Script for making movie
1122 // //---------------------------------
1123 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1126 levelSetFilter->Update();
1127 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1129 thresholder->SetInput( levelSetFilter->GetOutput() );
1130 thresholder->Update();
1131 statisticsFilter->SetInput(thresholder->GetOutput());
1132 statisticsFilter->Update();
1133 coverage=statisticsFilter->GetSum(0);
1135 // Compare the volumes
1136 if ( coverage >= m_ArgsInfo.fillingLevel_arg * volume/100 ) {
1137 if (m_Verbose) std::cout<<"Lungs filled for "<< m_ArgsInfo.fillingLevel_arg<<"% !"<<std::endl;
1140 if (m_Verbose) std::cout<<"After "<<totalNumberOfIterations<<" iterations, lungs are covered for "
1141 <<(double)coverage/volume*100.0<<"%..."<<std::endl;
1142 levelSetFilter->SetInput(levelSetFilter->GetOutput());
1143 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose);
1146 // std::ostringstream number_str;
1147 // number_str << number;
1148 // std::string param = number_str.str();
1149 // system((script+param).c_str());
1153 if ( (totalNumberOfIterations> 30000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS2_arg) ) break;
1157 if (m_Verbose) std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1158 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1159 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1162 thresholder->SetInput( levelSetFilter->GetOutput() );
1163 thresholder->Update();
1164 filledRibs=thresholder->GetOutput();
1165 // writeImage<InternalImageType>(filledRibs,"/home/jef/tmp/filled_ribs_before.mhd", m_Verbose);
1169 //---------------------------------
1170 // Write the filled ribs
1171 //---------------------------------
1172 if (m_ArgsInfo.writeFilledRibs_given) {
1173 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1174 caster->SetInput(filledRibs);
1175 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeFilledRibs_arg, m_Verbose);
1179 //----------------------------------------------------------------------------------------------------
1180 // Collapse to the lungs
1181 //----------------------------------------------------------------------------------------------------
1183 std::cout<<std::endl;
1184 std::cout<<"=========================================="<<std::endl;
1185 std::cout<<"|| Collapsing to the lung image ||"<<std::endl;
1186 std::cout<<"=========================================="<<std::endl;
1189 //---------------------------------
1190 // Make feature image air+bones
1191 //---------------------------------
1192 if (m_Verbose) std::cout<<"Making feature images..."<<std::endl;
1193 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
1194 setBackgroundFilter->SetInput(air);
1195 setBackgroundFilter->SetInput2(bones);
1196 setBackgroundFilter->SetMaskValue(0);
1197 setBackgroundFilter->SetOutsideValue(0);
1198 setBackgroundFilter->Update();
1199 typename InternalImageType::Pointer bones_air =setBackgroundFilter->GetOutput();
1201 typename SetBackgroundFilterType::Pointer setBackgroundFilter2 = SetBackgroundFilterType::New();
1202 setBackgroundFilter2->SetInput(bones_air);
1203 setBackgroundFilter2->SetInput2(lungs);
1204 setBackgroundFilter2->SetMaskValue(0);
1205 setBackgroundFilter2->SetOutsideValue(0);
1206 setBackgroundFilter2->Update();
1207 typename InternalImageType::Pointer lungs_bones_air =setBackgroundFilter2->GetOutput();
1209 //---------------------------------
1210 // Prepare previous solution
1211 //---------------------------------
1212 if (m_Verbose) std::cout<<"Resampling previous solution..."<<std::endl;
1213 typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
1214 typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
1215 typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
1216 resampler->SetOutputStartIndex(bones->GetLargestPossibleRegion().GetIndex());
1217 resampler->SetInput(filledRibs);
1218 resampler->SetInterpolator(interpolator);
1219 resampler->SetOutputParametersFromImage(bones);
1220 resampler->Update();
1221 filledRibs =resampler->GetOutput();
1223 if (m_Verbose) std::cout<<"Setting lungs to 1..."<<std::endl;
1224 typename SetBackgroundFilterType::Pointer setBackgroundFilter3 = SetBackgroundFilterType::New();
1225 setBackgroundFilter3->SetInput(filledRibs);
1226 setBackgroundFilter3->SetInput2(lungs);
1227 setBackgroundFilter3->SetMaskValue(0);
1228 setBackgroundFilter3->SetOutsideValue(1);
1229 setBackgroundFilter3->Update();
1230 filledRibs=setBackgroundFilter3->GetOutput();
1232 if (m_Verbose) std::cout<<"Removing overlap with bones..."<<std::endl;
1233 typename SetBackgroundFilterType::Pointer setBackgroundFilter4 = SetBackgroundFilterType::New();
1234 setBackgroundFilter4->SetInput(filledRibs);
1235 setBackgroundFilter4->SetInput2(bones);
1236 setBackgroundFilter4->SetMaskValue(0);
1237 setBackgroundFilter4->SetOutsideValue(0);
1238 setBackgroundFilter4->Update();
1239 filledRibs =setBackgroundFilter4->GetOutput();
1240 // writeImage<InternalImageType>(filledRibs,"/home/jef/tmp/filledRibs_pp.mhd");
1241 //---------------------------------
1242 // Calculate Distance Map
1243 //---------------------------------
1244 // typedef itk::ApproximateSignedDistanceMapImageFilter <InternalImageType, LevelSetImageType> DistanceMapImageFilterType2;
1245 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
1246 distanceMapImageFilter->SetInput(filledRibs);
1247 distanceMapImageFilter->SetInsideIsPositive(false);
1248 // distanceMapImageFilter->SetInsideValue(0);
1249 // distanceMapImageFilter->SetOutsideValue(1);
1250 if (m_Verbose) std::cout<<"Calculating distance map..."<<std::endl;
1251 distanceMapImageFilter->Update();
1253 //---------------------------------
1255 //---------------------------------
1256 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
1257 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
1258 levelSetFilter->SetFeatureImage( lungs_bones_air );
1259 levelSetFilter->SetPropagationScaling(m_ArgsInfo.prop3_arg);
1260 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve3_arg );
1261 levelSetFilter->SetAdvectionScaling( 0 );
1262 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS3_arg );
1263 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter3_arg );
1264 levelSetFilter->SetUseImageSpacing(true);
1267 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1270 if (m_Verbose) std::cout<<"Starting the levelset..."<<std::endl;
1271 unsigned int totalNumberOfIterations=0;
1273 levelSetFilter->Update();
1276 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1277 levelSetFilter->SetInput(levelSetFilter->GetOutput());
1278 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg);
1279 std::cout << "After "<<totalNumberOfIterations<<" iteration the RMS change is " << levelSetFilter->GetRMSChange() <<"..."<< std::endl;
1282 // std::ostringstream number_str;
1283 // number_str << number;
1284 // std::string param = number_str.str();
1285 // system((script+param).c_str());
1288 // stopping condition
1289 if ( (totalNumberOfIterations>= (unsigned int) m_ArgsInfo.maxIter3_arg) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS3_arg) ) break;
1290 levelSetFilter->SetNumberOfIterations( std::min ((unsigned int) m_ArgsInfo.iter3_arg, (unsigned int) m_ArgsInfo.maxIter3_arg-totalNumberOfIterations ) );
1295 std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1296 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1297 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1301 typedef itk::BinaryThresholdImageFilter< LevelSetImageType,InternalImageType> LevelSetBinarizeFilterType;
1302 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1303 thresholder->SetUpperThreshold( 0.0 );
1304 thresholder->SetOutsideValue( 0 );
1305 thresholder->SetInsideValue( 1 );
1306 thresholder->SetInput( levelSetFilter->GetOutput() );
1307 thresholder->Update();
1308 typename InternalImageType::Pointer output = thresholder->GetOutput();
1311 //----------------------------------------------------------------------------------------------------
1312 // Prepare the output
1313 //----------------------------------------------------------------------------------------------------
1315 //---------------------------------
1317 //---------------------------------
1318 if (m_Verbose) std::cout<<"Removing overlap mask with air..."<<std::endl;
1319 typename SetBackgroundFilterType::Pointer setBackgroundFilter5 = SetBackgroundFilterType::New();
1320 setBackgroundFilter5->SetInput(output);
1321 setBackgroundFilter5->SetInput2(air);
1322 setBackgroundFilter5->SetMaskValue(0);
1323 setBackgroundFilter5->SetOutsideValue(0);
1324 setBackgroundFilter5->Update();
1325 output=setBackgroundFilter5->GetOutput();
1327 //---------------------------------
1328 // Open & close to cleanup
1329 //---------------------------------
1330 if ( m_ArgsInfo.openClose_flag) {
1332 //---------------------------------
1333 // Structuring element
1334 //---------------------------------
1335 typedef itk::BinaryBallStructuringElement<InternalPixelType,Dimension > KernelType;
1336 KernelType structuringElement;
1337 structuringElement.SetRadius(1);
1338 structuringElement.CreateStructuringElement();
1340 //---------------------------------
1342 //---------------------------------
1343 typedef itk::BinaryMorphologicalOpeningImageFilter<InternalImageType, InternalImageType , KernelType> OpenFilterType;
1344 typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
1345 openFilter->SetInput(output);
1346 openFilter->SetBackgroundValue(0);
1347 openFilter->SetForegroundValue(1);
1348 openFilter->SetKernel(structuringElement);
1349 if(m_Verbose) std::cout<<"Opening the output image..."<<std::endl;
1351 //---------------------------------
1353 //---------------------------------
1354 typedef itk::BinaryMorphologicalClosingImageFilter<InternalImageType, InternalImageType , KernelType> CloseFilterType;
1355 typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
1356 closeFilter->SetInput(openFilter->GetOutput());
1357 closeFilter->SetSafeBorder(true);
1358 closeFilter->SetForegroundValue(1);
1359 closeFilter->SetKernel(structuringElement);
1360 if(m_Verbose) std::cout<<"Closing the output image..."<<std::endl;
1361 closeFilter->Update();
1362 output=closeFilter->GetOutput();
1365 //writeImage<InternalImageType>(output,"/home/jef/tmp/mm_double.mhd");
1367 // Extract the upper part
1368 typedef itk::CropImageFilter<InternalImageType, InternalImageType> CropImageFilterType;
1369 typename CropImageFilterType::Pointer cropFilter=CropImageFilterType::New();
1370 cropFilter->SetInput(output);
1371 typename InputImageType::SizeType lSize, uSize;
1374 lSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
1375 cropFilter->SetLowerBoundaryCropSize(lSize);
1376 cropFilter->SetUpperBoundaryCropSize(uSize);
1377 cropFilter->Update();
1380 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1381 caster->SetInput(cropFilter->GetOutput());
1382 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.output_arg, m_Verbose);
1388 #endif //#define clitkMotionMaskGenericFilter_txx