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 ===================================================*/
30 #include "itkLabelImageToShapeLabelMapFilter.h"
31 #include "itkShapeLabelObject.h"
36 //-------------------------------------------------------------------
37 // Update with the number of dimensions
38 //-------------------------------------------------------------------
39 template<unsigned int Dimension>
41 MotionMaskGenericFilter::UpdateWithDim(std::string PixelType)
43 if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<"..."<<std::endl;
45 if(PixelType == "short") {
46 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
47 UpdateWithDimAndPixelType<Dimension, signed short>();
49 // else if(PixelType == "unsigned_short"){
50 // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
51 // UpdateWithDimAndPixelType<Dimension, unsigned short>();
54 // else if (PixelType == "unsigned_char"){
55 // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
56 // UpdateWithDimAndPixelType<Dimension, unsigned char>();
59 // else if (PixelType == "char"){
60 // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
61 // UpdateWithDimAndPixelType<Dimension, signed char>();
64 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
65 UpdateWithDimAndPixelType<Dimension, float>();
69 //-------------------------------------------------------------------
71 //-------------------------------------------------------------------
72 template <unsigned int Dimension, class PixelType>
73 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
74 MotionMaskGenericFilter::GetAirImage(typename itk::Image<PixelType, Dimension>::Pointer input,
75 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer lungs)
79 typedef int InternalPixelType;
80 typedef itk::Image<PixelType, Dimension> InputImageType;
81 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
83 //----------------------------------------------------------------------------------------------------
84 // Typedef for Processing
85 //----------------------------------------------------------------------------------------------------
86 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
87 typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
88 typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> InversionFilterType;
89 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdFilterType;
90 typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
91 typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
92 typedef itk::MirrorPadImageFilter<InternalImageType, InternalImageType> MirrorPadImageFilterType;
95 typename InternalImageType::Pointer air;
96 if (m_ArgsInfo.featureAir_given) {
97 typename FeatureReaderType::Pointer featureReader =FeatureReaderType::New();
98 featureReader->SetFileName(m_ArgsInfo.featureAir_arg);
99 if (m_Verbose) std::cout<<"Reading the air feature image..."<<std::endl;
100 featureReader->Update();
101 air=featureReader->GetOutput();
103 //---------------------------------
105 //---------------------------------
106 if(m_ArgsInfo.pad_flag) {
107 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
108 IteratorType it(air, air->GetLargestPossibleRegion());
109 typename InternalImageType::IndexType index;
110 while(!it.IsAtEnd()) {
113 //Pad borders of all dimensions but 2 (cranio-caudal)
114 for (unsigned int i=0; i<Dimension; i++){
117 if(index[i]==air->GetLargestPossibleRegion().GetIndex()[i]
118 || index[i]==(unsigned int)air->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air->GetLargestPossibleRegion().GetSize()[i]-1)
125 if (m_Verbose) std::cout<<"Extracting the air feature image..."<<std::endl;
126 //---------------------------------
127 // Binarize the image
128 //---------------------------------
129 typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
130 binarizeFilter->SetInput(input);
131 binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdAir_arg));
132 binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdAir_arg));
133 if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<<m_ArgsInfo.lowerThresholdAir_arg
134 <<", "<<m_ArgsInfo.upperThresholdAir_arg<<"..."<<std::endl;
135 binarizeFilter->Update();
136 air = binarizeFilter->GetOutput();
138 //---------------------------------
140 //---------------------------------
141 if(m_ArgsInfo.pad_flag) {
142 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
143 IteratorType it(air, air->GetLargestPossibleRegion());
144 typename InternalImageType::IndexType index;
145 while(!it.IsAtEnd()) {
148 //Pad borders of all dimensions but 2 (cranio-caudal)
149 for (unsigned int i=0; i<Dimension; i++){
152 if(index[i]==air->GetLargestPossibleRegion().GetIndex()[i]
153 || index[i]==(unsigned int)air->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air->GetLargestPossibleRegion().GetSize()[i]-1)
154 it.Set(binarizeFilter->GetInsideValue());
160 //---------------------------------
161 // Remove lungs (in place)
162 //---------------------------------
163 typedef itk::ImageRegionIterator<InternalImageType> IteratorType;
164 IteratorType itAir(air, binarizeFilter->GetOutput()->GetLargestPossibleRegion());
165 IteratorType itLungs(lungs, binarizeFilter->GetOutput()->GetLargestPossibleRegion()); //lungs padded, used air region
168 while(!itAir.IsAtEnd()) {
169 if(!itLungs.Get()) // The lungs are set to 0 at that stage
175 //---------------------------------
176 // Label the connected components
177 //---------------------------------
178 typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
179 connectFilter->SetInput(binarizeFilter->GetOutput());
180 connectFilter->SetBackgroundValue(0);
181 connectFilter->SetFullyConnected(false);
182 if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
184 //---------------------------------
185 // Sort the labels according to size
186 //---------------------------------
187 typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
188 relabelFilter->SetInput(connectFilter->GetOutput());
189 if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
191 //---------------------------------
193 //---------------------------------
194 typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
195 thresholdFilter->SetInput(relabelFilter->GetOutput());
196 thresholdFilter->SetUpper(1);
197 if (m_Verbose) std::cout<<"Keeping the first label..."<<std::endl;
198 thresholdFilter->Update();
199 air=thresholdFilter->GetOutput();
202 //---------------------------------
204 //---------------------------------
205 typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
206 inversionFilter->SetInput(air);
207 inversionFilter->SetLowerThreshold(0);
208 inversionFilter->SetUpperThreshold(0);
209 inversionFilter->SetInsideValue(1);
210 inversionFilter->Update();
212 //---------------------------------
214 //---------------------------------
215 typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
216 mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
217 typename InternalImageType::SizeType padSize;
219 padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
220 mirrorPadImageFilter->SetPadLowerBound(padSize);
221 if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
222 mirrorPadImageFilter->Update();
223 air=mirrorPadImageFilter->GetOutput();
224 //writeImage<InternalImageType>(air,"/home/srit/tmp/air.mhd");
230 //-------------------------------------------------------------------
232 //-------------------------------------------------------------------
233 template <unsigned int Dimension, class PixelType>
234 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
235 MotionMaskGenericFilter::GetBonesImage(typename itk::Image<PixelType, Dimension>::Pointer input )
239 typedef int InternalPixelType;
240 typedef itk::Image<PixelType, Dimension> InputImageType;
241 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
243 //----------------------------------------------------------------------------------------------------
244 // Typedef for Processing
245 //----------------------------------------------------------------------------------------------------
246 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
247 typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
248 typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> InversionFilterType;
249 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdFilterType;
250 typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
251 typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
252 typedef itk::MirrorPadImageFilter<InternalImageType, InternalImageType> MirrorPadImageFilterType;
255 typename InternalImageType::Pointer bones;
256 if (m_ArgsInfo.featureBones_given) {
257 typename FeatureReaderType::Pointer featureReader =FeatureReaderType::New();
258 featureReader->SetFileName(m_ArgsInfo.featureBones_arg);
259 if (m_Verbose) std::cout<<"Reading the bones feature image..."<<std::endl;
260 featureReader->Update();
261 bones=featureReader->GetOutput();
263 if (m_Verbose) std::cout<<"Extracting the bones feature image..."<<std::endl;
264 //---------------------------------
265 // Binarize the image
266 //---------------------------------
267 typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
268 binarizeFilter->SetInput(input);
269 binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdBones_arg));
270 binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdBones_arg));
271 if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<<m_ArgsInfo.lowerThresholdBones_arg
272 <<", "<<m_ArgsInfo.upperThresholdBones_arg<<"..."<<std::endl;
274 //---------------------------------
275 // Label the connected components
276 //---------------------------------
277 typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
278 connectFilter->SetInput(binarizeFilter->GetOutput());
279 connectFilter->SetBackgroundValue(0);
280 connectFilter->SetFullyConnected(false);
281 if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
283 //---------------------------------
284 // Sort the labels according to size
285 //---------------------------------
286 typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
287 relabelFilter->SetInput(connectFilter->GetOutput());
288 if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
290 //---------------------------------
292 //---------------------------------
293 typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
294 thresholdFilter->SetInput(relabelFilter->GetOutput());
295 thresholdFilter->SetUpper(1);
296 if (m_Verbose) std::cout<<"Keeping the first label..."<<std::endl;
297 thresholdFilter->Update();
298 bones=thresholdFilter->GetOutput();
302 //---------------------------------
304 //---------------------------------
305 typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
306 inversionFilter->SetInput(bones);
307 inversionFilter->SetLowerThreshold(0);
308 inversionFilter->SetUpperThreshold(0);
309 inversionFilter->SetInsideValue(1);
310 inversionFilter->Update();
312 //---------------------------------
314 //---------------------------------
315 typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
316 mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
317 typename InternalImageType::SizeType padSize;
319 padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
320 mirrorPadImageFilter->SetPadLowerBound(padSize);
321 if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
322 mirrorPadImageFilter->Update();
323 bones=mirrorPadImageFilter->GetOutput();
324 // writeImage<InternalImageType>(bones,"/home/jef/tmp/bones.mhd");
332 //-------------------------------------------------------------------
334 //-------------------------------------------------------------------
335 template <unsigned int Dimension, class PixelType>
336 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
337 MotionMaskGenericFilter::GetLungsImage(typename itk::Image<PixelType, Dimension>::Pointer input )
340 typedef int InternalPixelType;
341 typedef itk::Image<PixelType, Dimension> InputImageType;
342 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
344 //----------------------------------------------------------------------------------------------------
345 // Typedef for Processing
346 //----------------------------------------------------------------------------------------------------
347 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
348 typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
349 typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> InversionFilterType;
350 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdFilterType;
351 typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
352 typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
353 typedef itk::MirrorPadImageFilter<InternalImageType, InternalImageType> MirrorPadImageFilterType;
355 typename InternalImageType::Pointer lungs;
356 if (m_ArgsInfo.featureLungs_given) {
357 typename FeatureReaderType::Pointer featureReader =FeatureReaderType::New();
358 featureReader->SetFileName(m_ArgsInfo.featureLungs_arg);
359 if (m_Verbose) std::cout<<"Reading the lungs feature image..."<<std::endl;
360 featureReader->Update();
361 lungs=featureReader->GetOutput();
363 if (m_Verbose) std::cout<<"Extracting the lungs feature image..."<<std::endl;
364 //---------------------------------
365 // Binarize the image
366 //---------------------------------
367 typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
368 binarizeFilter->SetInput(input);
369 binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdLungs_arg));
370 binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdLungs_arg));
371 if (m_Verbose) std::cout<<"Binarizing the image using a threshold "<<m_ArgsInfo.lowerThresholdLungs_arg
372 <<", "<<m_ArgsInfo.upperThresholdLungs_arg<<"..."<<std::endl;
374 //---------------------------------
375 // Label the connected components
376 //---------------------------------
377 typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
378 connectFilter->SetInput(binarizeFilter->GetOutput());
379 connectFilter->SetBackgroundValue(0);
380 connectFilter->SetFullyConnected(true);
381 if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
382 connectFilter->Update();
383 if (m_Verbose) std::cout<<"found "<< connectFilter->GetObjectCount() << std::endl;
385 //---------------------------------
386 // Sort the labels according to size
387 //---------------------------------
388 typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
389 relabelFilter->SetInput(connectFilter->GetOutput());
390 if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
391 // writeImage<InternalImageType> (relabelFilter->GetOutput(), "/home/vdelmon/tmp/labels.mhd");
393 //---------------------------------
395 //---------------------------------
396 typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
397 thresholdFilter->SetInput(relabelFilter->GetOutput());
398 thresholdFilter->SetLower(1);
399 thresholdFilter->SetUpper(2);
400 if (m_Verbose) std::cout<<"Keeping the first two labels..."<<std::endl;
401 thresholdFilter->Update();
402 lungs=thresholdFilter->GetOutput();
407 //---------------------------------
409 //---------------------------------
410 typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
411 inversionFilter->SetInput(lungs);
412 inversionFilter->SetLowerThreshold(0);
413 inversionFilter->SetUpperThreshold(0);
414 inversionFilter->SetInsideValue(1);
415 inversionFilter->Update();
417 //---------------------------------
419 //---------------------------------
420 typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
421 mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
422 typename InternalImageType::SizeType padSize;
424 padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
425 mirrorPadImageFilter->SetPadLowerBound(padSize);
426 if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
427 mirrorPadImageFilter->Update();
428 lungs=mirrorPadImageFilter->GetOutput();
429 // writeImage<InternalImageType>(lungs,"/home/jef/tmp/lungs.mhd");
434 //-------------------------------------------------------------------
436 //-------------------------------------------------------------------
437 template <unsigned int Dimension, class PixelType>
438 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
439 MotionMaskGenericFilter::Resample( typename itk::Image<InternalPixelType,Dimension>::Pointer input )
442 typedef int InternalPixelType;
443 typedef itk::Image<PixelType, Dimension> InputImageType;
444 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
446 //---------------------------------
447 // Resample to new spacing
448 //---------------------------------
449 typename InternalImageType::SizeType size_low;
450 typename InternalImageType::SpacingType spacing_low;
451 for (unsigned int i=0; i<Dimension; i++)
452 if (m_ArgsInfo.spacing_given==Dimension)
453 for (unsigned int i=0; i<Dimension; i++) {
454 spacing_low[i]=m_ArgsInfo.spacing_arg[i];
455 size_low[i]=ceil(input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i]/spacing_low[i]);
458 for (unsigned int i=0; i<Dimension; i++) {
459 spacing_low[i]=m_ArgsInfo.spacing_arg[0];
460 size_low[i]=ceil(input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i]/spacing_low[i]);
463 typename InternalImageType::PointType origin;
464 input->TransformIndexToPhysicalPoint(input->GetLargestPossibleRegion().GetIndex(), origin);
465 typedef itk::ResampleImageFilter<InternalImageType, InternalImageType> ResampleImageFilterType;
466 typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
467 typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
468 typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
469 resampler->SetInterpolator(interpolator);
470 resampler->SetOutputOrigin(origin);
471 resampler->SetOutputSpacing(spacing_low);
472 resampler->SetSize(size_low);
473 resampler->SetInput(input);
475 typename InternalImageType::Pointer output =resampler->GetOutput();
480 //-------------------------------------------------------------------
482 //-------------------------------------------------------------------
483 template <unsigned int Dimension, class PixelType>
484 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
485 MotionMaskGenericFilter::InitializeEllips( typename itk::Vector<double,Dimension> center, typename itk::Image<InternalPixelType,Dimension>::Pointer bones_low, typename itk::Image<InternalPixelType,Dimension>::Pointer lungs_low )
489 typedef int InternalPixelType;
490 typedef itk::Image<PixelType, Dimension> InputImageType;
491 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
492 typedef itk::Image< unsigned char , Dimension> OutputImageType;
493 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
494 typedef clitk::SetBackgroundImageFilter<InternalImageType,InternalImageType, InternalImageType> SetBackgroundFilterType;
495 typedef itk::LabelStatisticsImageFilter<InternalImageType,InternalImageType> LabelStatisticsImageFilterType;
496 typedef itk::CastImageFilter<InternalImageType,OutputImageType> CastImageFilterType;
499 typename InternalImageType::Pointer ellips;
500 if (m_ArgsInfo.ellips_given || m_ArgsInfo.grownEllips_given || m_ArgsInfo.filledRibs_given) {
501 if(m_ArgsInfo.ellips_given) {
502 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
503 typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
504 featureReader->SetFileName(m_ArgsInfo.ellips_arg);
505 featureReader->Update();
506 ellips=featureReader->GetOutput();
510 std::cout<<std::endl;
511 std::cout<<"=========================================="<<std::endl;
512 std::cout<<"|| Initializing ellipsoide ||"<<std::endl;
513 std::cout<<"=========================================="<<std::endl;
516 itk::Vector<double,Dimension> centerEllips;
517 typename itk::Vector<double, Dimension> axes;
518 if (m_ArgsInfo.ellipseAutoDetect_flag) {
520 std::cout<<"Auto-detecting intial ellipse..."<<std::endl;
523 typename InternalImageType::RegionType region, region_before = lungs_low->GetLargestPossibleRegion();
524 typename InternalImageType::SizeType size = region_before.GetSize();
526 region.SetSize(size);
527 lungs_low->SetRegions(region);
529 // compute statistics of the (mirroed) lung region in order to guess the params of the ellipse
530 typedef itk::LabelImageToShapeLabelMapFilter<InternalImageType> LabelImageToShapeLabelMapFilter;
531 typename LabelImageToShapeLabelMapFilter::Pointer label_filter = LabelImageToShapeLabelMapFilter::New();
532 label_filter->SetInput(lungs_low);
533 label_filter->Update();
535 typename InternalImageType::SpacingType spacing = lungs_low->GetSpacing();
536 typedef typename LabelImageToShapeLabelMapFilter::OutputImageType::LabelObjectType LabelType;
537 const LabelType* label = dynamic_cast<LabelType*>(label_filter->GetOutput()->GetNthLabelObject(0));
540 typename LabelType::CentroidType lung_centroid = label->GetCentroid();
542 // try to guess ideal ellipse axes from the lungs' bounding box and centroid
543 // with some hard-coded "magic" constants...
544 #if ITK_VERSION_MAJOR >= 4
545 typename LabelType::RegionType lung_bbox = label->GetBoundingBox();
547 typename LabelType::RegionType lung_bbox = label->GetRegion();
549 axes[0] = 0.95*lung_bbox.GetSize()[0]*spacing[0]/2;
550 axes[1] = 0.3*lung_bbox.GetSize()[1]*spacing[1]/2;
551 axes[2] = 1.25*fabs(lung_centroid[2] - center[2]);
553 // ellipse's center in XY is the center of the lungs' bounding box
554 typename InternalImageType::PointType origin = lungs_low->GetOrigin();
555 centerEllips[0] = origin[0] + (lung_bbox.GetIndex()[0] + lung_bbox.GetSize()[0]/2)*spacing[0];
556 centerEllips[1] = origin[1] + (lung_bbox.GetIndex()[1] + lung_bbox.GetSize()[1]/2)*spacing[1];
557 centerEllips[2] = center[2];
559 lungs_low->SetRegions(region_before);
562 std::cout << "Lungs'centroid at " << lung_centroid << std::endl;
563 std::cout << "Ellipse center at " << centerEllips << std::endl;
564 std::cout << "Ellipse axes as " << axes << std::endl;
565 std::cout << "Lung's bounding box (index,size) " << lung_bbox.GetIndex() << lung_bbox.GetSize() << std::endl;
569 //---------------------------------
570 // Offset from center
571 //---------------------------------
572 typename itk::Vector<double,Dimension> offset;
573 if(m_ArgsInfo.offset_given== Dimension) {
574 for(unsigned int i=0; i<Dimension; i++)
575 offset[i]=m_ArgsInfo.offset_arg[i];
580 centerEllips=center+offset;
582 std::cout<<"Placing the center of the initial ellipsoide at (mm) "<<std::endl;
583 std::cout<<centerEllips[0];
584 for (unsigned int i=1; i<Dimension; i++)
585 std::cout<<" "<<centerEllips[i];
586 std::cout<<std::endl;
590 if (m_ArgsInfo.axes_given==Dimension)
591 for(unsigned int i=0; i<Dimension; i++)
592 axes[i]=m_ArgsInfo.axes_arg[i];
600 //---------------------------------
602 //---------------------------------
603 ellips=InternalImageType::New();
604 ellips->SetRegions (bones_low->GetLargestPossibleRegion());
605 ellips->SetOrigin(bones_low->GetOrigin());
606 ellips->SetSpacing(bones_low->GetSpacing());
608 ellips->FillBuffer(0);
611 IteratorType itEllips(ellips, ellips->GetLargestPossibleRegion());
612 itEllips.GoToBegin();
613 typename InternalImageType::PointType point;
614 typename InternalImageType::IndexType index;
617 if (m_Verbose) std::cout<<"Drawing the initial ellipsoide..."<<std::endl;
618 while (!itEllips.IsAtEnd()) {
619 index=itEllips.GetIndex();
620 ellips->TransformIndexToPhysicalPoint(index, point);
622 for(unsigned int i=0; i<Dimension; i++)
623 distance+=powf( ( (centerEllips[i]-point[i])/axes[i] ), 2);
632 //---------------------------------
634 //---------------------------------
635 if (m_ArgsInfo.writeEllips_given) {
636 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
637 caster->SetInput(ellips);
638 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeEllips_arg, m_Verbose);
646 //-------------------------------------------------------------------
647 // Update with the number of dimensions and the pixeltype
648 //-------------------------------------------------------------------
649 template <unsigned int Dimension, class PixelType>
651 MotionMaskGenericFilter::UpdateWithDimAndPixelType()
655 typedef int InternalPixelType;
656 typedef itk::Image<PixelType, Dimension> InputImageType;
657 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
658 typedef itk::Image< float, Dimension> LevelSetImageType; //labeling doesn't work with unsigned char?
659 typedef itk::Image<unsigned char, Dimension> OutputImageType;
662 //----------------------------------------------------------------------------------------------------
663 // Typedef for Processing
664 //----------------------------------------------------------------------------------------------------
665 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
666 typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
667 typedef itk::BinaryThresholdImageFilter< LevelSetImageType,InternalImageType> LevelSetBinarizeFilterType;
668 typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> InversionFilterType;
669 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdFilterType;
670 typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
671 typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
672 typedef itk::MirrorPadImageFilter<InternalImageType, InternalImageType> MirrorPadImageFilterType;
673 typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
674 typedef itk::ResampleImageFilter<InternalImageType, InternalImageType> ResampleImageFilterType;
675 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
676 typedef itk::GeodesicActiveContourLevelSetImageFilter<LevelSetImageType, InternalImageType> LevelSetImageFilterType;
677 typedef itk::SignedDanielssonDistanceMapImageFilter<InternalImageType, LevelSetImageType> DistanceMapImageFilterType;
678 typedef clitk::SetBackgroundImageFilter<InternalImageType,InternalImageType, InternalImageType> SetBackgroundFilterType;
679 typedef itk::LabelStatisticsImageFilter<InternalImageType,InternalImageType> LabelStatisticsImageFilterType;
680 typedef itk::CastImageFilter<InternalImageType,OutputImageType> CastImageFilterType;
684 typedef itk::ImageFileReader<InputImageType> InputReaderType;
685 typename InputReaderType::Pointer reader = InputReaderType::New();
686 reader->SetFileName( m_InputFileName);
688 typename InputImageType::Pointer input= reader->GetOutput();
690 // // globals for avi
691 // unsigned int number=0;
695 std::cout<<std::endl;
696 std::cout<<"=========================================="<<std::endl;
697 std::cout<<"|| Making feature images ||"<<std::endl;
698 std::cout<<"=========================================="<<std::endl;
701 //--------------------------------------------------------------------------------
703 //-------------------------------------------------------------------------------
704 typename InternalImageType::Pointer lungs = this->GetLungsImage<Dimension, PixelType>(input);
706 //-------------------------------------------------------------------------------
708 //-------------------------------------------------------------------------------
709 typename InternalImageType::Pointer air = this->GetAirImage<Dimension, PixelType>(input, lungs);
711 //-------------------------------------------------------------------------------
713 //-------------------------------------------------------------------------------
714 typename InternalImageType::Pointer bones = this->GetBonesImage<Dimension, PixelType>(input);
716 //----------------------------------------------------------------------------------------------------
717 // Write feature images
718 //----------------------------------------------------------------------------------------------------
719 if(m_ArgsInfo.writeFeature_given) {
720 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
721 setBackgroundFilter->SetInput(air);
722 setBackgroundFilter->SetInput2(bones);
723 setBackgroundFilter->SetMaskValue(0);
724 setBackgroundFilter->SetOutsideValue(2);
725 setBackgroundFilter->Update();
726 typename InternalImageType::Pointer bones_air =setBackgroundFilter->GetOutput();
728 typename SetBackgroundFilterType::Pointer setBackgroundFilter2 = SetBackgroundFilterType::New();
729 setBackgroundFilter2->SetInput(bones_air);
730 setBackgroundFilter2->SetInput2(lungs);
731 setBackgroundFilter2->SetMaskValue(0);
732 setBackgroundFilter2->SetOutsideValue(3);
733 setBackgroundFilter2->Update();
734 typename InternalImageType::Pointer lungs_bones_air =setBackgroundFilter2->GetOutput();
736 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
737 caster->SetInput(lungs_bones_air);
738 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeFeature_arg, m_Verbose);
741 //----------------------------------------------------------------------------------------------------
742 // Low dimensional versions
743 //----------------------------------------------------------------------------------------------------
744 typename InternalImageType::Pointer bones_low =Resample<Dimension , PixelType>(bones);
745 typename InternalImageType::Pointer lungs_low =Resample<Dimension , PixelType>(lungs);
746 typename InternalImageType::Pointer air_low =Resample<Dimension , PixelType>(air);
748 //---------------------------------
750 //---------------------------------
751 if(m_ArgsInfo.pad_flag) {
752 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
753 IteratorType it(air_low, air_low->GetLargestPossibleRegion());
754 typename InternalImageType::IndexType index;
755 while(!it.IsAtEnd()) {
757 for (unsigned int i=0; i<Dimension; i++)
758 if(index[i]==air_low->GetLargestPossibleRegion().GetIndex()[i]
759 || index[i]==(unsigned int)air_low->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air_low->GetLargestPossibleRegion().GetSize()[i]-1)
765 //---------------------------------
767 //---------------------------------
768 typename itk::Vector<double,Dimension> center;
769 typedef itk::ImageMomentsCalculator<InternalImageType> MomentsCalculatorType;
770 typename MomentsCalculatorType::Pointer momentsCalculator=MomentsCalculatorType::New();
771 momentsCalculator->SetImage(air);
772 momentsCalculator->Compute();
773 center=momentsCalculator->GetCenterOfGravity();
775 std::cout<<"The center of gravity of the patient body is located at (mm) "<<std::endl;
776 std::cout<<center[0];
777 for (unsigned int i=1; i<Dimension; i++)
778 std::cout<<" "<<center[i];
779 std::cout<<std::endl;
782 //----------------------------------------------------------------------------------------------------
783 // Ellipsoide initialization
784 //----------------------------------------------------------------------------------------------------
785 typename InternalImageType::Pointer m_Ellips= InitializeEllips<Dimension, PixelType>(center,bones_low,lungs_low);
787 //----------------------------------------------------------------------------------------------------
789 //----------------------------------------------------------------------------------------------------
790 typename InternalImageType::Pointer grownEllips;
791 if (m_ArgsInfo.grownEllips_given || m_ArgsInfo.filledRibs_given) {
792 if (m_ArgsInfo.grownEllips_given) {
794 typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
795 featureReader->SetFileName(m_ArgsInfo.grownEllips_arg);
796 featureReader->Update();
797 grownEllips=featureReader->GetOutput();
802 std::cout<<std::endl;
803 std::cout<<"=========================================="<<std::endl;
804 std::cout<<"|| Growing ellipsoide ||"<<std::endl;
805 std::cout<<"=========================================="<<std::endl;
808 //---------------------------------
810 //---------------------------------
811 typename InternalImageType::PointType dPoint;
812 if (m_ArgsInfo.detectionPoint_given) {
813 for (unsigned int i=0; i<Dimension; i++)
814 dPoint[i]=m_ArgsInfo.detectionPoint_arg[i];
816 typename InternalImageType::RegionType searchRegion=air->GetLargestPossibleRegion();
817 typename InternalImageType::SizeType searchRegionSize = searchRegion.GetSize();
818 typename InternalImageType::IndexType searchRegionIndex = searchRegion.GetIndex();
819 searchRegionIndex[2]+=searchRegionSize[2]/2;
820 searchRegionSize[2]=1;
821 searchRegion.SetSize(searchRegionSize);
822 searchRegion.SetIndex(searchRegionIndex);
823 IteratorType itAbdomen(air, searchRegion);
824 itAbdomen.GoToBegin();
826 typename InternalImageType::PointType aPoint;
827 typename InternalImageType::IndexType aIndex;
829 if (m_Verbose) std::cout<<"Detecting the abdomen..."<<std::endl;
830 while (!itAbdomen.IsAtEnd()) {
831 if(itAbdomen.Get()) break;
834 aIndex=itAbdomen.GetIndex();
835 air->TransformIndexToPhysicalPoint(aIndex,aPoint);
836 if (m_Verbose) std::cout<<"Detected the abdomen at "<<aPoint<<".."<<std::endl;
839 //---------------------------------
840 // Detect abdomen in additional images?
841 //---------------------------------
842 if (m_ArgsInfo.detectionPairs_given) {
843 for (unsigned int i=0; i<m_ArgsInfo.detectionPairs_given; i++) {
844 typename InternalImageType::Pointer airAdd;
845 //---------------------------------
847 //--------------------------------
848 typedef itk::ImageFileReader<InputImageType> InputReaderType;
849 typename InputReaderType::Pointer reader = InputReaderType::New();
850 reader->SetFileName( m_ArgsInfo.detectionPairs_arg[i]);
852 typename InputImageType::Pointer additional= reader->GetOutput();
853 if (m_Verbose) std::cout<<"Extracting the air from additional image "<< i<<"..."<<std::endl;
855 //---------------------------------
856 // Binarize the image
857 //---------------------------------
858 typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
859 binarizeFilter->SetInput(additional);
860 binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdAir_arg));
861 binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdAir_arg));
862 if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<<m_ArgsInfo.lowerThresholdAir_arg
863 <<", "<<m_ArgsInfo.upperThresholdAir_arg<<"..."<<std::endl;
865 //---------------------------------
866 // Label the connected components
867 //---------------------------------
868 typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
869 connectFilter->SetInput(binarizeFilter->GetOutput());
870 connectFilter->SetBackgroundValue(0);
871 connectFilter->SetFullyConnected(false);
872 if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
874 //---------------------------------
875 // Sort the labels according to size
876 //---------------------------------
877 typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
878 relabelFilter->SetInput(connectFilter->GetOutput());
879 if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
881 //---------------------------------
883 //---------------------------------
884 typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
885 thresholdFilter->SetInput(relabelFilter->GetOutput());
886 thresholdFilter->SetUpper(1);
887 if (m_Verbose) std::cout<<"Keeping the first label..."<<std::endl;
888 thresholdFilter->Update();
889 airAdd=thresholdFilter->GetOutput();
892 //---------------------------------
894 //---------------------------------
895 typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
896 inversionFilter->SetInput(airAdd);
897 inversionFilter->SetLowerThreshold(0);
898 inversionFilter->SetUpperThreshold(0);
899 inversionFilter->SetInsideValue(1);
900 inversionFilter->Update();
902 //---------------------------------
904 //---------------------------------
905 typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
906 mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
907 typename InternalImageType::SizeType padSize;
909 padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
910 mirrorPadImageFilter->SetPadLowerBound(padSize);
911 if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
912 mirrorPadImageFilter->Update();
913 airAdd=mirrorPadImageFilter->GetOutput();
915 //---------------------------------
917 //---------------------------------
918 typename InternalImageType::RegionType searchRegion=airAdd->GetLargestPossibleRegion();
919 typename InternalImageType::SizeType searchRegionSize = searchRegion.GetSize();
920 typename InternalImageType::IndexType searchRegionIndex = searchRegion.GetIndex();
921 searchRegionIndex[2]+=searchRegionSize[2]/2;
922 searchRegionSize[2]=1;
923 searchRegion.SetSize(searchRegionSize);
924 searchRegion.SetIndex(searchRegionIndex);
925 IteratorType itAbdomen(airAdd, searchRegion);
926 itAbdomen.GoToBegin();
928 typename InternalImageType::PointType additionalPoint;
929 typename InternalImageType::IndexType aIndex;
931 if (m_Verbose) std::cout<<"Detecting the abdomen..."<<std::endl;
932 while (!itAbdomen.IsAtEnd()) {
933 if(itAbdomen.Get()) break;
936 aIndex=itAbdomen.GetIndex();
937 airAdd->TransformIndexToPhysicalPoint(aIndex,additionalPoint);
938 if (m_Verbose) std::cout<<"Detected the abdomen in the additional image at "<<additionalPoint<<".."<<std::endl;
940 if(additionalPoint[1]< aPoint[1]) {
941 aPoint=additionalPoint;
942 if (m_Verbose) std::cout<<"Modifying the detected abdomen to "<<aPoint<<".."<<std::endl;
949 // Determine the detection point
953 if(m_ArgsInfo.offsetDetect_given==Dimension)
954 for(unsigned int i=0; i <Dimension; i++)
955 dPoint[i]+=m_ArgsInfo.offsetDetect_arg[i];
960 if (m_Verbose) std::cout<<"Setting the detection point to "<<dPoint<<".."<<std::endl;
963 //---------------------------------
964 // Pad the rib image and ellips image
965 //---------------------------------
966 typename InternalImageType::Pointer padded_ellips;
967 typename InternalImageType::Pointer padded_bones_low;
969 // If detection point not inside the image: pad
970 typename InternalImageType::IndexType dIndex;
971 if (!bones_low->TransformPhysicalPointToIndex(dPoint, dIndex)) {
972 typename InternalImageType::SizeType padSize;
974 padSize[1]=abs(dIndex[1])+1;
975 if (m_Verbose) std::cout<<"Padding the images with padding size "<<padSize<<" to include the detection point..."<<std::endl;
977 typename MirrorPadImageFilterType::Pointer padBonesFilter=MirrorPadImageFilterType::New();
978 padBonesFilter->SetInput(bones_low);
979 padBonesFilter->SetPadLowerBound(padSize);
980 padBonesFilter->Update();
981 padded_bones_low=padBonesFilter->GetOutput();
983 typename MirrorPadImageFilterType::Pointer padEllipsFilter=MirrorPadImageFilterType::New();
984 padEllipsFilter->SetInput(m_Ellips);
985 padEllipsFilter->SetPadLowerBound(padSize);
986 padEllipsFilter->Update();
987 padded_ellips=padEllipsFilter->GetOutput();
991 padded_bones_low=bones_low;
992 padded_ellips=m_Ellips;
996 //---------------------------------
997 // Calculate distance map
998 //---------------------------------
999 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
1000 distanceMapImageFilter->SetInput(padded_ellips);
1001 distanceMapImageFilter->SetInsideIsPositive(false);
1002 if (m_Verbose) std::cout<<"Calculating the distance map..."<<std::endl;
1003 distanceMapImageFilter->Update();
1004 if (m_ArgsInfo.writeDistMap_given) {
1005 writeImage<LevelSetImageType>(distanceMapImageFilter->GetOutput(), m_ArgsInfo.writeDistMap_arg, m_Verbose);
1009 //---------------------------------
1010 // Grow while monitoring detection point
1011 //---------------------------------
1012 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
1013 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
1014 levelSetFilter->SetFeatureImage( padded_bones_low );
1015 levelSetFilter->SetPropagationScaling( 1 );
1016 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve1_arg );
1017 levelSetFilter->SetAdvectionScaling( 0 );
1018 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS1_arg );
1019 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter1_arg );
1020 levelSetFilter->SetUseImageSpacing(true);
1022 // //---------------------------------
1023 // // Script for making movie
1024 // //---------------------------------
1025 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1028 if (m_Verbose) std::cout<<"Starting level set segmentation..."<<std::endl;
1029 unsigned int totalNumberOfIterations=0;
1031 levelSetFilter->Update();
1032 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1034 if ( levelSetFilter->GetOutput()->GetPixel(dIndex) < 0 ) {
1035 if (m_Verbose) std::cout<<"Detection point reached!"<<std::endl;
1038 if (m_Verbose) std::cout<<"Detection point not reached after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1039 levelSetFilter->SetInput(levelSetFilter->GetOutput());
1040 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose);
1043 // std::ostringstream number_str;
1044 // number_str << number;
1045 // std::string param = number_str.str();
1046 // system((script+param).c_str());
1050 if ( (totalNumberOfIterations> 10000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS1_arg) ) break;
1054 if (m_Verbose) std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1055 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1056 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1059 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1060 thresholder->SetUpperThreshold( 0.0 );
1061 thresholder->SetOutsideValue( 0 );
1062 thresholder->SetInsideValue( 1 );
1063 thresholder->SetInput( levelSetFilter->GetOutput() );
1064 if (m_Verbose) std::cout<<"Thresholding the output level set..."<<std::endl;
1065 thresholder->Update();
1066 grownEllips=thresholder->GetOutput();
1069 //---------------------------------
1070 // Write the grown ellips
1071 //---------------------------------
1072 if (m_ArgsInfo.writeGrownEllips_given) {
1073 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1074 caster->SetInput(grownEllips);
1075 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeGrownEllips_arg, m_Verbose);
1079 //----------------------------------------------------------------------------------------------------
1081 //----------------------------------------------------------------------------------------------------
1082 typename InternalImageType::Pointer filledRibs;
1083 if (m_ArgsInfo.filledRibs_given) {
1084 typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
1085 featureReader->SetFileName(m_ArgsInfo.filledRibs_arg);
1086 if (m_Verbose) std::cout<<"Reading filled ribs mask..."<<std::endl;
1087 featureReader->Update();
1088 filledRibs=featureReader->GetOutput();
1091 std::cout<<std::endl;
1092 std::cout<<"=========================================="<<std::endl;
1093 std::cout<<"|| Filling the ribs image ||"<<std::endl;
1094 std::cout<<"=========================================="<<std::endl;
1096 //---------------------------------
1097 // Make feature image air+bones
1098 //---------------------------------
1099 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
1100 setBackgroundFilter->SetInput(air_low);
1101 setBackgroundFilter->SetInput2(bones_low);
1102 setBackgroundFilter->SetMaskValue(0);
1103 setBackgroundFilter->SetOutsideValue(0);
1104 setBackgroundFilter->Update();
1105 typename InternalImageType::Pointer bones_air_low =setBackgroundFilter->GetOutput();
1107 //---------------------------------
1108 // Resampling previous solution
1109 //---------------------------------
1110 if (m_Verbose) std::cout<<"Resampling previous solution..."<<std::endl;
1111 typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
1112 typename InternalImageType::PointType origin;
1113 bones_low->TransformIndexToPhysicalPoint(bones_low->GetLargestPossibleRegion().GetIndex(), origin);
1114 resampler->SetOutputOrigin(origin);
1115 resampler->SetOutputSpacing(bones_low->GetSpacing());
1116 typename InternalImageType::SizeType size_low= bones_low->GetLargestPossibleRegion().GetSize();
1117 resampler->SetSize(size_low);
1118 typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
1119 resampler->SetInterpolator(interpolator);
1120 resampler->SetInput(grownEllips);
1121 resampler->Update();
1122 typename InternalImageType::Pointer grownEllips =resampler->GetOutput();
1125 //---------------------------------
1126 // Calculate Distance Map
1127 //---------------------------------
1128 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
1129 distanceMapImageFilter->SetInput(grownEllips);
1130 distanceMapImageFilter->SetInsideIsPositive(false);
1131 if (m_Verbose) std::cout<<"Calculating distance map..."<<std::endl;
1132 distanceMapImageFilter->Update();
1133 if (m_ArgsInfo.writeDistMap_given) {
1134 writeImage<LevelSetImageType>(distanceMapImageFilter->GetOutput(), m_ArgsInfo.writeDistMap_arg, m_Verbose);
1138 //---------------------------------
1139 // Grow while monitoring lung volume coverage
1140 //---------------------------------
1141 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
1142 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
1143 levelSetFilter->SetFeatureImage( bones_air_low );
1144 levelSetFilter->SetPropagationScaling( 1 );
1145 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve2_arg );
1146 levelSetFilter->SetAdvectionScaling( 0 );
1147 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS2_arg );
1148 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter2_arg );
1149 levelSetFilter->SetUseImageSpacing(true);
1151 //---------------------------------
1153 //---------------------------------
1154 typename InversionFilterType::Pointer invertor= InversionFilterType::New();
1155 invertor->SetInput(lungs_low);
1156 invertor->SetLowerThreshold(0);
1157 invertor->SetUpperThreshold(0);
1158 invertor->SetInsideValue(1);
1160 typename InternalImageType::Pointer lungs_low_inv=invertor->GetOutput();
1162 //---------------------------------
1163 // Calculate lung volume
1164 //---------------------------------
1165 typename LabelStatisticsImageFilterType::Pointer statisticsFilter= LabelStatisticsImageFilterType::New();
1166 statisticsFilter->SetInput(lungs_low_inv);
1167 statisticsFilter->SetLabelInput(lungs_low);
1168 statisticsFilter->Update();
1169 unsigned int volume=statisticsFilter->GetSum(0);
1171 // Prepare threshold
1172 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1173 thresholder->SetUpperThreshold( 0.0 );
1174 thresholder->SetOutsideValue( 0 );
1175 thresholder->SetInsideValue( 1 );
1179 unsigned int totalNumberOfIterations=0;
1180 unsigned int coverage=0;
1183 // //---------------------------------
1184 // // Script for making movie
1185 // //---------------------------------
1186 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1189 levelSetFilter->Update();
1190 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1192 thresholder->SetInput( levelSetFilter->GetOutput() );
1193 thresholder->Update();
1194 statisticsFilter->SetInput(thresholder->GetOutput());
1195 statisticsFilter->Update();
1196 coverage=statisticsFilter->GetSum(0);
1198 // Compare the volumes
1199 if ( coverage >= m_ArgsInfo.fillingLevel_arg * volume/100 ) {
1200 if (m_Verbose) std::cout<<"Lungs filled for "<< m_ArgsInfo.fillingLevel_arg<<"% !"<<std::endl;
1203 if (m_Verbose) std::cout<<"After "<<totalNumberOfIterations<<" iterations, lungs are covered for "
1204 <<(double)coverage/volume*100.0<<"%..."<<std::endl;
1205 levelSetFilter->SetInput(levelSetFilter->GetOutput());
1206 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose);
1209 // std::ostringstream number_str;
1210 // number_str << number;
1211 // std::string param = number_str.str();
1212 // system((script+param).c_str());
1216 if ( (totalNumberOfIterations> 30000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS2_arg) ) break;
1220 if (m_Verbose) std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1221 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1222 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1225 thresholder->SetInput( levelSetFilter->GetOutput() );
1226 thresholder->Update();
1227 filledRibs=thresholder->GetOutput();
1228 // writeImage<InternalImageType>(filledRibs,"/home/jef/tmp/filled_ribs_before.mhd", m_Verbose);
1232 //---------------------------------
1233 // Write the filled ribs
1234 //---------------------------------
1235 if (m_ArgsInfo.writeFilledRibs_given) {
1236 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1237 caster->SetInput(filledRibs);
1238 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeFilledRibs_arg, m_Verbose);
1242 //----------------------------------------------------------------------------------------------------
1243 // Collapse to the lungs
1244 //----------------------------------------------------------------------------------------------------
1246 std::cout<<std::endl;
1247 std::cout<<"=========================================="<<std::endl;
1248 std::cout<<"|| Collapsing to the lung image ||"<<std::endl;
1249 std::cout<<"=========================================="<<std::endl;
1252 //---------------------------------
1253 // Make feature image air+bones
1254 //---------------------------------
1255 if (m_Verbose) std::cout<<"Making feature images..."<<std::endl;
1256 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
1257 setBackgroundFilter->SetInput(air);
1258 setBackgroundFilter->SetInput2(bones);
1259 setBackgroundFilter->SetMaskValue(0);
1260 setBackgroundFilter->SetOutsideValue(0);
1261 setBackgroundFilter->Update();
1262 typename InternalImageType::Pointer bones_air =setBackgroundFilter->GetOutput();
1264 typename SetBackgroundFilterType::Pointer setBackgroundFilter2 = SetBackgroundFilterType::New();
1265 setBackgroundFilter2->SetInput(bones_air);
1266 setBackgroundFilter2->SetInput2(lungs);
1267 setBackgroundFilter2->SetMaskValue(0);
1268 setBackgroundFilter2->SetOutsideValue(0);
1269 setBackgroundFilter2->Update();
1270 typename InternalImageType::Pointer lungs_bones_air =setBackgroundFilter2->GetOutput();
1272 //---------------------------------
1273 // Prepare previous solution
1274 //---------------------------------
1275 if (m_Verbose) std::cout<<"Resampling previous solution..."<<std::endl;
1276 typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
1277 typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
1278 typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
1279 resampler->SetOutputStartIndex(bones->GetLargestPossibleRegion().GetIndex());
1280 resampler->SetInput(filledRibs);
1281 resampler->SetInterpolator(interpolator);
1282 resampler->SetOutputParametersFromImage(bones);
1283 resampler->Update();
1284 filledRibs =resampler->GetOutput();
1286 if (m_Verbose) std::cout<<"Setting lungs to 1..."<<std::endl;
1287 typename SetBackgroundFilterType::Pointer setBackgroundFilter3 = SetBackgroundFilterType::New();
1288 setBackgroundFilter3->SetInput(filledRibs);
1289 setBackgroundFilter3->SetInput2(lungs);
1290 setBackgroundFilter3->SetMaskValue(0);
1291 setBackgroundFilter3->SetOutsideValue(1);
1292 setBackgroundFilter3->Update();
1293 filledRibs=setBackgroundFilter3->GetOutput();
1295 if (m_Verbose) std::cout<<"Removing overlap with bones..."<<std::endl;
1296 typename SetBackgroundFilterType::Pointer setBackgroundFilter4 = SetBackgroundFilterType::New();
1297 setBackgroundFilter4->SetInput(filledRibs);
1298 setBackgroundFilter4->SetInput2(bones);
1299 setBackgroundFilter4->SetMaskValue(0);
1300 setBackgroundFilter4->SetOutsideValue(0);
1301 setBackgroundFilter4->Update();
1302 filledRibs =setBackgroundFilter4->GetOutput();
1303 // writeImage<InternalImageType>(filledRibs,"/home/jef/tmp/filledRibs_pp.mhd");
1304 //---------------------------------
1305 // Calculate Distance Map
1306 //---------------------------------
1307 // typedef itk::ApproximateSignedDistanceMapImageFilter <InternalImageType, LevelSetImageType> DistanceMapImageFilterType2;
1308 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
1309 distanceMapImageFilter->SetInput(filledRibs);
1310 distanceMapImageFilter->SetInsideIsPositive(false);
1311 // distanceMapImageFilter->SetInsideValue(0);
1312 // distanceMapImageFilter->SetOutsideValue(1);
1313 if (m_Verbose) std::cout<<"Calculating distance map..."<<std::endl;
1314 distanceMapImageFilter->Update();
1316 //---------------------------------
1318 //---------------------------------
1319 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
1320 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
1321 levelSetFilter->SetFeatureImage( lungs_bones_air );
1322 levelSetFilter->SetPropagationScaling(m_ArgsInfo.prop3_arg);
1323 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve3_arg );
1324 levelSetFilter->SetAdvectionScaling( 0 );
1325 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS3_arg );
1326 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter3_arg );
1327 levelSetFilter->SetUseImageSpacing(true);
1330 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1333 if (m_Verbose) std::cout<<"Starting the levelset..."<<std::endl;
1334 unsigned int totalNumberOfIterations=0;
1336 levelSetFilter->Update();
1339 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1340 levelSetFilter->SetInput(levelSetFilter->GetOutput());
1341 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg);
1342 std::cout << "After "<<totalNumberOfIterations<<" iteration the RMS change is " << levelSetFilter->GetRMSChange() <<"..."<< std::endl;
1345 // std::ostringstream number_str;
1346 // number_str << number;
1347 // std::string param = number_str.str();
1348 // system((script+param).c_str());
1351 // stopping condition
1352 if ( (totalNumberOfIterations>= (unsigned int) m_ArgsInfo.maxIter3_arg) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS3_arg) ) break;
1353 levelSetFilter->SetNumberOfIterations( std::min ((unsigned int) m_ArgsInfo.iter3_arg, (unsigned int) m_ArgsInfo.maxIter3_arg-totalNumberOfIterations ) );
1358 std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1359 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1360 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1364 typedef itk::BinaryThresholdImageFilter< LevelSetImageType,InternalImageType> LevelSetBinarizeFilterType;
1365 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1366 thresholder->SetUpperThreshold( 0.0 );
1367 thresholder->SetOutsideValue( 0 );
1368 thresholder->SetInsideValue( 1 );
1369 thresholder->SetInput( levelSetFilter->GetOutput() );
1370 thresholder->Update();
1371 typename InternalImageType::Pointer output = thresholder->GetOutput();
1374 //----------------------------------------------------------------------------------------------------
1375 // Prepare the output
1376 //----------------------------------------------------------------------------------------------------
1378 //---------------------------------
1380 //---------------------------------
1381 if (m_Verbose) std::cout<<"Removing overlap mask with air..."<<std::endl;
1382 typename SetBackgroundFilterType::Pointer setBackgroundFilter5 = SetBackgroundFilterType::New();
1383 setBackgroundFilter5->SetInput(output);
1384 setBackgroundFilter5->SetInput2(air);
1385 setBackgroundFilter5->SetMaskValue(0);
1386 setBackgroundFilter5->SetOutsideValue(0);
1387 setBackgroundFilter5->Update();
1388 output=setBackgroundFilter5->GetOutput();
1390 //---------------------------------
1391 // Open & close to cleanup
1392 //---------------------------------
1393 if ( m_ArgsInfo.openClose_flag) {
1395 //---------------------------------
1396 // Structuring element
1397 //---------------------------------
1398 typedef itk::BinaryBallStructuringElement<InternalPixelType,Dimension > KernelType;
1399 KernelType structuringElement;
1400 structuringElement.SetRadius(1);
1401 structuringElement.CreateStructuringElement();
1403 //---------------------------------
1405 //---------------------------------
1406 typedef itk::BinaryMorphologicalOpeningImageFilter<InternalImageType, InternalImageType , KernelType> OpenFilterType;
1407 typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
1408 openFilter->SetInput(output);
1409 openFilter->SetBackgroundValue(0);
1410 openFilter->SetForegroundValue(1);
1411 openFilter->SetKernel(structuringElement);
1412 if(m_Verbose) std::cout<<"Opening the output image..."<<std::endl;
1414 //---------------------------------
1416 //---------------------------------
1417 typedef itk::BinaryMorphologicalClosingImageFilter<InternalImageType, InternalImageType , KernelType> CloseFilterType;
1418 typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
1419 closeFilter->SetInput(openFilter->GetOutput());
1420 closeFilter->SetSafeBorder(true);
1421 closeFilter->SetForegroundValue(1);
1422 closeFilter->SetKernel(structuringElement);
1423 if(m_Verbose) std::cout<<"Closing the output image..."<<std::endl;
1424 closeFilter->Update();
1425 output=closeFilter->GetOutput();
1428 //writeImage<InternalImageType>(output,"/home/jef/tmp/mm_double.mhd");
1430 // Extract the upper part
1431 typedef itk::CropImageFilter<InternalImageType, InternalImageType> CropImageFilterType;
1432 typename CropImageFilterType::Pointer cropFilter=CropImageFilterType::New();
1433 cropFilter->SetInput(output);
1434 typename InputImageType::SizeType lSize, uSize;
1437 lSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
1438 cropFilter->SetLowerBoundaryCropSize(lSize);
1439 cropFilter->SetUpperBoundaryCropSize(uSize);
1440 cropFilter->Update();
1443 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1444 caster->SetInput(cropFilter->GetOutput());
1445 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.output_arg, m_Verbose);
1451 #endif //#define clitkMotionMaskGenericFilter_txx