[Paraview] Volume positioning

padraig padraig.looney at gmail.com
Mon Nov 14 17:50:46 EST 2016


Hi Cory,


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. I thought initially the issue was with VTK but on 
emailing the devel group I got the following response


> 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
>>>>>>>
>>>>>>>
>>>
>
>



More information about the ParaView mailing list