[vtkusers] Announcement: New image rendering classes for VTK
Jothy
jothybasu at gmail.com
Wed Jul 6 07:23:07 EDT 2011
Could someone please contribute an example to VTK wiki on this issue. So it
will be easier to use and understand these methods.
Thank you
Jothy
2011/7/6 Dženan Zukić <dzenanz at gmail.com>
> 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
>>>>>> >>>> >
>>>>>> >>>> >
>>>>>> >>>
>>>>>> >>
>>>>>> >
>>>>>> >
>>>>>>
>>>>>>
>>>>>
>>>>
>>>
>>
>
> _______________________________________________
> 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
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.vtk.org/pipermail/vtkusers/attachments/20110706/e3e3c6de/attachment.htm>
More information about the vtkusers
mailing list