Volumes

Another fundamental class in highdicom is highdicom.Volume. Unlike many other highdicom classes, a volume is not a DICOM object at all. Instead it provides a means to work with regularly-sampled volumetric arrays conveniently. Not all DICOM images represent regularly-sampled volumetric data, but a significant proportion do, particularly many CT, MRI, and PET images. The highdicom.Volume class provides a means to work with these arrays, derive results (such as segmentations or parametric maps) from them, and store those results in new DICOM files with the correct spatial metadata without having to write code to handle spatial metadata.

A Volume object has two core components: an array of voxels with three or more dimensions, and an affine matrix that describes how that array is positioned in 3D space within the relevant frame-of-reference coordinate system (either the patient or slide coordinate system).

Creating a Basic Volume

To create a basic volume we need to provide the affine matrix and the array, as well as specifying the coordinate system being used.

import numpy as np
import highdicom as hd


array = np.zeros((32, 256, 256), dtype=np.uint8)
affine = np.array(
    [
        [1.0, 0.0, 0.0, 10.0],
        [0.0, 0.4, 0.0, 20.0],
        [0.0, 0.0, 0.4, -5.0],
        [0.0, 0.0, 0.0,  1.0],
    ]
)

vol = hd.Volume(
    array=array,
    affine=affine,
    coordinate_system="PATIENT",
)

# Once created, you can access the array and the affine via properties
assert np.array_equal(vol.array, array)
assert np.array_equal(vol.affine, affine)

# The datatype of the array, its shape, and the spatial shape (not
# including any channels, see below), may be accessed via
# properties
print(vol.dtype)
# uint8

print(vol.shape)
# (32, 256, 256)

print(vol.spatial_shape)
# (32, 256, 256)

Affine Matrices

The affine matrix is a 4x4 numpy array that can be used to map (zero-based) array indices into the position of the center of the indexed voxel in the frame-of-reference coordinate system.

Within highdicom, the affine matrices used within volumes are always stored using the “LPS” convention for describing the frame of reference, since this is the convention used within DICOM itself. Affine matrices defined in other conventions need to be converted before being used to describe a volume. Another common convention in the “RAS” convention used within NIfTI files. If you have an affine matrix generated elsewhere using a different convention, you can specify the from_reference_convention parameter to specify the convention used by the source affine, and highdicom will handle the conversion to “LPS” convention for you.

The top-left 3x3 matrix of the affine matrix consists of three orthogonal column vectors that give the vector travelled when incrementing the voxel indices in each o the three dimensions. The top three elements of the last column give the position of the voxel at index (0, 0, 0). The bottom row is always [0., 0., 0., 1.]. The highdicom.Volume.from_components() method allows you to construct a volume by specifying these components of the affine matrix. Furthermore, once constructed you can use properties and methods of the object to access these components individually.

This code snippet produces the same volume as in the previous that specified the full affine matrix.

import numpy as np
import highdicom as hd


array = np.zeros((32, 256, 256), dtype=np.uint8)
direction = np.eye(3)
spacing = np.array([1.0, 0.4, 0.4])
position = np.array([10.0, 20.0, -5.0])

vol = hd.Volume.from_components(
    array=array,
    direction=direction,
    spacing=spacing,
    position=position,
    coordinate_system="PATIENT",
)

print(vol.direction)
# [[1. 0. 0.]
#  [0. 1. 0.]
#  [0. 0. 1.]]

print(vol.spacing)
# (1.0, 0.4, 0.4)

print(vol.position)
# (10.0, 20.0, -5.0)

print(vol.center_position)
# (25.5, 71.0, 46.0)

print(vol.unit_vectors())
# (array([1., 0., 0.]), array([0., 1., 0.]), array([0., 0., 1.]))

Volume Geometries

Sometimes it is useful to work with the affine matrix of a volume without the full voxel array. The highdicom.VolumeGeometry fills this role. It has an API that is compatible with highdicom.Volume, except for any operations that require access to the pixel data.

import numpy as np
import highdicom as hd


direction = np.eye(3)
spacing = np.array([1.0, 0.4, 0.4])
position = np.array([10.0, 20.0, -5.0])

geometry = hd.VolumeGeometry.from_components(
    spatial_shape=(32, 256, 256),
    direction=direction,
    spacing=spacing,
    position=position,
    coordinate_system="PATIENT",
)

print(geometry.direction)
# [[1. 0. 0.]
#  [0. 1. 0.]
#  [0. 0. 1.]]

