[vtkusers] Announcement: New image rendering classes for VTK
Dženan Zukić
dzenanz at gmail.com
Wed Jul 6 07:04:43 EDT 2011
I just checked: vtkImageData origin is the same as itk::Image's origin.
But your earlier mail explained:
VTK -> scale by Spacing, translate by Origin, then rotate and translate by a
4x4 matrix
So SetPosition was needed (I realized that SetOrigin is not needed) to
counteract the built-in translation by origin, and then the 4x4 matrix does
its job.
If I knew earlier how VTK does it (exactly and concisely, like your
explanation above), I could have done the math (instead of stumbling onto
the solution by experimentation):
Mitk=Mscale*Mrot*Mtrans
Mvtk=Mscale*Mtrans*M4x4
Solving the matrix equation Mitk=Mvtk, we get M4x4=(Mtrans)^-1*Mrot*Mtrans.
Mrot*Mtrans is the 4x4 matrix I was already constructing (because scaling is
already handled by VTK). I didn't know that vtk also translates by the
origin before applying normal Prop3D transformations. And since
SetPosition(-origin) is the inverse of the translation matrix Mtrans, I was
in effect creating the above M4x4 matrix.
And after a note on notation:
http://www.gamedev.net/topic/321862-pre-or-post-multiplication-for-matrices/,
things are clear now. So finally we can have this (no usage of SetPosition):
VisualizingImageType::DirectionType d=logic->visualizing->GetDirection();
VisualizingImageType::PointType origin=logic->visualizing->GetOrigin();
vtkMatrix4x4 *mat=vtkMatrix4x4::New(); //start with identity matrix
for (int i=0; i<3; i++)
for (int k=0; k<3; k++)
mat->SetElement(i,k, d(i,k));
for (int i=0; i<3; i++)
mat->SetElement(i,3, origin[i]);
vtkMatrix4x4 *t=vtkMatrix4x4::New();
for (int i=0; i<3; i++)
t->SetElement(i,3, -origin[i]); //invert(Mtrans)
vtkMatrix4x4::Multiply4x4(mat,t,mat); //keep it all in UserMatrix
volume->SetUserMatrix(mat);
Dženan
2011/7/5 David Gobbi <david.gobbi at gmail.com>
> Of course the matrix is identity before you call SetOrigin()... you don't
> call SetUserMatrix() until a few lines after that.
>
> I suspect that the problem is that you are taking ITK's itkImage "Origin"
> and expecting that you can use it as VTK's vtkImageData "Origin". I don't
> think you can do this unless the Orientation is identity. I'm not 100% sure
> but I believe that ITK does image orientation like this:
>
> ITK: scale by itk::Image spacing, rotate by itk::Image direction, translate
> by itk::Image origin
>
> In VTK, image display works like this:
>
> VTK: scale by vtkImageData spacing, translate by vtkImageData origin,
> rotate (3x3 part of 4x4 matrix), translate (last column of 4x4 matrix)
>
> If the order of operations is different, then it is not possible to simply
> plug in the values like you do in your example. That is why I recommended
> setting the VTK image origin to zero (note that I said the vtkImageData
> origin, which is completely different from the vtkVolume origin even if they
> do happen to have the same name).
>
> If the itkImageToVTKImageFilter is setting the VTK image origin to the ITK
> image origin for images with a non-identity orientation then it is, to be
> blunt, screwing up. You will have to use vtkImageChangeInformation to fix
> the origin and move on from there.
>
> - David
>
>
>
> 2011/7/5 Dženan Zukić <dzenanz at gmail.com>
>
>> 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/20110706/441f88c0/attachment.htm>
More information about the vtkusers
mailing list