[Neuroimaging] DICOM Orientation and World<-->Image coordinate transformation

[Matthew Brett]

Hi,

On Tue, Sep 6, 2016 at 6:34 AM, Steve Pieper <pieper at isomics.com> wrote:
> Hi Athanasios -
>
> To get the scan direction you'll need to look at the relative
> ImagePositionPatient points from slice to slice.  Note that the scan
> direction is not always the cross product of the row and column orientations
> since the scan may go in the other direction from a right handed cross
> product or the slices can be sheared (or even at arbitrary locations)..
> There are lots of other things that can happen too, like irregular spacing,
> missing slices, etc, but usually just normalizing the vector between your
> origin and any slice in the scan will be what you want.
>
> This code will give you an idea:
>
> https://github.com/Slicer/Slicer/blob/master/Modules/Scripted/DICOMPlugins/DICOMScalarVolumePlugin.py#L195-L216
>

>From your code, you are missing a valid third column for your affine.
I believe that column will be all zeros from your code. This is what
the later part of the DICOM orientation page is talking about, and
what Steve is referring to as the "slice direction".

Steve is quite right that the slice direction need not be the
cross-product of the first two, and the DICOM information can often
tell you what that slice direction vector is, but assuming for a
moment that it is the cross product, and that you are looking at the
first slice of the volume, then you'd want something like:

"""
import numpy as np

ImageOrientationPatient = [0.999857, 0.00390641, 0.0164496,
                           -0.00741602, 0.975738, 0.218818]

ImagePositionPatient = [-127.773, -105.599, -94.5758]

PixelSpacing = [0.4688, 0.4688]

slice_spacing = 3.0  # ?

# Make F array from DICOM orientation page
F = np.fliplr(np.reshape(ImageOrientationPatient, (2, 3)).T)
rotations = np.eye(3)
rotations[:, :2] = F
# Third direction cosine from cross-product of first two
rotations[:, 2] = np.cross(F[:, 0], F[:, 1])
# Add the zooms
zooms = np.diag(PixelSpacing + [slice_spacing])

# Make the affine
affine = np.diag([0., 0, 0, 1])
affine[:3, :3] = rotations.dot(zooms)
affine[:3, 3] = ImagePositionPatient

np.set_printoptions(precision=4, suppress=True)
print(affine)
"""

But - Steve's suggestion is more general - this code is just to give
you an idea.

Best,

Matthew