[Insight-users] Re: Brain Tumor + Level Set Filter
Luis Ibanez
luis.ibanez at kitware.com
Thu Jan 6 22:33:54 EST 2005
Hi Neha,
Please read the Doxygen documentation of the filters that you use.
http://www.itk.org/Insight/Doxygen/html/classitk_1_1ShapeDetectionLevelSetImageFilter.html
Level set filters are typically templated over three types:
1) Image type of the input image (initial level set)
2) Image type of the feature image (used for the speed image)
3) Pixel type of the output image (final level set).
Note that the first two are ImageTypes and the last one is a PixelType.
The Third template parameter is set by default to "float" because it
is the type over which the LevelSet computations are performed. Note
that when you run LevelSets you are iteratively solving a partial
differential equations by iterations. In general you don't want to
do this using integer pixels types since they don't provide enough
precision.
To solve the problem in your code, simply instantiate the Binary
Threshold Image Filter using as input image type, the "OutputImageType"
from the ShapeDetectionLevelSetImageFilter.
Something like:
typedef
itk::BinaryThresholdImageFilter<
ShapeDetectionFilterType1::OutputImageType,
BinaryImageType
> ThresholdingFilterType1;
Enchaining the OutputImageType between consecutive filters is a safe
way of making sure that the image types will match between Output and
Input.
Regards,
Luis
=====================
neha k wrote:
> Hello Luis,
>
> Thanks for the reply. Currently, I am trying to use Shape Detection
> Level Set along with Conencted Threshold to delineate Tumour. With
> Shape Detection, it gives me tumor with inner boundary. To separate
> Tumor completely with its outer bounday, I am trying to give binary
> output to Shape Detection again. I have attached the program. In that,
> I receive error in the following lines of code -->
>
> typedef unsigned char BinaryPixelType; // For Dilate and Erode
> oparations using structuring element
> typedef itk::Image< BinaryPixelType, Dimension > BinaryImageType;
> typedef itk::ShapeDetectionLevelSetImageFilter< BinaryImageType,
> BinaryImageType> ShapeDetectionFilterType1;
> ShapeDetectionFilterType1::Pointer shapeDetection1 =
> ShapeDetectionFilterType1::New();
> typedef itk::SigmoidImageFilter< InputImageType, BinaryImageType >
> SigmoidFilterType1;
> SigmoidFilterType1::Pointer sigmoid1 = SigmoidFilterType1::New();
> *typedef itk::BinaryThresholdImageFilter< BinaryImageType,
> BinaryImageType> ThresholdingFilterType1;
> T*hresholdingFilterType1::Pointer thresholder1 =
> ThresholdingFilterType1::New();
>
> shapeDetection1->SetInput( binaryDilate->GetOutput() );
> shapeDetection1->SetFeatureImage( sigmoid1->GetOutput() );
> *thresholder1->SetInput( shapeDetection1->GetOutput() ); // ERROR
> *writer->SetInput( thresholder1->GetOutput() );
> Error is --> can not convert from <float,2> to <unsigned char,2>.
> Here, for shapeDetection1, all input is Binary. Same case for Thresholder.
> If original image is read as "float" and given to sigmoid filter, and
> feature image is binary (unsigned char), is above declaration of
> ShapeDetection correct? Is it ok for take output of Sigmoid as Binary?
>
> Thanks,
> Neha
>
>
>
>
>
>
> ------------------------------------------------------------------------
>
> int ShapeDetectionLevelSet()
> {
> int noOfIterations = 10, checkExtn = 0;
> char inputFile[255]="c:\\tumour57_smooth.png", // unsigned short type outputFile[255]="c:\\tumour57_smooth_ShapeDetection.png"; float propagationScaling = 0.9, curvatureScaling = 2.74;
> const double sigma = 1.2; const double alpha = -2.9; const double beta = 4.1;
>
> const unsigned int Dimension = 2;
> typedef itk::Image < InputPixelType, 2 > InputImageType; //float - as in guide
> typedef itk::Image < OutputPixelType, 2 > OutputImageType; //unsigned char - as in guide
>
> typedef itk::ImageFileReader< InputImageType > ReaderType; // read input image as float
> ReaderType::Pointer reader = ReaderType::New();
>
> typedef itk::Image<unsigned char, 2> WriteImageType;
> typedef itk::ImageFileWriter< WriteImageType > WriterType; // write output image as short - use rescaler
> WriterType::Pointer writer = WriterType::New();
>
> // declarations for "Float" type connectedThrshold Filter
> typedef itk::ConnectedThresholdImageFilter <InputImageType, InputImageType> ConnectedFilterType;
> ConnectedFilterType::Pointer connectedThreshold = ConnectedFilterType::New();
> // This has output image of float type, since this o/p goes to Geodesic Active Filter
>
> // Rescaler - To rescale From Float to Unsigned Short
> typedef itk::RescaleIntensityImageFilter< OutputImageType, WriteImageType > RescaleFilterType;
> RescaleFilterType::Pointer rescaler = RescaleFilterType::New();
> typedef itk::RescaleIntensityImageFilter< InputImageType, WriteImageType > RescaleFilterType1;
>
> typedef itk::GradientMagnitudeRecursiveGaussianImageFilter<InputImageType, InputImageType > FilterType1;
> typedef itk::GradientMagnitudeRecursiveGaussianImageFilter< InputImageType, InputImageType> GradientFilterType;
>
> typedef itk::SigmoidImageFilter< InputImageType, InputImageType > SigmoidFilterType;
>
> GradientFilterType::Pointer gradientMagnitude = GradientFilterType::New();
> FilterType1::Pointer recursiveGaussian = FilterType1::New();
>
> SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();
>
> typedef itk::ShapeDetectionLevelSetImageFilter< InputImageType,
> InputImageType > ShapeDetectionFilterType;
> ShapeDetectionFilterType::Pointer shapeDetection = ShapeDetectionFilterType::New();
>
> typedef itk::BinaryThresholdImageFilter< InputImageType, OutputImageType> ThresholdingFilterType;
> ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();
>
> typedef unsigned char BinaryPixelType; // For Dilate and Erode oparations
> typedef itk::Image< BinaryPixelType, Dimension > BinaryImageType;
>
> typedef itk::BinaryBallStructuringElement<BinaryPixelType,Dimension > StructuringElementType;
> // The structuring element type is then used along with the input and output image types for
> // instantiating the type of the filters.
> typedef itk::BinaryDilateImageFilter<BinaryImageType,BinaryImageType,StructuringElementType> DilateFilterType;
>
> BinaryPixelType background = 0; // Values on which dilate/erosion filters operate
> BinaryPixelType foreground = 255;
> DilateFilterType::Pointer binaryDilate = DilateFilterType::New();
>
> // Define and set index - for connected Threshold Seg. filter
> typedef OutputImageType::IndexType IndexType;
> OutputImageType::IndexType seedIndex;
>
> reader->SetFileName(inputFile);
> std::cout <<"\n\nReading A Raw File...\n";
> try {
> reader->Update();
> }
> catch( itk::ExceptionObject & e ) {
> std::cerr << "Exception caught during Raw file reading " << std::endl;
> std::cerr << e << std::endl;
> std::cout<< "\n\nMake Sure if the file name you specified already exists...\n";
> return -1;
> }
> std::cout << "File read succesfully ! " << std::endl;
>
> // NOTE - If you need to use smoothing, input image to this progarm should already be smoothed.
>
> connectedThreshold->SetInput( reader->GetOutput() );
> connectedThreshold->SetLower( 169 ); // 27
> connectedThreshold->SetUpper( 185 ); // 186
> connectedThreshold->SetReplaceValue( 255 );
> seedIndex[0] = 76; // first index on X // 84;//75;
> seedIndex[1] = 93; // first index on Y //116;//93;
> connectedThreshold->SetSeed( seedIndex );
>
> recursiveGaussian->SetSigma( sigma );
>
> // Set parameters for Sigmoid Filter
> sigmoid->SetOutputMinimum( 0.0 ); // What is this?? do I need to modify it??
> sigmoid->SetOutputMaximum( 1.0 ); // What is this?? do I need to modify it??
> sigmoid->SetAlpha( alpha );
> sigmoid->SetBeta( beta );
>
> // Set Paramters for Shape Detection Level Set
> shapeDetection->SetPropagationScaling( propagationScaling );
> shapeDetection->SetCurvatureScaling( curvatureScaling );
> shapeDetection->SetMaximumRMSError( 0.02 );
> shapeDetection->SetNumberOfIterations( 75 ); // 800
>
> recursiveGaussian->SetInput( reader->GetOutput() );
> sigmoid->SetInput( recursiveGaussian->GetOutput() );
>
> // Set Thresholder output
> thresholder->SetLowerThreshold( -5000.0 );
> thresholder->SetUpperThreshold( 0.0 );
> thresholder->SetOutsideValue( 255 );
> thresholder->SetInsideValue( 0 );
>
> // The filters are now connected in a pipeline
> shapeDetection->SetInput( connectedThreshold->GetOutput() );
> shapeDetection->SetFeatureImage( sigmoid->GetOutput() );
> thresholder->SetInput( shapeDetection->GetOutput() );
> writer->SetInput( thresholder->GetOutput() );
>
> // NOTE - Here we configure all the writers required to see the intermediate
> // outputs of the pipeline. This is added here for debugging purposes. Output of the final thresholding filter
> // is relevant. Observing intermediate output is helpful in the process of fine tuning the parameters
> // of filters in the pipeline.
>
> RescaleFilterType1::Pointer caster1 = RescaleFilterType1::New();
> RescaleFilterType1::Pointer caster2 = RescaleFilterType1::New();
> RescaleFilterType1::Pointer caster3 = RescaleFilterType1::New();
> RescaleFilterType1::Pointer caster4 = RescaleFilterType1::New();
>
> WriterType::Pointer writer1 = WriterType::New();
> WriterType::Pointer writer2 = WriterType::New();
> WriterType::Pointer writer3 = WriterType::New();
> WriterType::Pointer writer4 = WriterType::New();
>
> caster2->SetInput( recursiveGaussian->GetOutput() );
> writer2->SetInput( caster2->GetOutput() );
> writer2->SetFileName("GradientMagnitudeFilterOutput.png");
> caster2->SetOutputMinimum( 0 );
> caster2->SetOutputMaximum( 255 );
> writer2->Update();
>
> caster3->SetInput( sigmoid->GetOutput() );
> writer3->SetInput( caster3->GetOutput() );
> writer3->SetFileName("SigmoidImageFilterOutput.png");
> caster3->SetOutputMinimum( 0 );
> caster3->SetOutputMaximum( 255 );
> writer3->Update();
>
> caster4->SetInput( connectedThreshold->GetOutput() );
> writer4->SetInput( caster4->GetOutput() );
> writer4->SetFileName("ConnectedThresholdOutput.png");
> caster4->SetOutputMinimum( 0 );
> caster4->SetOutputMaximum( 255 );
> writer4->Update();
>
> rescaler->SetOutputMinimum( 0 );
> rescaler->SetOutputMaximum( 255 );
> rescaler->SetInput(thresholder->GetOutput());
>
> // Dilate the result and then apply another Shape Detection Step. Then Write output image
> StructuringElementType structuringElement;
> structuringElement.SetRadius( 1 ); // 3x3 structuring element
> structuringElement.CreateStructuringElement();
> binaryDilate->SetKernel( structuringElement );
>
>
> binaryDilate->SetDilateValue( foreground );
> binaryDilate->SetInput( rescaler->GetOutput() );
> writer->SetFileName( outputFile );
>
> typedef unsigned char BinaryPixelType; // For Dilate and Erode oparations using structuring element
> typedef itk::Image< BinaryPixelType, Dimension > BinaryImageType;
>
> typedef itk::ShapeDetectionLevelSetImageFilter< BinaryImageType,
> BinaryImageType> ShapeDetectionFilterType1;
> ShapeDetectionFilterType1::Pointer shapeDetection1 = ShapeDetectionFilterType1::New();
> typedef itk::SigmoidImageFilter< InputImageType, BinaryImageType > SigmoidFilterType1;
> SigmoidFilterType1::Pointer sigmoid1 = SigmoidFilterType1::New();
>
> typedef itk::BinaryThresholdImageFilter< BinaryImageType, BinaryImageType> ThresholdingFilterType1;
> ThresholdingFilterType1::Pointer thresholder1 = ThresholdingFilterType1::New();
>
> shapeDetection1->SetInput( binaryDilate->GetOutput() );
> shapeDetection1->SetFeatureImage( sigmoid1->GetOutput() );
> thresholder1->SetInput( shapeDetection1->GetOutput() );
> writer->SetInput( thresholder1->GetOutput() );
>
> // Execute the pipeline with Update() call.
> try {
> writer->Update();
> }
> catch( itk::ExceptionObject & excep ){
> std::cerr << "Exception caught !" << std::endl;
> std::cerr << excep << std::endl;
> }
> std::cout << "File succesfully writen ! " << std::endl;
>
> std::cout << std::endl;
> std::cout << "Max. no. iterations: " << shapeDetection->GetNumberOfIterations() << std::endl;
> std::cout << "Max. RMS error: " << shapeDetection->GetMaximumRMSError() << std::endl;
> std::cout << std::endl;
> std::cout << "No. elapsed iterations: " << shapeDetection->GetElapsedIterations() << std::endl;
> std::cout << "RMS change: " << shapeDetection->GetRMSChange() << std::endl;
>
> return 0;
> }
>
More information about the Insight-users
mailing list