[Insight-users] Fwd: translating image volume using affine transformation
michiel mentink
michael.mentink at st-hughs.ox.ac.uk
Mon Nov 30 04:42:21 EST 2009
Dear Bill,
I indeed didn't define the origin correctly. However, now there seems to be
a more
serious problem.
Since I'm not sure I can upload screenshots to this mailing list, I have put
a
ITKsnap screenshot of the original dcm volume view, together with output
after
the affine filter and code on this site:
http://sites.google.com/site/michielmentink/programming/itk/image-translate
You can see there that before and after the affine filter, using identity
transformation,
the filesize, number of pixels, and pixel sizes are identical.
There is a small bit of the original volume visible in the resulting dicom
volume, but
it contains only image data on a eight slices and the slices are deformed.
What I can make of the data, it seems that the filter has rotated the dicom
volume
so that the 'coronal' (seeing)from the front) view now shows the 'axial'
view (seeing
from above).
Also, the image is flipped in both X and Y axis.
If I change the Z origin to -69 instead of 69, most image slices actually
contain the
image data, but still the views have been changed completely, like mentioned
before.
What is going on?
What I really want to do is load in a dicom volume, translate it in the X,
Y, Z direction
and save it. Normally, the affine filter should be used for this, right?
For completeness, I've pasted the code again.
Kind regards, Michael
On Fri, Nov 27, 2009 at 4:05 PM, Bill Lorensen <bill.lorensen at gmail.com>wrote:
> What is the origin, spacing and direction of the input volume?
>
> On Fri, Nov 27, 2009 at 9:59 AM, michiel mentink
> <michael.mentink at st-hughs.ox.ac.uk> wrote:
> > Dear all,
> >
> > I'm trying to make a simple function that translates a dicom volume by a
> > certain amount.
> > It first opens all dicom slides of an mri scan.
> >
> > To do this, I pasted DicomSeriesReadImageWrite2.cxx and
> > ResampleImageFilter.cxx and adjusted
> > the latter so it accepts 3 dimensions instead of two and uses integer
> > instead of unsigned char.
> >
> > The function compiles and a .vtk volume is written. However, there's no
> > image data in the file.
> > It does display all pixels with value 100, as I set by:
> > filter->SetDefaultPixelValue( 100 );
> >
> > questions:
> > 1) am I too naive to think that the ResampleImageFilter can easily be
> > extended from 2 to 3 dimensions?
> > 2) if not, what did I do wrong?
> >
>
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#ifdef __BORLANDC__
#define ITK_LEAN_AND_MEAN
#endif
#include "itkOrientedImage.h"
#include "itkGDCMImageIO.h"
#include "itkGDCMSeriesFileNames.h"
#include "itkImageSeriesReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkAffineTransform.h"
#include "itkNearestNeighborInterpolateImageFunction.h"
// Software Guide : EndCodeSnippet
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 ) // If no optional series identifier
{
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;
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;
}
}
catch (itk::ExceptionObject &ex)
{
std::cout << ex << std::endl;
return EXIT_FAILURE;
}
// return EXIT_SUCCESS;
//////////////////////////////////////////////////////////////////////////////////////
// FROM HERE NEW CODE
/////////////////////////////////////////////////////////////////////////////////////
// const unsigned int Dimension = 3;
// typedef unsigned char InputPixelType;
// typedef unsigned char OutputPixelType;
typedef itk::Image< PixelType, Dimension > InputImageType;
typedef itk::Image< PixelType, Dimension > OutputImageType;
typedef itk::ImageFileWriter< ImageType > WriterType;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName( argv[2] );
typedef itk::ResampleImageFilter<ImageType,ImageType> FilterType;
FilterType::Pointer filter = FilterType::New();
typedef itk::AffineTransform< double, Dimension > TransformType;
TransformType::Pointer transform = TransformType::New();
filter->SetTransform( transform );
typedef itk::NearestNeighborInterpolateImageFunction<ImageType, double >
InterpolatorType;
InterpolatorType::Pointer interpolator = InterpolatorType::New();
filter->SetInterpolator( interpolator );
filter->SetDefaultPixelValue( 100 ); // if pixel in output image not
defined, what value?
double spacing[ Dimension ];
spacing[0] = 0.29; // pixel spacing in millimeters along X
spacing[1] = 0.29; // pixel spacing in millimeters along Y
spacing[2] = 3.3; // pixel spacing in millimeters along Z
filter->SetOutputSpacing( spacing );
double origin[ Dimension ];
origin[0] = -174.22; // X space coordinate of origin
origin[1] = -125.00; // Y space coordinate of origin
origin[2] = 69.37; // Z space coordinate of origin
filter->SetOutputOrigin( origin );
ImageType::DirectionType direction;
direction.SetIdentity();
filter->SetOutputDirection( direction );
ImageType::SizeType size;
size[0] = 560; // number of pixels along X
size[1] = 560; // number of pixels along Y
size[2] = 26;
filter->SetSize( size );
filter->SetInput( reader->GetOutput() );
writer->SetInput( filter->GetOutput() );
writer->Update();
std::cout << "Writing the image" << std::endl << std::endl;
/* (this code not executed)
TransformType::OutputVectorType translation;
translation[0] = 0; // X translation in millimeters
translation[1] = 0; // Y translation in millimeters
translation[2] = 0; // Z translation in millimeters
transform->Translate( translation );
// Software Guide : EndCodeSnippet
writer->SetInput( filter->GetOutput() );
std::cout << "Writing the image as " << std::endl << std::endl;
std::cout << argv[2] << std::endl << std::endl;
writer->Update();
*/
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20091130/42c48d29/attachment.htm>
More information about the Insight-users
mailing list