print(geometry.spacing)
# (1.0, 0.4, 0.4)

print(geometry.position)
# (10.0, 20.0, -5.0)

print(geometry.center_position)
# (25.5, 71.0, 46.0)

print(geometry.unit_vectors())
# (array([1., 0., 0.]), array([0., 1., 0.]), array([0., 0., 1.]))

Volumes From Images

Volumes are loaded from existing images more often than constructed directly. To load a volume from a single DICOM image (single frame or multi-frame), use the highdicom.Image.get_volume() method (see Images).

from pydicom.data import get_testdata_file

import highdicom as hd

# Load an enhanced (multiframe) CT image
im = hd.imread(get_testdata_file('eCT_Supplemental.dcm'))

geometry = im.get_volume_geometry()

assert geometry is not None

vol = im.get_volume()
print(vol.spatial_shape)
# (2, 512, 512)

print(vol.affine)
# [[   0.          0.         -0.388672   99.5     ]
#  [  -0.          0.388672    0.       -301.5     ]
#  [  10.          0.          0.       -159.      ]
#  [   0.          0.          0.          1.      ]]

Even if the image consists of a single plane, the resulting Volume will have three spatial dimensions and the singleton dimension is placed first.

Volumes From Image Series

In the case where the frames that make up a volume are stored across multiple, single-frame files from a series, the highdicom.get_volume_from_series() function may be used to create a volume.

import pydicom
from pydicom.data import get_testdata_file

import highdicom as hd

# Three test files from pydicom that form a volume
ct_files = [
    get_testdata_file('dicomdirtests/77654033/CT2/17136'),
    get_testdata_file('dicomdirtests/77654033/CT2/17196'),
    get_testdata_file('dicomdirtests/77654033/CT2/17166'),
]
ct_series = [pydicom.dcmread(f) for f in ct_files]

vol = get_volume_from_series(ct_series)

Array Manipulation

Since the volume’s array is just a NumPy array, it can be manipulated just like any other numpy array to process the image. However, any operation that changes the array’s shape is not allowed because changing the shape requires changing the affine matrix.

import numpy as np
import highdicom as hd


vol = hd.Volume.from_components(
    array=np.random.randint(0, 100, size=(32, 256, 256), dtype=np.uint8),
    direction=np.eye(3),
    spacing=[1., 0.4, 0.4],
    position=[10., 20., 30.],
    coordinate_system="PATIENT",
)

# OK
vol.array = vol.array + 10.0

# OK
vol.array /= 100

# OK
vol.array = np.exp(vol.array / 1000)

# Disallowed, changes shape
vol.array = vol.array[:10]

The above operations edit the volume in-place. If you want to create a new volume with a new array but the same geometry as an existing volume, use the highdicom.Volume.with_array() method.

Indexing

Volumes can be indexed along their spatial dimensions using square brackets in a largely similar way to any NumPy array. This operation crops the array and also updates the affine matrix to reflect the effect of the crop. However, there is one important change: spatial dimensions can be reduced to size one but never removed by indexing (volumes always have three spatial dimensions).

import numpy as np
import highdicom as hd


vol = hd.Volume.from_components(
    array=np.random.randint(0, 100, size=(32, 256, 256), dtype=np.uint8),
    direction=np.eye(3),
    spacing=[1., 0.4, 0.4],
    position=[10., 20., 30.],
    coordinate_system="PATIENT",
)

cropped = vol[:10]
print(cropped.shape)
# (10, 256, 256)

cropped = vol[10]
print(cropped.shape)
# (1, 256, 256)

cropped = vol[:, 20:100, -80:]
print(cropped.shape)
# (32, 80, 80)

cropped = vol[:, :, 200:120:-1]
print(cropped.shape)
# (32, 256, 80)

Spatial Operations

The highdicom.Volume class provides a number of spatial operations that manipulate the array and correctly update the affine matrix to reflect the change. Currently these only include operations that do not require resampling of the array:

Patient Orientation

For volumes in the patient frame-of-reference coordinate system, the “patient orientation” describes how the axes of the volume align with the axes of the patient coordinate system, which are defined from left-to-right, anterior-to-posterior, and foot-to-head, in that order. The axes of the volume do not need to be exactly aligned with the frame-of-reference axes to be described using a patient orientation; if they are not the closest match is used. For example, the patient orientation “FPL” means that the first axis of the volume is most closely aligned with the head-to-foot direction, the second axis of the volume is most closely aligned with the anterior-to-posterior direction, and the third axis is most closely aligned with the right-to-left direction.

