[Paraview] Volume positioning

Cory Quammen cory.quammen at kitware.com
Mon Nov 14 20:02:36 EST 2016 

> It doesn't seem right to me that the coordinates are not correctly applied
> in Paraview and that a user should have to apply a transformation.

Yes, it is unfortunate. We wish ParaView could support every feature
one may find useful. Contributions from the community to support new
features are welcome!

- Cory

>
>
>> Hi Padraig,
>>
>> What you are seeing is one of the big (very big) differences between ITK
>> images and VTK images.  It comes about because in ITK, the direction cosines
>> are stored in the image class, but in VTK, they are not, and they must
>> instead be stored in a separate matrix.  You're probably wondering why this
>> matters at all, since your image has direction cosines of
>> ((1,0,0),(0,1,0),(0,0,1)).  Hopefully I'll be able to provide an explanation
>> that isn't too muddy.
>>
>> In ITK, the transformation from the (I,J,K) indices of each voxel to the
>> (x,y,z) coordinates is as follows:
>>
>>   let q = (I,J,K)
>>   let p = (x,y,z)
>>
>>   p = R*S*q + o
>>
>> where R is ITK's Direction matrix, S is the "scale" matrix (a diagonal
>> matrix made from the Spacing), and "o" is ITK's Origin.
>>
>> Since vtkImageData does not have a Direction matrix, but most medical
>> image formats do have a Direction, it is necessary for VTK medical image
>> readers to produce a vtkMatrix4x4 in addition to the vtkImageData.  In VTK,
>> the equivalent equation is:
>>
>>   p = R*(S*q + o)
>>
>> Or we can expand this,
>>
>>   p = R*S*q + R*o
>>
>> See how the ITK Origin and the VTK Origin mean different things?  The ITK
>> origin is rotated with respect to the VTK origin.  In fact, the way that VTK
>> defines the image "Origin" is incompatible with NIFTI (and with DICOM, too).
>>
>> So, the vtkNIFTIImageReader always sets the VTK image Origin to (0,0,0).
>> Then, it stores an ITK-style Origin as the 4th column of the SFormMatrix and
>> the QFormMatrix.  That's an advantage of using a 4x4 matrix: it can store
>> both a rotation and an offset.  When we introduce this matrix "M" and set
>> the VTK image Origin to (0,0,0),
>>
>>   p = M*S*q
>>
>> Of course, now we have to make p and q into (I, J, K, 1) and (x, y, z, 1)
>> in order to use the 4x4 matrix.  I'll let you work out the rest of the
>> details.
>>
>> In order to properly compute the bounds of an oriented image in VTK, you
>> have to compute the coordinates of the four corners and then multiply them
>> by the matrix M (i.e. by the SFormMatrix).  In your case, since the
>> "orientation" is identity, you could simply add the offset.
>>
>> If you're wondering why the reader doesn't set the VTK image Origin to be
>> equal to the offset for axial images, it's to avoid inconsistency between
>> the way the reader treats axial images as compared to sagittal or coronal.
>>
>> For Paraview, you'll have to find a way to get Paraview to use the
>> SFormMatrix (or the QFormMatrix) that is provided by the reader.
>> Unfortunately, I don't know enough about Paraview to help with that...
>> hopefully one of the Paraview users around here can help.
>>
> Thanks
>
> Pádraig
>
>
>
> On 14/11/16 22:34, Cory Quammen wrote:
>>
>> On Mon, Nov 14, 2016 at 5:25 PM, padraig <padraig.looney at gmail.com> wrote:
>>>
>>> I have been using the analyseniftiio plugin to open .nii files in
>>> paraview.
>>
>> Ah, I see. I had forgotten about that one.
>>
>>> In Slicer 4.6.2 they do not occupy the same physical space. The direction
>>> is
>>> different although the origin is the same.
>>
>> I looked again in Slicer and you are right. Looking at
>> vtkNIfTIReader.cxx, the direction information seems to be ignored by
>> ParaView. Which makes sense because VTK and ParaView do not support
>> the concept of direction cosines, just origin and spacing. In
>> ParaView, you'll have to transform one of the volumes using the
>> Transform filter and scaling the Z direction by -1. Because the x and
>> y direction diagonals are also -1, you may want to scale those axes by
>> -1 as well.
>>
>> HTH,
>> Cory
>>
>>> Thanks
>>>
>>>
>>>
>>> On 14/11/16 22:11, Cory Quammen wrote:
>>>>
>>>> Are you somehow opening these files in ParaView? I was not able to
>>>> either gzipped or not, which is what I expect because ParaView does
>>>> not have a Nifti reader to my knowledge.
>>>>
>>>> In Slicer, which can open these files, the volumes are reported to
>>>> occupy the same space with an origin at (102.536, 106.413, 30.149).
>>>>
>>>> Please provide more details, step-by-step, about how you are viewing
>>>> these files in ParaView.
>>>>
>>>> Thanks,
>>>> Cory
>>>>
>>>> On Mon, Nov 14, 2016 at 10:29 AM, padraig <padraig.looney at gmail.com>
>>>> wrote:
>>>>>
>>>>> Attached are two volumes that have an overlap in paraview but the ITK
>>>>> volume
>>>>> I find using the code below means there should be no overlap
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> On 14/11/16 15:09, Cory Quammen wrote:
>>>>>>
>>>>>> MHA and NIFTII definitely contain position information that ParaView
>>>>>> should read. Do you have any small-ish representative volumes you can
>>>>>> share (privately with me if needed).
>>>>>>
>>>>>> On Mon, Nov 14, 2016 at 10:04 AM, padraig <padraig.looney at gmail.com>
>>>>>> wrote:
>>>>>>>
>>>>>>> I have used MHA and NIFTII. I have converted the MHA into NIFTII
>>>>>>> using
>>>>>>> both
>>>>>>> c3d and ITK.
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> On 14/11/16 14:59, Cory Quammen wrote:
>>>>>>>>
>>>>>>>> What file format are you using to load the volumes into ParaView? A
>>>>>>>> number of formats support volume positioning, so this should be
>>>>>>>> possible, unless you are loading a series of TIFF images, for
>>>>>>>> example.
>>>>>>>>
>>>>>>>> Thanks,
>>>>>>>> Cory
>>>>>>>>
>>>>>>>> On Mon, Nov 14, 2016 at 5:32 AM, padraig <padraig.looney at gmail.com>
>>>>>>>> wrote:
>>>>>>>>>
>>>>>>>>> Dear list,
>>>>>>>>>
>>>>>>>>> I have been having problems with the positioning of volumes using
>>>>>>>>> Paraview.
>>>>>>>>> ITK tells me that, using,
>>>>>>>>>
>>>>>>>>>         IteratorType  it2( img_input,
>>>>>>>>> img_input->GetLargestPossibleRegion()
>>>>>>>>> );
>>>>>>>>>
>>>>>>>>>         it2.GoToBegin();
>>>>>>>>>         ImageType::IndexType begin = it2.GetIndex();
>>>>>>>>>
>>>>>>>>> img_input->TransformIndexToPhysicalPoint(it2.GetIndex(),p0);
>>>>>>>>>         it2.GoToEnd();
>>>>>>>>>         --it2;
>>>>>>>>>
>>>>>>>>> img_input->TransformIndexToPhysicalPoint(it2.GetIndex(),p1);
>>>>>>>>>         std::cout << p0 <<  p1 << std::endl;
>>>>>>>>>
>>>>>>>>>      two volumes I have have the positions
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> [-102.536, -106.413, 30.1491][102.512, 106.414, 177.564]
>>>>>>>>>
>>>>>>>>> and
>>>>>>>>>
>>>>>>>>> [-102.536, -106.413, 30.1491][102.512, 106.414, -117.265]
>>>>>>>>>
>>>>>>>>> When I load these into Paraview they occupy the same volume. In
>>>>>>>>> Slicer3D
>>>>>>>>> the
>>>>>>>>> volumes are distinct as I expect above.
>>>>>>>>>
>>>>>>>>> Thanks
>>>>>>>>> Pádraig
>>>>>>>>> _______________________________________________
>>>>>>>>> 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 ParaView Wiki at:
>>>>>>>>> http://paraview.org/Wiki/ParaView
>>>>>>>>>
>>>>>>>>> Search the list archives at: http://markmail.org/search/?q=ParaView
>>>>>>>>>
>>>>>>>>> Follow this link to subscribe/unsubscribe:
>>>>>>>>> http://public.kitware.com/mailman/listinfo/paraview
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>
>>
>>
>



-- 
Cory Quammen
Staff R&D Engineer
Kitware, Inc.

More information about the ParaView mailing list