[Insight-users] how to write *.vtk file
SUTRISNO SUTRISNO
sutrisno_link at yahoo.com
Tue Apr 13 06:09:58 EDT 2010
Hi Luis,
I tried to segmented DICOM files and save them in vtk format.
I used DicomSeriesReadImageWrite2.cxx, DeformableModel1.cxx, DeformableModel2.cxx to build my program.
the problem is the result does not match what I expected.
when I opened bone.vtk (result file) with file editor the content just like this:
# vtk DataFile Version 3.0
VTK File Generated by Insight Segmentation and Registration Toolkit (ITK)
ASCII
DATASET STRUCTURED_POINTS
DIMENSIONS 124 124 45
SPACING 1.0000000000000000e+00 1.0000000000000000e+00 1.0000000000000000e+00
ORIGIN 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
POINT_DATA 691920
SCALARS scalars short 1
LOOKUP_TABLE default
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
....
and when I opened this file with vtk reader (PolyDataReader) there was error :
ERROR: In /build/buildd/vtk-5.0.3/IO/vtkPolyDataReader.cxx, line 142
vtkPolyDataReader (0x8765510): Cannot read dataset type: structured_points
This is my source code list :
/*=========================================================================
Name : 3D Segmentation of Bone Structure on CT Image
=========================================================================*/
// Include the required header
#include "itkOrientedImage.h"
#include "itkGDCMImageIO.h"
#include "itkGDCMSeriesFileNames.h"
#include "itkImageSeriesReader.h"
#include "itkImageFileWriter.h"
#include "itkMesh.h"
#include "itkDeformableMesh3DFilter.h"
#include "itkGradientRecursiveGaussianImageFilter.h"
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkImage.h"
#include "itkCovariantVector.h"
#include "itkSphereMeshSource.h"
#include "itkImageFileReader.h"
#include "itkPointSetToImageFilter.h"
#include "itkVTKImageIO.h"
int main(int argc, char* argv[] )
{
if( argc < 3 )
{
std::cerr << "Usage: " << std::endl;
std::cerr << argv[0] << " DicomDirectory outputFileName [seriesName]"
<< std::endl;
return EXIT_FAILURE;
}
typedef signed short PixelType;
const unsigned int Dimension = 3;
typedef itk::OrientedImage< PixelType, Dimension > ImageType;
typedef itk::ImageSeriesReader< ImageType > ReaderType;
ReaderType::Pointer reader = ReaderType::New();
typedef itk::GDCMImageIO ImageIOType;
ImageIOType::Pointer dicomIO = ImageIOType::New();
reader->SetImageIO( dicomIO );
typedef itk::GDCMSeriesFileNames NamesGeneratorType;
NamesGeneratorType::Pointer nameGenerator = NamesGeneratorType::New();
nameGenerator->SetUseSeriesDetails( true );
nameGenerator->AddSeriesRestriction("0008|0021" );
nameGenerator->SetDirectory( argv[1] );
try
{
std::cout << std::endl << "The directory: " << std::endl;
std::cout << std::endl << argv[1] << std::endl << std::endl;
std::cout << "Contains the following DICOM Series: ";
std::cout << std::endl << std::endl;
typedef std::vector< std::string > SeriesIdContainer;
const SeriesIdContainer & seriesUID = nameGenerator->GetSeriesUIDs();
SeriesIdContainer::const_iterator seriesItr = seriesUID.begin();
SeriesIdContainer::const_iterator seriesEnd = seriesUID.end();
while( seriesItr != seriesEnd )
{
std::cout << seriesItr->c_str() << std::endl;
seriesItr++;
}
std::string seriesIdentifier;
if( argc > 3 )
{
seriesIdentifier = argv[3];
}
else
{
seriesIdentifier = seriesUID.begin()->c_str();
}
std::cout << std::endl << std::endl;
std::cout << "Now reading series: " << std::endl << std::endl;
std::cout << seriesIdentifier << std::endl;
std::cout << std::endl << std::endl;
//read DICOM
typedef std::vector< std::string > FileNamesContainer;
FileNamesContainer fileNames;
fileNames = nameGenerator->GetFileNames( seriesIdentifier );
reader->SetFileNames( fileNames );
try
{
reader->Update();
}
catch (itk::ExceptionObject &ex)
{
std::cout << ex << std::endl;
return EXIT_FAILURE;
}
//segmentation with DeformableModel
const float sigma = 1.0 ;
const int numberOfIterations = 100 ;
const float timeStep = 0.1 ;
const float externalForceScale = 10 ;
const float stiffness = 0.1 ;
typedef double MeshPixelType; typedef itk::Mesh< MeshPixelType > MeshType;
typedef itk::GradientRecursiveGaussianImageFilter<ImageType> GradientFilterType; typedef itk::GradientMagnitudeRecursiveGaussianImageFilter<ImageType,ImageType> GradientMagnitudeFilterType;
typedef itk::DeformableMesh3DFilter<MeshType,MeshType> DeformableFilterType; GradientMagnitudeFilterType::Pointer gradientMagnitudeFilter = GradientMagnitudeFilterType::New();
ImageType::ConstPointer inputImage = reader->GetOutput(); gradientMagnitudeFilter->SetInput( inputImage ); gradientMagnitudeFilter->SetSigma( sigma ); GradientFilterType::Pointer gradientMapFilter = GradientFilterType::New();
gradientMapFilter->SetInput( gradientMagnitudeFilter->GetOutput()); gradientMapFilter->SetSigma( sigma ); gradientMapFilter->Update();
DeformableFilterType::Pointer deformableModelFilter = DeformableFilterType::New();
typedef itk::SphereMeshSource< MeshType > MeshSourceType; MeshSourceType::Pointer meshSource = MeshSourceType::New();
// Set the initial sphere in the center of the image
const ImageType::SpacingType spacing = inputImage->GetSpacing();
ImageType::PointType origin = inputImage->GetOrigin(); ImageType::SizeType size = inputImage->GetBufferedRegion().GetSize();
MeshType::PointType center;
center[0] = origin[0] + spacing[0] * size[0] / 2.0; center[1] = origin[1] + spacing[1] * size[1] / 2.0; center[2] = origin[2] + spacing[2] * size[2] / 2.0;
meshSource->SetCenter( center );
MeshType::PointType radius; radius[0] = spacing[0] * size[0] / 4.0;
radius[1] = spacing[1] * size[1] / 4.0; radius[2] = spacing[2] * size[2] / 4.0;
meshSource->SetScale( radius );
meshSource->SetResolutionX( 50 ); meshSource->SetResolutionY( 50 ); meshSource->Update();
deformableModelFilter->SetInput( meshSource->GetOutput() ); deformableModelFilter->SetGradient( gradientMapFilter->GetOutput() );
typedef itk::CovariantVector<double, 2> StiffnessType;
StiffnessType stiffnessVector; stiffnessVector[0] = stiffness; stiffnessVector[1] = stiffness;
deformableModelFilter->SetTimeStep( timeStep ); deformableModelFilter->SetStiffness( stiffnessVector );
deformableModelFilter->SetStepThreshold( numberOfIterations ); deformableModelFilter->SetGradientMagnitude( externalForceScale );
try {
deformableModelFilter->Update();
}
catch( itk::ExceptionObject & excep ) { std::cerr << "Exception Caught !" << std::endl; std::cerr << excep << std::endl; }
//writer to *.vtk file
typedef itk::PointSetToImageFilter<MeshType,ImageType> MeshFilterType;
MeshFilterType::Pointer meshFilter = MeshFilterType::New();
meshFilter->SetInput(meshSource->GetOutput());
//(I confused in there !!!!! or )
try
{
meshFilter->Update();
}
catch( itk::ExceptionObject & excep )
{
std::cerr << "Exception Caught !" << std::endl;
std::cerr << excep << std::endl;
}
typedef itk::ImageFileWriter< ImageType > WriterType;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName( argv[2] );
writer->SetInput( meshFilter->GetOutput() );
std::cout << "Writing the image as " << std::endl << std::endl;
std::cout << argv[2] << std::endl << std::endl;
typedef itk::VTKImageIO ImageIOType;
ImageIOType::Pointer vtkIO = ImageIOType::New();
vtkIO->SetFileTypeToASCII();
writer->SetImageIO( vtkIO );
try
{
writer->Update();
}
catch (itk::ExceptionObject &ex)
{
std::cout << ex << std::endl;
return EXIT_FAILURE;
}
}
catch (itk::ExceptionObject &ex)
{
std::cout << ex << std::endl;
return EXIT_FAILURE;
}
}
//end list=========================================================================
Regards,
Sutrisno
__________________________________________________
Apakah Anda Yahoo!?
Lelah menerima spam? Surat Yahoo! memiliki perlindungan terbaik terhadap spam
http://id.mail.yahoo.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20100413/e76604d8/attachment-0001.htm>
More information about the Insight-users
mailing list