- //----------------------------------------------------------------------------------------------------
- // Ellipsoide initialization
- //----------------------------------------------------------------------------------------------------
- typename InternalImageType::Pointer m_Ellips= InitializeEllips<Dimension, PixelType>(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<<std::endl;
- std::cout<<"=========================================="<<std::endl;
- std::cout<<"|| Growing ellipsoide ||"<<std::endl;
- std::cout<<"=========================================="<<std::endl;
- }
-
- //---------------------------------
- // Detect abdomen
- //---------------------------------
- typename InternalImageType::PointType dPoint;
- if (m_ArgsInfo.detectionPoint_given)
- {
- for (unsigned int i=0;i<Dimension;i++)
- dPoint[i]=m_ArgsInfo.detectionPoint_arg[i];
- }
- else
- {
- typename InternalImageType::RegionType searchRegion=air->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..."<<std::endl;
- while (!itAbdomen.IsAtEnd())
- {
- if(itAbdomen.Get()) break;
- ++itAbdomen;
- }
- aIndex=itAbdomen.GetIndex();
- air->TransformIndexToPhysicalPoint(aIndex,aPoint);
- if (m_Verbose) std::cout<<"Detected the abdomen at "<<aPoint<<".."<<std::endl;
-
-
- //---------------------------------
- // Detect abdomen in additional images?
- //---------------------------------
- if (m_ArgsInfo.detectionPairs_given)
- {
- for (unsigned int i=0; i<m_ArgsInfo.detectionPairs_given; i++)
- {
- typename InternalImageType::Pointer airAdd;
- //---------------------------------
- // Read the input
- //--------------------------------
- typedef itk::ImageFileReader<InputImageType> 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<<"..."<<std::endl;
-
- //---------------------------------
- // Binarize the image
- //---------------------------------
- typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
- binarizeFilter->SetInput(additional);
- binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdAir_arg));
- binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdAir_arg));
- if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<<m_ArgsInfo.lowerThresholdAir_arg
- <<", "<<m_ArgsInfo.upperThresholdAir_arg<<"..."<<std::endl;
-
- //---------------------------------
- // 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..."<<std::endl;
-
- //---------------------------------
- // Sort the labels according to size
- //---------------------------------
- typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
- relabelFilter->SetInput(connectFilter->GetOutput());
- if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
-
- //---------------------------------
- // Keep the label
- //---------------------------------
- typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
- thresholdFilter->SetInput(relabelFilter->GetOutput());
- thresholdFilter->SetUpper(1);
- if (m_Verbose) std::cout<<"Keeping the first label..."<<std::endl;
- thresholdFilter->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..."<<std::endl;
- mirrorPadImageFilter->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..."<<std::endl;
- while (!itAbdomen.IsAtEnd())
- {
- if(itAbdomen.Get()) break;
- ++itAbdomen;
- }
- aIndex=itAbdomen.GetIndex();
- airAdd->TransformIndexToPhysicalPoint(aIndex,additionalPoint);
- if (m_Verbose) std::cout<<"Detected the abdomen in the additional image at "<<additionalPoint<<".."<<std::endl;
-
- if(additionalPoint[1]< aPoint[1])
- {
- aPoint=additionalPoint;
- if (m_Verbose) std::cout<<"Modifying the detected abdomen to "<<aPoint<<".."<<std::endl;
-
- }
- }
- }
-
-
- // Determine the detection point
- dPoint.Fill(0.0);
- dPoint+=center;
- dPoint[1]=aPoint[1];
- if(m_ArgsInfo.offsetDetect_given==Dimension)
- for(unsigned int i=0; i <Dimension; i++)
- dPoint[i]+=m_ArgsInfo.offsetDetect_arg[i];
- else
- dPoint[1]+=-10;
-
- }
- if (m_Verbose) std::cout<<"Setting the detection point to "<<dPoint<<".."<<std::endl;
-
-
- //---------------------------------
- // Pad the rib image and ellips image
- //---------------------------------
- typename InternalImageType::Pointer padded_ellips;
- typename InternalImageType::Pointer padded_bones_low;
-
- // If detection point not inside the image: pad
- typename InternalImageType::IndexType dIndex;
- if (!bones_low->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 "<<padSize<<" to include the detection point..."<<std::endl;
-
- typename MirrorPadImageFilterType::Pointer padBonesFilter=MirrorPadImageFilterType::New();
- padBonesFilter->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..."<<std::endl;
- distanceMapImageFilter->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..."<<std::endl;
- unsigned int totalNumberOfIterations=0;
- while (true)
- {
- levelSetFilter->Update();
- totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
-
- if ( levelSetFilter->GetOutput()->GetPixel(dIndex) < 0 )
- {
- if (m_Verbose) std::cout<<"Detection point reached!"<<std::endl;
- break;
- }
- else
- {
- if (m_Verbose) std::cout<<"Detection point not reached after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
- levelSetFilter->SetInput(levelSetFilter->GetOutput());
- if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(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 "<<totalNumberOfIterations<<" iterations..."<<std::endl;
- std::cout << "Max. RMS error: " << levelSetFilter->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..."<<std::endl;
- thresholder->Update();
- grownEllips=thresholder->GetOutput();
- }