[vtk-developers] vtkNIFTIImageReader coordinates
David Gobbi
david.gobbi at gmail.com
Sun Nov 13 16:26:21 EST 2016
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.
- David
On Sun, Nov 13, 2016 at 9:50 AM, padraig <padraig.looney at gmail.com> wrote:
> Hi David,
>
>
> I am still confused by this. Volumes that ITK tells me have the same
> bounds, VTK displays in Paraview as not having the same bounds.
>
> Consider this header from
>
> vtkNIFTIImageHeader (0x221e6b0)
> Debug: Off
> Modified Time: 85
> Reference Count: 1
> Registered Events: (none)
> DimInfo: 0x0
> Dim: 3 224 153 52 1 1 1 1
> PixDim: 1 0.9195 1.40018 2.89048 0 0 0 0
> VoxOffset:352
> IntentP1: 0
> IntentP2: 0
> IntentP3: 0
> IntentCode: 0
> DataType: 2
> BitPix: 8
> SliceStart: 0
> SclSlope: 1
> SclInter: 0
> SliceEnd: 0
> SliceCode: 0
> XYZTUnits: 0x2
> CalMax: 0
> CalMin: 0
> SliceDuration: 0
> TOffset: 0
> Descrip: ""
> AuxFile: ""
> QFormCode: 2
> SFormCode: 1
> QuaternB: 0
> QuaternC: 0
> QuaternD: 1
> QOffsetX: 102.536
> QOffsetY: 106.413
> QOffsetZ: 30.1491
> SRowX: -0.9195 0 0 102.536
> SRowY: 0 -1.40018 0 106.413
> SRowZ: 0 0 2.89048 30.1491
> IntentName: ""
> Magic: "n+1"
>
> That is the same as the affine from nibabel
>
>
> array([[ -0.91949999, 0. , 0. , 102.53600311],
> [ 0. , -1.40017998, 0. , 106.41300201],
> [ 0. , 0. , 2.89048004, 30.14909935],
> [ 0. , 0. , 0. , 1. ]])
>
>
> If you look at the SForm matrix below shouldn't the bounds be
>
> X: (224-1)*-0.9195 + 102.536 = -102.51249999999999
>
> Y: (153-1)*-1.40018 + 106.413 =--106.41436
> Z: (52-1)*2.89048 + 30.1491 = 177.56358
>
>
> In ITK I have checked the physical coordinates of the image and found the
> bounds to be, as expected from the above,
>
> [-102.536, -106.413, 30.1491][102.512, 106.414, 177.564]
>
> In VTK I have checked the centre and origin with
>
> vtkNew<vtkNIFTIImageReader> niiReader;
> niiReader->SetFileName(str);
> niiReader->Update();
> int size[3];
> double center[3], spacing[3], bounds[6];
> niiReader->GetOutput()->GetDimensions(size);
> niiReader->GetOutput()->GetCenter(center);
> niiReader->GetOutput()->GetBounds(bounds);
> std::vector<double> bounds_vec(&bounds[0],&bounds[6]);
> std::cout << bounds[0] << " " << bounds[1] << " " << bounds[2] << " "
> << bounds[3] << " " << bounds[4] << " " << bounds[5] << std::endl;
> niiReader->GetOutput()->GetSpacing(spacing);
> double center1[3] = { center[0], center[1], center[2] };
> double center2[3] = { center[0], center[1], center[2] };
> std::cout << center[0] << " " << center[1] << " " << center[2] <<
> std::endl;
> std::cout << size[0] << " " << size[1] << " " << size[2] << std::endl;
>
> and get
>
> 0 205.048 0 212.827 0 147.414
> 102.524 106.414 73.7072
> 224 153 52
>
>
> On 12/11/16 15:39, David Gobbi wrote:
>
> Hi Padraig,
>
> The "direction" or "orientation" is the 3x3 portion of the 4x4 matrix,
> after dividing out the scale.
>
> So if ITK is giving a direction like the following, it means that ITK
> allows the direction to include a mirror-flip (in this case, it's a flip
> along the z-axis):
>
> 1 0 0
> 0 1 0
> 0 0 -1
>
> The vtkNIFTIImageReader doesn't allow a mirror flip in the QFormMatrix, so
> instead it re-orders the nifti slices in memory to undo the flip. When it
> does this, it also has to adjust the 4th column of the QFormMatrix since
> the 4th column now indicates the position of the last slice in the file,
> instead of the first slice in the file. The difference in position between
> the first slice and the last slice is (n - 1)*sz, where n is the number of
> slices and sz is the slice spacing. For your image, we can add (n - 1)*sz
> to -117.265 in order to get 30.149:
>
> -117.265 + (52 - 1)*2.89048 = 30.149
>
> This accounts for the discrepancy.
>
> - David
>
> On Sat, Nov 12, 2016 at 7:42 AM, padraig <padraig.looney at gmail.com> wrote:
>
>> Hi David,
>>
>>
>> vtkNIFTIImageHeader (0x3346680)
>> Debug: Off
>> Modified Time: 85
>> Reference Count: 1
>> Registered Events: (none)
>> DimInfo: 0x0
>> Dim: 3 224 153 52 1 1 1 1
>> PixDim: -1 0.9195 1.40018 2.89048 0 0 0 0
>> VoxOffset:352
>> IntentP1: 0
>> IntentP2: 0
>> IntentP3: 0
>> IntentCode: 0
>> DataType: 2
>> BitPix: 8
>> SliceStart: 0
>> SclSlope: 1
>> SclInter: 0
>> SliceEnd: 0
>> SliceCode: 0
>> XYZTUnits: 0x2
>> CalMax: 0
>> CalMin: 0
>> SliceDuration: 0
>> TOffset: 0
>> Descrip: ""
>> AuxFile: ""
>> QFormCode: 2
>> SFormCode: 1
>> QuaternB: 0
>> QuaternC: 0
>> QuaternD: 1
>> QOffsetX: 102.536
>> QOffsetY: 106.413
>> QOffsetZ: 30.1491
>> SRowX: -0.9195 0 -0 102.536
>> SRowY: 0 -1.40018 -0 106.413
>> SRowZ: 0 0 -2.89048 30.1491
>> IntentName: ""
>> Magic: "n+1"
>>
>> The direction in ITK is
>>
>>
>> 1 0 0
>> 0 1 0
>> 0 0 -1
>> I don't know how to see the direction in nibabel
>>
>>
>> On 12/11/16 14:10, David Gobbi wrote:
>>
>> Hi Padraig,
>>
>> Can you show us the raw qform and sform fields from the NIFTI header
>> itself? In VTK, you can do this as follows:
>>
>> >>> print(reader.GetNIFTIHeader())
>>
>> The difference you see between nibabel, ITK, and VTK has to do with
>> differences in the way these packages deal with the pixdim and the qfac of
>> the nifti file. If you print the Direction from the ITK image you might
>> find that it differs from nibabel, too.
>>
>> For VTK, you can read all the gory details here:
>>
>> http://gitlab.kitware.com/vtk/vtk/blob/v7.0.0/IO/Image/vtkNI
>> FTIImageReader.cxx#L738
>>
>> - David
>>
>>
>> On Sat, Nov 12, 2016 at 6:10 AM, padraig <padraig.looney at gmail.com>
>> wrote:
>>
>>> I have noticed some inconsistencies in how NIFTI images are handled in
>>> VTK, ITK, and NiBabel. ITK and NiBabel are in agreement with each other and
>>> differ with VTK on the position of the volume. The VTK output for a file is
>>> below.
>>>
>>> vtkNIFTIImageReader (0x1de67d0)
>>> Debug: Off
>>> Modified Time: 91
>>> Reference Count: 2
>>> Registered Events: (none)
>>> Executive: 0x1de7720
>>> ErrorCode: Success
>>> Information: 0x1de6ad0
>>> AbortExecute: Off
>>> Progress: 1
>>> Progress Text: (None)
>>> FileName: /home/padraig/Desktop/59674/Seeding/I0000003_59674_sn.nii
>>> FileNames: 0
>>> FilePrefix: (none)
>>> FilePattern: %s.%d
>>> FileNameSliceOffset: 0
>>> FileNameSliceSpacing: 1
>>> DataScalarType: unsigned char
>>> NumberOfScalarComponents: 1
>>> File Dimensionality: 3
>>> File Lower Left: On
>>> Swap Bytes: Off
>>> DataIncrements: (1, 1)
>>> DataExtent: (0, 223, 0, 152, 0, 51)
>>> DataSpacing: (0.9195, 1.40018, 2.89048)
>>> DataOrigin: (0, 0, 0)
>>> HeaderSize: 352
>>> Internal File Name: (none)
>>> TimeAsVector: Off
>>> TimeDimension: 1
>>> TimeSpacing: 1
>>> RescaleSlope: 1
>>> RescaleIntercept: 0
>>> QFac: -1
>>> QFormMatrix: -1 0 0 102.536 0 -1 0 106.413 0 0 1 -117.265 0 0 0 1
>>> SFormMatrix: -1 0 -0 102.536 0 -1 -0 106.413 0 0 1 -117.265 0 0 0 1
>>> NIFTIHeader:
>>> PlanarRGB: Off
>>>
>>> the qform given by nibabel is
>>>
>>> >>> ni_ob.get_qform()
>>> array([[ -0.91949999, 0. , 0. , 102.53600311],
>>> [ 0. , -1.40017998, 0. , 106.41300201],
>>> [ 0. , 0. , 2.89048004, 30.14909935],
>>> [ 0. , 0. , 0. , 1. ]])
>>>
>>>
>>> ITK has the origin at
>>>
>>> [-102.536, -106.413, 30.1491]
>>>
>>
>>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/vtk-developers/attachments/20161113/6201df4e/attachment-0001.html>
More information about the vtk-developers
mailing list