[Insight-users] applying Shape Detection filter successively, causes exception

neha k itkneha1 at yahoo.com
Wed Jan 12 12:15:35 EST 2005


Skipped content of type multipart/alternative-------------- next part --------------
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkGradientAnisotropicDiffusionImageFilter.h"
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkConnectedThresholdImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkCastImageFilter.h"
// Filters for Level Set Filters
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkSigmoidImageFilter.h"
#include "itkFastMarchingImageFilter.h"
#include "itkShapeDetectionLevelSetImageFilter.h"
#include "itkBinaryThresholdImageFilter.h"
// Morphology Filters
#include "itkBinaryErodeImageFilter.h"
#include "itkBinaryDilateImageFilter.h"
// Header File Structuring Element for Morphology Filters
#include "itkBinaryBallStructuringElement.h"

int ShapeDetectionLevelSet();

int main()
{
	ShapeDetectionLevelSet();
}

int ShapeDetectionLevelSet()
{
	typedef float InputPixelType;
	typedef unsigned char OutputPixelType;
	typedef unsigned char WritePixelType;

	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
	typedef itk::Image  < OutputPixelType, 2 >   OutputImageType; //unsigned char
	
	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 using structuring element
	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;
	//typedef itk::BinaryErodeImageFilter<BinaryImageType,BinaryImageType,StructuringElementType> ErodeFilterType;

	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;
	}  
	
	connectedThreshold->SetInput( reader->GetOutput() );		
	connectedThreshold->SetLower( 169 ); 
	connectedThreshold->SetUpper( 185 ); 
	connectedThreshold->SetReplaceValue( 255 );			
	seedIndex[0] = 76; // first index on X 
	seedIndex[1] = 93; // first index on Y	
	connectedThreshold->SetSeed( seedIndex );
	
	// Set parameters for Gradient Magnitude Filter 
	
	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< InputImageType, 
                              InputImageType>    ShapeDetectionFilterType1;
	ShapeDetectionFilterType1::Pointer shapeDetection1 = ShapeDetectionFilterType1::New();

	typedef   itk::SigmoidImageFilter< InputImageType, InputImageType> SigmoidFilterType1;	
	SigmoidFilterType1::Pointer sigmoid1 = SigmoidFilterType1::New();

	typedef itk::BinaryThresholdImageFilter< InputImageType, BinaryImageType> ThresholdingFilterType1;  
	ThresholdingFilterType1::Pointer thresholder1 = ThresholdingFilterType1::New(); 
	
	typedef itk::CastImageFilter<BinaryImageType, InputImageType> CastToShortFilterType;
	CastToShortFilterType::Pointer uchar_to_float = CastToShortFilterType::New();			

	uchar_to_float->SetInput( binaryDilate->GetOutput() );	
	shapeDetection1->SetInput( uchar_to_float->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 << 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