[ITK-users] segmentation a 3D image takes a long time , any help to speed it ?
Jim Miller
millerjv at gmail.com
Tue Mar 25 07:39:31 EDT 2014
Are you building Release of Debug?
Jim
> On Mar 25, 2014, at 7:37 AM, alaamegawer <alaamegawer at yahoo.com> wrote:
>
>
> 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 Insight-users
mailing list