Images in physical space in Python

I remember reading about limitations about reallocating a numpy array that is mapped from a SimpleITK image. For example, if I create a SimpleITK image, get it as a numpy array, and modify the numpy array (reshape and resize the array and change some voxel values), then is the SimpleITK image going to remain valid?

My main concern is that it seems that nibabel seems to use uses nifti image object to store images. Nifti format has limitations on which dimensions of an image may contain (see for example this discussion) and it does not support arbitrary metadata fields. Does nibabel.nifti1.Nifti1Image class have these limitations, too?

Hi,

Yes, that’s true - we do concentrate on NIfTI. We have thought of
extending NIfTI with JSON metadata in the header, and we could revisit
that.

The restrictions on the image axes are, I agree, pretty severe. I
wonder whether something like NRRD would be a good basis for a more
general image format. It’s a pretty simple format - and I never get
tired of the name.

Cheers,

Matthew

2 Likes

SimpleITK has a conservative model for interfacing with numpy with a focus on easy to use methods that are safe. There are GetImageFromArray and GetArrayFromImage methods which do a deep-copy, and a GetArrayViewFromImage method which provides “read-only” access to the itk Image buffer. If you have specific code of the issue please share, but the SimpleITK interface is designed to be safe and protected from this type of aliasing.

Hi @lassoan!

A great question, and one that we haven’t yet resolved at scikit-image (see here for a sample discussion). The feedback from our roadmap post last year was that continuing to support bare NumPy arrays as the first-class citizen was essential (and the core team agrees), and that people were generally OK with scikit-image being the low-level voxel-space library on top of which physical-space libraries could be built. We might revisit this question periodically, but I expect that this is a post-1.0 question.

In the meantime, from what I can tell, the most flexible, “SciPythonic” library to work with physical space arrays is xarray [1], and I suspect that if/when we start supporting physical space in scikit-image, it’ll be through xarray or some close derivative.

We are certainly amenable to further discussion, either here or on any of the communication channels listed on our home page [2]! (You will reach more of our community there — thanks @ctrueden for bringing this post to our attention!)

…[1, 2]: Discourse annoyingly only allows two links in posts for new users, so links upcoming in next post

5 Likes

[1] http://xarray.pydata.org/en/stable/
[2] https://scikit-image.org/

3 Likes

Thanks for the note @jni – this limit has been increased to 10.

I heard MINC has full HDF5 backing for proper slicing and chunking, plus extensive metadata support and compression all designed in from the start.

I keep wondering why people want to “extend” existing formats which have missing features instead of implementing a solution that already has a complete set.

nibabel doesn’t use libminc for example. It has a hacky HDF5 wrapper which doesn’t enable any of the good bits.

The itk Python package is gaining support for NumPy arrays and NumPy-like arrays. And, NumPy is adding support for additional NumPy-like arrays, “duck arrays”, like Dask arrays, Zarr arrays, xarray, CuPy arrays, … we have an opportunity to improve the image processing experience across the scientific Python ecosystem by following the principle of compatible, NumPy-based interfaces.

We can add critical spatial support by extending the ndarray interface with three spatial image attributes:

  • origin: the physical distance in the coordinate system to the center of the “lower left” pixel.
  • spacing: the N-dimensional physical spacing between pixels.
  • direction: an orientation matrix consisting of direction cosines.

This would make

  • Resampling
  • Registration
  • Multi-resolution

image processing elegant and effective in Python.

Soon, itk.Image will be a partial duck array. And, ITK image filters will support passing and returning NumPy duck arrays. When NumPy duck arrays are passed, if hasattr(image, 'spacing'): [...] can be used to check and use spacing information, if present. The returned object from an itk filter will have origin, spacing, and direction attributes set accordingly, whether an itk.Image type or a spatial image type, i.e. has origin, spacing, and direction attributes, derived from numpy.ndarray.

After a few issues are implemented, we will have the basics covered, and we can create a pre-release to test and refine.

It is absolutely wonderful to have folks from the broader community like @matthew.brett, @jni, @lassoan, @ctrueden, @pieper participating in the discussion. Ideally, additional tools in the community incrementally support a NumPy-based spatial image interface – an interface approach is nice because it is not tied to an implementation, can be adopted incrementally, and is optional. Ideally, scikit-image would utilize spatial attributes, if present, as more spatial functionality is added to scikit-image, SimpleITK Python starts to support these attributes, pyimagej transfers these attributes when converting an ImgLib2 image to NumPy with imglyb, 3D Slicer supports these attributes when working with volume nodes, and nibabel provides the interface. We also discussed the issue of preserving spatial metatdata in the imageio package and its interface.

5 Likes

This sounds like a really great an exciting approach. It would be great if these libraries define the attributes the same way so they are interoperable. That I believe would be the critical thing ensuring that all libraries that use these attribute define them the same way. I see some critical points:

  • What coordinations system LPS? ( from right towards left, …) ITK uses LPS while VTK/Slicer uses RAS.
  • Where is the origin defined? ITK uses the center of the 0-indexed voxel, other libraries use the corner.
  • Specifically defining the mapping from index space to physical space.
3 Likes

I think most medical software has standardized on LPS.

In Slicer, we are slowly moving towards LPS, too. Input/output files are already assumed to be in LPS (with a few exceptions that we are going to address soon by storing coordinate system information in the file header). Within 1-2 years we’ll probably switch the internal coordinate system to LPS, too.

DICOM, ITK, VTK, nrrd, metaimage, etc. all uses center of the voxel as physical position within a 3D volume.

Overall, for medical image computing, nowadays things are quite well standardized. However, other application areas may use conventions. For example, communities that only work with 2D images may find pixel origin in the corner a more natural choice.

The main focus of this discussion is about how to represent image data in memory in Python (it should be possible to read/write the image data to/from any file format). How to keep the number of research file formats at the minimum would be an interesting discussion, too - you could start a new topic on that.

This sounds awesome.

Great! Is this planned for SimpleITK, too? (Slicer uses SimpleITK’s Python wrapping)

2 Likes

It looks like it is still under development in NumPy with NEP 30, so we’ll follow the development behind the tail blazers until it’s reader for prime time. SimpleITK is committed to continual improving the interface with NumPy and interoperability.

1 Like

Another approach I’ve heard used is array + affine transform matrix, which gives you a bit of extra flexibility.

But, in general, awesome idea to standardise on a minimum amount of extra information.

2 Likes

:+1: :+1: for establishing this metadata standard across image libraries in Python.

Would it then make sense to create some validation code, which implementors can use to check their image objects for appropriate structure? Is there a solid Python library for doing “duck-type conformance checking” like this? Maybe this SO post is helpful?

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?

1 Like

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.
3 Likes

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. :duck: :microphone:

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.

1 Like

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

2 Likes

Fully agree.

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.

2 Likes

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,

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

1 Like

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.

1 Like