API

AbstractDetectorGeometry{T} <: AbstractProjectionGeometry{T}

Geometries of a source (or ray direction) and a flat detector per view; the projection data has size (nu, nv, nviews). Besides rays they describe the detector cells, which footprint methods need.

Interface (in addition to that of AbstractProjectionGeometry):

• nviews(g), detectorsize(g) = (nu, nv)
• subset(g, views): the geometry of a subset of the views

AbstractProjectionGeometry{T}

Supertype of acquisition geometries: where the rays run and the shape of the projection data. T is the element type of the stored coordinates.

Interface:

• datasize(g): size of the projection data, e.g. (nu, nv, nviews)
• ray(g, i) -> (a, b): ray i (linear index into the projection data) as two points SVector{3,T}; callable inside kernels
• raytype(g): RaySegment() if the ray is the segment from a to b, RayLine() if it is the whole line through a and b

Geometries are adapted (moved and converted) with Adapt.adapt, e.g. adapt(CuArray{Float32}, g).

RayLine()

raytype of geometries whose ray i is the whole line through the two points of ray(g, i), e.g. parallel beams. Only the part inside the volume contributes.

RaySegment()

raytype of geometries whose ray i is the segment between the two points of ray(g, i), e.g. from the source to a detector pixel.

datasize(g::AbstractProjectionGeometry)

Size of the projection data.

detectorsize(g::AbstractDetectorGeometry) -> (nu, nv)

nrays(g::AbstractProjectionGeometry)

Number of rays, prod(datasize(g)).

nviews(g::AbstractDetectorGeometry)

projection_geometry(g)

The projection geometry of CTOperators.jl that describes g: g itself for the geometries of this package; package extensions convert the geometries of other packages, e.g. those of Sinograms.jl (SinoPar, SinoFanFlat, SinoFanArc, CtPar, CtFanFlat, CtFanArc). The projectors call it, so they accept these geometries directly.

ray(g::AbstractProjectionGeometry, i) -> (a, b)

Ray i (linear index into the projection data) as two points SVector{3,T}, see raytype. Callable inside kernels.

raytype(g::AbstractProjectionGeometry)

RaySegment() or RayLine().

subset(g::AbstractDetectorGeometry, views)

The geometry of the views views (indices or a range) of g, for ordered-subset methods.

RayList(starts, ends; datasize=(size(starts, 2),))

Arbitrary rays given by their start and end points (columns of the 3 × nrays matrices starts and ends), as in parallelproj. Ray i is the segment from starts[:, i] to ends[:, i]. datasize gives the projection data a shape (column-major order of the rays), e.g. a sinogram.

RayList(g::AbstractProjectionGeometry)

The rays of a geometry with RaySegment rays, computed on the device of g.

ConeBeam(source, centre, u, v; nu, nv)

Cone beam with a flat detector of nu × nv cells and an arbitrary trajectory, given per view k (columns of 3 × nviews matrices):

• source[:, k]: position of the source
• centre[:, k]: centre of the detector
• u[:, k], v[:, k]: vectors from one detector cell to the next along the rows (u) and the columns (v); their lengths are the cell sizes

The centre of cell (iu, iv) is centre + (iu − (nu + 1)/2) u + (iv − (nv + 1)/2) v, and ray (iu, iv, k) runs from the source to this point. The projection data has size (nu, nv, nviews).

circular_cone_beam and helical_cone_beam construct the common trajectories.

ArcConeBeam(source, centre, u, v; nu, nv)

Cone beam with a curved detector of nu × nv cells on a cylinder about the source, as in clinical CT scanners (third-generation geometry), and an arbitrary trajectory, given per view k (columns of 3 × nviews matrices):

• source[:, k]: position of the source, on the axis of the detector cylinder (the axis runs along v)
• centre[:, k]: centre of the detector; its distance from the axis is the radius of the cylinder
• u[:, k]: the direction of the detector rows at the centre (tangent to the arc); its length is the cell size along the arc
• v[:, k]: vector from one detector row to the next along the cylinder axis; its length is the cell size

Cell (iu, iv) sits at the fan angle γ = (iu − (nu + 1)/2) |u| / r from the centre, with r the radius, and at the height (iv − (nv + 1)/2) v from it; ray (iu, iv, k) runs from the source to it. The projection data has size (nu, nv, nviews). u must be orthogonal to v and to centre − source.

