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 typename LabelType::RegionType lung_bbox = label->GetBoundingBox();
546 axes[0] = 0.95*lung_bbox.GetSize()[0]*spacing[0]/2;
547 axes[1] = 0.3*lung_bbox.GetSize()[1]*spacing[1]/2;
548 axes[2] = 1.25*fabs(lung_centroid[2] - center[2]);
550 // ellipse's center in XY is the center of the lungs' bounding box
551 typename InternalImageType::PointType origin = lungs_low->GetOrigin();
552 centerEllips[0] = origin[0] + (lung_bbox.GetIndex()[0] + lung_bbox.GetSize()[0]/2)*spacing[0];
553 centerEllips[1] = origin[1] + (lung_bbox.GetIndex()[1] + lung_bbox.GetSize()[1]/2)*spacing[1];
554 centerEllips[2] = center[2];
556 lungs_low->SetRegions(region_before);
559 std::cout << "Lungs'centroid at " << lung_centroid << std::endl;
560 std::cout << "Ellipse center at " << centerEllips << std::endl;
561 std::cout << "Ellipse axes as " << axes << std::endl;
562 std::cout << "Lung's bounding box (index,size) " << lung_bbox.GetIndex() << lung_bbox.GetSize() << std::endl;
566 //---------------------------------
567 // Offset from center
568 //---------------------------------
569 typename itk::Vector<double,Dimension> offset;
570 if(m_ArgsInfo.offset_given== Dimension) {
571 for(unsigned int i=0; i<Dimension; i++)
572 offset[i]=m_ArgsInfo.offset_arg[i];
577 centerEllips=center+offset;
579 std::cout<<"Placing the center of the initial ellipsoide at (mm) "<<std::endl;
580 std::cout<<centerEllips[0];
581 for (unsigned int i=1; i<Dimension; i++)
582 std::cout<<" "<<centerEllips[i];
583 std::cout<<std::endl;
587 if (m_ArgsInfo.axes_given==Dimension)
588 for(unsigned int i=0; i<Dimension; i++)
589 axes[i]=m_ArgsInfo.axes_arg[i];
597 //---------------------------------
599 //---------------------------------
600 ellips=InternalImageType::New();
601 ellips->SetRegions (bones_low->GetLargestPossibleRegion());
602 ellips->SetOrigin(bones_low->GetOrigin());
603 ellips->SetSpacing(bones_low->GetSpacing());
605 ellips->FillBuffer(0);
608 IteratorType itEllips(ellips, ellips->GetLargestPossibleRegion());
609 itEllips.GoToBegin();
610 typename InternalImageType::PointType point;
611 typename InternalImageType::IndexType index;
614 if (m_Verbose) std::cout<<"Drawing the initial ellipsoide..."<<std::endl;
615 while (!itEllips.IsAtEnd()) {
616 index=itEllips.GetIndex();
617 ellips->TransformIndexToPhysicalPoint(index, point);
619 for(unsigned int i=0; i<Dimension; i++)
620 distance+=powf( ( (centerEllips[i]-point[i])/axes[i] ), 2);
629 //---------------------------------
631 //---------------------------------
632 if (m_ArgsInfo.writeEllips_given) {
633 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
634 caster->SetInput(ellips);
635 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeEllips_arg, m_Verbose);
643 //-------------------------------------------------------------------
644 // Update with the number of dimensions and the pixeltype
645 //-------------------------------------------------------------------
646 template <unsigned int Dimension, class PixelType>
648 MotionMaskGenericFilter::UpdateWithDimAndPixelType()
652 typedef int InternalPixelType;
653 typedef itk::Image<PixelType, Dimension> InputImageType;
654 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
655 typedef itk::Image< float, Dimension> LevelSetImageType; //labeling doesn't work with unsigned char?
656 typedef itk::Image<unsigned char, Dimension> OutputImageType;
659 //----------------------------------------------------------------------------------------------------
660 // Typedef for Processing
661 //----------------------------------------------------------------------------------------------------
662 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
663 typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
664 typedef itk::BinaryThresholdImageFilter< LevelSetImageType,InternalImageType> LevelSetBinarizeFilterType;
665 typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> InversionFilterType;
666 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdFilterType;
667 typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
668 typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
669 typedef itk::MirrorPadImageFilter<InternalImageType, InternalImageType> MirrorPadImageFilterType;
670 typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
671 typedef itk::ResampleImageFilter<InternalImageType, InternalImageType> ResampleImageFilterType;
672 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
673 typedef itk::GeodesicActiveContourLevelSetImageFilter<LevelSetImageType, InternalImageType> LevelSetImageFilterType;
674 typedef itk::SignedDanielssonDistanceMapImageFilter<InternalImageType, LevelSetImageType> DistanceMapImageFilterType;
675 typedef clitk::SetBackgroundImageFilter<InternalImageType,InternalImageType, InternalImageType> SetBackgroundFilterType;
676 typedef itk::LabelStatisticsImageFilter<InternalImageType,InternalImageType> LabelStatisticsImageFilterType;
677 typedef itk::CastImageFilter<InternalImageType,OutputImageType> CastImageFilterType;
681 typedef itk::ImageFileReader<InputImageType> InputReaderType;
682 typename InputReaderType::Pointer reader = InputReaderType::New();
683 reader->SetFileName( m_InputFileName);
685 typename InputImageType::Pointer input= reader->GetOutput();
687 // // globals for avi
688 // unsigned int number=0;
692 std::cout<<std::endl;
693 std::cout<<"=========================================="<<std::endl;
694 std::cout<<"|| Making feature images ||"<<std::endl;
695 std::cout<<"=========================================="<<std::endl;
698 //--------------------------------------------------------------------------------
700 //-------------------------------------------------------------------------------
701 typename InternalImageType::Pointer lungs = this->GetLungsImage<Dimension, PixelType>(input);
703 //-------------------------------------------------------------------------------
705 //-------------------------------------------------------------------------------
706 typename InternalImageType::Pointer air = this->GetAirImage<Dimension, PixelType>(input, lungs);
708 //-------------------------------------------------------------------------------
710 //-------------------------------------------------------------------------------
711 typename InternalImageType::Pointer bones = this->GetBonesImage<Dimension, PixelType>(input);
713 //----------------------------------------------------------------------------------------------------
714 // Write feature images
715 //----------------------------------------------------------------------------------------------------
716 if(m_ArgsInfo.writeFeature_given) {
717 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
718 setBackgroundFilter->SetInput(air);
719 setBackgroundFilter->SetInput2(bones);
720 setBackgroundFilter->SetMaskValue(0);
721 setBackgroundFilter->SetOutsideValue(2);
722 setBackgroundFilter->Update();
723 typename InternalImageType::Pointer bones_air =setBackgroundFilter->GetOutput();
725 typename SetBackgroundFilterType::Pointer setBackgroundFilter2 = SetBackgroundFilterType::New();
726 setBackgroundFilter2->SetInput(bones_air);
727 setBackgroundFilter2->SetInput2(lungs);
728 setBackgroundFilter2->SetMaskValue(0);
729 setBackgroundFilter2->SetOutsideValue(3);
730 setBackgroundFilter2->Update();
731 typename InternalImageType::Pointer lungs_bones_air =setBackgroundFilter2->GetOutput();
733 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
734 caster->SetInput(lungs_bones_air);
735 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeFeature_arg, m_Verbose);
738 //----------------------------------------------------------------------------------------------------
739 // Low dimensional versions
740 //----------------------------------------------------------------------------------------------------
741 typename InternalImageType::Pointer bones_low =Resample<Dimension , PixelType>(bones);
742 typename InternalImageType::Pointer lungs_low =Resample<Dimension , PixelType>(lungs);
743 typename InternalImageType::Pointer air_low =Resample<Dimension , PixelType>(air);
745 //---------------------------------
747 //---------------------------------
748 if(m_ArgsInfo.pad_flag) {
749 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
750 IteratorType it(air_low, air_low->GetLargestPossibleRegion());
751 typename InternalImageType::IndexType index;
752 while(!it.IsAtEnd()) {
754 for (unsigned int i=0; i<Dimension; i++)
755 if(index[i]==air_low->GetLargestPossibleRegion().GetIndex()[i]
756 || index[i]==(unsigned int)air_low->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air_low->GetLargestPossibleRegion().GetSize()[i]-1)
762 //---------------------------------
764 //---------------------------------
765 typename itk::Vector<double,Dimension> center;
766 typedef itk::ImageMomentsCalculator<InternalImageType> MomentsCalculatorType;
767 typename MomentsCalculatorType::Pointer momentsCalculator=MomentsCalculatorType::New();
768 momentsCalculator->SetImage(air);
769 momentsCalculator->Compute();
770 center=momentsCalculator->GetCenterOfGravity();
772 std::cout<<"The center of gravity of the patient body is located at (mm) "<<std::endl;
773 std::cout<<center[0];
774 for (unsigned int i=1; i<Dimension; i++)
775 std::cout<<" "<<center[i];
776 std::cout<<std::endl;
779 //----------------------------------------------------------------------------------------------------
780 // Ellipsoide initialization
781 //----------------------------------------------------------------------------------------------------
782 typename InternalImageType::Pointer m_Ellips= InitializeEllips<Dimension, PixelType>(center,bones_low,lungs_low);
784 //----------------------------------------------------------------------------------------------------
786 //----------------------------------------------------------------------------------------------------
787 typename InternalImageType::Pointer grownEllips;
788 if (m_ArgsInfo.grownEllips_given || m_ArgsInfo.filledRibs_given) {
789 if (m_ArgsInfo.grownEllips_given) {
791 typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
792 featureReader->SetFileName(m_ArgsInfo.grownEllips_arg);
793 featureReader->Update();
794 grownEllips=featureReader->GetOutput();
799 std::cout<<std::endl;
800 std::cout<<"=========================================="<<std::endl;
801 std::cout<<"|| Growing ellipsoide ||"<<std::endl;
802 std::cout<<"=========================================="<<std::endl;
805 //---------------------------------
807 //---------------------------------
808 typename InternalImageType::PointType dPoint;
809 if (m_ArgsInfo.detectionPoint_given) {
810 for (unsigned int i=0; i<Dimension; i++)
811 dPoint[i]=m_ArgsInfo.detectionPoint_arg[i];
813 typename InternalImageType::RegionType searchRegion=air->GetLargestPossibleRegion();
814 typename InternalImageType::SizeType searchRegionSize = searchRegion.GetSize();
815 typename InternalImageType::IndexType searchRegionIndex = searchRegion.GetIndex();
816 searchRegionIndex[2]+=searchRegionSize[2]/2;
817 searchRegionSize[2]=1;
818 searchRegion.SetSize(searchRegionSize);
819 searchRegion.SetIndex(searchRegionIndex);
820 IteratorType itAbdomen(air, searchRegion);
821 itAbdomen.GoToBegin();
823 typename InternalImageType::PointType aPoint;
824 typename InternalImageType::IndexType aIndex;
826 if (m_Verbose) std::cout<<"Detecting the abdomen..."<<std::endl;
827 while (!itAbdomen.IsAtEnd()) {
828 if(itAbdomen.Get()) break;
831 aIndex=itAbdomen.GetIndex();
832 air->TransformIndexToPhysicalPoint(aIndex,aPoint);
833 if (m_Verbose) std::cout<<"Detected the abdomen at "<<aPoint<<".."<<std::endl;
836 //---------------------------------
837 // Detect abdomen in additional images?
838 //---------------------------------
839 if (m_ArgsInfo.detectionPairs_given) {
840 for (unsigned int i=0; i<m_ArgsInfo.detectionPairs_given; i++) {
841 typename InternalImageType::Pointer airAdd;
842 //---------------------------------
844 //--------------------------------
845 typedef itk::ImageFileReader<InputImageType> InputReaderType;
846 typename InputReaderType::Pointer reader = InputReaderType::New();
847 reader->SetFileName( m_ArgsInfo.detectionPairs_arg[i]);
849 typename InputImageType::Pointer additional= reader->GetOutput();
850 if (m_Verbose) std::cout<<"Extracting the air from additional image "<< i<<"..."<<std::endl;
852 //---------------------------------
853 // Binarize the image
854 //---------------------------------
855 typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
856 binarizeFilter->SetInput(additional);
857 binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdAir_arg));
858 binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdAir_arg));
859 if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<<m_ArgsInfo.lowerThresholdAir_arg
860 <<", "<<m_ArgsInfo.upperThresholdAir_arg<<"..."<<std::endl;
862 //---------------------------------
863 // Label the connected components
864 //---------------------------------
865 typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
866 connectFilter->SetInput(binarizeFilter->GetOutput());
867 connectFilter->SetBackgroundValue(0);
868 connectFilter->SetFullyConnected(false);
869 if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
871 //---------------------------------
872 // Sort the labels according to size
873 //---------------------------------
874 typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
875 relabelFilter->SetInput(connectFilter->GetOutput());
876 if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
878 //---------------------------------
880 //---------------------------------
881 typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
882 thresholdFilter->SetInput(relabelFilter->GetOutput());
883 thresholdFilter->SetUpper(1);
884 if (m_Verbose) std::cout<<"Keeping the first label..."<<std::endl;
885 thresholdFilter->Update();
886 airAdd=thresholdFilter->GetOutput();
889 //---------------------------------
891 //---------------------------------
892 typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
893 inversionFilter->SetInput(airAdd);
894 inversionFilter->SetLowerThreshold(0);
895 inversionFilter->SetUpperThreshold(0);
896 inversionFilter->SetInsideValue(1);
897 inversionFilter->Update();
899 //---------------------------------
901 //---------------------------------
902 typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
903 mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
904 typename InternalImageType::SizeType padSize;
906 padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
907 mirrorPadImageFilter->SetPadLowerBound(padSize);
908 if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
909 mirrorPadImageFilter->Update();
910 airAdd=mirrorPadImageFilter->GetOutput();
912 //---------------------------------
914 //---------------------------------
915 typename InternalImageType::RegionType searchRegion=airAdd->GetLargestPossibleRegion();
916 typename InternalImageType::SizeType searchRegionSize = searchRegion.GetSize();
917 typename InternalImageType::IndexType searchRegionIndex = searchRegion.GetIndex();
918 searchRegionIndex[2]+=searchRegionSize[2]/2;
919 searchRegionSize[2]=1;
920 searchRegion.SetSize(searchRegionSize);
921 searchRegion.SetIndex(searchRegionIndex);
922 IteratorType itAbdomen(airAdd, searchRegion);
923 itAbdomen.GoToBegin();
925 typename InternalImageType::PointType additionalPoint;
926 typename InternalImageType::IndexType aIndex;
928 if (m_Verbose) std::cout<<"Detecting the abdomen..."<<std::endl;
929 while (!itAbdomen.IsAtEnd()) {
930 if(itAbdomen.Get()) break;
933 aIndex=itAbdomen.GetIndex();
934 airAdd->TransformIndexToPhysicalPoint(aIndex,additionalPoint);
935 if (m_Verbose) std::cout<<"Detected the abdomen in the additional image at "<<additionalPoint<<".."<<std::endl;
937 if(additionalPoint[1]< aPoint[1]) {
938 aPoint=additionalPoint;
939 if (m_Verbose) std::cout<<"Modifying the detected abdomen to "<<aPoint<<".."<<std::endl;
946 // Determine the detection point
950 if(m_ArgsInfo.offsetDetect_given==Dimension)
951 for(unsigned int i=0; i <Dimension; i++)
952 dPoint[i]+=m_ArgsInfo.offsetDetect_arg[i];
957 if (m_Verbose) std::cout<<"Setting the detection point to "<<dPoint<<".."<<std::endl;
960 //---------------------------------
961 // Pad the rib image and ellips image
962 //---------------------------------
963 typename InternalImageType::Pointer padded_ellips;
964 typename InternalImageType::Pointer padded_bones_low;
966 // If detection point not inside the image: pad
967 typename InternalImageType::IndexType dIndex;
968 if (!bones_low->TransformPhysicalPointToIndex(dPoint, dIndex)) {
969 typename InternalImageType::SizeType padSize;
971 padSize[1]=abs(dIndex[1])+1;
972 if (m_Verbose) std::cout<<"Padding the images with padding size "<<padSize<<" to include the detection point..."<<std::endl;
974 typename MirrorPadImageFilterType::Pointer padBonesFilter=MirrorPadImageFilterType::New();
975 padBonesFilter->SetInput(bones_low);
976 padBonesFilter->SetPadLowerBound(padSize);
977 padBonesFilter->Update();
978 padded_bones_low=padBonesFilter->GetOutput();
980 typename MirrorPadImageFilterType::Pointer padEllipsFilter=MirrorPadImageFilterType::New();
981 padEllipsFilter->SetInput(m_Ellips);
982 padEllipsFilter->SetPadLowerBound(padSize);
983 padEllipsFilter->Update();
984 padded_ellips=padEllipsFilter->GetOutput();
988 padded_bones_low=bones_low;
989 padded_ellips=m_Ellips;
993 //---------------------------------
994 // Calculate distance map
995 //---------------------------------
996 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
997 distanceMapImageFilter->SetInput(padded_ellips);
998 distanceMapImageFilter->SetInsideIsPositive(false);
999 if (m_Verbose) std::cout<<"Calculating the distance map..."<<std::endl;
1000 distanceMapImageFilter->Update();
1001 if (m_ArgsInfo.writeDistMap_given) {
1002 writeImage<LevelSetImageType>(distanceMapImageFilter->GetOutput(), m_ArgsInfo.writeDistMap_arg, m_Verbose);
1006 //---------------------------------
1007 // Grow while monitoring detection point
1008 //---------------------------------
1009 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
1010 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
1011 levelSetFilter->SetFeatureImage( padded_bones_low );
1012 levelSetFilter->SetPropagationScaling( 1 );
1013 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve1_arg );
1014 levelSetFilter->SetAdvectionScaling( 0 );
1015 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS1_arg );
1016 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter1_arg );
1017 levelSetFilter->SetUseImageSpacing(true);
1019 // //---------------------------------
1020 // // Script for making movie
1021 // //---------------------------------
1022 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1025 if (m_Verbose) std::cout<<"Starting level set segmentation..."<<std::endl;
1026 unsigned int totalNumberOfIterations=0;
1028 levelSetFilter->Update();
1029 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1031 if ( levelSetFilter->GetOutput()->GetPixel(dIndex) < 0 ) {
1032 if (m_Verbose) std::cout<<"Detection point reached!"<<std::endl;
1035 if (m_Verbose) std::cout<<"Detection point not reached after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1036 levelSetFilter->SetInput(levelSetFilter->GetOutput());
1037 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose);
1040 // std::ostringstream number_str;
1041 // number_str << number;
1042 // std::string param = number_str.str();
1043 // system((script+param).c_str());
1047 if ( (totalNumberOfIterations> 10000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS1_arg) ) break;
1051 if (m_Verbose) std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1052 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1053 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1056 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1057 thresholder->SetUpperThreshold( 0.0 );
1058 thresholder->SetOutsideValue( 0 );
1059 thresholder->SetInsideValue( 1 );
1060 thresholder->SetInput( levelSetFilter->GetOutput() );
1061 if (m_Verbose) std::cout<<"Thresholding the output level set..."<<std::endl;
1062 thresholder->Update();
1063 grownEllips=thresholder->GetOutput();
1066 //---------------------------------
1067 // Write the grown ellips
1068 //---------------------------------
1069 if (m_ArgsInfo.writeGrownEllips_given) {
1070 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1071 caster->SetInput(grownEllips);
1072 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeGrownEllips_arg, m_Verbose);
1076 //----------------------------------------------------------------------------------------------------
1078 //----------------------------------------------------------------------------------------------------
1079 typename InternalImageType::Pointer filledRibs;
1080 if (m_ArgsInfo.filledRibs_given) {
1081 typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
1082 featureReader->SetFileName(m_ArgsInfo.filledRibs_arg);
1083 if (m_Verbose) std::cout<<"Reading filled ribs mask..."<<std::endl;
1084 featureReader->Update();
1085 filledRibs=featureReader->GetOutput();
1088 std::cout<<std::endl;
1089 std::cout<<"=========================================="<<std::endl;
1090 std::cout<<"|| Filling the ribs image ||"<<std::endl;
1091 std::cout<<"=========================================="<<std::endl;
1093 //---------------------------------
1094 // Make feature image air+bones
1095 //---------------------------------
1096 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
1097 setBackgroundFilter->SetInput(air_low);
1098 setBackgroundFilter->SetInput2(bones_low);
1099 setBackgroundFilter->SetMaskValue(0);
1100 setBackgroundFilter->SetOutsideValue(0);
1101 setBackgroundFilter->Update();
1102 typename InternalImageType::Pointer bones_air_low =setBackgroundFilter->GetOutput();
1104 //---------------------------------
1105 // Resampling previous solution
1106 //---------------------------------
1107 if (m_Verbose) std::cout<<"Resampling previous solution..."<<std::endl;
1108 typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
1109 typename InternalImageType::PointType origin;
1110 bones_low->TransformIndexToPhysicalPoint(bones_low->GetLargestPossibleRegion().GetIndex(), origin);
1111 resampler->SetOutputOrigin(origin);
1112 resampler->SetOutputSpacing(bones_low->GetSpacing());
1113 typename InternalImageType::SizeType size_low= bones_low->GetLargestPossibleRegion().GetSize();
1114 resampler->SetSize(size_low);
1115 typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
1116 resampler->SetInterpolator(interpolator);
1117 resampler->SetInput(grownEllips);
1118 resampler->Update();
1119 typename InternalImageType::Pointer grownEllips =resampler->GetOutput();
1122 //---------------------------------
1123 // Calculate Distance Map
1124 //---------------------------------
1125 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
1126 distanceMapImageFilter->SetInput(grownEllips);
1127 distanceMapImageFilter->SetInsideIsPositive(false);
1128 if (m_Verbose) std::cout<<"Calculating distance map..."<<std::endl;
1129 distanceMapImageFilter->Update();
1130 if (m_ArgsInfo.writeDistMap_given) {
1131 writeImage<LevelSetImageType>(distanceMapImageFilter->GetOutput(), m_ArgsInfo.writeDistMap_arg, m_Verbose);
1135 //---------------------------------
1136 // Grow while monitoring lung volume coverage
1137 //---------------------------------
1138 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
1139 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
1140 levelSetFilter->SetFeatureImage( bones_air_low );
1141 levelSetFilter->SetPropagationScaling( 1 );
1142 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve2_arg );
1143 levelSetFilter->SetAdvectionScaling( 0 );
1144 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS2_arg );
1145 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter2_arg );
1146 levelSetFilter->SetUseImageSpacing(true);
1148 //---------------------------------
1150 //---------------------------------
1151 typename InversionFilterType::Pointer invertor= InversionFilterType::New();
1152 invertor->SetInput(lungs_low);
1153 invertor->SetLowerThreshold(0);
1154 invertor->SetUpperThreshold(0);
1155 invertor->SetInsideValue(1);
1157 typename InternalImageType::Pointer lungs_low_inv=invertor->GetOutput();
1159 //---------------------------------
1160 // Calculate lung volume
1161 //---------------------------------
1162 typename LabelStatisticsImageFilterType::Pointer statisticsFilter= LabelStatisticsImageFilterType::New();
1163 statisticsFilter->SetInput(lungs_low_inv);
1164 statisticsFilter->SetLabelInput(lungs_low);
1165 statisticsFilter->Update();
1166 unsigned int volume=statisticsFilter->GetSum(0);
1168 // Prepare threshold
1169 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1170 thresholder->SetUpperThreshold( 0.0 );
1171 thresholder->SetOutsideValue( 0 );
1172 thresholder->SetInsideValue( 1 );
1176 unsigned int totalNumberOfIterations=0;
1177 unsigned int coverage=0;
1180 // //---------------------------------
1181 // // Script for making movie
1182 // //---------------------------------
1183 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1186 levelSetFilter->Update();
1187 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1189 thresholder->SetInput( levelSetFilter->GetOutput() );
1190 thresholder->Update();
1191 statisticsFilter->SetInput(thresholder->GetOutput());
1192 statisticsFilter->Update();
1193 coverage=statisticsFilter->GetSum(0);
1195 // Compare the volumes
1196 if ( coverage >= m_ArgsInfo.fillingLevel_arg * volume/100 ) {
1197 if (m_Verbose) std::cout<<"Lungs filled for "<< m_ArgsInfo.fillingLevel_arg<<"% !"<<std::endl;
1200 if (m_Verbose) std::cout<<"After "<<totalNumberOfIterations<<" iterations, lungs are covered for "
1201 <<(double)coverage/volume*100.0<<"%..."<<std::endl;
1202 levelSetFilter->SetInput(levelSetFilter->GetOutput());
1203 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose);
1206 // std::ostringstream number_str;
1207 // number_str << number;
1208 // std::string param = number_str.str();
1209 // system((script+param).c_str());
1213 if ( (totalNumberOfIterations> 30000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS2_arg) ) break;
1217 if (m_Verbose) std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1218 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1219 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1222 thresholder->SetInput( levelSetFilter->GetOutput() );
1223 thresholder->Update();
1224 filledRibs=thresholder->GetOutput();
1225 // writeImage<InternalImageType>(filledRibs,"/home/jef/tmp/filled_ribs_before.mhd", m_Verbose);
1229 //---------------------------------
1230 // Write the filled ribs
1231 //---------------------------------
1232 if (m_ArgsInfo.writeFilledRibs_given) {
1233 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1234 caster->SetInput(filledRibs);
1235 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeFilledRibs_arg, m_Verbose);
1239 //----------------------------------------------------------------------------------------------------
1240 // Collapse to the lungs
1241 //----------------------------------------------------------------------------------------------------
1243 std::cout<<std::endl;
1244 std::cout<<"=========================================="<<std::endl;
1245 std::cout<<"|| Collapsing to the lung image ||"<<std::endl;
1246 std::cout<<"=========================================="<<std::endl;
1249 //---------------------------------
1250 // Make feature image air+bones
1251 //---------------------------------
1252 if (m_Verbose) std::cout<<"Making feature images..."<<std::endl;
1253 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
1254 setBackgroundFilter->SetInput(air);
1255 setBackgroundFilter->SetInput2(bones);
1256 setBackgroundFilter->SetMaskValue(0);
1257 setBackgroundFilter->SetOutsideValue(0);
1258 setBackgroundFilter->Update();
1259 typename InternalImageType::Pointer bones_air =setBackgroundFilter->GetOutput();
1261 typename SetBackgroundFilterType::Pointer setBackgroundFilter2 = SetBackgroundFilterType::New();
1262 setBackgroundFilter2->SetInput(bones_air);
1263 setBackgroundFilter2->SetInput2(lungs);
1264 setBackgroundFilter2->SetMaskValue(0);
1265 setBackgroundFilter2->SetOutsideValue(0);
1266 setBackgroundFilter2->Update();
1267 typename InternalImageType::Pointer lungs_bones_air =setBackgroundFilter2->GetOutput();
1269 //---------------------------------
1270 // Prepare previous solution
1271 //---------------------------------
1272 if (m_Verbose) std::cout<<"Resampling previous solution..."<<std::endl;
1273 typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
1274 typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
1275 typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
1276 resampler->SetOutputStartIndex(bones->GetLargestPossibleRegion().GetIndex());
1277 resampler->SetInput(filledRibs);
1278 resampler->SetInterpolator(interpolator);
1279 resampler->SetOutputParametersFromImage(bones);
1280 resampler->Update();
1281 filledRibs =resampler->GetOutput();
1283 if (m_Verbose) std::cout<<"Setting lungs to 1..."<<std::endl;
1284 typename SetBackgroundFilterType::Pointer setBackgroundFilter3 = SetBackgroundFilterType::New();
1285 setBackgroundFilter3->SetInput(filledRibs);
1286 setBackgroundFilter3->SetInput2(lungs);
1287 setBackgroundFilter3->SetMaskValue(0);
1288 setBackgroundFilter3->SetOutsideValue(1);
1289 setBackgroundFilter3->Update();
1290 filledRibs=setBackgroundFilter3->GetOutput();
1292 if (m_Verbose) std::cout<<"Removing overlap with bones..."<<std::endl;
1293 typename SetBackgroundFilterType::Pointer setBackgroundFilter4 = SetBackgroundFilterType::New();
1294 setBackgroundFilter4->SetInput(filledRibs);
1295 setBackgroundFilter4->SetInput2(bones);
1296 setBackgroundFilter4->SetMaskValue(0);
1297 setBackgroundFilter4->SetOutsideValue(0);
1298 setBackgroundFilter4->Update();
1299 filledRibs =setBackgroundFilter4->GetOutput();
1300 // writeImage<InternalImageType>(filledRibs,"/home/jef/tmp/filledRibs_pp.mhd");
1301 //---------------------------------
1302 // Calculate Distance Map
1303 //---------------------------------
1304 // typedef itk::ApproximateSignedDistanceMapImageFilter <InternalImageType, LevelSetImageType> DistanceMapImageFilterType2;
1305 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
1306 distanceMapImageFilter->SetInput(filledRibs);
1307 distanceMapImageFilter->SetInsideIsPositive(false);
1308 // distanceMapImageFilter->SetInsideValue(0);
1309 // distanceMapImageFilter->SetOutsideValue(1);
1310 if (m_Verbose) std::cout<<"Calculating distance map..."<<std::endl;
1311 distanceMapImageFilter->Update();
1313 //---------------------------------
1315 //---------------------------------
1316 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
1317 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
1318 levelSetFilter->SetFeatureImage( lungs_bones_air );
1319 levelSetFilter->SetPropagationScaling(m_ArgsInfo.prop3_arg);
1320 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve3_arg );
1321 levelSetFilter->SetAdvectionScaling( 0 );
1322 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS3_arg );
1323 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter3_arg );
1324 levelSetFilter->SetUseImageSpacing(true);
1327 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1330 if (m_Verbose) std::cout<<"Starting the levelset..."<<std::endl;
1331 unsigned int totalNumberOfIterations=0;
1333 levelSetFilter->Update();
1336 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1337 levelSetFilter->SetInput(levelSetFilter->GetOutput());
1338 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg);
1339 std::cout << "After "<<totalNumberOfIterations<<" iteration the RMS change is " << levelSetFilter->GetRMSChange() <<"..."<< std::endl;
1342 // std::ostringstream number_str;
1343 // number_str << number;
1344 // std::string param = number_str.str();
1345 // system((script+param).c_str());
1348 // stopping condition
1349 if ( (totalNumberOfIterations>= (unsigned int) m_ArgsInfo.maxIter3_arg) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS3_arg) ) break;
1350 levelSetFilter->SetNumberOfIterations( std::min ((unsigned int) m_ArgsInfo.iter3_arg, (unsigned int) m_ArgsInfo.maxIter3_arg-totalNumberOfIterations ) );
1355 std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1356 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1357 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1361 typedef itk::BinaryThresholdImageFilter< LevelSetImageType,InternalImageType> LevelSetBinarizeFilterType;
1362 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1363 thresholder->SetUpperThreshold( 0.0 );
1364 thresholder->SetOutsideValue( 0 );
1365 thresholder->SetInsideValue( 1 );
1366 thresholder->SetInput( levelSetFilter->GetOutput() );
1367 thresholder->Update();
1368 typename InternalImageType::Pointer output = thresholder->GetOutput();
1371 //----------------------------------------------------------------------------------------------------
1372 // Prepare the output
1373 //----------------------------------------------------------------------------------------------------
1375 //---------------------------------
1377 //---------------------------------
1378 if (m_Verbose) std::cout<<"Removing overlap mask with air..."<<std::endl;
1379 typename SetBackgroundFilterType::Pointer setBackgroundFilter5 = SetBackgroundFilterType::New();
1380 setBackgroundFilter5->SetInput(output);
1381 setBackgroundFilter5->SetInput2(air);
1382 setBackgroundFilter5->SetMaskValue(0);
1383 setBackgroundFilter5->SetOutsideValue(0);
1384 setBackgroundFilter5->Update();
1385 output=setBackgroundFilter5->GetOutput();
1387 //---------------------------------
1388 // Open & close to cleanup
1389 //---------------------------------
1390 if ( m_ArgsInfo.openClose_flag) {
1392 //---------------------------------
1393 // Structuring element
1394 //---------------------------------
1395 typedef itk::BinaryBallStructuringElement<InternalPixelType,Dimension > KernelType;
1396 KernelType structuringElement;
1397 structuringElement.SetRadius(1);
1398 structuringElement.CreateStructuringElement();
1400 //---------------------------------
1402 //---------------------------------
1403 typedef itk::BinaryMorphologicalOpeningImageFilter<InternalImageType, InternalImageType , KernelType> OpenFilterType;
1404 typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
1405 openFilter->SetInput(output);
1406 openFilter->SetBackgroundValue(0);
1407 openFilter->SetForegroundValue(1);
1408 openFilter->SetKernel(structuringElement);
1409 if(m_Verbose) std::cout<<"Opening the output image..."<<std::endl;
1411 //---------------------------------
1413 //---------------------------------
1414 typedef itk::BinaryMorphologicalClosingImageFilter<InternalImageType, InternalImageType , KernelType> CloseFilterType;
1415 typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
1416 closeFilter->SetInput(openFilter->GetOutput());
1417 closeFilter->SetSafeBorder(true);
1418 closeFilter->SetForegroundValue(1);
1419 closeFilter->SetKernel(structuringElement);
1420 if(m_Verbose) std::cout<<"Closing the output image..."<<std::endl;
1421 closeFilter->Update();
1422 output=closeFilter->GetOutput();
1425 //writeImage<InternalImageType>(output,"/home/jef/tmp/mm_double.mhd");
1427 // Extract the upper part
1428 typedef itk::CropImageFilter<InternalImageType, InternalImageType> CropImageFilterType;
1429 typename CropImageFilterType::Pointer cropFilter=CropImageFilterType::New();
1430 cropFilter->SetInput(output);
1431 typename InputImageType::SizeType lSize, uSize;
1434 lSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
1435 cropFilter->SetLowerBoundaryCropSize(lSize);
1436 cropFilter->SetUpperBoundaryCropSize(uSize);
1437 cropFilter->Update();
1440 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1441 caster->SetInput(cropFilter->GetOutput());
1442 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.output_arg, m_Verbose);
1448 #endif //#define clitkMotionMaskGenericFilter_txx