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 typename LabelType::RegionType lung_bbox = label->GetRegion();
536 axes[0] = 0.95*lung_bbox.GetSize()[0]*spacing[0]/2;
537 axes[1] = 0.3*lung_bbox.GetSize()[1]*spacing[1]/2;
538 axes[2] = 0.95*lung_bbox.GetSize()[2]*spacing[2]/2;
540 // ellipse's center is the center of the lungs' bounding box
541 typename InternalImageType::PointType origin = lungs_low->GetOrigin();
542 centerEllips[0] = origin[0] + (lung_bbox.GetIndex()[0] + lung_bbox.GetSize()[0]/2)*spacing[0];
543 centerEllips[1] = origin[1] + (lung_bbox.GetIndex()[1] + lung_bbox.GetSize()[1]/2)*spacing[1];
544 centerEllips[2] = origin[2] + (lung_bbox.GetIndex()[2] + lung_bbox.GetSize()[2]/2)*spacing[2];
547 std::cout << "Ellipse center at " << centerEllips << std::endl;
548 std::cout << "Ellipse axes as " << axes << std::endl;
549 std::cout << "Lung's bounding box (index,size) " << lung_bbox.GetIndex() << lung_bbox.GetSize() << std::endl;
553 //---------------------------------
554 // Offset from center
555 //---------------------------------
556 typename itk::Vector<double,Dimension> offset;
557 if(m_ArgsInfo.offset_given== Dimension) {
558 for(unsigned int i=0; i<Dimension; i++)
559 offset[i]=m_ArgsInfo.offset_arg[i];
564 centerEllips=center+offset;
566 std::cout<<"Placing the center of the initial ellipsoide at (mm) "<<std::endl;
567 std::cout<<centerEllips[0];
568 for (unsigned int i=1; i<Dimension; i++)
569 std::cout<<" "<<centerEllips[i];
570 std::cout<<std::endl;
574 if (m_ArgsInfo.axes_given==Dimension)
575 for(unsigned int i=0; i<Dimension; i++)
576 axes[i]=m_ArgsInfo.axes_arg[i];
584 //---------------------------------
586 //---------------------------------
587 ellips=InternalImageType::New();
588 ellips->SetRegions (bones_low->GetLargestPossibleRegion());
589 ellips->SetOrigin(bones_low->GetOrigin());
590 ellips->SetSpacing(bones_low->GetSpacing());
592 ellips->FillBuffer(0);
595 IteratorType itEllips(ellips, ellips->GetLargestPossibleRegion());
596 itEllips.GoToBegin();
597 typename InternalImageType::PointType point;
598 typename InternalImageType::IndexType index;
601 if (m_Verbose) std::cout<<"Drawing the initial ellipsoide..."<<std::endl;
602 while (!itEllips.IsAtEnd()) {
603 index=itEllips.GetIndex();
604 ellips->TransformIndexToPhysicalPoint(index, point);
606 for(unsigned int i=0; i<Dimension; i++)
607 distance+=powf( ( (centerEllips[i]-point[i])/axes[i] ), 2);
616 //---------------------------------
618 //---------------------------------
619 if (m_ArgsInfo.writeEllips_given) {
620 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
621 caster->SetInput(ellips);
622 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeEllips_arg, m_Verbose);
630 //-------------------------------------------------------------------
631 // Update with the number of dimensions and the pixeltype
632 //-------------------------------------------------------------------
633 template <unsigned int Dimension, class PixelType>
635 MotionMaskGenericFilter::UpdateWithDimAndPixelType()
639 typedef int InternalPixelType;
640 typedef itk::Image<PixelType, Dimension> InputImageType;
641 typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
642 typedef itk::Image< float, Dimension> LevelSetImageType; //labeling doesn't work with unsigned char?
643 typedef itk::Image<unsigned char, Dimension> OutputImageType;
646 //----------------------------------------------------------------------------------------------------
647 // Typedef for Processing
648 //----------------------------------------------------------------------------------------------------
649 typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
650 typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
651 typedef itk::BinaryThresholdImageFilter< LevelSetImageType,InternalImageType> LevelSetBinarizeFilterType;
652 typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> InversionFilterType;
653 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdFilterType;
654 typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
655 typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
656 typedef itk::MirrorPadImageFilter<InternalImageType, InternalImageType> MirrorPadImageFilterType;
657 typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
658 typedef itk::ResampleImageFilter<InternalImageType, InternalImageType> ResampleImageFilterType;
659 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
660 typedef itk::GeodesicActiveContourLevelSetImageFilter<LevelSetImageType, InternalImageType> LevelSetImageFilterType;
661 typedef itk::SignedDanielssonDistanceMapImageFilter<InternalImageType, LevelSetImageType> DistanceMapImageFilterType;
662 typedef clitk::SetBackgroundImageFilter<InternalImageType,InternalImageType, InternalImageType> SetBackgroundFilterType;
663 typedef itk::LabelStatisticsImageFilter<InternalImageType,InternalImageType> LabelStatisticsImageFilterType;
664 typedef itk::CastImageFilter<InternalImageType,OutputImageType> CastImageFilterType;
668 typedef itk::ImageFileReader<InputImageType> InputReaderType;
669 typename InputReaderType::Pointer reader = InputReaderType::New();
670 reader->SetFileName( m_InputFileName);
672 typename InputImageType::Pointer input= reader->GetOutput();
674 // // globals for avi
675 // unsigned int number=0;
679 std::cout<<std::endl;
680 std::cout<<"=========================================="<<std::endl;
681 std::cout<<"|| Making feature images ||"<<std::endl;
682 std::cout<<"=========================================="<<std::endl;
685 //--------------------------------------------------------------------------------
687 //-------------------------------------------------------------------------------
688 typename InternalImageType::Pointer lungs = this->GetLungsImage<Dimension, PixelType>(input);
690 //-------------------------------------------------------------------------------
692 //-------------------------------------------------------------------------------
693 typename InternalImageType::Pointer air = this->GetAirImage<Dimension, PixelType>(input, lungs);
695 //-------------------------------------------------------------------------------
697 //-------------------------------------------------------------------------------
698 typename InternalImageType::Pointer bones = this->GetBonesImage<Dimension, PixelType>(input);
700 //----------------------------------------------------------------------------------------------------
701 // Write feature images
702 //----------------------------------------------------------------------------------------------------
703 if(m_ArgsInfo.writeFeature_given) {
704 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
705 setBackgroundFilter->SetInput(air);
706 setBackgroundFilter->SetInput2(bones);
707 setBackgroundFilter->SetMaskValue(0);
708 setBackgroundFilter->SetOutsideValue(2);
709 setBackgroundFilter->Update();
710 typename InternalImageType::Pointer bones_air =setBackgroundFilter->GetOutput();
712 typename SetBackgroundFilterType::Pointer setBackgroundFilter2 = SetBackgroundFilterType::New();
713 setBackgroundFilter2->SetInput(bones_air);
714 setBackgroundFilter2->SetInput2(lungs);
715 setBackgroundFilter2->SetMaskValue(0);
716 setBackgroundFilter2->SetOutsideValue(3);
717 setBackgroundFilter2->Update();
718 typename InternalImageType::Pointer lungs_bones_air =setBackgroundFilter2->GetOutput();
720 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
721 caster->SetInput(lungs_bones_air);
722 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeFeature_arg, m_Verbose);
725 //----------------------------------------------------------------------------------------------------
726 // Low dimensional versions
727 //----------------------------------------------------------------------------------------------------
728 typename InternalImageType::Pointer bones_low =Resample<Dimension , PixelType>(bones);
729 typename InternalImageType::Pointer lungs_low =Resample<Dimension , PixelType>(lungs);
730 typename InternalImageType::Pointer air_low =Resample<Dimension , PixelType>(air);
732 //---------------------------------
734 //---------------------------------
735 if(m_ArgsInfo.pad_flag) {
736 typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
737 IteratorType it(air_low, air_low->GetLargestPossibleRegion());
738 typename InternalImageType::IndexType index;
739 while(!it.IsAtEnd()) {
741 for (unsigned int i=0; i<Dimension; i++)
742 if(index[i]==air_low->GetLargestPossibleRegion().GetIndex()[i]
743 || index[i]==(unsigned int)air_low->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air_low->GetLargestPossibleRegion().GetSize()[i]-1)
749 //---------------------------------
751 //---------------------------------
752 typename itk::Vector<double,Dimension> center;
753 typedef itk::ImageMomentsCalculator<InternalImageType> MomentsCalculatorType;
754 typename MomentsCalculatorType::Pointer momentsCalculator=MomentsCalculatorType::New();
755 momentsCalculator->SetImage(air);
756 momentsCalculator->Compute();
757 center=momentsCalculator->GetCenterOfGravity();
759 std::cout<<"The center of gravity of the patient body is located at (mm) "<<std::endl;
760 std::cout<<center[0];
761 for (unsigned int i=1; i<Dimension; i++)
762 std::cout<<" "<<center[i];
763 std::cout<<std::endl;
766 //----------------------------------------------------------------------------------------------------
767 // Ellipsoide initialization
768 //----------------------------------------------------------------------------------------------------
769 typename InternalImageType::Pointer m_Ellips= InitializeEllips<Dimension, PixelType>(center,bones_low,lungs_low);
771 //----------------------------------------------------------------------------------------------------
773 //----------------------------------------------------------------------------------------------------
774 typename InternalImageType::Pointer grownEllips;
775 if (m_ArgsInfo.grownEllips_given || m_ArgsInfo.filledRibs_given) {
776 if (m_ArgsInfo.grownEllips_given) {
778 typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
779 featureReader->SetFileName(m_ArgsInfo.grownEllips_arg);
780 featureReader->Update();
781 grownEllips=featureReader->GetOutput();
786 std::cout<<std::endl;
787 std::cout<<"=========================================="<<std::endl;
788 std::cout<<"|| Growing ellipsoide ||"<<std::endl;
789 std::cout<<"=========================================="<<std::endl;
792 //---------------------------------
794 //---------------------------------
795 typename InternalImageType::PointType dPoint;
796 if (m_ArgsInfo.detectionPoint_given) {
797 for (unsigned int i=0; i<Dimension; i++)
798 dPoint[i]=m_ArgsInfo.detectionPoint_arg[i];
800 typename InternalImageType::RegionType searchRegion=air->GetLargestPossibleRegion();
801 typename InternalImageType::SizeType searchRegionSize = searchRegion.GetSize();
802 typename InternalImageType::IndexType searchRegionIndex = searchRegion.GetIndex();
803 searchRegionIndex[2]+=searchRegionSize[2]/2;
804 searchRegionSize[2]=1;
805 searchRegion.SetSize(searchRegionSize);
806 searchRegion.SetIndex(searchRegionIndex);
807 IteratorType itAbdomen(air, searchRegion);
808 itAbdomen.GoToBegin();
810 typename InternalImageType::PointType aPoint;
811 typename InternalImageType::IndexType aIndex;
813 if (m_Verbose) std::cout<<"Detecting the abdomen..."<<std::endl;
814 while (!itAbdomen.IsAtEnd()) {
815 if(itAbdomen.Get()) break;
818 aIndex=itAbdomen.GetIndex();
819 air->TransformIndexToPhysicalPoint(aIndex,aPoint);
820 if (m_Verbose) std::cout<<"Detected the abdomen at "<<aPoint<<".."<<std::endl;
823 //---------------------------------
824 // Detect abdomen in additional images?
825 //---------------------------------
826 if (m_ArgsInfo.detectionPairs_given) {
827 for (unsigned int i=0; i<m_ArgsInfo.detectionPairs_given; i++) {
828 typename InternalImageType::Pointer airAdd;
829 //---------------------------------
831 //--------------------------------
832 typedef itk::ImageFileReader<InputImageType> InputReaderType;
833 typename InputReaderType::Pointer reader = InputReaderType::New();
834 reader->SetFileName( m_ArgsInfo.detectionPairs_arg[i]);
836 typename InputImageType::Pointer additional= reader->GetOutput();
837 if (m_Verbose) std::cout<<"Extracting the air from additional image "<< i<<"..."<<std::endl;
839 //---------------------------------
840 // Binarize the image
841 //---------------------------------
842 typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
843 binarizeFilter->SetInput(additional);
844 binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdAir_arg));
845 binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdAir_arg));
846 if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<<m_ArgsInfo.lowerThresholdAir_arg
847 <<", "<<m_ArgsInfo.upperThresholdAir_arg<<"..."<<std::endl;
849 //---------------------------------
850 // Label the connected components
851 //---------------------------------
852 typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
853 connectFilter->SetInput(binarizeFilter->GetOutput());
854 connectFilter->SetBackgroundValue(0);
855 connectFilter->SetFullyConnected(false);
856 if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
858 //---------------------------------
859 // Sort the labels according to size
860 //---------------------------------
861 typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
862 relabelFilter->SetInput(connectFilter->GetOutput());
863 if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
865 //---------------------------------
867 //---------------------------------
868 typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
869 thresholdFilter->SetInput(relabelFilter->GetOutput());
870 thresholdFilter->SetUpper(1);
871 if (m_Verbose) std::cout<<"Keeping the first label..."<<std::endl;
872 thresholdFilter->Update();
873 airAdd=thresholdFilter->GetOutput();
876 //---------------------------------
878 //---------------------------------
879 typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
880 inversionFilter->SetInput(airAdd);
881 inversionFilter->SetLowerThreshold(0);
882 inversionFilter->SetUpperThreshold(0);
883 inversionFilter->SetInsideValue(1);
884 inversionFilter->Update();
886 //---------------------------------
888 //---------------------------------
889 typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
890 mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
891 typename InternalImageType::SizeType padSize;
893 padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
894 mirrorPadImageFilter->SetPadLowerBound(padSize);
895 if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
896 mirrorPadImageFilter->Update();
897 airAdd=mirrorPadImageFilter->GetOutput();
899 //---------------------------------
901 //---------------------------------
902 typename InternalImageType::RegionType searchRegion=airAdd->GetLargestPossibleRegion();
903 typename InternalImageType::SizeType searchRegionSize = searchRegion.GetSize();
904 typename InternalImageType::IndexType searchRegionIndex = searchRegion.GetIndex();
905 searchRegionIndex[2]+=searchRegionSize[2]/2;
906 searchRegionSize[2]=1;
907 searchRegion.SetSize(searchRegionSize);
908 searchRegion.SetIndex(searchRegionIndex);
909 IteratorType itAbdomen(airAdd, searchRegion);
910 itAbdomen.GoToBegin();
912 typename InternalImageType::PointType additionalPoint;
913 typename InternalImageType::IndexType aIndex;
915 if (m_Verbose) std::cout<<"Detecting the abdomen..."<<std::endl;
916 while (!itAbdomen.IsAtEnd()) {
917 if(itAbdomen.Get()) break;
920 aIndex=itAbdomen.GetIndex();
921 airAdd->TransformIndexToPhysicalPoint(aIndex,additionalPoint);
922 if (m_Verbose) std::cout<<"Detected the abdomen in the additional image at "<<additionalPoint<<".."<<std::endl;
924 if(additionalPoint[1]< aPoint[1]) {
925 aPoint=additionalPoint;
926 if (m_Verbose) std::cout<<"Modifying the detected abdomen to "<<aPoint<<".."<<std::endl;
933 // Determine the detection point
937 if(m_ArgsInfo.offsetDetect_given==Dimension)
938 for(unsigned int i=0; i <Dimension; i++)
939 dPoint[i]+=m_ArgsInfo.offsetDetect_arg[i];
944 if (m_Verbose) std::cout<<"Setting the detection point to "<<dPoint<<".."<<std::endl;
947 //---------------------------------
948 // Pad the rib image and ellips image
949 //---------------------------------
950 typename InternalImageType::Pointer padded_ellips;
951 typename InternalImageType::Pointer padded_bones_low;
953 // If detection point not inside the image: pad
954 typename InternalImageType::IndexType dIndex;
955 if (!bones_low->TransformPhysicalPointToIndex(dPoint, dIndex)) {
956 typename InternalImageType::SizeType padSize;
958 padSize[1]=abs(dIndex[1])+1;
959 if (m_Verbose) std::cout<<"Padding the images with padding size "<<padSize<<" to include the detection point..."<<std::endl;
961 typename MirrorPadImageFilterType::Pointer padBonesFilter=MirrorPadImageFilterType::New();
962 padBonesFilter->SetInput(bones_low);
963 padBonesFilter->SetPadLowerBound(padSize);
964 padBonesFilter->Update();
965 padded_bones_low=padBonesFilter->GetOutput();
967 typename MirrorPadImageFilterType::Pointer padEllipsFilter=MirrorPadImageFilterType::New();
968 padEllipsFilter->SetInput(m_Ellips);
969 padEllipsFilter->SetPadLowerBound(padSize);
970 padEllipsFilter->Update();
971 padded_ellips=padEllipsFilter->GetOutput();
975 padded_bones_low=bones_low;
976 padded_ellips=m_Ellips;
980 //---------------------------------
981 // Calculate distance map
982 //---------------------------------
983 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
984 distanceMapImageFilter->SetInput(padded_ellips);
985 distanceMapImageFilter->SetInsideIsPositive(false);
986 if (m_Verbose) std::cout<<"Calculating the distance map..."<<std::endl;
987 distanceMapImageFilter->Update();
988 if (m_ArgsInfo.writeDistMap_given) {
989 writeImage<LevelSetImageType>(distanceMapImageFilter->GetOutput(), m_ArgsInfo.writeDistMap_arg, m_Verbose);
993 //---------------------------------
994 // Grow while monitoring detection point
995 //---------------------------------
996 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
997 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
998 levelSetFilter->SetFeatureImage( padded_bones_low );
999 levelSetFilter->SetPropagationScaling( 1 );
1000 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve1_arg );
1001 levelSetFilter->SetAdvectionScaling( 0 );
1002 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS1_arg );
1003 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter1_arg );
1004 levelSetFilter->SetUseImageSpacing(true);
1006 // //---------------------------------
1007 // // Script for making movie
1008 // //---------------------------------
1009 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1012 if (m_Verbose) std::cout<<"Starting level set segmentation..."<<std::endl;
1013 unsigned int totalNumberOfIterations=0;
1015 levelSetFilter->Update();
1016 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1018 if ( levelSetFilter->GetOutput()->GetPixel(dIndex) < 0 ) {
1019 if (m_Verbose) std::cout<<"Detection point reached!"<<std::endl;
1022 if (m_Verbose) std::cout<<"Detection point not reached after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1023 levelSetFilter->SetInput(levelSetFilter->GetOutput());
1024 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose);
1027 // std::ostringstream number_str;
1028 // number_str << number;
1029 // std::string param = number_str.str();
1030 // system((script+param).c_str());
1034 if ( (totalNumberOfIterations> 10000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS1_arg) ) break;
1038 if (m_Verbose) std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1039 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1040 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1043 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1044 thresholder->SetUpperThreshold( 0.0 );
1045 thresholder->SetOutsideValue( 0 );
1046 thresholder->SetInsideValue( 1 );
1047 thresholder->SetInput( levelSetFilter->GetOutput() );
1048 if (m_Verbose) std::cout<<"Thresholding the output level set..."<<std::endl;
1049 thresholder->Update();
1050 grownEllips=thresholder->GetOutput();
1053 //---------------------------------
1054 // Write the grown ellips
1055 //---------------------------------
1056 if (m_ArgsInfo.writeGrownEllips_given) {
1057 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1058 caster->SetInput(grownEllips);
1059 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeGrownEllips_arg, m_Verbose);
1063 //----------------------------------------------------------------------------------------------------
1065 //----------------------------------------------------------------------------------------------------
1066 typename InternalImageType::Pointer filledRibs;
1067 if (m_ArgsInfo.filledRibs_given) {
1068 typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
1069 featureReader->SetFileName(m_ArgsInfo.filledRibs_arg);
1070 if (m_Verbose) std::cout<<"Reading filled ribs mask..."<<std::endl;
1071 featureReader->Update();
1072 filledRibs=featureReader->GetOutput();
1075 std::cout<<std::endl;
1076 std::cout<<"=========================================="<<std::endl;
1077 std::cout<<"|| Filling the ribs image ||"<<std::endl;
1078 std::cout<<"=========================================="<<std::endl;
1080 //---------------------------------
1081 // Make feature image air+bones
1082 //---------------------------------
1083 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
1084 setBackgroundFilter->SetInput(air_low);
1085 setBackgroundFilter->SetInput2(bones_low);
1086 setBackgroundFilter->SetMaskValue(0);
1087 setBackgroundFilter->SetOutsideValue(0);
1088 setBackgroundFilter->Update();
1089 typename InternalImageType::Pointer bones_air_low =setBackgroundFilter->GetOutput();
1091 //---------------------------------
1092 // Resampling previous solution
1093 //---------------------------------
1094 if (m_Verbose) std::cout<<"Resampling previous solution..."<<std::endl;
1095 typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
1096 typename InternalImageType::PointType origin;
1097 bones_low->TransformIndexToPhysicalPoint(bones_low->GetLargestPossibleRegion().GetIndex(), origin);
1098 resampler->SetOutputOrigin(origin);
1099 resampler->SetOutputSpacing(bones_low->GetSpacing());
1100 typename InternalImageType::SizeType size_low= bones_low->GetLargestPossibleRegion().GetSize();
1101 resampler->SetSize(size_low);
1102 typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
1103 resampler->SetInterpolator(interpolator);
1104 resampler->SetInput(grownEllips);
1105 resampler->Update();
1106 typename InternalImageType::Pointer grownEllips =resampler->GetOutput();
1109 //---------------------------------
1110 // Calculate Distance Map
1111 //---------------------------------
1112 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
1113 distanceMapImageFilter->SetInput(grownEllips);
1114 distanceMapImageFilter->SetInsideIsPositive(false);
1115 if (m_Verbose) std::cout<<"Calculating distance map..."<<std::endl;
1116 distanceMapImageFilter->Update();
1117 if (m_ArgsInfo.writeDistMap_given) {
1118 writeImage<LevelSetImageType>(distanceMapImageFilter->GetOutput(), m_ArgsInfo.writeDistMap_arg, m_Verbose);
1122 //---------------------------------
1123 // Grow while monitoring lung volume coverage
1124 //---------------------------------
1125 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
1126 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
1127 levelSetFilter->SetFeatureImage( bones_air_low );
1128 levelSetFilter->SetPropagationScaling( 1 );
1129 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve2_arg );
1130 levelSetFilter->SetAdvectionScaling( 0 );
1131 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS2_arg );
1132 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter2_arg );
1133 levelSetFilter->SetUseImageSpacing(true);
1135 //---------------------------------
1137 //---------------------------------
1138 typename InversionFilterType::Pointer invertor= InversionFilterType::New();
1139 invertor->SetInput(lungs_low);
1140 invertor->SetLowerThreshold(0);
1141 invertor->SetUpperThreshold(0);
1142 invertor->SetInsideValue(1);
1144 typename InternalImageType::Pointer lungs_low_inv=invertor->GetOutput();
1146 //---------------------------------
1147 // Calculate lung volume
1148 //---------------------------------
1149 typename LabelStatisticsImageFilterType::Pointer statisticsFilter= LabelStatisticsImageFilterType::New();
1150 statisticsFilter->SetInput(lungs_low_inv);
1151 statisticsFilter->SetLabelInput(lungs_low);
1152 statisticsFilter->Update();
1153 unsigned int volume=statisticsFilter->GetSum(0);
1155 // Prepare threshold
1156 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1157 thresholder->SetUpperThreshold( 0.0 );
1158 thresholder->SetOutsideValue( 0 );
1159 thresholder->SetInsideValue( 1 );
1163 unsigned int totalNumberOfIterations=0;
1164 unsigned int coverage=0;
1167 // //---------------------------------
1168 // // Script for making movie
1169 // //---------------------------------
1170 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1173 levelSetFilter->Update();
1174 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1176 thresholder->SetInput( levelSetFilter->GetOutput() );
1177 thresholder->Update();
1178 statisticsFilter->SetInput(thresholder->GetOutput());
1179 statisticsFilter->Update();
1180 coverage=statisticsFilter->GetSum(0);
1182 // Compare the volumes
1183 if ( coverage >= m_ArgsInfo.fillingLevel_arg * volume/100 ) {
1184 if (m_Verbose) std::cout<<"Lungs filled for "<< m_ArgsInfo.fillingLevel_arg<<"% !"<<std::endl;
1187 if (m_Verbose) std::cout<<"After "<<totalNumberOfIterations<<" iterations, lungs are covered for "
1188 <<(double)coverage/volume*100.0<<"%..."<<std::endl;
1189 levelSetFilter->SetInput(levelSetFilter->GetOutput());
1190 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose);
1193 // std::ostringstream number_str;
1194 // number_str << number;
1195 // std::string param = number_str.str();
1196 // system((script+param).c_str());
1200 if ( (totalNumberOfIterations> 30000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS2_arg) ) break;
1204 if (m_Verbose) std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1205 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1206 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1209 thresholder->SetInput( levelSetFilter->GetOutput() );
1210 thresholder->Update();
1211 filledRibs=thresholder->GetOutput();
1212 // writeImage<InternalImageType>(filledRibs,"/home/jef/tmp/filled_ribs_before.mhd", m_Verbose);
1216 //---------------------------------
1217 // Write the filled ribs
1218 //---------------------------------
1219 if (m_ArgsInfo.writeFilledRibs_given) {
1220 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1221 caster->SetInput(filledRibs);
1222 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeFilledRibs_arg, m_Verbose);
1226 //----------------------------------------------------------------------------------------------------
1227 // Collapse to the lungs
1228 //----------------------------------------------------------------------------------------------------
1230 std::cout<<std::endl;
1231 std::cout<<"=========================================="<<std::endl;
1232 std::cout<<"|| Collapsing to the lung image ||"<<std::endl;
1233 std::cout<<"=========================================="<<std::endl;
1236 //---------------------------------
1237 // Make feature image air+bones
1238 //---------------------------------
1239 if (m_Verbose) std::cout<<"Making feature images..."<<std::endl;
1240 typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
1241 setBackgroundFilter->SetInput(air);
1242 setBackgroundFilter->SetInput2(bones);
1243 setBackgroundFilter->SetMaskValue(0);
1244 setBackgroundFilter->SetOutsideValue(0);
1245 setBackgroundFilter->Update();
1246 typename InternalImageType::Pointer bones_air =setBackgroundFilter->GetOutput();
1248 typename SetBackgroundFilterType::Pointer setBackgroundFilter2 = SetBackgroundFilterType::New();
1249 setBackgroundFilter2->SetInput(bones_air);
1250 setBackgroundFilter2->SetInput2(lungs);
1251 setBackgroundFilter2->SetMaskValue(0);
1252 setBackgroundFilter2->SetOutsideValue(0);
1253 setBackgroundFilter2->Update();
1254 typename InternalImageType::Pointer lungs_bones_air =setBackgroundFilter2->GetOutput();
1256 //---------------------------------
1257 // Prepare previous solution
1258 //---------------------------------
1259 if (m_Verbose) std::cout<<"Resampling previous solution..."<<std::endl;
1260 typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
1261 typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
1262 typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
1263 resampler->SetOutputStartIndex(bones->GetLargestPossibleRegion().GetIndex());
1264 resampler->SetInput(filledRibs);
1265 resampler->SetInterpolator(interpolator);
1266 resampler->SetOutputParametersFromImage(bones);
1267 resampler->Update();
1268 filledRibs =resampler->GetOutput();
1270 if (m_Verbose) std::cout<<"Setting lungs to 1..."<<std::endl;
1271 typename SetBackgroundFilterType::Pointer setBackgroundFilter3 = SetBackgroundFilterType::New();
1272 setBackgroundFilter3->SetInput(filledRibs);
1273 setBackgroundFilter3->SetInput2(lungs);
1274 setBackgroundFilter3->SetMaskValue(0);
1275 setBackgroundFilter3->SetOutsideValue(1);
1276 setBackgroundFilter3->Update();
1277 filledRibs=setBackgroundFilter3->GetOutput();
1279 if (m_Verbose) std::cout<<"Removing overlap with bones..."<<std::endl;
1280 typename SetBackgroundFilterType::Pointer setBackgroundFilter4 = SetBackgroundFilterType::New();
1281 setBackgroundFilter4->SetInput(filledRibs);
1282 setBackgroundFilter4->SetInput2(bones);
1283 setBackgroundFilter4->SetMaskValue(0);
1284 setBackgroundFilter4->SetOutsideValue(0);
1285 setBackgroundFilter4->Update();
1286 filledRibs =setBackgroundFilter4->GetOutput();
1287 // writeImage<InternalImageType>(filledRibs,"/home/jef/tmp/filledRibs_pp.mhd");
1288 //---------------------------------
1289 // Calculate Distance Map
1290 //---------------------------------
1291 // typedef itk::ApproximateSignedDistanceMapImageFilter <InternalImageType, LevelSetImageType> DistanceMapImageFilterType2;
1292 typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
1293 distanceMapImageFilter->SetInput(filledRibs);
1294 distanceMapImageFilter->SetInsideIsPositive(false);
1295 // distanceMapImageFilter->SetInsideValue(0);
1296 // distanceMapImageFilter->SetOutsideValue(1);
1297 if (m_Verbose) std::cout<<"Calculating distance map..."<<std::endl;
1298 distanceMapImageFilter->Update();
1300 //---------------------------------
1302 //---------------------------------
1303 typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
1304 levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() );
1305 levelSetFilter->SetFeatureImage( lungs_bones_air );
1306 levelSetFilter->SetPropagationScaling(m_ArgsInfo.prop3_arg);
1307 levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve3_arg );
1308 levelSetFilter->SetAdvectionScaling( 0 );
1309 levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS3_arg );
1310 levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter3_arg );
1311 levelSetFilter->SetUseImageSpacing(true);
1314 // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1317 if (m_Verbose) std::cout<<"Starting the levelset..."<<std::endl;
1318 unsigned int totalNumberOfIterations=0;
1320 levelSetFilter->Update();
1323 totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1324 levelSetFilter->SetInput(levelSetFilter->GetOutput());
1325 if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg);
1326 std::cout << "After "<<totalNumberOfIterations<<" iteration the RMS change is " << levelSetFilter->GetRMSChange() <<"..."<< std::endl;
1329 // std::ostringstream number_str;
1330 // number_str << number;
1331 // std::string param = number_str.str();
1332 // system((script+param).c_str());
1335 // stopping condition
1336 if ( (totalNumberOfIterations>= (unsigned int) m_ArgsInfo.maxIter3_arg) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS3_arg) ) break;
1337 levelSetFilter->SetNumberOfIterations( std::min ((unsigned int) m_ArgsInfo.iter3_arg, (unsigned int) m_ArgsInfo.maxIter3_arg-totalNumberOfIterations ) );
1342 std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1343 std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1344 std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1348 typedef itk::BinaryThresholdImageFilter< LevelSetImageType,InternalImageType> LevelSetBinarizeFilterType;
1349 typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1350 thresholder->SetUpperThreshold( 0.0 );
1351 thresholder->SetOutsideValue( 0 );
1352 thresholder->SetInsideValue( 1 );
1353 thresholder->SetInput( levelSetFilter->GetOutput() );
1354 thresholder->Update();
1355 typename InternalImageType::Pointer output = thresholder->GetOutput();
1358 //----------------------------------------------------------------------------------------------------
1359 // Prepare the output
1360 //----------------------------------------------------------------------------------------------------
1362 //---------------------------------
1364 //---------------------------------
1365 if (m_Verbose) std::cout<<"Removing overlap mask with air..."<<std::endl;
1366 typename SetBackgroundFilterType::Pointer setBackgroundFilter5 = SetBackgroundFilterType::New();
1367 setBackgroundFilter5->SetInput(output);
1368 setBackgroundFilter5->SetInput2(air);
1369 setBackgroundFilter5->SetMaskValue(0);
1370 setBackgroundFilter5->SetOutsideValue(0);
1371 setBackgroundFilter5->Update();
1372 output=setBackgroundFilter5->GetOutput();
1374 //---------------------------------
1375 // Open & close to cleanup
1376 //---------------------------------
1377 if ( m_ArgsInfo.openClose_flag) {
1379 //---------------------------------
1380 // Structuring element
1381 //---------------------------------
1382 typedef itk::BinaryBallStructuringElement<InternalPixelType,Dimension > KernelType;
1383 KernelType structuringElement;
1384 structuringElement.SetRadius(1);
1385 structuringElement.CreateStructuringElement();
1387 //---------------------------------
1389 //---------------------------------
1390 typedef itk::BinaryMorphologicalOpeningImageFilter<InternalImageType, InternalImageType , KernelType> OpenFilterType;
1391 typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
1392 openFilter->SetInput(output);
1393 openFilter->SetBackgroundValue(0);
1394 openFilter->SetForegroundValue(1);
1395 openFilter->SetKernel(structuringElement);
1396 if(m_Verbose) std::cout<<"Opening the output image..."<<std::endl;
1398 //---------------------------------
1400 //---------------------------------
1401 typedef itk::BinaryMorphologicalClosingImageFilter<InternalImageType, InternalImageType , KernelType> CloseFilterType;
1402 typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
1403 closeFilter->SetInput(openFilter->GetOutput());
1404 closeFilter->SetSafeBorder(true);
1405 closeFilter->SetForegroundValue(1);
1406 closeFilter->SetKernel(structuringElement);
1407 if(m_Verbose) std::cout<<"Closing the output image..."<<std::endl;
1408 closeFilter->Update();
1409 output=closeFilter->GetOutput();
1412 //writeImage<InternalImageType>(output,"/home/jef/tmp/mm_double.mhd");
1414 // Extract the upper part
1415 typedef itk::CropImageFilter<InternalImageType, InternalImageType> CropImageFilterType;
1416 typename CropImageFilterType::Pointer cropFilter=CropImageFilterType::New();
1417 cropFilter->SetInput(output);
1418 typename InputImageType::SizeType lSize, uSize;
1421 lSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
1422 cropFilter->SetLowerBoundaryCropSize(lSize);
1423 cropFilter->SetUpperBoundaryCropSize(uSize);
1424 cropFilter->Update();
1427 typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1428 caster->SetInput(cropFilter->GetOutput());
1429 writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.output_arg, m_Verbose);
1435 #endif //#define clitkMotionMaskGenericFilter_txx