Create a 4D image with python bindings

Hello everybody,

I have a 4D numpy array and would like to write it as a 4D non-vector MHA image. Unfortunately, GetImageFromArray apparently cannot do that. I can edit the file header manually after writing to disk; it achieves the result I’m looking for, but I’m wondering if there is a cleaner way to do that, using SimpleITK ideally, but itk python bindings would be OK too.

Right now, the header looks like this:

ObjectType = Image
NDims = 3
BinaryData = True
BinaryDataByteOrderMSB = False
CompressedData = False
TransformMatrix = 1 0 0 0 1 0 0 0 1
Offset = 0 0 0
CenterOfRotation = 0 0 0
AnatomicalOrientation = RAI
ElementSpacing = 1 1 1
DimSize = 256 42 8
ElementNumberOfChannels = 256
ElementType = MET_SHORT
ElementDataFile = LOCAL

but I would need something like this:

ObjectType = Image
NDims = 4
BinaryData = True
BinaryDataByteOrderMSB = False
CompressedData = False
TransformMatrix = 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1
Offset = 0 0 0 0
CenterOfRotation = 0 0 0 0
AnatomicalOrientation = RAI
ElementSpacing = 1 1 1 1
DimSize = 256 256 42 8
ElementType = MET_SHORT
ElementDataFile = LOCAL

Thanks for you help.

SimpleITK in Python has a function GetImageFromArray that theoretically ought to convert a 4d numpy array to a SimpleITK image. That ought to solve your problem, but I tried it, and doesn’t work correctly for 4d image. This is a problem that needs to be rectified, obviously.

Until we can do so, here is an example that takes a 4d numpy array, extracts each 3d slice as a SimpleITK image, and then joins the slices into one 4D image.

#! /usr/bin/env python

import numpy as np
import SimpleITK as sitk

np_array = np.zeros( (10,10,10,10) )

tdim = np_array.shape[3]

slices = []

for t in range(tdim):
    slices.append( sitk.GetImageFromArray( np_array[t], False ) )

img = sitk.JoinSeries(slices)

sitk.WriteImage(img, "test.mha")
2 Likes

sitk.JoinSeries is exactly what I was after. No more ugly script to edit the MHA header!

Thanks!

Thank you for reporting this issue! There is a fix for this issue.

This enables the following:

nda = np.arange(210, dtype=np.float32).reshape([3,5,7,2])
img = sitk.GetImageFromArray(nda, isVector=True)
print(img.GetSize()) # (2,7,5,3)

More improvements to the GetImageFromArray are planned…

2 Likes