[Insight-users] Convert image's direction cosines to transform matrix

[Dženan Zukić]

Hi Luis!

Thanks for the provided code, but in the end it didn't really help me (aside
from encouraging me that transformations are not that complicated). I
re-read 3D matrix operations, did a lot of trial and error, and ended up
with this code (it might be useful to someone):

osg::ref_ptr<osg::Group> stuff=segmentInterestingStructures(...); //vertices
have physical image coordinates
osg::ref_ptr<osg::Group> isosurface=marchingCubesGPU(...); //accounts for
spacing, so scaling is handled
VisualizingImageType::DirectionType d=visualizing->GetDirection();
d=d.GetTranspose(); //ITK matrix is in column-major format, OSG matrix is
row major (trial and error)

osg::Matrix matrix(d(0,0), d(0,1), d(0,2), 0,
   d(1,0), d(1,1), d(1,2), 0,
   d(2,0), d(2,1), d(2,2), 0,
   0,   0,   0,   1);
osg::Matrix t;
VisualizingImageType::PointType origin=visualizing->GetOrigin();
t.makeTranslate(origin[0], origin[1], origin[2]);
matrix.postMult(t);
osg::ref_ptr<osg::MatrixTransform> mt=new osg::MatrixTransform(matrix);

osg::ref_ptr<Group> grp=new osg::Group;
mt->addChild(isosurface);
grp->addChild(mt);
grp->addChild(stuff);
mainForm.vis->setSceneData(grp);

Regards,
Dženan

2010/4/28 Luis Ibanez <luis.ibanez at kitware.com>

>
> Hi Dženan
>
> This is indeed a common problem,
> so, it is great that you bring it up.
>
>
> Following the numbering of your questions:
>
> a) No, there is no direct way of forcing the direction
>     matrix to be an identity. (and if there was such
>     method we will strongly discourage you from using it)
>
> b)  You can change the direction matrix of the image,
>      by using the ChangeImageInformationFilter.
>      (but we strongly discourage you from doing so).
>
> c)  The RIGHT thing to do is to use the direction
>      matrix to correct the orientation of the isosurface
>      after you extract it.  In that way, the iso-surface
>      will occupy the same coordinate system as the
>      original image.
>
>
> Options (a) and (b) are VERY dangerous if you are
> dealing with real patient data.
>
>
> Please find attached an example of how to do this
> with ITK and VTK. You could directly use this example,
> or you can probably follow the math and apply the
> equivalent processing using OSG.
>
>
> You will find the full source code at:
>
>                InsightApplications/Auxiliary/vtk
>
>
>
>
>       Regards,
>
>
>              Luis
>
>
> ---------------------------------------------------------------------
> 2010/4/27 Dženan Zukić <dzenanz at gmail.com>
>
>> Hi everybody,
>>
>> I guess someone has dealt with a problem similar to mine before. I am
>> doing segmentation in a magnetic resonance image. The result of my
>> segmentation is a polygonal model. Problem: How to visualize the resulting
>> polygonal model along with scalar-field image?
>>
>> I use Qt, ITK and OpenSceneGraph. First I tried using osgVolume
>> (ray-casting renderer) together with polygonal data in the same scene, but
>> with that setup the depth information is lost. I resorted to using
>> isosurface of the image instead of the image itself. It works, but I have to
>> manage physical coordinates by myself - which is both a reinvention of the
>> wheel and error-prone.
>>
>> So I decided to use ITK's physical coordinate managing routines
>> (TransformPhysicalPointToIndex), but I ran into a problem. GetDirection does
>> not return identity matrix for my images (at least not for sagittal ones),
>> therefore my segmented parts end up displayed *outside* of the image.
>> a) Is there a way to force image reading routines to load images into
>> memory in such a way to have an identity matrix for directions?
>> b) If not, is there a way to transform image data (using a filter) to get
>> an image with identity matrix as "directions". Is that
>> filter ChangeInformationImageFilter?
>> c) What is the proper way to construct a transformation matrix to be
>> applied to the isosurface of the image (which was generated assuming
>> identity matrix directions) in order to bring it to same coordinates I get
>> when using TransformIndexToPhysicalPoint? Just using matrix returned
>> by GetDirection and adding translation to account for origin does not work.
>>
>> If anyone knows a better way to visualize a polygonal model alongside 3D
>> image, I would like to hear it.
>> Off topic question: is combining polygons with scalar fields so simple in
>> VTK to justify a switch from OSG to VTK?
>>
>> Thanks to anyone who takes interest in this.
>>
>> Best regards,
>> Dženan
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20100428/8c9ea160/attachment.htm>