dicomstack

A thin layer over pydicom to load DICOM data (folders, zip files, multi-frame files), explore and query it by field values, sort it, and build 3d volumes with their geometry.

At a glance:

from dicomstack import DicomStack, DICOM

stack = DicomStack("data")              # load a folder (or a file, a zip file...)
print(stack.describe())                 # list its studies and series

stack.unique("EchoTime")                # -> [5, 10, 20]: values of a field
t2 = stack(SeriesNumber=2)              # filter by field values
first_echo = t2(DICOM.EchoTime < 15)    # ... or with queries

volume = first_echo.as_volume()         # numpy array (x, y, z) with its geometry
volume.spacing                          # -> (1.484375, 1.484375, 2.0)
echo_times, volumes = t2.as_volume(by="EchoTime")   # one volume per echo time

The examples below use a folder data/ with two series (a T1 series of 4 slices, and a T2 series with 2 echoes of 4 slices) and a non-DICOM file, the same data zipped in data.zip, and a multi-frame file multiframe.dcm.

Installation

pip install dicomstack
pip install "dicomstack[pandas]"    # optional: stack.to_dataframe()

Requires Python ≥ 3.10, pydicom ≥ 3 and numpy.

Loading

from dicomstack import DicomStack

stack = DicomStack("data")                    # a folder and its sub-folders
stack = DicomStack("data/series1/IM1.dcm")    # a single file
stack = DicomStack(["data/series1", "data/series2/IM1_1.dcm"])  # several paths
stack = DicomStack("data/series*")            # glob pattern
stack = DicomStack("data.zip")                # a zip file
stack = DicomStack("data.zip/series1")        # a folder inside a zip file
stack = DicomStack(filenames=["data/series1/IM1.dcm", "data/series1/IM2.dcm"])

stack = DicomStack("data", include="*.dcm")   # file name patterns (fnmatch)
stack = DicomStack("data", exclude=["*.txt", "series2/*"], recursive=True)

stack = DicomStack("data")
stack.root                                    # folder of the files
stack.filenames[:2]                           # -> ['series1/IM1.dcm', 'series1/IM2.dcm'] (relative to root)
stack.non_dicom                               # -> ['notes.txt']: files which are not DICOM

# multi-frame (enhanced) files: one frame per image
multi = DicomStack("multiframe.dcm")
print(multi)                                  # -> DICOM(90 frames | 1 study | 1 series)

Loading reads the files' headers only; field values are converted when accessed, and pixels are read when needed. The values of private fields are not kept in memory: they are read again from the file when accessed (except for zip members). A progress bar is shown while loading when the output is a terminal (progress=False to hide it, True to always show it, or a function progress(done, total)).

Duplicates. Frames are identified by their SOPInstanceUID (and frame index): the same image loaded twice (two copies of a file, overlapping paths) appears once. duplicates="error" raises DuplicatedFramesError instead, duplicates="ignore" keeps them. dedup="content" compares whole file contents instead of UIDs (slower: whole files are read); files without SOPInstanceUID are never merged.

stack = DicomStack(["data/series1", "data/series1"])       # same files twice
len(stack)                                    # -> 4
stack = DicomStack("data", duplicates="error", dedup="content")

Header cache. Re-opening the same data parses all headers again. cache=True stores the converted values in a file (.dicomstack-cache.json.gz in the data folder, or in the user cache directory if the folder is not writable), used on the next loads for the files that are unchanged (same size and modification time). The first load with cache=True is slower (all values are converted); the next ones are much faster. cache="refresh" rebuilds it, cache="path/to/file.json.gz" sets its location.

stack = DicomStack("data", cache=True)        # writes data/.dicomstack-cache.json.gz
stack = DicomStack("data", cache=True)        # unchanged files: not parsed

Threads. By default, 8 threads read the files ahead (workers=8, also for as_volume and iter_pixels), while headers and pixels are decoded one by one: much faster on network storage, about the same on a local disk. workers=None reads the files sequentially.

In-memory datasets.

import pydicom

datasets = [pydicom.dcmread("data/series1/IM1.dcm"), pydicom.dcmread("data/series1/IM2.dcm")]
stack = DicomStack.from_datasets(datasets)

Exploring

stack = DicomStack("data")
print(stack.describe())

0: (non-DICOM data)
	0: notes.txt
1: Legs
	1: T1 axial (n=4)
	2: T2 multi-echo (n=8)

Studies are listed by StudyID (with a suffix if several studies share it), series by SeriesNumber. stack.dicomtree returns the same structure as a dictionary.

