[vtkusers] using itk transform in vtk

Thu Nov 18 16:25:49 EST 2010

Mathieu + vtk community,
Thanks for replying.  I looked into gdcmrelice.cxx example today.  I was
unsuccessful in using that approach to achieve results.  In case it matters
I have gdcm 2.0.16 and vtk 5.4.2.

If I understand correctly the following process is applied:
1)Load image with vtkGDCMImageReader.  Due to differences in origin the
reader shifts the origin along the y-axis (based on the directional cosines)
and reverses the order of the slices along the y direction.

2)This reversal of order along the y-axis also results in the x-axis
appearing flipped, so the vtkImageFlip filter is used to undo that effect.

3)vtkImageReslice is then used to reorientate the data to lie along the
actual dicom x,y,z axes instead of the row, column, slice sampling of the
data set.  This doesn't seem the best approach because it results in
obliquely resampling the entire data set.
Using this approach does not work correctly.  The data does appear with RAH
orientation but does not get placed correctly.

My code goes something like this:
vtkGDCMImageReader *reader = vtkGDCMImageReader::New();
reader->SetFileName( filename );

vtkImageFlip *flip = vtkImageFlip::New();
flip->SetInput(reader ->GetOutput());
flip->SetFilteredAxis(0);
flip->Update();

vtkImageReslice *reslice = vtkImageReslice::New();
reslice->SetInput(flip->GetOutput());
vtkMatrix4x4 *invert = vtkMatrix4x4::New();
invert->DeepCopy(reader->GetDirectionCosines() );
invert->Invert();
reslice->SetResliceAxes(invert);
reslice->Update();

vtkMarchingCubes* march = vtkMarchingCubes::New();
march->SetInputConnection(reslice->GetOutputPort());

// create transform matrix
vtkMatrix4x4* matrix = vtkMatrix4x4::New(); ...
vtkTransform* matrixTransform = vtkTransform::New();
matrixTransform->SetMatrix(matrix);

vtkTransformPolyDataFilter* transform = vtkTransformPolyDataFilter::New();
transform->SetTranform(matrixTransform);
transform->SetInputConnection(march->GetOutputPort());

vtkPolyDataMapper* mapper = vtkPolyDataMapper::New();
mapper->SetInputConnection(transform->GetOutputPort());

vtkActor* actor = vtkActor::New();
actor->SetMapper(mapper);

This is basically done twice.  Once with the original data and the transform
from ITK, and once with the data set that has already been transformed by
ITK and an identity transform.  The two objects do not appear in the same
location.  The dicom images do start out having different origins and
directional cosines.  But if both are mapped to physical space correctly
they should overlap exactly.

Thanks

-----Original Message-----
From: Mathieu Malaterre [mailto:mathieu.malaterre at gmail.com] 
Sent: Thursday, November 18, 2010 7:04 AM
To: bulkmailaddress at gmail.com
Subject: Re: [vtkusers] using itk transform in vtk

Have you tried the gdcmreslice.cxx example in GDCM ? Is it working for yuo ?

Thx

