From b830398b27ac2c14957bdb67f6924148b065305a Mon Sep 17 00:00:00 2001 From: Niklas Hackelberg Date: Mon, 2 Nov 2026 18:05:19 +0100 Subject: [PATCH 1/6] Benchmark: GPU selection with --device; guide: choose the GPU bench.jl --device N runs the GPU libraries on the N-th GPU (KernelAbstractions' device! for CTOperators, CuPy for parallelproj), e.g. a free one on a shared machine. The guide "Run on a GPU" explains selecting the GPU before creating a projector. Co-Authored-By: Claude Opus 5.5 --- benchmark/README.md | 3 ++- benchmark/bench.jl | 7 ++++--- benchmark/worker/adapters/ctoperators-cuda.jl | 3 +++ benchmark/worker/adapters/ctoperators-rocm.jl | 3 +++ benchmark/worker/adapters/parallelproj-cuda.jl | 3 ++- benchmark/worker/worker.jl | 2 ++ docs/literate/howto/gpu.jl | 13 +++++++++++++ 7 files changed, 29 insertions(+), 5 deletions(-) diff --git a/benchmark/README.md b/benchmark/README.md index 37397a0..cd950de 100644 --- a/benchmark/README.md +++ b/benchmark/README.md @@ -19,7 +19,8 @@ GPU libraries (`*-rocm`: AMD GPUs through AMDGPU.jl; `*-cuda`: NVIDIA GPUs throu `julia --project=benchmark/envs/parallelproj-cuda benchmark/check_parallelproj_cuda.jl` checks that parallelproj computes on CuPy arrays directly (otherwise its timings would include transfers). -Options: `--libs` (comma separated), `--threads` (comma separated, default all cores), `--budget` +Options: `--libs` (comma separated), `--threads` (comma separated, default all cores), `--device` +(the GPU of GPU libraries, counted from 1, default 1; e.g. a free one on a shared machine), `--budget` (seconds of repeated calls per measurement, default 2). Measure on a quiet machine. ## How it works diff --git a/benchmark/bench.jl b/benchmark/bench.jl index f0cec6b..283163f 100644 --- a/benchmark/bench.jl +++ b/benchmark/bench.jl @@ -38,7 +38,7 @@ function instantiate(lib) run(`$(Base.julia_cmd()) --project=$env -e "using Pkg; Pkg.instantiate()"`) end -function run_suite(suite, libs, threads, budget) +function run_suite(suite, libs, threads, budget, device) haskey(SUITES, suite) || error("unknown suite $suite; see `julia benchmark/bench.jl list`") mkpath(joinpath(ROOT, "out")) out = joinpath(ROOT, "out", Dates.format(now(), "yyyymmdd-HHMMSS") * "-$suite.csv") @@ -48,7 +48,8 @@ function run_suite(suite, libs, threads, budget) for nt in threads env = environment(lib) cmd = `$(Base.julia_cmd()) --project=$env --threads=$nt $worker $lib $out $budget $(SUITES[suite])` - run(addenv(cmd, "OMP_NUM_THREADS" => string(nt))) # OpenMP libraries (parallelproj) + # OpenMP libraries (parallelproj); the GPU of the GPU libraries + run(addenv(cmd, "OMP_NUM_THREADS" => string(nt), "BENCH_DEVICE" => device)) end end println("results: ", out) @@ -63,7 +64,7 @@ function main(args) elseif args[1] == "run" libs = split(option(args, "libs", join(LIBS, ",")), ",") threads = parse.(Int, split(option(args, "threads", string(Sys.CPU_THREADS)), ",")) - run_suite(args[2], libs, threads, option(args, "budget", "2")) + run_suite(args[2], libs, threads, option(args, "budget", "2"), option(args, "device", "1")) else error("unknown command $(args[1])") end diff --git a/benchmark/worker/adapters/ctoperators-cuda.jl b/benchmark/worker/adapters/ctoperators-cuda.jl index 9c6f0b4..b03e4dd 100644 --- a/benchmark/worker/adapters/ctoperators-cuda.jl +++ b/benchmark/worker/adapters/ctoperators-cuda.jl @@ -4,6 +4,9 @@ using CUDA using LinearAlgebra: mul! +# the GPU to run on, before any array is created on it +CTOperators.KA.device!(CUDABackend(), DEVICE) + library_version() = string(pkgversion(CTOperators)) method_name() = string(nameof(PROJECTOR))[1:end-length("Projector")] device_name() = CUDA.name(CUDA.device()) diff --git a/benchmark/worker/adapters/ctoperators-rocm.jl b/benchmark/worker/adapters/ctoperators-rocm.jl index bb8d20c..6982075 100644 --- a/benchmark/worker/adapters/ctoperators-rocm.jl +++ b/benchmark/worker/adapters/ctoperators-rocm.jl @@ -4,6 +4,9 @@ using AMDGPU using LinearAlgebra: mul! +# the GPU to run on, before any array is created on it +CTOperators.KA.device!(ROCBackend(), DEVICE) + library_version() = string(pkgversion(CTOperators)) method_name() = string(nameof(PROJECTOR))[1:end-length("Projector")] device_name() = AMDGPU.HIP.name(AMDGPU.device()) diff --git a/benchmark/worker/adapters/parallelproj-cuda.jl b/benchmark/worker/adapters/parallelproj-cuda.jl index 2030332..17847b6 100644 --- a/benchmark/worker/adapters/parallelproj-cuda.jl +++ b/benchmark/worker/adapters/parallelproj-cuda.jl @@ -11,10 +11,11 @@ using LinearAlgebra: normalize const np = pyimport("numpy") const cp = pyimport("cupy") const ppc = pyimport("parallelproj_core") +cp.cuda.Device(DEVICE - 1).use() # the GPU to run on (CuPy counts from 0) library_version() = pyconvert(String, ppc.__version__) method_name() = "Joseph" -device_name() = pyconvert(String, cp.cuda.runtime.getDeviceProperties(0)["name"].decode()) +device_name() = pyconvert(String, cp.cuda.runtime.getDeviceProperties(DEVICE - 1)["name"].decode()) supports(p) = true measures_error(p) = true diff --git a/benchmark/worker/worker.jl b/benchmark/worker/worker.jl index 6762970..ce859d2 100644 --- a/benchmark/worker/worker.jl +++ b/benchmark/worker/worker.jl @@ -7,6 +7,8 @@ using BenchmarkTools, Printf const LIB, OUT, BUDGET = ARGS[1], ARGS[2], parse(Float64, ARGS[3]) +# the GPU of GPU libraries, counted from 1 (bench.jl --device) +const DEVICE = parse(Int, get(ENV, "BENCH_DEVICE", "1")) const PROBLEM_NAMES = ARGS[4:end] # Defaults for CPU libraries; adapters for devices override them: the device in the results, moving diff --git a/docs/literate/howto/gpu.jl b/docs/literate/howto/gpu.jl index 040c79b..d539dc8 100644 --- a/docs/literate/howto/gpu.jl +++ b/docs/literate/howto/gpu.jl @@ -22,6 +22,19 @@ # For AMD GPUs use `ROCArray` (AMDGPU.jl), for Apple GPUs `MtlArray` (Metal.jl); Metal has no # `Float64`. The test suite runs on CPUs and on NVIDIA and AMD GPUs. # +# ## Choose the GPU +# +# On a machine with several GPUs, select one with `device!` of KernelAbstractions.jl (counted from 1) +# *before* creating the projector: the projector stores its geometry on the GPU that is active when +# it is created, and all data passed to it must live there. +# +# ```julia +# using KernelAbstractions: device!, ndevices +# ndevices(CUDABackend()) # number of GPUs +# device!(CUDABackend(), 2) # the second GPU, for this task +# op = JosephProjector(image, geometry; storage=CuArray) +# ``` + # ## Apply it to data on the device # # Image and projection data must live on the same device as the projector: From fb2c8db446721490e4df8d196b7879f7ea07e68d Mon Sep 17 00:00:00 2001 From: Niklas Hackelberg Date: Mon, 2 Nov 2026 18:10:39 +0100 Subject: [PATCH 2/6] Benchmark: script to tune the launch configuration on a GPU tune_launch.jl times a projector on one GPU (CUDA or ROCm, device selectable) for workgroup sizes and tiles of detector cells, the device's defaults first, and lists the fastest configurations. Co-Authored-By: Claude Opus 5.5 --- benchmark/README.md | 5 ++- benchmark/tune_launch.jl | 66 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 70 insertions(+), 1 deletion(-) create mode 100644 benchmark/tune_launch.jl diff --git a/benchmark/README.md b/benchmark/README.md index cd950de..f8e73fa 100644 --- a/benchmark/README.md +++ b/benchmark/README.md @@ -17,7 +17,10 @@ GPU libraries (`*-rocm`: AMD GPUs through AMDGPU.jl; `*-cuda`: NVIDIA GPUs throu `parallelproj-cuda`, the CUDA build of parallelproj with CuPy arrays in its own environment `envs/parallelproj-cuda`) run only when named with `--libs`. Before a first CUDA run, `julia --project=benchmark/envs/parallelproj-cuda benchmark/check_parallelproj_cuda.jl` checks that -parallelproj computes on CuPy arrays directly (otherwise its timings would include transfers). +parallelproj computes on CuPy arrays directly (otherwise its timings would include transfers). To choose +a `LaunchConfig` for a GPU, `julia --project=benchmark/envs/ctoperators benchmark/tune_launch.jl + [device] [problem] [projector]` times workgroup sizes and tiles of detector cells on one +problem. Options: `--libs` (comma separated), `--threads` (comma separated, default all cores), `--device` (the GPU of GPU libraries, counted from 1, default 1; e.g. a free one on a shared machine), `--budget` diff --git a/benchmark/tune_launch.jl b/benchmark/tune_launch.jl new file mode 100644 index 0000000..42398b6 --- /dev/null +++ b/benchmark/tune_launch.jl @@ -0,0 +1,66 @@ +# Times the launch configurations of a projector on one GPU, to choose `LaunchConfig` for a device: +# julia --project=benchmark/envs/ctoperators benchmark/tune_launch.jl [device] [problem] [projector] +# device: the GPU, counted from 1 (default 1); problem: one of benchmark/src/problems.jl (default +# cone256); projector: Siddon or Joseph (default Joseph). Prints forward and back projection times +# (minimum over repeated calls) for workgroup sizes and tiles of detector cells, the device's +# defaults first, and the fastest configurations last. + +using CTOperators, BenchmarkTools, Printf +using LinearAlgebra: mul! +include(joinpath(@__DIR__, "src", "problems.jl")) + +const BACKEND_NAME = get(ARGS, 1, "cuda") +const DEVICE = parse(Int, get(ARGS, 2, "1")) +const PROBLEM = get(ARGS, 3, "cone256") +const PROJECTOR = Dict("Siddon" => SiddonProjector, "Joseph" => JosephProjector)[get(ARGS, 4, "Joseph")] + +if BACKEND_NAME == "cuda" + using CUDA + const ArrayType, backend, sync = CuArray, CUDABackend(), CUDA.synchronize + const devicename() = CUDA.name(CUDA.device()) +elseif BACKEND_NAME == "rocm" + using AMDGPU + const ArrayType, backend, sync = ROCArray, ROCBackend(), AMDGPU.synchronize + const devicename() = AMDGPU.HIP.name(AMDGPU.device()) +else + error("backend must be cuda or rocm, got $BACKEND_NAME") +end +CTOperators.KA.device!(backend, DEVICE) + +p = PROBLEMS[PROBLEM]() +x = ArrayType(Float32.(phantom_image(p))) +y = ArrayType(rand(Float32, datasize(p.geometry)...)) +xb, yf = similar(x), similar(y) + +# minimum time in ms of f, which launches kernels, including the wait for the device +time_ms(f) = (f(); sync(); minimum(@benchmark ($f(); $sync()) seconds = 2 evals = 1).time / 1e6) + +function measure(launch) + op = PROJECTOR(p.image, p.geometry; storage=ArrayType, launch) + time_ms(() -> mul!(yf, op, x)), time_ms(() -> mul!(xb, adjoint(op), y)) +end + +configs = Any[nothing] # the defaults of the device +for workgroupsize in (nothing, 32, 64, 128, 256), tile in (nothing, (8, 4), (8, 8), (16, 4), (16, 8), (16, 16)) + push!(configs, LaunchConfig(; workgroupsize, tile)) +end + +label(::Nothing) = "defaults of the device" +label(c::LaunchConfig) = "workgroupsize=$(something(c.workgroupsize, "backend")), tile=$(something(c.tile, "none"))" + +@printf("%s on %s (device %d), %s\n\n", nameof(PROJECTOR), devicename(), DEVICE, PROBLEM) +results = [] +for c in configs + tf, tb = measure(c) + push!(results, (c, tf, tb)) + @printf("%-44s forward %8.2f ms back %8.2f ms\n", label(c), tf, tb) +end + +println("\nfastest forward:") +for (c, tf, tb) in sort(results; by=r -> r[2])[1:3] + @printf(" %-44s forward %8.2f ms back %8.2f ms\n", label(c), tf, tb) +end +println("fastest back:") +for (c, tf, tb) in sort(results; by=r -> r[3])[1:3] + @printf(" %-44s forward %8.2f ms back %8.2f ms\n", label(c), tf, tb) +end From b694b6bbbac0ebf1f43eea91bbf5936789be9e3e Mon Sep 17 00:00:00 2001 From: Niklas Hackelberg Date: Mon, 2 Nov 2026 19:35:42 +0100 Subject: [PATCH 3/6] Docs: explanation of the parameters of a scan, step by step A new explanation page builds a scan up from one ray to a helical cone-beam scan and introduces the parameters of the trajectory functions where they appear (angles, dso, dsd, nu, nv, du, dv, pitch, z0), with a figure per step. Getting Started lists the parameters and links to it. Co-Authored-By: Claude Opus 5.5 --- docs/literate/examples/getting_started.jl | 27 ++- docs/literate/explanation/conventions.jl | 3 + docs/literate/explanation/scan_parameters.jl | 101 +++++++++++ docs/make.jl | 1 + docs/src/index.md | 3 +- docs/utils/scan_steps.jl | 174 +++++++++++++++++++ 6 files changed, 301 insertions(+), 8 deletions(-) create mode 100644 docs/literate/explanation/scan_parameters.jl create mode 100644 docs/utils/scan_steps.jl diff --git a/docs/literate/examples/getting_started.jl b/docs/literate/examples/getting_started.jl index 4ae5ac4..3bd65d9 100644 --- a/docs/literate/examples/getting_started.jl +++ b/docs/literate/examples/getting_started.jl @@ -39,13 +39,26 @@ x_true = Float32.(phantom(axis(image, 1), axis(image, 2), axis(image, 3), object # and detector rotate about the object, and in our scan the object also moves along the rotation # axis. # -# The geometry describes this: per view the position of the source and of the detector. The source -# is 128 mm from the rotation axis and 256 mm from the detector, so the object is magnified by 2 on -# the detector. The detector has 192 × 96 cells of 1 mm; at the rotation axis a cell covers 0.5 mm, -# so two rays pass through each voxel per direction. (With one ray per voxel, the rays of a view -# cover the voxels unevenly; see [Siddon's Method](@ref).) The object moves by 32 mm, half the -# volume height, per turn. In the object's frame, source and detector therefore move on a helix: two -# turns, 240 views. +# The geometry describes this: per view the position of the source and of the detector. The +# functions that build the common scans take these parameters: +# +# - `angles`: the angle β of the source in every view, counterclockwise from the `y` axis (radians); +# - `dso`: the distance from the source to the rotation axis, the `z` axis; +# - `dsd`: the distance from the source to the detector; +# - `nu`, `nv`: the number of detector cells along `u` (a detector row) and along `v` (a column, +# parallel to the rotation axis), `du`, `dv` their sizes; +# - for helical scans, `pitch`: how far source and detector move along `z` per turn, and `z0`: the +# height of the first view. +# +# All lengths are in the unit of the image (mm here). [The Parameters of a Scan](@ref) introduces +# them step by step, with figures. + +# In our scan the source is 128 mm from the rotation axis and 256 mm from the detector, so the +# object is magnified by 2 on the detector. The detector has 192 × 96 cells of 1 mm; at the rotation +# axis a cell covers 0.5 mm, so two rays pass through each voxel per direction. (With one ray per +# voxel, the rays of a view cover the voxels unevenly; see [Siddon's Method](@ref).) The object moves +# by 32 mm, half the volume height, per turn. In the object's frame, source and detector therefore +# move on a helix: two turns, 240 views, starting 32 mm below the center. geometry = helical_cone_beam(; angles=range(0, 4pi, 240), pitch=n / 2, z0=-n / 2, dso=2n, dsd=4n, nu=192, nv=96, du=1.0, dv=1.0); diff --git a/docs/literate/explanation/conventions.jl b/docs/literate/explanation/conventions.jl index 3b5f695..5f79b39 100644 --- a/docs/literate/explanation/conventions.jl +++ b/docs/literate/explanation/conventions.jl @@ -32,6 +32,9 @@ include(joinpath(pkgdir(CTOperators), "docs", "utils", "convention_figures.jl")) # # ## Cone beam and its trajectories # +# [The Parameters of a Scan](@ref) introduces the parameters of the trajectory functions (`dso`, +# `dsd`, `nu`, `nv`, `du`, `dv`, `pitch`, `z0`) step by step; this section defines them exactly. +# # A view is given by the position of the source, the center of the flat detector and two vectors # `u` and `v` from one detector cell to the next, along the rows and the columns of the detector # ([`ConeBeam`](@ref)). The ray of cell `(iu, iv)` runs from the source to the center of the cell. diff --git a/docs/literate/explanation/scan_parameters.jl b/docs/literate/explanation/scan_parameters.jl new file mode 100644 index 0000000..04be746 --- /dev/null +++ b/docs/literate/explanation/scan_parameters.jl @@ -0,0 +1,101 @@ +# # The Parameters of a Scan +# +# The functions that build the common scans, such as [`circular_cone_beam`](@ref) and +# [`helical_cone_beam`](@ref), take a handful of parameters: `angles`, `dso`, `dsd`, `nu`, `nv`, +# `du`, `dv`, and for helical scans `pitch` and `z0`. This page builds a scan up step by step, from +# a single ray to a helical cone-beam scan, and introduces each parameter where it first appears. + +using CTOperators, CairoMakie +include(joinpath(pkgdir(CTOperators), "docs", "utils", "scan_steps.jl")) #hide +include(joinpath(pkgdir(CTOperators), "docs", "utils", "tutorial_figures.jl")) #hide + +# ## One ray +# +# An X-ray source sends a ray through the object to a detector cell. The cell measures how much the +# ray was attenuated; after taking the logarithm, this is the line integral of the attenuation +# coefficient along the ray: the sum of the image values along the ray, weighted by the length of +# the ray in each voxel. Computing these line integrals for all rays is the forward projection. + +fig = step_ray_figure() #hide +save("steps_ray.png", fig) #hide +nothing #hide +# ![One ray](steps_ray.png) + +# ## A detector row +# +# A detector has many cells side by side: `nu` cells of size `du` along the direction `u`. Every +# cell gets its own ray from the source, so a single position of source and detector measures a +# fan of line integrals. `dsd` is the distance from the source to the detector. + +fig = step_row_figure() #hide +save("steps_row.png", fig) #hide +nothing #hide +# ![A detector row](steps_row.png) + +# ## Views +# +# Source and detector rotate together about the rotation axis, the `z` axis through the origin. The +# source moves on a circle of radius `dso`, the distance from the source to the axis. Each position +# is a view; its angle β is measured counterclockwise from the `y` axis. `angles` lists the angles of +# all views, in radians. + +fig = step_rotation_figure() #hide +save("steps_rotation.png", fig) #hide +nothing #hide +# ![Views](steps_rotation.png) + +# ## Cone beam +# +# A cone-beam detector has `nv` rows of cells along the direction `v`, parallel to the rotation axis, +# each of size `dv`. The rays from the source to all cells form a cone, and a single view measures a +# 2D image of line integrals, `nu × nv` values. The cells are drawn as a checkerboard: + +fig = step_cone_figure() #hide +save("steps_cone.png", fig) #hide +nothing #hide +# ![Cone beam: the detector cells as a checkerboard](steps_cone.png) + +# Together these are the parameters of a circular cone-beam scan. With 64 views over one turn: + +image = ImageGeom(; dims=(64, 64, 64), deltas=(1.0, 1.0, 1.0)) +circle = circular_cone_beam(; angles=range(0, 2pi, 65)[1:64], dso=128, dsd=256, nu=96, nv=64, du=2.0, dv=2.0) +datasize(circle) + +# The data has one detector image per view, `(nu, nv, number of views)`. View 1 in 3D: the source +# (orange), the cone of rays (blue edges, a few rays orange) through the volume (gray box) to the +# detector; the orange circle is the path of the source. + +fig = cone_scan_figure(circle, image; title="circular cone beam, view 1", height=false) #hide +save("steps_circle.png", fig) #hide +nothing #hide +# ![View 1 of the circular scan](steps_circle.png) + +# ## Helical scan +# +# In a helical scan source and detector also move along the rotation axis: by `pitch` per turn, +# starting at the height `z0` in the first view. Their height therefore rises linearly with the +# angle. The same geometry describes an object that moves along the axis during a circular scan, +# for example on a table. + +fig = step_helix_figure() #hide +save("steps_helix.png", fig) #hide +nothing #hide +# ![Height of the source in a helical scan](steps_helix.png) + +# A helical scan takes the parameters of the circular one plus `pitch` and `z0`; here two turns, 24 +# mm per turn, starting 24 mm below the center: + +helix = helical_cone_beam(; angles=range(0, 4pi, 128), pitch=24, z0=-24, dso=128, dsd=256, + nu=96, nv=64, du=2.0, dv=2.0) + +fig = cone_scan_figure(helix, image; title="helical cone beam, view 1") #hide +save("steps_helical.png", fig) #hide +nothing #hide +# ![View 1 of the helical scan and the height of the source](steps_helical.png) + +# ## More +# +# All lengths are in the unit of the image (mm here). [Coordinates and Conventions](@ref) defines +# the parameters exactly, including curved detectors (`detector=:arc`) and detector offsets +# (`offset_u`, `offset_v`), and [Set Up a Geometry](@ref) shows parallel-beam scans, arbitrary +# trajectories and ray lists. diff --git a/docs/make.jl b/docs/make.jl index 6fc4333..e7087bb 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -45,6 +45,7 @@ makedocs(; "Implement a CT Operator" => "generated/howto/implement_operator.md", ], "Explanation" => [ + "The Parameters of a Scan" => "generated/explanation/scan_parameters.md", "Coordinates and Conventions" => "generated/explanation/conventions.md", "Siddon's Method" => "generated/explanation/siddon.md", "Joseph's Method" => "generated/explanation/joseph.md", diff --git a/docs/src/index.md b/docs/src/index.md index 0474e53..5632d7a 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -54,7 +54,8 @@ y = op * rand(Float32, 128, 128, 128) # projections, size (256, 128, 360) [differentiate through a projector](@ref "Differentiate through a Projector"), how to [use geometries of Sinograms.jl](@ref "Use Geometries of Sinograms.jl") and how to [implement a CT operator](@ref "Implement a CT Operator"). -- The explanation pages describe the [coordinates and conventions](@ref "Coordinates and Conventions"), +- The explanation pages introduce [the parameters of a scan](@ref "The Parameters of a Scan") step by step, + describe the [coordinates and conventions](@ref "Coordinates and Conventions"), the projection methods, [Siddon's Method](@ref) and [Joseph's Method](@ref), and [workgroups and ray order on GPUs](@ref "Workgroups and Ray Order on GPUs"). - The [API](@ref) lists all types and functions. diff --git a/docs/utils/scan_steps.jl b/docs/utils/scan_steps.jl new file mode 100644 index 0000000..6aa0bd2 --- /dev/null +++ b/docs/utils/scan_steps.jl @@ -0,0 +1,174 @@ +# Figures of the page "The Parameters of a Scan" (CairoMakie): a scan built up step by step, from one +# ray to a helical cone-beam scan. Drawn from geometries of the package, so they follow its +# conventions. Included invisibly by the page. + +using CTOperators, CairoMakie +using LinearAlgebra: normalize, dot + +const STEP_COLORS = (source=:darkorange, ray=(:darkorange, 0.6), cells=(:steelblue, :lightsteelblue), + image=(:gray60, 0.25), frame=:gray50) + +# a dimension line from a to b, shifted by `shift`, with ticks at its ends and `label` in the middle +# (rotated along the line, never upside down), on the side of the shift +function dimension!(ax, a, b, shift, label; labelshift=0.2, fontsize=15) + p, q = Point2f(a + shift), Point2f(b + shift) + d = normalize(q - p) + n = Point2f(-d[2], d[1]) + lines!(ax, [p, q]; color=:black, linewidth=1) + for e in (p, q) + lines!(ax, [e - 0.08n, e + 0.08n]; color=:black, linewidth=1) + end + angle = atan(d[2], d[1]) + angle = angle > pi / 2 ? angle - pi : angle <= -pi / 2 ? angle + pi : angle + side = dot(n, shift) >= 0 ? 1 : -1 + text!(ax, ((p + q) / 2 + side * labelshift * n)...; text=label, align=(:center, :center), rotation=angle, fontsize) +end + +# the cells of a detector row between the edge points `edges`, alternately colored +function cells!(ax, edges) + for k in 1:length(edges)-1 + lines!(ax, [Point2f(edges[k]), Point2f(edges[k+1])]; color=STEP_COLORS.cells[isodd(k) ? 1 : 2], linewidth=9) + end +end + +function step_axis(fig, pos; title="", limits=(-3.4, 3.4, -3.2, 3.6)) + ax = Axis(fig[pos...]; aspect=DataAspect(), title, limits) + hidedecorations!(ax); hidespines!(ax) + ax +end + +image_square!(ax; h=0.9) = poly!(ax, Rect2f(-h, -h, 2h, 2h); color=STEP_COLORS.image, strokecolor=STEP_COLORS.frame, strokewidth=1) + +# view 1 of a fan of nu cells (the central detector row of a cone beam) from above: source, the edges +# of the cells and the center of the detector, in x, y +function topview_geometry(; dso=3.0, dsd=5.0, nu=6, du=0.55, beta=0.0) + g = circular_cone_beam(; angles=[beta], dso, dsd, nu, nv=1, du) + s, c, u = g.source[1:2, 1], g.center[1:2, 1], g.u[1:2, 1] + s, c, u, [c + (k - nu / 2) * u for k in 0:nu] +end + +# 1. one ray from the source to one detector cell through the image +function step_ray_figure() + fig = Figure(size=(520, 430)) + ax = step_axis(fig, (1, 1); title="one ray") + s, c, u, edges = topview_geometry(; nu=1, du=0.55) + image_square!(ax) + cells!(ax, edges) + lines!(ax, [Point2f(s), Point2f(c)]; color=STEP_COLORS.ray, linewidth=2) + scatter!(ax, [Point2f(s)]; color=STEP_COLORS.source, markersize=16) + text!(ax, s...; text="source", align=(:left, :center), offset=(12, 0), fontsize=15) + text!(ax, edges[end]...; text="detector cell", align=(:left, :center), offset=(10, 0), fontsize=15) + text!(ax, 0, 0; text="image", align=(:right, :bottom), offset=(-6, 4), fontsize=14, color=:gray30) + text!(ax, 0.12, 1.6; text="ray: the cell measures the\nline integral of the image\nalong it", align=(:left, :center), fontsize=13) + fig +end + +# 2. a detector row: nu cells of size du along u, the rays to their centers, the distance dsd +function step_row_figure(; nu=6, du=0.55) + fig = Figure(size=(520, 430)) + ax = step_axis(fig, (1, 1); title="a detector row") + s, c, u, edges = topview_geometry(; nu, du) + image_square!(ax) + for k in 1:nu + lines!(ax, [Point2f(s), Point2f((edges[k] + edges[k+1]) / 2)]; color=STEP_COLORS.ray, linewidth=1) + end + cells!(ax, edges) + scatter!(ax, [Point2f(s)]; color=STEP_COLORS.source, markersize=16) + text!(ax, s...; text="source", align=(:left, :center), offset=(12, 0), fontsize=15) + ray = normalize(c - s) + dimension!(ax, edges[1], edges[2], -0.3 * ray, "du") + arrows2d!(ax, [Point2f(edges[end] + 0.1 * normalize(u))], [Vec2f(0.7 * normalize(u))]; color=STEP_COLORS.cells[1]) + text!(ax, (edges[end] + 0.85 * normalize(u))...; text="u", align=(:left, :center), offset=(4, 0), fontsize=15) + text!(ax, (c - 0.75 * ray)...; text="nu cells", align=(:center, :center), fontsize=15) + dimension!(ax, s, c, [-2.3, 0.0], "dsd") + fig +end + +# 3. rotation: the source on a circle of radius dso about the rotation axis, at the angle β of the view +function step_rotation_figure(; dso=3.0, dsd=5.0, nu=6, du=0.55, beta=pi / 5) + fig = Figure(size=(520, 470)) + ax = step_axis(fig, (1, 1); title="views: the scanner rotates", limits=(-3.6, 3.6, -3.4, 3.8)) + image_square!(ax) + lines!(ax, [Point2f(dso * cos(t), dso * sin(t)) for t in range(0, 2pi, 200)]; color=(STEP_COLORS.source, 0.4), linestyle=:dash) + # an earlier view, faint: the first one, at β = 0 + s0, _, _, edges0 = topview_geometry(; dso, dsd, nu, du, beta=0.0) + lines!(ax, [Point2f(edges0[1]), Point2f(edges0[end])]; color=(STEP_COLORS.cells[2], 0.6), linewidth=6) + scatter!(ax, [Point2f(s0)]; color=(STEP_COLORS.source, 0.4), markersize=12) + text!(ax, s0...; text="first view", align=(:left, :bottom), offset=(8, 4), fontsize=12, color=:gray30) + # the view at β + s, c, u, edges = topview_geometry(; dso, dsd, nu, du, beta) + for e in edges[[1, end]] + lines!(ax, [Point2f(s), Point2f(e)]; color=STEP_COLORS.ray, linewidth=1) + end + cells!(ax, edges) + scatter!(ax, [Point2f(s)]; color=STEP_COLORS.source, markersize=16) + scatter!(ax, [Point2f(0, 0)]; color=:black, markersize=8) + text!(ax, 0, 0; text="rotation axis", align=(:left, :top), offset=(6, -6), fontsize=13) + dimension!(ax, s, [0.0, 0.0], -0.35 * normalize([s[2], -s[1]]), "dso") + arc = [Point2f(-1.4 * sin(t), 1.4 * cos(t)) for t in range(0, beta, 30)] + lines!(ax, [Point2f(0, 0), Point2f(0, 1.9)]; color=:gray40, linewidth=0.8, linestyle=:dot) + lines!(ax, arc; color=:black) + text!(ax, arc[end÷2]...; text="β", align=(:left, :bottom), offset=(-4, 6), fontsize=16) + fig +end + +# 4. cone beam in 3D, from the angle of the other 3D views (`scene_axis` of geometry_figures.jl): the +# detector has nu × nv cells of size du × dv, drawn as a checkerboard, along u and v +function step_cone_figure(; nu=8, nv=6, du=20.0, dv=20.0, dso=128.0, dsd=256.0, beta=0.0) + g = circular_cone_beam(; angles=[beta], dso, dsd, nu, nv, du, dv) + image = ImageGeom(; dims=(64, 64, 64), deltas=(1.0, 1.0, 1.0)) + s = Point3f(g.source[:, 1]) + c, u, v = g.center[:, 1], g.u[:, 1], g.v[:, 1] + corner(a, b) = Point3f(c + (a - nu / 2) * u + (b - nv / 2) * v) # a in 0:nu, b in 0:nv + fig = Figure(size=(620, 480)) + ax = scene_axis(fig, (1, 1); title="cone beam: nu × nv detector cells") + draw_volume!(ax, image) + # the cells as a checkerboard + for a in 0:nu-1, b in 0:nv-1 + quad = [corner(a, b), corner(a + 1, b), corner(a + 1, b + 1), corner(a, b + 1)] + mesh!(ax, quad, [1 2 3; 1 3 4]; color=STEP_COLORS.cells[isodd(a + b) ? 1 : 2], shading=NoShading) + end + # the cone: the rays to the four corners of the detector, and a few rays to cells + for (a, b) in ((0, 0), (nu, 0), (nu, nv), (0, nv)) + lines!(ax, [s, corner(a, b)]; color=:steelblue, linewidth=1) + end + for (a, b) in ((2, 2), (5, 3), (3, 4)) + lines!(ax, [s, corner(a + 0.5, b + 0.5)]; color=STEP_COLORS.ray, linewidth=1) + end + scatter!(ax, [s]; color=STEP_COLORS.source, markersize=14) + text!(ax, s; text="source", align=(:right, :bottom), offset=(-6, 6)) + # u and v from the first cell, and the counts + o = corner(0, 0) + for (vec, name) in ((u, "u"), (v, "v")) + tip = o + Point3f(60 * normalize(vec)) + lines!(ax, [o, tip]; color=:black, linewidth=2) + text!(ax, tip; text=name, fontsize=16, align=(:left, :bottom)) + end + text!(ax, corner(nu / 2, nv); text="nu cells", align=(:center, :bottom), offset=(0, 8), fontsize=14) + text!(ax, corner(0, nv / 2); text="nv rows", align=(:right, :center), offset=(-8, 0), fontsize=14) + # one cell outlined: its size du × dv + cell = [corner(0, nv - 1), corner(1, nv - 1), corner(1, nv), corner(0, nv), corner(0, nv - 1)] + lines!(ax, cell; color=STEP_COLORS.source, linewidth=2.5) + text!(ax, corner(0, nv); text="du × dv", align=(:right, :bottom), offset=(-4, 4), fontsize=13) + fig +end + +# 5. helical scan: the height of the source over the angle, rising by pitch per turn from z0 +function step_helix_figure(; pitch=24.0, z0=-32.0, turns=2.5) + angles = range(0, 2pi * turns, 200) + g = helical_cone_beam(; angles, pitch, z0, dso=128, dsd=256, nu=8, nv=4, du=1.0) + fig = Figure(size=(560, 400)) + ax = Axis(fig[1, 1]; xlabel="angle β of the view (turns)", ylabel="height z of source and detector", + title="helical scan") + t = angles ./ 2pi + lines!(ax, t, g.source[3, :]; color=STEP_COLORS.source, linewidth=2) + hlines!(ax, [0]; color=:gray60, linestyle=:dot) + scatter!(ax, [0], [z0]; color=STEP_COLORS.source, markersize=12) + text!(ax, 0, z0; text="first view: z0", align=(:left, :top), offset=(8, -4), fontsize=14) + # pitch: the rise in one turn + lines!(ax, [1, 2], [z0 + pitch, z0 + pitch]; color=:gray40, linestyle=:dash) + lines!(ax, [2, 2], [z0 + pitch, z0 + 2pitch]; color=:black) + text!(ax, 2, z0 + 1.5pitch; text="pitch\n(one turn)", align=(:left, :center), offset=(8, 0), fontsize=14) + ylims!(ax, z0 - 12, z0 + pitch * turns + 6) + fig +end From a6663e1f5b56330174b810621d8fcc39e30c9147 Mon Sep 17 00:00:00 2001 From: Niklas Hackelberg Date: Mon, 2 Nov 2026 20:58:49 +0100 Subject: [PATCH 4/6] Comments: state what the code does, not its history Kernel and launch comments keep the reasons and invariants; how they were found and the measured numbers are left to the commit log and benchmark/README.md. Co-Authored-By: Claude Opus 5.5 --- src/operators/launch.jl | 7 +++---- src/operators/siddon/kernels.jl | 6 +++--- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/src/operators/launch.jl b/src/operators/launch.jl index ff211bc..a1b3efe 100644 --- a/src/operators/launch.jl +++ b/src/operators/launch.jl @@ -1,9 +1,8 @@ # How the projection kernels are launched: the workgroup size and the order in which the threads # take the rays. Thread j (1-based) traces ray `ray_index(order, j)`, or none if that is 0 -# (padding). Measured on an RX 7900 XT (benchmark/README.md): workgroups of 64 instead of the -# backend's default and tiles of 8 × 8 detector cells make forward projections at 256^3 up to -# 2.2 times faster; on the CPU neither changes much. The best values depend on the device, hence -# the configuration. +# (padding). Both affect the speed on GPUs (rays of neighboring detector cells read neighboring +# voxels), and the best values depend on the device, hence the configuration; the defaults are in +# `default_launch`. """ LaunchConfig(; workgroupsize=nothing, tile=nothing) diff --git a/src/operators/siddon/kernels.jl b/src/operators/siddon/kernels.jl index be9d5bc..6164397 100644 --- a/src/operators/siddon/kernels.jl +++ b/src/operators/siddon/kernels.jl @@ -9,12 +9,12 @@ # length, step that axis's voxel index (and the linear index) by +-1. Leaving the volume on the # stepped axis ends the ray. The planes are derived from the voxel, never computed independently # of it: in Float32 an entry point within rounding of a plane could otherwise put voxel and plane -# counter one apart, shifting all later voxels of that axis (found by comparing with Float64). +# counter one apart, shifting all later voxels of that axis. # Plane crossings are computed from the plane index rather than accumulated, so rounding errors do # not add up along the ray. Simultaneous crossings give zero-length segments. # # No closures in this code: a closure capturing a variable that is reassigned (the loop state) is -# boxed and the kernel no longer compiles for GPUs (docs/memory/ka-lessons.md of the project). +# boxed and the kernel no longer compiles for GPUs. @kernel inbounds = true function siddon_forward_kernel!(y, @Const(image), grid, geometry, order) i = ray_index(order, @index(Global, Linear)) @@ -63,7 +63,7 @@ end offset = map(crossing_offset, p, d) alpha = map(crossing, plane, invd, offset, d) # The loop state in scalars, one branch per axis: indexing the small vectors with a run-time axis - # would keep them in memory instead of registers (measured: 9.0 ns per step with SVectors). + # would keep them in memory instead of registers, several times slower per step. nx, ny, nz = dims ix, iy, iz = vox[1], vox[2], vox[3] sx, sy, sz = step[1], step[2], step[3] From 6322d6139886c82e4810c7e22840740cc639117f Mon Sep 17 00:00:00 2001 From: Niklas Hackelberg Date: Mon, 2 Nov 2026 20:58:49 +0100 Subject: [PATCH 5/6] Operators: image and ray masks Projectors apply M_rays A M_image and its exact adjoint. The generic layer applies the masks around forward!/backward!, so every operator gets these semantics; Siddon and Joseph skip masked rays in their kernels (trait handles_raymask) and clip rays to the bounding box of the image mask (VoxelGrid lo/hi). subset and adapt keep the masks. Re-exports the mask types of ImageGeoms. Co-Authored-By: Claude Opus 5.5 --- src/CTOperators.jl | 4 +- src/image/voxelgrid.jl | 34 +++++++++++- src/operators/joseph/joseph.jl | 26 ++++++--- src/operators/joseph/kernels.jl | 20 +++---- src/operators/launch.jl | 5 +- src/operators/masks.jl | 90 +++++++++++++++++++++++++++++++ src/operators/operators.jl | 22 +++++--- src/operators/siddon/kernels.jl | 10 ++-- src/operators/siddon/siddon.jl | 26 ++++++--- test/operators/masks.jl | 93 +++++++++++++++++++++++++++++++++ 10 files changed, 289 insertions(+), 41 deletions(-) create mode 100644 src/operators/masks.jl create mode 100644 test/operators/masks.jl diff --git a/src/CTOperators.jl b/src/CTOperators.jl index ac244e0..c692d69 100644 --- a/src/CTOperators.jl +++ b/src/CTOperators.jl @@ -5,11 +5,11 @@ using StaticArrays using Adapt using KernelAbstractions using Atomix -using ImageGeoms: ImageGeom, axis +using ImageGeoms: ImageGeom, axis, MaskCircle, MaskEllipse, MaskAllButEdge const KA = KernelAbstractions export AbstractCTOperator -export ImageGeom +export ImageGeom, MaskCircle, MaskEllipse, MaskAllButEdge export AbstractProjectionGeometry, AbstractDetectorGeometry, RaySegment, RayLine export RayList, ConeBeam, ArcConeBeam, ParallelBeam export circular_cone_beam, helical_cone_beam, parallel_beam diff --git a/src/image/voxelgrid.jl b/src/image/voxelgrid.jl index 8ed5948..e7602c6 100644 --- a/src/image/voxelgrid.jl +++ b/src/image/voxelgrid.jl @@ -3,16 +3,30 @@ 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`, centered at `z = 0`), element type `T`, -no units, no mask. `first` is the world coordinate of the center of voxel `(1, 1, 1)`; voxel +no units. `first` is the world coordinate of the center of voxel `(1, 1, 1)`; voxel `(i, j, k)` has its center at `first + ((i, j, k) .- 1) .* deltas`. Negative `deltas` (flipped axes, allowed by ImageGeoms) are kept. +`lo` and `hi` are the corners of the bounding box of the mask of the `ImageGeom` (the smallest box +containing its voxels), in voxel coordinates (see [`to_index`](@ref CTOperators.to_index)): +`[1/2, n + 1/2]` per dimension without a mask. Projectors clip their rays to this box (Joseph widens +it by half a voxel across the ray, for the interpolation), so that they do not traverse voxels that +the mask sets to zero. + Units of the `ImageGeom` are stripped; they must be the units of the projection geometry. """ struct VoxelGrid{T} dims::NTuple{3,Int} deltas::SVector{3,T} first::SVector{3,T} + lo::SVector{3,T} + hi::SVector{3,T} +end + +# the whole volume as the box +function VoxelGrid{T}(dims::NTuple{3,Int}, deltas::SVector{3,T}, first::SVector{3,T}) where {T} + half = T(1 // 2) + VoxelGrid{T}(dims, deltas, first, SVector(half, half, half), SVector(map(n -> T(n) + half, dims))) end function VoxelGrid{T}(ig::ImageGeom{D}) where {T<:AbstractFloat,D} @@ -20,7 +34,23 @@ function VoxelGrid{T}(ig::ImageGeom{D}) where {T<:AbstractFloat,D} dims = ntuple(d -> d <= D ? ig.dims[d] : 1, 3) deltas = SVector(ntuple(d -> d <= D ? T(strip_unit(ig.deltas[d])) : one(T), 3)) first = SVector(ntuple(d -> d <= D ? T(strip_unit(axis(ig, d)[1])) : zero(T), 3)) - VoxelGrid{T}(dims, deltas, first) + lo, hi = mask_box(ig.mask, T) + VoxelGrid{T}(dims, deltas, first, lo, hi) +end + +# The box around the true voxels of a mask, in voxel coordinates (3D; a 2D mask covers z in +# [1/2, 3/2]). An empty mask gives lo > hi: rays miss it. +function mask_box(mask::AbstractArray{Bool,D}, ::Type{T}) where {D,T} + half = T(1 // 2) + lo, hi = ones(T, 3) .- half, ones(T, 3) .+ half + for d in 1:D + hits = findall(vec(any(mask; dims=Tuple(e for e in 1:D if e != d)))) + if isempty(hits) + return SVector{3,T}(1, 1, 1), SVector{3,T}(0, 0, 0) + end + lo[d], hi[d] = first(hits) - half, last(hits) + half + end + SVector{3,T}(lo), SVector{3,T}(hi) end # a length without its unit (Unitful quantities), plain numbers unchanged diff --git a/src/operators/joseph/joseph.jl b/src/operators/joseph/joseph.jl index f878871..4821fb5 100644 --- a/src/operators/joseph/joseph.jl +++ b/src/operators/joseph/joseph.jl @@ -1,5 +1,5 @@ """ - JosephProjector(image::ImageGeom, geometry; eltype=Float32, storage=Array, launch=nothing) + JosephProjector(image::ImageGeom, geometry; eltype=Float32, storage=Array, launch=nothing, raymask=nothing) Projector with Joseph's method (Joseph 1982): the image is interpolated bilinearly within each layer of voxels perpendicular to the axis the ray runs most along (its principal axis), and the @@ -25,6 +25,14 @@ Keywords: and projection data passed to the operator must live on that device. * `launch` is a [`LaunchConfig`](@ref) to tune the kernel launches for a device, `nothing` for the defaults of the device. +* `raymask` is `nothing` or an array with one entry per ray, of the data size + (`raymask[u, v, view]`) or a vector, on any device: the rays that are used (`true`). It is stored + as `Bool` on the device of the operator. Masked rays give 0 in the forward projection and + contribute nothing to the back projection. + +The mask of `image` (e.g. `ImageGeom(MaskCircle(); dims, deltas)`, the circular field of view) +selects the voxels: voxels outside count as 0 in the forward projection and get 0 in the back +projection, and the rays only traverse the bounding box of the mask. See [Use Masks](@ref). # Example @@ -42,31 +50,35 @@ Reference: * P. M. Joseph, "An improved algorithm for reprojecting rays through pixel images", IEEE Trans. Med. Imaging 1(3), 192–196, 1982. """ -struct JosephProjector{T,G<:AbstractProjectionGeometry{T},B<:KA.Backend,I<:ImageGeom,L<:Union{Nothing,LaunchConfig}} <: AbstractProjector{T,G} +struct JosephProjector{T,G<:AbstractProjectionGeometry{T},B<:KA.Backend,I<:ImageGeom,L<:Union{Nothing,LaunchConfig},M<:ProjectorMasks} <: AbstractProjector{T,G} image::I grid::VoxelGrid{T} geometry::G backend::B launch::L + masks::M end function JosephProjector(image::ImageGeom, geometry::AbstractProjectionGeometry; eltype::Type{T}=Float32, - storage=Array, launch::Union{Nothing,LaunchConfig}=nothing) where {T<:AbstractFloat} - joseph_projector(image, device_geometry(geometry, T, storage), launch) + storage=Array, launch::Union{Nothing,LaunchConfig}=nothing, raymask=nothing) where {T<:AbstractFloat} + joseph_projector(image, device_geometry(geometry, T, storage), launch, raymask) end JosephProjector(image::ImageGeom, geometry; kwargs...) = JosephProjector(image, projection_geometry(geometry); kwargs...) # from a geometry already in the operator's element type and on its device -function joseph_projector(image::ImageGeom, geometry::AbstractProjectionGeometry{T}, launch) where {T} +function joseph_projector(image::ImageGeom, geometry::AbstractProjectionGeometry{T}, launch, raymask) where {T} backend = KA.get_backend(storage_array(geometry)) - JosephProjector(image, VoxelGrid{T}(image), geometry, backend, launch) + JosephProjector(image, VoxelGrid{T}(image), geometry, backend, launch, ProjectorMasks(image, geometry, raymask)) end geometry(op::JosephProjector) = op.geometry image_geometry(op::JosephProjector) = op.image -similar_projector(op::JosephProjector, geometry::AbstractProjectionGeometry) = joseph_projector(op.image, geometry, op.launch) +similar_projector(op::JosephProjector, geometry::AbstractProjectionGeometry; raymask=ray_mask(op)) = + joseph_projector(op.image, geometry, op.launch, raymask) +projector_masks(op::JosephProjector) = op.masks +handles_raymask(::JosephProjector) = true KA.get_backend(op::JosephProjector) = op.backend forward!(y, op::JosephProjector, x) = launch_rays!(joseph_forward_kernel!, op, y, vec(x)) diff --git a/src/operators/joseph/kernels.jl b/src/operators/joseph/kernels.jl index fe3635a..7b568a8 100644 --- a/src/operators/joseph/kernels.jl +++ b/src/operators/joseph/kernels.jl @@ -18,16 +18,16 @@ # # No closures in this code (see siddon/kernels.jl). -@kernel inbounds = true function joseph_forward_kernel!(y, @Const(image), grid, geometry, order) +@kernel inbounds = true function joseph_forward_kernel!(y, @Const(image), grid, geometry, order, raymask) i = ray_index(order, @index(Global, Linear)) if i > 0 - y[i] = joseph_ray(Gather(image), grid, geometry, i) + y[i] = is_masked(raymask, i) ? zero(eltype(y)) : joseph_ray(Gather(image), grid, geometry, i) end end -@kernel inbounds = true function joseph_backward_kernel!(image, @Const(y), grid, geometry, order) +@kernel inbounds = true function joseph_backward_kernel!(image, @Const(y), grid, geometry, order, raymask) i = ray_index(order, @index(Global, Linear)) - if i > 0 + if i > 0 && !is_masked(raymask, i) v = y[i] if !iszero(v) joseph_ray(Scatter(image, v), grid, geometry, i) @@ -64,10 +64,12 @@ end end # the region where the image can be nonzero for rays with principal axis k: the slabs of the -# layers along k, and up to the next voxel center outside the volume across k -@inline function joseph_support(dims, k, ::Type{T}) where {T} - lo = SVector(ntuple(m -> m == k ? T(1 // 2) : zero(T), 3)) - hi = SVector(ntuple(m -> m == k ? T(dims[m]) + T(1 // 2) : T(dims[m] + 1), 3)) +# layers along k within the box of the image mask (the volume without a mask), and across k up to +# the next voxel center outside the box +@inline function joseph_support(grid::VoxelGrid{T}, k) where {T} + half = T(1 // 2) + lo = SVector(ntuple(m -> m == k ? grid.lo[m] : grid.lo[m] - half, 3)) + hi = SVector(ntuple(m -> m == k ? grid.hi[m] : grid.hi[m] + half, 3)) lo, hi end @@ -78,7 +80,7 @@ end d = (b - a) ./ grid.deltas k = principal_axis(d) iszero(pick(d, k)) && return zero(T) # a point, not a ray - lo, hi = joseph_support(grid.dims, k, T) + lo, hi = joseph_support(grid, k) amin, amax = clip(p, d, lo, hi, raytype(geometry)) amin < amax || return zero(T) smin, smax = segment_range(raytype(geometry), T) diff --git a/src/operators/launch.jl b/src/operators/launch.jl index a1b3efe..7ff02fd 100644 --- a/src/operators/launch.jl +++ b/src/operators/launch.jl @@ -87,7 +87,8 @@ nthreads(o::TileOrder, geometry) = cld(o.nu, o.tu) * o.tu * cld(o.nv, o.tv) * o. ifelse(iu < o.nu && iv < o.nv, 1 + iu + o.nu * (iv + o.nv * k), 0) end -# Launches `kernel(out, in, grid, geometry, order)` of the projector `op` over its rays. +# Launches `kernel(out, in, grid, geometry, order, raymask)` of the projector `op` over its rays; +# `raymask` is `nothing` or the ray mask on the device (the kernels skip masked rays). function launch_rays!(kernel, op, out, in) g = geometry(op) nrays(g) == 0 && return out @@ -96,6 +97,6 @@ function launch_rays!(kernel, op, out, in) order = ray_order(config, g) wg = config.workgroupsize k = isnothing(wg) ? kernel(backend) : kernel(backend, wg) - k(out, in, op.grid, g, order; ndrange=nthreads(order, g)) + k(out, in, op.grid, g, order, ray_mask(op); ndrange=nthreads(order, g)) out end diff --git a/src/operators/masks.jl b/src/operators/masks.jl new file mode 100644 index 0000000..9ea6080 --- /dev/null +++ b/src/operators/masks.jl @@ -0,0 +1,90 @@ +# Masks: an image mask (the voxels that are reconstructed, from the mask of the `ImageGeom`) and a +# ray mask (the rays that are used, the projector keyword `raymask`); `true` means used, `false` +# masked. A projector with masks applies +# +# forward: y = M_rays A M_image x back: x = M_image Aᵀ M_rays y +# +# (the exact adjoint), with M the diagonal 0/1 matrices of the masks. The generic layer (`mul!`) +# applies the masks around `forward!`/`backward!`, so every operator gets these semantics; operators +# that apply a mask themselves say so with a trait (`handles_raymask`), and the generic layer leaves +# that mask to them. + +""" + ProjectorMasks(image::ImageGeom, geometry, raymask) + +The masks of a projector on the device of `geometry`: the image mask of `image` (`nothing` if all +voxels are used) and the ray mask `raymask` (`nothing`, or a `Bool` array of the data size, `true` +for the rays that are used). For implementations of projectors. +""" +struct ProjectorMasks{I,R} + image::I + rays::R +end + +const NO_MASKS = ProjectorMasks(nothing, nothing) + +function ProjectorMasks(image::ImageGeom, geometry::AbstractProjectionGeometry, raymask) + backend = KA.get_backend(storage_array(geometry)) + imagemask = all(image.mask) ? nothing : device_mask(backend, image.mask) + if !isnothing(raymask) + length(raymask) == nrays(geometry) || + throw(DimensionMismatch("raymask has $(length(raymask)) entries for $(nrays(geometry)) rays $(datasize(geometry))")) + raymask = reshape(device_mask(backend, raymask), datasize(geometry)) + end + ProjectorMasks(imagemask, raymask) +end + +# a Bool array on the device of `backend` with the values of `mask` (from any device) +device_mask(backend, mask) = copyto!(KA.allocate(backend, Bool, size(mask)), Array{Bool}(Array(mask))) + +""" + projector_masks(op::AbstractCTOperator) -> ProjectorMasks + +The masks of `op`. Projectors of this package store them. Default for an `AbstractProjector`: the +mask of its `ImageGeom`, copied to the device at every call, and no ray mask. Default for any other +`AbstractCTOperator`: no masks. +""" +projector_masks(op::AbstractCTOperator) = NO_MASKS +projector_masks(op::AbstractProjector) = ProjectorMasks(image_geometry(op), geometry(op), nothing) + +""" + image_mask(op) -> Union{Nothing,AbstractArray{Bool}} + ray_mask(op) -> Union{Nothing,AbstractArray{Bool}} + +The image mask (the voxels that are reconstructed, of the image size) and the ray mask (the rays +that are used, of the data size) of `op` on its device, `nothing` without a mask. +""" +image_mask(op::AbstractCTOperator) = projector_masks(op).image +ray_mask(op::AbstractCTOperator) = projector_masks(op).rays +@doc (@doc image_mask) ray_mask + +""" + handles_raymask(op) -> Bool + +Whether `op` applies its ray mask itself (masked rays give 0 in the forward projection and add +nothing in the back projection); otherwise the generic layer multiplies by the mask. Default +`false`; the projectors of this package define it as `true` and skip masked rays in their kernels. +""" +handles_raymask(op::AbstractCTOperator) = false + +# the operator with its masks, as `mul!` applies it +apply!(y, op::AbstractCTOperator, x) = masked_forward!(y, op, x) +apply!(x, a::AdjointCTOperator, y) = masked_backward!(x, parent(a), y) + +function masked_forward!(y, op, x) + masks = projector_masks(op) + forward!(y, op, isnothing(masks.image) ? x : x .* masks.image) + isnothing(masks.rays) || handles_raymask(op) || (y .*= masks.rays) + y +end + +function masked_backward!(x, op, y) + masks = projector_masks(op) + backward!(x, op, isnothing(masks.rays) || handles_raymask(op) ? y : y .* masks.rays) + isnothing(masks.image) || (x .*= masks.image) + x +end + +# is ray i masked out? Inside kernels; without a mask the check compiles away +@inline is_masked(::Nothing, i) = false +@inline is_masked(mask, i) = @inbounds !mask[i] diff --git a/src/operators/operators.jl b/src/operators/operators.jl index 85b4f34..bd93d15 100644 --- a/src/operators/operators.jl +++ b/src/operators/operators.jl @@ -8,7 +8,8 @@ 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`) + options, for another geometry (used by `subset` and `adapt`); a projector that stores a ray mask + also accepts the keyword `raymask` (default `ray_mask(op)`) * the methods of [`AbstractCTOperator`](@ref) """ abstract type AbstractProjector{T,G<:AbstractProjectionGeometry} <: AbstractCTOperator{T} end @@ -28,9 +29,10 @@ The `ImageGeom` the operator was built for. function image_geometry end """ - similar_projector(op::AbstractProjector, geometry) + similar_projector(op::AbstractProjector, geometry; raymask=ray_mask(op)) A projector of the same kind and with the same image geometry and options as `op`, for `geometry`. +The keyword `raymask` exists only for projectors that store a ray mask. """ function similar_projector end @@ -85,11 +87,12 @@ KA.get_backend(a::AdjointCTOperator) = KA.get_backend(a.parent) forward!(x, a::AdjointCTOperator, y) = backward!(x, a.parent, y) backward!(y, a::AdjointCTOperator, x) = forward!(y, a.parent, x) -# Shaped arrays and flat vectors (RegularizedLeastSquares works on vectors) both work. +# Shaped arrays and flat vectors (RegularizedLeastSquares works on vectors) both work. `apply!` +# (masks.jl) applies the operator with its masks. function LinearAlgebra.mul!(y::AbstractArray, op::AbstractCTOperator, x::AbstractArray) check_length(x, imagesize(op), "input") check_length(y, datasize(op), "output") - forward!(reshape(y, datasize(op)), op, reshape(x, imagesize(op))) + apply!(reshape(y, datasize(op)), op, reshape(x, imagesize(op))) y end @@ -141,10 +144,14 @@ KA.get_backend(op::AbstractProjector) = KA.get_backend(storage_array(geometry(op subset(op::AbstractProjector, views) The projector of the views `views` of `op` (detector geometries; a vector or range of view -indices), for ordered-subset methods. It shares the image geometry and options with `op`. +indices), for ordered-subset methods. It shares the image geometry and options with `op`, and the +views of its ray mask. """ -subset(op::AbstractProjector{T,<:AbstractDetectorGeometry}, views) where {T} = - similar_projector(op, subset(geometry(op), views)) +function subset(op::AbstractProjector{T,<:AbstractDetectorGeometry}, views) where {T} + g = subset(geometry(op), views) + rays = ray_mask(op) + isnothing(rays) ? similar_projector(op, g) : similar_projector(op, g; raymask=Array(rays)[:, :, views]) +end Adapt.adapt_structure(to, op::AbstractProjector) = similar_projector(op, adapt(to, geometry(op))) @@ -158,6 +165,7 @@ with their keywords `eltype` and `storage`. device_geometry(geometry, ::Type{T}, storage) where {T} = adapt(storage isa UnionAll ? storage{T} : storage, geometry) +include("masks.jl") include("launch.jl") include("siddon/siddon.jl") include("joseph/joseph.jl") diff --git a/src/operators/siddon/kernels.jl b/src/operators/siddon/kernels.jl index 6164397..926d127 100644 --- a/src/operators/siddon/kernels.jl +++ b/src/operators/siddon/kernels.jl @@ -16,16 +16,16 @@ # No closures in this code: a closure capturing a variable that is reassigned (the loop state) is # boxed and the kernel no longer compiles for GPUs. -@kernel inbounds = true function siddon_forward_kernel!(y, @Const(image), grid, geometry, order) +@kernel inbounds = true function siddon_forward_kernel!(y, @Const(image), grid, geometry, order, raymask) i = ray_index(order, @index(Global, Linear)) if i > 0 - y[i] = siddon_ray(Gather(image), grid, geometry, i) + y[i] = is_masked(raymask, i) ? zero(eltype(y)) : siddon_ray(Gather(image), grid, geometry, i) end end -@kernel inbounds = true function siddon_backward_kernel!(image, @Const(y), grid, geometry, order) +@kernel inbounds = true function siddon_backward_kernel!(image, @Const(y), grid, geometry, order, raymask) i = ray_index(order, @index(Global, Linear)) - if i > 0 + if i > 0 && !is_masked(raymask, i) v = y[i] if !iszero(v) siddon_ray(Scatter(image, v), grid, geometry, i) @@ -38,7 +38,7 @@ end a, b = ray(geometry, i) p = to_index(grid, a) d = (b - a) ./ grid.deltas - amin, amax = clip(p, d, grid.dims, raytype(geometry)) + amin, amax = clip(p, d, grid.lo, grid.hi, raytype(geometry)) # the box of the image mask amin < amax || return zero(T) siddon_traverse(action, grid.dims, p, d, norm(b - a), amin, amax) end diff --git a/src/operators/siddon/siddon.jl b/src/operators/siddon/siddon.jl index ef07b43..a4b5c54 100644 --- a/src/operators/siddon/siddon.jl +++ b/src/operators/siddon/siddon.jl @@ -1,5 +1,5 @@ """ - SiddonProjector(image::ImageGeom, geometry; eltype=Float32, storage=Array, launch=nothing) + SiddonProjector(image::ImageGeom, geometry; eltype=Float32, storage=Array, launch=nothing, raymask=nothing) 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 @@ -12,6 +12,14 @@ Keywords: and projection data passed to the operator must live on that device. * `launch` is a [`LaunchConfig`](@ref) to tune the kernel launches for a device, `nothing` for the defaults of the device. +* `raymask` is `nothing` or an array with one entry per ray, of the data size + (`raymask[u, v, view]`) or a vector, on any device: the rays that are used (`true`). It is stored + as `Bool` on the device of the operator. Masked rays give 0 in the forward projection and + contribute nothing to the back projection. + +The mask of `image` (e.g. `ImageGeom(MaskCircle(); dims, deltas)`, the circular field of view) +selects the voxels: voxels outside count as 0 in the forward projection and get 0 in the back +projection, and the rays only traverse the bounding box of the mask. See [Use Masks](@ref). # Example @@ -32,31 +40,35 @@ References: calculate the exact radiological path through a pixel or voxel space", J. Comput. Inf. Technol. 6(1), 89–94, 1998. """ -struct SiddonProjector{T,G<:AbstractProjectionGeometry{T},B<:KA.Backend,I<:ImageGeom,L<:Union{Nothing,LaunchConfig}} <: AbstractProjector{T,G} +struct SiddonProjector{T,G<:AbstractProjectionGeometry{T},B<:KA.Backend,I<:ImageGeom,L<:Union{Nothing,LaunchConfig},M<:ProjectorMasks} <: AbstractProjector{T,G} image::I grid::VoxelGrid{T} geometry::G backend::B launch::L + masks::M end function SiddonProjector(image::ImageGeom, geometry::AbstractProjectionGeometry; eltype::Type{T}=Float32, - storage=Array, launch::Union{Nothing,LaunchConfig}=nothing) where {T<:AbstractFloat} - siddon_projector(image, device_geometry(geometry, T, storage), launch) + storage=Array, launch::Union{Nothing,LaunchConfig}=nothing, raymask=nothing) where {T<:AbstractFloat} + siddon_projector(image, device_geometry(geometry, T, storage), launch, raymask) end SiddonProjector(image::ImageGeom, geometry; kwargs...) = SiddonProjector(image, projection_geometry(geometry); kwargs...) # from a geometry already in the operator's element type and on its device -function siddon_projector(image::ImageGeom, geometry::AbstractProjectionGeometry{T}, launch) where {T} +function siddon_projector(image::ImageGeom, geometry::AbstractProjectionGeometry{T}, launch, raymask) where {T} backend = KA.get_backend(storage_array(geometry)) - SiddonProjector(image, VoxelGrid{T}(image), geometry, backend, launch) + SiddonProjector(image, VoxelGrid{T}(image), geometry, backend, launch, ProjectorMasks(image, geometry, raymask)) end geometry(op::SiddonProjector) = op.geometry image_geometry(op::SiddonProjector) = op.image -similar_projector(op::SiddonProjector, geometry::AbstractProjectionGeometry) = siddon_projector(op.image, geometry, op.launch) +similar_projector(op::SiddonProjector, geometry::AbstractProjectionGeometry; raymask=ray_mask(op)) = + siddon_projector(op.image, geometry, op.launch, raymask) +projector_masks(op::SiddonProjector) = op.masks +handles_raymask(::SiddonProjector) = true KA.get_backend(op::SiddonProjector) = op.backend forward!(y, op::SiddonProjector, x) = launch_rays!(siddon_forward_kernel!, op, y, vec(x)) diff --git a/test/operators/masks.jl b/test/operators/masks.jl new file mode 100644 index 0000000..d181826 --- /dev/null +++ b/test/operators/masks.jl @@ -0,0 +1,93 @@ +# Masks (src/operators/masks.jl): a projector with an image mask (from its ImageGeom) and a ray +# mask (keyword `raymask`) applies M_rays A M_image, the exact adjoint M_image Aᵀ M_rays, on every +# device. The projectors of the package skip masked rays in their kernels and clip the rays to the +# box of the image mask; the results must equal those of the unmasked projector applied to masked +# inputs. + +@testsnippet MaskProblem begin + using LinearAlgebra, Adapt + using CTOperators: image_mask, ray_mask, handles_raymask + + # a 3D image with the circular field of view plus a box-shaped region of interest, a helical + # cone beam and a random ray mask with about 20 % masked rays + function mask_problem() + dims = (10, 9, 7) + circle = ImageGeom(MaskCircle(); dims, deltas=(1.0, 1.25, 0.8), offsets=(0.37, 0.11, -0.23)) + roi = falses(dims...) + roi[3:7, 2:6, 2:5] .= true + box = ImageGeom(; dims, deltas=(1.0, 1.25, 0.8), offsets=(0.37, 0.11, -0.23), mask=roi) + g = helical_cone_beam(; angles=range(0, 3pi, 7), pitch=4, dso=30, dsd=60, nu=13, nv=9, du=1.17, dv=0.93) + raymask = rand(datasize(g)...) .> 0.2 + circle, box, g, raymask + end +end + +@testitem "masks: semantics and adjoint" tags = [:operators, :gpu] setup = [Backends, MaskProblem] begin + using Random + Random.seed!(3) + circle, box, g, raymask = mask_problem() + for P in (SiddonProjector, JosephProjector), image in (circle, box), arrayType in arrayTypes, T in floatTypes + @testset "$P, mask of $(count(image.mask)) voxels, $arrayType{$T}" begin + op = P(image, g; eltype=T, storage=arrayType, raymask) + plain = P(ImageGeom(; dims=image.dims, deltas=image.deltas, offsets=image.offsets), g; eltype=Float64) + m = Array(image.mask) + x, y = rand(image.dims...), rand(datasize(g)...) + fx = Array(op * adapt(arrayType{T}, x)) + by = Array(adjoint(op) * adapt(arrayType{T}, y)) + # equal to the unmasked projector on masked inputs (also checks the clipping to the mask's box) + @test isapprox(fx, raymask .* (plain * (m .* x)); rtol=sqrt(eps(T))) + @test isapprox(by, m .* (adjoint(plain) * (raymask .* y)); rtol=sqrt(eps(T))) + # exact zeros outside the masks + @test all(iszero, fx[.!raymask]) && all(iszero, by[.!m]) + # exact adjoint + xd, yd = adapt(arrayType{T}, rand(T, image.dims...)), adapt(arrayType{T}, rand(T, datasize(g)...)) + @test isapprox(sum(vec(op * xd) .* vec(yd)), sum(vec(xd) .* vec(adjoint(op) * yd)); rtol=100eps(T)) + end + end +end + +@testitem "masks: interface, subset, adapt" tags = [:operators] setup = [MaskProblem] begin + circle, box, g, raymask = mask_problem() + op = JosephProjector(circle, g; eltype=Float64, raymask) + @test image_mask(op) == circle.mask && ray_mask(op) == raymask && handles_raymask(op) + plain = JosephProjector(ImageGeom(; dims=circle.dims), g) + @test isnothing(image_mask(plain)) && isnothing(ray_mask(plain)) + # subset takes the views of the ray mask; adapt and similar projectors keep the masks + sub = subset(op, [2, 5]) + @test ray_mask(sub) == raymask[:, :, [2, 5]] + x = rand(circle.dims...) + @test sub * x ≈ (op * x)[:, :, [2, 5]] + op32 = adapt(Array{Float32}, op) + @test ray_mask(op32) == raymask && image_mask(op32) == circle.mask + # the ray mask also as a vector, and with the wrong size + @test ray_mask(JosephProjector(circle, g; raymask=vec(raymask))) == raymask + @test_throws DimensionMismatch JosephProjector(circle, g; raymask=trues(3)) + # an empty image mask: no voxel, all projections zero + none = ImageGeom(; dims=(6, 5, 4), mask=falses(6, 5, 4)) + @test all(iszero, SiddonProjector(none, g; eltype=Float64) * rand(6, 5, 4)) +end + +# A projector of another package that knows nothing about masks gets the semantics of the image mask +# from the generic layer (the default `projector_masks` takes the mask of its ImageGeom). +@testitem "masks: generic layer for other projectors" tags = [:operators] setup = [MaskProblem] begin + using CTOperators: AbstractProjector + import CTOperators: geometry, image_geometry, forward!, backward! + struct WrappedProjector{T,G,P} <: AbstractProjector{T,G} + image::ImageGeom + inner::P # a projector without mask, doing the work + end + WrappedProjector(image, g) = + WrappedProjector{Float64,typeof(g),Any}(image, SiddonProjector(ImageGeom(; dims=image.dims), g; eltype=Float64)) + geometry(op::WrappedProjector) = geometry(op.inner) + image_geometry(op::WrappedProjector) = op.image + CTOperators.KA.get_backend(op::WrappedProjector) = CTOperators.KA.get_backend(op.inner) + forward!(y, op::WrappedProjector, x) = forward!(y, op.inner, x) + backward!(x, op::WrappedProjector, y) = backward!(x, op.inner, y) + + circle, box, g, raymask = mask_problem() + op = WrappedProjector(box, g) + @test !handles_raymask(op) && image_mask(op) == box.mask + x, y = rand(box.dims...), rand(datasize(g)...) + @test op * x ≈ op.inner * (box.mask .* x) + @test adjoint(op) * y ≈ box.mask .* (adjoint(op.inner) * y) +end From 927fa7988ece3bb1f348b9cf114f65e33a028d7d Mon Sep 17 00:00:00 2001 From: Niklas Hackelberg Date: Mon, 2 Nov 2026 20:58:49 +0100 Subject: [PATCH 6/6] Docs: masks how-to and explanation, mask benchmarks How-to "Use Masks" (circular field of view, defective detector cells, region of interest) and explanation "Masks" (semantics, adjointness, where masks save time, masks for other operators), a section on masks in the implementation guide, and GPU timings with masks in the benchmark README. The navigation collapses to the first level. Co-Authored-By: Claude Opus 5.5 --- benchmark/README.md | 20 ++++ docs/literate/explanation/masks.jl | 74 +++++++++++++ docs/literate/howto/implement_operator.jl | 14 +++ docs/literate/howto/masks.jl | 120 ++++++++++++++++++++++ docs/make.jl | 3 + docs/src/index.md | 1 + docs/utils/mask_figures.jl | 49 +++++++++ 7 files changed, 281 insertions(+) create mode 100644 docs/literate/explanation/masks.jl create mode 100644 docs/literate/howto/masks.jl create mode 100644 docs/utils/mask_figures.jl diff --git a/benchmark/README.md b/benchmark/README.md index f8e73fa..fa26950 100644 --- a/benchmark/README.md +++ b/benchmark/README.md @@ -137,3 +137,23 @@ Observations: against 102 ms), slower at 128³ (20 ms against 4.9 ms). - Joseph's back projection takes 4–5× its forward projection (four atomic additions per sample), Siddon's about 2.3× (one per voxel): the next target. + +### Masks (GPU) + +AMD Radeon RX 7900 XT, Float32, cone256 (256³ volume, 256² detector, 360 views), minimum time, +2026-11-02, branch `nh/masks` (the GPU was also the desktop GPU in use: approximate). Image masks: +the circular field of view (`MaskCircle`) and a box of 64³ voxels in the center; ray mask: 20 % of +the rays masked at random. + +| image mask | ray mask | Siddon forward / back | Joseph forward / back | +|---|---|---|---| +| none | none | 111 / 264 ms | 113 / 537 ms | +| none | 20 % masked | 96 / 212 ms | 108 / 438 ms | +| circle | none | 111 / 257 ms | 115 / 542 ms | +| box 64³ | none | 2.6 / 6.7 ms | 3.1 / 18.0 ms | +| box 64³ | 20 % masked | 2.2 / 5.5 ms | 2.9 / 15.0 ms | + +Rays are clipped to the bounding box of the image mask: a small region of interest makes the +projections 30 to 40× faster; the circle fills its box and gains nothing. +Masked rays are not traced: the back projection gets faster by about the masked fraction. + diff --git a/docs/literate/explanation/masks.jl b/docs/literate/explanation/masks.jl new file mode 100644 index 0000000..2dd326a --- /dev/null +++ b/docs/literate/explanation/masks.jl @@ -0,0 +1,74 @@ +# # Masks +# +# A mask removes voxels or rays from the problem. This page explains what a projector with masks +# computes, why it stays the exact adjoint, where it saves time, and how operators of other packages +# take part. [Use Masks](@ref) shows how to set the masks. + +using CTOperators +include(joinpath(pkgdir(CTOperators), "docs", "utils", "mask_figures.jl")) #hide + +# ## What a mask means +# +# A mask is a `Bool` array: `true` means "used", `false` means "masked" (excluded). +# The image mask (the mask of the `ImageGeom`) and the ray mask (the keyword `raymask`) are diagonal +# matrices with 0 and 1 on the diagonal, ``M_\text{image}`` and ``M_\text{rays}``. A projector ``A`` +# with masks computes +# +# ```math +# y = M_\text{rays} \, A \, M_\text{image} \, x, \qquad x = M_\text{image} \, A^T M_\text{rays} \, y +# ``` +# +# for the forward and the back projection. Since the masks are diagonal, the second is the transpose +# of the first: the pair stays exactly adjoint, as iterative solvers and automatic differentiation +# need. It is the same as removing the masked columns (voxels) and rows (rays) from the system +# matrix, as the [ASTRA Toolbox](https://astra-toolbox.com) does, except that images and data keep +# their full shape, with zeros at the masked entries. +# +# For a solver this means: masked voxels are not unknowns (they stay zero), and masked rays are not +# measurements (their values do not matter). The row and column sums, the normalizations of SART and +# OS-EM, are those of the masked matrix. +# +# ## Where masks save time +# +# A projector that applied the masks only by multiplying would compute everything and then throw +# parts away. The projectors of this package use the masks to compute less: +# +# * **Masked rays are not traced.** The kernel checks the ray mask before it sets up a ray: a masked +# ray writes 0 in the forward projection and returns at once in the back projection. With 20 % of +# the rays masked, the back projection, whose atomic additions dominate, gets about 20 % faster. +# * **Rays are clipped to the bounding box of the image mask.** Every ray is clipped to a box before +# it is traced; without a mask the box is the volume, with a mask it is the smallest box that +# contains the mask (for Joseph widened by half a voxel across the ray, for the interpolation). +# Outside the box the image is zero, so nothing is lost. In the figure, rays are solid where they +# are traced and dashed where they are not; the bottom ray misses the box and is not traced at all: + +fig = mask_box_figure() #hide +save("mask_box.png", fig) #hide +nothing #hide +# ![Rays traced only inside the bounding box of the mask](mask_box.png) + +# For a small region of interest this saves most of the work: a mask of 64³ voxels in a 256³ volume +# made the projections 30 to 40 times faster in our +# [benchmarks](https://github.com/JuliaImageRecon/CTOperators.jl/blob/main/benchmark/README.md). The +# circular field of view (`MaskCircle`) fills its bounding box, the whole square, and gains nothing. +# +# Inside the box, masked voxels are still visited, as by the top ray in the figure, which crosses a +# corner of the box outside the ellipse: the forward projection reads the image multiplied by the +# mask, and the back projection writes them and sets them to zero afterwards. +# +# ## Masks for every operator +# +# The masks are applied in the generic layer of the operators (`mul!`, which calls the operator's +# `forward!` and `backward!`), so every operator gets the semantics above, also one implemented +# outside this package (see [Implement a CT Operator](@ref)): +# +# * `CTOperators.image_mask(op)` and `CTOperators.ray_mask(op)` return the masks of `op` on its +# device, or `nothing`. A subtype of `AbstractProjector` from another package has the mask of its +# `ImageGeom` and no ray mask (there is no way to give it one from outside); other operators have +# no masks. +# * `CTOperators.handles_raymask(op)` states whether `op` applies its ray mask itself, as the +# projectors of this package do in their kernels. If it does not, the generic layer multiplies the +# data by the ray mask. +# +# A projector of another package that ignores masks is therefore still correct; it only computes the +# masked parts. [Implement a CT Operator](@ref) shows how a projector uses the masks itself. diff --git a/docs/literate/howto/implement_operator.jl b/docs/literate/howto/implement_operator.jl index 3fdcb33..7b69d8c 100644 --- a/docs/literate/howto/implement_operator.jl +++ b/docs/literate/howto/implement_operator.jl @@ -217,6 +217,20 @@ y_joseph = JosephProjector(image, geometry) * x_true # Further tests check what the method promises, for example exact line integrals for images the # method represents exactly, compared with the analytic projections of ImagePhantoms.jl. # +# ## Support for masks +# +# The generic layer applies the masks of an operator (see [Masks](@ref)). The sampling projector +# above gets the image mask of its `ImageGeom` this way, applied by multiplication, and has no ray +# mask. A projector can do more in two independent steps: +# +# * **Clip the rays to the bounding box of the image mask:** clip with `grid.lo` and `grid.hi` of the +# `VoxelGrid` instead of the volume. +# * **Store a ray mask and skip masked rays:** store a `CTOperators.ProjectorMasks` built from the +# keyword `raymask`, return it from `CTOperators.projector_masks`, pass `CTOperators.ray_mask(op)` +# to the kernels (`CTOperators.is_masked(raymask, i)`), define +# `CTOperators.handles_raymask(::MyProjector) = true`, and accept the keyword +# `raymask=CTOperators.ray_mask(op)` in `similar_projector`, which `subset` and `adapt` call. +# # ## Wrapping another package # # A projector of another package, with its own geometry, becomes an `AbstractCTOperator` with five diff --git a/docs/literate/howto/masks.jl b/docs/literate/howto/masks.jl new file mode 100644 index 0000000..3e5bddc --- /dev/null +++ b/docs/literate/howto/masks.jl @@ -0,0 +1,120 @@ +# # Use Masks +# +# A projector has two masks, both `Bool` arrays in which `true` means "used" and `false` means +# "masked" (excluded). The **image mask** selects the voxels that are reconstructed; it is the mask +# of the `ImageGeom`. The **ray mask** selects the rays whose data is used; it is the keyword +# `raymask` of the projector. Masked voxels are zero in the image and in the back projection; masked +# rays are zero in the data and add nothing to the back projection. The projectors skip masked rays +# and trace rays only through the bounding box of the image mask. [Masks](@ref) explains what +# happens inside. + +using CTOperators, ImagePhantoms, RegularizedLeastSquares, LinearOperatorCollection, CairoMakie +using ImageGeoms: axis +include(joinpath(pkgdir(CTOperators), "docs", "utils", "tutorial_figures.jl")) #hide + +# The examples use the 2D Shepp-Logan phantom on 128 × 128 pixels and a parallel-beam scan with 180 +# views, reconstructed with 30 CGNR iterations: + +n = 128 +x_true = Float32.(phantom(axis(ImageGeom(; dims=(n, n)), 1), axis(ImageGeom(; dims=(n, n)), 2), + shepp_logan(SheppLoganToft(); fovs=(n, n)), 2)) +geometry = parallel_beam(; angles=range(0, pi, 181)[1:180], nu=192, nv=1, du=1.0) + +function reconstruct(op, y) + solver = createLinearSolver(CGNR, op; AHA=normalOperator(op), iterations=30) + reshape(solve!(solver, vec(y)), n, n) +end; + +# ## The circular field of view +# +# A rotating scanner sees the full object only within the circle that all views cover. The image +# mask `MaskCircle()` of ImageGeoms.jl restricts the image to this circle (a cylinder in 3D): + +image = ImageGeom(MaskCircle(); dims=(n, n), deltas=(1.0, 1.0)) +op = JosephProjector(image, geometry) +x = reconstruct(op, op * x_true); + +# The reconstruction is zero outside the circle; inside, it is the reconstruction of the data +# restricted to the circle. The phantom lies inside the circle, so its reconstruction is the same as +# without the mask; the mask matters for objects that extend past the field of view, whose parts +# outside the circle are seen only in some views: + +fig = images_figure([("image mask", Float32.(image.mask), (0, 1)), ("reconstruction", x)]; #hide + colorrange=(0, 0.4)) #hide +save("masks_circle.png", fig) #hide +nothing #hide +# ![The circular image mask and the reconstruction](masks_circle.png) + +# Any `Bool` array of the image size is a mask: `ImageGeom(; dims, deltas, mask)`. ImageGeoms.jl also +# offers `MaskEllipse()` and `MaskAllButEdge()`. +# +# ## Defective detector cells +# +# A defective cell measures nothing useful in any view. Here three cells of the detector row report +# zero: + +defective = [70, 71, 120] +image = ImageGeom(; dims=(n, n), deltas=(1.0, 1.0)) +op = JosephProjector(image, geometry) +y = op * x_true +y[defective, :, :] .= 0; + +# A defective cell lies at the same distance from the rotation axis in every view, so its wrong +# values form a circle in the reconstruction: a ring artifact. The ray mask excludes these rays. It +# has the size of the data, `(nu, nv, nviews)` (see `datasize`), so `raymask[u, v, view]` belongs to +# detector cell `(u, v)` in view `view`; a vector with one entry per ray works too. + +raymask = trues(size(y)) +raymask[defective, :, :] .= false +op_masked = JosephProjector(image, geometry; raymask) +x_rings, x_masked = reconstruct(op, y), reconstruct(op_masked, y); + +#- + +fig = images_figure([("phantom", x_true), ("defective cells", x_rings), ("ray mask", x_masked)]; #hide + colorrange=(0, 0.4), truth=x_true) #hide +save("masks_rings.png", fig) #hide +nothing #hide +# ![Defective detector cells, without and with the ray mask](masks_rings.png) + +# The figure shows the reconstructions in the color window 0 to 0.4, where the soft tissue of the +# phantom is visible, and the relative error ``\|x - x_\text{true}\| / \|x_\text{true}\|``. With the +# ray mask the rings almost disappear: the other views measure the same lines with healthy cells. +# The ray mask works the same way for any rays whose data cannot be used, e.g. truncated or +# saturated ones. +# +# ## A region of interest +# +# If only part of the object matters, the image mask can be a box around it, here deliberately off +# center. The rays then only traverse this box, which makes the projections much faster for a small +# region: in a 256³ volume, a mask of 64³ voxels made the projections 30 to 40 times faster in our +# [benchmarks](https://github.com/JuliaImageRecon/CTOperators.jl/blob/main/benchmark/README.md). + +roi = falses(n, n) +roi[33:96, 49:112] .= true +image = ImageGeom(; dims=(n, n), deltas=(1.0, 1.0), mask=roi) +op = JosephProjector(image, geometry) +x = reconstruct(op, op * x_true); + +#- + +fig = images_figure([("image mask", Float32.(roi), (0, 1)), ("reconstruction", x)]; colorrange=(0, 0.4)) #hide +save("masks_roi.png", fig) #hide +nothing #hide +# ![A region of interest](masks_roi.png) + +# Here the data were simulated with the masked operator, so they contain nothing from outside the +# region. With measured data of the whole object, the region of interest misses the contributions of +# the object outside it: reconstructing a region alone from such data (interior tomography) needs +# more than a mask. +# +# ## Masks and other functions +# +# * `subset(op, views)` keeps the part of the ray mask for these views; `adapt` moves both masks to +# the device. +# * The masks are `Bool` arrays on the device of the operator, whatever its element type; a ray mask +# from any device is copied there. The ray mask takes one byte per ray, and the forward projection +# with an image mask allocates the masked image on every call. +# * `rowsums` and `colsums` include the masks, and so do the normalizations of ordered-subset methods +# built on them. +# * The masks are part of the operator in every solver and in automatic differentiation. diff --git a/docs/make.jl b/docs/make.jl index e7087bb..2809b5b 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -22,6 +22,7 @@ makedocs(; prettyurls=get(ENV, "CI", "false") == "true", # also in utils/bonito_page.jl canonical="https://JuliaImageRecon.github.io/CTOperators.jl", edit_link="main", + collapselevel=1, # the navigation shows the sections, collapsed ), # Diátaxis: tutorials to learn, how-to guides for tasks, explanation for understanding, # reference for looking things up @@ -40,6 +41,7 @@ makedocs(; "Set Up a Geometry" => "generated/howto/geometries.md", "Run on a GPU" => "generated/howto/gpu.md", "Ordered Subsets and Normalization" => "generated/howto/subsets.md", + "Use Masks" => "generated/howto/masks.md", "Differentiate through a Projector" => "generated/howto/differentiation.md", "Use Geometries of Sinograms.jl" => "generated/howto/sinograms.md", "Implement a CT Operator" => "generated/howto/implement_operator.md", @@ -49,6 +51,7 @@ makedocs(; "Coordinates and Conventions" => "generated/explanation/conventions.md", "Siddon's Method" => "generated/explanation/siddon.md", "Joseph's Method" => "generated/explanation/joseph.md", + "Masks" => "generated/explanation/masks.md", "Workgroups and Ray Order on GPUs" => "generated/explanation/launch.md", ], "Reference" => [ diff --git a/docs/src/index.md b/docs/src/index.md index 5632d7a..a547a3e 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -51,6 +51,7 @@ y = op * rand(Float32, 128, 128, 128) # projections, size (256, 128, 360) - The how-to guides show how to [set up a geometry](@ref "Set Up a Geometry"), including an object that moves during the scan, how to [run on a GPU](@ref "Run on a GPU"), how to build [ordered-subset methods](@ref "Ordered Subsets and Normalization"), how to + [use masks](@ref "Use Masks") for the image and the rays, how to [differentiate through a projector](@ref "Differentiate through a Projector"), how to [use geometries of Sinograms.jl](@ref "Use Geometries of Sinograms.jl") and how to [implement a CT operator](@ref "Implement a CT Operator"). diff --git a/docs/utils/mask_figures.jl b/docs/utils/mask_figures.jl new file mode 100644 index 0000000..4ae1f77 --- /dev/null +++ b/docs/utils/mask_figures.jl @@ -0,0 +1,49 @@ +# Figure of the explanation page "Masks" (CairoMakie): rays through a 2D image with a mask; each ray +# is traced only inside the bounding box of the mask. Included invisibly by the page. + +using CairoMakie + +function mask_box_figure(; n=12) + # an elliptical mask; its box is computed from it as in the projectors + mask = [((i - 6.5) / 3.2)^2 + ((j - 8) / 2.6)^2 <= 1 for i in 1:n, j in 1:n] + rows, cols = findall(vec(any(mask; dims=2))), findall(vec(any(mask; dims=1))) + fig = Figure(size=(520, 460)) + ax = Axis(fig[1, 1]; aspect=DataAspect(), title="rays: solid where traced, dashed where not", + limits=(-0.5, n + 1.5, -0.5, n + 1.5)) + hidedecorations!(ax); hidespines!(ax) + heatmap!(ax, 1:n, 1:n, Float32.(mask); colormap=[:white, (:steelblue, 0.45)], colorrange=(0, 1)) + for w in 0.5:1:n+0.5 + lines!(ax, [w, w], [0.5, n + 0.5]; color=:gray80, linewidth=0.6) + lines!(ax, [0.5, n + 0.5], [w, w]; color=:gray80, linewidth=0.6) + end + # the box around the mask + lo, hi = (first(rows) - 0.5, first(cols) - 0.5), (last(rows) + 0.5, last(cols) + 0.5) + poly!(ax, Rect2f(lo..., (hi .- lo)...); color=:transparent, strokecolor=:steelblue, strokewidth=2) + text!(ax, hi[1], hi[2]; text="bounding box of the mask", align=(:left, :bottom), offset=(4, 2), fontsize=13, color=:steelblue) + # rays: dashed where they are not traced, solid inside the box + for (a, b) in (((-0.5, 2.0), (13.5, 9.0)), ((-0.5, 7.5), (13.5, 11.0)), ((2.0, -0.5), (11.0, 13.5)), ((-0.5, 1.0), (13.5, 2.5))) + lines!(ax, [Point2f(a), Point2f(b)]; color=(:darkorange, 0.5), linestyle=:dash) + t0, t1 = clip_segment(a, b, lo, hi) + if t0 < t1 + p0, p1 = a .+ t0 .* (b .- a), a .+ t1 .* (b .- a) + lines!(ax, [Point2f(p0), Point2f(p1)]; color=:darkorange, linewidth=3) + end + end + text!(ax, 0.5, 0.2; text="grid: image voxels; blue: image mask", align=(:left, :top), fontsize=12, color=:gray30) + fig +end + +# the parameters t in [0, 1] of the segment a -> b inside the box [lo, hi] +function clip_segment(a, b, lo, hi) + t0, t1 = 0.0, 1.0 + for k in 1:2 + d = b[k] - a[k] + if iszero(d) + (lo[k] <= a[k] <= hi[k]) || return 1.0, 0.0 + else + s0, s1 = (lo[k] - a[k]) / d, (hi[k] - a[k]) / d + t0, t1 = max(t0, min(s0, s1)), min(t1, max(s0, s1)) + end + end + t0, t1 +end