Patient orientations may be used to describe a volume, and the highdicom.Volume.to_patient_orientation() is used to manipulate a volume to align with the given patient orientation as well as possible via permutations and flips.

Patient orientations may be represented as strings or as tuples of the highdicom.PatientOrientationValuesBiped class.

from pydicom.data import get_testdata_file

import highdicom as hd

# Load an enhanced (multiframe) CT image
im = hd.imread(get_testdata_file('eCT_Supplemental.dcm'))

vol = im.get_volume()

print(vol.get_closest_patient_orientation())
# (<PatientOrientationValuesBiped.H: 'H'>, <PatientOrientationValuesBiped.P: 'P'>, <PatientOrientationValuesBiped.R: 'R'>)

vol = vol.to_patient_orientation("LAF")

print(vol.get_closest_patient_orientation())
# (<PatientOrientationValuesBiped.L: 'L'>, <PatientOrientationValuesBiped.A: 'A'>, <PatientOrientationValuesBiped.F: 'F'>)

Channels

In addition to the three spatial dimensions, a volume may have further non-spatial dimensions that are referred to as “channels”. Channel dimensions are stacked after the spatial dimensions in the volume’s pixel array. The meaning of each channel is explicitly described in the volume. Common uses for channels include RGB channels in color images, optical paths in microscopy images, or contrast phases in radiology images.

The highdicom.ChannelDescriptor class is used to describe the meaning of a single channel dimension. Where possible, it is recommended to use DICOM attributes to describe channels. A DICOM keyword or the corresponding tag value may be passed to the highdicom.ChannelDescriptor constructor.

When using a DICOM attribute, each channel of the volume is associated with a particular value for that attribute. For example, if the descriptor uses the “OpticalPathIdentifier” attribute, each channel will be associated with a string. Alternatively if an integer-valued attribute like “SegmentNumber” is used, each channel will be associated with an integer. We refer to this type as the descriptor’s “value type”.

This code snippet creates channel descriptors using some DICOM attribute, and checks the corresponding value types:

import highdicom as hd


# Channel descriptor using the "OpticalPathIdentifier"
optical_path_descriptor = hd.ChannelDescriptor('OpticalPathIdentifier')

# Using the hexcode for the attribute is equivalent
optical_path_descriptor = hd.ChannelDescriptor(0x0048_0106)

# Channel descriptor using the "DiffusionBValue"
bvalue_descriptor = hd.ChannelDescriptor('DiffusionBValue')

# Check that the value types are as expected
print(optical_path_descriptor.value_type)
# <class 'str'>

print(bvalue_descriptor.value_type)
# <class 'float'>

Alternatively, it is possible to define custom identifiers that do not use a DICOM attribute. In this case, you must specify the value type yourself. The value type must be either int, str, or float (or a sub-type of one of these types), or an enumerated type derived from the Python standard library enum.Enum.

from enum import Enum
import highdicom as hd

# A custom descriptor using integer values
custom_int_descriptor = hd.ChannelDescriptor(
    'my_int_descriptor',
    is_custom=True,
    value_type=int,
)

# A custom descriptor using an enumerated type
class MyEnum(Enum):
    VALUE1 = "VALUE1"
    VALUE2 = "VALUE2"

custom_enum_descriptor = hd.ChannelDescriptor(
    'my_enum_descriptor',
    is_custom=True,
    value_type=MyEnum,
)

One very common channel descriptor that does not correspond to a DICOM attribute is RGB color channels. The enum highdicom.RGBColorChannels is used as the value type for volumes with color channels, and the descriptor for this channel is provided as a constant in highdicom.RGB_COLOR_CHANNEL_DESCRIPTOR.

To create a volume with channels, you must provide a dictionary that contains, for each channel dimension, the channel descriptor and the values of each channel along that dimension:

import numpy as np
import highdicom as hd

# Array with three spatial dimensions plus 3 color channels and 4 optical
# paths
array = np.random.randint(0, 10, size=(1, 50, 50, 3, 4))

# Names of the 4 optical paths
path_names = ['path1', 'path2', 'path3', 'path4']

