[vtkusers] using itk transform in vtk

[David Gobbi]

Hi Karl,

I meant to reply to you this morning but was distracted by other tasks...

The DICOM code that I have (which is actually written in Python) does
things like this:

        matrix = vtk.vtkMatrix4x4()
        for j in range(3):
            matrix.SetElement(j,0,rowOrientation[j])
            matrix.SetElement(j,1,columnOrientation[j])
            matrix.SetElement(j,2,sliceOrientation[j])

The above code builds the direction cosines from the DICOM
orientation, the slice orientation is the cross product of the other
two orientations.  Then, to convert the DICOM position into the
ImageData origin, I use the following.  This code is so old it uses
Python Numeric... but I'll give it to you raw, and you can figure out
how to do the same thing with vtkMatrix4x4.

        origin = Numeric.array(position)
        dircos = Numeric.transpose(Numeric.array([rowOrientation,
columnOrientation, sliceOrientation]))
        origin = list(LinearAlgebra.solve_linear_equations(dircos,origin))

I apologize for the python code, but as you can see, your own code
more-or-less duplicates what I do.  So all that I can suggest is the
following:

1) Check the value that reader->GetDataOrigin() returns.  Make sure
that it really is the same as the ImagePositionPatient value in the
DICOM header.

2) Check the math, just to make sure that matrices aren't accidentally
transposed or any of the other mistakes that we often make.

  David


On Fri, Nov 19, 2010 at 12:00 PM, Karl <bulkmailaddress at gmail.com> wrote:
> David,
> Thanks for walking me through another way.  I have attempted to implement
> what you propose, but am not doing something quite right.  Would you mind
> reviewing my procedure below and pointing out where I go wrong.  I suspect
> it is with the origin.  The following code is called twice, once for the
> non-transformed data, and once for the transformed data, hopefully they show
> up in the same spot:)
>
> // load data
> vtkGDCMImageReader *reader = vtkGDCMImageReader::New();
> reader->SetFileName( filename );
> reader->SetFileLowerLeft(1); // do not recalculate origin or y-flip data
>
> // get directional cosines
> double directionalCosines[3][3];
> image->GetImageOrientationPatient(directionalCosines[0][0],directionalCosine
> s[0][1],
>
> directionalCosines[0][2],directionalCosines[1][0],directionalCosines[1][1],
>    directionalCosines[1][2]);
> directionalCosines[2][0] =
> directionalCosines[0][1]*directionalCosines[1][2]-directionalCosines[0][2]*d
> irectionalCosines[1][1];
> directionalCosines[2][1] =
> directionalCosines[0][2]*directionalCosines[1][0]-directionalCosines[0][0]*d
> irectionalCosines[1][2];
> directionalCosines[2][2] =
> directionalCosines[0][0]*directionalCosines[1][1]-directionalCosines[0][1]*d
> irectionalCosines[1][0];
>
> // create directional cosine matrix, transforms from raw data space to
> patient space
> vtkMatrix4x4* dc = vtkMatrix4x4::New();
> dc->Identity();
> dc->SetElement(0,0,directionalCosines[0][0]);
> dc->SetElement(0,1,directionalCosines[0][1]);
> dc->SetElement(0,2,directionalCosines[0][2]);
> dc->SetElement(1,0,directionalCosines[1][0]);
> dc->SetElement(1,1,directionalCosines[1][1]);
> dc->SetElement(1,2,directionalCosines[1][2]);
> dc->SetElement(2,0,directionalCosines[2][0]);
> dc->SetElement(2,1,directionalCosines[2][1]);
> dc->SetElement(2,2,directionalCosines[2][2]);
>
> // inverse directional cosine matrix, transform from patient space to raw
> data space
> vtkMatrix4x4* dcInverse = vtkMatrix4x4::New();
> dc->Invert(dc,dcInverse);
>
> // calculate new origin, transform from patient space to raw data space
> double origin[4];
> origin[3]=1;
> reader->GetDataOrigin(origin);
> dcInverse->MultiplyPoint(origin,origin);
>
> // apply new origin
> vtkImageChangeInformation* information = vtkImageChangeInformation::New();
> information->SetOutputOrigin(origin[0], origin[1], origin[2]);
> information->SetInputConnection(reader->GetOutputPort());
>
> // create polydata
> vtkMarchingCubes* march = vtkMarchingCubes::New();
> march->SetInputConnection(information->GetOutputPort());
>
> // transform polydata from raw data space to patient space
> vtkTransformPolyDataFilter* transformToPatientSpace =
> vtkTransformPolyDataFilter::New();
> vtkTransform* matrixTransform = vtkTransform::New();
> matrixTransform->SetMatrix(dc);
> transformToPatientSpace->SetTransform(matrixTransform);
> transformToPatientSpace->SetInputConnection(march->GetOutputPort());
>
> //*************************************
> // this section is only called for the non-transformed data
> // load and apply itk transform
> vtkMatrix4x4* itkMatrix =
> [param0 param1 param2 param9]
> [param3 param4 param5 param10]
> [param6 param7 param8 param11]
> [0 0 0 1] // from itk writing out matrix
> itkMatrix->Invert(itkMatrix,itkMatrix); // itk actually provides inverse
> transforms
>
> vtkTransform* matrixTransform = vtkTransform::New();
> matrixTransform->SetMatrix(itkMatrix);
>
> vtkTransformPolyDataFilter* itkTransform =
> vtkTransformPolyDataFilter::New();
> itkTransform->SetTranform(matrixTransform);
> itkTransform->SetInputConnection(transformToPatientSpace->GetOutputPort());
> **************************************//
>
> // display polydata
> vtkPolyDataMapper* mapper = vtkPolyDataMapper::New();
> mapper->SetInputConnection(transformToPatientSpace ->GetOutputPort());
>
> vtkActor* actor = vtkActor::New();
> actor->SetMapper(mapper);
>
>
> I believe I am accomplishing the following.
> Load the data already transformed by itk (data A).  Apply transform to
> change it from raw data space to patient space.  Display data.   VTK space
> now corresponds to this patient space.
>
> Load untransformed data (data B).  Apply transform to change it from its raw
> data space to its patient space.  Apply inverse itk transform to change it
> from its patient space to data A's patient space.  Display data.
>
> Thanks
>
>
> -----Original Message-----
> From: vtkusers-bounces at vtk.org [mailto:vtkusers-bounces at vtk.org] On Behalf
> Of David Gobbi
> Sent: Thursday, November 18, 2010 4:55 PM
> To: vtkusers at vtk.org
> Subject: Re: [vtkusers] using itk transform in vtk
>
> 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
>>
> _______________________________________________
> 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
>
>