On Thu, Nov 18, 2010 at 6:20 AM, Karl <bulkmailaddress at gmail.com> wrote:
> Vtk folks and David,
>
> Thanks for responding.  I have yet to have success, but I believe I have
> addressed the issues you point out.
> I am using vtkGDCMImageReader with the FileLowerLeft flag set to 0.  If I
> understand correctly this importer correctly flips the Y and adjusts the
> origin.
>
> ITK transform is about the 0,0,0 physical coordinate.  When I load the
image
> in itk the origin information should be used by this importer to properly
> position the image so that the 0,0,0 physical coordinate corresponds to
the
> 0,0,0 vtk location.  The inverse of the itk transform can then be directly
> applied to the imported image.
>
> This still does not work however.  I need to somehow get the vtk
coordinate
> system to accurately represent the physical coordinate system so that
> objects from different dicom images, with different origins and
directional
> cosines, appear in proper relation to each other.
>
> Any suggestions on how to proceed?
>
> Thanks
>
> -----Original Message-----
> From: David Gobbi [mailto:david.gobbi at gmail.com]
> Sent: Tuesday, November 16, 2010 4:31 PM
> To: bulkmailaddress at gmail.com
> Cc: vtkusers at vtk.org
> Subject: Re: [vtkusers] using itk transform in vtk
>
> Hi Karl,
>
> The VTK transform operates on the physical origin, which can be at any
> desired location in the image.  The relationship between DICOM
> coordinates and VTK coordinates is explained thusly:
>
> DICOM has three fields related to the positions of the image voxels in
> physical (or patient) space:
> - ImagePositionPatient (position of upper right corner)
> - ImageOrientationPatient (row and column directions)
> - PixelSpacing (X and Y)
>
> VTK image data has two fields (no orientation):
> - Origin (position of lower right corner)
> - Spacing (X, Y, and Z)
>
> Note that the image data "Origin" is a misnomer, it is not the
> physical origin of the coordinate system.  It is just the position of
> the voxel in the lower left corner.  It is the responsibility of the
> VTK image reader that you are using to set the Origin and Spacing of
> the data correctly.  (And again, I will repeat that what VTK calls the
> image "Origin" is not the physical image origin, it is just a
> position).
>
> The fundamental difference between these two coordinate systems is not
> the location of the origin, it is a vertical flip since DICOM images
> are top-to-bottom while VTK images are bottom-to-top.  When dealing
> with DICOM, I set the VTK image "origin" to the DICOM
> ImagePositionPatient, and then apply a view transform that flips the
> image by using a rotation of 180 degrees about the X axis.
>
>  - David
>
> On Tue, Nov 16, 2010 at 2:02 PM, Karl <bulkmailaddress at gmail.com> wrote:
>> Hi,
>> I am trying to use an itk transform in vtk.  Here is a high level
overview
>> of what I do:
>>
>> 1) In itk use itk::ImageRegistrationMethod with itk::AffineTransform to
> find
>> a transform.
>>
>> 2) I use itk::ResampleImageFilter to create an image transformed using
the
>> itk::AffineTransform.
>>
>> 3) I use the GetLastTransformParameters() function to get the affine
>> transformation matrix to use in vtk as.
>> [param0 param1 param2 param9]
>> [param3 param4 param5 param10]
>> [param6 param7 param8 param11]
>> [0 0 0 1]
>>
>> 4) Now in vtk I try to display both the resampled image and the original
>> image with the transform applied.  They should show as exactly the same
if
> I
>> can apply the transform in vtk the same as it is applied in itk.
>>
>> 5) I load both images in vtk.  The images are binary so I use marching
> cubes
>> to create meshes of them.  I display the transformed image polydata as
is.
>>
>> 6) I apply a vtkTransform to the original image polydata.  I do this by
>> concatenating 3 vtkMatrix4x4. Two to adjust for different origins in itk
> and
>> vtk and one for the itk transform.
>> Translate_Back * itkTransform * Translate_itk_origin_to_vtk_origin
>>
>> Where
>> itkTransform is as shown in 3 above.
>>
>> Translate_itk_origin_to_vtk_origin is
>> [1 0 0 -originX]
>> [0 1 0 -originY]
>> [0 0 1 -originZ]
>> [0 0 0 1]
>> Where the origin is as found in the dicom header for the original image
>>
>> Translate_Back
>> [1 0 0 originX]
>> [0 1 0 originY]
>> [0 0 1 originZ]
>> [0 0 0 1]
>>
>> 7) Set vtkTransform inverse flag and transform and display the original
>> image.
>>
>> 8) The two objects show up at different places in the scene instead of at
>> the same location.
>>
>> If I understand this correctly the itk transform is to be applied about
> the
>> physical origin.  Not the location of the corner pixel or the center of
> the
>> image.
>>
>> The vtk transform uses an origin at the corner of the image.
>>
>> By translating by the origin defined in the dicom header I move the
> physical
>> origin to the corner pixel of the image.  This corrects for the
>> discrepancies between the point with itk and vtk apply their transforms.
>>
>> Also the itk transform goes from output to input, so the inverse
transform
>> needs to be applied in vtk.
>>
>> Can anyone see what I am missing?
>> Thanks
>
> _______________________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
>
> Please keep messages on-topic and check the VTK FAQ at:
http://www.vtk.org/Wiki/VTK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.vtk.org/mailman/listinfo/vtkusers
>

-- 
Mathieu