[Insight-users] Help-ERROR: VTKImageImport
firat.sarialtun at boun.edu.tr
firat.sarialtun at boun.edu.tr
Wed Mar 25 08:54:54 EDT 2009
Hi. I am a new user and I have been trying to implement
GeodesicActiveContourLevelSetImageFilter in
ITK for the purpose of Liver Segmentation. I modified the code in order to take
a VTK image as input, convert it into ITK image, process it with the filter,
then reconvert the segmented image into VTK image again. But i have faced with
the problem below. by the way I use the type float for all input variables.
[frat at orion build]$ ./liver
/home/frat/Desktop/liverdata/liverDicom/GE-ser002img?00150?.dcm 1 56
/home/frat/Desktop/outtt.jpg 650 300 5.0 1.0 -0.5 3.0 2.0 250 500
terminate called after throwing an instance of 'itk::ExceptionObject'
what(): /usr/local/ITK/ITK/Code/BasicFilters/itkVTKImageImport.txx:251:
itk::ERROR: VTKImageImport(0xa4206f0): Input scalar type is double but should be
float
Aborted
my Code is :
#include <iostream>
#include <fstream>
#include "vtkImageActor.h"
#include "vtkImageMapToColors.h"
#include "vtkColorTransferFunction.h"
#include "vtkRenderWindow.h"
#include "vtkRenderer.h"
#include "vtkRenderWindowInteractor.h"
#include "itkImage.h"
#include
"/home/frat/Desktop/liverdata/svn/src/VavLoader/itkImageToVTKImageFilter.h"
#include
"/home/frat/Desktop/liverdata/svn/src/VavLoader/itkVTKImageToImageFilter.h"
#include "/home/frat/Desktop/liverdata/VavDataLoader/VavDataLoader.h"
#include "/home/frat/Desktop/liverdata/VavDICOMLoader/VavDICOMLoader.h"
#include "itkGeodesicActiveContourLevelSetImageFilter.h"
#include
"/usr/local/ITK/ITK/Code/BasicFilters/itkCurvatureAnisotropicDiffusionImageFilter.h"
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkSigmoidImageFilter.h"
#include "itkFastMarchingImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#ifdef __BORLANDC__
#define ITK_LEAN_AND_MEAN
#endif
using namespace std;
void render(vtkImageData *data, int viewSlice, int width, int height){
}
int main( int argc, char * argv[] )
{
// argv[1] = first file
// argv[2] = modality
// argv[3] = number of slices
VavDataLoader *loader;
loader = new VavDataLoader( );
int modality = atoi( argv[2] );
int slices = atoi( argv[3] );
VavScalar* data_scalar;
data_scalar = new VavScalar();
VavCT* CT = new VavCT();
loader -> load( argv[1], CT, NULL, modality, slices );
//VTK -> ITK Conversion
typedef itk::Image< float, 2 > InputImageType;
typedef itk::VTKImageToImageFilter< InputImageType >
VTK2ITKConnectorFilterType;
VTK2ITKConnectorFilterType::Pointer VTK2ITKconnector =
VTK2ITKConnectorFilterType::New();
VTK2ITKconnector->SetInput( CT -> dataScalar -> getImageArray() );
//insert your ITK filters here
if( argc < 12 )
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " inputImage outputImage";
std::cerr << " seedX seedY InitialDistance";
std::cerr << " Sigma SigmoidAlpha SigmoidBeta";
std::cerr << " PropagationScaling" << std::endl;
return 1;
}
typedef float InternalPixelType;
const unsigned int Dimension = 2;
typedef itk::Image< InternalPixelType, Dimension > InternalImageType;
// Software Guide : EndCodeSnippet
// The following lines instantiate the thresholding filter that will
// process the final level set at the output of the
// GeodesicActiveContourLevelSetImageFilter.
//
typedef unsigned char OutputPixelType;
typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
typedef itk::BinaryThresholdImageFilter<
InternalImageType,
OutputImageType > ThresholdingFilterType;
ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();
thresholder->SetLowerThreshold( -1000.0 );
thresholder->SetUpperThreshold( 0.0 );
thresholder->SetOutsideValue( 0 );
thresholder->SetInsideValue( 255 );
// We instantiate reader and writer types in the following lines.
//
typedef itk::ImageFileReader< InternalImageType > ReaderType;
typedef itk::ImageFileWriter< OutputImageType > WriterType;
ReaderType::Pointer reader = ReaderType::New();
WriterType::Pointer writer = WriterType::New();
reader->SetFileName( argv[1] );
writer->SetFileName( argv[4] );
// The RescaleIntensityImageFilter type is declared below. This filter will
// renormalize image before sending them to writers.
//
typedef itk::RescaleIntensityImageFilter<
InternalImageType,
OutputImageType > CastFilterType;
// The \doxygen{CurvatureAnisotropicDiffusionImageFilter} type is
// instantiated using the internal image type.
//
typedef itk::CurvatureAnisotropicDiffusionImageFilter<
InternalImageType,
InternalImageType > SmoothingFilterType;
SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();
// The types of the
// GradientMagnitudeRecursiveGaussianImageFilter and
// SigmoidImageFilter are instantiated using the internal image
// type.
//
typedef itk::GradientMagnitudeRecursiveGaussianImageFilter<
InternalImageType,
InternalImageType > GradientFilterType;
typedef itk::SigmoidImageFilter<
InternalImageType,
InternalImageType > SigmoidFilterType;
GradientFilterType::Pointer gradientMagnitude = GradientFilterType::New();
SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();
// The minimum and maximum values of the SigmoidImageFilter output
// are defined with the methods \code{SetOutputMinimum()} and
// \code{SetOutputMaximum()}. In our case, we want these two values to be
// $0.0$ and $1.0$ respectively in order to get a nice speed image to feed
// the \code{FastMarchingImageFilter}. Additional details on the user of the
// \doxygen{SigmoidImageFilter} are presented in
// section~\ref{sec:IntensityNonLinearMapping}.
sigmoid->SetOutputMinimum( 0.0 );
sigmoid->SetOutputMaximum( 1.0 );
// We declare now the type of the FastMarchingImageFilter that
// will be used to generate the initial level set in the form of a distance
// map.
//
typedef itk::FastMarchingImageFilter<
InternalImageType,
InternalImageType > FastMarchingFilterType;
// Next we construct one filter of this class using the \code{New()}
// method.
//
FastMarchingFilterType::Pointer fastMarching =
FastMarchingFilterType::New();
// Software Guide : BeginLatex
//
// In the following lines we instantiate the type of the
// GeodesicActiveContourLevelSetImageFilter and create an object of this
// type using the \code{New()} method.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::GeodesicActiveContourLevelSetImageFilter< InternalImageType,
InternalImageType > GeodesicActiveContourFilterType;
GeodesicActiveContourFilterType::Pointer geodesicActiveContour =
GeodesicActiveContourFilterType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// For the GeodesicActiveContourLevelSetImageFilter, scaling parameters
// are used to trade off between the propagation (inflation), the
// curvature (smoothing) and the advection terms. These parameters are set
// using methods \code{SetPropagationScaling()},
// \code{SetCurvatureScaling()} and \code{SetAdvectionScaling()}. In this
// example, we will set the curvature and advection scales to one and let
// the propagation scale be a command-line argument.
//
//
\index{itk::Geodesic\-Active\-Contour\-LevelSet\-Image\-Filter!SetPropagationScaling()}
//
\index{itk::Segmentation\-Level\-Set\-Image\-Filter!SetPropagationScaling()}
//
\index{itk::Geodesic\-Active\-Contour\-LevelSet\-Image\-Filter!SetCurvatureScaling()}
//
\index{itk::Segmentation\-Level\-Set\-Image\-Filter!SetCurvatureScaling()}
//
\index{itk::Geodesic\-Active\-Contour\-LevelSet\-Image\-Filter!SetAdvectionScaling()}
//
\index{itk::Segmentation\-Level\-Set\-Image\-Filter!SetAdvectionScaling()}
//
// Software Guide : EndLatex
const double propagationScaling = atof( argv[11] );
// Software Guide : BeginCodeSnippet
geodesicActiveContour->SetPropagationScaling( propagationScaling );
geodesicActiveContour->SetCurvatureScaling( 1.0 );
geodesicActiveContour->SetAdvectionScaling( 1.0 );
// Software Guide : EndCodeSnippet
// Once activiated the level set evolution will stop if the convergence
// criteria or if the maximum number of iterations is reached. The
// convergence criteria is defined in terms of the root mean squared (RMS)
// change in the level set function. The evolution is said to have
// converged if the RMS change is below a user specified threshold. In a
// real application is desirable to couple the evolution of the zero set
// to a visualization module allowing the user to follow the evolution of
// the zero set. With this feedback, the user may decide when to stop the
// algorithm before the zero set leaks through the regions of low gradient
// in the contour of the anatomical structure to be segmented.
geodesicActiveContour->SetMaximumRMSError( 0.02 );
geodesicActiveContour->SetNumberOfIterations( 800 );
// Software Guide : BeginLatex
//
// The filters are now connected in a pipeline indicated in
// Figure~\ref{fig:GeodesicActiveContoursCollaborationDiagram} using the
// following lines:
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
// smoothing->SetInput( reader->GetOutput()
);///////////////////////////////////////////// inputtttt
smoothing->SetInput( VTK2ITKconnector->GetOutput() );
gradientMagnitude->SetInput( smoothing->GetOutput() );
sigmoid->SetInput( gradientMagnitude->GetOutput() );
geodesicActiveContour->SetInput( fastMarching->GetOutput() );
geodesicActiveContour->SetFeatureImage( sigmoid->GetOutput() );
thresholder->SetInput( geodesicActiveContour->GetOutput() );
writer->SetInput( thresholder->GetOutput() );
// Software Guide : EndCodeSnippet
// The CurvatureAnisotropicDiffusionImageFilter requires a couple of
// parameter to be defined. The following are typical values for $2D$
// images. However they may have to be adjusted depending on the amount of
// noise present in the input image. This filter has been discussed in
// section~\ref{sec:GradientAnisotropicDiffusionImageFilter}.
smoothing->SetTimeStep( 0.125 );
smoothing->SetNumberOfIterations( 5 );
smoothing->SetConductanceParameter( 9.0 );
// The GradientMagnitudeRecursiveGaussianImageFilter performs the
// equivalent of a convolution with a Gaussian kernel, followed by a
// derivative operator. The sigma of this Gaussian can be used to control
// the range of influence of the image edges. This filter has been discussed
// in Section~\ref{sec:GradientMagnitudeRecursiveGaussianImageFilter}
const double sigma = atof( argv[8] );
gradientMagnitude->SetSigma( sigma );
// The SigmoidImageFilter requires two parameters that define the linear
// transformation to be applied to the sigmoid argument. This parameters
// have been discussed in Sections~\ref{sec:IntensityNonLinearMapping} and
// \ref{sec:FastMarchingImageFilter}.
const double alpha = atof( argv[9] );
const double beta = atof( argv[10] );
sigmoid->SetAlpha( alpha );
sigmoid->SetBeta( beta );
// The FastMarchingImageFilter requires the user to provide a seed
// point from which the level set will be generated. The user can actually
// pass not only one seed point but a set of them. Note the the
// FastMarchingImageFilter is used here only as a helper in the
// determination of an initial level set. We could have used the
// \doxygen{DanielssonDistanceMapImageFilter} in the same way.
//
// The seeds are passed stored in a container. The type of this
// container is defined as \code{NodeContainer} among the
// FastMarchingImageFilter traits.
//
typedef FastMarchingFilterType::NodeContainer NodeContainer;
typedef FastMarchingFilterType::NodeType NodeType;
NodeContainer::Pointer seeds = NodeContainer::New();
InternalImageType::IndexType seedPosition;
seedPosition[0] = atoi( argv[5] );
seedPosition[1] = atoi( argv[6] );
// Nodes are created as stack variables and initialized with a value and an
// \doxygen{Index} position. Note that here we assign the value of minus the
// user-provided distance to the unique node of the seeds passed to the
// FastMarchingImageFilter. In this way, the value will increment
// as the front is propagated, until it reaches the zero value corresponding
// to the contour. After this, the front will continue propagating until it
// fills up the entire image. The initial distance is taken here from the
// command line arguments. The rule of thumb for the user is to select this
// value as the distance from the seed points at which she want the initial
// contour to be.
const double initialDistance = atof( argv[7] );
NodeType node;
const double seedValue = - initialDistance;
node.SetValue( seedValue );
node.SetIndex( seedPosition );
// The list of nodes is initialized and then every node is inserted using
// the \code{InsertElement()}.
seeds->Initialize();
seeds->InsertElement( 0, node );
// modification by firat~~~~~~~~~~~~~~~~~~~~~~~More Seed Points
/*
NodeType node2;
seedPosition[0] = atoi( argv[10] );;
seedPosition[1] = atoi( argv[11] );;
node2.SetValue( seedValue );
node2.SetIndex( seedPosition );
seeds->InsertElement( 1, node2 );*/
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// The set of seed nodes is passed now to the
// FastMarchingImageFilter with the method
// \code{SetTrialPoints()}.
//
fastMarching->SetTrialPoints( seeds );
// Since the FastMarchingImageFilter is used here just as a
// Distance Map generator. It does not require a speed image as input.
// Instead the constant value $1.0$ is passed using the
// \code{SetSpeedConstant()} method.
//
fastMarching->SetSpeedConstant( 1.0 );
// Here we configure all the writers required to see the intermediate
// outputs of the pipeline. This is added here only for
// pedagogical/debugging purposes. These intermediate output are normaly not
// required. Only the output of the final thresholding filter should be
// relevant. Observing intermediate output is helpful in the process of
// fine tuning the parameters of filters in the pipeline.
//
CastFilterType::Pointer caster1 = CastFilterType::New();
CastFilterType::Pointer caster2 = CastFilterType::New();
CastFilterType::Pointer caster3 = CastFilterType::New();
CastFilterType::Pointer caster4 = CastFilterType::New();
WriterType::Pointer writer1 = WriterType::New();
WriterType::Pointer writer2 = WriterType::New();
WriterType::Pointer writer3 = WriterType::New();
WriterType::Pointer writer4 = WriterType::New();
caster1->SetInput( smoothing->GetOutput() );
writer1->SetInput( caster1->GetOutput() );
writer1->SetFileName("GeodesicActiveContourImageFilterOutput1.png");
caster1->SetOutputMinimum( 0 );
caster1->SetOutputMaximum( 255 );
writer1->Update();
caster2->SetInput( gradientMagnitude->GetOutput() );
writer2->SetInput( caster2->GetOutput() );
writer2->SetFileName("GeodesicActiveContourImageFilterOutput2.png");
caster2->SetOutputMinimum( 0 );
caster2->SetOutputMaximum( 255 );
writer2->Update();
caster3->SetInput( sigmoid->GetOutput() );
writer3->SetInput( caster3->GetOutput() );
writer3->SetFileName("GeodesicActiveContourImageFilterOutput3.png");
caster3->SetOutputMinimum( 0 );
caster3->SetOutputMaximum( 255 );
writer3->Update();
caster4->SetInput( fastMarching->GetOutput() );
writer4->SetInput( caster4->GetOutput() );
writer4->SetFileName("GeodesicActiveContourImageFilterOutput4.png");
caster4->SetOutputMinimum( 0 );
caster4->SetOutputMaximum( 255 );
// The FastMarchingImageFilter requires the user to specify the
// size of the image to be produced as output. This is done using the
// \code{SetOutputSize()}. Note that the size is obtained here from the
// output image of the smoothing filter. The size of this image is valid
// only after the \code{Update()} methods of this filter has been called
// directly or indirectly.
//
fastMarching->SetOutputSize(
reader->GetOutput()->GetBufferedRegion().GetSize() );
// Software Guide : BeginLatex
//
// The invocation of the \code{Update()} method on the writer triggers the
// execution of the pipeline. As usual, the call is placed in a
// \code{try/catch} block should any errors occur or exceptions be thrown.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
try
{
writer->Update();
}
catch( itk::ExceptionObject & excep )
{
std::cerr << "Exception caught !" << std::endl;
std::cerr << excep << std::endl;
}
// Software Guide : EndCodeSnippet
// Print out some useful information
std::cout << std::endl;
std::cout << "Max. no. iterations: " <<
geodesicActiveContour->GetNumberOfIterations() << std::endl;
std::cout << "Max. RMS error: " << geodesicActiveContour->GetMaximumRMSError()
<< std::endl;
std::cout << std::endl;
std::cout << "No. elpased iterations: " <<
geodesicActiveContour->GetElapsedIterations() << std::endl;
std::cout << "RMS change: " << geodesicActiveContour->GetRMSChange() <<
std::endl;
writer4->Update();
// The following writer type is used to save the output of the time-crossing
// map in a file with apropiate pixel representation. The advantage of saving
// this image in native format is that it can be used with a viewer to help
// determine an appropriate threshold to be used on the output of the
// fastmarching filter.
//
typedef itk::ImageFileWriter< InternalImageType > InternalWriterType;
InternalWriterType::Pointer mapWriter = InternalWriterType::New();
mapWriter->SetInput( fastMarching->GetOutput() );
mapWriter->SetFileName("GeodesicActiveContourImageFilterOutput4.mha");
mapWriter->Update();
InternalWriterType::Pointer speedWriter = InternalWriterType::New();
speedWriter->SetInput( sigmoid->GetOutput() );
speedWriter->SetFileName("GeodesicActiveContourImageFilterOutput3.mha");
speedWriter->Update();
InternalWriterType::Pointer gradientWriter = InternalWriterType::New();
gradientWriter->SetInput( gradientMagnitude->GetOutput() );
gradientWriter->SetFileName("GeodesicActiveContourImageFilterOutput2.mha");
gradientWriter->Update();
//ITK -> VTK Conversion
typedef itk::ImageToVTKImageFilter< InputImageType >
ITK2VTKConnectorFilterType;
ITK2VTKConnectorFilterType::Pointer ITK2VTKconnector =
ITK2VTKConnectorFilterType::New();
ITK2VTKconnector->GetExporter()->SetInput(VTK2ITKconnector->GetOutput());
ITK2VTKconnector->GetImporter()->Update();
//display resulting image
//render(ITK2VTKconnector->GetImporter()->GetOutput(), 1, 512, 512);
cout << "Done!" << endl;
return 0;
}
More information about the Insight-users
mailing list