[vtkusers] using itk transform in vtk

[David Gobbi]

Hi Karl,

I'm glad to hear that you got it working.  I have the cosines going down the
columns of the matrix, so that the matrix transforms "raw data" coordinates
into patient coordinates.  In other words, option (1).

I was trying to remember why I chose that option... and I think this is the
reason: the Z direction cosine vector is not always orthogonal to the X and
Y cosine vectors.  In a tilted-gantry CT scan, the volume is rhombus-shaped,
and the Z direction cosine is the direction in which the gantry moves
between slices (i.e. not the cross product of the X and Y vectors).  In this
case, you will get the correct transform if you choose option (1), but will
get the wrong transform if you choose option (2).

Most DICOM scans provide an orthogonal volume, and for these it makes no
difference.  It is only with the rare non-orthogonal scan volumes that the
difference becomes apparent.

  David

On Mon, Nov 22, 2010 at 12:35 PM, Karl <bulkmailaddress at gmail.com> wrote:

> David,
>
> Thanks for the help.  I have gotten it to work.  It was an issue with a
> matrix being transposed when it shouldn't.  I have one final question.
>  When
> I load the directional cosines from the header should they go into the rows
> or columns of the directional cosine matrix.
>
> For example I get the same results if I:
> 1) load the cosines into columns of matrix, multiply origin by inverse, and
> multiply data set by cosine matrix
> 2) load the cosines into rows of matrix, multiply origin by matrix,
> multiply
> data set by inverse of cosine matrix
>
> I know these are the same thing mathematically, but which one is correct
> based on convention.  I believe the first one is correct.
>
> Thanks for all the help getting this to work.
>
> Karl
>
> -----Original Message-----
> From: vtkusers-bounces at vtk.org [mailto:vtkusers-bounces at vtk.org] On Behalf
> Of David Gobbi
> Sent: Friday, November 19, 2010 4:22 PM
> To: vtkusers at vtk.org
> Subject: Re: [vtkusers] using itk transform in vtk
>
> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.vtk.org/pipermail/vtkusers/attachments/20101122/65dc001c/attachment.htm>