Images in physical space in Python

Engineering
1

ITK developers and long-time users understand importance of performing image processing operations in physical space, but many Python users are not aware of this and might have the temptation to work in voxel space instead.

Numpy arrays do not know anything about physical space and it seems that even scipy.ndimage and scikit-image do not store image geometry in the image object itself. Instead, there are ad-hoc solutions for passing some extra parameters to some processing functions.

Is there a widely used Python image class that can store image origin, spacing, and axis directions? Should we promote usage of that instead of (or at least in addition to) plain numpy arrays?

Should we try to convince Python image processing communities (scipy, scikit-image, …) to add support for physical mapping to their packages? Or, should we just keep ITK as the only (or few) Python image processing package that can do this properly?

5 Likes
2

In pyradiomics we use SimpleITK images, which is a bit heavy-weight but works well.

Another thing I know is nibabel which handles many neuroimaging conventions well.

3 Likes
3

My (perhaps cynical) view is that probably about 80% of the users asking to do this type of operation are doing it in preparation for “Deep Learning”. This trend has made it much harder to argue for the importance of physical space.

4 Likes
4

I think the percentage is much higher. All groups I know now uses “classic” image processing methods only for pre-processing for AI - regardless if it is actually works better or not.

Anyway, spatial normalization/denormalization steps that are usually performed on the input/output of the network, data augmentation, proper visualization, validation, etc. all require physical space information. You can surely manage with ad-hoc solutions but it would be so much simpler and safer to store the information properly where it belongs.

SimpleITK could be a good option, but it might be an issue that it is not a native Python object and pixel data is not stored in a numpy array. I know that it can be mapped to a numpy array but there are (or at least there were) limitations about what operations are allowed on the array (if the memory layout changes then things may break).

nibabel image looks quite good. Unfortunately, the package is quite neuroimaging-focused, which makes things simpler for neuroimaging people but imposes limitations for other specialties.

1 Like
5

Hi,

I’m one of the Nibabel maintainers. We’d be happy to hear
suggestions how Nibabel could be more useful for not-neuroimaging
cases …

Cheers,

Matthew

2 Likes
6

If you are aware of any limitation please post and issue on SimpleITK’s GitHub to create a bug report or feature request. We would be happy to address these issues.

7

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?

8

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?

9

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
10

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.

11

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
12

[1] Xarray documentation
[2] https://scikit-image.org/

3 Likes
13

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

14

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.

15

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.

6 Likes
16

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
17

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
18

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
19

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
20

:+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?

2 Likes