Images
Contents
Images#
Surfa provides utilities for working with 3D and 2D multi-frame images. Image array data is structured in one of the following classes, based on spatial dimensionality:
3D volumetric images are represented by Volume
objects.
2D image slices are represented by Slice
objects.
These structures also encompass an inherent geometry, which defines the linear mapping of voxel coordinates on the image grid to Cartesian coordinates in 3D world space (sometimes referred to as scanner space).
Framed Image Basics#
From a high level, Volume
and Slice
are universal array objects that store an image buffer along with user-defined metadata. This array buffer, with a defined shape, also supports an additional non-spatial dimension to encode multiple image frames, for example, across time-series data or diffusion directions. Generally, this non-spatial dimension is represented by the last data axis, or the frame axis.
At the level of the underlying array, we can access data properties, such shape, data type, and number of frames. For example, a single-frame structural 3D MRI is loaded below as a Volume
object.
>>> import surfa as sf
>>> image = sf.load_volume('structural.nii.gz')
>>> image.shape
(256, 256, 256)
>>> image.dtype
dtype('uint8')
>>> image.nframes
1
Alternatively, a functional MRI aquistion might contain a series of timepoints stacked along the last axis. While the shape
property includes this extra dimension, baseshape
defines only the shape of the spatial domain.
>>> image = sf.load_volume('functional.nii.gz')
>>> image.shape
(256, 256, 256, 16)
>>> image.baseshape
(256, 256, 256)
>>> image.nframes
16
This same concept applies to 2D Slice
images:
>>> image = sf.load_slice('slice.nii.gz')
>>> image.shape
(256, 256, 16)
>>> image.baseshape
(256, 256)
>>> image.basedim
2
The basedim
member indicates the number of spatial dimension in the array. A Slice
always has a basedim
of 2 and a Volume
always has a basedim
of 3.
Array Representation#
At their core, Volume
and Slice
framed images operate like numpy arrays and can be constructed with such:
data = np.random.rand(128, 128, 128)
image = sf.Volume(data)
In fact, the image buffer is internally stored as a np.ndarray
instance, which can be retrieved directly via the data
member. Note that when constructing a framed array, the input numpy array is not copied by default, and so the framed array will share its memory.
Framed images implement numpy’s array-like functionality, meaning they can be implicitly converted to the np.ndarray
class when necessary, for example, when used as input to a numpy function:
mean = np.mean(image)
Framed images also emulate numpy-style operations, in-place modifications, and indexing.
image += 10
mask = image < 20
image[mask] = 0
image = image1 * image2
Array-specific operations can also be called directly as class methods.
imax = image.max()
mean = image.mean()
perc = image.percentile(99)
This is only a handful of the array utilities. Visit the Volume API reference for more detailed class documentation.
Geometry#
An important component of these framed images is the mapping of their voxel coordinates to locations in a world, or scanner, coordinate system. This relationship is defined by the geom
image geometry attribute. This ImageGeometry
object represents position in two ways: 1) using a set of linear parameters, including voxel size (resolution), image center, and rotation, and 2) using a singular vox2world
affine matrix that transforms voxel coordinates to world coordinates. When either of the parameters are updated, the affine is recomputed accordingly, and vice versa.
>>> image.geom.voxsize
array([1., 1., 1.])
>>> image.geom.center
array([0., 0., 0.])
>>> image.geom.rotation
array([[ 1., 0., 0.],
[ 0., 1., 0.],
[ 0., 0., 1.]])
>>> image.geom.vox2world
sf.Affine([
[ 1., 0., 0., -127.],
[ 0., 1., 0., -128.],
[ 0., 0., 1., -127.],
[ 0., 0., 0., 1.]])
Image Cropping#
Numpy-style indexing can be used to crop Volume
and Slice
objects while correctly accounting for the underlying geometry. When an image is cropped, as shown below, the image geometry center
and vox2world
will be updated.
>>> cropped = image[60:180, 60:180, 60:180]
>>> cropped.geom.center
array([-8., -8., -8.])
When the dimensionality of a Volume
object is reduced by 1 during cropping, the resulting image is cast as a Slice
object.
>>> image
sf.Volume(shape=(256, 256, 256), dtype=uint8)
>>> sliced = image[:, :, 110]
>>> sliced
sf.Slice(shape=(256, 256), dtype=uint8)
Reorienting#
The anatomical orientation of the voxel array is implicitly defined by the rotation matrix of the image geometry. In surfa, the default world coordinate system corresponds to Right-Anterior-Superior (RAS) anatomical ordering of the x, y, z axes. So, the image geometry orientation
can be coded by a tag that maps each anatomical axis (Left/Right, Inferior/Superior, Anterior/Posterior) to each image axis.
>>> image.geom.orientation
'RAS'
To re-order the image data to a particular orientation while preserving the correct world-space position, use the reorient
function.
>>> image = image.reorient('LIA')
Transforming and Resampling#
Other functions can be used to modify the underlying image data while correctly updating (or preserving) image geometry. To resample the image data to a particular voxel size:
>>> resized = image.resize((2, 2, 2), method='linear')
>>> resized.geom.voxsize
array([2., 2., 2.])
To fit the image to a particular shape by padding or cropping around the center:
>>> reshaped = image.reshape((180, 256, 180))
>>> reshaped.shape
(180, 256, 180)
Images can be conformed to a particular set of geometric requirements all at once using the conform
function:
>>> conformed = image.conform(voxsize=2, orientation='LIA')
To resample an image with a particular target geometry, use the resample_like
function.
>>> resampled = conformed.resample_like(image)
Affine Transformation#
To linear resample the image data with an affine transform matrix, use the transform
function.
>>> affine = sf.transform.compose_affine(rotation=(10, 0, 0))
>>> rotated = image.transform(affine, method='linear', rotation='center')
Deformation#
Voxel displacement fields can also be applied to the image to transform it nonlinearly.
>>> disp = sf.Volume(np.random.randn(*image.shape, 3))
>>> deformed = image.transform(disp, method='linear')