[ITK] [ITK-users] segmentation a 3D image takes a long time , any help to speed it ?
alaamegawer
alaamegawer at yahoo.com
Tue Mar 25 07:37:29 EDT 2014
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.
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php
Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ
Follow this link to subscribe/unsubscribe:
http://www.itk.org/mailman/listinfo/insight-users
More information about the Community
mailing list