[ITK-users] segmentation a 3D image takes a long time , any help to speed it ?

[alaamegawer]

Hi All 

I'm trying to segment the lung from 3D Ct image using this algorithm 
1- determine the optimal threshold 
2- using connected region to delete the unwanted region 
3- filling holes by using closing operation 

but the algorithm is very slow it takes more than 10 min  

so please i want to ask if any one can guide me to speed it .

and what's the difference between segment the lung slice per slice or as a
whole 

image information

NDims = 3
DimSize = 512  512  148
ElementSpacing =  0.69531  0.69531  2.00001
Position =  -17.76523  -42.56523  -6.5
ElementByteOrderMSB = False
ElementType = MET_USHORT
ElementDataFile = ImageSet_5.img



my code is 


#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkImageRegionIterator.h"
#include "itkImageConstIterator.h"
#include "itkNeighborhoodConnectedImageFilter.h"
#include "itkMaskNegatedImageFilter.h"
#include "itkMaskImageFilter.h"
#include "itkVotingBinaryIterativeHoleFillingImageFilter.h"
#include "itkBinaryErodeImageFilter.h"
#include "itkBinaryDilateImageFilter.h"
#include "itkBinaryBallStructuringElement.h"
#include "itkRegionOfInterestImageFilter.h"

#define IMAGE_DIMENSIONS 3
#define IMAGE_DATA_EXT short
#define IMAGE_DATA_INT unsigned short
#define INITIAL_THRESHOLD_VALUE 128


typedef itk::Image< IMAGE_DATA_EXT, IMAGE_DIMENSIONS > ImageTypeExt;
typedef itk::Image< IMAGE_DATA_INT, IMAGE_DIMENSIONS > ImageTypeInt;
typedef itk::ImageFileReader< ImageTypeInt > ReaderType;
typedef itk::ImageFileWriter< ImageTypeInt > WriterType;
typedef itk::RescaleIntensityImageFilter< ImageTypeInt, ImageTypeInt >
RescaleType;
typedef itk::BinaryThresholdImageFilter < ImageTypeInt, ImageTypeInt > 
ThresholdFilterType;
typedef itk::ImageRegionConstIterator< ImageTypeInt > ConstIteratorType;
typedef itk::ImageRegionIterator< ImageTypeInt > IteratorType;
typedef itk::NeighborhoodConnectedImageFilter< ImageTypeInt, ImageTypeInt >
ConnectedFilterType;
typedef itk::MaskNegatedImageFilter<ImageTypeInt, ImageTypeInt,
ImageTypeInt> MaskNegatedFilterType;
typedef itk::MaskImageFilter<ImageTypeInt, ImageTypeInt> MaskFilterType;
typedef itk::VotingBinaryIterativeHoleFillingImageFilter< ImageTypeInt >
HoleFillingType;
typedef itk::BinaryBallStructuringElement < IMAGE_DATA_INT, IMAGE_DIMENSIONS
> StructuringElementType;
typedef itk::BinaryErodeImageFilter < ImageTypeInt, ImageTypeInt,
StructuringElementType> ErodeFilterType;
typedef itk::BinaryDilateImageFilter < ImageTypeInt, ImageTypeInt,
StructuringElementType> DilateFilterType;
typedef itk::RegionOfInterestImageFilter< ImageTypeInt, ImageTypeInt >
ROIFilterType;


void invertBinaryValues(ImageTypeInt *image) {
	IteratorType iterator(image, image->GetLargestPossibleRegion());
	for(iterator.GoToBegin(); !iterator.IsAtEnd(); ++iterator) {
		iterator.Set(1 - iterator.Value());
	}
}

