RasmahRasmah

DICOM: Scalar-Volume Imaging

Overview

DICOM (Digital Imaging and Communications in Medicine) is the medical-imaging exchange format: a series of 2D grayscale slices — a CT or MRI scan — that together define a 3D scalar field (e.g. Hounsfield density). It is the one format in Rasmah that is not a mesh or a boundary representation, so it enters through the implicit/SDF path: read a scan into a scalar volume, threshold it into a shape, and mesh it.

Rasmah reads and writes DICOM natively (no external library) with four entry points — read_dicom / read_dicom_series and write_dicom / write_dicom_series — around the DicomVolume type.

The scalar-volume representation

A DicomVolume is a regular 3D grid of values plus the voxel geometry: values is a rows × columns × slices array indexed [row, column, slice], spacing is the [x, y, z] voxel size in millimetres, and origin is the patient position of the first voxel’s centre.

The volume is callable: it trilinearly interpolates the field at any point, and it has a bounds:

using Rasmah

# a tiny 3×4×2 volume, unit spacing, first voxel centred at the origin
values = zeros(3, 4, 2)
for k in 1:2, j in 1:4, i in 1:3
    values[i, j, k] = i + 10j + 100k
end
vol = DicomVolume(values, [1.0, 1.0, 1.0], [0.0, 0.0, 0.0])

println(vol([0.0, 0.0, 0.0]))          # exactly voxel (1,1,1)
println(vol([0.5, 0.0, 0.0]))          # halfway to the next column
println(bounds(vol))
111.0
111.5
([-0.5, -0.5, -0.5], [2.5, 3.5, 1.5])

Building a synthetic volume

Real scans are large; for this tutorial we synthesise one. A CT image stores Hounsfield units — roughly +1000 for bone and −1000 for air — so we model a spherical “bone” of density 1000 inside air of −1000. sample_grid samples any callable field onto a regular grid and returns a DicomVolume:

density(x) = field(sphere(6.0), x) < 0 ? 1000.0 : -1000.0

vol = sample_grid(density, [-12.0, -12.0, -12.0], [12.0, 12.0, 12.0], (48, 48, 48))

println(size(vol.values))
println("spacing: ", round.(vol.spacing; digits=3))
println("centre density: ", vol([0.0, 0.0, 0.0]))
println("corner density: ", vol([12.0, 12.0, 12.0]))
(48, 48, 48)
spacing: [0.511, 0.511, 0.511]
centre density: 1000.0
corner density: -1000.0

Writing DICOM

write_dicom serialises the volume as a single Explicit-VR-Little-Endian file (MONOCHROME2, signed 16-bit pixels). Because DICOM stores integer pixels, the values are rounded to the nearest integer — exact here, since our densities are already whole numbers:

mktempdir() do dir
    path = joinpath(dir, "bone.dcm")
    write_dicom(vol, path)
    println("wrote ", filesize(path), " bytes")
end
wrote 221982 bytes

Reading DICOM back

read_dicom parses the file (single- or multi-frame) and reconstructs the DicomVolume, recovering the values, the voxel spacing (PixelSpacing + SliceThickness), and the origin (ImagePositionPatient):

mktempdir() do dir
    path = joinpath(dir, "bone.dcm")
    write_dicom(vol, path)

    vol2 = read_dicom(path)
    println("values identical:  ", vol2.values == vol.values)
    println("spacing identical: ", vol2.spacing ≈ vol.spacing)
    println("origin identical:  ", vol2.origin ≈ vol.origin)
end
values identical:  true
spacing identical: true
origin identical:  true

Segmenting and meshing

segment(vol, threshold) turns the volume into a level-set ImplicitGeometry, φ(x) = threshold − vol(x), whose zero set is the vol = threshold isosurface. Because it is an ImplicitGeometry, the whole implicit toolkit applies — bounds, surface_mesh, marching_tetrahedra, gradient. Segmenting at 0 picks out the bone/air boundary and meshes it:

seg = segment(vol, 0.0)

println("inside the bone:  ", seg([0.0, 0.0, 0.0]))    # negative
println("in the air:       ", seg([11.0, 0.0, 0.0]))   # positive

m = surface_mesh(seg; n=64)
println(size(m.vertices, 2), " vertices, ", size(m.faces, 2), " triangles")

vtk(m)
inside the bone:  -1000.0
in the air:       1000.0
3822 vertices, 7640 triangles

The slice-series round-trip

write_dicom_series writes one single-frame file per slice (the classic CT layout), and read_dicom_series reassembles them — ordered by patient position, with the between-slice spacing recovered from the median slice distance:

mktempdir() do dir
    series = joinpath(dir, "series")
    write_dicom_series(vol, series)

    files = sort(readdir(series))
    println(length(files), " slice files: ", files[1], " … ", files[end])

    vol3 = read_dicom_series(series)
    println("series values identical:  ", vol3.values == vol.values)
    println("series spacing identical: ", vol3.spacing ≈ vol.spacing)
end
48 slice files: slice_0001.dcm … slice_0048.dcm
series values identical:  true
series spacing identical: true

NIfTI (the research-side counterpart)

NIfTI-1 (.nii) is the neuroimaging counterpart to DICOM: a single-file scalar volume (a 352-byte header followed by the voxel data). Rasmah reads and writes the axis-aligned subset, mapping a NIfTI volume onto the same DicomVolume representation — and, unlike DICOM’s integer pixels, it preserves floating-point values exactly:

mktempdir() do dir
    p = joinpath(dir, "bone.nii")
    write_nifti(vol, p)
    println("wrote ", filesize(p), " bytes")

    vol4 = read_nifti(p)
    println("values identical:  ", vol4.values == vol.values)
    println("spacing identical: ", vol4.spacing ≈ vol.spacing)
    println("origin identical:  ", vol4.origin ≈ vol.origin)
end
wrote 885088 bytes
values identical:  true
spacing identical: true
origin identical:  true

Next steps