[prev in list] [next in list] [prev in thread] [next in thread]
List: insight-users
Subject: [Insight-users] Fwd: translating image volume using affine
From: michiel mentink <michael.mentink () st-hughs ! ox ! ac ! uk>
Date: 2009-11-30 9:42:21
Message-ID: 5488ad6a0911300142p39543a39y266311c7ef53994 () mail ! gmail ! com
[Download RAW message or body]
[Attachment #2 (multipart/alternative)]
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@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@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();
*/
}
[Attachment #5 (text/html)]
<div class="gmail_quote">Dear Bill,<br><br>I indeed didn't define the origin \
correctly. However, now there seems to be a more<br>serious problem.<br><br>Since \
I'm not sure I can upload screenshots to this mailing list, I have put a <br>
ITKsnap screenshot of the original dcm volume view, together with output after<br>
the affine filter and code on this site:<br><br><a \
href="http://sites.google.com/site/michielmentink/programming/itk/image-translate" \
target="_blank">http://sites.google.com/site/michielmentink/programming/itk/image-translate</a><br>
<br>
You can see there that before and after the affine filter, using identity \
transformation,<br>the filesize, number of pixels, and pixel sizes are \
identical.<br><br>There is a small bit of the original volume visible in the \
resulting dicom volume, but<br>
it contains only image data on a eight slices and the slices are deformed. \
<br><br>What I can make of the data, it seems that the filter has rotated the dicom \
volume<br>so that the 'coronal' (seeing)from the front) view now shows the \
'axial' view (seeing<br>
from above).<br>Also, the image is flipped in both X and Y axis.<br><br>If I change \
the Z origin to -69 instead of 69, most image slices actually contain the<br>image \
data, but still the views have been changed completely, like mentioned before.<br>
<br>What is going on? <br>What I really want to do is load in a dicom volume, \
translate it in the X, Y, Z direction<br>and save it. Normally, the affine filter \
should be used for this, right?<br> <br>For completeness, I've pasted the code \
again.<br><br>Kind regards, Michael<div class="im"><br><br><br><br><br><div \
class="gmail_quote">On Fri, Nov 27, 2009 at 4:05 PM, Bill Lorensen <span \
dir="ltr"><<a href="mailto:bill.lorensen@gmail.com" \
target="_blank">bill.lorensen@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); \
margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">What is the origin, spacing and \
direction of the input volume?<br> <div><div></div><div><br>
On Fri, Nov 27, 2009 at 9:59 AM, michiel mentink<br>
<<a href="mailto:michael.mentink@st-hughs.ox.ac.uk" \
target="_blank">michael.mentink@st-hughs.ox.ac.uk</a>> wrote:<br> > Dear \
all,<br> ><br>
> I'm trying to make a simple function that translates a dicom volume by a<br>
> certain amount.<br>
> It first opens all dicom slides of an mri scan.<br>
><br>
> To do this, I pasted DicomSeriesReadImageWrite2.cxx and<br>
> ResampleImageFilter.cxx and adjusted<br>
> the latter so it accepts 3 dimensions instead of two and uses integer<br>
> instead of unsigned char.<br>
><br>
> The function compiles and a .vtk volume is written. However, there's no<br>
> image data in the file.<br>
> It does display all pixels with value 100, as I set by:<br>
> filter->SetDefaultPixelValue( 100 );<br>
><br>
> questions:<br>
> 1) am I too naive to think that the ResampleImageFilter can easily be<br>
> extended from 2 to 3 dimensions?<br>
> 2) if not, what did I do wrong?<br>
><br>
</div></div></blockquote></div><br><br></div><div><div></div><div class="h5">#if \
defined(_MSC_VER)<br>#pragma warning ( disable : 4786 )<br>#endif<br><br>#ifdef \
__BORLANDC__<br>#define ITK_LEAN_AND_MEAN<br>#endif<br><br>
#include "itkOrientedImage.h"<br>
#include "itkGDCMImageIO.h"<br>#include \
"itkGDCMSeriesFileNames.h"<br>#include \
"itkImageSeriesReader.h"<br>#include \
"itkImageFileWriter.h"<br>#include "itkResampleImageFilter.h"<br>
<br>#include "itkAffineTransform.h"<br>#include \
"itkNearestNeighborInterpolateImageFunction.h"<br>// Software Guide : \
EndCodeSnippet<br><br>int main( int argc, char* argv[] )<br>{<br><br> if( argc < \
3 )<br>
{<br> std::cerr << "Usage: " << std::endl;<br> \
std::cerr << argv[0] << " DicomDirectory outputFileName \
[seriesName]" <br> << std::endl;<br> return \
EXIT_FAILURE;<br>
}<br><br> typedef signed short PixelType;<br> const unsigned int \
Dimension = 3;<br><br> typedef itk::OrientedImage< PixelType, Dimension > \
ImageType;<br><br> typedef itk::ImageSeriesReader< ImageType > \
ReaderType;<br>
ReaderType::Pointer reader = ReaderType::New();<br><br> typedef itk::GDCMImageIO \
ImageIOType;<br> ImageIOType::Pointer dicomIO = ImageIOType::New();<br> <br> \
reader->SetImageIO( dicomIO );<br><br> typedef itk::GDCMSeriesFileNames \
NamesGeneratorType;<br>
NamesGeneratorType::Pointer nameGenerator = NamesGeneratorType::New();<br><br> \
nameGenerator->SetUseSeriesDetails( true );<br> \
nameGenerator->AddSeriesRestriction("0008|0021" );<br><br> \
nameGenerator->SetDirectory( argv[1] );<br>
<br> try<br> {<br> std::cout << std::endl << "The directory: \
" << std::endl;<br> std::cout << std::endl << argv[1] \
<< std::endl << std::endl;<br> std::cout << "Contains the \
following DICOM Series: ";<br>
std::cout << std::endl << std::endl;<br><br> typedef \
std::vector< std::string > SeriesIdContainer;<br> <br> const \
SeriesIdContainer & seriesUID = nameGenerator->GetSeriesUIDs();<br>
<br> SeriesIdContainer::const_iterator seriesItr = seriesUID.begin();<br> \
SeriesIdContainer::const_iterator seriesEnd = seriesUID.end();<br> while( \
seriesItr != seriesEnd )<br> {<br> std::cout << seriesItr->c_str() \
<< std::endl;<br>
seriesItr++;<br> }<br><br> std::string seriesIdentifier;<br><br> if( \
argc > 3 ) // If no optional series identifier<br> {<br> \
seriesIdentifier = argv[3];<br> }<br> else<br> {<br>
seriesIdentifier = seriesUID.begin()->c_str();<br> }<br><br> \
std::cout << std::endl << std::endl;<br> std::cout << "Now \
reading series: " << std::endl << std::endl;<br>
std::cout << seriesIdentifier << std::endl;<br> std::cout << \
std::endl << std::endl;<br><br> typedef std::vector< std::string > \
FileNamesContainer;<br> FileNamesContainer fileNames;<br>
<br> fileNames = nameGenerator->GetFileNames( seriesIdentifier );<br> \
reader->SetFileNames( fileNames );<br> try<br> {<br> \
reader->Update();<br> }<br> catch (itk::ExceptionObject &ex)<br>
{<br> std::cout << ex << std::endl;<br> return \
EXIT_FAILURE;<br> }<br> }<br> catch (itk::ExceptionObject &ex)<br> \
{<br> std::cout << ex << std::endl;<br> return EXIT_FAILURE;<br>
}<br>// return EXIT_SUCCESS;<br><br>//////////////////////////////////////////////////////////////////////////////////////<br></div></div>// \
FROM HERE NEW CODE<div \
class="im"><br>/////////////////////////////////////////////////////////////////////////////////////<br>
<br>// const unsigned int Dimension = 3;<br>// typedef unsigned char \
InputPixelType;<br>// typedef unsigned char OutputPixelType;<br> typedef \
itk::Image< PixelType, Dimension > InputImageType;<br>
</div><div class="im">
typedef itk::Image< PixelType, Dimension > OutputImageType;<br></div><div \
class="im"> typedef itk::ImageFileWriter< ImageType > WriterType;<br> \
WriterType::Pointer writer = WriterType::New();<br><br> writer->SetFileName( \
argv[2] );<br>
<br></div><div class="im"> typedef \
itk::ResampleImageFilter<ImageType,ImageType> FilterType;<br> \
FilterType::Pointer filter = FilterType::New();<br><br> typedef \
itk::AffineTransform< double, Dimension > TransformType;<br>
TransformType::Pointer transform = TransformType::New();<br>
filter->SetTransform( transform );<br><br><br> typedef \
itk::NearestNeighborInterpolateImageFunction<ImageType, double > \
InterpolatorType;<br> InterpolatorType::Pointer interpolator = \
InterpolatorType::New();<br>
filter->SetInterpolator( interpolator );<br> <br> \
filter->SetDefaultPixelValue( 100 ); // if pixel in output image not defined, what \
value?<br><br> double spacing[ Dimension ];<br> spacing[0] = 0.29; // pixel spacing \
in millimeters along X<br>
spacing[1] = 0.29; // pixel spacing in millimeters along Y<br></div> spacing[2] = \
3.3; // pixel spacing in millimeters along Z<div class="im"><br><br> \
filter->SetOutputSpacing( spacing );<br><br> double origin[ Dimension ];<br>
</div> origin[0] = -174.22; // X space coordinate of origin<br>
origin[1] = -125.00; // Y space coordinate of origin<br> origin[2] = 69.37; // \
Z space coordinate of origin<div class="im"><br><br> filter->SetOutputOrigin( \
origin );<br> ImageType::DirectionType direction;<br>
direction.SetIdentity();<br>
filter->SetOutputDirection( direction );<br> ImageType::SizeType \
size;<br><br> size[0] = 560; // number of pixels along X<br> size[1] = 560; // \
number of pixels along Y<br></div> size[2] = 26;<div class="im">
<br><br> filter->SetSize( size );<br>
filter->SetInput( reader->GetOutput() );<br><br> writer->SetInput( \
filter->GetOutput() );<br> writer->Update();<br><br></div> std::cout << \
"Writing the image" << std::endl << std::endl;<br>
<br><br> /* (this code not executed)<div class="im"><br> \
TransformType::OutputVectorType translation;<br> translation[0] = 0; // X \
translation in millimeters<br> translation[1] = 0; // Y translation in \
millimeters<br>
</div> translation[2] = 0; // Z translation in millimeters<div class="im"><br>
<br> transform->Translate( translation );<br> // Software Guide : \
EndCodeSnippet<br><br> writer->SetInput( filter->GetOutput() );<br><br> \
std::cout << "Writing the image as " << std::endl << \
std::endl;<br>
std::cout << argv[2] << std::endl << std::endl;<br> \
writer->Update();<br></div> */<br><br>}<br> </div><br>
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.html
Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ
Follow this link to subscribe/unsubscribe:
http://www.itk.org/mailman/listinfo/insight-users
[prev in list] [next in list] [prev in thread] [next in thread]
Configure |
About |
News |
Add a list |
Sponsored by KoreLogic