Nh/masks

benchmark/README.md

@@ -17,9 +17,13 @@ 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
+<cuda|rocm> [device] [problem] [projector]` times workgroup sizes and tiles of detector cells on one
+problem.
 
-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
@@ -133,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.
+

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

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 <cuda|rocm> [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

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())

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())

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
 

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

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);

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.

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.

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.