vol = hd.Volume.from_components(
    direction=np.eye(3),
    center_position=[98.1, 78.4, 23.1],
    spacing=[2.0, 0.5, 0.5],
    coordinate_system="SLIDE",
    array=array,
    channels={
        hd.RGB_COLOR_CHANNEL_DESCRIPTOR: ['R', 'G', 'B'],
        'OpticalPathIdentifier': path_names
    },
)

# The total shape of the volume includes the channel dimensions
assert vol.shape == (1, 50, 50, 3, 4)

# But the spatial shape excludes them
assert vol.spatial_shape == (1, 50, 50)

# The channel shape includes only the channel dimensions, not the spatial
# dimensions
assert vol.channel_shape == (3, 4)
assert vol.number_of_channel_dimensions == 2

# You can access the descriptors like this
assert vol.channel_descriptors == (
    hd.RGB_COLOR_CHANNEL_DESCRIPTOR,
    hd.ChannelDescriptor('OpticalPathIdentifier'),
)

The order of the items in the dictionary is significant and must match the order of the channel dimensions in the array.

For most purposes, a volume with channels can be treated just like one without. All spatial operations (including indexing) only alter the array along the spatial dimensions and leave the channel dimensions unchanged. A separate set of methods are used to alter the channel dimensions:

This snippet, using the same volume as above, demonstrates how to use these methods:

import numpy as np
import highdicom as hd

# Array with three spatial dimensions plus 3 color channels and 4 optical
# paths
array = np.random.randint(0, 10, size=(1, 50, 50, 3, 4))

# Names of the 4 optical paths
path_names = ['path1', 'path2', 'path3', 'path4']

vol = hd.Volume.from_components(
    direction=np.eye(3),
    center_position=[98.1, 78.4, 23.1],
    spacing=[2.0, 0.5, 0.5],
    coordinate_system="SLIDE",
    array=array,
    channels={
        hd.RGB_COLOR_CHANNEL_DESCRIPTOR: ['R', 'G', 'B'],
        'OpticalPathIdentifier': path_names
    },
)

assert (
    vol.get_channel_values('OpticalPathIdentifier') ==
    path_names
)

# Get a new volume containing just optical path 'path2'
path_2_vol = vol.get_channel(OpticalPathIdentifier='path2')

# Swap the two channel axes by descriptor
permuted_vol = vol.permute_channel_axes(
    ['OpticalPathIdentifier', 'RGBColorChannel']
)

# Swap the two channel axes by index
permuted_vol = vol.permute_channel_axes_by_index([1, 0])

Full Example

This full example presents a typical workflow of how volumes are used within highdicom. First, a volume is extracted from an existing image. Then it is manipulated to prepare it for some automated analysis tool (a simple example segmentation in this case). The tool’s output is placed back into a volume, which is then passed to the constructor of a highdicom class to ensure that the spatial metadata in the output object is correct.

import numpy as np

from pydicom.sr.codedict import codes
from pydicom import pixel_array
from pydicom.data import get_testdata_file
from pydicom.uid import JPEGLSLossless

import highdicom as hd


def complex_segmentation_tool(arr: np.ndarray) -> np.ndarray:
    """This is a stand-in for a generic segmentation tool.

    We assume that the tool has certain requirements on the input array, in
    this case that it has patient orientation "FLP" and a shape of (2, 400,
    400).

    Further, we assume that the tool takes in a numpy array and returns a
    binary segmentation that is pixel-for-pixel aligned with its input array
    (i.e. the tool itself does not do any further spatial manipulation.

    """
    # Basic thresholding as a simple example
    return arr > 0

# Load an enhanced (multiframe) CT image
im = hd.imread(get_testdata_file('eCT_Supplemental.dcm'))

# Load the input volume
original_volume = im.get_volume()

# Manipulate the original volume to give a suitable input for the tool
input_volume = (
    original_volume
    .to_patient_orientation("FLP")
    .crop_to_spatial_shape((2, 400, 400))
)

# Run the "complex segmentation tool"
seg_array = complex_segmentation_tool(input_volume.array)

# Since the seg array shares its geometry with the inupt array, we can combine
# the two to create a volume of the segmentation array
seg_volume = input_volume.with_array(seg_array)

algorithm_identification = hd.AlgorithmIdentificationSequence(
    name='Complex Segmentation Tool',
    version='v1.0',
    family=codes.cid7162.ArtificialIntelligence
)

# metadata needed for a segmentation
brain_description = hd.seg.SegmentDescription(
    segment_number=1,
    segment_label='brain',
    segmented_property_category=codes.SCT.Organ,
    segmented_property_type=codes.SCT.Brain,
    algorithm_type=hd.seg.SegmentAlgorithmTypeValues.AUTOMATIC,
    algorithm_identification=algorithm_identification,
)