circular_cone_beam and helical_cone_beam with detector=:arc construct the common trajectories.

ParallelBeam(direction, centre, u, v; nu, nv)

Parallel beam with a flat detector of nu × nv cells, given per view k (columns of 3 × nviews matrices):

• direction[:, k]: direction of the rays (any nonzero length)
• centre[:, k]: a point of the detector plane, the centre of the detector
• u[:, k], v[:, k]: vectors from one detector cell to the next along the rows (u) and the columns (v); their lengths are the cell sizes

Ray (iu, iv, k) is the whole line through the centre of cell (iu, iv) (see ConeBeam) along direction[:, k] (RayLine). The projection data has size (nu, nv, nviews).

parallel_beam constructs the common rotation about the z axis.

circular_cone_beam(; angles, dso, dsd, nu, nv, du, dv=du, offset_u=0, offset_v=0, source_z=0,
                   detector=:flat)

Cone-beam geometry of a source on a circle of radius dso about the z axis and a detector at distance dsd from the source, nu × nv cells of size du × dv, one view per angle in angles (radians, source counter-clockwise from the y axis).

detector=:flat gives a flat detector (ConeBeam), detector=:arc a detector curved on a cylinder about the source (ArcConeBeam), with du the cell size along the arc at distance dsd (an angle of du / dsd per cell) and offset_u shifting the cells along the arc.

source_z is the height of source and detector centre, a number or one value per view; with source_z varying the trajectory is no longer circular (see helical_cone_beam).

helical_cone_beam(; angles, pitch, z0=0, kwargs...)

ConeBeam geometry of a helical trajectory: as circular_cone_beam (same keywords), with source and detector moving along z by pitch (a length) per turn, z = z0 + pitch (β − β₁)/2π for view angle β (β₁ = angles[1]).

An object that translates along the rotation axis with speed s during a circular scan sees this geometry with pitch = −s · (time per turn).

parallel_beam(; angles, nu, nv, du, dv=du, offset_u=0, offset_v=0)

ParallelBeam geometry rotating about the z axis, nu × nv cells of size du × dv, one view per angle in angles (radians). At angle ϕ the rays run along (−sin ϕ, cos ϕ, 0) and the detector rows along (cos ϕ, sin ϕ, 0): the coordinates (u, v, ϕ) of ImagePhantoms.jl's radon with polar angle 0.

AbstractCTOperator{T}

Supertype of the CT projection operators of CTOperators.jl. T is the real element type of image and projection data. An operator holds its geometry on the device and the workspace of its kernels, created once and reused by every application.

AbstractProjector{T,G<:AbstractProjectionGeometry} <: AbstractCTOperator{T}

Supertype of the projectors of this package: operators between an image on an ImageGeom and projection data of a geometry of type G (stored with its arrays on the device).

Interface:

• geometry(op)::G, image_geometry(op)::ImageGeom
• similar_projector(op, geometry): the same kind of projector, with the same image geometry and options, for another geometry (used by subset and adapt)
• the methods of AbstractCTOperator

AdjointCTOperator(op)

The adjoint (back projection) of op, created by adjoint(op); applying it calls backward!(x, op, y).

colsums(op::AbstractCTOperator)

adjoint(op) * ones: for projectors the sensitivity image, the summed weights of all rays per voxel. A normalisation of SART/SIRT/OS-EM-type methods.

geometry(op::AbstractProjector)

The projection geometry of op, with its arrays on the device of the operator.

image_geometry(op::AbstractProjector)

The ImageGeom the operator was built for.

imagesize(op)

Size of the image the operator takes (2D or 3D, as its ImageGeom).

rowsums(op::AbstractCTOperator)

op * ones: for projectors the length of each ray inside the image. A normalisation of SART/SIRT-type methods.

subset(op::AbstractProjector, views)

The projector of the views views of op (detector geometries), for ordered-subset methods. It shares the image geometry and options with op.

SiddonProjector(image::ImageGeom, geometry; eltype=Float32, storage=Array)

Projector with Siddon's method: the weight of a voxel for a ray is the exact length of the ray inside the voxel (Siddon 1985), computed incrementally as in Jacobs et al. 1998. Works with every AbstractProjectionGeometry. The back projection is the exact adjoint.

eltype is the element type of the computation and of the stored geometry; storage the array type that selects the device (Array, CuArray, ROCArray, …). Image and projection data passed to the operator must live on that device.

Example

