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 // compute statistics of the (mirroed) lung region in order to guess the params of the ellipse
524 typedef itk::LabelImageToShapeLabelMapFilter<InternalImageType> LabelImageToShapeLabelMapFilter;
525 typename LabelImageToShapeLabelMapFilter::Pointer label_filter = LabelImageToShapeLabelMapFilter::New();
526 label_filter->SetInput(lungs_low);
527 label_filter->Update();
529 typename InternalImageType::SpacingType spacing = lungs_low->GetSpacing();
530 typedef typename LabelImageToShapeLabelMapFilter::OutputImageType::LabelObjectType LabelType;
531 const LabelType* label = dynamic_cast<LabelType*>(label_filter->GetOutput()->GetNthLabelObject(0));
534 // try to guess ideal ellipse axes from the lungs' bounding box
535 #if ITK_VERSION_MAJOR >= 4
536 typename LabelType::RegionType lung_bbox = label->GetBoundingBox();
538 typename LabelType::RegionType lung_bbox = label->GetRegion();
540 axes[0] = 0.95*lung_bbox.GetSize()[0]*spacing[0]/2;
541 axes[1] = 0.3*lung_bbox.GetSize()[1]*spacing[1]/2;
542 axes[2] = 0.95*lung_bbox.GetSize()[2]*spacing[2]/2;
544 // ellipse's center is the center of the lungs' bounding box
545 typename InternalImageType::PointType origin = lungs_low->GetOrigin();
546 centerEllips[0] = origin[0] + (lung_bbox.GetIndex()[0] + lung_bbox.GetSize()[0]/2)*spacing[0];
547 centerEllips[1] = origin[1] + (lung_bbox.GetIndex()[1] + lung_bbox.GetSize()[1]/2)*spacing[1];
548 centerEllips[2] = origin[2] + (lung_bbox.GetIndex()[2] + lung_bbox.GetSize()[2]/2)*spacing[2];
551 std::cout << "Ellipse center at " << centerEllips << std::endl;
552 std::cout << "Ellipse axes as " << axes << std::endl;
553 std::cout << "Lung's bounding box (index,size) " << lung_bbox.GetIndex() << lung_bbox.GetSize() << std::endl;
557 //---------------------------------
558 // Offset from center
559 //---------------------------------
560 typename itk::Vector<double,Dimension> offset;
561 if(m_ArgsInfo.offset_given== Dimension) {
562 for(unsigned int i=0; i<Dimension; i++)
563 offset[i]=m_ArgsInfo.offset_arg[i];
568 centerEllips=center+offset;
570 std::cout<<"Placing the center of the initial ellipsoide at (mm) "<<std::endl;
571 std::cout<<centerEllips[0];
572 for (unsigned int i=1; i<Dimension; i++)
573 std::cout<<" "<<centerEllips[i];
574 std::cout<<std::endl;
578 if (m_ArgsInfo.axes_given==Dimension)
579 for(unsigned int i=0; i<Dimension; i++)
580 axes[i]=m_ArgsInfo.axes_arg[i];
588 //---------------------------------
590 //---------------------------------
591 ellips=InternalImageType::New();
592 ellips->SetRegions (bones_low->GetLargestPossibleRegion());
593 ellips->SetOrigin(bones_low->GetOrigin());
594 ellips->SetSpacing(bones_low->GetSpacing());
596 ellips->FillBuffer(0);
599 IteratorType itEllips(ellips, ellips->GetLargestPossibleRegion());
600 itEllips.GoToBegin();
601 typename InternalImageType::PointType point;
602 typename InternalImageType::IndexType index;
605 if (m_Verbose) std::cout<<"Drawing the initial ellipsoide..."<<std::endl;
606 while (!itEllips.IsAtEnd()) {
607 index=itEllips.GetIndex();
608 ellips->TransformIndexToPhysicalPoint(index, point);
610 for(unsigned int i=0; i<Dimension; i++)
611 distance+=powf( ( (centerEllips[i]-point[i])/axes[i] ), 2);
620 //---------------------------------
622 //---------------------------------
623 if (m_ArgsInfo.writeEllips_given) {
624 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
625 caster->SetInput(ellips);
626 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeEllips_arg, m_Verbose);
634 //-------------------------------------------------------------------
635 // Update with the number of dimensions and the pixeltype
636 //-------------------------------------------------------------------
637 template <unsigned int Dimension, class PixelType>
639 MotionMaskGenericFilter::UpdateWithDimAndPixelType()
643 typedef int InternalPixelType;
644 typedef itk::Image<PixelType, Dimension> InputImageType;
645 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
646 typedef itk::Image< float, Dimension> LevelSetImageType; //labeling doesn't work with unsigned char?
647 typedef itk::Image<unsigned char, Dimension> OutputImageType;
650 //----------------------------------------------------------------------------------------------------
651 // Typedef for Processing
652 //----------------------------------------------------------------------------------------------------
653 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
654 typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
655 typedef itk::BinaryThresholdImageFilter< LevelSetImageType,InternalImageType> LevelSetBinarizeFilterType;
656 typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> InversionFilterType;
657 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdFilterType;
658 typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
659 typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
660 typedef itk::MirrorPadImageFilter<InternalImageType, InternalImageType> MirrorPadImageFilterType;
661 typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
662 typedef itk::ResampleImageFilter<InternalImageType, InternalImageType> ResampleImageFilterType;
663 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
664 typedef itk::GeodesicActiveContourLevelSetImageFilter<LevelSetImageType, InternalImageType> LevelSetImageFilterType;
665 typedef itk::SignedDanielssonDistanceMapImageFilter<InternalImageType, LevelSetImageType> DistanceMapImageFilterType;
666 typedef clitk::SetBackgroundImageFilter<InternalImageType,InternalImageType, InternalImageType> SetBackgroundFilterType;
667 typedef itk::LabelStatisticsImageFilter<InternalImageType,InternalImageType> LabelStatisticsImageFilterType;
668 typedef itk::CastImageFilter<InternalImageType,OutputImageType> CastImageFilterType;
672 typedef itk::ImageFileReader<InputImageType> InputReaderType;
673 typename InputReaderType::Pointer reader = InputReaderType::New();
674 reader->SetFileName( m_InputFileName);
676 typename InputImageType::Pointer input= reader->GetOutput();
678 // // globals for avi
679 // unsigned int number=0;
683 std::cout<<std::endl;
684 std::cout<<"=========================================="<<std::endl;
685 std::cout<<"|| Making feature images ||"<<std::endl;
686 std::cout<<"=========================================="<<std::endl;
689 //--------------------------------------------------------------------------------
691 //-------------------------------------------------------------------------------
692 typename InternalImageType::Pointer lungs = this->GetLungsImage<Dimension, PixelType>(input);
694 //-------------------------------------------------------------------------------
696 //-------------------------------------------------------------------------------
697 typename InternalImageType::Pointer air = this->GetAirImage<Dimension, PixelType>(input, lungs);
699 //-------------------------------------------------------------------------------
701 //-------------------------------------------------------------------------------
702 typename InternalImageType::Pointer bones = this->GetBonesImage<Dimension, PixelType>(input);
704 //----------------------------------------------------------------------------------------------------
705 // Write feature images
706 //----------------------------------------------------------------------------------------------------
707 if(m_ArgsInfo.writeFeature_given) {
708 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
709 setBackgroundFilter->SetInput(air);
710 setBackgroundFilter->SetInput2(bones);
711 setBackgroundFilter->SetMaskValue(0);
712 setBackgroundFilter->SetOutsideValue(2);
713 setBackgroundFilter->Update();
714 typename InternalImageType::Pointer bones_air =setBackgroundFilter->GetOutput();
716 typename SetBackgroundFilterType::Pointer setBackgroundFilter2 = SetBackgroundFilterType::New();
717 setBackgroundFilter2->SetInput(bones_air);
718 setBackgroundFilter2->SetInput2(lungs);
719 setBackgroundFilter2->SetMaskValue(0);
720 setBackgroundFilter2->SetOutsideValue(3);
721 setBackgroundFilter2->Update();
722 typename InternalImageType::Pointer lungs_bones_air =setBackgroundFilter2->GetOutput();
724 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
725 caster->SetInput(lungs_bones_air);
726 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeFeature_arg, m_Verbose);
729 //----------------------------------------------------------------------------------------------------
730 // Low dimensional versions
731 //----------------------------------------------------------------------------------------------------
732 typename InternalImageType::Pointer bones_low =Resample<Dimension , PixelType>(bones);
733 typename InternalImageType::Pointer lungs_low =Resample<Dimension , PixelType>(lungs);
734 typename InternalImageType::Pointer air_low =Resample<Dimension , PixelType>(air);
736 //---------------------------------
738 //---------------------------------
739 if(m_ArgsInfo.pad_flag) {
740 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
741 IteratorType it(air_low, air_low->GetLargestPossibleRegion());
742 typename InternalImageType::IndexType index;
743 while(!it.IsAtEnd()) {
745 for (unsigned int i=0; i<Dimension; i++)
746 if(index[i]==air_low->GetLargestPossibleRegion().GetIndex()[i]
747 || index[i]==(unsigned int)air_low->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air_low->GetLargestPossibleRegion().GetSize()[i]-1)
753 //---------------------------------
755 //---------------------------------
756 typename itk::Vector<double,Dimension> center;
757 typedef itk::ImageMomentsCalculator<InternalImageType> MomentsCalculatorType;
758 typename MomentsCalculatorType::Pointer momentsCalculator=MomentsCalculatorType::New();
759 momentsCalculator->SetImage(air);
760 momentsCalculator->Compute();
761 center=momentsCalculator->GetCenterOfGravity();
763 std::cout<<"The center of gravity of the patient body is located at (mm) "<<std::endl;
764 std::cout<<center[0];
765 for (unsigned int i=1; i<Dimension; i++)
766 std::cout<<" "<<center[i];
767 std::cout<<std::endl;
770 //----------------------------------------------------------------------------------------------------
771 // Ellipsoide initialization
772 //----------------------------------------------------------------------------------------------------
773 typename InternalImageType::Pointer m_Ellips= InitializeEllips<Dimension, PixelType>(center,bones_low,lungs_low);
775 //----------------------------------------------------------------------------------------------------
777 //----------------------------------------------------------------------------------------------------
778 typename InternalImageType::Pointer grownEllips;
779 if (m_ArgsInfo.grownEllips_given || m_ArgsInfo.filledRibs_given) {
780 if (m_ArgsInfo.grownEllips_given) {
782 typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
783 featureReader->SetFileName(m_ArgsInfo.grownEllips_arg);
784 featureReader->Update();
785 grownEllips=featureReader->GetOutput();
790 std::cout<<std::endl;
791 std::cout<<"=========================================="<<std::endl;
792 std::cout<<"|| Growing ellipsoide ||"<<std::endl;
793 std::cout<<"=========================================="<<std::endl;
796 //---------------------------------
798 //---------------------------------
799 typename InternalImageType::PointType dPoint;
800 if (m_ArgsInfo.detectionPoint_given) {
801 for (unsigned int i=0; i<Dimension; i++)
802 dPoint[i]=m_ArgsInfo.detectionPoint_arg[i];
804 typename InternalImageType::RegionType searchRegion=air->GetLargestPossibleRegion();
805 typename InternalImageType::SizeType searchRegionSize = searchRegion.GetSize();
806 typename InternalImageType::IndexType searchRegionIndex = searchRegion.GetIndex();
807 searchRegionIndex[2]+=searchRegionSize[2]/2;
808 searchRegionSize[2]=1;
809 searchRegion.SetSize(searchRegionSize);
810 searchRegion.SetIndex(searchRegionIndex);
811 IteratorType itAbdomen(air, searchRegion);
812 itAbdomen.GoToBegin();
814 typename InternalImageType::PointType aPoint;
815 typename InternalImageType::IndexType aIndex;
817 if (m_Verbose) std::cout<<"Detecting the abdomen..."<<std::endl;
818 while (!itAbdomen.IsAtEnd()) {
819 if(itAbdomen.Get()) break;
822 aIndex=itAbdomen.GetIndex();
823 air->TransformIndexToPhysicalPoint(aIndex,aPoint);
824 if (m_Verbose) std::cout<<"Detected the abdomen at "<<aPoint<<".."<<std::endl;
827 //---------------------------------
828 // Detect abdomen in additional images?
829 //---------------------------------
830 if (m_ArgsInfo.detectionPairs_given) {
831 for (unsigned int i=0; i<m_ArgsInfo.detectionPairs_given; i++) {
832 typename InternalImageType::Pointer airAdd;
833 //---------------------------------
835 //--------------------------------
836 typedef itk::ImageFileReader<InputImageType> InputReaderType;
837 typename InputReaderType::Pointer reader = InputReaderType::New();
838 reader->SetFileName( m_ArgsInfo.detectionPairs_arg[i]);
840 typename InputImageType::Pointer additional= reader->GetOutput();
841 if (m_Verbose) std::cout<<"Extracting the air from additional image "<< i<<"..."<<std::endl;
843 //---------------------------------
844 // Binarize the image
845 //---------------------------------
846 typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
847 binarizeFilter->SetInput(additional);
848 binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdAir_arg));
849 binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdAir_arg));
850 if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<<m_ArgsInfo.lowerThresholdAir_arg
851 <<", "<<m_ArgsInfo.upperThresholdAir_arg<<"..."<<std::endl;
853 //---------------------------------
854 // Label the connected components
855 //---------------------------------
856 typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
857 connectFilter->SetInput(binarizeFilter->GetOutput());
858 connectFilter->SetBackgroundValue(0);
859 connectFilter->SetFullyConnected(false);
860 if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
862 //---------------------------------
863 // Sort the labels according to size
864 //---------------------------------
865 typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
866 relabelFilter->SetInput(connectFilter->GetOutput());
867 if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
869 //---------------------------------
871 //---------------------------------
872 typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
873 thresholdFilter->SetInput(relabelFilter->GetOutput());
874 thresholdFilter->SetUpper(1);
875 if (m_Verbose) std::cout<<"Keeping the first label..."<<std::endl;
876 thresholdFilter->Update();
877 airAdd=thresholdFilter->GetOutput();
880 //---------------------------------
882 //---------------------------------
883 typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
884 inversionFilter->SetInput(airAdd);
885 inversionFilter->SetLowerThreshold(0);
886 inversionFilter->SetUpperThreshold(0);
887 inversionFilter->SetInsideValue(1);
888 inversionFilter->Update();
890 //---------------------------------
892 //---------------------------------
893 typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
894 mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
895 typename InternalImageType::SizeType padSize;
897 padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
898 mirrorPadImageFilter->SetPadLowerBound(padSize);
899 if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
900 mirrorPadImageFilter->Update();
901 airAdd=mirrorPadImageFilter->GetOutput();
903 //---------------------------------
905 //---------------------------------
906 typename InternalImageType::RegionType searchRegion=airAdd->GetLargestPossibleRegion();
907 typename InternalImageType::SizeType searchRegionSize = searchRegion.GetSize();
908 typename InternalImageType::IndexType searchRegionIndex = searchRegion.GetIndex();
909 searchRegionIndex[2]+=searchRegionSize[2]/2;
910 searchRegionSize[2]=1;
911 searchRegion.SetSize(searchRegionSize);
912 searchRegion.SetIndex(searchRegionIndex);
913 IteratorType itAbdomen(airAdd, searchRegion);
914 itAbdomen.GoToBegin();
916 typename InternalImageType::PointType additionalPoint;
917 typename InternalImageType::IndexType aIndex;
919 if (m_Verbose) std::cout<<"Detecting the abdomen..."<<std::endl;
920 while (!itAbdomen.IsAtEnd()) {
921 if(itAbdomen.Get()) break;
924 aIndex=itAbdomen.GetIndex();
925 airAdd->TransformIndexToPhysicalPoint(aIndex,additionalPoint);
926 if (m_Verbose) std::cout<<"Detected the abdomen in the additional image at "<<additionalPoint<<".."<<std::endl;
928 if(additionalPoint[1]< aPoint[1]) {
929 aPoint=additionalPoint;
930 if (m_Verbose) std::cout<<"Modifying the detected abdomen to "<<aPoint<<".."<<std::endl;
937 // Determine the detection point
941 if(m_ArgsInfo.offsetDetect_given==Dimension)
942 for(unsigned int i=0; i <Dimension; i++)
943 dPoint[i]+=m_ArgsInfo.offsetDetect_arg[i];
948 if (m_Verbose) std::cout<<"Setting the detection point to "<<dPoint<<".."<<std::endl;
951 //---------------------------------
952 // Pad the rib image and ellips image
953 //---------------------------------
954 typename InternalImageType::Pointer padded_ellips;
955 typename InternalImageType::Pointer padded_bones_low;
957 // If detection point not inside the image: pad
958 typename InternalImageType::IndexType dIndex;
959 if (!bones_low->TransformPhysicalPointToIndex(dPoint, dIndex)) {
960 typename InternalImageType::SizeType padSize;
962 padSize[1]=abs(dIndex[1])+1;
963 if (m_Verbose) std::cout<<"Padding the images with padding size "<<padSize<<" to include the detection point..."<<std::endl;
965 typename MirrorPadImageFilterType::Pointer padBonesFilter=MirrorPadImageFilterType::New();
966 padBonesFilter->SetInput(bones_low);
967 padBonesFilter->SetPadLowerBound(padSize);
968 padBonesFilter->Update();
969 padded_bones_low=padBonesFilter->GetOutput();
971 typename MirrorPadImageFilterType::Pointer padEllipsFilter=MirrorPadImageFilterType::New();
972 padEllipsFilter->SetInput(m_Ellips);
973 padEllipsFilter->SetPadLowerBound(padSize);
974 padEllipsFilter->Update();
975 padded_ellips=padEllipsFilter->GetOutput();
979 padded_bones_low=bones_low;
980 padded_ellips=m_Ellips;
984 //---------------------------------
985 // Calculate distance map
986 //---------------------------------
987 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
988 distanceMapImageFilter->SetInput(padded_ellips);
989 distanceMapImageFilter->SetInsideIsPositive(false);
990 if (m_Verbose) std::cout<<"Calculating the distance map..."<<std::endl;
991 distanceMapImageFilter->Update();
992 if (m_ArgsInfo.writeDistMap_given) {
993 writeImage<LevelSetImageType>(distanceMapImageFilter->GetOutput(), m_ArgsInfo.writeDistMap_arg, m_Verbose);
997 //---------------------------------
998 // Grow while monitoring detection point
999 //---------------------------------
1000 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
1001 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
1002 levelSetFilter->SetFeatureImage( padded_bones_low );
1003 levelSetFilter->SetPropagationScaling( 1 );
1004 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve1_arg );
1005 levelSetFilter->SetAdvectionScaling( 0 );
1006 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS1_arg );
1007 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter1_arg );
1008 levelSetFilter->SetUseImageSpacing(true);
1010 // //---------------------------------
1011 // // Script for making movie
1012 // //---------------------------------
1013 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1016 if (m_Verbose) std::cout<<"Starting level set segmentation..."<<std::endl;
1017 unsigned int totalNumberOfIterations=0;
1019 levelSetFilter->Update();
1020 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1022 if ( levelSetFilter->GetOutput()->GetPixel(dIndex) < 0 ) {
1023 if (m_Verbose) std::cout<<"Detection point reached!"<<std::endl;
1026 if (m_Verbose) std::cout<<"Detection point not reached after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1027 levelSetFilter->SetInput(levelSetFilter->GetOutput());
1028 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose);
1031 // std::ostringstream number_str;
1032 // number_str << number;
1033 // std::string param = number_str.str();
1034 // system((script+param).c_str());
1038 if ( (totalNumberOfIterations> 10000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS1_arg) ) break;
1042 if (m_Verbose) std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1043 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1044 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1047 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1048 thresholder->SetUpperThreshold( 0.0 );
1049 thresholder->SetOutsideValue( 0 );
1050 thresholder->SetInsideValue( 1 );
1051 thresholder->SetInput( levelSetFilter->GetOutput() );
1052 if (m_Verbose) std::cout<<"Thresholding the output level set..."<<std::endl;
1053 thresholder->Update();
1054 grownEllips=thresholder->GetOutput();
1057 //---------------------------------
1058 // Write the grown ellips
1059 //---------------------------------
1060 if (m_ArgsInfo.writeGrownEllips_given) {
1061 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1062 caster->SetInput(grownEllips);
1063 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeGrownEllips_arg, m_Verbose);
1067 //----------------------------------------------------------------------------------------------------
1069 //----------------------------------------------------------------------------------------------------
1070 typename InternalImageType::Pointer filledRibs;
1071 if (m_ArgsInfo.filledRibs_given) {
1072 typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
1073 featureReader->SetFileName(m_ArgsInfo.filledRibs_arg);
1074 if (m_Verbose) std::cout<<"Reading filled ribs mask..."<<std::endl;
1075 featureReader->Update();
1076 filledRibs=featureReader->GetOutput();
1079 std::cout<<std::endl;
1080 std::cout<<"=========================================="<<std::endl;
1081 std::cout<<"|| Filling the ribs image ||"<<std::endl;
1082 std::cout<<"=========================================="<<std::endl;
1084 //---------------------------------
1085 // Make feature image air+bones
1086 //---------------------------------
1087 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
1088 setBackgroundFilter->SetInput(air_low);
1089 setBackgroundFilter->SetInput2(bones_low);
1090 setBackgroundFilter->SetMaskValue(0);
1091 setBackgroundFilter->SetOutsideValue(0);
1092 setBackgroundFilter->Update();
1093 typename InternalImageType::Pointer bones_air_low =setBackgroundFilter->GetOutput();
1095 //---------------------------------
1096 // Resampling previous solution
1097 //---------------------------------
1098 if (m_Verbose) std::cout<<"Resampling previous solution..."<<std::endl;
1099 typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
1100 typename InternalImageType::PointType origin;
1101 bones_low->TransformIndexToPhysicalPoint(bones_low->GetLargestPossibleRegion().GetIndex(), origin);
1102 resampler->SetOutputOrigin(origin);
1103 resampler->SetOutputSpacing(bones_low->GetSpacing());
1104 typename InternalImageType::SizeType size_low= bones_low->GetLargestPossibleRegion().GetSize();
1105 resampler->SetSize(size_low);
1106 typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
1107 resampler->SetInterpolator(interpolator);
1108 resampler->SetInput(grownEllips);
1109 resampler->Update();
1110 typename InternalImageType::Pointer grownEllips =resampler->GetOutput();
1113 //---------------------------------
1114 // Calculate Distance Map
1115 //---------------------------------
1116 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
1117 distanceMapImageFilter->SetInput(grownEllips);
1118 distanceMapImageFilter->SetInsideIsPositive(false);
1119 if (m_Verbose) std::cout<<"Calculating distance map..."<<std::endl;
1120 distanceMapImageFilter->Update();
1121 if (m_ArgsInfo.writeDistMap_given) {
1122 writeImage<LevelSetImageType>(distanceMapImageFilter->GetOutput(), m_ArgsInfo.writeDistMap_arg, m_Verbose);
1126 //---------------------------------
1127 // Grow while monitoring lung volume coverage
1128 //---------------------------------
1129 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
1130 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
1131 levelSetFilter->SetFeatureImage( bones_air_low );
1132 levelSetFilter->SetPropagationScaling( 1 );
1133 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve2_arg );
1134 levelSetFilter->SetAdvectionScaling( 0 );
1135 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS2_arg );
1136 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter2_arg );
1137 levelSetFilter->SetUseImageSpacing(true);
1139 //---------------------------------
1141 //---------------------------------
1142 typename InversionFilterType::Pointer invertor= InversionFilterType::New();
1143 invertor->SetInput(lungs_low);
1144 invertor->SetLowerThreshold(0);
1145 invertor->SetUpperThreshold(0);
1146 invertor->SetInsideValue(1);
1148 typename InternalImageType::Pointer lungs_low_inv=invertor->GetOutput();
1150 //---------------------------------
1151 // Calculate lung volume
1152 //---------------------------------
1153 typename LabelStatisticsImageFilterType::Pointer statisticsFilter= LabelStatisticsImageFilterType::New();
1154 statisticsFilter->SetInput(lungs_low_inv);
1155 statisticsFilter->SetLabelInput(lungs_low);
1156 statisticsFilter->Update();
1157 unsigned int volume=statisticsFilter->GetSum(0);
1159 // Prepare threshold
1160 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1161 thresholder->SetUpperThreshold( 0.0 );
1162 thresholder->SetOutsideValue( 0 );
1163 thresholder->SetInsideValue( 1 );
1167 unsigned int totalNumberOfIterations=0;
1168 unsigned int coverage=0;
1171 // //---------------------------------
1172 // // Script for making movie
1173 // //---------------------------------
1174 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1177 levelSetFilter->Update();
1178 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1180 thresholder->SetInput( levelSetFilter->GetOutput() );
1181 thresholder->Update();
1182 statisticsFilter->SetInput(thresholder->GetOutput());
1183 statisticsFilter->Update();
1184 coverage=statisticsFilter->GetSum(0);
1186 // Compare the volumes
1187 if ( coverage >= m_ArgsInfo.fillingLevel_arg * volume/100 ) {
1188 if (m_Verbose) std::cout<<"Lungs filled for "<< m_ArgsInfo.fillingLevel_arg<<"% !"<<std::endl;
1191 if (m_Verbose) std::cout<<"After "<<totalNumberOfIterations<<" iterations, lungs are covered for "
1192 <<(double)coverage/volume*100.0<<"%..."<<std::endl;
1193 levelSetFilter->SetInput(levelSetFilter->GetOutput());
1194 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose);
1197 // std::ostringstream number_str;
1198 // number_str << number;
1199 // std::string param = number_str.str();
1200 // system((script+param).c_str());
1204 if ( (totalNumberOfIterations> 30000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS2_arg) ) break;
1208 if (m_Verbose) std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1209 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1210 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1213 thresholder->SetInput( levelSetFilter->GetOutput() );
1214 thresholder->Update();
1215 filledRibs=thresholder->GetOutput();
1216 // writeImage<InternalImageType>(filledRibs,"/home/jef/tmp/filled_ribs_before.mhd", m_Verbose);
1220 //---------------------------------
1221 // Write the filled ribs
1222 //---------------------------------
1223 if (m_ArgsInfo.writeFilledRibs_given) {
1224 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1225 caster->SetInput(filledRibs);
1226 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeFilledRibs_arg, m_Verbose);
1230 //----------------------------------------------------------------------------------------------------
1231 // Collapse to the lungs
1232 //----------------------------------------------------------------------------------------------------
1234 std::cout<<std::endl;
1235 std::cout<<"=========================================="<<std::endl;
1236 std::cout<<"|| Collapsing to the lung image ||"<<std::endl;
1237 std::cout<<"=========================================="<<std::endl;
1240 //---------------------------------
1241 // Make feature image air+bones
1242 //---------------------------------
1243 if (m_Verbose) std::cout<<"Making feature images..."<<std::endl;
1244 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
1245 setBackgroundFilter->SetInput(air);
1246 setBackgroundFilter->SetInput2(bones);
1247 setBackgroundFilter->SetMaskValue(0);
1248 setBackgroundFilter->SetOutsideValue(0);
1249 setBackgroundFilter->Update();
1250 typename InternalImageType::Pointer bones_air =setBackgroundFilter->GetOutput();
1252 typename SetBackgroundFilterType::Pointer setBackgroundFilter2 = SetBackgroundFilterType::New();
1253 setBackgroundFilter2->SetInput(bones_air);
1254 setBackgroundFilter2->SetInput2(lungs);
1255 setBackgroundFilter2->SetMaskValue(0);
1256 setBackgroundFilter2->SetOutsideValue(0);
1257 setBackgroundFilter2->Update();
1258 typename InternalImageType::Pointer lungs_bones_air =setBackgroundFilter2->GetOutput();
1260 //---------------------------------
1261 // Prepare previous solution
1262 //---------------------------------
1263 if (m_Verbose) std::cout<<"Resampling previous solution..."<<std::endl;
1264 typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
1265 typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
1266 typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
1267 resampler->SetOutputStartIndex(bones->GetLargestPossibleRegion().GetIndex());
1268 resampler->SetInput(filledRibs);
1269 resampler->SetInterpolator(interpolator);
1270 resampler->SetOutputParametersFromImage(bones);
1271 resampler->Update();
1272 filledRibs =resampler->GetOutput();
1274 if (m_Verbose) std::cout<<"Setting lungs to 1..."<<std::endl;
1275 typename SetBackgroundFilterType::Pointer setBackgroundFilter3 = SetBackgroundFilterType::New();
1276 setBackgroundFilter3->SetInput(filledRibs);
1277 setBackgroundFilter3->SetInput2(lungs);
1278 setBackgroundFilter3->SetMaskValue(0);
1279 setBackgroundFilter3->SetOutsideValue(1);
1280 setBackgroundFilter3->Update();
1281 filledRibs=setBackgroundFilter3->GetOutput();
1283 if (m_Verbose) std::cout<<"Removing overlap with bones..."<<std::endl;
1284 typename SetBackgroundFilterType::Pointer setBackgroundFilter4 = SetBackgroundFilterType::New();
1285 setBackgroundFilter4->SetInput(filledRibs);
1286 setBackgroundFilter4->SetInput2(bones);
1287 setBackgroundFilter4->SetMaskValue(0);
1288 setBackgroundFilter4->SetOutsideValue(0);
1289 setBackgroundFilter4->Update();
1290 filledRibs =setBackgroundFilter4->GetOutput();
1291 // writeImage<InternalImageType>(filledRibs,"/home/jef/tmp/filledRibs_pp.mhd");
1292 //---------------------------------
1293 // Calculate Distance Map
1294 //---------------------------------
1295 // typedef itk::ApproximateSignedDistanceMapImageFilter <InternalImageType, LevelSetImageType> DistanceMapImageFilterType2;
1296 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
1297 distanceMapImageFilter->SetInput(filledRibs);
1298 distanceMapImageFilter->SetInsideIsPositive(false);
1299 // distanceMapImageFilter->SetInsideValue(0);
1300 // distanceMapImageFilter->SetOutsideValue(1);
1301 if (m_Verbose) std::cout<<"Calculating distance map..."<<std::endl;
1302 distanceMapImageFilter->Update();
1304 //---------------------------------
1306 //---------------------------------
1307 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
1308 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
1309 levelSetFilter->SetFeatureImage( lungs_bones_air );
1310 levelSetFilter->SetPropagationScaling(m_ArgsInfo.prop3_arg);
1311 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve3_arg );
1312 levelSetFilter->SetAdvectionScaling( 0 );
1313 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS3_arg );
1314 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter3_arg );
1315 levelSetFilter->SetUseImageSpacing(true);
1318 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1321 if (m_Verbose) std::cout<<"Starting the levelset..."<<std::endl;
1322 unsigned int totalNumberOfIterations=0;
1324 levelSetFilter->Update();
1327 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1328 levelSetFilter->SetInput(levelSetFilter->GetOutput());
1329 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg);
1330 std::cout << "After "<<totalNumberOfIterations<<" iteration the RMS change is " << levelSetFilter->GetRMSChange() <<"..."<< std::endl;
1333 // std::ostringstream number_str;
1334 // number_str << number;
1335 // std::string param = number_str.str();
1336 // system((script+param).c_str());
1339 // stopping condition
1340 if ( (totalNumberOfIterations>= (unsigned int) m_ArgsInfo.maxIter3_arg) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS3_arg) ) break;
1341 levelSetFilter->SetNumberOfIterations( std::min ((unsigned int) m_ArgsInfo.iter3_arg, (unsigned int) m_ArgsInfo.maxIter3_arg-totalNumberOfIterations ) );
1346 std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1347 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1348 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1352 typedef itk::BinaryThresholdImageFilter< LevelSetImageType,InternalImageType> LevelSetBinarizeFilterType;
1353 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1354 thresholder->SetUpperThreshold( 0.0 );
1355 thresholder->SetOutsideValue( 0 );
1356 thresholder->SetInsideValue( 1 );
1357 thresholder->SetInput( levelSetFilter->GetOutput() );
1358 thresholder->Update();
1359 typename InternalImageType::Pointer output = thresholder->GetOutput();
1362 //----------------------------------------------------------------------------------------------------
1363 // Prepare the output
1364 //----------------------------------------------------------------------------------------------------
1366 //---------------------------------
1368 //---------------------------------
1369 if (m_Verbose) std::cout<<"Removing overlap mask with air..."<<std::endl;
1370 typename SetBackgroundFilterType::Pointer setBackgroundFilter5 = SetBackgroundFilterType::New();
1371 setBackgroundFilter5->SetInput(output);
1372 setBackgroundFilter5->SetInput2(air);
1373 setBackgroundFilter5->SetMaskValue(0);
1374 setBackgroundFilter5->SetOutsideValue(0);
1375 setBackgroundFilter5->Update();
1376 output=setBackgroundFilter5->GetOutput();
1378 //---------------------------------
1379 // Open & close to cleanup
1380 //---------------------------------
1381 if ( m_ArgsInfo.openClose_flag) {
1383 //---------------------------------
1384 // Structuring element
1385 //---------------------------------
1386 typedef itk::BinaryBallStructuringElement<InternalPixelType,Dimension > KernelType;
1387 KernelType structuringElement;
1388 structuringElement.SetRadius(1);
1389 structuringElement.CreateStructuringElement();
1391 //---------------------------------
1393 //---------------------------------
1394 typedef itk::BinaryMorphologicalOpeningImageFilter<InternalImageType, InternalImageType , KernelType> OpenFilterType;
1395 typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
1396 openFilter->SetInput(output);
1397 openFilter->SetBackgroundValue(0);
1398 openFilter->SetForegroundValue(1);
1399 openFilter->SetKernel(structuringElement);
1400 if(m_Verbose) std::cout<<"Opening the output image..."<<std::endl;
1402 //---------------------------------
1404 //---------------------------------
1405 typedef itk::BinaryMorphologicalClosingImageFilter<InternalImageType, InternalImageType , KernelType> CloseFilterType;
1406 typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
1407 closeFilter->SetInput(openFilter->GetOutput());
1408 closeFilter->SetSafeBorder(true);
1409 closeFilter->SetForegroundValue(1);
1410 closeFilter->SetKernel(structuringElement);
1411 if(m_Verbose) std::cout<<"Closing the output image..."<<std::endl;
1412 closeFilter->Update();
1413 output=closeFilter->GetOutput();
1416 //writeImage<InternalImageType>(output,"/home/jef/tmp/mm_double.mhd");
1418 // Extract the upper part
1419 typedef itk::CropImageFilter<InternalImageType, InternalImageType> CropImageFilterType;
1420 typename CropImageFilterType::Pointer cropFilter=CropImageFilterType::New();
1421 cropFilter->SetInput(output);
1422 typename InputImageType::SizeType lSize, uSize;
1425 lSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
1426 cropFilter->SetLowerBoundaryCropSize(lSize);
1427 cropFilter->SetUpperBoundaryCropSize(uSize);
1428 cropFilter->Update();
1431 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1432 caster->SetInput(cropFilter->GetOutput());
1433 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.output_arg, m_Verbose);
1439 #endif //#define clitkMotionMaskGenericFilter_txx