# Use the segmentation volume as input to create a DICOM Segmentation
seg_dataset = hd.seg.Segmentation(
    pixel_array=seg_volume,
    source_images=[im],
    segmentation_type=hd.seg.SegmentationTypeValues.LABELMAP,
    segment_descriptions=[brain_description],
    series_instance_uid=hd.UID(),
    series_number=1,
    sop_instance_uid=hd.UID(),
    instance_number=1,
    manufacturer='Complex Segmentations Plc.',
    manufacturer_model_name='Complex Segmentation Tool',
    software_versions='0.0.1',
    device_serial_number='1234567890',
    transfer_syntax_uid=JPEGLSLossless,
    series_description='Example Segmentation of CT',
)

seg_dataset.save_as('segmentation.dcm')

# Alternatively, it may be desirable to match the geometry of the output
# segmentation image to that of the input image. This will "undo" the
# cropping and axis permutation operations done to the image volume above.
seg_volume_matched = seg_volume.match_geometry(original_volume)

# Use the segmentation volume as input to create a DICOM Segmentation
seg_dataset_matched = hd.seg.Segmentation(
    pixel_array=seg_volume_matched,
    source_images=[im],
    segmentation_type=hd.seg.SegmentationTypeValues.LABELMAP,
    segment_descriptions=[brain_description],
    series_instance_uid=hd.UID(),
    series_number=1,
    sop_instance_uid=hd.UID(),
    instance_number=1,
    manufacturer='Complex Segmentations Plc.',
    manufacturer_model_name='Complex Segmentation Tool',
    software_versions='0.0.1',
    device_serial_number='1234567890',
    transfer_syntax_uid=JPEGLSLossless,
    series_description='Example Segmentation of CT',
)

seg_dataset_matched.save_as('segmentation_matched.dcm')

Volumes To/From NIfTI Files

NIfTI is a file format used to store volumetric imaging data. It arose from neuro-imaging but is now used in other areas of radiology and beyond. When converting between highdicom Volumes and NIfTI files, it is critical to remember to account for the difference in convention used to specify the frame-of-reference coordinate system: highdicom (and DICOM) uses “LPS” convention, NIfTI uses “RAS” convention.

We plan to add tools to handle this conversion in the near future, but for now these snippets should correctly handle simple situations converting to and from NIfTI using the nibabel package.

Reading a volume from a NIfTI:

import nibabel as nib
import highdicom as hd


nifti_path = '/path/to/nifti.nii'  # or .nii.gz
nifti = nib.load(nifti_path)

vol = hd.Volume(
    array=nifti.get_fdata(),
    affine=nifti.affine,
    coordinate_system="PATIENT",
    from_reference_convention='RAS',
)

Writing a volume to a NIfTI file:

import nibabel
import highdicom as hd


vol = hd.Volume(...)

nifti = nib.Nifti1Image(
    vol.array,
    vol.get_affine('RAS'),
)

nifti_path = '/path/to/nifti.nii'  # or .nii.gz
nib.save(nifti, nifti_path)

Volumes To/From ITK Images

ITK is a widely-used library for volumetric image processing. Its Image class shares many similarities with our highdicom.Volume class. Like highdicom, ITK uses the “LPS” convention. However, when converting to and from NumPy arrays, ITK reverses the order of dimensions. It is important to account for this when performing conversions.

We plan to add tools to handle this conversion in the near future, but for now these snippets should correctly handle simple situations converting to and from ITK Images.

Creating a volume from an ITK Image:

import itk
import numpy as np
import highdicom as hd


im = itk.image(...)

# Reverse array dimension order
array = itk.array_from_image(im).transpose([2, 1, 0])

vol2 = hd.Volume.from_components(
    array=array,
    direction=np.asarray(im.GetDirection()),
    spacing=np.asarray(im.GetSpacing()),
    position=np.asarray(im.GetOrigin()),
    coordinate_system="PATIENT"
)

Creating an ITK Image from a Volume:

import itk
import highdicom as hd


vol = hd.Volume(...)

# Reverse array dimension order
array = vol.array.transpose([2, 1, 0])

im = itk.image_from_array(array)
im.SetOrigin(vol.position)
im.SetDirection(vol.direction)
im.SetSpacing(vol.spacing)