<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>