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();
949 if (m_ArgsInfo.writeDistMap_given) {
950 writeImage<LevelSetImageType>(distanceMapImageFilter->GetOutput(), m_ArgsInfo.writeDistMap_arg, m_Verbose);
954 //---------------------------------
955 // Grow while monitoring detection point
956 //---------------------------------
957 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
958 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
959 levelSetFilter->SetFeatureImage( padded_bones_low );
960 levelSetFilter->SetPropagationScaling( 1 );
961 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve1_arg );
962 levelSetFilter->SetAdvectionScaling( 0 );
963 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS1_arg );
964 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter1_arg );
965 levelSetFilter->SetUseImageSpacing(true);
967 // //---------------------------------
968 // // Script for making movie
969 // //---------------------------------
970 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
973 if (m_Verbose) std::cout<<"Starting level set segmentation..."<<std::endl;
974 unsigned int totalNumberOfIterations=0;
976 levelSetFilter->Update();
977 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
979 if ( levelSetFilter->GetOutput()->GetPixel(dIndex) < 0 ) {
980 if (m_Verbose) std::cout<<"Detection point reached!"<<std::endl;
983 if (m_Verbose) std::cout<<"Detection point not reached after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
984 levelSetFilter->SetInput(levelSetFilter->GetOutput());
985 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose);
988 // std::ostringstream number_str;
989 // number_str << number;
990 // std::string param = number_str.str();
991 // system((script+param).c_str());
995 if ( (totalNumberOfIterations> 10000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS1_arg) ) break;
999 if (m_Verbose) std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1000 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1001 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1004 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1005 thresholder->SetUpperThreshold( 0.0 );
1006 thresholder->SetOutsideValue( 0 );
1007 thresholder->SetInsideValue( 1 );
1008 thresholder->SetInput( levelSetFilter->GetOutput() );
1009 if (m_Verbose) std::cout<<"Thresholding the output level set..."<<std::endl;
1010 thresholder->Update();
1011 grownEllips=thresholder->GetOutput();
1014 //---------------------------------
1015 // Write the grown ellips
1016 //---------------------------------
1017 if (m_ArgsInfo.writeGrownEllips_given) {
1018 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1019 caster->SetInput(grownEllips);
1020 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeGrownEllips_arg, m_Verbose);
1024 //----------------------------------------------------------------------------------------------------
1026 //----------------------------------------------------------------------------------------------------
1027 typename InternalImageType::Pointer filledRibs;
1028 if (m_ArgsInfo.filledRibs_given) {
1029 typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
1030 featureReader->SetFileName(m_ArgsInfo.filledRibs_arg);
1031 if (m_Verbose) std::cout<<"Reading filled ribs mask..."<<std::endl;
1032 featureReader->Update();
1033 filledRibs=featureReader->GetOutput();
1036 std::cout<<std::endl;
1037 std::cout<<"=========================================="<<std::endl;
1038 std::cout<<"|| Filling the ribs image ||"<<std::endl;
1039 std::cout<<"=========================================="<<std::endl;
1041 //---------------------------------
1042 // Make feature image air+bones
1043 //---------------------------------
1044 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
1045 setBackgroundFilter->SetInput(air_low);
1046 setBackgroundFilter->SetInput2(bones_low);
1047 setBackgroundFilter->SetMaskValue(0);
1048 setBackgroundFilter->SetOutsideValue(0);
1049 setBackgroundFilter->Update();
1050 typename InternalImageType::Pointer bones_air_low =setBackgroundFilter->GetOutput();
1052 //---------------------------------
1053 // Resampling previous solution
1054 //---------------------------------
1055 if (m_Verbose) std::cout<<"Resampling previous solution..."<<std::endl;
1056 typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
1057 typename InternalImageType::PointType origin;
1058 bones_low->TransformIndexToPhysicalPoint(bones_low->GetLargestPossibleRegion().GetIndex(), origin);
1059 resampler->SetOutputOrigin(origin);
1060 resampler->SetOutputSpacing(bones_low->GetSpacing());
1061 typename InternalImageType::SizeType size_low= bones_low->GetLargestPossibleRegion().GetSize();
1062 resampler->SetSize(size_low);
1063 typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
1064 resampler->SetInterpolator(interpolator);
1065 resampler->SetInput(grownEllips);
1066 resampler->Update();
1067 typename InternalImageType::Pointer grownEllips =resampler->GetOutput();
1070 //---------------------------------
1071 // Calculate Distance Map
1072 //---------------------------------
1073 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
1074 distanceMapImageFilter->SetInput(grownEllips);
1075 distanceMapImageFilter->SetInsideIsPositive(false);
1076 if (m_Verbose) std::cout<<"Calculating distance map..."<<std::endl;
1077 distanceMapImageFilter->Update();
1079 //---------------------------------
1080 // Grow while monitoring lung volume coverage
1081 //---------------------------------
1082 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
1083 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
1084 levelSetFilter->SetFeatureImage( bones_air_low );
1085 levelSetFilter->SetPropagationScaling( 1 );
1086 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve2_arg );
1087 levelSetFilter->SetAdvectionScaling( 0 );
1088 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS2_arg );
1089 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter2_arg );
1090 levelSetFilter->SetUseImageSpacing(true);
1092 //---------------------------------
1094 //---------------------------------
1095 typename InversionFilterType::Pointer invertor= InversionFilterType::New();
1096 invertor->SetInput(lungs_low);
1097 invertor->SetLowerThreshold(0);
1098 invertor->SetUpperThreshold(0);
1099 invertor->SetInsideValue(1);
1101 typename InternalImageType::Pointer lungs_low_inv=invertor->GetOutput();
1103 //---------------------------------
1104 // Calculate lung volume
1105 //---------------------------------
1106 typename LabelStatisticsImageFilterType::Pointer statisticsFilter= LabelStatisticsImageFilterType::New();
1107 statisticsFilter->SetInput(lungs_low_inv);
1108 statisticsFilter->SetLabelInput(lungs_low);
1109 statisticsFilter->Update();
1110 unsigned int volume=statisticsFilter->GetSum(0);
1112 // Prepare threshold
1113 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1114 thresholder->SetUpperThreshold( 0.0 );
1115 thresholder->SetOutsideValue( 0 );
1116 thresholder->SetInsideValue( 1 );
1120 unsigned int totalNumberOfIterations=0;
1121 unsigned int coverage=0;
1124 // //---------------------------------
1125 // // Script for making movie
1126 // //---------------------------------
1127 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1130 levelSetFilter->Update();
1131 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1133 thresholder->SetInput( levelSetFilter->GetOutput() );
1134 thresholder->Update();
1135 statisticsFilter->SetInput(thresholder->GetOutput());
1136 statisticsFilter->Update();
1137 coverage=statisticsFilter->GetSum(0);
1139 // Compare the volumes
1140 if ( coverage >= m_ArgsInfo.fillingLevel_arg * volume/100 ) {
1141 if (m_Verbose) std::cout<<"Lungs filled for "<< m_ArgsInfo.fillingLevel_arg<<"% !"<<std::endl;
1144 if (m_Verbose) std::cout<<"After "<<totalNumberOfIterations<<" iterations, lungs are covered for "
1145 <<(double)coverage/volume*100.0<<"%..."<<std::endl;
1146 levelSetFilter->SetInput(levelSetFilter->GetOutput());
1147 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose);
1150 // std::ostringstream number_str;
1151 // number_str << number;
1152 // std::string param = number_str.str();
1153 // system((script+param).c_str());
1157 if ( (totalNumberOfIterations> 30000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS2_arg) ) break;
1161 if (m_Verbose) std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1162 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1163 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1166 thresholder->SetInput( levelSetFilter->GetOutput() );
1167 thresholder->Update();
1168 filledRibs=thresholder->GetOutput();
1169 // writeImage<InternalImageType>(filledRibs,"/home/jef/tmp/filled_ribs_before.mhd", m_Verbose);
1173 //---------------------------------
1174 // Write the filled ribs
1175 //---------------------------------
1176 if (m_ArgsInfo.writeFilledRibs_given) {
1177 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1178 caster->SetInput(filledRibs);
1179 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeFilledRibs_arg, m_Verbose);
1183 //----------------------------------------------------------------------------------------------------
1184 // Collapse to the lungs
1185 //----------------------------------------------------------------------------------------------------
1187 std::cout<<std::endl;
1188 std::cout<<"=========================================="<<std::endl;
1189 std::cout<<"|| Collapsing to the lung image ||"<<std::endl;
1190 std::cout<<"=========================================="<<std::endl;
1193 //---------------------------------
1194 // Make feature image air+bones
1195 //---------------------------------
1196 if (m_Verbose) std::cout<<"Making feature images..."<<std::endl;
1197 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
1198 setBackgroundFilter->SetInput(air);
1199 setBackgroundFilter->SetInput2(bones);
1200 setBackgroundFilter->SetMaskValue(0);
1201 setBackgroundFilter->SetOutsideValue(0);
1202 setBackgroundFilter->Update();
1203 typename InternalImageType::Pointer bones_air =setBackgroundFilter->GetOutput();
1205 typename SetBackgroundFilterType::Pointer setBackgroundFilter2 = SetBackgroundFilterType::New();
1206 setBackgroundFilter2->SetInput(bones_air);
1207 setBackgroundFilter2->SetInput2(lungs);
1208 setBackgroundFilter2->SetMaskValue(0);
1209 setBackgroundFilter2->SetOutsideValue(0);
1210 setBackgroundFilter2->Update();
1211 typename InternalImageType::Pointer lungs_bones_air =setBackgroundFilter2->GetOutput();
1213 //---------------------------------
1214 // Prepare previous solution
1215 //---------------------------------
1216 if (m_Verbose) std::cout<<"Resampling previous solution..."<<std::endl;
1217 typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
1218 typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
1219 typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
1220 resampler->SetOutputStartIndex(bones->GetLargestPossibleRegion().GetIndex());
1221 resampler->SetInput(filledRibs);
1222 resampler->SetInterpolator(interpolator);
1223 resampler->SetOutputParametersFromImage(bones);
1224 resampler->Update();
1225 filledRibs =resampler->GetOutput();
1227 if (m_Verbose) std::cout<<"Setting lungs to 1..."<<std::endl;
1228 typename SetBackgroundFilterType::Pointer setBackgroundFilter3 = SetBackgroundFilterType::New();
1229 setBackgroundFilter3->SetInput(filledRibs);
1230 setBackgroundFilter3->SetInput2(lungs);
1231 setBackgroundFilter3->SetMaskValue(0);
1232 setBackgroundFilter3->SetOutsideValue(1);
1233 setBackgroundFilter3->Update();
1234 filledRibs=setBackgroundFilter3->GetOutput();
1236 if (m_Verbose) std::cout<<"Removing overlap with bones..."<<std::endl;
1237 typename SetBackgroundFilterType::Pointer setBackgroundFilter4 = SetBackgroundFilterType::New();
1238 setBackgroundFilter4->SetInput(filledRibs);
1239 setBackgroundFilter4->SetInput2(bones);
1240 setBackgroundFilter4->SetMaskValue(0);
1241 setBackgroundFilter4->SetOutsideValue(0);
1242 setBackgroundFilter4->Update();
1243 filledRibs =setBackgroundFilter4->GetOutput();
1244 // writeImage<InternalImageType>(filledRibs,"/home/jef/tmp/filledRibs_pp.mhd");
1245 //---------------------------------
1246 // Calculate Distance Map
1247 //---------------------------------
1248 // typedef itk::ApproximateSignedDistanceMapImageFilter <InternalImageType, LevelSetImageType> DistanceMapImageFilterType2;
1249 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
1250 distanceMapImageFilter->SetInput(filledRibs);
1251 distanceMapImageFilter->SetInsideIsPositive(false);
1252 // distanceMapImageFilter->SetInsideValue(0);
1253 // distanceMapImageFilter->SetOutsideValue(1);
1254 if (m_Verbose) std::cout<<"Calculating distance map..."<<std::endl;
1255 distanceMapImageFilter->Update();
1257 //---------------------------------
1259 //---------------------------------
1260 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
1261 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
1262 levelSetFilter->SetFeatureImage( lungs_bones_air );
1263 levelSetFilter->SetPropagationScaling(m_ArgsInfo.prop3_arg);
1264 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve3_arg );
1265 levelSetFilter->SetAdvectionScaling( 0 );
1266 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS3_arg );
1267 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter3_arg );
1268 levelSetFilter->SetUseImageSpacing(true);
1271 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1274 if (m_Verbose) std::cout<<"Starting the levelset..."<<std::endl;
1275 unsigned int totalNumberOfIterations=0;
1277 levelSetFilter->Update();
1280 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1281 levelSetFilter->SetInput(levelSetFilter->GetOutput());
1282 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg);
1283 std::cout << "After "<<totalNumberOfIterations<<" iteration the RMS change is " << levelSetFilter->GetRMSChange() <<"..."<< std::endl;
1286 // std::ostringstream number_str;
1287 // number_str << number;
1288 // std::string param = number_str.str();
1289 // system((script+param).c_str());
1292 // stopping condition
1293 if ( (totalNumberOfIterations>= (unsigned int) m_ArgsInfo.maxIter3_arg) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS3_arg) ) break;
1294 levelSetFilter->SetNumberOfIterations( std::min ((unsigned int) m_ArgsInfo.iter3_arg, (unsigned int) m_ArgsInfo.maxIter3_arg-totalNumberOfIterations ) );
1299 std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1300 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1301 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1305 typedef itk::BinaryThresholdImageFilter< LevelSetImageType,InternalImageType> LevelSetBinarizeFilterType;
1306 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1307 thresholder->SetUpperThreshold( 0.0 );
1308 thresholder->SetOutsideValue( 0 );
1309 thresholder->SetInsideValue( 1 );
1310 thresholder->SetInput( levelSetFilter->GetOutput() );
1311 thresholder->Update();
1312 typename InternalImageType::Pointer output = thresholder->GetOutput();
1315 //----------------------------------------------------------------------------------------------------
1316 // Prepare the output
1317 //----------------------------------------------------------------------------------------------------
1319 //---------------------------------
1321 //---------------------------------
1322 if (m_Verbose) std::cout<<"Removing overlap mask with air..."<<std::endl;
1323 typename SetBackgroundFilterType::Pointer setBackgroundFilter5 = SetBackgroundFilterType::New();
1324 setBackgroundFilter5->SetInput(output);
1325 setBackgroundFilter5->SetInput2(air);
1326 setBackgroundFilter5->SetMaskValue(0);
1327 setBackgroundFilter5->SetOutsideValue(0);
1328 setBackgroundFilter5->Update();
1329 output=setBackgroundFilter5->GetOutput();
1331 //---------------------------------
1332 // Open & close to cleanup
1333 //---------------------------------
1334 if ( m_ArgsInfo.openClose_flag) {
1336 //---------------------------------
1337 // Structuring element
1338 //---------------------------------
1339 typedef itk::BinaryBallStructuringElement<InternalPixelType,Dimension > KernelType;
1340 KernelType structuringElement;
1341 structuringElement.SetRadius(1);
1342 structuringElement.CreateStructuringElement();
1344 //---------------------------------
1346 //---------------------------------
1347 typedef itk::BinaryMorphologicalOpeningImageFilter<InternalImageType, InternalImageType , KernelType> OpenFilterType;
1348 typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
1349 openFilter->SetInput(output);
1350 openFilter->SetBackgroundValue(0);
1351 openFilter->SetForegroundValue(1);
1352 openFilter->SetKernel(structuringElement);
1353 if(m_Verbose) std::cout<<"Opening the output image..."<<std::endl;
1355 //---------------------------------
1357 //---------------------------------
1358 typedef itk::BinaryMorphologicalClosingImageFilter<InternalImageType, InternalImageType , KernelType> CloseFilterType;
1359 typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
1360 closeFilter->SetInput(openFilter->GetOutput());
1361 closeFilter->SetSafeBorder(true);
1362 closeFilter->SetForegroundValue(1);
1363 closeFilter->SetKernel(structuringElement);
1364 if(m_Verbose) std::cout<<"Closing the output image..."<<std::endl;
1365 closeFilter->Update();
1366 output=closeFilter->GetOutput();
1369 //writeImage<InternalImageType>(output,"/home/jef/tmp/mm_double.mhd");
1371 // Extract the upper part
1372 typedef itk::CropImageFilter<InternalImageType, InternalImageType> CropImageFilterType;
1373 typename CropImageFilterType::Pointer cropFilter=CropImageFilterType::New();
1374 cropFilter->SetInput(output);
1375 typename InputImageType::SizeType lSize, uSize;
1378 lSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
1379 cropFilter->SetLowerBoundaryCropSize(lSize);
1380 cropFilter->SetUpperBoundaryCropSize(uSize);
1381 cropFilter->Update();
1384 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1385 caster->SetInput(cropFilter->GetOutput());
1386 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.output_arg, m_Verbose);
1392 #endif //#define clitkMotionMaskGenericFilter_txx