int main (int argc, char *argv[])
{
	const char * inputFilename  ="C:/Users/Desktop/lung images/ImageSet_5.mha";
	const char * outputFilename ="C:/Users/Desktop/lung
images/ImageSettest.mha";

	ReaderType::Pointer reader = ReaderType::New();
	reader->SetFileName(inputFilename);
	WriterType::Pointer writer = WriterType::New();
	RescaleType::Pointer rescalerFilter = RescaleType::New();
	rescalerFilter->SetInput(reader->GetOutput());
	rescalerFilter->SetOutputMinimum(0);
	rescalerFilter->SetOutputMaximum(255);
	
	IMAGE_DATA_INT thresholdValue = INITIAL_THRESHOLD_VALUE;
	IMAGE_DATA_INT thresholdValueOLD = INITIAL_THRESHOLD_VALUE + 1;
	ImageTypeInt::IndexType index;
	ThresholdFilterType::Pointer thresholdFilter = ThresholdFilterType::New();
	thresholdFilter->SetOutsideValue(1); 
	thresholdFilter->SetInsideValue(0); 
	thresholdFilter->SetLowerThreshold(0);
	thresholdFilter->SetInput(rescalerFilter->GetOutput());
	
	while(thresholdValue != thresholdValueOLD) {
		unsigned long long averageAbove, averageBelow;
		unsigned int quantityAbove, quantityBelow;
		
		averageAbove = averageBelow = quantityAbove = quantityBelow = 0;
		
		thresholdFilter->SetUpperThreshold(thresholdValue);
		thresholdFilter->Update();

		ConstIteratorType threshold_it(thresholdFilter->GetOutput(),
thresholdFilter->GetOutput()->GetLargestPossibleRegion());
		ConstIteratorType reader_it(rescalerFilter->GetOutput(),
rescalerFilter->GetOutput()->GetLargestPossibleRegion());
		for(threshold_it.GoToBegin(), reader_it.GoToBegin();
!threshold_it.IsAtEnd(); ++threshold_it, ++reader_it) {
			if(threshold_it.Get()) {
				averageAbove += reader_it.Get();
				quantityAbove++;
			}
			else {
				averageBelow += reader_it.Get();
				quantityBelow++;
			}
		}
		if(quantityAbove==0)quantityAbove=1;
		if(quantityBelow==0)quantityBelow=1;

		averageAbove /= quantityAbove;
		averageBelow /= quantityBelow;
		thresholdValueOLD = thresholdValue;
		thresholdValue = (averageAbove + averageBelow)/2;
	}
	invertBinaryValues(thresholdFilter->GetOutput());

	ConnectedFilterType::Pointer neighborhoodConnected =
ConnectedFilterType::New();
 	neighborhoodConnected->SetInput( thresholdFilter->GetOutput() );
	neighborhoodConnected->SetLower(1);
	neighborhoodConnected->SetUpper(1);
	neighborhoodConnected->SetReplaceValue(1);//TODO: magic number
	//Seeds
	index[0] = 0;
	index[1] = 0;
	index[2] = 0;
	neighborhoodConnected->SetSeed(index);
	{
		unsigned int x =
thresholdFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[0];
		unsigned int y =
thresholdFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[1];
		unsigned int z =
thresholdFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[2];
		int i;
		//for(int k =0 ; k<z ; k++){
		for (i = 1; i < x; i++) {
			index[0] = i;
			index[1] = 0;
			index[2] = 0;
			neighborhoodConnected->AddSeed(index);
		}

		//for(int k =0 ; k<z ; k++){
		for (i = 0; i < x; i++) {
			index[0] = i;
			index[1] = y-1;
			index[2] = 0;
			neighborhoodConnected->AddSeed(index);
		}

	}
	ImageTypeInt::SizeType radius;
 	radius[0] = 0;
 	radius[1] = 0;
	radius[2] = 0;
 	neighborhoodConnected->SetRadius(radius);
	
	MaskNegatedFilterType::Pointer maskNegatedFilter =
MaskNegatedFilterType::New();
	maskNegatedFilter->SetInput1(thresholdFilter->GetOutput());
	maskNegatedFilter->SetInput2(neighborhoodConnected->GetOutput());

	HoleFillingType::Pointer holeFillingFilter = HoleFillingType::New();
	radius[0] = 2;
 	radius[1] = 2;
	radius[2] = 2;
 	holeFillingFilter->SetRadius(radius);
	holeFillingFilter->SetInput(maskNegatedFilter->GetOutput());
	holeFillingFilter->SetBackgroundValue( 1 );
	holeFillingFilter->SetForegroundValue( 0 );
	holeFillingFilter->SetMajorityThreshold( 2 );
	holeFillingFilter->SetMaximumNumberOfIterations( 4 );
	
	DilateFilterType::Pointer binaryDilateFilter = DilateFilterType::New();   
 	ErodeFilterType::Pointer binaryErodeFilter = ErodeFilterType::New();
 	
 	StructuringElementType structuringElement;
 	structuringElement.SetRadius( 6 ); 
 	structuringElement.CreateStructuringElement();
 	binaryDilateFilter->SetKernel( structuringElement );
 	binaryErodeFilter->SetKernel( structuringElement );
 	
 	binaryDilateFilter->SetInput( holeFillingFilter->GetOutput() );
 	binaryErodeFilter->SetInput( binaryDilateFilter->GetOutput() );
 	binaryDilateFilter->SetDilateValue( 1 );
 	binaryErodeFilter->SetErodeValue( 1 );


	ImageTypeInt *image;
	binaryErodeFilter->Update();
	image = binaryErodeFilter->GetOutput();
	

	writer->SetFileName( outputFilename );
	writer->SetInput(image);
	writer->Update();





--
View this message in context: http://itk-insight-users.2283740.n2.nabble.com/segmentation-a-3D-image-takes-a-long-time-any-help-to-speed-it-tp7585238.html
Sent from the ITK Insight Users mailing list archive at Nabble.com.