X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=segmentation%2FclitkMotionMaskGenericFilter.txx;h=fd21b4629d4da99590b741e4fe7fa9ff634cfb50;hb=HEAD;hp=0b8860f61a51fa6b3533ac55c24c49a97e3b11dd;hpb=fa1e89be40a1eb3279719490b736a118723fa314;p=clitk.git diff --git a/segmentation/clitkMotionMaskGenericFilter.txx b/segmentation/clitkMotionMaskGenericFilter.txx old mode 100755 new mode 100644 index 0b8860f..fd21b46 --- a/segmentation/clitkMotionMaskGenericFilter.txx +++ b/segmentation/clitkMotionMaskGenericFilter.txx @@ -1,9 +1,9 @@ /*========================================================================= Program: vv http://www.creatis.insa-lyon.fr/rio/vv - Authors belong to: + Authors belong to: - University of LYON http://www.universite-lyon.fr/ - - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr + - Léon Bérard cancer center http://www.centreleonberard.fr - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr This software is distributed WITHOUT ANY WARRANTY; without even @@ -14,1391 +14,1435 @@ - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html -======================================================================-====*/ +===========================================================================**/ #ifndef clitkMotionMaskGenericFilter_txx #define clitkMotionMaskGenericFilter_txx /* ================================================= * @file clitkMotionMaskGenericFilter.txx - * @author - * @date - * - * @brief - * + * @author + * @date + * + * @brief + * ===================================================*/ +#include "itkLabelImageToShapeLabelMapFilter.h" +#include "itkShapeLabelObject.h" namespace clitk { - //------------------------------------------------------------------- - // Update with the number of dimensions - //------------------------------------------------------------------- - template - void - MotionMaskGenericFilter::UpdateWithDim(std::string PixelType) - { - if (m_Verbose) std::cout << "Image was detected to be "<(); - } - // else if(PixelType == "unsigned_short"){ - // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl; - // UpdateWithDimAndPixelType(); - // } - - // else if (PixelType == "unsigned_char"){ - // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl; - // UpdateWithDimAndPixelType(); - // } - - // else if (PixelType == "char"){ - // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl; - // UpdateWithDimAndPixelType(); - // } - else { - if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl; - UpdateWithDimAndPixelType(); - } +//------------------------------------------------------------------- +// Update with the number of dimensions +//------------------------------------------------------------------- +template +void +MotionMaskGenericFilter::UpdateWithDim(std::string PixelType) +{ + if (m_Verbose) std::cout << "Image was detected to be "<(); + } + // else if(PixelType == "unsigned_short"){ + // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl; + // UpdateWithDimAndPixelType(); + // } + + // else if (PixelType == "unsigned_char"){ + // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl; + // UpdateWithDimAndPixelType(); + // } + + // else if (PixelType == "char"){ + // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl; + // UpdateWithDimAndPixelType(); + // } + else { + if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl; + UpdateWithDimAndPixelType(); } +} + +//------------------------------------------------------------------- +// Air +//------------------------------------------------------------------- +template +typename itk::Image::Pointer +MotionMaskGenericFilter::GetAirImage(typename itk::Image::Pointer input, + typename itk::Image::Pointer lungs) +{ + + // ImageTypes + typedef int InternalPixelType; + typedef itk::Image InputImageType; + typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char? + + //---------------------------------------------------------------------------------------------------- + // Typedef for Processing + //---------------------------------------------------------------------------------------------------- + typedef itk::ImageFileReader FeatureReaderType; + typedef itk::BinaryThresholdImageFilter InputBinarizeFilterType; + typedef itk::BinaryThresholdImageFilter InversionFilterType; + typedef itk::ThresholdImageFilter ThresholdFilterType; + typedef itk::ConnectedComponentImageFilter ConnectFilterType; + typedef itk::RelabelComponentImageFilter RelabelFilterType; + typedef itk::MirrorPadImageFilter MirrorPadImageFilterType; + + + typename InternalImageType::Pointer air; + if (m_ArgsInfo.featureAir_given) { + typename FeatureReaderType::Pointer featureReader =FeatureReaderType::New(); + featureReader->SetFileName(m_ArgsInfo.featureAir_arg); + if (m_Verbose) std::cout<<"Reading the air feature image..."<Update(); + air=featureReader->GetOutput(); - //------------------------------------------------------------------- - // Air - //------------------------------------------------------------------- - template - typename itk::Image::Pointer - MotionMaskGenericFilter::GetAirImage(typename itk::Image::Pointer input ) - { - - // ImageTypes - typedef int InternalPixelType; - typedef itk::Image InputImageType; - typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char? - - //---------------------------------------------------------------------------------------------------- - // Typedef for Processing - //---------------------------------------------------------------------------------------------------- - typedef itk::ImageFileReader FeatureReaderType; - typedef itk::BinaryThresholdImageFilter InputBinarizeFilterType; - typedef itk::BinaryThresholdImageFilter InversionFilterType; - typedef itk::ThresholdImageFilter ThresholdFilterType; - typedef itk::ConnectedComponentImageFilter ConnectFilterType; - typedef itk::RelabelComponentImageFilter RelabelFilterType; - typedef itk::MirrorPadImageFilter MirrorPadImageFilterType; - - - typename InternalImageType::Pointer air; - if (m_ArgsInfo.featureAir_given) - { - typename FeatureReaderType::Pointer featureReader =FeatureReaderType::New(); - featureReader->SetFileName(m_ArgsInfo.featureAir_arg); - if (m_Verbose) std::cout<<"Reading the air feature image..."<Update(); - air=featureReader->GetOutput(); - } - else - { - if (m_Verbose) std::cout<<"Extracting the air feature image..."<SetInput(input); - binarizeFilter->SetLowerThreshold(static_cast(m_ArgsInfo.lowerThresholdAir_arg)); - binarizeFilter->SetUpperThreshold(static_cast(m_ArgsInfo.upperThresholdAir_arg)); - if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<SetInput(binarizeFilter->GetOutput()); - connectFilter->SetBackgroundValue(0); - connectFilter->SetFullyConnected(false); - if (m_Verbose) std::cout<<"Labeling the connected components..."<SetInput(connectFilter->GetOutput()); - if (m_Verbose) std::cout<<"Sorting the labels..."<SetInput(relabelFilter->GetOutput()); - thresholdFilter->SetUpper(1); - if (m_Verbose) std::cout<<"Keeping the first label..."<Update(); - air=thresholdFilter->GetOutput(); - } - //--------------------------------- - // Invert + // Pad //--------------------------------- - typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New(); - inversionFilter->SetInput(air); - inversionFilter->SetLowerThreshold(0); - inversionFilter->SetUpperThreshold(0); - inversionFilter->SetInsideValue(1); - inversionFilter->Update(); - + if(m_ArgsInfo.pad_flag) { + typedef itk::ImageRegionIteratorWithIndex IteratorType; + IteratorType it(air, air->GetLargestPossibleRegion()); + typename InternalImageType::IndexType index; + while(!it.IsAtEnd()) { + index=it.GetIndex(); + + //Pad borders of all dimensions but 2 (cranio-caudal) + for (unsigned int i=0; iGetLargestPossibleRegion().GetIndex()[i] + || index[i]==(unsigned int)air->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air->GetLargestPossibleRegion().GetSize()[i]-1) + it.Set(1); + } + ++it; + } + } + } else { + if (m_Verbose) std::cout<<"Extracting the air feature image..."<SetInput(inversionFilter->GetOutput()); - typename InternalImageType::SizeType padSize; - padSize.Fill(0); - padSize[2]=input->GetLargestPossibleRegion().GetSize()[2]; - mirrorPadImageFilter->SetPadLowerBound(padSize); - if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<Update(); - air=mirrorPadImageFilter->GetOutput(); - // writeImage(air,"/home/jef/tmp/air.mhd"); + typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New(); + binarizeFilter->SetInput(input); + binarizeFilter->SetLowerThreshold(static_cast(m_ArgsInfo.lowerThresholdAir_arg)); + binarizeFilter->SetUpperThreshold(static_cast(m_ArgsInfo.upperThresholdAir_arg)); + if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<Update(); + air = binarizeFilter->GetOutput(); //--------------------------------- // Pad //--------------------------------- - if(m_ArgsInfo.pad_flag) - { - typedef itk::ImageRegionIteratorWithIndex IteratorType; - IteratorType it(air, air->GetLargestPossibleRegion()); - typename InternalImageType::IndexType index; - while(!it.IsAtEnd()) - { - index=it.GetIndex(); - for (unsigned int i=0;iGetLargestPossibleRegion().GetIndex()[i] - || index[i]==(unsigned int)air->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air->GetLargestPossibleRegion().GetSize()[i]-1) - it.Set(0); - ++it; - } + if(m_ArgsInfo.pad_flag) { + typedef itk::ImageRegionIteratorWithIndex IteratorType; + IteratorType it(air, air->GetLargestPossibleRegion()); + typename InternalImageType::IndexType index; + while(!it.IsAtEnd()) { + index=it.GetIndex(); + + //Pad borders of all dimensions but 2 (cranio-caudal) + for (unsigned int i=0; iGetLargestPossibleRegion().GetIndex()[i] + || index[i]==(unsigned int)air->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air->GetLargestPossibleRegion().GetSize()[i]-1) + it.Set(binarizeFilter->GetInsideValue()); + } + ++it; } - - return air; - } + } - - //------------------------------------------------------------------- - // Bones - //------------------------------------------------------------------- - template - typename itk::Image::Pointer - MotionMaskGenericFilter::GetBonesImage(typename itk::Image::Pointer input ) - { - - // ImageTypes - typedef int InternalPixelType; - typedef itk::Image InputImageType; - typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char? - - //---------------------------------------------------------------------------------------------------- - // Typedef for Processing - //---------------------------------------------------------------------------------------------------- - typedef itk::ImageFileReader FeatureReaderType; - typedef itk::BinaryThresholdImageFilter InputBinarizeFilterType; - typedef itk::BinaryThresholdImageFilter InversionFilterType; - typedef itk::ThresholdImageFilter ThresholdFilterType; - typedef itk::ConnectedComponentImageFilter ConnectFilterType; - typedef itk::RelabelComponentImageFilter RelabelFilterType; - typedef itk::MirrorPadImageFilter MirrorPadImageFilterType; - - - typename InternalImageType::Pointer bones; - if (m_ArgsInfo.featureBones_given) - { - typename FeatureReaderType::Pointer featureReader =FeatureReaderType::New(); - featureReader->SetFileName(m_ArgsInfo.featureBones_arg); - if (m_Verbose) std::cout<<"Reading the bones feature image..."<Update(); - bones=featureReader->GetOutput(); - } - else - { - if (m_Verbose) std::cout<<"Extracting the bones feature image..."<SetInput(input); - binarizeFilter->SetLowerThreshold(static_cast(m_ArgsInfo.lowerThresholdBones_arg)); - binarizeFilter->SetUpperThreshold(static_cast(m_ArgsInfo.upperThresholdBones_arg)); - if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<SetInput(binarizeFilter->GetOutput()); - connectFilter->SetBackgroundValue(0); - connectFilter->SetFullyConnected(false); - if (m_Verbose) std::cout<<"Labeling the connected components..."<SetInput(connectFilter->GetOutput()); - if (m_Verbose) std::cout<<"Sorting the labels..."<SetInput(relabelFilter->GetOutput()); - thresholdFilter->SetUpper(1); - if (m_Verbose) std::cout<<"Keeping the first label..."<Update(); - bones=thresholdFilter->GetOutput(); - - } - //--------------------------------- - // Invert + // Remove lungs (in place) + //--------------------------------- + typedef itk::ImageRegionIterator IteratorType; + IteratorType itAir(air, binarizeFilter->GetOutput()->GetLargestPossibleRegion()); + IteratorType itLungs(lungs, binarizeFilter->GetOutput()->GetLargestPossibleRegion()); //lungs padded, used air region + itAir.GoToBegin(); + itLungs.GoToBegin(); + while(!itAir.IsAtEnd()) { + if(!itLungs.Get()) // The lungs are set to 0 at that stage + itAir.Set(0); + ++itAir; + ++itLungs; + } + + //--------------------------------- + // Label the connected components + //--------------------------------- + typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New(); + connectFilter->SetInput(binarizeFilter->GetOutput()); + connectFilter->SetBackgroundValue(0); + connectFilter->SetFullyConnected(false); + if (m_Verbose) std::cout<<"Labeling the connected components..."<SetInput(bones); - inversionFilter->SetLowerThreshold(0); - inversionFilter->SetUpperThreshold(0); - inversionFilter->SetInsideValue(1); - inversionFilter->Update(); - + typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New(); + relabelFilter->SetInput(connectFilter->GetOutput()); + if (m_Verbose) std::cout<<"Sorting the labels..."<SetInput(inversionFilter->GetOutput()); - typename InternalImageType::SizeType padSize; - padSize.Fill(0); - padSize[2]=input->GetLargestPossibleRegion().GetSize()[2]; - mirrorPadImageFilter->SetPadLowerBound(padSize); - if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<Update(); - bones=mirrorPadImageFilter->GetOutput(); - // writeImage(bones,"/home/jef/tmp/bones.mhd"); - - return bones; + typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New(); + thresholdFilter->SetInput(relabelFilter->GetOutput()); + thresholdFilter->SetUpper(1); + if (m_Verbose) std::cout<<"Keeping the first label..."<Update(); + air=thresholdFilter->GetOutput(); } + //--------------------------------- + // Invert + //--------------------------------- + typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New(); + inversionFilter->SetInput(air); + inversionFilter->SetLowerThreshold(0); + inversionFilter->SetUpperThreshold(0); + inversionFilter->SetInsideValue(1); + inversionFilter->Update(); + + //--------------------------------- + // Mirror + //--------------------------------- + typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New(); + mirrorPadImageFilter->SetInput(inversionFilter->GetOutput()); + typename InternalImageType::SizeType padSize; + padSize.Fill(0); + padSize[2]=input->GetLargestPossibleRegion().GetSize()[2]; + mirrorPadImageFilter->SetPadLowerBound(padSize); + if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<Update(); + air=mirrorPadImageFilter->GetOutput(); + //writeImage(air,"/home/srit/tmp/air.mhd"); + + return air; +} + + +//------------------------------------------------------------------- +// Bones +//------------------------------------------------------------------- +template +typename itk::Image::Pointer +MotionMaskGenericFilter::GetBonesImage(typename itk::Image::Pointer input ) +{ + // ImageTypes + typedef int InternalPixelType; + typedef itk::Image InputImageType; + typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char? + + //---------------------------------------------------------------------------------------------------- + // Typedef for Processing + //---------------------------------------------------------------------------------------------------- + typedef itk::ImageFileReader FeatureReaderType; + typedef itk::BinaryThresholdImageFilter InputBinarizeFilterType; + typedef itk::BinaryThresholdImageFilter InversionFilterType; + typedef itk::ThresholdImageFilter ThresholdFilterType; + typedef itk::ConnectedComponentImageFilter ConnectFilterType; + typedef itk::RelabelComponentImageFilter RelabelFilterType; + typedef itk::MirrorPadImageFilter MirrorPadImageFilterType; + + + typename InternalImageType::Pointer bones; + if (m_ArgsInfo.featureBones_given) { + typename FeatureReaderType::Pointer featureReader =FeatureReaderType::New(); + featureReader->SetFileName(m_ArgsInfo.featureBones_arg); + if (m_Verbose) std::cout<<"Reading the bones feature image..."<Update(); + bones=featureReader->GetOutput(); + } else { + if (m_Verbose) std::cout<<"Extracting the bones feature image..."<SetInput(input); + binarizeFilter->SetLowerThreshold(static_cast(m_ArgsInfo.lowerThresholdBones_arg)); + binarizeFilter->SetUpperThreshold(static_cast(m_ArgsInfo.upperThresholdBones_arg)); + if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<SetInput(binarizeFilter->GetOutput()); + connectFilter->SetBackgroundValue(0); + connectFilter->SetFullyConnected(false); + if (m_Verbose) std::cout<<"Labeling the connected components..."< - typename itk::Image::Pointer - MotionMaskGenericFilter::GetLungsImage(typename itk::Image::Pointer input ) - { - // ImageTypes - typedef int InternalPixelType; - typedef itk::Image InputImageType; - typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char? - - //---------------------------------------------------------------------------------------------------- - // Typedef for Processing - //---------------------------------------------------------------------------------------------------- - typedef itk::ImageFileReader FeatureReaderType; - typedef itk::BinaryThresholdImageFilter InputBinarizeFilterType; - typedef itk::BinaryThresholdImageFilter InversionFilterType; - typedef itk::ThresholdImageFilter ThresholdFilterType; - typedef itk::ConnectedComponentImageFilter ConnectFilterType; - typedef itk::RelabelComponentImageFilter RelabelFilterType; - typedef itk::MirrorPadImageFilter MirrorPadImageFilterType; - - typename InternalImageType::Pointer lungs; - if (m_ArgsInfo.featureLungs_given) - { - typename FeatureReaderType::Pointer featureReader =FeatureReaderType::New(); - featureReader->SetFileName(m_ArgsInfo.featureLungs_arg); - if (m_Verbose) std::cout<<"Reading the lungs feature image..."<Update(); - lungs=featureReader->GetOutput(); - } - else - { - if (m_Verbose) std::cout<<"Extracting the lungs feature image..."<SetInput(input); - binarizeFilter->SetLowerThreshold(static_cast(m_ArgsInfo.lowerThresholdLungs_arg)); - binarizeFilter->SetUpperThreshold(static_cast(m_ArgsInfo.upperThresholdLungs_arg)); - if (m_Verbose) std::cout<<"Binarizing the image using a threshold "<SetInput(binarizeFilter->GetOutput()); - connectFilter->SetBackgroundValue(0); - connectFilter->SetFullyConnected(true); - if (m_Verbose) std::cout<<"Labeling the connected components..."<SetInput(connectFilter->GetOutput()); - if (m_Verbose) std::cout<<"Sorting the labels..."< (relabelFilter->GetOutput(), "/home/jef/tmp/labels.mhd"); - - //--------------------------------- - // Keep the label - //--------------------------------- - typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New(); - thresholdFilter->SetInput(relabelFilter->GetOutput()); - thresholdFilter->SetLower(1); - thresholdFilter->SetUpper(2); - if (m_Verbose) std::cout<<"Keeping the first two labels..."<Update(); - lungs=thresholdFilter->GetOutput(); - - } - - //--------------------------------- - // Invert + // Sort the labels according to size //--------------------------------- - typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New(); - inversionFilter->SetInput(lungs); - inversionFilter->SetLowerThreshold(0); - inversionFilter->SetUpperThreshold(0); - inversionFilter->SetInsideValue(1); - inversionFilter->Update(); - + typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New(); + relabelFilter->SetInput(connectFilter->GetOutput()); + if (m_Verbose) std::cout<<"Sorting the labels..."<SetInput(inversionFilter->GetOutput()); - typename InternalImageType::SizeType padSize; - padSize.Fill(0); - padSize[2]=input->GetLargestPossibleRegion().GetSize()[2]; - mirrorPadImageFilter->SetPadLowerBound(padSize); - if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<Update(); - lungs=mirrorPadImageFilter->GetOutput(); - // writeImage(lungs,"/home/jef/tmp/lungs.mhd"); - - return lungs; + typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New(); + thresholdFilter->SetInput(relabelFilter->GetOutput()); + thresholdFilter->SetUpper(1); + if (m_Verbose) std::cout<<"Keeping the first label..."<Update(); + bones=thresholdFilter->GetOutput(); + } - //------------------------------------------------------------------- - // Resample - //------------------------------------------------------------------- - template - typename itk::Image::Pointer - MotionMaskGenericFilter::Resample( typename itk::Image::Pointer input ) - { + //--------------------------------- + // Invert + //--------------------------------- + typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New(); + inversionFilter->SetInput(bones); + inversionFilter->SetLowerThreshold(0); + inversionFilter->SetUpperThreshold(0); + inversionFilter->SetInsideValue(1); + inversionFilter->Update(); + + //--------------------------------- + // Mirror + //--------------------------------- + typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New(); + mirrorPadImageFilter->SetInput(inversionFilter->GetOutput()); + typename InternalImageType::SizeType padSize; + padSize.Fill(0); + padSize[2]=input->GetLargestPossibleRegion().GetSize()[2]; + mirrorPadImageFilter->SetPadLowerBound(padSize); + if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<Update(); + bones=mirrorPadImageFilter->GetOutput(); + // writeImage(bones,"/home/jef/tmp/bones.mhd"); + + return bones; +} + + + + +//------------------------------------------------------------------- +// Lungs +//------------------------------------------------------------------- +template +typename itk::Image::Pointer +MotionMaskGenericFilter::GetLungsImage(typename itk::Image::Pointer input ) +{ + // ImageTypes + typedef int InternalPixelType; + typedef itk::Image InputImageType; + typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char? + + //---------------------------------------------------------------------------------------------------- + // Typedef for Processing + //---------------------------------------------------------------------------------------------------- + typedef itk::ImageFileReader FeatureReaderType; + typedef itk::BinaryThresholdImageFilter InputBinarizeFilterType; + typedef itk::BinaryThresholdImageFilter InversionFilterType; + typedef itk::ThresholdImageFilter ThresholdFilterType; + typedef itk::ConnectedComponentImageFilter ConnectFilterType; + typedef itk::RelabelComponentImageFilter RelabelFilterType; + typedef itk::MirrorPadImageFilter MirrorPadImageFilterType; + + typename InternalImageType::Pointer lungs; + if (m_ArgsInfo.featureLungs_given) { + typename FeatureReaderType::Pointer featureReader =FeatureReaderType::New(); + featureReader->SetFileName(m_ArgsInfo.featureLungs_arg); + if (m_Verbose) std::cout<<"Reading the lungs feature image..."<Update(); + lungs=featureReader->GetOutput(); + } else { + if (m_Verbose) std::cout<<"Extracting the lungs feature image..."<SetInput(input); + binarizeFilter->SetLowerThreshold(static_cast(m_ArgsInfo.lowerThresholdLungs_arg)); + binarizeFilter->SetUpperThreshold(static_cast(m_ArgsInfo.upperThresholdLungs_arg)); + if (m_Verbose) std::cout<<"Binarizing the image using a threshold "< InputImageType; - typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char? + //--------------------------------- + // Label the connected components + //--------------------------------- + typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New(); + connectFilter->SetInput(binarizeFilter->GetOutput()); + connectFilter->SetBackgroundValue(0); + connectFilter->SetFullyConnected(true); + if (m_Verbose) std::cout<<"Labeling the connected components..."<Update(); + if (m_Verbose) std::cout<<"found "<< connectFilter->GetObjectCount() << std::endl; //--------------------------------- - // Resample to new spacing + // Sort the labels according to size //--------------------------------- - typename InternalImageType::SizeType size_low; - typename InternalImageType::SpacingType spacing_low; - for (unsigned int i=0; iGetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i]/spacing_low[i]); - } - else - for (unsigned int i=0; iGetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i]/spacing_low[i]); - } - - typename InternalImageType::PointType origin; - input->TransformIndexToPhysicalPoint(input->GetLargestPossibleRegion().GetIndex(), origin); - typedef itk::ResampleImageFilter ResampleImageFilterType; - typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New(); - typedef itk::NearestNeighborInterpolateImageFunction InterpolatorType; - typename InterpolatorType::Pointer interpolator=InterpolatorType::New(); - resampler->SetInterpolator(interpolator); - resampler->SetOutputOrigin(origin); - resampler->SetOutputSpacing(spacing_low); - resampler->SetSize(size_low); - resampler->SetInput(input); - resampler->Update(); - typename InternalImageType::Pointer output =resampler->GetOutput(); - return output; + typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New(); + relabelFilter->SetInput(connectFilter->GetOutput()); + if (m_Verbose) std::cout<<"Sorting the labels..."< (relabelFilter->GetOutput(), "/home/vdelmon/tmp/labels.mhd"); + + //--------------------------------- + // Keep the label + //--------------------------------- + typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New(); + thresholdFilter->SetInput(relabelFilter->GetOutput()); + thresholdFilter->SetLower(1); + thresholdFilter->SetUpper(2); + if (m_Verbose) std::cout<<"Keeping the first two labels..."<Update(); + lungs=thresholdFilter->GetOutput(); + } - //------------------------------------------------------------------- - // Initialize ellips - //------------------------------------------------------------------- - template - typename itk::Image::Pointer - MotionMaskGenericFilter::InitializeEllips( typename itk::Vector center, typename itk::Image::Pointer bones_low ) - { - - // ImageTypes - typedef int InternalPixelType; - typedef itk::Image InputImageType; - typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char? - typedef itk::Image< unsigned char , Dimension> OutputImageType; - typedef itk::ImageRegionIteratorWithIndex IteratorType; - typedef clitk::SetBackgroundImageFilter SetBackgroundFilterType; - typedef itk::LabelStatisticsImageFilter LabelStatisticsImageFilterType; - typedef itk::CastImageFilter CastImageFilterType; - - - typename InternalImageType::Pointer ellips; - if (m_ArgsInfo.ellips_given || m_ArgsInfo.grownEllips_given || m_ArgsInfo.filledRibs_given) - { - if(m_ArgsInfo.ellips_given) - { - typedef itk::ImageFileReader FeatureReaderType; - typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New(); - featureReader->SetFileName(m_ArgsInfo.ellips_arg); - featureReader->Update(); - ellips=featureReader->GetOutput(); - } + //--------------------------------- + // Invert + //--------------------------------- + typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New(); + inversionFilter->SetInput(lungs); + inversionFilter->SetLowerThreshold(0); + inversionFilter->SetUpperThreshold(0); + inversionFilter->SetInsideValue(1); + inversionFilter->Update(); + + //--------------------------------- + // Mirror + //--------------------------------- + typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New(); + mirrorPadImageFilter->SetInput(inversionFilter->GetOutput()); + typename InternalImageType::SizeType padSize; + padSize.Fill(0); + padSize[2]=input->GetLargestPossibleRegion().GetSize()[2]; + mirrorPadImageFilter->SetPadLowerBound(padSize); + if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<Update(); + lungs=mirrorPadImageFilter->GetOutput(); + // writeImage(lungs,"/home/jef/tmp/lungs.mhd"); + + return lungs; +} + +//------------------------------------------------------------------- +// Resample +//------------------------------------------------------------------- +template +typename itk::Image::Pointer +MotionMaskGenericFilter::Resample( typename itk::Image::Pointer input ) +{ + + typedef int InternalPixelType; + typedef itk::Image InputImageType; + typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char? + + //--------------------------------- + // Resample to new spacing + //--------------------------------- + typename InternalImageType::SizeType size_low; + typename InternalImageType::SpacingType spacing_low; + for (unsigned int i=0; iGetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i]/spacing_low[i]); } else - { - if(m_Verbose) - { - std::cout< offset; - if(m_ArgsInfo.offset_given== Dimension) - { - for(unsigned int i=0; i centerEllips=center+offset; - if (m_Verbose) - { - std::cout<<"Placing the center of the initial ellipsoide at (mm) "<SetRegions (bones_low->GetLargestPossibleRegion()); - ellips->SetOrigin(bones_low->GetOrigin()); - ellips->SetSpacing(bones_low->GetSpacing()); - ellips->Allocate(); - ellips->FillBuffer(0); - - // Axes - typename itk::Vector axes; - if (m_ArgsInfo.axes_given==Dimension) - for(unsigned int i=0; iGetLargestPossibleRegion()); - itEllips.GoToBegin(); - typename InternalImageType::PointType point; - typename InternalImageType::IndexType index; - double distance; - - if (m_Verbose) std::cout<<"Drawing the initial ellipsoide..."<TransformIndexToPhysicalPoint(index, point); - distance=0.0; - for(unsigned int i=0; iGetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i]/spacing_low[i]); } - - + + typename InternalImageType::PointType origin; + input->TransformIndexToPhysicalPoint(input->GetLargestPossibleRegion().GetIndex(), origin); + typedef itk::ResampleImageFilter ResampleImageFilterType; + typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New(); + typedef itk::NearestNeighborInterpolateImageFunction InterpolatorType; + typename InterpolatorType::Pointer interpolator=InterpolatorType::New(); + resampler->SetInterpolator(interpolator); + resampler->SetOutputOrigin(origin); + resampler->SetOutputSpacing(spacing_low); + resampler->SetSize(size_low); + resampler->SetInput(input); + resampler->Update(); + typename InternalImageType::Pointer output =resampler->GetOutput(); + return output; +} + + +//------------------------------------------------------------------- +// Initialize ellips +//------------------------------------------------------------------- +template +typename itk::Image::Pointer +MotionMaskGenericFilter::InitializeEllips( typename itk::Vector center, typename itk::Image::Pointer bones_low, typename itk::Image::Pointer lungs_low ) +{ + + // ImageTypes + typedef int InternalPixelType; + typedef itk::Image InputImageType; + typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char? + typedef itk::Image< unsigned char , Dimension> OutputImageType; + typedef itk::ImageRegionIteratorWithIndex IteratorType; + typedef clitk::SetBackgroundImageFilter SetBackgroundFilterType; + typedef itk::LabelStatisticsImageFilter LabelStatisticsImageFilterType; + typedef itk::CastImageFilter CastImageFilterType; + + + typename InternalImageType::Pointer ellips; + if (m_ArgsInfo.ellips_given || m_ArgsInfo.grownEllips_given || m_ArgsInfo.filledRibs_given) { + if(m_ArgsInfo.ellips_given) { + typedef itk::ImageFileReader FeatureReaderType; + typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New(); + featureReader->SetFileName(m_ArgsInfo.ellips_arg); + featureReader->Update(); + ellips=featureReader->GetOutput(); + } + } else { + if(m_Verbose) { + std::cout< centerEllips; + typename itk::Vector axes; + if (m_ArgsInfo.ellipseAutoDetect_flag) { + if(m_Verbose) { + std::cout<<"Auto-detecting intial ellipse..."<GetLargestPossibleRegion(); + typename InternalImageType::SizeType size = region_before.GetSize(); + size[2] /= 2; + region.SetSize(size); + lungs_low->SetRegions(region); + + // compute statistics of the (mirroed) lung region in order to guess the params of the ellipse + typedef itk::LabelImageToShapeLabelMapFilter LabelImageToShapeLabelMapFilter; + typename LabelImageToShapeLabelMapFilter::Pointer label_filter = LabelImageToShapeLabelMapFilter::New(); + label_filter->SetInput(lungs_low); + label_filter->Update(); + + typename InternalImageType::SpacingType spacing = lungs_low->GetSpacing(); + typedef typename LabelImageToShapeLabelMapFilter::OutputImageType::LabelObjectType LabelType; + const LabelType* label = dynamic_cast(label_filter->GetOutput()->GetNthLabelObject(0)); + + // lungs' barycenter + typename LabelType::CentroidType lung_centroid = label->GetCentroid(); + + // try to guess ideal ellipse axes from the lungs' bounding box and centroid + // with some hard-coded "magic" constants... + typename LabelType::RegionType lung_bbox = label->GetBoundingBox(); + + axes[0] = 0.95*lung_bbox.GetSize()[0]*spacing[0]/2; + axes[1] = 0.3*lung_bbox.GetSize()[1]*spacing[1]/2; + axes[2] = 1.25*fabs(lung_centroid[2] - center[2]); + + // ellipse's center in XY is the center of the lungs' bounding box + typename InternalImageType::PointType origin = lungs_low->GetOrigin(); + centerEllips[0] = origin[0] + (lung_bbox.GetIndex()[0] + lung_bbox.GetSize()[0]/2)*spacing[0]; + centerEllips[1] = origin[1] + (lung_bbox.GetIndex()[1] + lung_bbox.GetSize()[1]/2)*spacing[1]; + centerEllips[2] = center[2]; + + lungs_low->SetRegions(region_before); + + if(m_Verbose) { + std::cout << "Lungs'centroid at " << lung_centroid << std::endl; + std::cout << "Ellipse center at " << centerEllips << std::endl; + std::cout << "Ellipse axes as " << axes << std::endl; + std::cout << "Lung's bounding box (index,size) " << lung_bbox.GetIndex() << lung_bbox.GetSize() << std::endl; + } + } + else { + //--------------------------------- + // Offset from center + //--------------------------------- + typename itk::Vector offset; + if(m_ArgsInfo.offset_given== Dimension) { + for(unsigned int i=0; iSetInput(ellips); - writeImage(caster->GetOutput(), m_ArgsInfo.writeEllips_arg, m_Verbose); - } + ellips=InternalImageType::New(); + ellips->SetRegions (bones_low->GetLargestPossibleRegion()); + ellips->SetOrigin(bones_low->GetOrigin()); + ellips->SetSpacing(bones_low->GetSpacing()); + ellips->Allocate(); + ellips->FillBuffer(0); + + // Draw an ellips + IteratorType itEllips(ellips, ellips->GetLargestPossibleRegion()); + itEllips.GoToBegin(); + typename InternalImageType::PointType point; + typename InternalImageType::IndexType index; + double distance; + + if (m_Verbose) std::cout<<"Drawing the initial ellipsoide..."<TransformIndexToPhysicalPoint(index, point); + distance=0.0; + for(unsigned int i=0; iSetInput(ellips); + writeImage(caster->GetOutput(), m_ArgsInfo.writeEllips_arg, m_Verbose); } + return ellips; + +} + + +//------------------------------------------------------------------- +// Update with the number of dimensions and the pixeltype +//------------------------------------------------------------------- +template +void +MotionMaskGenericFilter::UpdateWithDimAndPixelType() +{ + + // ImageTypes + typedef int InternalPixelType; + typedef itk::Image InputImageType; + typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char? + typedef itk::Image< float, Dimension> LevelSetImageType; //labeling doesn't work with unsigned char? + typedef itk::Image OutputImageType; + + + //---------------------------------------------------------------------------------------------------- + // Typedef for Processing + //---------------------------------------------------------------------------------------------------- + typedef itk::ImageFileReader FeatureReaderType; + typedef itk::BinaryThresholdImageFilter InputBinarizeFilterType; + typedef itk::BinaryThresholdImageFilter< LevelSetImageType,InternalImageType> LevelSetBinarizeFilterType; + typedef itk::BinaryThresholdImageFilter InversionFilterType; + typedef itk::ThresholdImageFilter ThresholdFilterType; + typedef itk::ConnectedComponentImageFilter ConnectFilterType; + typedef itk::RelabelComponentImageFilter RelabelFilterType; + typedef itk::MirrorPadImageFilter MirrorPadImageFilterType; + typedef itk::NearestNeighborInterpolateImageFunction InterpolatorType; + typedef itk::ResampleImageFilter ResampleImageFilterType; + typedef itk::ImageRegionIteratorWithIndex IteratorType; + typedef itk::GeodesicActiveContourLevelSetImageFilter LevelSetImageFilterType; + typedef itk::SignedDanielssonDistanceMapImageFilter DistanceMapImageFilterType; + typedef clitk::SetBackgroundImageFilter SetBackgroundFilterType; + typedef itk::LabelStatisticsImageFilter LabelStatisticsImageFilterType; + typedef itk::CastImageFilter CastImageFilterType; + + + // Read the input + typedef itk::ImageFileReader InputReaderType; + typename InputReaderType::Pointer reader = InputReaderType::New(); + reader->SetFileName( m_InputFileName); + reader->Update(); + typename InputImageType::Pointer input= reader->GetOutput(); + + // // globals for avi + // unsigned int number=0; + + + if(m_Verbose) { + std::cout<GetLungsImage(input); + + //------------------------------------------------------------------------------- + // Air + //------------------------------------------------------------------------------- + typename InternalImageType::Pointer air = this->GetAirImage(input, lungs); + + //------------------------------------------------------------------------------- + // Bones + //------------------------------------------------------------------------------- + typename InternalImageType::Pointer bones = this->GetBonesImage(input); + + //---------------------------------------------------------------------------------------------------- + // Write feature images + //---------------------------------------------------------------------------------------------------- + if(m_ArgsInfo.writeFeature_given) { + typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New(); + setBackgroundFilter->SetInput(air); + setBackgroundFilter->SetInput2(bones); + setBackgroundFilter->SetMaskValue(0); + setBackgroundFilter->SetOutsideValue(2); + setBackgroundFilter->Update(); + typename InternalImageType::Pointer bones_air =setBackgroundFilter->GetOutput(); + + typename SetBackgroundFilterType::Pointer setBackgroundFilter2 = SetBackgroundFilterType::New(); + setBackgroundFilter2->SetInput(bones_air); + setBackgroundFilter2->SetInput2(lungs); + setBackgroundFilter2->SetMaskValue(0); + setBackgroundFilter2->SetOutsideValue(3); + setBackgroundFilter2->Update(); + typename InternalImageType::Pointer lungs_bones_air =setBackgroundFilter2->GetOutput(); + + typename CastImageFilterType::Pointer caster=CastImageFilterType::New(); + caster->SetInput(lungs_bones_air); + writeImage(caster->GetOutput(), m_ArgsInfo.writeFeature_arg, m_Verbose); + } - //------------------------------------------------------------------- - // Update with the number of dimensions and the pixeltype - //------------------------------------------------------------------- - template - void - MotionMaskGenericFilter::UpdateWithDimAndPixelType() - { - - // ImageTypes - typedef int InternalPixelType; - typedef itk::Image InputImageType; - typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char? - typedef itk::Image< float, Dimension> LevelSetImageType; //labeling doesn't work with unsigned char? - typedef itk::Image OutputImageType; - - - //---------------------------------------------------------------------------------------------------- - // Typedef for Processing - //---------------------------------------------------------------------------------------------------- - typedef itk::ImageFileReader FeatureReaderType; - typedef itk::BinaryThresholdImageFilter InputBinarizeFilterType; - typedef itk::BinaryThresholdImageFilter< LevelSetImageType,InternalImageType> LevelSetBinarizeFilterType; - typedef itk::BinaryThresholdImageFilter InversionFilterType; - typedef itk::ThresholdImageFilter ThresholdFilterType; - typedef itk::ConnectedComponentImageFilter ConnectFilterType; - typedef itk::RelabelComponentImageFilter RelabelFilterType; - typedef itk::MirrorPadImageFilter MirrorPadImageFilterType; - typedef itk::NearestNeighborInterpolateImageFunction InterpolatorType; - typedef itk::ResampleImageFilter ResampleImageFilterType; + //---------------------------------------------------------------------------------------------------- + // Low dimensional versions + //---------------------------------------------------------------------------------------------------- + typename InternalImageType::Pointer bones_low =Resample(bones); + typename InternalImageType::Pointer lungs_low =Resample(lungs); + typename InternalImageType::Pointer air_low =Resample(air); + + //--------------------------------- + // Pad + //--------------------------------- + if(m_ArgsInfo.pad_flag) { typedef itk::ImageRegionIteratorWithIndex IteratorType; - typedef itk::GeodesicActiveContourLevelSetImageFilter LevelSetImageFilterType; - typedef itk::SignedDanielssonDistanceMapImageFilter DistanceMapImageFilterType; - typedef clitk::SetBackgroundImageFilter SetBackgroundFilterType; - typedef itk::LabelStatisticsImageFilter LabelStatisticsImageFilterType; - typedef itk::CastImageFilter CastImageFilterType; - - - // Read the input - typedef itk::ImageFileReader InputReaderType; - typename InputReaderType::Pointer reader = InputReaderType::New(); - reader->SetFileName( m_InputFileName); - reader->Update(); - typename InputImageType::Pointer input= reader->GetOutput(); - - // // globals for avi - // unsigned int number=0; - - - if(m_Verbose) - { - std::cout<GetLargestPossibleRegion()); + typename InternalImageType::IndexType index; + while(!it.IsAtEnd()) { + index=it.GetIndex(); + for (unsigned int i=0; iGetLargestPossibleRegion().GetIndex()[i] + || index[i]==(unsigned int)air_low->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air_low->GetLargestPossibleRegion().GetSize()[i]-1) + it.Set(0); + ++it; + } + } - //------------------------------------------------------------------------------- - // Air - //------------------------------------------------------------------------------- - typename InternalImageType::Pointer air = this->GetAirImage(input); - - //------------------------------------------------------------------------------- - // Bones - //------------------------------------------------------------------------------- - typename InternalImageType::Pointer bones = this->GetBonesImage(input); - - //-------------------------------------------------------------------------------- - // Lungs - //------------------------------------------------------------------------------- - typename InternalImageType::Pointer lungs = this->GetLungsImage(input); - - //---------------------------------------------------------------------------------------------------- - // Write feature images - //---------------------------------------------------------------------------------------------------- - if(m_ArgsInfo.writeFeature_given) - { - typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New(); - setBackgroundFilter->SetInput(air); - setBackgroundFilter->SetInput2(bones); - setBackgroundFilter->SetMaskValue(0); - setBackgroundFilter->SetOutsideValue(2); - setBackgroundFilter->Update(); - typename InternalImageType::Pointer bones_air =setBackgroundFilter->GetOutput(); - - typename SetBackgroundFilterType::Pointer setBackgroundFilter2 = SetBackgroundFilterType::New(); - setBackgroundFilter2->SetInput(bones_air); - setBackgroundFilter2->SetInput2(lungs); - setBackgroundFilter2->SetMaskValue(0); - setBackgroundFilter2->SetOutsideValue(3); - setBackgroundFilter2->Update(); - typename InternalImageType::Pointer lungs_bones_air =setBackgroundFilter2->GetOutput(); - - typename CastImageFilterType::Pointer caster=CastImageFilterType::New(); - caster->SetInput(lungs_bones_air); - writeImage(caster->GetOutput(), m_ArgsInfo.writeFeature_arg, m_Verbose); - } + //--------------------------------- + // Center + //--------------------------------- + typename itk::Vector center; + typedef itk::ImageMomentsCalculator MomentsCalculatorType; + typename MomentsCalculatorType::Pointer momentsCalculator=MomentsCalculatorType::New(); + momentsCalculator->SetImage(air); + momentsCalculator->Compute(); + center=momentsCalculator->GetCenterOfGravity(); + if (m_Verbose) { + std::cout<<"The center of gravity of the patient body is located at (mm) "<(center,bones_low,lungs_low); + + //---------------------------------------------------------------------------------------------------- + // Grow ellips + //---------------------------------------------------------------------------------------------------- + typename InternalImageType::Pointer grownEllips; + if (m_ArgsInfo.grownEllips_given || m_ArgsInfo.filledRibs_given) { + if (m_ArgsInfo.grownEllips_given) { + + typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New(); + featureReader->SetFileName(m_ArgsInfo.grownEllips_arg); + featureReader->Update(); + grownEllips=featureReader->GetOutput(); + } + } else { + + if(m_Verbose) { + std::cout<(bones); - typename InternalImageType::Pointer lungs_low =Resample(lungs); - typename InternalImageType::Pointer air_low =Resample(air); - //--------------------------------- - // Pad + // Detect abdomen //--------------------------------- - if(m_ArgsInfo.pad_flag) - { - typedef itk::ImageRegionIteratorWithIndex IteratorType; - IteratorType it(air_low, air_low->GetLargestPossibleRegion()); - typename InternalImageType::IndexType index; - while(!it.IsAtEnd()) - { - index=it.GetIndex(); - for (unsigned int i=0;iGetLargestPossibleRegion().GetIndex()[i] - || index[i]==(unsigned int)air_low->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air_low->GetLargestPossibleRegion().GetSize()[i]-1) - it.Set(0); - ++it; - } + typename InternalImageType::PointType dPoint; + if (m_ArgsInfo.detectionPoint_given) { + for (unsigned int i=0; iGetLargestPossibleRegion(); + typename InternalImageType::SizeType searchRegionSize = searchRegion.GetSize(); + typename InternalImageType::IndexType searchRegionIndex = searchRegion.GetIndex(); + searchRegionIndex[2]+=searchRegionSize[2]/2; + searchRegionSize[2]=1; + searchRegion.SetSize(searchRegionSize); + searchRegion.SetIndex(searchRegionIndex); + IteratorType itAbdomen(air, searchRegion); + itAbdomen.GoToBegin(); + + typename InternalImageType::PointType aPoint; + typename InternalImageType::IndexType aIndex; + + if (m_Verbose) std::cout<<"Detecting the abdomen..."<TransformIndexToPhysicalPoint(aIndex,aPoint); + if (m_Verbose) std::cout<<"Detected the abdomen at "< InputReaderType; + typename InputReaderType::Pointer reader = InputReaderType::New(); + reader->SetFileName( m_ArgsInfo.detectionPairs_arg[i]); + reader->Update(); + typename InputImageType::Pointer additional= reader->GetOutput(); + if (m_Verbose) std::cout<<"Extracting the air from additional image "<< i<<"..."<SetInput(additional); + binarizeFilter->SetLowerThreshold(static_cast(m_ArgsInfo.lowerThresholdAir_arg)); + binarizeFilter->SetUpperThreshold(static_cast(m_ArgsInfo.upperThresholdAir_arg)); + if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<SetInput(binarizeFilter->GetOutput()); + connectFilter->SetBackgroundValue(0); + connectFilter->SetFullyConnected(false); + if (m_Verbose) std::cout<<"Labeling the connected components..."<SetInput(connectFilter->GetOutput()); + if (m_Verbose) std::cout<<"Sorting the labels..."<SetInput(relabelFilter->GetOutput()); + thresholdFilter->SetUpper(1); + if (m_Verbose) std::cout<<"Keeping the first label..."<Update(); + airAdd=thresholdFilter->GetOutput(); + + + //--------------------------------- + // Invert + //--------------------------------- + typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New(); + inversionFilter->SetInput(airAdd); + inversionFilter->SetLowerThreshold(0); + inversionFilter->SetUpperThreshold(0); + inversionFilter->SetInsideValue(1); + inversionFilter->Update(); + + //--------------------------------- + // Mirror + //--------------------------------- + typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New(); + mirrorPadImageFilter->SetInput(inversionFilter->GetOutput()); + typename InternalImageType::SizeType padSize; + padSize.Fill(0); + padSize[2]=input->GetLargestPossibleRegion().GetSize()[2]; + mirrorPadImageFilter->SetPadLowerBound(padSize); + if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<Update(); + airAdd=mirrorPadImageFilter->GetOutput(); + + //--------------------------------- + // Detect abdomen + //--------------------------------- + typename InternalImageType::RegionType searchRegion=airAdd->GetLargestPossibleRegion(); + typename InternalImageType::SizeType searchRegionSize = searchRegion.GetSize(); + typename InternalImageType::IndexType searchRegionIndex = searchRegion.GetIndex(); + searchRegionIndex[2]+=searchRegionSize[2]/2; + searchRegionSize[2]=1; + searchRegion.SetSize(searchRegionSize); + searchRegion.SetIndex(searchRegionIndex); + IteratorType itAbdomen(airAdd, searchRegion); + itAbdomen.GoToBegin(); + + typename InternalImageType::PointType additionalPoint; + typename InternalImageType::IndexType aIndex; + + if (m_Verbose) std::cout<<"Detecting the abdomen..."<TransformIndexToPhysicalPoint(aIndex,additionalPoint); + if (m_Verbose) std::cout<<"Detected the abdomen in the additional image at "< center; - typedef itk::ImageMomentsCalculator MomentsCalculatorType; - typename MomentsCalculatorType::Pointer momentsCalculator=MomentsCalculatorType::New(); - momentsCalculator->SetImage(air); - momentsCalculator->Compute(); - center=momentsCalculator->GetCenterOfGravity(); - if (m_Verbose) - { - std::cout<<"The center of gravity of the patient body is located at (mm) "<TransformPhysicalPointToIndex(dPoint, dIndex)) { + typename InternalImageType::SizeType padSize; + padSize.Fill(0); + padSize[1]=abs(dIndex[1])+1; + if (m_Verbose) std::cout<<"Padding the images with padding size "<SetInput(bones_low); + padBonesFilter->SetPadLowerBound(padSize); + padBonesFilter->Update(); + padded_bones_low=padBonesFilter->GetOutput(); + + typename MirrorPadImageFilterType::Pointer padEllipsFilter=MirrorPadImageFilterType::New(); + padEllipsFilter->SetInput(m_Ellips); + padEllipsFilter->SetPadLowerBound(padSize); + padEllipsFilter->Update(); + padded_ellips=padEllipsFilter->GetOutput(); + } + + else { + padded_bones_low=bones_low; + padded_ellips=m_Ellips; + } - //---------------------------------------------------------------------------------------------------- - // Ellipsoide initialization - //---------------------------------------------------------------------------------------------------- - typename InternalImageType::Pointer m_Ellips= InitializeEllips(center,bones_low); - - //---------------------------------------------------------------------------------------------------- - // Grow ellips - //---------------------------------------------------------------------------------------------------- - typename InternalImageType::Pointer grownEllips; - if (m_ArgsInfo.grownEllips_given || m_ArgsInfo.filledRibs_given) - { - if (m_ArgsInfo.grownEllips_given) - { - - typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New(); - featureReader->SetFileName(m_ArgsInfo.grownEllips_arg); - featureReader->Update(); - grownEllips=featureReader->GetOutput(); - } - } - else - { - - if(m_Verbose) - { - std::cout<GetLargestPossibleRegion(); - typename InternalImageType::SizeType searchRegionSize = searchRegion.GetSize(); - typename InternalImageType::IndexType searchRegionIndex = searchRegion.GetIndex(); - searchRegionIndex[2]+=searchRegionSize[2]/2; - searchRegionSize[2]=1; - searchRegion.SetSize(searchRegionSize); - searchRegion.SetIndex(searchRegionIndex); - IteratorType itAbdomen(air, searchRegion); - itAbdomen.GoToBegin(); - - typename InternalImageType::PointType aPoint; - typename InternalImageType::IndexType aIndex; - - if (m_Verbose) std::cout<<"Detecting the abdomen..."<TransformIndexToPhysicalPoint(aIndex,aPoint); - if (m_Verbose) std::cout<<"Detected the abdomen at "< InputReaderType; - typename InputReaderType::Pointer reader = InputReaderType::New(); - reader->SetFileName( m_ArgsInfo.detectionPairs_arg[i]); - reader->Update(); - typename InputImageType::Pointer additional= reader->GetOutput(); - if (m_Verbose) std::cout<<"Extracting the air from additional image "<< i<<"..."<SetInput(additional); - binarizeFilter->SetLowerThreshold(static_cast(m_ArgsInfo.lowerThresholdAir_arg)); - binarizeFilter->SetUpperThreshold(static_cast(m_ArgsInfo.upperThresholdAir_arg)); - if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<SetInput(binarizeFilter->GetOutput()); - connectFilter->SetBackgroundValue(0); - connectFilter->SetFullyConnected(false); - if (m_Verbose) std::cout<<"Labeling the connected components..."<SetInput(connectFilter->GetOutput()); - if (m_Verbose) std::cout<<"Sorting the labels..."<SetInput(relabelFilter->GetOutput()); - thresholdFilter->SetUpper(1); - if (m_Verbose) std::cout<<"Keeping the first label..."<Update(); - airAdd=thresholdFilter->GetOutput(); - - - //--------------------------------- - // Invert - //--------------------------------- - typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New(); - inversionFilter->SetInput(airAdd); - inversionFilter->SetLowerThreshold(0); - inversionFilter->SetUpperThreshold(0); - inversionFilter->SetInsideValue(1); - inversionFilter->Update(); - - //--------------------------------- - // Mirror - //--------------------------------- - typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New(); - mirrorPadImageFilter->SetInput(inversionFilter->GetOutput()); - typename InternalImageType::SizeType padSize; - padSize.Fill(0); - padSize[2]=input->GetLargestPossibleRegion().GetSize()[2]; - mirrorPadImageFilter->SetPadLowerBound(padSize); - if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<Update(); - airAdd=mirrorPadImageFilter->GetOutput(); - - //--------------------------------- - // Detect abdomen - //--------------------------------- - typename InternalImageType::RegionType searchRegion=airAdd->GetLargestPossibleRegion(); - typename InternalImageType::SizeType searchRegionSize = searchRegion.GetSize(); - typename InternalImageType::IndexType searchRegionIndex = searchRegion.GetIndex(); - searchRegionIndex[2]+=searchRegionSize[2]/2; - searchRegionSize[2]=1; - searchRegion.SetSize(searchRegionSize); - searchRegion.SetIndex(searchRegionIndex); - IteratorType itAbdomen(airAdd, searchRegion); - itAbdomen.GoToBegin(); - - typename InternalImageType::PointType additionalPoint; - typename InternalImageType::IndexType aIndex; - - if (m_Verbose) std::cout<<"Detecting the abdomen..."<TransformIndexToPhysicalPoint(aIndex,additionalPoint); - if (m_Verbose) std::cout<<"Detected the abdomen in the additional image at "<TransformPhysicalPointToIndex(dPoint, dIndex)) - { - typename InternalImageType::SizeType padSize; - padSize.Fill(0); - padSize[1]=abs(dIndex[1])+1; - if (m_Verbose) std::cout<<"Padding the images with padding size "<SetInput(bones_low); - padBonesFilter->SetPadLowerBound(padSize); - padBonesFilter->Update(); - padded_bones_low=padBonesFilter->GetOutput(); - - typename MirrorPadImageFilterType::Pointer padEllipsFilter=MirrorPadImageFilterType::New(); - padEllipsFilter->SetInput(m_Ellips); - padEllipsFilter->SetPadLowerBound(padSize); - padEllipsFilter->Update(); - padded_ellips=padEllipsFilter->GetOutput(); - } - - else - { - padded_bones_low=bones_low; - padded_ellips=m_Ellips; - } - - - //--------------------------------- - // Calculate distance map - //--------------------------------- - typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New(); - distanceMapImageFilter->SetInput(padded_ellips); - distanceMapImageFilter->SetInsideIsPositive(false); - if (m_Verbose) std::cout<<"Calculating the distance map..."<Update(); - - //--------------------------------- - // Grow while monitoring detection point - //--------------------------------- - typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New(); - levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() ); - levelSetFilter->SetFeatureImage( padded_bones_low ); - levelSetFilter->SetPropagationScaling( 1 ); - levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve1_arg ); - levelSetFilter->SetAdvectionScaling( 0 ); - levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS1_arg ); - levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter1_arg ); - levelSetFilter->SetUseImageSpacing(true); - - // //--------------------------------- - // // Script for making movie - // //--------------------------------- - // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh "; - - - if (m_Verbose) std::cout<<"Starting level set segmentation..."<Update(); - totalNumberOfIterations+=levelSetFilter->GetElapsedIterations(); - - if ( levelSetFilter->GetOutput()->GetPixel(dIndex) < 0 ) - { - if (m_Verbose) std::cout<<"Detection point reached!"<SetInput(levelSetFilter->GetOutput()); - if(m_ArgsInfo.monitor_given) writeImage(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose); - - // // script - // std::ostringstream number_str; - // number_str << number; - // std::string param = number_str.str(); - // system((script+param).c_str()); - // number+=1; - - } - if ( (totalNumberOfIterations> 10000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS1_arg) ) break; - } - - // Output values - if (m_Verbose) std::cout<<"Level set segmentation stopped after "<GetMaximumRMSError() << std::endl; - std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl; - - // Threshold - typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New(); - thresholder->SetUpperThreshold( 0.0 ); - thresholder->SetOutsideValue( 0 ); - thresholder->SetInsideValue( 1 ); - thresholder->SetInput( levelSetFilter->GetOutput() ); - if (m_Verbose) std::cout<<"Thresholding the output level set..."<Update(); - grownEllips=thresholder->GetOutput(); - } //--------------------------------- - // Write the grown ellips + // Calculate distance map //--------------------------------- - if (m_ArgsInfo.writeGrownEllips_given) - { - typename CastImageFilterType::Pointer caster=CastImageFilterType::New(); - caster->SetInput(grownEllips); - writeImage(caster->GetOutput(), m_ArgsInfo.writeGrownEllips_arg, m_Verbose); - } - - - //---------------------------------------------------------------------------------------------------- - // Grow inside ribs - //---------------------------------------------------------------------------------------------------- - typename InternalImageType::Pointer filledRibs; - if (m_ArgsInfo.filledRibs_given) - { - typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New(); - featureReader->SetFileName(m_ArgsInfo.filledRibs_arg); - if (m_Verbose) std::cout<<"Reading filled ribs mask..."<Update(); - filledRibs=featureReader->GetOutput(); - } - else - { - if(m_Verbose) - { - std::cout<SetInput(air_low); - setBackgroundFilter->SetInput2(bones_low); - setBackgroundFilter->SetMaskValue(0); - setBackgroundFilter->SetOutsideValue(0); - setBackgroundFilter->Update(); - typename InternalImageType::Pointer bones_air_low =setBackgroundFilter->GetOutput(); - - //--------------------------------- - // Resampling previous solution - //--------------------------------- - if (m_Verbose) std::cout<<"Resampling previous solution..."<TransformIndexToPhysicalPoint(bones_low->GetLargestPossibleRegion().GetIndex(), origin); - resampler->SetOutputOrigin(origin); - resampler->SetOutputSpacing(bones_low->GetSpacing()); - typename InternalImageType::SizeType size_low= bones_low->GetLargestPossibleRegion().GetSize(); - resampler->SetSize(size_low); - typename InterpolatorType::Pointer interpolator=InterpolatorType::New(); - resampler->SetInterpolator(interpolator); - resampler->SetInput(grownEllips); - resampler->Update(); - typename InternalImageType::Pointer grownEllips =resampler->GetOutput(); - - - //--------------------------------- - // Calculate Distance Map - //--------------------------------- - typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New(); - distanceMapImageFilter->SetInput(grownEllips); - distanceMapImageFilter->SetInsideIsPositive(false); - if (m_Verbose) std::cout<<"Calculating distance map..."<Update(); - - //--------------------------------- - // Grow while monitoring lung volume coverage - //--------------------------------- - typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New(); - levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() ); - levelSetFilter->SetFeatureImage( bones_air_low ); - levelSetFilter->SetPropagationScaling( 1 ); - levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve2_arg ); - levelSetFilter->SetAdvectionScaling( 0 ); - levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS2_arg ); - levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter2_arg ); - levelSetFilter->SetUseImageSpacing(true); - - //--------------------------------- - // Invert lung mask - //--------------------------------- - typename InversionFilterType::Pointer invertor= InversionFilterType::New(); - invertor->SetInput(lungs_low); - invertor->SetLowerThreshold(0); - invertor->SetUpperThreshold(0); - invertor->SetInsideValue(1); - invertor->Update(); - typename InternalImageType::Pointer lungs_low_inv=invertor->GetOutput(); - - //--------------------------------- - // Calculate lung volume - //--------------------------------- - typename LabelStatisticsImageFilterType::Pointer statisticsFilter= LabelStatisticsImageFilterType::New(); - statisticsFilter->SetInput(lungs_low_inv); - statisticsFilter->SetLabelInput(lungs_low); - statisticsFilter->Update(); - unsigned int volume=statisticsFilter->GetSum(0); - - // Prepare threshold - typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New(); - thresholder->SetUpperThreshold( 0.0 ); - thresholder->SetOutsideValue( 0 ); - thresholder->SetInsideValue( 1 ); - - - // Start monitoring - unsigned int totalNumberOfIterations=0; - unsigned int coverage=0; - - - // //--------------------------------- - // // Script for making movie - // //--------------------------------- - // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh "; - - while (true) - { - levelSetFilter->Update(); - totalNumberOfIterations+=levelSetFilter->GetElapsedIterations(); - - thresholder->SetInput( levelSetFilter->GetOutput() ); - thresholder->Update(); - statisticsFilter->SetInput(thresholder->GetOutput()); - statisticsFilter->Update(); - coverage=statisticsFilter->GetSum(0); - - // Compare the volumes - if ( coverage >= m_ArgsInfo.fillingLevel_arg * volume/100 ) - { - if (m_Verbose) std::cout<<"Lungs filled for "<< m_ArgsInfo.fillingLevel_arg<<"% !"<SetInput(levelSetFilter->GetOutput()); - if(m_ArgsInfo.monitor_given) writeImage(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose); - - // // script - // std::ostringstream number_str; - // number_str << number; - // std::string param = number_str.str(); - // system((script+param).c_str()); - // number+=1; - - } - if ( (totalNumberOfIterations> 30000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS2_arg) ) break; - } - - // Output values - if (m_Verbose) std::cout<<"Level set segmentation stopped after "<GetMaximumRMSError() << std::endl; - std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl; - - // Threshold - thresholder->SetInput( levelSetFilter->GetOutput() ); - thresholder->Update(); - filledRibs=thresholder->GetOutput(); - // writeImage(filledRibs,"/home/jef/tmp/filled_ribs_before.mhd", m_Verbose); - - } + typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New(); + distanceMapImageFilter->SetInput(padded_ellips); + distanceMapImageFilter->SetInsideIsPositive(false); + if (m_Verbose) std::cout<<"Calculating the distance map..."<Update(); + if (m_ArgsInfo.writeDistMap_given) { + writeImage(distanceMapImageFilter->GetOutput(), m_ArgsInfo.writeDistMap_arg, m_Verbose); + + } //--------------------------------- - // Write the filled ribs + // Grow while monitoring detection point //--------------------------------- - if (m_ArgsInfo.writeFilledRibs_given) - { - typename CastImageFilterType::Pointer caster=CastImageFilterType::New(); - caster->SetInput(filledRibs); - writeImage(caster->GetOutput(), m_ArgsInfo.writeFilledRibs_arg, m_Verbose); - } + typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New(); + levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() ); + levelSetFilter->SetFeatureImage( padded_bones_low ); + levelSetFilter->SetPropagationScaling( 1 ); + levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve1_arg ); + levelSetFilter->SetAdvectionScaling( 0 ); + levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS1_arg ); + levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter1_arg ); + levelSetFilter->SetUseImageSpacing(true); + + // //--------------------------------- + // // Script for making movie + // //--------------------------------- + // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh "; + + + if (m_Verbose) std::cout<<"Starting level set segmentation..."<Update(); + totalNumberOfIterations+=levelSetFilter->GetElapsedIterations(); + + if ( levelSetFilter->GetOutput()->GetPixel(dIndex) < 0 ) { + if (m_Verbose) std::cout<<"Detection point reached!"<SetInput(levelSetFilter->GetOutput()); + if(m_ArgsInfo.monitor_given) writeImage(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose); + + // // script + // std::ostringstream number_str; + // number_str << number; + // std::string param = number_str.str(); + // system((script+param).c_str()); + // number+=1; - - //---------------------------------------------------------------------------------------------------- - // Collapse to the lungs - //---------------------------------------------------------------------------------------------------- - if(m_Verbose) - { - std::cout< 10000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS1_arg) ) break; + } + + // Output values + if (m_Verbose) std::cout<<"Level set segmentation stopped after "<GetMaximumRMSError() << std::endl; + std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl; + + // Threshold + typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New(); + thresholder->SetUpperThreshold( 0.0 ); + thresholder->SetOutsideValue( 0 ); + thresholder->SetInsideValue( 1 ); + thresholder->SetInput( levelSetFilter->GetOutput() ); + if (m_Verbose) std::cout<<"Thresholding the output level set..."<Update(); + grownEllips=thresholder->GetOutput(); + } + + //--------------------------------- + // Write the grown ellips + //--------------------------------- + if (m_ArgsInfo.writeGrownEllips_given) { + typename CastImageFilterType::Pointer caster=CastImageFilterType::New(); + caster->SetInput(grownEllips); + writeImage(caster->GetOutput(), m_ArgsInfo.writeGrownEllips_arg, m_Verbose); + } + + //---------------------------------------------------------------------------------------------------- + // Grow inside ribs + //---------------------------------------------------------------------------------------------------- + typename InternalImageType::Pointer filledRibs; + if (m_ArgsInfo.filledRibs_given) { + typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New(); + featureReader->SetFileName(m_ArgsInfo.filledRibs_arg); + if (m_Verbose) std::cout<<"Reading filled ribs mask..."<Update(); + filledRibs=featureReader->GetOutput(); + } else { + if(m_Verbose) { + std::cout<SetInput(air); - setBackgroundFilter->SetInput2(bones); + typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New(); + setBackgroundFilter->SetInput(air_low); + setBackgroundFilter->SetInput2(bones_low); setBackgroundFilter->SetMaskValue(0); setBackgroundFilter->SetOutsideValue(0); setBackgroundFilter->Update(); - typename InternalImageType::Pointer bones_air =setBackgroundFilter->GetOutput(); - - typename SetBackgroundFilterType::Pointer setBackgroundFilter2 = SetBackgroundFilterType::New(); - setBackgroundFilter2->SetInput(bones_air); - setBackgroundFilter2->SetInput2(lungs); - setBackgroundFilter2->SetMaskValue(0); - setBackgroundFilter2->SetOutsideValue(0); - setBackgroundFilter2->Update(); - typename InternalImageType::Pointer lungs_bones_air =setBackgroundFilter2->GetOutput(); + typename InternalImageType::Pointer bones_air_low =setBackgroundFilter->GetOutput(); //--------------------------------- - // Prepare previous solution + // Resampling previous solution //--------------------------------- if (m_Verbose) std::cout<<"Resampling previous solution..."< InterpolatorType; + typename InternalImageType::PointType origin; + bones_low->TransformIndexToPhysicalPoint(bones_low->GetLargestPossibleRegion().GetIndex(), origin); + resampler->SetOutputOrigin(origin); + resampler->SetOutputSpacing(bones_low->GetSpacing()); + typename InternalImageType::SizeType size_low= bones_low->GetLargestPossibleRegion().GetSize(); + resampler->SetSize(size_low); typename InterpolatorType::Pointer interpolator=InterpolatorType::New(); - resampler->SetOutputStartIndex(bones->GetLargestPossibleRegion().GetIndex()); - resampler->SetInput(filledRibs); resampler->SetInterpolator(interpolator); - resampler->SetOutputParametersFromImage(bones); + resampler->SetInput(grownEllips); resampler->Update(); - filledRibs =resampler->GetOutput(); - - if (m_Verbose) std::cout<<"Setting lungs to 1..."<SetInput(filledRibs); - setBackgroundFilter3->SetInput2(lungs); - setBackgroundFilter3->SetMaskValue(0); - setBackgroundFilter3->SetOutsideValue(1); - setBackgroundFilter3->Update(); - filledRibs=setBackgroundFilter3->GetOutput(); - - if (m_Verbose) std::cout<<"Removing overlap with bones..."<SetInput(filledRibs); - setBackgroundFilter4->SetInput2(bones); - setBackgroundFilter4->SetMaskValue(0); - setBackgroundFilter4->SetOutsideValue(0); - setBackgroundFilter4->Update(); - filledRibs =setBackgroundFilter4->GetOutput(); - // writeImage(filledRibs,"/home/jef/tmp/filledRibs_pp.mhd"); + typename InternalImageType::Pointer grownEllips =resampler->GetOutput(); + + //--------------------------------- // Calculate Distance Map //--------------------------------- - // typedef itk::ApproximateSignedDistanceMapImageFilter DistanceMapImageFilterType2; typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New(); - distanceMapImageFilter->SetInput(filledRibs); + distanceMapImageFilter->SetInput(grownEllips); distanceMapImageFilter->SetInsideIsPositive(false); - // distanceMapImageFilter->SetInsideValue(0); -// distanceMapImageFilter->SetOutsideValue(1); if (m_Verbose) std::cout<<"Calculating distance map..."<Update(); - + if (m_ArgsInfo.writeDistMap_given) { + writeImage(distanceMapImageFilter->GetOutput(), m_ArgsInfo.writeDistMap_arg, m_Verbose); + + } + //--------------------------------- - // Collapse + // Grow while monitoring lung volume coverage //--------------------------------- typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New(); levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() ); - levelSetFilter->SetFeatureImage( lungs_bones_air ); - levelSetFilter->SetPropagationScaling(m_ArgsInfo.prop3_arg); - levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve3_arg ); + levelSetFilter->SetFeatureImage( bones_air_low ); + levelSetFilter->SetPropagationScaling( 1 ); + levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve2_arg ); levelSetFilter->SetAdvectionScaling( 0 ); - levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS3_arg ); - levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter3_arg ); + levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS2_arg ); + levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter2_arg ); levelSetFilter->SetUseImageSpacing(true); - - // //script - // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh "; - + + //--------------------------------- + // Invert lung mask + //--------------------------------- + typename InversionFilterType::Pointer invertor= InversionFilterType::New(); + invertor->SetInput(lungs_low); + invertor->SetLowerThreshold(0); + invertor->SetUpperThreshold(0); + invertor->SetInsideValue(1); + invertor->Update(); + typename InternalImageType::Pointer lungs_low_inv=invertor->GetOutput(); + + //--------------------------------- + // Calculate lung volume + //--------------------------------- + typename LabelStatisticsImageFilterType::Pointer statisticsFilter= LabelStatisticsImageFilterType::New(); + statisticsFilter->SetInput(lungs_low_inv); + statisticsFilter->SetLabelInput(lungs_low); + statisticsFilter->Update(); + unsigned int volume=statisticsFilter->GetSum(0); + + // Prepare threshold + typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New(); + thresholder->SetUpperThreshold( 0.0 ); + thresholder->SetOutsideValue( 0 ); + thresholder->SetInsideValue( 1 ); + + // Start monitoring - if (m_Verbose) std::cout<<"Starting the levelset..."<Update(); - - // monitor state - totalNumberOfIterations+=levelSetFilter->GetElapsedIterations(); - levelSetFilter->SetInput(levelSetFilter->GetOutput()); - if(m_ArgsInfo.monitor_given) writeImage(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg); - std::cout << "After "<GetRMSChange() <<"..."<< std::endl; - - // // script - // std::ostringstream number_str; - // number_str << number; - // std::string param = number_str.str(); - // system((script+param).c_str()); - // number+=1; - - // stopping condition - if ( (totalNumberOfIterations>= (unsigned int) m_ArgsInfo.maxIter3_arg) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS3_arg) ) break; - levelSetFilter->SetNumberOfIterations( std::min ((unsigned int) m_ArgsInfo.iter3_arg, (unsigned int) m_ArgsInfo.maxIter3_arg-totalNumberOfIterations ) ); + unsigned int coverage=0; + + + // //--------------------------------- + // // Script for making movie + // //--------------------------------- + // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh "; + + while (true) { + levelSetFilter->Update(); + totalNumberOfIterations+=levelSetFilter->GetElapsedIterations(); + + thresholder->SetInput( levelSetFilter->GetOutput() ); + thresholder->Update(); + statisticsFilter->SetInput(thresholder->GetOutput()); + statisticsFilter->Update(); + coverage=statisticsFilter->GetSum(0); + + // Compare the volumes + if ( coverage >= m_ArgsInfo.fillingLevel_arg * volume/100 ) { + if (m_Verbose) std::cout<<"Lungs filled for "<< m_ArgsInfo.fillingLevel_arg<<"% !"<SetInput(levelSetFilter->GetOutput()); + if(m_ArgsInfo.monitor_given) writeImage(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose); + + // // script + // std::ostringstream number_str; + // number_str << number; + // std::string param = number_str.str(); + // system((script+param).c_str()); + // number+=1; + } - + if ( (totalNumberOfIterations> 30000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS2_arg) ) break; + } + // Output values - if (m_Verbose) - { - std::cout<<"Level set segmentation stopped after "<GetMaximumRMSError() << std::endl; - std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl; - } + if (m_Verbose) std::cout<<"Level set segmentation stopped after "<GetMaximumRMSError() << std::endl; + std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl; // Threshold - typedef itk::BinaryThresholdImageFilter< LevelSetImageType,InternalImageType> LevelSetBinarizeFilterType; - typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New(); - thresholder->SetUpperThreshold( 0.0 ); - thresholder->SetOutsideValue( 0 ); - thresholder->SetInsideValue( 1 ); thresholder->SetInput( levelSetFilter->GetOutput() ); thresholder->Update(); - typename InternalImageType::Pointer output = thresholder->GetOutput(); + filledRibs=thresholder->GetOutput(); + // writeImage(filledRibs,"/home/jef/tmp/filled_ribs_before.mhd", m_Verbose); + } + + //--------------------------------- + // Write the filled ribs + //--------------------------------- + if (m_ArgsInfo.writeFilledRibs_given) { + typename CastImageFilterType::Pointer caster=CastImageFilterType::New(); + caster->SetInput(filledRibs); + writeImage(caster->GetOutput(), m_ArgsInfo.writeFilledRibs_arg, m_Verbose); + } + + + //---------------------------------------------------------------------------------------------------- + // Collapse to the lungs + //---------------------------------------------------------------------------------------------------- + if(m_Verbose) { + std::cout<SetInput(air); + setBackgroundFilter->SetInput2(bones); + setBackgroundFilter->SetMaskValue(0); + setBackgroundFilter->SetOutsideValue(0); + setBackgroundFilter->Update(); + typename InternalImageType::Pointer bones_air =setBackgroundFilter->GetOutput(); + + typename SetBackgroundFilterType::Pointer setBackgroundFilter2 = SetBackgroundFilterType::New(); + setBackgroundFilter2->SetInput(bones_air); + setBackgroundFilter2->SetInput2(lungs); + setBackgroundFilter2->SetMaskValue(0); + setBackgroundFilter2->SetOutsideValue(0); + setBackgroundFilter2->Update(); + typename InternalImageType::Pointer lungs_bones_air =setBackgroundFilter2->GetOutput(); + + //--------------------------------- + // Prepare previous solution + //--------------------------------- + if (m_Verbose) std::cout<<"Resampling previous solution..."< InterpolatorType; + typename InterpolatorType::Pointer interpolator=InterpolatorType::New(); + resampler->SetOutputStartIndex(bones->GetLargestPossibleRegion().GetIndex()); + resampler->SetInput(filledRibs); + resampler->SetInterpolator(interpolator); + resampler->SetOutputParametersFromImage(bones); + resampler->Update(); + filledRibs =resampler->GetOutput(); + + if (m_Verbose) std::cout<<"Setting lungs to 1..."<SetInput(filledRibs); + setBackgroundFilter3->SetInput2(lungs); + setBackgroundFilter3->SetMaskValue(0); + setBackgroundFilter3->SetOutsideValue(1); + setBackgroundFilter3->Update(); + filledRibs=setBackgroundFilter3->GetOutput(); + + if (m_Verbose) std::cout<<"Removing overlap with bones..."<SetInput(filledRibs); + setBackgroundFilter4->SetInput2(bones); + setBackgroundFilter4->SetMaskValue(0); + setBackgroundFilter4->SetOutsideValue(0); + setBackgroundFilter4->Update(); + filledRibs =setBackgroundFilter4->GetOutput(); + // writeImage(filledRibs,"/home/jef/tmp/filledRibs_pp.mhd"); + //--------------------------------- + // Calculate Distance Map + //--------------------------------- + // typedef itk::ApproximateSignedDistanceMapImageFilter DistanceMapImageFilterType2; + typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New(); + distanceMapImageFilter->SetInput(filledRibs); + distanceMapImageFilter->SetInsideIsPositive(false); + // distanceMapImageFilter->SetInsideValue(0); +// distanceMapImageFilter->SetOutsideValue(1); + if (m_Verbose) std::cout<<"Calculating distance map..."<Update(); + + //--------------------------------- + // Collapse + //--------------------------------- + typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New(); + levelSetFilter->SetInput( distanceMapImageFilter->GetOutput() ); + levelSetFilter->SetFeatureImage( lungs_bones_air ); + levelSetFilter->SetPropagationScaling(m_ArgsInfo.prop3_arg); + levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve3_arg ); + levelSetFilter->SetAdvectionScaling( 0 ); + levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS3_arg ); + levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter3_arg ); + levelSetFilter->SetUseImageSpacing(true); + + // //script + // std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh "; + + // Start monitoring + if (m_Verbose) std::cout<<"Starting the levelset..."<Update(); + + // monitor state + totalNumberOfIterations+=levelSetFilter->GetElapsedIterations(); + levelSetFilter->SetInput(levelSetFilter->GetOutput()); + if(m_ArgsInfo.monitor_given) writeImage(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg); + std::cout << "After "<GetRMSChange() <<"..."<< std::endl; + + // // script + // std::ostringstream number_str; + // number_str << number; + // std::string param = number_str.str(); + // system((script+param).c_str()); + // number+=1; + + // stopping condition + if ( (totalNumberOfIterations>= (unsigned int) m_ArgsInfo.maxIter3_arg) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS3_arg) ) break; + levelSetFilter->SetNumberOfIterations( std::min ((unsigned int) m_ArgsInfo.iter3_arg, (unsigned int) m_ArgsInfo.maxIter3_arg-totalNumberOfIterations ) ); + } + + // Output values + if (m_Verbose) { + std::cout<<"Level set segmentation stopped after "<GetMaximumRMSError() << std::endl; + std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl; + } + + // Threshold + typedef itk::BinaryThresholdImageFilter< LevelSetImageType,InternalImageType> LevelSetBinarizeFilterType; + typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New(); + thresholder->SetUpperThreshold( 0.0 ); + thresholder->SetOutsideValue( 0 ); + thresholder->SetInsideValue( 1 ); + thresholder->SetInput( levelSetFilter->GetOutput() ); + thresholder->Update(); + typename InternalImageType::Pointer output = thresholder->GetOutput(); + + + //---------------------------------------------------------------------------------------------------- + // Prepare the output + //---------------------------------------------------------------------------------------------------- + + //--------------------------------- + // Set Air to zero + //--------------------------------- + if (m_Verbose) std::cout<<"Removing overlap mask with air..."<SetInput(output); + setBackgroundFilter5->SetInput2(air); + setBackgroundFilter5->SetMaskValue(0); + setBackgroundFilter5->SetOutsideValue(0); + setBackgroundFilter5->Update(); + output=setBackgroundFilter5->GetOutput(); + + //--------------------------------- + // Open & close to cleanup + //--------------------------------- + if ( m_ArgsInfo.openClose_flag) { - //---------------------------------------------------------------------------------------------------- - // Prepare the output - //---------------------------------------------------------------------------------------------------- - //--------------------------------- - // Set Air to zero + // Structuring element //--------------------------------- - if (m_Verbose) std::cout<<"Removing overlap mask with air..."<SetInput(output); - setBackgroundFilter5->SetInput2(air); - setBackgroundFilter5->SetMaskValue(0); - setBackgroundFilter5->SetOutsideValue(0); - setBackgroundFilter5->Update(); - output=setBackgroundFilter5->GetOutput(); - + typedef itk::BinaryBallStructuringElement KernelType; + KernelType structuringElement; + structuringElement.SetRadius(1); + structuringElement.CreateStructuringElement(); + //--------------------------------- - // Open & close to cleanup + // Open //--------------------------------- - if ( m_ArgsInfo.openClose_flag) - { - - //--------------------------------- - // Structuring element - //--------------------------------- - typedef itk::BinaryBallStructuringElement KernelType; - KernelType structuringElement; - structuringElement.SetRadius(1); - structuringElement.CreateStructuringElement(); - - //--------------------------------- - // Open - //--------------------------------- - typedef itk::BinaryMorphologicalOpeningImageFilter OpenFilterType; - typename OpenFilterType::Pointer openFilter = OpenFilterType::New(); - openFilter->SetInput(output); - openFilter->SetBackgroundValue(0); - openFilter->SetForegroundValue(1); - openFilter->SetKernel(structuringElement); - if(m_Verbose) std::cout<<"Opening the output image..."< CloseFilterType; - typename CloseFilterType::Pointer closeFilter = CloseFilterType::New(); - closeFilter->SetInput(openFilter->GetOutput()); - closeFilter->SetSafeBorder(true); - closeFilter->SetForegroundValue(1); - closeFilter->SetKernel(structuringElement); - if(m_Verbose) std::cout<<"Closing the output image..."<Update(); - output=closeFilter->GetOutput(); + typedef itk::BinaryMorphologicalOpeningImageFilter OpenFilterType; + typename OpenFilterType::Pointer openFilter = OpenFilterType::New(); + openFilter->SetInput(output); + openFilter->SetBackgroundValue(0); + openFilter->SetForegroundValue(1); + openFilter->SetKernel(structuringElement); + if(m_Verbose) std::cout<<"Opening the output image..."< CloseFilterType; + typename CloseFilterType::Pointer closeFilter = CloseFilterType::New(); + closeFilter->SetInput(openFilter->GetOutput()); + closeFilter->SetSafeBorder(true); + closeFilter->SetForegroundValue(1); + closeFilter->SetKernel(structuringElement); + if(m_Verbose) std::cout<<"Closing the output image..."<Update(); + output=closeFilter->GetOutput(); - } - //writeImage(output,"/home/jef/tmp/mm_double.mhd"); - - // Extract the upper part - typedef itk::CropImageFilter CropImageFilterType; - typename CropImageFilterType::Pointer cropFilter=CropImageFilterType::New(); - cropFilter->SetInput(output); - typename InputImageType::SizeType lSize, uSize; - uSize.Fill(0); - lSize.Fill(0); - lSize[2]=input->GetLargestPossibleRegion().GetSize()[2]; - cropFilter->SetLowerBoundaryCropSize(lSize); - cropFilter->SetUpperBoundaryCropSize(uSize); - cropFilter->Update(); - - // Output - typename CastImageFilterType::Pointer caster=CastImageFilterType::New(); - caster->SetInput(cropFilter->GetOutput()); - writeImage(caster->GetOutput(), m_ArgsInfo.output_arg, m_Verbose); - } + //writeImage(output,"/home/jef/tmp/mm_double.mhd"); + + // Extract the upper part + typedef itk::CropImageFilter CropImageFilterType; + typename CropImageFilterType::Pointer cropFilter=CropImageFilterType::New(); + cropFilter->SetInput(output); + typename InputImageType::SizeType lSize, uSize; + uSize.Fill(0); + lSize.Fill(0); + lSize[2]=input->GetLargestPossibleRegion().GetSize()[2]; + cropFilter->SetLowerBoundaryCropSize(lSize); + cropFilter->SetUpperBoundaryCropSize(uSize); + cropFilter->Update(); + + // Output + typename CastImageFilterType::Pointer caster=CastImageFilterType::New(); + caster->SetInput(cropFilter->GetOutput()); + writeImage(caster->GetOutput(), m_ArgsInfo.output_arg, m_Verbose); + +} }//end clitk - + #endif //#define clitkMotionMaskGenericFilter_txx