Images in physical space in Python

[matt.mccormick]

@blowekamp

It would be great if these libraries define the attributes the same way so they are interoperable.

Yes, lessons from the past emphasize the need for explicit, specific definitions.

To take another pass, and borrow from the ITK Software Guide:

• origin: N-dimensional sequence of floating-point values that specifies the location of the center of the 0-indexed pixel.
• spacing: N-dimensional sequence of floating-point values that specifies the physical distance between pixel centers in each direction.
• direction: NxN matrix of floating-point values that specifies the orientation of the image. Columns of the matrix are orthonormal direction cosines. For medical images, direction corresponds to an LPS (DICOM) coordinate system, i.e. the right-to-left, anterior-to-posterior, inferior-to-superior directions of the imaged person.

For a given index \vec{I} in 3D, the physical location \vec{P} is calculated as:

\begin{equation} \begin{pmatrix} P_1\\\ P_2\\\ P_3 \end{pmatrix} = \begin{pmatrix} O_1\\\ O_2\\\ O_3 \end{pmatrix} + \begin{pmatrix} D_{1,1} & D_{1,2} & D_{1,3}\\\ D_{2,1} & D_{2,2} & D_{2,3}\\\ D_{3,1} & D_{3,2} & D_{3,3} \end{pmatrix} * \begin{pmatrix} S_1 & 0 & 0 \\\ 0 & S_2 & 0 \\\ 0 & 0 & S_3 \end{pmatrix} * \begin{pmatrix} I_1\\\ I_2\\\ I_3 \end{pmatrix} \end{equation}

Where:

• \vec{I}: pixel space index.
• \vec{P}: resulting physical space position of the image index \vec{I}.
• \vec{O}, origin: physical space origin of the first image index.
• \mathcal{D}, direction: direction cosines matrix (orthonormal). It represents the orientation relationship between the image and the physical space coordinate system.
• \vec{S}, spacing: physical spacing between pixels of the same axis.

[matt.mccormick]

@ctrueden

Would it then make sense to create some validation code, which implementors can use to check their image objects for appropriate structure?

Good point. We can create a reference, pure NumPy-based implementation, create tests for the attributes and methods. We can also provide those tests in a way that any implementation that wants to be spatial-image-array-like duck can test if its interface quacks like a spatial-image-array-like duck.

@ctrueden

If so, where would it live? Should it be part of scipy? Who do we invite to the conversation? Or file an issue on GitHub, perhaps?

Perhaps we should create an independent package for now, e.g. spatialimage. This package could only depend on numpy, so it is lightweight dependency. We can move rapidly – after it matures, it could move to a different home.

[blowekamp (Bradley Lowekamp)]

@matt.mccormick That is looking like a good start.

There was a recent discussion of how to handle time. Similarly there are a couple different ways to order color. i.e. CXYZ, XYZTC etc. Simularly tensor images too. How can non-physical dimensions be handled? ITK uses a templated pixel type to separate a physical dimension form “pixel components”

There is also the definition and order of an “index”. We mean (i,j,k) where i is the fasted changing index, but then ndarrays are indexed via [k,j,i].

I expect others from different domains will have additional concerns and comments to improve the specification too. It may be worthwhile to go through a more formal proposal and comment process to aid in getting more eyes on the spec, and interest.

Would the LPS direction imply that 2-D it just LP with the origin being in the top right?

[lassoan (Andras Lasso)]

@blowekamp

It may be worth wild to go through a more formal proposal and comment process to aid in getting more eyes on the spec, and interest.

Fully agree.

@blowekamp

Would the LPS direction imply that 2-D it just LP with the origin being in the top right?

In most medical imaging applications, 2D images are still 3D images (just happen to have a single slice), with their position, spacing (pixel size and slice thickness), and axis directions defined in 3D space. In these images, axes would be still LPS.

If the 2D image is not defined in 3D space (microscopy, optical camera, etc.) then I don’t think LPS coordinate system would be applicable.

[matt.mccormick]

@blowekamp

There was a recent discussion of how to handle time. Similarly there are a couple different ways to order color. i.e. CXYZ, XYZTC etc.

@blowekamp

There is also the definition and order of an “index”. We mean (i,j,k) where i is the fasted changing index, but then ndarrays are indexed via [k,j,i] .

It is reasonable to apply the conventions of scikit-image for indexing, space, channel, and time order. These are compatible with how data is natively organized in ITK and NumPy array interfaces to itk.Image. That is, [k,j,i] indexing,

@blowekamp

I expect others from different domains will have additional concerns and comments to improve the specification too.

We can add documentation to a repository and use pull requests and line-based comments for additional discussion.

[jni (Juan Nunez-Iglesias)]

@matt.mccormick

We can create a reference, pure NumPy-based implementation, create tests for the attributes and methods.

To reiterate, I really think it might be better to build a spec on top of xarray rather than numpy. They exactly solve the problem of “arrays with annotations”, and we can leverage their extensive efforts to provide compatibility with dask, distributed computing, etc. Matt Rocklin specifically pulled me aside at SciPy to implore me to look long and hard at xarray, which he considers to be ready for production, before developing our own metadata handling at scikit-image.