[Insight-users] translating image volume using affine transformation
michiel mentink
michael.mentink at st-hughs.ox.ac.uk
Mon Nov 30 06:18:33 EST 2009
FIXED IT!
found the problem explained before (sorry, can't find that particular post
back, since I closed the window and don't
recall the search term.
The problem was that the direction of the image wasn't copied to the filter,
instead the filter got a identity direction.
I solved it by adding the lines:
ImageType::DirectionType direction2 = reader->GetOutput()->GetDirection();
and subsituting direction by direction2 in
filter->SetOutputDirection(direction2 );
Now it works!
greets, Michael
On Mon, Nov 30, 2009 at 10:10 AM, michiel mentink <
michael.mentink at st-hughs.ox.ac.uk> wrote:
> To check this, I also followed Luiz Ibanez' suggestion to John Drozd in the
> message "problem with orientation of
> dicom output from segmentation.
>
> I pasted in the following lines into the codes:
>
> std::cout << reader->GetOutput()->GetDirection() << std::endl;
> std::cout << filter->GetOutput()->GetDirection() << std::endl;
>
> and got the following result:
>
> 0.997256 -0.0188216 -0.0716039
> 0.0740362 0.253523 0.964492
> -8.34627e-10 -0.967146 0.254221
>
> 1 0 0
> 0 1 0
> 0 0 1
>
> Does this mean for certain that the filter changes the axis' orientation?
>
> Kind regards,
>
> Michael
>
>
>
> On Mon, Nov 30, 2009 at 9:53 AM, michiel mentink <
> michael.mentink at st-hughs.ox.ac.uk> wrote:
>
>> might this be the same problem as in
>> problem with orientation of dicom output from segmentation where the
>> dicom 3d volume is changed in orientation?
>>
>> Kind regards, Michael
>>
>>
>> On Mon, Nov 30, 2009 at 9:42 AM, michiel mentink <
>> michael.mentink at st-hughs.ox.ac.uk> wrote:
>>
>>> 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/747af69b/attachment.htm>
More information about the Insight-users
mailing list