[vtkusers] Announcement: New image rendering classes for VTK
Dženan Zukić
dzenanz at gmail.com
Tue Jul 5 16:35:12 EDT 2011
Take a look at the example:
http://vtk.org/Wiki/VTK/Examples/Cxx/VolumeRendering/itkVtkImageConvert
If I call GetMatrix before the line "volume->SetOrigin(...)", it returns
identity. If I call it after SetPosition, it contains the translation part.
I usually convert DICOM data to MetaImage format which is more convenient,
but still preserves the crucial metadata. For the image in the screenshot I
manually edited the image and set an orientation 45 degrees rotated to
coordinate frame, because the orientations usually do not visibly deviate
from coordinate axes. The point of ITK is to deliver these values to us.
When I get to university tomorrow I will send you the image in question with
some accompanying polydata so you can also experiment yourself. The data
size is above the list limit, so I will not send it to the list.
2011/7/5 David Gobbi <david.gobbi at gmail.com>
> What did you set the image Origin to? As I said, it should be set to zero.
> So by calling volume->SetOrigin(-origin), my guess is that you were
> essentially cancelling the origin you had set in the image.
>
> When you say GetMatrix() returns an identity matrix, do you mean that you
> get an identity matrix after doing the following:
>
> [ set the elements of a 4x4 matrix ]
> volume->SetUserMatrix(matrix);
> newmatrix = volume->GetMatrix();
>
> the "newmatrix" is identity, even though "matrix" was not? Because that is
> not what I see myself... I see a "newmatrix" that is the result of
> compositing my "matrix" with whatever other transformations have been set
> for the volume (i.e. via SetOrigin(), SetPosition(), SetOrientation()).
> I.e. the final transform matrix used by the volume is your usermatrix
> _plus_ the other transformations.
>
> In any case, if you have to call SetPosition(), SetOrigin() on the volume,
> then there must be something you have not accounted for in the matrix.
> Either that, or you have set the Origin of your vtkImageData when you
> shouldn't have.
>
> Also, please note that my discussion assumes that you retrieved the DICOM
> Spacing, Orientation, and Position directly from the DICOM header. I'm
> assuming that ITK is not manipulating these values in any way before
> presenting them to you.
>
> - David
>
>
> 2011/7/5 Dženan Zukić <dzenanz at gmail.com>
>
>> I built matrix myself and used SetUserMatrix. But that did not work. In
>> addition to that I had to use SetPosition and SetOrigin. My question was why
>> did I have to do that in addition to the manually constructed matrix which
>> also contained translation part.
>>
>> Without any manual trickery GetMatrix returns an identity matrix. The
>> spacing which is taken from vtkImage is not in that matrix, but somewhere
>> else. The other information (orientation and position) is completely
>> ignored.
>>
>> Regards,
>> Dženan
>>
>> 2011/7/5 David Gobbi <david.gobbi at gmail.com>
>>
>>> 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/f687e666/attachment.htm>
More information about the vtkusers
mailing list