[vtkusers] using itk transform in vtk

[David Gobbi]

I haven't looked into gdcmreslice myself, but the procedure that you
have outlined is very different from what I have used with DICOM
files.  What I do is as follows:

1) Read the image as-is, i.e. do not allow the reader to flip it or do
any other adjustments.  Computing the vtkImageData Origin is tricky,
this Origin is in VTK's "data coordinate space" for the vtkImageData,
while DICOM's ImagePositionPatient is in DICOM patient coordinates.
So I have to use the direction cosines matrix (see step 3) to convert
ImagePositionPatient to the vtkImageData Origin that I use.

2) Compute the Z spacing by measuring the distance between adjacent
DICOM slices, and checking the direction of the cross product of x, y
direction cosines from ImageOrientationPatient to make sure the slices
are in the right order

3) Use the direction cosines to create an "image orientation matrix".
I do not re-order the axes of the image with vtkImageReslice (as far
as I'm concerned, doing so just complicates the whole process, and
cannot account for oblique orientations)

So at this point, I have a DICOM image with X as the row direction
(left-to-right), Y as the column direction (top-to-bottom), and Z as
the slice direction.  I also have a vtkMatrix4x4 that is a
transformation matrix from this "raw" coordinate system to the DICOM
patient coordinate system.  Then, I identify one of my DICOM data sets
as my "primary" data set, and identify its DICOM patient coordinate
system as being identical to my VTK world coordinate system.

When doing registrations, I use a vtkTransform that converts from the
DICOM patient coordinates of one image, to the DICOM patient
coordinates of another image.  Each of these images also has a
vtkMatrix4x4 that contains its direction cosines, which relate the
vtkImageData data coordinates to the DICOM patient coordinates of that
image.

By doing the above, I keep the various coordinate spaces as
compartmentalizes as possible.  I avoid flipping the image or
resampling it along a new set of axes because that tends to make
things more confusing, I prefer to just store a direction cosines
matrix for each image and leave the data alone.

  David


On Thu, Nov 18, 2010 at 2:25 PM, Karl <bulkmailaddress at gmail.com> wrote:
> 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
>
> _______________________________________________
> 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
>