<font size="2"><font face="verdana,sans-serif">I just checked: vtkImageData origin is the same as itk::Image's origin.</font></font><div><font size="2"><font face="verdana,sans-serif"><br></font></font></div><div><font size="2"><font face="verdana,sans-serif">But your earlier mail explained:<br>
</font></font></div><div><span style="border-collapse:collapse;font-family:arial, sans-serif;font-size:13px">VTK -> scale by Spacing, translate by Origin, then rotate and translate by a 4x4 matrix</span></div>
<div><font face="arial, sans-serif"><span style="border-collapse:collapse"><br></span></font></div><div><font face="arial, sans-serif"><span style="border-collapse:collapse">So SetPosition was needed </span></font><span style="border-collapse:collapse;font-family:arial, sans-serif">(I realized that SetOrigin is not needed</span><span style="border-collapse:collapse;font-family:arial, sans-serif">) to counteract the built-in translation by origin, and then the 4x4 matrix does its job.</span></div>
<div><font face="arial, sans-serif"><span style="border-collapse:collapse"><br></span></font></div><div><font face="arial, sans-serif"><span style="border-collapse:collapse">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):</span></font></div>
<div><font face="arial, sans-serif"><span style="border-collapse:collapse"><br></span></font></div><div><font face="arial, sans-serif"><span style="border-collapse:collapse">Mitk=Mscale*Mrot*Mtrans</span></font></div>
<div><font face="arial, sans-serif"><span style="border-collapse:collapse">Mvtk=Mscale*Mtrans*M4x4</span></font></div><div><font face="arial, sans-serif"><span style="border-collapse:collapse"><br>
</span></font></div><div><font face="arial, sans-serif"><span style="border-collapse:collapse">Solving the matrix equation Mitk=Mvtk, we get M4x4=(Mtrans)^-1*Mrot*Mtrans.</span></font></div>
<div><font face="arial, sans-serif"><span style="border-collapse:collapse"><br></span></font></div><div><font face="arial, sans-serif"><span style="border-collapse:collapse">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.</span></font></div>
<div><font face="arial, sans-serif"><span style="border-collapse:collapse"><br></span></font></div><div><font face="arial, sans-serif"><span style="border-collapse:collapse">And after a note on notation: </span></font><a href="http://www.gamedev.net/topic/321862-pre-or-post-multiplication-for-matrices/" target="_blank">http://www.gamedev.net/topic/321862-pre-or-post-multiplication-for-matrices/</a>, things are clear now. So finally we can have this (no usage of SetPosition):</div>
<div><br></div><div><div><div>VisualizingImageType::DirectionType d=logic->visualizing->GetDirection();</div><div>VisualizingImageType::PointType origin=logic->visualizing->GetOrigin();</div><div>vtkMatrix4x4 *mat=vtkMatrix4x4::New(); //start with identity matrix</div>
<div>for (int i=0; i<3; i++)</div><div> for (int k=0; k<3; k++)</div><div> mat->SetElement(i,k, d(i,k));</div><div>for (int i=0; i<3; i++)</div><div> mat->SetElement(i,3, origin[i]);</div><div>
vtkMatrix4x4 *t=vtkMatrix4x4::New();</div><div>for (int i=0; i<3; i++)</div><div> t->SetElement(i,3, -origin[i]); //invert(Mtrans)</div><div>vtkMatrix4x4::Multiply4x4(mat,t,mat); //keep it all in UserMatrix</div>
<div>volume->SetUserMatrix(mat);</div></div></div><div><br></div><div>Dženan</div>
<div><font face="arial, sans-serif"><span style="border-collapse:collapse"><br></span></font><div class="gmail_quote">2011/7/5 David Gobbi <span dir="ltr"><<a href="mailto:david.gobbi@gmail.com" target="_blank">david.gobbi@gmail.com</a>></span><br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div>Of course the matrix is identity before you call SetOrigin()... you don't call SetUserMatrix() until a few lines after that.</div>
<div><br></div><div>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:</div>
<div><br></div><div>ITK: scale by itk::Image spacing, rotate by itk::Image direction, translate by itk::Image origin</div><div><br></div><div>In VTK, image display works like this:</div><div><br></div><div>VTK: scale by vtkImageData spacing, translate by vtkImageData origin, rotate (3x3 part of 4x4 matrix), translate (last column of 4x4 matrix)</div>
<div><br></div><div>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).</div>
<div><br></div><div>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.</div>
<div><div></div><div>
<div><br></div><div> - David</div><div><br></div><div><br><br><div class="gmail_quote">2011/7/5 Dženan Zukić <span dir="ltr"><<a href="mailto:dzenanz@gmail.com" target="_blank">dzenanz@gmail.com</a>></span><br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<font size="2"><font face="verdana,sans-serif">Take a look at the example: </font></font><font size="2"><font face="verdana,sans-serif"><a href="http://vtk.org/Wiki/VTK/Examples/Cxx/VolumeRendering/itkVtkImageConvert" target="_blank">http://vtk.org/Wiki/VTK/Examples/Cxx/VolumeRendering/itkVtkImageConvert</a><br>
</font></font><div><font face="verdana, sans-serif"><br></font></div>If I call GetMatrix before the line "volume->SetOrigin(...)", it returns identity. If I call it after SetPosition, it contains the translation part.<div>
<br></div><div>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.</div>
<div><br></div><div>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.<div>
<div></div><div><br>
<div><font face="verdana, sans-serif"><br></font><div class="gmail_quote">2011/7/5 David Gobbi <span dir="ltr"><<a href="mailto:david.gobbi@gmail.com" target="_blank">david.gobbi@gmail.com</a>></span><br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
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.<div><br></div>
<div>When you say GetMatrix() returns an identity matrix, do you mean that you get an identity matrix after doing the following:</div><div><br></div><div>[ set the elements of a 4x4 matrix ]</div><div>volume->SetUserMatrix(matrix);</div>
<div>newmatrix = volume->GetMatrix();</div><div><br></div><div>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.</div>
<div><br></div><div>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.</div>
<div><br></div><div>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.</div>
<div><div></div><div>
<div><br></div><div> - David<br><div><br></div><div><br><div class="gmail_quote">2011/7/5 Dženan Zukić <span dir="ltr"><<a href="mailto:dzenanz@gmail.com" target="_blank">dzenanz@gmail.com</a>></span><br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<font size="2"><font face="verdana,sans-serif">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.<br>
</font></font><div><font face="verdana, sans-serif"><br></font></div><div><font face="verdana, sans-serif">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.</font></div>
<div><font face="verdana, sans-serif"><br></font></div><div><font face="verdana, sans-serif">Regards,</font></div><div><font face="verdana, sans-serif">Dženan</font></div><div><div></div><div>
<div><font face="verdana, sans-serif"><br></font><div class="gmail_quote">2011/7/5 David Gobbi <span dir="ltr"><<a href="mailto:david.gobbi@gmail.com" target="_blank">david.gobbi@gmail.com</a>></span><br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
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.<br><br>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.<br>
<br>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".<br>
<br>In summary:<br>DICOM -> scale by Spacing, rotate by Orientation, translate by Position<br>VTK -> scale by Spacing, translate by Origin, then rotate and translate by a 4x4 matrix<br><br>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.<br>
<br>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.<br>
<br>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.<div>
<br></div><div>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.</div>
<div><br></div><div> - David<br><br><br><br><br>2011/7/5 Dženan Zukić <<a href="mailto:dzenanz@gmail.com" target="_blank">dzenanz@gmail.com</a>>:<div><div></div><div><br>> I solved this (see my other mail). But I still don't know why I had to set<br>
> origin and position:<br>> volume->SetOrigin(-origin[0], -origin[1], -origin[2]);<br>> volume->SetPosition(-origin[0], -origin[1], -origin[2]);<br>> Can you explain it? Also, a link to some text explaining vtk's coordinate<br>
> transformation system for future reference would be very appreciated (if<br>> such a text exists). I had been searching but could not find a source which<br>> explains in which order are the transformation matrices applied, in which<br>
> space are voxels of a vtkVolume etc.<br>> Regards,<br>> Dženan<br>> 2011/7/1 Dženan Zukić <<a href="mailto:dzenanz@gmail.com" target="_blank">dzenanz@gmail.com</a>><br>>><br>>> I did the math, and it should be as simple as this:<br>
>><br>>> VisualizingImageType::DirectionType d=logic->visualizing->GetDirection();<br>>> vtkMatrix4x4 *mat=vtkMatrix4x4::New(); //identity matrix<br>>> for (int i=0; i<3; i++)<br>>> for (int k=0; k<3; k++)<br>
>> mat->SetElement(i,k, d(i,k));<br>>> VisualizingImageType::PointType origin=logic->visualizing->GetOrigin();<br>>> for (int i=0; i<3; i++)<br>>> mat->SetElement(i,3, origin[i]);<br>
>> //mat->Invert(); //inverse produces even wrong orientation<br>>> volume->SetUserMatrix(mat);<br>>> However this does not work. This is a simple transformation from scaled<br>>> index space (vtkVolume handles spacing) to DICOM physical space. If no one<br>
>> has any suggestion, I will simply have to create a system of linear<br>>> equations from which to construct the VTK's transform matrix. And equations<br>>> I will have to obtain by passing a number of points through ITK's<br>
>> transformation and add some simple objects to the vtk scene at that spot and<br>>> on index space spot but with VTK transformation applied.<br>>> Regards,<br>>> Dženan<br>>> 2011/6/29 Dženan Zukić <<a href="mailto:dzenanz@gmail.com" target="_blank">dzenanz@gmail.com</a>><br>
>>><br>>>> I tried trial and error: premultiply translation, postmultiply, divide by<br>>>> and multiply by spacing, all that unsuccessfully. Now I need to work out the<br>>>> math :D<br>
>>><br>>>> 2011/6/29 David Gobbi <<a href="mailto:david.gobbi@gmail.com" target="_blank">david.gobbi@gmail.com</a>><br>>>>><br>>>>> Hi Dzenan,<br>>>>><br>>>>> For your code, my guess is that the "translation" part of the 4x4<br>
>>>> matrix should be set to something like the origin multiplied by the<br>>>>> rotation matrix, but there might be other little details that you need<br>>>>> to take care of. My only suggestion is that you don't try to fix it<br>
>>>> through trial-end-error. Use pencil-and-paper to work through the<br>>>>> math.<br>>>>><br>>>>> - David<br>>>>><br>>>>> 2011/6/29 Dženan Zukić <<a href="mailto:dzenanz@gmail.com" target="_blank">dzenanz@gmail.com</a>>:<br>
>>>> > Hi David,<br>>>>> ><br>>>>> > with this I try to do what you suggested, i.e. consider VTK's world<br>>>>> > coordinate system as DICOM patient coordinate system:<br>
>>>> > //red is the transformation-related code<br>>>>> > void MainWindow::updateVisualization()<br>>>>> > {<br>>>>> > typedef itk::ImageToVTKImageFilter<VisualizingImageType><br>
>>>> > itkVtkConverter;<br>>>>> > itkVtkConverter::Pointer conv=itkVtkConverter::New();<br>>>>> > conv->SetInput(logic->visualizing);<br>>>>> > vtkGPUVolumeRayCastMapper *mapper =<br>
>>>> > vtkGPUVolumeRayCastMapper::New();<br>>>>> > mapper->SetInput(conv->GetOutput());<br>>>>> > if (volume)<br>>>>> > volume->Delete();<br>
>>>> > volume=vtkVolume::New();<br>>>>> > volume->SetProperty( myTransferFunction );<br>>>>> > volume->SetMapper( mapper );<br>>>>> > mapper->SetBlendModeToComposite();<br>
>>>> > VisualizingImageType::DirectionType<br>>>>> > d=logic->visualizing->GetDirection();<br>>>>> > vtkMatrix4x4 *mat=vtkMatrix4x4::New();<br>>>>> > for (int i=0; i<3; i++)<br>
>>>> > for (int k=0; k<3; k++)<br>>>>> > mat->SetElement(i,k, d(i,k));<br>>>>> > mat->SetElement(3,3, 1);<br>>>>> > VisualizingImageType::SpacingType sp =<br>
>>>> > logic->visualizing->GetSpacing();<br>>>>> > VisualizingImageType::PointType<br>>>>> > origin=logic->visualizing->GetOrigin();<br>>>>> > for (int i=0; i<3; i++)<br>
>>>> > mat->SetElement(i,3, origin[i]/sp[i]);<br>>>>> > volume->SetUserMatrix(mat); //orientation and size OK, position<br>>>>> > wrong<br>>>>> > vtkRenderer *renderer =<br>
>>>> > vis->GetRenderWindow()->GetRenderers()->GetFirstRenderer();<br>>>>> > renderer->RemoveAllViewProps();<br>>>>> > renderer->AddVolume( volume );<br>
>>>> > defaultCameraPos(); //centers view on the volume, looking at it<br>>>>> > from<br>>>>> > left<br>>>>> > vis->GetRenderWindow()->Render();<br>
>>>> ><br>>>>> ><br>>>>> > vis->GetRenderWindow()->GetRenderers()->GetFirstRenderer()->ResetCamera();<br>>>>> > }<br>>>>> > The polygonal data is in DICOM's patient coordinate system (vertex<br>
>>>> > positions<br>>>>> > from itk::Image->TransformIndexToPhysicalPoint() etc). The orientation<br>>>>> > is<br>>>>> > correct, but I can't get the position correctly. I have tried<br>
>>>> > -origin[i]/sp[i], origin[i]*sp[i], origin[i], -origin[i], -2*origin[i]<br>>>>> > and<br>>>>> > similar combinations but none of them worked. Any suggestion to what<br>>>>> > am I<br>
>>>> > doing wrong?<br>>>>> > Regards,<br>>>>> > Dženan<br>>>>> > 2011/6/29 David Gobbi <<a href="mailto:david.gobbi@gmail.com" target="_blank">david.gobbi@gmail.com</a>><br>
>>>> >><br>>>>> >> On Tue, Jun 28, 2011 at 9:40 AM, Dženan Zukić <<a href="mailto:dzenanz@gmail.com" target="_blank">dzenanz@gmail.com</a>><br>>>>> >> wrote:<br>
>>>> >> > Hi David,<br>
>>>> >> ><br>>>>> >> > I am interested in this too:<br>>>>> >> > Does the above approach affect the polygonal actors in the scene<br>>>>> >> > (segmented<br>
>>>> >> > parts of the image)?<br>>>>> >> > Regards,<br>>>>> >> > Dženan<br>>>>> >><br>>>>> >> When the DICOM patient coordinate system is used as the VTK world<br>
>>>> >> coordinate system, the assumption is that everything (images AND<br>>>>> >> polydata) would have a either already be in the patient coordinate<br>>>>> >> system, or if not, you would to have a 4x4 transform to set as the<br>
>>>> >> actor's UserMatrix (or UserTransform). I.e. a transform to bring the<br>>>>> >> data from its original coordinate system into the patient coordinate<br>>>>> >> system.<br>
>>>> >><br>>>>> >> - David<br>>>>> ><br>>>>> ><br>>>><br>>><br>><br>><br><br></div></div></div>
</blockquote></div><br></div>
</div></div></blockquote></div><br></div></div>
</div></div></blockquote></div><br></div>
</div></div></div>
</blockquote></div><br></div>
</div></div></blockquote></div><br></div>