[vtkusers] Announcement: New image rendering classes for VTK

[David Gobbi]

The best way to work out coordinate transformation is to consider what
happens to the position of a single pixel (I,J,K) in the image.

DICOM transformations go like this: first, there is a stack of images where
each pixel identified by column and row, i.e. by (I,J,0) indices (for now,
I'm just considering a single slice so the Z index is zero).  These are then
multiplied by the column and row spacing, to get (I*sx, J*sy, 0).  Then they
are multiplied by the Orientation matrix, which is constructed by creating a
3x3 matrix where the first two columns are from the ImageOrientationPatient
tag in the header, and the third column (usually) the cross product of the
first two columns.  Finally, the last step in creating (x,y,z) patient
coordinates is to add the ImagePositionPatient to the stretched and rotated
(I,J,0) pixel indices.  So, that is how patient coordinates are computed for
DICOM images.

In VTK, the (I,J,K) indices of the voxels are called the "structured coords"
of the image.  To convert these to VTK's "data coords", i.e. what VTK uses
for polydata and for the mappers, the IJK indices have to be multiplied by
the image Spacing (just like they are multiplied by the spacing for DICOM)
and then the Origin is added to these values to get VTK "data coordinates."
 Finally, the data coordinates are multiplied by a 4x4 matrix (the
vtkProp3D's matrix, i.e. the matrix that is composed of the position and
rotation of the vtkVolume).  This results in VTK "world coordinates".

In summary:
DICOM -> scale by Spacing, rotate by Orientation, translate by Position
VTK -> scale by Spacing, translate by Origin, then rotate and translate by a
4x4 matrix

There is nothing in DICOM that corresponds to the Origin of a VTK image
data, i.e. in DICOM the scale operation is followed directly by the
rotation.  So when dealing with pure DICOM patient coords in VTK, the image
Origin should generally be set to (0,0,0).  Furthermore, the DICOM
Orientation and Position become subsumed into a singe 4x4 VTK matrix
transformation.  The 4x4 matrix is constructed from the 3x3 DICOM
orientation matrix, and the last column is set to the DICOM Position.

Now what I said about the Origin always being set to zero... that represents
the "purest" way of expressing DICOM patient coords in VTK.  But it is not
the only way.  Another choice would be to take the DICOM Position, multiply
it by the inverse of the orientation matrix, and then set the Origin to
that.  In that case, the "translation" part of the VTK 4x4 matrix would be
zero.

Finally, getting to your question about calling SetOrigin() and
SetPosition() on the vtkVolume.  The vtkVolume, as I discussed, has an
internal 4x4 matrix (that you can get by calling volume->GetMatrix()), and
this matrix is constructed from the Origin, Orientation, and Position of the
vtkVolume.  The order of operations here is that the 4x4 matrix is
constructed by translating by the _negative_ of the origin, rotating by the
orientation, and translating once again by the position.  But the
Orientation is a nasty issue here... there is no easy way to figure out the
angles you would need to use for SetOrientation to achieve the correct
ImageOrientationPatient for your DICOM.  However, it is straightforward to
directly construct a 4x4 matrix as I discussed above.  So my recommendation
is to _never, ever_ use the SetOrigin(), SetOrientation(), and SetPosition()
methods when visualizing DICOM data.  Instead, build the 4x4
rotation+position matrix yourself, and use SetUserMatrix() to attach it to
your volume.

The reason that SetOrigin() and SetPosition() worked for you is that, by
setting both to the same value, you end up changing the center of rotation.
 But heed my advice: if you want to do things quantitatively, the best way
to do so is by building the 4x4 transformation matrix yourself.  Letting VTK
build it for you is IMHO just adding a lot of extra unknowns into the mix.

 - David




2011/7/5 Dženan Zukić <dzenanz at gmail.com>:
> I solved this (see my other mail). But I still don't know why I had to set
> origin and position:
> volume->SetOrigin(-origin[0], -origin[1], -origin[2]);
> volume->SetPosition(-origin[0], -origin[1], -origin[2]);
> Can you explain it? Also, a link to some text explaining vtk's coordinate
> transformation system for future reference would be very appreciated (if
> such a text exists). I had been searching but could not find a source
which
> explains in which order are the transformation matrices applied, in which
> space are voxels of a vtkVolume etc.
> Regards,
> Dženan
> 2011/7/1 Dženan Zukić <dzenanz at gmail.com>
>>
>> I did the math, and it should be as simple as this:
>>
>> VisualizingImageType::DirectionType d=logic->visualizing->GetDirection();
>> vtkMatrix4x4 *mat=vtkMatrix4x4::New(); //identity matrix
>> for (int i=0; i<3; i++)
>>     for (int k=0; k<3; k++)
>>         mat->SetElement(i,k, d(i,k));
>> VisualizingImageType::PointType origin=logic->visualizing->GetOrigin();
>> for (int i=0; i<3; i++)
>>     mat->SetElement(i,3, origin[i]);
>> //mat->Invert(); //inverse produces even wrong orientation
>> volume->SetUserMatrix(mat);
>> However this does not work. This is a simple transformation from scaled
>> index space (vtkVolume handles spacing) to DICOM physical space. If no
one
>> has any suggestion, I will simply have to create a system of linear
>> equations from which to construct the VTK's transform matrix. And
equations
>> I will have to obtain by passing a number of points through ITK's
>> transformation and add some simple objects to the vtk scene at that spot
and
>> on index space spot but with VTK transformation applied.
>> Regards,
>> Dženan
>> 2011/6/29 Dženan Zukić <dzenanz at gmail.com>
>>>
>>> I tried trial and error: premultiply translation,
postmultiply, divide by
>>> and multiply by spacing, all that unsuccessfully. Now I need to work out
the
>>> math :D
>>>
>>> 2011/6/29 David Gobbi <david.gobbi at gmail.com>
>>>>
>>>> Hi Dzenan,
>>>>
>>>> For your code, my guess is that the "translation" part of the 4x4
>>>> matrix should be set to something like the origin multiplied by the
>>>> rotation matrix, but there might be other little details that you need
>>>> to take care of.  My only suggestion is that you don't try to fix it
>>>> through trial-end-error.  Use pencil-and-paper to work through the
>>>> math.
>>>>
>>>>  - David
>>>>
>>>> 2011/6/29 Dženan Zukić <dzenanz at gmail.com>:
>>>> > Hi David,
>>>> >
>>>> > with this I try to do what you suggested, i.e. consider VTK's world
>>>> > coordinate system as DICOM patient coordinate system:
>>>> > //red is the transformation-related code
>>>> > void MainWindow::updateVisualization()
>>>> > {
>>>> >     typedef itk::ImageToVTKImageFilter<VisualizingImageType>
>>>> > itkVtkConverter;
>>>> >     itkVtkConverter::Pointer conv=itkVtkConverter::New();
>>>> >     conv->SetInput(logic->visualizing);
>>>> >     vtkGPUVolumeRayCastMapper *mapper =
>>>> > vtkGPUVolumeRayCastMapper::New();
>>>> >     mapper->SetInput(conv->GetOutput());
>>>> >     if (volume)
>>>> >         volume->Delete();
>>>> >     volume=vtkVolume::New();
>>>> >     volume->SetProperty( myTransferFunction );
>>>> >     volume->SetMapper( mapper );
>>>> >     mapper->SetBlendModeToComposite();
>>>> >     VisualizingImageType::DirectionType
>>>> > d=logic->visualizing->GetDirection();
>>>> >     vtkMatrix4x4 *mat=vtkMatrix4x4::New();
>>>> >     for (int i=0; i<3; i++)
>>>> >         for (int k=0; k<3; k++)
>>>> >             mat->SetElement(i,k, d(i,k));
>>>> >     mat->SetElement(3,3, 1);
>>>> >     VisualizingImageType::SpacingType sp =
>>>> > logic->visualizing->GetSpacing();
>>>> >     VisualizingImageType::PointType
>>>> > origin=logic->visualizing->GetOrigin();
>>>> >     for (int i=0; i<3; i++)
>>>> >         mat->SetElement(i,3, origin[i]/sp[i]);
>>>> >     volume->SetUserMatrix(mat); //orientation and size OK, position
>>>> > wrong
>>>> >     vtkRenderer *renderer =
>>>> > vis->GetRenderWindow()->GetRenderers()->GetFirstRenderer();
>>>> >     renderer->RemoveAllViewProps();
>>>> >     renderer->AddVolume( volume );
>>>> >     defaultCameraPos(); //centers view on the volume, looking at it
>>>> > from
>>>> > left
>>>> >     vis->GetRenderWindow()->Render();
>>>> >
>>>> >
>>>> >
vis->GetRenderWindow()->GetRenderers()->GetFirstRenderer()->ResetCamera();
>>>> > }
>>>> > The polygonal data is in DICOM's patient coordinate system (vertex
>>>> > positions
>>>> > from itk::Image->TransformIndexToPhysicalPoint() etc). The
orientation
>>>> > is
>>>> > correct, but I can't get the position correctly. I have tried
>>>> > -origin[i]/sp[i], origin[i]*sp[i], origin[i], -origin[i],
-2*origin[i]
>>>> > and
>>>> > similar combinations but none of them worked. Any suggestion to what
>>>> > am I
>>>> > doing wrong?
>>>> > Regards,
>>>> > Dženan
>>>> > 2011/6/29 David Gobbi <david.gobbi at gmail.com>
>>>> >>
>>>> >> On Tue, Jun 28, 2011 at 9:40 AM, Dženan Zukić <dzenanz at gmail.com>
>>>> >> wrote:
>>>> >> > Hi David,
>>>> >> >
>>>> >> > I am interested in this too:
>>>> >> > Does the above approach affect the polygonal actors in the scene
>>>> >> > (segmented
>>>> >> > parts of the image)?
>>>> >> > Regards,
>>>> >> > Dženan
>>>> >>
>>>> >> When the DICOM patient coordinate system is used as the VTK world
>>>> >> coordinate system, the assumption is that everything (images AND
>>>> >> polydata) would have a either already be in the patient coordinate
>>>> >> system, or if not, you would to have a 4x4 transform to set as the
>>>> >> actor's UserMatrix (or UserTransform).  I.e. a transform to bring
the
>>>> >> data from its original coordinate system into the patient coordinate
>>>> >> system.
>>>> >>
>>>> >>  - David
>>>> >
>>>> >
>>>
>>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.vtk.org/pipermail/vtkusers/attachments/20110705/18b91713/attachment.htm>