stack.dicomtree.keys()                        # -> dict_keys(['0', '1']): studies ("0": non-DICOM files)
len(stack.fields())                           # -> 116: keywords present in all frames
stack.fields("any")                           # keywords present in any frame
stack.varying()                               # -> ['EchoTime', 'SeriesDescription', 'SeriesInstanceUID', 'SeriesNumber']
stack.varying(within="SeriesNumber")          # -> ['EchoTime']: fields varying inside a series
stack.table("SeriesNumber", "EchoTime")[:2]   # -> [{'SeriesNumber': 1, 'EchoTime': 5}, {'SeriesNumber': 1, 'EchoTime': 5}]

varying() ignores per-image fields (SOPInstanceUID, InstanceNumber, positions, window values...: see dicomstack.dicomstack.VARYING_EXCLUDE; exclude=() to include them).

df = stack.to_dataframe("SeriesNumber", "EchoTime", "SliceThickness")  # requires pandas

Field access

from dicomstack import DICOM

stack["EchoTime"]                             # one value per frame
stack["SeriesNumber", "EchoTime"][:2]         # -> [(1, 5), (1, 5)]
frame = stack[0]                              # a frame (stack[1:3]: a sub-stack)
frame["EchoTime"]                             # -> 5
frame.get("StudyComments")                    # -> None (missing field)

Field names:

stack["ImageType_1"][0]                       # -> 'ORIGINAL'
stack[DICOM.ImageType[0]][0]                  # -> 'ORIGINAL'
stack["IconImageSequence.1.Rows"][0]          # -> 64
stack["[CSAImageHeaderType]"][0]              # -> 'IMAGE NUM 4'
stack[(0x0018, 0x0081)][0]                    # -> 5

Values are converted: dates and times to ISO strings ('2014-11-24', '14:11:39.046000'), decimal strings to int when integral (5) or float, multi-valued fields to tuples (PixelSpacing → (1.484375, 1.484375)), empty values to None, sequences to lists of datasets. frame.dataset returns all the converted elements of a frame (frame.dataset.elements, elements_by_keyword, ...).

stack.unique("SeriesNumber")                  # -> [1, 2] (sorted, missing values ignored)
stack.unique("SeriesNumber", "EchoTime")      # -> [(1, 5), (2, 10), (2, 20)]
stack(SeriesNumber=2).single("SeriesDescription")  # -> 'T2 multi-echo'
stack.has_field("EchoTime")                   # -> True (all frames; how="any")
"EchoTime" in stack                           # -> True

single() raises MultipleValuesError if the values differ and NoValueError if there is none (both are ValueError), unless a default is given:

from dicomstack import MultipleValuesError

try:
    stack.single("EchoTime")
except MultipleValuesError as exc:
    print(exc)                                # -> Multiple values found for ('EchoTime',)
stack.single("StudyComments", default=None)   # -> None

Filtering and queries

stack(SeriesNumber=2)                         # frames with SeriesNumber == 2
stack(SeriesNumber=2, EchoTime=[10, 20])      # several fields, list of accepted values

Queries are built from selectors (DICOM.Keyword, DICOM["[Private Name]"], DICOM[(group, element)], or stack.Keyword):

stack(DICOM.EchoTime > 5)
stack(DICOM.ImageType[1] == "PRIMARY")
stack(DICOM.SeriesNumber.isin([1, 2]))
stack(DICOM.SeriesDescription.startswith("T2"))
stack((DICOM.SeriesNumber == 2) & ~(DICOM.EchoTime == 20))   # &, |, ~
stack(stack.EchoTime.between(5, 10))

A comparison which fails (e.g. a missing value compared with <, or a number compared with a string) excludes the frame with a warning; strict=True raises the error instead:

stack(DICOM.SeriesNumber > "1")               # warning: no frame
stack(DICOM.SeriesNumber > "1", strict=True)  # TypeError

Sorting and grouping

stack.sort("SeriesNumber", "InstanceNumber")  # frames with missing values are dropped
stack.sort("StudyComments", missing="last")   # ... or kept at the end ("error": NoValueError)
stack(SeriesNumber=1).reorder()               # sorted along the slice direction

for (series, echo_time), substack in stack.groupby("SeriesNumber", "EchoTime"):
    print(series, echo_time, len(substack))   # -> 1 5 4, 2 10 4, 2 20 4

for uid, series in stack.series():            # groupby("SeriesInstanceUID")
    print(series.single("SeriesDescription"))

stack[0], stack[-1], stack[2:6]               # frame, frame, sub-stack
t1, t2 = stack(SeriesNumber=1), stack(SeriesNumber=2)
len(t1 + t2), len(stack - t1), len(stack & t2)    # -> (12, 8, 8): union, difference, intersection