image = ImageGeom(; dims=(128, 128, 64), deltas=(1, 1, 1))
geometry = circular_cone_beam(; angles=range(0, 2pi, 361)[1:360], dso=500, dsd=1000,
                              nu=256, nv=128, du=1, dv=1)
op = SiddonProjector(image, geometry)
y = op * x                 # forward projection, size (256, 128, 360)
z = adjoint(op) * y        # back projection

References:

• R. L. Siddon, "Fast calculation of the exact radiological path for a three-dimensional CT array", Med. Phys. 12(2), 252–255, 1985.
• F. Jacobs, E. Sundermann, B. De Sutter, M. Christiaens, I. Lemahieu, "A fast algorithm to calculate the exact radiological path through a pixel or voxel space", J. Comput. Inf. Technol. 6(1), 89–94, 1998.

JosephProjector(image::ImageGeom, geometry; eltype=Float32, storage=Array)

Projector with Joseph's method: the image is interpolated linearly between the voxel centres in every layer of voxels perpendicular to the axis the ray runs most along, and the line integral is the sum over these layers of the interpolated value times the length of the ray in the layer (Joseph 1982). Compared with SiddonProjector the image model is smooth: the projection changes continuously with the rays, the rays of a view cover the voxels more evenly (except near 45°, where both methods coincide), and the projection is exact for images that are linear functions of the position rather than for voxel images. Works with every AbstractProjectionGeometry. The back projection is the exact adjoint.

Outside the volume the image is zero; across the ray it fades to zero within one voxel outside the outermost voxel centres. Segments that start or end inside the volume take the part of the first and last layer they cover.

eltype is the element type of the computation and of the stored geometry; storage the array type that selects the device (Array, CuArray, ROCArray, …). Image and projection data passed to the operator must live on that device.

Example

image = ImageGeom(; dims=(128, 128, 64), deltas=(1, 1, 1))
geometry = circular_cone_beam(; angles=range(0, 2pi, 361)[1:360], dso=500, dsd=1000,
                              nu=256, nv=128, du=1, dv=1)
op = JosephProjector(image, geometry)
y = op * x                 # forward projection, size (256, 128, 360)
z = adjoint(op) * y        # back projection

Reference:

• P. M. Joseph, "An improved algorithm for reprojecting rays through pixel images", IEEE Trans. Med. Imaging 1(3), 192–196, 1982.

VoxelGrid{T}(ig::ImageGeom)

The voxel grid of an ImageGeom in the form the kernels use: three dimensions (a 2D ImageGeom becomes a volume with nz = 1 and voxel size 1 along z, centred at z = 0), element type T, no units, no mask. first is the world coordinate of the centre of voxel (1, 1, 1); voxel (i, j, k) has its centre at first + ((i, j, k) .- 1) .* deltas. Negative deltas (flipped axes, allowed by ImageGeoms) are kept.

Units of the ImageGeom are stripped; they must be the units of the projection geometry.

backward!(x, op::AbstractCTOperator, y)

Back projection (the adjoint) of the projection data y into the image x, overwriting x. Implemented by every operator; users call mul!(x, adjoint(op), y) or adjoint(op) * y.

bounds(g::VoxelGrid) -> (lo, hi)

Corners of the box covered by the voxels (outer faces, not centres), lo .<= hi.

clip(p, d, dims, ::RaySegment | ::RayLine) -> (amin, amax)

Parameters of the part of the ray p + α d (voxel coordinates) inside the volume of size dims: the box [1/2, n + 1/2] per dimension, intersected with α ∈ [0, 1] for segments. The ray misses the volume if amin >= amax.

clip(p, d, lo, hi, ::RaySegment | ::RayLine) -> (amin, amax)

Parameters of the part of the ray p + α d inside the box [lo, hi], intersected with α ∈ [0, 1] for segments.

forward!(y, op::AbstractCTOperator, x)

Forward projection of the image x (of size imagesize(op)) into y (of size datasize(op)), overwriting y. Implemented by every operator; users call mul!(y, op, x) or op * x.

similar_projector(op::AbstractProjector, geometry)

A projector of the same kind and with the same image geometry and options as op, for geometry.

to_index(g::VoxelGrid, p)

Continuous voxel coordinates of the world point p: voxel (i, j, k) covers [i − 1/2, i + 1/2) × …, its centre has the integer coordinates (i, j, k). In these coordinates the voxels are unit cubes with positive orientation, whatever the sign of deltas.