Volumes

volume = stack(SeriesNumber=1).as_volume()
volume.shape                                  # -> (256, 192, 4): (x, y, z) = (columns, rows, slices)
volume.origin                                 # -> (-190, -142.5, 2): position of voxel (0, 0, 0)
volume.spacing                                # -> (1.484375, 1.484375, 2.0)
volume.transform                              # -> ((1, 0, 0), (0, 1, 0), (0.0, 0.0, 1.0))
volume.affine                                 # 4x4: voxel indices (i, j, k, 1) -> coordinates

Volumes are numpy arrays (DicomVolume) with their geometry, in LPS coordinates (as DICOM and ITK): transform[i] is the direction of array axis i and spacing[i] the voxel size along it. Slices are sorted by position (reorder=False keeps the stack's order), pixel values are rescaled with RescaleSlope/RescaleIntercept (rescale=False to keep the stored values), the data type is float64 (dtype="float32" halves the memory). A warning is shown if the slices are not evenly spaced (missing slice), strict=True raises a ValueError instead.

Slicing a volume keeps a consistent geometry (volume[:, :, 2:3].origin is the position of the third slice). To write volumes to image files, see With tomio.

echo_times, volumes = stack(SeriesNumber=2).as_volume(by="EchoTime")
echo_times                                    # -> [10, 20]: one volume per echo time
echo_times, volume = stack(SeriesNumber=2).as_volume(by="EchoTime", stack=True)
volume.shape                                  # -> (256, 192, 4, 2): echo times on the last axis
volume.metadata                               # -> {'EchoTime': '[10, 20]'}

for pixels in stack.iter_pixels():            # frames' pixel arrays (rows, columns)
    pass

With tomio

Volumes have the attributes of tomio's Volume (spacing, origin, transform, metadata): tomio writes them to image files (NIfTI, NRRD, MetaImage...) with their geometry.

import tomio

echo_times, volume = stack(SeriesNumber=2).as_volume(by="EchoTime", stack=True)
tomio.write_volume("multiecho.nii.gz", volume)      # 4d NIfTI (x, y, z, echo)
tomio.write_volume("multiecho.nrrd", volume)        # NRRD also keeps the metadata (echo times)

vol = tomio.read_volume("multiecho.nrrd")
vol.spacing                                   # -> (1.484375, 1.484375, 2.0)
vol.metadata["EchoTime"]                      # -> '[10, 20]'
vol = tomio.asvolume(volume)                  # in memory: tomio Volume

Export and anonymisation

stack(SeriesNumber=2).export("anonymous")     # copies of the stack's files, anonymized
stack.export("copy", anonymize=False, remove_private_tags=True)

from dicomstack import utils
utils.export_stack("data", "anonymous")       # all DICOM files of a folder
utils.export_file("data/series1/IM1.dcm", "anonymous/IM1.dcm")

Anonymisation empties the patient's name, ID, birth date and sex, the accession number and the referring physician's name (utils.FIELDS_TYPE_2), removes the fields of utils.FIELDS_TYPE_3 (institution, station, study description, patient's age, size, weight, comments, device serial number, protocol name...), and empties all person names. Curves (groups 50xx) are removed. Values can be set with a mapping table (mapper, mapkey: see the command line).

Command line

dicomstack list data                          # List DICOM studies, series and non-DICOM files.
dicomstack values data EchoTime --series 2    # Show unique values of one or more DICOM elements.
dicomstack elements data -i 1 -f SeriesNumber=2   # List the DICOM elements of one frame.
dicomstack export data anonymous              # Export a modified/anonymous copy of a DICOM stack.
dicomstack export -t table.csv -k PatientID data anonymous   # with a mapping table

--cache uses the header cache (see Loading). dicomstack export -h describes the mapping tables (CSV or JSON).

Performance

Indicative numbers (2,000 single-frame MR files, 350 MB, Windows laptop, warm disk cache; benchmarks/):

Loading is dominated by pydicom's header parsing, which threads do not speed up; the default reading threads (workers=8) pay off on slow (network) storage. Uncompressed pixel data is read directly from the files; other transfer syntaxes are decoded by pydicom.

Development

pip install -e ".[test]"
python -m pytest

• test/make_golden.py writes a snapshot of the outputs (values, types, volumes) of the reference version in test/golden/ (not versioned); test/test_golden.py compares the current code with it (skipped without snapshot).
• benchmarks/make_data.py generates the benchmark data (benchmarks/data/, not versioned), benchmarks/bench.py runs the benchmark scenarios.
• Every Python example of this README is run by test/test_readme.py.