From c1a428906ac0382bfccb159d02403ec7ccaeadb5 Mon Sep 17 00:00:00 2001 From: Niklas Hackelberg Date: Tue, 3 Nov 2026 21:46:33 +0100 Subject: [PATCH 01/11] Geometry: rays at continuous positions on the detector cell_ray(g, cu, cv, k) gives the ray of view k that ends at position (cu, cv) on the detector, in cell units (integers are cell centers, half-integers cell edges), for cone beams with flat and arc detectors and parallel beams; ray(g, i) is the case of the cell center. Co-Authored-By: Claude Opus 5.5 --- docs/src/{api.md => api/geometries.md} | 0 src/geometry/arcconebeam.jl | 9 ++++++-- src/geometry/conebeam.jl | 6 ++++- src/geometry/geometry.jl | 15 +++++++++--- src/geometry/parallelbeam.jl | 7 +++++- test/geometry/geometry.jl | 32 ++++++++++++++++++++++++++ 6 files changed, 62 insertions(+), 7 deletions(-) rename docs/src/{api.md => api/geometries.md} (100%) diff --git a/docs/src/api.md b/docs/src/api/geometries.md similarity index 100% rename from docs/src/api.md rename to docs/src/api/geometries.md diff --git a/src/geometry/arcconebeam.jl b/src/geometry/arcconebeam.jl index 18b73f4..c09440e 100644 --- a/src/geometry/arcconebeam.jl +++ b/src/geometry/arcconebeam.jl @@ -64,6 +64,11 @@ storage_array(g::ArcConeBeam) = g.source @inline function ray(g::ArcConeBeam, i) iu, iv, k = pixel_index(i, g.nu, g.nv) + T = eltype(g) + cell_ray(g, T(iu), T(iv), k) +end + +@inline function cell_ray(g::ArcConeBeam, cu, cv, k) s = column(g.source, k) w = column(g.center, k) - s u = column(g.u, k) @@ -74,8 +79,8 @@ storage_array(g::ArcConeBeam) = g.source radial = w - h r = norm(radial) du = norm(u) - sg, cg = sincos((T(iu) - T(g.nu + 1) / 2) * du / r) - b = s + h + cg * radial + (sg * r / du) * u + (T(iv) - T(g.nv + 1) / 2) * v + sg, cg = sincos((cu - T(g.nu + 1) / 2) * du / r) + b = s + h + cg * radial + (sg * r / du) * u + (cv - T(g.nv + 1) / 2) * v s, b end diff --git a/src/geometry/conebeam.jl b/src/geometry/conebeam.jl index 5fc297f..7f5838c 100644 --- a/src/geometry/conebeam.jl +++ b/src/geometry/conebeam.jl @@ -44,8 +44,12 @@ storage_array(g::ConeBeam) = g.source @inline function ray(g::ConeBeam, i) iu, iv, k = pixel_index(i, g.nu, g.nv) - column(g.source, k), pixel_position(column(g.center, k), column(g.u, k), column(g.v, k), iu, iv, g.nu, g.nv) + T = eltype(g) + cell_ray(g, T(iu), T(iv), k) end +@inline cell_ray(g::ConeBeam, cu, cv, k) = + column(g.source, k), cell_position(column(g.center, k), column(g.u, k), column(g.v, k), cu, cv, g.nu, g.nv) + Adapt.adapt_structure(to, g::ConeBeam) = ConeBeam(adapt(to, g.source), adapt(to, g.center), adapt(to, g.u), adapt(to, g.v); nu=g.nu, nv=g.nv) diff --git a/src/geometry/geometry.jl b/src/geometry/geometry.jl index ae1a360..9a81026 100644 --- a/src/geometry/geometry.jl +++ b/src/geometry/geometry.jl @@ -113,12 +113,21 @@ datasize(g::AbstractDetectorGeometry) = (detectorsize(g)..., nviews(g)) ju + 1, jv + 1, k + 1 end -# the position of pixel (iu, iv) on a detector with center c and cell vectors u, v -@inline function pixel_position(c, u, v, iu, iv, nu, nv) +# the point at position (cu, cv) on a flat detector with center c and cell vectors u, v; positions +# are in cell units: integers are cell centers, half-integers cell edges +@inline function cell_position(c, u, v, cu, cv, nu, nv) T = eltype(c) - c + (T(iu) - T(nu + 1) / 2) * u + (T(iv) - T(nv + 1) / 2) * v + c + (cu - T(nu + 1) / 2) * u + (cv - T(nv + 1) / 2) * v end +# cell_ray(g, cu, cv, k) -> (a, b): the ray of view k of the detector geometry g that ends at +# position (cu, cv) on the detector, as two points like `ray` (cone beams: source and detector +# point; parallel beams: detector point and that point plus the direction). Positions are in cell +# units: integers are cell centers, so `cell_ray(g, iu, iv, k)` is the ray of cell (iu, iv), equal to +# `ray`; half-integers are cell edges, e.g. `cell_ray(g, iu + 1/2, iv, k)` runs through the midpoint +# of the edge between cells iu and iu + 1. Callable inside kernels. +function cell_ray end + function check_views(name, ms...) for m in ms size(m, 1) == 3 || throw(DimensionMismatch("$name: per-view vectors must be 3 × nviews, got $(size(m))")) diff --git a/src/geometry/parallelbeam.jl b/src/geometry/parallelbeam.jl index 21cb863..50d2a98 100644 --- a/src/geometry/parallelbeam.jl +++ b/src/geometry/parallelbeam.jl @@ -44,7 +44,12 @@ storage_array(g::ParallelBeam) = g.direction @inline function ray(g::ParallelBeam, i) iu, iv, k = pixel_index(i, g.nu, g.nv) - a = pixel_position(column(g.center, k), column(g.u, k), column(g.v, k), iu, iv, g.nu, g.nv) + T = eltype(g) + cell_ray(g, T(iu), T(iv), k) +end + +@inline function cell_ray(g::ParallelBeam, cu, cv, k) + a = cell_position(column(g.center, k), column(g.u, k), column(g.v, k), cu, cv, g.nu, g.nv) a, a + column(g.direction, k) end diff --git a/test/geometry/geometry.jl b/test/geometry/geometry.jl index 4944f1d..35903e4 100644 --- a/test/geometry/geometry.jl +++ b/test/geometry/geometry.jl @@ -106,6 +106,38 @@ end end end +# cell_ray gives rays at continuous cell coordinates (centers at integers, edges at half-integers): +# at the cell centers it equals ray; on flat detectors the end point is linear in the coordinates, +# on arc detectors the angle about the source is +@testitem "cell_ray" tags = [:geometry] setup = [Rays] begin + using CTOperators: cell_ray + using LinearAlgebra: norm, dot + angles = range(0, 2pi, 5)[1:4] + nu, nv = 6, 4 + flat = (circular_cone_beam(; angles, dso=500, dsd=900, nu, nv, du=1.5, dv=2.0, offset_u=0.25), + parallel_beam(; angles, nu, nv, du=0.8, dv=1.2)) + arc = circular_cone_beam(; angles, dso=500, dsd=900, nu, nv, du=1.5, dv=2.0, detector=:arc) + for g in (flat..., arc) + li = LinearIndices(datasize(g)) + for k in 1:4, iv in 1:nv, iu in 1:nu + @test all(cell_ray(g, Float64(iu), Float64(iv), k) .≈ ray(g, li[iu, iv, k])) + end + end + for g in flat, k in 1:4, (p, q, r) in (((2.0, 3.0), (3.0, 3.0), (2.25, 3.0)), ((2.0, 3.0), (2.0, 4.0), (2.0, 3.75))) + (a1, b1), (a2, b2), (a, b) = cell_ray(g, p..., k), cell_ray(g, q..., k), cell_ray(g, r..., k) + w = norm(collect(r .- p)) + @test isapprox(a, (1 - w) * a1 + w * a2) && isapprox(b, (1 - w) * b1 + w * b2) + end + for k in 1:4 + s, b1 = cell_ray(arc, 2.0, 3.0, k) + _, b2 = cell_ray(arc, 3.0, 3.0, k) + _, b = cell_ray(arc, 2.5, 3.0, k) + # same distance from the source in the plane of the arc, half the angle + @test norm(b - s) ≈ norm(b1 - s) ≈ norm(b2 - s) + @test dot(b - s, b1 - s) ≈ dot(b - s, b2 - s) + end +end + @testitem "subset" tags = [:geometry] setup = [Rays] begin views = [2, 5, 6] for g in (helical_cone_beam(; angles=range(0, 2pi, 8), pitch=10, dso=500, dsd=900, nu=4, nv=3, du=1.0), From c3c371a0ac55a2fce3a395b30064d4bc3019d7e0 Mon Sep 17 00:00:00 2001 From: Niklas Hackelberg Date: Tue, 3 Nov 2026 21:46:33 +0100 Subject: [PATCH 02/11] Docs: API reference split by domain One page each for images, projection geometries, operators (with masks, launch configuration and extensions), projectors and internals, instead of a single page. Co-Authored-By: Claude Opus 5.5 --- docs/make.jl | 6 +++- docs/src/api/geometries.md | 64 ++++++-------------------------------- docs/src/api/images.md | 13 ++++++++ docs/src/api/internals.md | 54 ++++++++++++++++++++++++++++++++ docs/src/api/operators.md | 46 +++++++++++++++++++++++++++ docs/src/api/projectors.md | 21 +++++++++++++ docs/src/index.md | 2 +- 7 files changed, 149 insertions(+), 57 deletions(-) create mode 100644 docs/src/api/images.md create mode 100644 docs/src/api/internals.md create mode 100644 docs/src/api/operators.md create mode 100644 docs/src/api/projectors.md diff --git a/docs/make.jl b/docs/make.jl index 2809b5b..38fe68b 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -55,7 +55,11 @@ makedocs(; "Workgroups and Ray Order on GPUs" => "generated/explanation/launch.md", ], "Reference" => [ - "API" => "api.md", + "Images" => "api/images.md", + "Projection Geometries" => "api/geometries.md", + "Operators" => "api/operators.md", + "Projectors" => "api/projectors.md", + "Internals" => "api/internals.md", ], ], ) diff --git a/docs/src/api/geometries.md b/docs/src/api/geometries.md index 74b9c1f..7e13c9d 100644 --- a/docs/src/api/geometries.md +++ b/docs/src/api/geometries.md @@ -1,18 +1,10 @@ -# API +# Projection Geometries -The sections follow the parts of the package: the image geometry, the projection geometries and -their rays, the generic operator interface, and one section per projector. +A projection geometry describes where the rays run and the shape of the projection data. Detector +geometries (cone beam with a flat or an arc detector, parallel beam) also describe the detector +cells; ray lists hold arbitrary rays. See [Set Up a Geometry](@ref) for examples. -## Image geometry - -Images are described by `ImageGeom` from [ImageGeoms.jl](https://github.com/JuliaImageRecon/ImageGeoms.jl), -re-exported by CTOperators.jl. Voxel `i` has its center at `((i - 1) - (n - 1)/2 - offset) * Δ` per -dimension; a 2D image is treated as a slice at `z = 0` that is 1 length unit thick. Lengths have -no fixed unit: image and projection geometry must use the same one. - -## Projection geometries - -### Interface +## Interface ```@autodocs Modules = [CTOperators] @@ -20,7 +12,7 @@ Pages = ["geometry/geometry.jl"] Private = false ``` -### Geometry types +## Geometry types ```@autodocs Modules = [CTOperators] @@ -28,50 +20,12 @@ Pages = ["geometry/raylist.jl", "geometry/conebeam.jl", "geometry/arcconebeam.jl Private = false ``` -### Trajectories - -```@autodocs -Modules = [CTOperators] -Pages = ["geometry/trajectories.jl"] -Private = false -``` - -## Operators +## Trajectories -### Generic interface - -```@docs -AbstractCTOperator -``` +Constructors for the common scans. ```@autodocs Modules = [CTOperators] -Pages = ["operators/operators.jl", "operators/launch.jl"] -Private = false -``` - -### Siddon - -```@autodocs -Modules = [CTOperators] -Pages = ["operators/siddon/siddon.jl", "operators/siddon/kernels.jl"] -Private = false -``` - -### Joseph - -```@autodocs -Modules = [CTOperators] -Pages = ["operators/joseph/joseph.jl", "operators/joseph/kernels.jl"] +Pages = ["geometry/trajectories.jl"] Private = false ``` - -## Internals - -For developers of operators: the kernel form of the image grid, clipping of rays and the -functions each operator implements. - -```@autodocs -Modules = [CTOperators] -Public = false -``` diff --git a/docs/src/api/images.md b/docs/src/api/images.md new file mode 100644 index 0000000..77ff0e6 --- /dev/null +++ b/docs/src/api/images.md @@ -0,0 +1,13 @@ +# Images + +Images are described by `ImageGeom` from [ImageGeoms.jl](https://github.com/JuliaImageRecon/ImageGeoms.jl), +re-exported by CTOperators.jl together with the mask types `MaskCircle`, `MaskEllipse` and +`MaskAllButEdge`. Voxel `i` has its center at `((i - 1) - (n - 1)/2 - offset) * Δ` per dimension; a +2D image is treated as a slice at `z = 0` that is 1 length unit thick. Lengths have no fixed unit: +image and projection geometry must use the same one. + +The mask of the `ImageGeom` is the image mask of a projector: the voxels that are reconstructed +(see [Use Masks](@ref)). + +How the projectors use the image geometry internally (`VoxelGrid`) is described under +[Internals](@ref). diff --git a/docs/src/api/internals.md b/docs/src/api/internals.md new file mode 100644 index 0000000..d9129c4 --- /dev/null +++ b/docs/src/api/internals.md @@ -0,0 +1,54 @@ +# Internals + +For developers of operators (see [Implement a CT Operator](@ref)): the functions each operator +implements, the masks a projector stores, the kernel form of the image grid and the clipping of +rays. None of these is exported. + +## Operator implementation + +```@autodocs +Modules = [CTOperators] +Pages = ["operators/operators.jl", "operators/launch.jl"] +Public = false +``` + +## Masks of a projector + +```@docs +CTOperators.ProjectorMasks +CTOperators.projector_masks +``` + +## Geometries + +```@autodocs +Modules = [CTOperators] +Pages = ["geometry/geometry.jl", "geometry/raylist.jl", "geometry/conebeam.jl", "geometry/arcconebeam.jl", + "geometry/parallelbeam.jl", "geometry/trajectories.jl"] +Public = false +``` + +## Image grid + +```@autodocs +Modules = [CTOperators] +Pages = ["image/image.jl", "image/voxelgrid.jl"] +Public = false +``` + +## Ray tracing + +```@autodocs +Modules = [CTOperators] +Pages = ["raytracing/raytracing.jl", "raytracing/clip.jl", "raytracing/actions.jl"] +Public = false +``` + +## Siddon and Joseph + +```@autodocs +Modules = [CTOperators] +Pages = ["operators/siddon/siddon.jl", "operators/siddon/kernels.jl", "operators/joseph/joseph.jl", + "operators/joseph/kernels.jl"] +Public = false +``` diff --git a/docs/src/api/operators.md b/docs/src/api/operators.md new file mode 100644 index 0000000..cfb7d0f --- /dev/null +++ b/docs/src/api/operators.md @@ -0,0 +1,46 @@ +# Operators + +The generic interface of all CT operators: application with `*` and `mul!`, adjoints, sizes, +row and column sums, subsets of views, moving to a device with `adapt`, and the launch +configuration of the kernels. + +## Generic interface + +```@docs +AbstractCTOperator +``` + +```@autodocs +Modules = [CTOperators] +Pages = ["operators/operators.jl"] +Private = false +``` + +## Image and ray masks + +The image mask comes from the `ImageGeom`, the ray mask from the keyword `raymask` of a projector +(see [Use Masks](@ref) and [Masks](@ref)). These functions are not exported. + +```@docs +CTOperators.image_mask +CTOperators.ray_mask +CTOperators.handles_raymask +``` + +## Launch configuration + +```@autodocs +Modules = [CTOperators] +Pages = ["operators/launch.jl"] +Private = false +``` + +## Extensions + +* With LinearOperators.jl loaded, `normalOperator(op)` of LinearOperatorCollection.jl gives the + normal operator `AᵀA` of a CT operator, e.g. for `AHA` of the solvers of + RegularizedLeastSquares.jl. +* With ChainRulesCore.jl loaded, the operators are differentiable (see + [Differentiate through a Projector](@ref)). +* With Sinograms.jl loaded, its geometries are converted by `projection_geometry` (see + [Use Geometries of Sinograms.jl](@ref)). diff --git a/docs/src/api/projectors.md b/docs/src/api/projectors.md new file mode 100644 index 0000000..d1e12d7 --- /dev/null +++ b/docs/src/api/projectors.md @@ -0,0 +1,21 @@ +# Projectors + +The projectors of the package. Both trace one ray per detector cell through the voxel grid and +have an exact adjoint; [Siddon's Method](@ref) and [Joseph's Method](@ref) explain how they +weight the voxels. + +## Siddon + +```@autodocs +Modules = [CTOperators] +Pages = ["operators/siddon/siddon.jl", "operators/siddon/kernels.jl"] +Private = false +``` + +## Joseph + +```@autodocs +Modules = [CTOperators] +Pages = ["operators/joseph/joseph.jl", "operators/joseph/kernels.jl"] +Private = false +``` diff --git a/docs/src/index.md b/docs/src/index.md index a547a3e..a35152c 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -59,7 +59,7 @@ y = op * rand(Float32, 128, 128, 128) # projections, size (256, 128, 360) 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. +- The reference lists all types and functions by domain: [Images](api/images.md), [Projection Geometries](api/geometries.md), [Operators](api/operators.md), [Projectors](api/projectors.md) and, for implementers, [Internals](api/internals.md). ## Related packages From 23a5d3d3a8807f7efabd4f6fc55fa35806942ec3 Mon Sep 17 00:00:00 2001 From: Niklas Hackelberg Date: Tue, 3 Nov 2026 22:39:02 +0100 Subject: [PATCH 03/11] Geometry: several rays per detector cell (MultiRay) MultiRay wraps a detector geometry with a point strategy (RegularPoints, JitteredPoints, RandomPoints, Points) and a weighting strategy (UniformWeights, GaussianWeights); Siddon and Joseph average the line integrals of the rays of each cell with normalized weights, with the exact adjoint. The number of rays per cell is a type parameter, so one ray per cell compiles to the plain projector. Random positions come from the hash pcg4d of Jarzynski & Olano (ported under its MIT license) of cell, view, ray and seed: the same in forward and back projection, on every device and in subsets. For 2D images the projectors keep the rays in the image plane. Co-Authored-By: Claude Opus 5.5 --- docs/src/api/geometries.md | 11 ++ docs/src/api/internals.md | 2 +- src/CTOperators.jl | 1 + src/geometry/geometry.jl | 2 + src/geometry/multiray.jl | 263 ++++++++++++++++++++++++++++++++ src/geometry/pcg4d.jl | 49 ++++++ src/operators/joseph/joseph.jl | 8 +- src/operators/joseph/kernels.jl | 13 +- src/operators/launch.jl | 6 +- src/operators/operators.jl | 6 +- src/operators/siddon/kernels.jl | 11 +- src/operators/siddon/siddon.jl | 8 +- src/raytracing/raytracing.jl | 1 + src/raytracing/subrays.jl | 32 ++++ test/geometry/multiray.jl | 56 +++++++ test/operators/multiray.jl | 119 +++++++++++++++ 16 files changed, 566 insertions(+), 22 deletions(-) create mode 100644 src/geometry/multiray.jl create mode 100644 src/geometry/pcg4d.jl create mode 100644 src/raytracing/subrays.jl create mode 100644 test/geometry/multiray.jl create mode 100644 test/operators/multiray.jl diff --git a/docs/src/api/geometries.md b/docs/src/api/geometries.md index 7e13c9d..2b0865c 100644 --- a/docs/src/api/geometries.md +++ b/docs/src/api/geometries.md @@ -20,6 +20,17 @@ Pages = ["geometry/raylist.jl", "geometry/conebeam.jl", "geometry/arcconebeam.jl Private = false ``` +## Several rays per cell + +A detector geometry wrapped with a point strategy (where in each cell the rays end) and a +weighting strategy. + +```@autodocs +Modules = [CTOperators] +Pages = ["geometry/multiray.jl"] +Private = false +``` + ## Trajectories Constructors for the common scans. diff --git a/docs/src/api/internals.md b/docs/src/api/internals.md index d9129c4..5534c2b 100644 --- a/docs/src/api/internals.md +++ b/docs/src/api/internals.md @@ -24,7 +24,7 @@ CTOperators.projector_masks ```@autodocs Modules = [CTOperators] Pages = ["geometry/geometry.jl", "geometry/raylist.jl", "geometry/conebeam.jl", "geometry/arcconebeam.jl", - "geometry/parallelbeam.jl", "geometry/trajectories.jl"] + "geometry/parallelbeam.jl", "geometry/trajectories.jl", "geometry/multiray.jl"] Public = false ``` diff --git a/src/CTOperators.jl b/src/CTOperators.jl index c692d69..15759d8 100644 --- a/src/CTOperators.jl +++ b/src/CTOperators.jl @@ -12,6 +12,7 @@ export AbstractCTOperator export ImageGeom, MaskCircle, MaskEllipse, MaskAllButEdge export AbstractProjectionGeometry, AbstractDetectorGeometry, RaySegment, RayLine export RayList, ConeBeam, ArcConeBeam, ParallelBeam +export MultiRay, RegularPoints, JitteredPoints, RandomPoints, Points, UniformWeights, GaussianWeights export circular_cone_beam, helical_cone_beam, parallel_beam export datasize, nrays, ray, raytype, nviews, detectorsize, subset, projection_geometry export AbstractProjector, AdjointCTOperator, SiddonProjector, JosephProjector, LaunchConfig diff --git a/src/geometry/geometry.jl b/src/geometry/geometry.jl index 9a81026..f4918d2 100644 --- a/src/geometry/geometry.jl +++ b/src/geometry/geometry.jl @@ -170,3 +170,5 @@ include("conebeam.jl") include("arcconebeam.jl") include("parallelbeam.jl") include("trajectories.jl") +include("pcg4d.jl") +include("multiray.jl") diff --git a/src/geometry/multiray.jl b/src/geometry/multiray.jl new file mode 100644 index 0000000..46d48ad --- /dev/null +++ b/src/geometry/multiray.jl @@ -0,0 +1,263 @@ +# Several rays per detector cell: a detector geometry wrapped with a point strategy (where in the +# cell the subrays end) and a weighting strategy. Offsets are in cell units: (0, 0) is the cell +# center, +-1/2 the cell edges. The number of subrays K is a type parameter of the point strategy, +# so the loop over the subrays has a length known to the compiler. + +""" + RegularPoints(su, sv=su) + +The centers of `su × sv` subcells of each detector cell: `su` along `u`, `sv` along `v`. Point +strategy of [`MultiRay`](@ref). +""" +struct RegularPoints{K} + su::Int + sv::Int +end + +""" + JitteredPoints(su, sv=su; seed=0) + +The points of [`RegularPoints`](@ref), each moved to a random position within its subcell +("jittered sampling", R. L. Cook, "Stochastic sampling in computer graphics", ACM Transactions on +Graphics 5(1), 51-72, 1986). The positions are fixed by `seed`, the detector cell, the view and the +subray, so they are the same in the forward and the back projection, on every device and in +`subset`. Point strategy of [`MultiRay`](@ref). +""" +struct JitteredPoints{K} + su::Int + sv::Int + seed::UInt32 +end + +""" + RandomPoints(n; seed=0) + +`n` random points in each detector cell, fixed like those of [`JitteredPoints`](@ref). Point +strategy of [`MultiRay`](@ref). +""" +struct RandomPoints{K} + seed::UInt32 +end + +""" + Points(offsets) + +Given points, the same in every detector cell: column `k` of the `2 × K` matrix `offsets` is the +position of subray `k` in cell units, `(0, 0)` the cell center and `±1/2` the cell edges. Points +outside the cell are allowed. Point strategy of [`MultiRay`](@ref). +""" +struct Points{K,T} + offsets::SVector{K,SVector{2,T}} +end + +function check_subcells(su, sv) + su > 0 && sv > 0 || throw(ArgumentError("numbers of subcells must be positive, got ($su, $sv)")) +end + +function RegularPoints(su::Integer, sv::Integer=su) + check_subcells(su, sv) + RegularPoints{Int(su * sv)}(su, sv) +end + +function JitteredPoints(su::Integer, sv::Integer=su; seed::Integer=0) + check_subcells(su, sv) + JitteredPoints{Int(su * sv)}(su, sv, UInt32(seed)) +end + +function RandomPoints(n::Integer; seed::Integer=0) + n > 0 || throw(ArgumentError("number of points must be positive, got $n")) + RandomPoints{Int(n)}(UInt32(seed)) +end + +function Points(offsets::AbstractMatrix{<:Real}) + size(offsets, 1) == 2 && size(offsets, 2) > 0 || + throw(DimensionMismatch("offsets must be 2 × K with K > 0, got $(size(offsets))")) + K = size(offsets, 2) + T = float(eltype(offsets)) + Points{K,T}(SVector{K}(ntuple(k -> SVector{2,T}(offsets[1, k], offsets[2, k]), K))) +end + +npoints(::Union{RegularPoints{K},JitteredPoints{K},RandomPoints{K},Points{K}}) where {K} = K + +# the points in element type T (only given points store numbers) +points_eltype(p, ::Type{T}) where {T} = p +points_eltype(p::Points{K}, ::Type{T}) where {K,T} = Points{K,T}(map(o -> SVector{2,T}(o), p.offsets)) + +# The offset (du, dv) of subray k in the cell (iu, iv) of view `view` (numbered as in the full +# geometry, so that subsets draw the same random positions). +@inline function point_offset(p::RegularPoints, ::Type{T}, iu, iv, view, k) where {T} + ju, jv = rem(k - 1, p.su), div(k - 1, p.su) + half = T(1 // 2) + (ju + half) / p.su - half, (jv + half) / p.sv - half +end + +@inline function point_offset(p::JitteredPoints, ::Type{T}, iu, iv, view, k) where {T} + ju, jv = rem(k - 1, p.su), div(k - 1, p.su) + ru, rv = random_pair(T, p.seed, iu, iv, view, k) + half = T(1 // 2) + (ju + ru) / p.su - half, (jv + rv) / p.sv - half +end + +@inline function point_offset(p::RandomPoints, ::Type{T}, iu, iv, view, k) where {T} + ru, rv = random_pair(T, p.seed, iu, iv, view, k) + half = T(1 // 2) + ru - half, rv - half +end + +@inline function point_offset(p::Points, ::Type{T}, iu, iv, view, k) where {T} + o = @inbounds p.offsets[k] + T(o[1]), T(o[2]) +end + +# two numbers in [0, 1) from the hash of (cell, view, subray, seed); the cell as one number +# (iu + 65536 iv) keeps the four inputs distinct for detectors of up to 65536 cells per row +@inline function random_pair(::Type{T}, seed, iu, iv, view, k) where {T} + cell = (iu % UInt32) + ((iv % UInt32) << 16) + hu, hv, _, _ = pcg4d(cell, view % UInt32, k % UInt32, seed) + unit_interval(T, hu), unit_interval(T, hv) +end + +""" + UniformWeights() + +All subrays of a detector cell count equally. Weighting strategy of [`MultiRay`](@ref). +""" +struct UniformWeights end + +""" + GaussianWeights(sigma_u, sigma_v=sigma_u) + +Subrays weighted by a Gaussian of the offset from the cell center, with standard deviations +`sigma_u`, `sigma_v` in cell units: a detector response that falls off towards the cell edges. +Weighting strategy of [`MultiRay`](@ref). +""" +struct GaussianWeights{T} + sigma_u::T + sigma_v::T + function GaussianWeights{T}(sigma_u, sigma_v) where {T} + sigma_u > 0 && sigma_v > 0 || throw(ArgumentError("standard deviations must be positive, got ($sigma_u, $sigma_v)")) + new{T}(sigma_u, sigma_v) + end +end + +GaussianWeights(sigma_u::Real, sigma_v::Real=sigma_u) = + GaussianWeights{float(promote_type(typeof(sigma_u), typeof(sigma_v)))}(sigma_u, sigma_v) + +weights_eltype(w, ::Type{T}) where {T} = w +weights_eltype(w::GaussianWeights, ::Type{T}) where {T} = GaussianWeights{T}(w.sigma_u, w.sigma_v) + +# the (unnormalized) weight of a subray at offset (du, dv) +@inline subray_weight(::UniformWeights, du, dv) = one(du) +@inline subray_weight(w::GaussianWeights, du, dv) = exp(-((du / w.sigma_u)^2 + (dv / w.sigma_v)^2) / 2) + +""" + MultiRay(geometry, points; weights=UniformWeights()) + +The detector geometry `geometry` (`ConeBeam`, `ArcConeBeam`, `ParallelBeam`) with several rays per +detector cell. The rays of a cell end at the positions given by the point strategy `points` +([`RegularPoints`](@ref), [`JitteredPoints`](@ref), [`RandomPoints`](@ref), [`Points`](@ref)) and +are averaged with the weights of `weights` ([`UniformWeights`](@ref), +[`GaussianWeights`](@ref)), normalized to sum 1 per cell. A ray-driven projector built on it +computes + + y[i] = sum_k w_k * (line integral along subray k of data element i) + +and the exact adjoint. This averages the line integrals over the cell (closer to what a detector +cell measures than one line through its center) and removes the gaps between rays that occur when +detector cells are large against the voxels. The data has the size of `geometry`. + +For 2D images (one voxel layer along `z`) the projectors keep the subrays in the image plane: they +are spread along `u` only, and point strategies with several subcells along `v` are rejected. + +A `MultiRay` has no single ray per data element (no [`ray`](@ref)); projectors that trace only one +ray per cell reject it. `subset` and `adapt` keep the points and weights; random positions stay the +same for the views of a subset. + +# Example +```julia +g = MultiRay(circular_cone_beam(; angles, dso, dsd, nu, nv, du), JitteredPoints(2, 2; seed=1)) +op = JosephProjector(image, g) +``` +""" +# Whether the subrays of a cell move along both detector axes or only along u. For 2D images (one +# voxel layer along z) the rays must stay in the image plane, which the projectors ensure with +# `in_plane`. +struct FullCell end +struct InPlane end + +@inline restrict(::FullCell, du, dv) = (du, dv) +@inline restrict(::InPlane, du, dv) = (du, zero(dv)) + +struct MultiRay{T,G<:AbstractDetectorGeometry{T},P,W,V,S} <: AbstractProjectionGeometry{T} + geometry::G + points::P + weights::W + # nothing, or the numbers of the views in the geometry this one is a subset of, stored like the + # other arrays of the geometry (element type T: exact for up to 2^24 views in Float32) + views::V + span::S # FullCell() or InPlane() +end + +function MultiRay(geometry, points; weights=UniformWeights()) + g = projection_geometry(geometry) + g isa AbstractDetectorGeometry || + throw(ArgumentError("MultiRay needs a detector geometry (cells to sample), got $(nameof(typeof(g)))")) + multiray(g, points, weights, nothing, FullCell()) +end + +function multiray(g::AbstractDetectorGeometry{T}, points, weights, views, span) where {T} + P, W = points_eltype(points, T), weights_eltype(weights, T) + MultiRay{T,typeof(g),typeof(P),typeof(W),typeof(views),typeof(span)}(g, P, W, views, span) +end + +nsubrays(g::MultiRay) = npoints(g.points) +datasize(g::MultiRay) = datasize(g.geometry) +raytype(g::MultiRay) = raytype(g.geometry) +nviews(g::MultiRay) = nviews(g.geometry) +detectorsize(g::MultiRay) = detectorsize(g.geometry) +storage_array(g::MultiRay) = storage_array(g.geometry) + +function subset(g::MultiRay{T}, views) where {T} + numbers = isnothing(g.views) ? collect(views) : Array(g.views)[views] + inner = subset(g.geometry, views) + multiray(inner, g.points, g.weights, copyto!(geometry_array(inner, length(numbers)), T.(numbers)), g.span) +end + +Adapt.adapt_structure(to, g::MultiRay) = + multiray(adapt(to, g.geometry), g.points, g.weights, adapt(to, g.views), g.span) + +# the geometry with its subrays in the plane v = 0, for 2D images (plain geometries are unchanged) +in_plane(g) = g +function in_plane(g::MultiRay) + check_in_plane(g.points) + multiray(g.geometry, g.points, g.weights, g.views, InPlane()) +end + +# subcells along v would all fall onto v = 0 and trace the same rays several times +check_in_plane(p::Union{RegularPoints,JitteredPoints}) = + p.sv == 1 || throw(ArgumentError("a 2D image has one row of rays per view: use $(nameof(typeof(p)))(su, 1)")) +check_in_plane(p) = nothing + +function Base.show(io::IO, ::MIME"text/plain", g::MultiRay) + print(io, "MultiRay with ", nsubrays(g), " rays per cell (", nameof(typeof(g.points)), ", ", + nameof(typeof(g.weights)), ") of ") + show(io, MIME("text/plain"), g.geometry) +end + +# the view number in the full geometry +@inline full_view(::Nothing, view) = view +@inline full_view(views, view) = @inbounds unsafe_trunc(Int, views[view]) + +""" + subray(g::MultiRay, i, k) -> (a, b, w) + +Subray `k` of data element `i` as two points like [`ray`](@ref), and its weight (not normalized). +Callable inside kernels. +""" +@inline function subray(g::MultiRay{T}, i, k) where {T} + nu, nv = detectorsize(g.geometry) + iu, iv, view = pixel_index(i, nu, nv) + du, dv = restrict(g.span, point_offset(g.points, T, iu, iv, full_view(g.views, view), k)...) + a, b = cell_ray(g.geometry, T(iu) + du, T(iv) + dv, view) + a, b, subray_weight(g.weights, du, dv) +end diff --git a/src/geometry/pcg4d.jl b/src/geometry/pcg4d.jl new file mode 100644 index 0000000..62e8624 --- /dev/null +++ b/src/geometry/pcg4d.jl @@ -0,0 +1,49 @@ +# The hash pcg4d of M. Jarzynski and M. Olano, "Hash Functions for GPU Rendering", Journal of +# Computer Graphics Techniques 9(3), 21-38, 2020, ported from their HLSL +# (https://github.com/markjarzynski/PCG3D, pcg4d.hlsl), which carries this notice: +# +# Copyright 2020 Mark Jarzynski & Marc Olano +# +# Permission is hereby granted, free of charge, to any person obtaining a copy of this software and +# associated documentation files (the "Software"), to deal in the Software without restriction, +# including without limitation the rights to use, copy, modify, merge, publish, distribute, +# sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all copies or +# substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT +# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, +# DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT +# OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + +# Four UInt32 inputs to four UInt32 outputs that look independent even for neighboring inputs; +# random numbers computed from indices, the same on every device. UInt32 arithmetic wraps. +@inline function pcg4d(x::UInt32, y::UInt32, z::UInt32, w::UInt32) + x = x * 0x0019660d + 0x3c6ef35f + y = y * 0x0019660d + 0x3c6ef35f + z = z * 0x0019660d + 0x3c6ef35f + w = w * 0x0019660d + 0x3c6ef35f + + x += y * w + y += z * x + z += x * y + w += y * z + + x = xor(x, x >> 16) + y = xor(y, y >> 16) + z = xor(z, z >> 16) + w = xor(w, w >> 16) + + x += y * w + y += z * x + z += x * y + w += y * z + + x, y, z, w +end + +# a number in [0, 1) from the upper 24 bits of a hash value (exact in Float32) +@inline unit_interval(::Type{T}, h::UInt32) where {T} = T(h >> 8) * T(0x1p-24) diff --git a/src/operators/joseph/joseph.jl b/src/operators/joseph/joseph.jl index 4821fb5..417a4ab 100644 --- a/src/operators/joseph/joseph.jl +++ b/src/operators/joseph/joseph.jl @@ -67,10 +67,12 @@ 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 +# from a geometry already in the operator's element type and on its device; for 2D images the +# subrays of a MultiRay stay in the image plane 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, ProjectorMasks(image, geometry, raymask)) + g = length(image.dims) == 2 ? in_plane(geometry) : geometry + backend = KA.get_backend(storage_array(g)) + JosephProjector(image, VoxelGrid{T}(image), g, backend, launch, ProjectorMasks(image, g, raymask)) end geometry(op::JosephProjector) = op.geometry diff --git a/src/operators/joseph/kernels.jl b/src/operators/joseph/kernels.jl index 7b568a8..90011dc 100644 --- a/src/operators/joseph/kernels.jl +++ b/src/operators/joseph/kernels.jl @@ -73,17 +73,20 @@ end lo, hi end -# The ray i of the geometry through the grid; returns the accumulated value of the action. -@inline function joseph_ray(action, grid::VoxelGrid{T}, geometry, i) where {T} - a, b = ray(geometry, i) +# The rays of data element i of the geometry through the grid (one, or the subrays of a MultiRay); +# returns the accumulated value of the action. +@inline joseph_ray(action, grid, geometry, i) = trace_element(joseph_segment, action, grid, geometry, i) + +# one ray from a to b of type rt +@inline function joseph_segment(action, grid::VoxelGrid{T}, a, b, rt) where {T} p = to_index(grid, a) 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, k) - amin, amax = clip(p, d, lo, hi, raytype(geometry)) + amin, amax = clip(p, d, lo, hi, rt) amin < amax || return zero(T) - smin, smax = segment_range(raytype(geometry), T) + smin, smax = segment_range(rt, T) joseph_traverse(action, grid.dims, p, d, k, norm(b - a), amin, amax, smin, smax) end diff --git a/src/operators/launch.jl b/src/operators/launch.jl index 7ff02fd..ba70653 100644 --- a/src/operators/launch.jl +++ b/src/operators/launch.jl @@ -52,7 +52,7 @@ launch_config(::Nothing, geometry, backend) = default_launch(geometry, backend) # (KA.GPU is no help: since KernelAbstractions 0.10 it is the same as KA.Backend, the CPU included) function default_launch(geometry, backend) backend isa KA.CPU && return LaunchConfig() - LaunchConfig(; workgroupsize=64, tile=geometry isa AbstractDetectorGeometry ? (8, 8) : nothing) + LaunchConfig(; workgroupsize=64, tile=geometry isa Union{AbstractDetectorGeometry,MultiRay} ? (8, 8) : nothing) end # the rays in data order @@ -67,9 +67,9 @@ struct TileOrder tv::Int end -# tiles only for detector geometries; ray lists keep their order +# tiles only for detector geometries (also with several rays per cell); ray lists keep their order ray_order(launch::LaunchConfig, geometry) = DataOrder() -ray_order(launch::LaunchConfig, geometry::AbstractDetectorGeometry) = +ray_order(launch::LaunchConfig, geometry::Union{AbstractDetectorGeometry,MultiRay}) = isnothing(launch.tile) ? DataOrder() : TileOrder(detectorsize(geometry)..., launch.tile...) nthreads(::DataOrder, geometry) = nrays(geometry) diff --git a/src/operators/operators.jl b/src/operators/operators.jl index bd93d15..1c56e6e 100644 --- a/src/operators/operators.jl +++ b/src/operators/operators.jl @@ -143,11 +143,11 @@ 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`, and the +The projector of the views `views` of `op` (detector geometries, also with several rays per cell; +a vector or range of view indices), for ordered-subset methods. It shares the image geometry and options with `op`, and the views of its ray mask. """ -function subset(op::AbstractProjector{T,<:AbstractDetectorGeometry}, views) where {T} +function subset(op::AbstractProjector{T,<:Union{AbstractDetectorGeometry,MultiRay}}, 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]) diff --git a/src/operators/siddon/kernels.jl b/src/operators/siddon/kernels.jl index 926d127..ad8dd87 100644 --- a/src/operators/siddon/kernels.jl +++ b/src/operators/siddon/kernels.jl @@ -33,12 +33,15 @@ end end end -# The ray i of the geometry through the grid; returns the accumulated value of the action. -@inline function siddon_ray(action, grid::VoxelGrid{T}, geometry, i) where {T} - a, b = ray(geometry, i) +# The rays of data element i of the geometry through the grid (one, or the subrays of a MultiRay); +# returns the accumulated value of the action. +@inline siddon_ray(action, grid, geometry, i) = trace_element(siddon_segment, action, grid, geometry, i) + +# one ray from a to b of type rt +@inline function siddon_segment(action, grid::VoxelGrid{T}, a, b, rt) where {T} p = to_index(grid, a) d = (b - a) ./ grid.deltas - amin, amax = clip(p, d, grid.lo, grid.hi, raytype(geometry)) # the box of the image mask + amin, amax = clip(p, d, grid.lo, grid.hi, rt) # 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 a4b5c54..a46560c 100644 --- a/src/operators/siddon/siddon.jl +++ b/src/operators/siddon/siddon.jl @@ -57,10 +57,12 @@ 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 +# from a geometry already in the operator's element type and on its device; for 2D images the +# subrays of a MultiRay stay in the image plane 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, ProjectorMasks(image, geometry, raymask)) + g = length(image.dims) == 2 ? in_plane(geometry) : geometry + backend = KA.get_backend(storage_array(g)) + SiddonProjector(image, VoxelGrid{T}(image), g, backend, launch, ProjectorMasks(image, g, raymask)) end geometry(op::SiddonProjector) = op.geometry diff --git a/src/raytracing/raytracing.jl b/src/raytracing/raytracing.jl index 2049454..c08bf1c 100644 --- a/src/raytracing/raytracing.jl +++ b/src/raytracing/raytracing.jl @@ -1,2 +1,3 @@ include("clip.jl") include("actions.jl") +include("subrays.jl") diff --git a/src/raytracing/subrays.jl b/src/raytracing/subrays.jl new file mode 100644 index 0000000..f863b3d --- /dev/null +++ b/src/raytracing/subrays.jl @@ -0,0 +1,32 @@ +# The rays of data element i for the ray-driven methods: one ray for plain geometries, the weighted +# subrays of a MultiRay. `segment(action, grid, a, b, raytype)` traces one ray from a to b and +# returns the value of the action (Gather: the weighted sum of the image). + +@inline trace_element(segment::F, action, grid, geometry, i) where {F} = + segment(action, grid, ray(geometry, i)..., raytype(geometry)) + +# forward: the weighted mean of the subrays' line integrals +@inline function trace_element(segment::F, action::Gather, grid::VoxelGrid{T}, g::MultiRay, i) where {F,T} + acc, wsum = zero(T), zero(T) + for k in 1:nsubrays(g) + a, b, w = subray(g, i, k) + acc = muladd(w, segment(action, grid, a, b, raytype(g)), acc) + wsum += w + end + acc / wsum +end + +# back: every subray adds its share w / sum(w) of the value; the weights are summed first (cheap +# next to tracing), so that the shares are known before any subray is traced +@inline function trace_element(segment::F, action::Scatter, grid::VoxelGrid{T}, g::MultiRay, i) where {F,T} + wsum = zero(T) + for k in 1:nsubrays(g) + _, _, w = subray(g, i, k) + wsum += w + end + for k in 1:nsubrays(g) + a, b, w = subray(g, i, k) + segment(Scatter(action.image, action.value * (w / wsum)), grid, a, b, raytype(g)) + end + zero(T) +end diff --git a/test/geometry/multiray.jl b/test/geometry/multiray.jl new file mode 100644 index 0000000..7b0397f --- /dev/null +++ b/test/geometry/multiray.jl @@ -0,0 +1,56 @@ +# MultiRay (src/geometry/multiray.jl): the subrays of a detector cell end inside the cell at the +# positions of the point strategy; random positions depend only on seed, cell, view (numbered as in +# the full geometry) and subray, so subsets and other element types draw the same ones. + +@testitem "MultiRay: points and subrays" tags = [:geometry, :multiray] begin + using CTOperators: subray, nsubrays, point_offset, pixel_index + g = circular_cone_beam(; angles=range(0, 2pi, 6)[1:5], dso=50, dsd=100, nu=7, nv=3, du=1.5, dv=2.0) + # one subray at the cell center is the ray of the cell, with weight 1 + m = MultiRay(g, RegularPoints(1)) + @test nsubrays(m) == 1 && datasize(m) == datasize(g) + @test all(subray(m, i, 1) == (ray(g, i)..., 1.0) for i in 1:nrays(g)) + # offsets: regular points at the subcell centers, jittered ones inside their subcell, random ones + # inside the cell + su, sv = 3, 2 + for i in 1:nrays(g) + iu, iv, view = pixel_index(i, 7, 3) + for k in 1:su*sv + ju, jv = rem(k - 1, su), div(k - 1, su) + @test all(isapprox.(point_offset(RegularPoints(su, sv), Float64, iu, iv, view, k), ((ju + 0.5) / su - 0.5, (jv + 0.5) / sv - 0.5))) + du, dv = point_offset(JitteredPoints(su, sv; seed=7), Float64, iu, iv, view, k) + @test ju / su - 0.5 <= du < (ju + 1) / su - 0.5 && jv / sv - 0.5 <= dv < (jv + 1) / sv - 0.5 + du, dv = point_offset(RandomPoints(su * sv; seed=7), Float64, iu, iv, view, k) + @test -0.5 <= du < 0.5 && -0.5 <= dv < 0.5 + end + end + # different seeds, cells and subrays give different positions + offsets = [point_offset(RandomPoints(4; seed=s), Float64, iu, 1, 1, k) for s in 0:1, iu in 1:7, k in 1:4] + @test allunique(offsets) + # given points; Gaussian weights fall off from the center + m = MultiRay(g, Points([0.0 0.25; 0.0 0.0]); weights=GaussianWeights(0.5)) + @test subray(m, 1, 1)[3] == 1 && isapprox(subray(m, 1, 2)[3], exp(-0.125)) +end + +@testitem "MultiRay: subset, adapt, errors" tags = [:geometry, :multiray] begin + using Adapt + using CTOperators: subray + g = helical_cone_beam(; angles=range(0, 4pi, 9), pitch=10, dso=50, dsd=100, nu=5, nv=4, du=1.5) + m = MultiRay(g, JitteredPoints(2, 3; seed=11); weights=GaussianWeights(0.3, 0.6)) + # the subrays of a subset (also of a subset of a subset) are those of the same views of m + views = [2, 5, 9] + s = subset(subset(m, [1, 2, 5, 8, 9]), [2, 3, 5]) + li, ls = LinearIndices(datasize(m)), LinearIndices(datasize(s)) + @test all(all(isapprox.(subray(s, ls[iu, iv, j], k), subray(m, li[iu, iv, views[j]], k))) + for iu in 1:5, iv in 1:4, j in 1:3, k in 1:6) + # adapt to Float32 keeps points, weights and the view numbers + s32 = adapt(Array{Float32}, s) + @test eltype(s32) == Float32 && s32.weights isa GaussianWeights{Float32} + @test all(isapprox.(subray(s32, 7, 4), subray(s, 7, 4); rtol=1e-5)) + # ray lists have no cells; a MultiRay has no single ray per element + @test_throws ArgumentError MultiRay(RayList(g), RegularPoints(2)) + @test_throws MethodError ray(m, 1) + @test_throws ArgumentError RegularPoints(0, 2) + @test_throws ArgumentError RandomPoints(0) + @test_throws ArgumentError GaussianWeights(0.5, -1.0) + @test_throws DimensionMismatch Points(zeros(3, 2)) +end diff --git a/test/operators/multiray.jl b/test/operators/multiray.jl new file mode 100644 index 0000000..c630d66 --- /dev/null +++ b/test/operators/multiray.jl @@ -0,0 +1,119 @@ +# Projectors on a MultiRay average the line integrals of the subrays of each detector cell with +# normalized weights, with the exact adjoint, on every device (random positions from a hash, so the +# same everywhere). + +@testsnippet MultiRayProblem begin + using LinearAlgebra + # a small volume, detector cells larger than the voxels + image = ImageGeom(; dims=(12, 10, 6), deltas=(1.0, 1.1, 0.9), offsets=(0.2, -0.1, 0.3)) + cone(; nu=6, nv=3, du=2.6, dv=2.4) = + circular_cone_beam(; angles=range(0, 2pi, 6)[1:5], dso=40, dsd=70, nu, nv, du, dv) + parallel(; nu=6, nv=3, du=2.2, dv=2.0) = parallel_beam(; angles=range(0, pi, 6)[1:5], nu, nv, du, dv) + arc(; nu=6, nv=3, du=2.6, dv=2.4) = + circular_cone_beam(; angles=range(0, 2pi, 6)[1:5], dso=40, dsd=70, nu, nv, du, dv, detector=:arc) + # the mean over blocks of su x sv cells of data of a detector of su nu x sv nv cells + function block_mean(y, su, sv) + nu, nv, nk = div(size(y, 1), su), div(size(y, 2), sv), size(y, 3) + dropdims(sum(reshape(y, su, nu, sv, nv, nk); dims=(1, 3)); dims=(1, 3)) ./ (su * sv) + end +end + +@testitem "MultiRay: regular points are a finer detector" tags = [:operators, :multiray] setup = [MultiRayProblem] begin + # the subcell centers of RegularPoints(su, sv) are the cell centers of a detector with su x sv + # times the cells, so with uniform weights the projection is the block mean of that detector's + # (on the arc detector the angle about the source is linear in the position, so also there) + su, sv = 3, 2 + for (coarse, fine) in ((cone(), cone(; nu=18, nv=6, du=2.6 / 3, dv=1.2)), + (arc(), arc(; nu=18, nv=6, du=2.6 / 3, dv=1.2)), + (parallel(), parallel(; nu=18, nv=6, du=2.2 / 3, dv=1.0))), + P in (SiddonProjector, JosephProjector) + op = P(image, MultiRay(coarse, RegularPoints(su, sv)); eltype=Float64) + x = rand(image.dims...) + @test isapprox(op * x, block_mean(P(image, fine; eltype=Float64) * x, su, sv)) + end +end + +@testitem "MultiRay: weights" tags = [:operators, :multiray] setup = [MultiRayProblem] begin + x = rand(image.dims...) + for P in (SiddonProjector, JosephProjector) + plain = P(image, cone(); eltype=Float64) + # one point at the center (any weight) is the plain projector + @test isapprox(P(image, MultiRay(cone(), Points([0.0; 0.0;;]); weights=GaussianWeights(0.2)); eltype=Float64) * x, plain * x) + # weights are normalized: equal points with Gaussian weights give the plain projection + @test isapprox(P(image, MultiRay(cone(), Points([0.0 0.0; 0.0 0.0]); weights=GaussianWeights(0.3)); eltype=Float64) * x, plain * x) + # Gaussian weights with a large width are uniform + uniform = P(image, MultiRay(cone(), JitteredPoints(2, 2; seed=5)); eltype=Float64) + wide = P(image, MultiRay(cone(), JitteredPoints(2, 2; seed=5); weights=GaussianWeights(1e4)); eltype=Float64) + @test isapprox(wide * x, uniform * x) + end +end + +@testitem "MultiRay: adjoint and devices" tags = [:operators, :multiray, :gpu] setup = [Backends, MultiRayProblem] begin + using Adapt + for P in (SiddonProjector, JosephProjector), + (g, pts) in ((cone(), JitteredPoints(2, 3; seed=1)), (arc(), RegularPoints(2)), (parallel(), RandomPoints(5; seed=2))), + arrayType in arrayTypes, T in floatTypes + @testset "$P, $(nameof(typeof(pts))), $arrayType{$T}" begin + m = MultiRay(g, pts; weights=GaussianWeights(0.4)) + op = P(image, m; eltype=T, storage=arrayType) + ref = P(image, m; eltype=Float64) + x, y = rand(image.dims...), rand(datasize(g)...) + # the same subrays on every device and in both precisions + @test isapprox(Array(op * adapt(arrayType{T}, x)), ref * x; rtol=sqrt(eps(T))) + @test isapprox(Array(adjoint(op) * adapt(arrayType{T}, y)), adjoint(ref) * y; rtol=sqrt(eps(T))) + 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 "MultiRay: subset and ray mask" tags = [:operators, :multiray] setup = [MultiRayProblem] begin + using CTOperators: ray_mask + g = cone() + raymask = rand(datasize(g)...) .> 0.3 + op = JosephProjector(image, MultiRay(g, JitteredPoints(2; seed=9)); eltype=Float64, raymask) + x = rand(image.dims...) + y = op * x + # a subset of views gives the same rows (same random positions) + @test isapprox(subset(op, [2, 4]) * x, y[:, :, [2, 4]]) + # masked data elements skip all their subrays + @test all(iszero, y[.!raymask]) + # the default launch tiles the detector also with several rays per cell; results do not depend on it + tiled = JosephProjector(image, MultiRay(g, JitteredPoints(2; seed=9)); eltype=Float64, raymask, + launch=LaunchConfig(; tile=(4, 2))) + @test isapprox(tiled * x, y) +end + +# A 2D image is one voxel layer along z: the projectors keep the subrays in the image plane (offsets +# along u only), for flat and arc detectors; 3D images keep offsets along v +@testitem "MultiRay: 2D images" tags = [:operators, :multiray] setup = [MultiRayProblem] begin + using CTOperators: subray, cell_ray, pixel_index, point_offset, InPlane, FullCell + image2 = ImageGeom(; dims=(24, 20), deltas=(1.0, 1.0)) + fan(detector) = + circular_cone_beam(; angles=range(0, 2pi, 9)[1:8], dso=40, dsd=70, nu=10, nv=1, du=3.0, dv=3.0, detector) + pts = JitteredPoints(4, 1; seed=3) + for g in (fan(:flat), fan(:arc), parallel_beam(; angles=range(0, pi, 9)[1:8], nu=10, nv=1, du=3.0, dv=3.0)) + m = MultiRay(g, pts) + s = geometry(SiddonProjector(image2, m; eltype=Float64)) + @test s.span isa InPlane && geometry(SiddonProjector(image, m; eltype=Float64)).span isa FullCell + # every subray ends at its jittered position along u, in the row of the cell centers + for i in 1:nrays(g), k in 1:4 + iu, iv, view = pixel_index(i, 10, 1) + du, _ = point_offset(pts, Float64, iu, iv, view, k) + a, b, _ = subray(s, i, k) + @test all(isapprox.((a, b), cell_ray(g, iu + du, 1.0, view))) + end + end + # the projection is the mean of the projections of the subrays as single rays (a ray list): the + # in-plane subrays see the image like any ray in the plane + x = rand(24, 20) + for g in (fan(:flat), fan(:arc)), P in (SiddonProjector, JosephProjector) + s = geometry(P(image2, MultiRay(g, pts); eltype=Float64)) + rays = [subray(s, i, k) for k in 1:4, i in 1:nrays(g)] + list = RayList(stack(vec([r[1] for r in rays])), stack(vec([r[2] for r in rays])); datasize=(4, nrays(g))) + expected = reshape(sum(P(image2, list; eltype=Float64) * x; dims=1) ./ 4, datasize(g)) + @test isapprox(P(image2, MultiRay(g, pts); eltype=Float64) * x, expected) + end + # subcells along v would collapse onto the plane + @test_throws ArgumentError SiddonProjector(image2, MultiRay(fan(:flat), RegularPoints(2, 2))) +end From 75e8b474331d8db827ffaf7b895ec847a84dff80 Mon Sep 17 00:00:00 2001 From: Niklas Hackelberg Date: Tue, 3 Nov 2026 22:42:34 +0100 Subject: [PATCH 04/11] Geometry: MultiRay docstring before its definition The types for the in-plane rule sat between the docstring and the struct, so the docstring documented FullCell instead of MultiRay. Co-Authored-By: Claude Opus 5.5 --- src/geometry/multiray.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/geometry/multiray.jl b/src/geometry/multiray.jl index 46d48ad..855696b 100644 --- a/src/geometry/multiray.jl +++ b/src/geometry/multiray.jl @@ -150,6 +150,15 @@ weights_eltype(w::GaussianWeights, ::Type{T}) where {T} = GaussianWeights{T}(w.s @inline subray_weight(::UniformWeights, du, dv) = one(du) @inline subray_weight(w::GaussianWeights, du, dv) = exp(-((du / w.sigma_u)^2 + (dv / w.sigma_v)^2) / 2) +# Whether the subrays of a cell move along both detector axes or only along u. For 2D images (one +# voxel layer along z) the rays must stay in the image plane, which the projectors ensure with +# `in_plane`. +struct FullCell end +struct InPlane end + +@inline restrict(::FullCell, du, dv) = (du, dv) +@inline restrict(::InPlane, du, dv) = (du, zero(dv)) + """ MultiRay(geometry, points; weights=UniformWeights()) @@ -179,15 +188,6 @@ g = MultiRay(circular_cone_beam(; angles, dso, dsd, nu, nv, du), JitteredPoints( op = JosephProjector(image, g) ``` """ -# Whether the subrays of a cell move along both detector axes or only along u. For 2D images (one -# voxel layer along z) the rays must stay in the image plane, which the projectors ensure with -# `in_plane`. -struct FullCell end -struct InPlane end - -@inline restrict(::FullCell, du, dv) = (du, dv) -@inline restrict(::InPlane, du, dv) = (du, zero(dv)) - struct MultiRay{T,G<:AbstractDetectorGeometry{T},P,W,V,S} <: AbstractProjectionGeometry{T} geometry::G points::P From 7dff5f2f900d694c6a2c379f25b38e20a8d671f2 Mon Sep 17 00:00:00 2001 From: Niklas Hackelberg Date: Tue, 3 Nov 2026 22:42:34 +0100 Subject: [PATCH 05/11] Docs: guides for several rays per cell How-to "Use Several Rays per Cell" (the sampling problem of one ray per cell on a 2D example, positions on the detector, MultiRay with regular and jittered points) and explanation "Several Rays per Cell" (model, point and weighting strategies, random positions from a hash, 2D images, implementation), linked from the home page and Siddon's page; GPU timings in the benchmark README. Co-Authored-By: Claude Opus 5.5 --- benchmark/README.md | 17 ++++ docs/literate/explanation/multiray.jl | 105 +++++++++++++++++++++ docs/literate/explanation/siddon.jl | 3 +- docs/literate/howto/multiray.jl | 129 ++++++++++++++++++++++++++ docs/make.jl | 2 + docs/src/index.md | 4 +- docs/utils/multiray_figures.jl | 58 ++++++++++++ 7 files changed, 316 insertions(+), 2 deletions(-) create mode 100644 docs/literate/explanation/multiray.jl create mode 100644 docs/literate/howto/multiray.jl create mode 100644 docs/utils/multiray_figures.jl diff --git a/benchmark/README.md b/benchmark/README.md index fa26950..9eab0fe 100644 --- a/benchmark/README.md +++ b/benchmark/README.md @@ -157,3 +157,20 @@ Rays are clipped to the bounding box of the image mask: a small region of intere 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. + +### Several rays per cell (GPU) + +AMD Radeon RX 7900 XT, Float32, 128³ volume, 128² detector, 180 views (circular cone beam), +minimum of 5 runs, 2026-11-03, branch `nh/multiray` (the GPU was also the desktop GPU in use: +approximate). + +| geometry | Siddon forward / back | Joseph forward / back | +|---|---|---| +| plain | 6.0 / 13.7 ms | 6.8 / 31.0 ms | +| `MultiRay`, `RegularPoints(1)` | 5.8 / 13.5 ms | 6.8 / 31.0 ms | +| `MultiRay`, `RegularPoints(2)` | 24.7 / 53.3 ms | 29.9 / 123.7 ms | +| `MultiRay`, `JitteredPoints(2)` | 25.8 / 54.9 ms | 30.2 / 125.4 ms | +| `MultiRay`, `JitteredPoints(2)`, `GaussianWeights(0.4)` | 24.1 / 54.9 ms | 29.8 / 125.7 ms | + +One ray per cell costs the same as the plain geometry; 4 rays per cell take about 4 times as long. +Random positions (a hash per ray) and Gaussian weights are not measurable. diff --git a/docs/literate/explanation/multiray.jl b/docs/literate/explanation/multiray.jl new file mode 100644 index 0000000..c98bddd --- /dev/null +++ b/docs/literate/explanation/multiray.jl @@ -0,0 +1,105 @@ +# # Several Rays per Cell +# +# This page explains what a projector computes on a `MultiRay` geometry, how the rays are placed and +# weighted, where the random positions come from, and how it is implemented. +# [Use Several Rays per Cell](@ref) shows how to use it. + +using CTOperators +include(joinpath(pkgdir(CTOperators), "docs", "utils", "multiray_figures.jl")) #hide + +# ## What is computed +# +# A ray-driven projector such as [Siddon's Method](@ref) or [Joseph's Method](@ref) treats a +# detector cell as a point: one line from the source to the center of the cell. When the cells are +# large compared to the voxels, the lines miss voxels between them, and the back projection covers +# the voxels unevenly (De Man & Basu 2004). A real cell integrates over its area. `MultiRay(geometry, points; weights)` replaces the one line by `K` rays +# per cell, ending at the positions given by `points`, and averages their line integrals: +# +# ```math +# y_i = \sum_{k=1}^K \frac{w_{ik}}{\sum_l w_{il}} \int_{\text{ray } (i, k)} x +# ``` +# +# The back projection is the exact adjoint: every ray adds its share ``w_{ik} / \sum_l w_{il}`` of +# ``y_i``. The rays of a cell approximate the integral over the cell, in the classification of +# system models by Miao (2013) the "area integral" model, by sampling. +# +# The weights average **line integrals**, which keeps the operator linear. A detector averages +# intensities, ``\exp(-\int x)``, which for a cell that sees an edge of the object is not the same +# (the partial-volume effect); modeling that would make the operator nonlinear. +# +# ## Positions on the detector +# +# Positions are counted in cells: cell `(iu, iv)` has its center at `(iu, iv)` and its edges at +# half-integers. Every detector geometry knows the ray that ends at any position +# (`CTOperators.cell_ray`): on a flat detector a point on the plane, on an arc detector a point on +# the cylinder around the source at the angle that belongs to the position, on a parallel-beam +# detector a line through that point. So the rays of a cell sample its actual shape, also on curved +# detectors. +# +# The point strategies give the offsets of the rays from the cell center, between -1/2 and 1/2; +# the weighting strategies give a weight for every offset. Marker areas show the weights: + +fig = points_figure([("RegularPoints(3, 3)", RegularPoints(3, 3), UniformWeights()), #hide + ("JitteredPoints(3, 3)", JitteredPoints(3, 3; seed=1), UniformWeights()), #hide + ("RandomPoints(9)", RandomPoints(9; seed=1), UniformWeights()), #hide + ("RegularPoints(3, 3),\nGaussianWeights(0.3)", RegularPoints(3, 3), GaussianWeights(0.3))]) #hide +save("multiray_points.png", fig) #hide +nothing #hide +# ![Point strategies in one cell](multiray_points.png) + +# * **Regular points** sit at the centers of `su × sv` subcells: the midpoint rule for the integral +# over the cell. For data that are cell averages they are the most accurate with few points. +# * **Jittered points** move each regular point to a random position within its subcell ("jittered +# sampling", Cook 1986, also called stratified sampling). They keep the even spread, but neighboring +# cells no longer share the same pattern, so regular structures of the rays do not line up with +# the voxel grid. +# * **Random points** are spread anywhere in the cell. With few points their mean can differ +# noticeably from the cell average. +# * **Given points** (`Points`) are the same in every cell, also outside it. +# * **Gaussian weights** model a detector response that falls off towards the edges of a cell. +# +# For 2D images, which have one layer of voxels along `z`, the projectors keep all rays in the image +# plane: the offsets along `v` are set to 0, and strategies with several subcells along `v` are +# rejected. +# +# ## Random positions without a random number generator +# +# The random positions must be the same in the forward and in the back projection, in every call +# and on every device: otherwise the back projection would not be the adjoint of the forward +# projection, and an iterative solver would see a different matrix in every iteration. A random +# number generator with state cannot give that in parallel kernels. Instead the positions are +# computed from the indices with a hash function, as GPU renderers do: `pcg4d` of Jarzynski and +# Olano (2020) maps four numbers, here the detector cell, the view, the number of the ray in the cell +# and the seed, to four numbers that look independent even for neighboring inputs; two of them +# give the offsets. +# +# The view is numbered as in the full scan: `subset(op, views)` keeps these numbers, so the rays of +# a subset are those of the same views of the full operator, as ordered-subset methods require. A +# different seed gives a different, again fixed, set of positions: `MultiRay(geometry, +# JitteredPoints(4; seed=2))`. +# +# ## Implementation +# +# * The number of rays per cell `K` is a type parameter of the point strategy. The loop over the +# rays of a cell has a length known to the compiler, and with `K = 1` the projector is as fast as +# on the plain geometry. +# * The kernels are those of the plain projectors: each thread handles one data element. For a +# `MultiRay` it traces the `K` rays one after the other and adds them up (forward), or first sums +# the weights and then lets every ray add its share (back). No memory per ray is needed. +# * The time is about `K` times that of one ray; computing the positions and weights is negligible +# (see the [benchmarks](https://github.com/JuliaImageRecon/CTOperators.jl/blob/main/benchmark/README.md)). +# * The ray mask applies per data element (a masked element skips all its rays), the clipping to the +# bounding box of the image mask per ray. +# * A `MultiRay` has no single ray per data element. A projector that traces only one ray per cell +# does not accept it, instead of silently ignoring the extra rays. +# +# ## References +# +# - R. L. Cook, "Stochastic sampling in computer graphics", ACM Transactions on Graphics 5(1), +# 51–72, 1986. +# - M. Jarzynski, M. Olano, "Hash functions for GPU rendering", Journal of Computer Graphics +# Techniques 9(3), 21–38, 2020. +# - C. Miao, "Comparative studies of different system models for iterative CT image +# reconstruction", master's thesis, Wake Forest University, 2013. +# - B. De Man, S. Basu, "Distance-driven projection and backprojection in three dimensions", Phys. +# Med. Biol. 49(11), 2463–2475, 2004. diff --git a/docs/literate/explanation/siddon.jl b/docs/literate/explanation/siddon.jl index 32003b3..61c5021 100644 --- a/docs/literate/explanation/siddon.jl +++ b/docs/literate/explanation/siddon.jl @@ -74,7 +74,8 @@ nothing #hide # ![Coverage of the pixels by the rays of one view](coverage.png) # In an iterative reconstruction this pattern can appear as a checkerboard in the image (De Man & -# Basu 2004 describe the effect). Remedies are more rays than voxels (finer detector sampling), +# Basu 2004 describe the effect). Remedies are more rays than voxels (finer detector sampling, or +# several rays per cell, see [Use Several Rays per Cell](@ref)), # [Joseph's Method](@ref), whose interpolation evens out the coverage at most angles, or projectors # that integrate over the detector cells, such as distance-driven or separable-footprint projectors. # diff --git a/docs/literate/howto/multiray.jl b/docs/literate/howto/multiray.jl new file mode 100644 index 0000000..7476595 --- /dev/null +++ b/docs/literate/howto/multiray.jl @@ -0,0 +1,129 @@ +# # Use Several Rays per Cell +# +# A detector cell measures the radiation that falls onto its whole area. The projectors trace one +# ray per cell, through its center. When the cells are large compared to the voxels, these rays +# miss voxels between them, and the reconstruction shows patterns that come from the sampling, not +# from the object. This guide shows the problem on a small example and how `MultiRay` traces several +# rays per cell instead. [Several Rays per Cell](@ref) explains how the rays are placed and weighted. + +using CTOperators, ImagePhantoms, RegularizedLeastSquares, LinearOperatorCollection, CairoMakie +using ImageGeoms: axis +include(joinpath(pkgdir(CTOperators), "docs", "utils", "tutorial_figures.jl")) #hide +include(joinpath(pkgdir(CTOperators), "docs", "utils", "multiray_figures.jl")) #hide + +# ## The example +# +# A 2D Shepp-Logan phantom of 128 × 128 pixels and a parallel-beam scan with 180 views. The detector +# cells are 2 pixels wide: + +n = 128 +image = ImageGeom(; dims=(n, n), deltas=(1.0, 1.0)) +x_true = Float32.(phantom(axis(image, 1), axis(image, 2), shepp_logan(SheppLoganToft(); fovs=(n, n)), 2)) +geometry = parallel_beam(; angles=range(0, pi, 181)[1:180], nu=96, nv=1, du=2.0); + +# The data should come from cells that average over their width, as a real detector does. They are +# simulated on a 4 times finer image with 16 rays per cell, so that the data do not come from any of +# the models compared below: + +fine = ImageGeom(; dims=(4n, 4n), deltas=(0.25, 0.25)) +x_fine = Float32.(phantom(axis(fine, 1), axis(fine, 2), shepp_logan(SheppLoganToft(); fovs=(n, n)), 2)) +y = JosephProjector(fine, MultiRay(geometry, RegularPoints(16, 1))) * x_fine; + +# The reconstructions use 40 CGNR iterations: + +function reconstruct(op) + solver = createLinearSolver(CGNR, op; AHA=normalOperator(op), iterations=40) + reshape(solve!(solver, vec(y)), n, n) +end; + +# ## One ray per cell +# +# So far each data element was one ray, through the center of its cell. Positions on the detector +# are counted in cells: cell 2 has its center at position 2 and its edges at 1.5 and 2.5. With cells +# twice as wide as the pixels, the rays of a view pass through every other column of pixels (left): + +fig = rays_figure() #hide +save("multiray_rays.png", fig) #hide +nothing #hide +# ![Rays through the cell centers and four rays per cell](multiray_rays.png) + +# Which pixels the rays miss changes from view to view, but not evenly. The column sums of the +# projector (`colsums`, the summed weights of all rays in every pixel) show it: with even coverage +# they would be smooth, here they vary by about 10 % in a fine pattern. + +op = SiddonProjector(image, geometry) +x_center = reconstruct(op); + +# ## Several rays per cell +# +# `MultiRay` wraps the geometry and traces several rays per cell, here 4 at the centers of 4 +# subcells along `u` (right in the figure above). The projection of a cell is the mean of its rays' +# line integrals; the data keep their size. + +op_multi = SiddonProjector(image, MultiRay(geometry, RegularPoints(4, 1))) +x_multi = reconstruct(op_multi); + +#- + +c, c_multi = colsums(op), colsums(op_multi) #hide +m = sum(c[49:80, 49:80]) / 32^2 #hide +fig = images_figure([("column sums, 1 ray per cell", c), ("column sums, 4 rays per cell", c_multi)]; #hide + colorrange=(0.8m, 1.2m)) #hide +save("multiray_colsums.png", fig) #hide +nothing #hide +# ![Column sums with one and with four rays per cell](multiray_colsums.png) + +# Both images are shown in the window of ±20 % around the mean. With 4 rays per cell the column sums +# vary by less than 1 % inside the field of view. The +# reconstructions, in the color window 0 to 0.4 with the relative error +# ``\|x - x_\text{true}\| / \|x_\text{true}\|``: + +fig = images_figure([("phantom", x_true), ("1 ray per cell", x_center), ("4 rays per cell", x_multi)]; #hide + colorrange=(0, 0.4), truth=x_true) #hide +save("multiray_reconstructions.png", fig) #hide +nothing #hide +# ![Reconstructions with one and with four rays per cell](multiray_reconstructions.png) + +# The model with several rays per cell is closer to what the cells measure, and the error drops from +# about 0.33 to 0.27. `JosephProjector` takes a `MultiRay` the same way. In 2D the rays stay in the +# image plane, so the points are spread along `u` only (`RegularPoints(4, 1)`); for a 3D image, +# `RegularPoints(4)` places 4 × 4 rays per cell. +# +# ## Other placements and weights +# +# Besides the regular points, a cell can be sampled at random positions, fixed by a seed: +# +# * `JitteredPoints(4, 1; seed=1)`: each of the regular points moved to a random position within its +# subcell; +# * `RandomPoints(4; seed=1)`: 4 random points anywhere in the cell; +# * `Points(offsets)`: given positions, the same in every cell; +# +# and the rays can be weighted, e.g. by a detector response that falls off towards the cell edges: + +g = MultiRay(geometry, JitteredPoints(4, 1; seed=1); weights=GaussianWeights(0.4)) +op_jittered = SiddonProjector(image, g) +x_jittered = reconstruct(op_jittered); + +#- + +fig = images_figure([("4 regular rays", x_multi), ("4 jittered rays,\nGaussian weights", x_jittered)]; #hide + colorrange=(0, 0.4), truth=x_true) #hide +save("multiray_jittered.png", fig) #hide +nothing #hide +# ![Regular and jittered rays](multiray_jittered.png) + +# With data that are exact cell averages, as here, the regular points fit best: they are the +# midpoint rule for the average over the cell. The jittered points with Gaussian weights describe a +# different detector, one whose response falls off towards the cell edges, and the error is larger. +# Random positions break up the regular pattern of the rays; which placement and weights fit a +# scanner depends on its detector. +# +# ## Cost and devices +# +# * The time grows with the number of rays per cell: 4 rays take about 4 times as long; the +# placement and the weights cost nothing measurable (see the +# [benchmarks](https://github.com/JuliaImageRecon/CTOperators.jl/blob/main/benchmark/README.md)). +# No extra memory is needed. +# * `MultiRay` works on every device (`storage=CuArray` etc.), and the random positions are the same +# on all of them, in the forward and the back projection, and for `subset(op, views)`. +# * Masks work as before: a masked data element skips all its rays. diff --git a/docs/make.jl b/docs/make.jl index 38fe68b..ebbe9d5 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -42,6 +42,7 @@ makedocs(; "Run on a GPU" => "generated/howto/gpu.md", "Ordered Subsets and Normalization" => "generated/howto/subsets.md", "Use Masks" => "generated/howto/masks.md", + "Use Several Rays per Cell" => "generated/howto/multiray.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", @@ -52,6 +53,7 @@ makedocs(; "Siddon's Method" => "generated/explanation/siddon.md", "Joseph's Method" => "generated/explanation/joseph.md", "Masks" => "generated/explanation/masks.md", + "Several Rays per Cell" => "generated/explanation/multiray.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 a35152c..592c182 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -52,12 +52,14 @@ y = op * rand(Float32, 128, 128, 128) # projections, size (256, 128, 360) 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 + [use several rays per detector cell](@ref "Use Several Rays per Cell"), 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"). - 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 + the projection methods, [Siddon's Method](@ref) and [Joseph's Method](@ref), + [several rays per cell](@ref "Several Rays per Cell"), and [workgroups and ray order on GPUs](@ref "Workgroups and Ray Order on GPUs"). - The reference lists all types and functions by domain: [Images](api/images.md), [Projection Geometries](api/geometries.md), [Operators](api/operators.md), [Projectors](api/projectors.md) and, for implementers, [Internals](api/internals.md). diff --git a/docs/utils/multiray_figures.jl b/docs/utils/multiray_figures.jl new file mode 100644 index 0000000..e20e8c4 --- /dev/null +++ b/docs/utils/multiray_figures.jl @@ -0,0 +1,58 @@ +# Figures of the pages on several rays per detector cell (CairoMakie): rays of one parallel-beam +# view through a pixel grid, and the positions of the point strategies in one cell. Included +# invisibly by the pages. + +using CairoMakie +using CTOperators: point_offset, subray_weight + +# One parallel-beam view along y through an n × n pixel grid, with detector cells of width du +# (pixels) below it: the rays through the cell centers (left) and `k` regular rays per cell (right). +# Positions on the detector are in cell units; integers are the cell centers, half-integers the +# edges. +function rays_figure(; n=8, du=2.0, k=4, size=(820, 470)) + ncells = round(Int, n / du) + 2 + fig = Figure(; size) + for (col, (title, offsets)) in enumerate((("one ray per cell, at its center", [0.0]), + ("$k rays per cell", [(j + 0.5) / k - 0.5 for j in 0:k-1]))) + ax = Axis(fig[1, col]; aspect=DataAspect(), title, limits=(-0.5, n + 0.5, -2.6, n + 0.5)) + hidedecorations!(ax); hidespines!(ax) + for w in 0:n + lines!(ax, [w, w], [0, n]; color=:gray75, linewidth=0.8) + lines!(ax, [0, n], [w, w]; color=:gray75, linewidth=0.8) + end + # the detector: cells c = 1, ..., ncells with centers at x = (c - (ncells + 1)/2) du + n/2 + x(c) = (c - (ncells + 1) / 2) * du + n / 2 + for c in 1:ncells + poly!(ax, Rect2f(x(c - 0.5), -1.2, du, 0.5); color=isodd(c) ? (:steelblue, 0.25) : (:steelblue, 0.45), + strokecolor=:steelblue, strokewidth=0.8) + for o in offsets + xr = x(c + o) + 0 <= xr <= n && lines!(ax, [xr, xr], [-0.7, n]; color=:darkorange, linewidth=1.5) + end + end + if col == 1 + for c in 2:3 + text!(ax, x(c), -1.35; text="$c", align=(:center, :top), fontsize=13, color=:steelblue) + text!(ax, x(c + 0.5), -1.35; text="$(c + 0.5)", align=(:center, :top), fontsize=11, color=:gray40) + end + text!(ax, n / 2, -2.0; text="positions on the detector, in cells", align=(:center, :top), fontsize=12, color=:gray30) + end + end + fig +end + +# The points of several strategies in one cell (-1/2 to 1/2 in u and v), the marker area +# proportional to the weight. `strategies` holds (title, points, weights) tuples. +function points_figure(strategies; size=(250 * length(strategies), 290)) + fig = Figure(; size) + for (col, (title, points, weights)) in enumerate(strategies) + ax = Axis(fig[1, col]; aspect=DataAspect(), title, limits=(-0.6, 0.6, -0.6, 0.6)) + hidedecorations!(ax) + poly!(ax, Rect2f(-0.5, -0.5, 1, 1); color=(:steelblue, 0.12), strokecolor=:steelblue, strokewidth=1.5) + K = typeof(points).parameters[1] + offs = [point_offset(points, Float64, 3, 2, 1, k) for k in 1:K] + w = [subray_weight(weights, o...) for o in offs] + scatter!(ax, first.(offs), last.(offs); markersize=14 .* sqrt.(w ./ maximum(w)), color=:darkorange) + end + fig +end From 4f80eee3cac6be37d7a61926af1d58d13ba2750f Mon Sep 17 00:00:00 2001 From: Niklas Hackelberg Date: Wed, 4 Nov 2026 18:37:58 +0100 Subject: [PATCH 06/11] Geometry: GaussianPoints n random points per cell, normally distributed around the cell center, all rays counting equally. The cell of a data element is computed once per element instead of once per subray. Co-Authored-By: Claude Opus 5.5 --- docs/literate/explanation/multiray.jl | 9 +++- docs/literate/howto/multiray.jl | 2 + docs/utils/multiray_figures.jl | 4 +- src/CTOperators.jl | 3 +- src/geometry/multiray.jl | 74 +++++++++++++++++++++++---- src/raytracing/subrays.jl | 13 ++--- test/geometry/multiray.jl | 26 ++++++++++ test/operators/multiray.jl | 2 +- 8 files changed, 111 insertions(+), 22 deletions(-) diff --git a/docs/literate/explanation/multiray.jl b/docs/literate/explanation/multiray.jl index c98bddd..cc492e4 100644 --- a/docs/literate/explanation/multiray.jl +++ b/docs/literate/explanation/multiray.jl @@ -41,7 +41,8 @@ include(joinpath(pkgdir(CTOperators), "docs", "utils", "multiray_figures.jl")) # fig = points_figure([("RegularPoints(3, 3)", RegularPoints(3, 3), UniformWeights()), #hide ("JitteredPoints(3, 3)", JitteredPoints(3, 3; seed=1), UniformWeights()), #hide - ("RandomPoints(9)", RandomPoints(9; seed=1), UniformWeights()), #hide + ("RandomPoints(9)", RandomPoints(9; seed=3), UniformWeights()), #hide + ("GaussianPoints(9, 0.3)", GaussianPoints(9, 0.3; seed=1), UniformWeights()), #hide ("RegularPoints(3, 3),\nGaussianWeights(0.3)", RegularPoints(3, 3), GaussianWeights(0.3))]) #hide save("multiray_points.png", fig) #hide nothing #hide @@ -55,8 +56,12 @@ nothing #hide # the voxel grid. # * **Random points** are spread anywhere in the cell. With few points their mean can differ # noticeably from the cell average. +# * **Gaussian points** are normally distributed around the cell center, and all rays count equally: +# the distribution of the positions is the weighting. Points may land outside the cell, which +# models blur across neighboring cells. # * **Given points** (`Points`) are the same in every cell, also outside it. -# * **Gaussian weights** model a detector response that falls off towards the edges of a cell. +# * **Gaussian weights** describe a similar response with fixed positions: the points stay where +# they are, and the rays near the edges count less. # # For 2D images, which have one layer of voxels along `z`, the projectors keep all rays in the image # plane: the offsets along `v` are set to 0, and strategies with several subcells along `v` are diff --git a/docs/literate/howto/multiray.jl b/docs/literate/howto/multiray.jl index 7476595..8848859 100644 --- a/docs/literate/howto/multiray.jl +++ b/docs/literate/howto/multiray.jl @@ -96,6 +96,8 @@ nothing #hide # * `JitteredPoints(4, 1; seed=1)`: each of the regular points moved to a random position within its # subcell; # * `RandomPoints(4; seed=1)`: 4 random points anywhere in the cell; +# * `GaussianPoints(4, 0.3; seed=1)`: 4 random points, normally distributed around the cell center, +# which can also land in neighboring cells (blur); # * `Points(offsets)`: given positions, the same in every cell; # # and the rays can be weighted, e.g. by a detector response that falls off towards the cell edges: diff --git a/docs/utils/multiray_figures.jl b/docs/utils/multiray_figures.jl index e20e8c4..df7f6db 100644 --- a/docs/utils/multiray_figures.jl +++ b/docs/utils/multiray_figures.jl @@ -43,10 +43,10 @@ end # The points of several strategies in one cell (-1/2 to 1/2 in u and v), the marker area # proportional to the weight. `strategies` holds (title, points, weights) tuples. -function points_figure(strategies; size=(250 * length(strategies), 290)) +function points_figure(strategies; size=(220 * length(strategies), 270)) fig = Figure(; size) for (col, (title, points, weights)) in enumerate(strategies) - ax = Axis(fig[1, col]; aspect=DataAspect(), title, limits=(-0.6, 0.6, -0.6, 0.6)) + ax = Axis(fig[1, col]; aspect=DataAspect(), title, limits=(-0.9, 0.9, -0.9, 0.9)) hidedecorations!(ax) poly!(ax, Rect2f(-0.5, -0.5, 1, 1); color=(:steelblue, 0.12), strokecolor=:steelblue, strokewidth=1.5) K = typeof(points).parameters[1] diff --git a/src/CTOperators.jl b/src/CTOperators.jl index 15759d8..f4dc0af 100644 --- a/src/CTOperators.jl +++ b/src/CTOperators.jl @@ -12,7 +12,8 @@ export AbstractCTOperator export ImageGeom, MaskCircle, MaskEllipse, MaskAllButEdge export AbstractProjectionGeometry, AbstractDetectorGeometry, RaySegment, RayLine export RayList, ConeBeam, ArcConeBeam, ParallelBeam -export MultiRay, RegularPoints, JitteredPoints, RandomPoints, Points, UniformWeights, GaussianWeights +export MultiRay, RegularPoints, JitteredPoints, RandomPoints, GaussianPoints, Points +export UniformWeights, GaussianWeights export circular_cone_beam, helical_cone_beam, parallel_beam export datasize, nrays, ray, raytype, nviews, detectorsize, subset, projection_geometry export AbstractProjector, AdjointCTOperator, SiddonProjector, JosephProjector, LaunchConfig diff --git a/src/geometry/multiray.jl b/src/geometry/multiray.jl index 855696b..3fe0aec 100644 --- a/src/geometry/multiray.jl +++ b/src/geometry/multiray.jl @@ -39,6 +39,21 @@ struct RandomPoints{K} seed::UInt32 end +""" + GaussianPoints(n, sigma_u, sigma_v=sigma_u; seed=0) + +`n` random points per detector cell, normally distributed around the cell center with standard +deviations `sigma_u`, `sigma_v` in cell units: a detector response set by where the rays land, with +all rays counting equally. Points can land outside the cell, which models blur across neighboring +cells. The positions are fixed like those of [`JitteredPoints`](@ref). Point strategy of +[`MultiRay`](@ref). +""" +struct GaussianPoints{K,T} + sigma_u::T + sigma_v::T + seed::UInt32 +end + """ Points(offsets) @@ -69,6 +84,13 @@ function RandomPoints(n::Integer; seed::Integer=0) RandomPoints{Int(n)}(UInt32(seed)) end +function GaussianPoints(n::Integer, sigma_u::Real, sigma_v::Real=sigma_u; seed::Integer=0) + n > 0 || throw(ArgumentError("number of points must be positive, got $n")) + sigma_u > 0 && sigma_v > 0 || throw(ArgumentError("standard deviations must be positive, got ($sigma_u, $sigma_v)")) + T = float(promote_type(typeof(sigma_u), typeof(sigma_v))) + GaussianPoints{Int(n),T}(sigma_u, sigma_v, UInt32(seed)) +end + function Points(offsets::AbstractMatrix{<:Real}) size(offsets, 1) == 2 && size(offsets, 2) > 0 || throw(DimensionMismatch("offsets must be 2 × K with K > 0, got $(size(offsets))")) @@ -77,10 +99,11 @@ function Points(offsets::AbstractMatrix{<:Real}) Points{K,T}(SVector{K}(ntuple(k -> SVector{2,T}(offsets[1, k], offsets[2, k]), K))) end -npoints(::Union{RegularPoints{K},JitteredPoints{K},RandomPoints{K},Points{K}}) where {K} = K +npoints(::Union{RegularPoints{K},JitteredPoints{K},RandomPoints{K},GaussianPoints{K},Points{K}}) where {K} = K -# the points in element type T (only given points store numbers) +# the points in element type T (given and Gaussian points store numbers) points_eltype(p, ::Type{T}) where {T} = p +points_eltype(p::GaussianPoints{K}, ::Type{T}) where {K,T} = GaussianPoints{K,T}(p.sigma_u, p.sigma_v, p.seed) points_eltype(p::Points{K}, ::Type{T}) where {K,T} = Points{K,T}(map(o -> SVector{2,T}(o), p.offsets)) # The offset (du, dv) of subray k in the cell (iu, iv) of view `view` (numbered as in the full @@ -104,6 +127,14 @@ end ru - half, rv - half end +# normally distributed from two uniform numbers (Box-Muller); 1 - ru lies in (0, 1] +@inline function point_offset(p::GaussianPoints, ::Type{T}, iu, iv, view, k) where {T} + ru, rv = random_pair(T, p.seed, iu, iv, view, k) + r = sqrt(-2 * log(1 - ru)) + s, c = sincos(2 * T(pi) * rv) + T(p.sigma_u) * r * c, T(p.sigma_v) * r * s +end + @inline function point_offset(p::Points, ::Type{T}, iu, iv, view, k) where {T} o = @inbounds p.offsets[k] T(o[1]), T(o[2]) @@ -164,7 +195,8 @@ struct InPlane end The detector geometry `geometry` (`ConeBeam`, `ArcConeBeam`, `ParallelBeam`) with several rays per detector cell. The rays of a cell end at the positions given by the point strategy `points` -([`RegularPoints`](@ref), [`JitteredPoints`](@ref), [`RandomPoints`](@ref), [`Points`](@ref)) and +([`RegularPoints`](@ref), [`JitteredPoints`](@ref), [`RandomPoints`](@ref), +[`GaussianPoints`](@ref), [`Points`](@ref)) and are averaged with the weights of `weights` ([`UniformWeights`](@ref), [`GaussianWeights`](@ref)), normalized to sum 1 per cell. A ray-driven projector built on it computes @@ -248,16 +280,38 @@ end @inline full_view(::Nothing, view) = view @inline full_view(views, view) = @inbounds unsafe_trunc(Int, views[view]) +# the cell of data element i, once per element: its detector indices, the view and the view number +# in the full geometry +@inline function element_cell(g::MultiRay, i) + nu, nv = detectorsize(g.geometry) + iu, iv, view = pixel_index(i, nu, nv) + iu, iv, view, full_view(g.views, view) +end + +# the offset of subray k in the cell +@inline function subray_offset(g::MultiRay{T}, cell, k) where {T} + iu, iv, _, number = cell + restrict(g.span, point_offset(g.points, T, iu, iv, number, k)...) +end + +# the weight (not normalized) of subray k in the cell +@inline function subray_weight(g::MultiRay, cell, k) + du, dv = subray_offset(g, cell, k) + subray_weight(g.weights, du, dv) +end + +# subray k in the cell +@inline function cell_subray(g::MultiRay{T}, cell, k) where {T} + iu, iv, view, _ = cell + du, dv = subray_offset(g, cell, k) + a, b = cell_ray(g.geometry, T(iu) + du, T(iv) + dv, view) + a, b, subray_weight(g.weights, du, dv) +end + """ subray(g::MultiRay, i, k) -> (a, b, w) Subray `k` of data element `i` as two points like [`ray`](@ref), and its weight (not normalized). Callable inside kernels. """ -@inline function subray(g::MultiRay{T}, i, k) where {T} - nu, nv = detectorsize(g.geometry) - iu, iv, view = pixel_index(i, nu, nv) - du, dv = restrict(g.span, point_offset(g.points, T, iu, iv, full_view(g.views, view), k)...) - a, b = cell_ray(g.geometry, T(iu) + du, T(iv) + dv, view) - a, b, subray_weight(g.weights, du, dv) -end +@inline subray(g::MultiRay, i, k) = cell_subray(g, element_cell(g, i), k) diff --git a/src/raytracing/subrays.jl b/src/raytracing/subrays.jl index f863b3d..79d8228 100644 --- a/src/raytracing/subrays.jl +++ b/src/raytracing/subrays.jl @@ -7,25 +7,26 @@ # forward: the weighted mean of the subrays' line integrals @inline function trace_element(segment::F, action::Gather, grid::VoxelGrid{T}, g::MultiRay, i) where {F,T} + cell = element_cell(g, i) acc, wsum = zero(T), zero(T) for k in 1:nsubrays(g) - a, b, w = subray(g, i, k) + a, b, w = cell_subray(g, cell, k) acc = muladd(w, segment(action, grid, a, b, raytype(g)), acc) wsum += w end acc / wsum end -# back: every subray adds its share w / sum(w) of the value; the weights are summed first (cheap -# next to tracing), so that the shares are known before any subray is traced +# back: every subray adds its share w / sum(w) of the value; the weights are summed first, so that +# the shares are known before any subray is traced @inline function trace_element(segment::F, action::Scatter, grid::VoxelGrid{T}, g::MultiRay, i) where {F,T} + cell = element_cell(g, i) wsum = zero(T) for k in 1:nsubrays(g) - _, _, w = subray(g, i, k) - wsum += w + wsum += subray_weight(g, cell, k) end for k in 1:nsubrays(g) - a, b, w = subray(g, i, k) + a, b, w = cell_subray(g, cell, k) segment(Scatter(action.image, action.value * (w / wsum)), grid, a, b, raytype(g)) end zero(T) diff --git a/test/geometry/multiray.jl b/test/geometry/multiray.jl index 7b0397f..bfc1eed 100644 --- a/test/geometry/multiray.jl +++ b/test/geometry/multiray.jl @@ -31,6 +31,32 @@ @test subray(m, 1, 1)[3] == 1 && isapprox(subray(m, 1, 2)[3], exp(-0.125)) end +# Gaussian points: normally distributed around the cell center with the given widths, identical +# rays counting equally; the seed fixes them +@testitem "MultiRay: Gaussian points" tags = [:geometry, :multiray] begin + using Adapt + using CTOperators: point_offset + mean(x) = sum(x) / length(x) + std(x) = sqrt(mean((x .- mean(x)) .^ 2)) + sigma_u, sigma_v = 0.3, 0.6 + p = GaussianPoints(16, sigma_u, sigma_v; seed=4) + offsets = [point_offset(p, Float64, iu, iv, 1, k) for iu in 1:200, iv in 1:200, k in 1:16] + du, dv = vec(first.(offsets)), vec(last.(offsets)) + N = length(du) # 640000 samples: tolerances of 5 standard errors + @test abs(mean(du)) < 5 * sigma_u / sqrt(N) && abs(mean(dv)) < 5 * sigma_v / sqrt(N) + @test isapprox(std(du), sigma_u; rtol=5 / sqrt(2N)) && isapprox(std(dv), sigma_v; rtol=5 / sqrt(2N)) + # fraction of a normal distribution within one standard deviation + @test isapprox(count(abs.(du) .< sigma_u) / N, 0.6827; atol=0.002) + @test abs(mean(du .* dv)) < 5 * sigma_u * sigma_v / sqrt(N) # the two axes are independent + @test all(isfinite, du) && all(isfinite, dv) + @test point_offset(p, Float64, 3, 2, 1, 5) == point_offset(GaussianPoints(16, sigma_u, sigma_v; seed=4), Float64, 3, 2, 1, 5) + @test point_offset(p, Float64, 3, 2, 1, 5) != point_offset(GaussianPoints(16, sigma_u, sigma_v; seed=5), Float64, 3, 2, 1, 5) + g = circular_cone_beam(; angles=range(0, 2pi, 5)[1:4], dso=50, dsd=100, nu=5, nv=3, du=1.5) + @test adapt(Array{Float32}, MultiRay(g, GaussianPoints(4, 0.3))).points isa GaussianPoints{4,Float32} + @test_throws ArgumentError GaussianPoints(4, -0.1) + @test_throws ArgumentError GaussianPoints(0, 0.1) +end + @testitem "MultiRay: subset, adapt, errors" tags = [:geometry, :multiray] begin using Adapt using CTOperators: subray diff --git a/test/operators/multiray.jl b/test/operators/multiray.jl index c630d66..bd9833b 100644 --- a/test/operators/multiray.jl +++ b/test/operators/multiray.jl @@ -51,7 +51,7 @@ end @testitem "MultiRay: adjoint and devices" tags = [:operators, :multiray, :gpu] setup = [Backends, MultiRayProblem] begin using Adapt for P in (SiddonProjector, JosephProjector), - (g, pts) in ((cone(), JitteredPoints(2, 3; seed=1)), (arc(), RegularPoints(2)), (parallel(), RandomPoints(5; seed=2))), + (g, pts) in ((cone(), JitteredPoints(2, 3; seed=1)), (arc(), RegularPoints(2)), (parallel(), RandomPoints(5; seed=2)), (cone(), GaussianPoints(4, 0.3; seed=6))), arrayType in arrayTypes, T in floatTypes @testset "$P, $(nameof(typeof(pts))), $arrayType{$T}" begin m = MultiRay(g, pts; weights=GaussianWeights(0.4)) From 710410bc160bb8e44b2cafcd87aab529eba03a07 Mon Sep 17 00:00:00 2001 From: Niklas Hackelberg Date: Wed, 4 Nov 2026 20:29:21 +0100 Subject: [PATCH 07/11] Geometry: open point and weighting strategies AbstractPoints and AbstractWeights replace the closed lists of strategies, so that types of other packages fit in with one method each. The functions such types and custom projectors use are declared public. Co-Authored-By: Claude Opus 5.5 --- ...plement_operator.jl => write_projector.jl} | 0 src/CTOperators.jl | 10 +- src/geometry/multiray.jl | 97 +++++++++++++++---- src/raytracing/actions.jl | 18 +++- src/raytracing/subrays.jl | 13 ++- test/operators/multiray.jl | 24 +++++ 6 files changed, 133 insertions(+), 29 deletions(-) rename docs/literate/howto/{implement_operator.jl => write_projector.jl} (100%) diff --git a/docs/literate/howto/implement_operator.jl b/docs/literate/howto/write_projector.jl similarity index 100% rename from docs/literate/howto/implement_operator.jl rename to docs/literate/howto/write_projector.jl diff --git a/src/CTOperators.jl b/src/CTOperators.jl index f4dc0af..3a6125d 100644 --- a/src/CTOperators.jl +++ b/src/CTOperators.jl @@ -12,13 +12,19 @@ export AbstractCTOperator export ImageGeom, MaskCircle, MaskEllipse, MaskAllButEdge export AbstractProjectionGeometry, AbstractDetectorGeometry, RaySegment, RayLine export RayList, ConeBeam, ArcConeBeam, ParallelBeam -export MultiRay, RegularPoints, JitteredPoints, RandomPoints, GaussianPoints, Points -export UniformWeights, GaussianWeights +export MultiRay, AbstractPoints, RegularPoints, JitteredPoints, RandomPoints, GaussianPoints, Points +export AbstractWeights, UniformWeights, GaussianWeights export circular_cone_beam, helical_cone_beam, parallel_beam export datasize, nrays, ray, raytype, nviews, detectorsize, subset, projection_geometry export AbstractProjector, AdjointCTOperator, SiddonProjector, JosephProjector, LaunchConfig export geometry, image_geometry, imagesize, rowsums, colsums +# not exported, but stable: the interface for types of other packages +public forward!, backward!, similar_projector, storage_array, device_geometry +public point_offset, random_pair, points_eltype, subray_weight, weights_eltype +public trace_element, in_plane, visit, Gather, Scatter, to_index, clip, VoxelGrid +public image_mask, ray_mask, handles_raymask, projector_masks, ProjectorMasks + """ AbstractCTOperator{T} diff --git a/src/geometry/multiray.jl b/src/geometry/multiray.jl index 3fe0aec..381903a 100644 --- a/src/geometry/multiray.jl +++ b/src/geometry/multiray.jl @@ -3,13 +3,32 @@ # center, +-1/2 the cell edges. The number of subrays K is a type parameter of the point strategy, # so the loop over the subrays has a length known to the compiler. +""" + AbstractPoints{K} + +Supertype of the point strategies of [`MultiRay`](@ref): where the `K` rays of a detector cell end. +A subtype implements [`CTOperators.point_offset`](@ref). If it stores floating-point numbers, it +also implements [`CTOperators.points_eltype`](@ref), so that it follows the element type of the +projector. The type can live in any package; see [Write a Sampling Strategy](@ref). +""" +abstract type AbstractPoints{K} end + +""" + AbstractWeights + +Supertype of the weighting strategies of [`MultiRay`](@ref). A subtype implements +[`CTOperators.subray_weight`](@ref), and [`CTOperators.weights_eltype`](@ref) if it stores +floating-point numbers. The type can live in any package; see [Write a Sampling Strategy](@ref). +""" +abstract type AbstractWeights end + """ RegularPoints(su, sv=su) The centers of `su × sv` subcells of each detector cell: `su` along `u`, `sv` along `v`. Point strategy of [`MultiRay`](@ref). """ -struct RegularPoints{K} +struct RegularPoints{K} <: AbstractPoints{K} su::Int sv::Int end @@ -23,7 +42,7 @@ Graphics 5(1), 51-72, 1986). The positions are fixed by `seed`, the detector cel subray, so they are the same in the forward and the back projection, on every device and in `subset`. Point strategy of [`MultiRay`](@ref). """ -struct JitteredPoints{K} +struct JitteredPoints{K} <: AbstractPoints{K} su::Int sv::Int seed::UInt32 @@ -35,7 +54,7 @@ end `n` random points in each detector cell, fixed like those of [`JitteredPoints`](@ref). Point strategy of [`MultiRay`](@ref). """ -struct RandomPoints{K} +struct RandomPoints{K} <: AbstractPoints{K} seed::UInt32 end @@ -48,7 +67,7 @@ all rays counting equally. Points can land outside the cell, which models blur a cells. The positions are fixed like those of [`JitteredPoints`](@ref). Point strategy of [`MultiRay`](@ref). """ -struct GaussianPoints{K,T} +struct GaussianPoints{K,T} <: AbstractPoints{K} sigma_u::T sigma_v::T seed::UInt32 @@ -61,7 +80,7 @@ Given points, the same in every detector cell: column `k` of the `2 × K` matrix position of subray `k` in cell units, `(0, 0)` the cell center and `±1/2` the cell edges. Points outside the cell are allowed. Point strategy of [`MultiRay`](@ref). """ -struct Points{K,T} +struct Points{K,T} <: AbstractPoints{K} offsets::SVector{K,SVector{2,T}} end @@ -99,15 +118,27 @@ function Points(offsets::AbstractMatrix{<:Real}) Points{K,T}(SVector{K}(ntuple(k -> SVector{2,T}(offsets[1, k], offsets[2, k]), K))) end -npoints(::Union{RegularPoints{K},JitteredPoints{K},RandomPoints{K},GaussianPoints{K},Points{K}}) where {K} = K +npoints(::AbstractPoints{K}) where {K} = K + +""" + points_eltype(points, T) -# the points in element type T (given and Gaussian points store numbers) -points_eltype(p, ::Type{T}) where {T} = p +`points` with the floating-point numbers it stores in element type `T`. The default returns `points` +unchanged, for strategies that store none. +""" +points_eltype(p::AbstractPoints, ::Type{T}) where {T} = p points_eltype(p::GaussianPoints{K}, ::Type{T}) where {K,T} = GaussianPoints{K,T}(p.sigma_u, p.sigma_v, p.seed) points_eltype(p::Points{K}, ::Type{T}) where {K,T} = Points{K,T}(map(o -> SVector{2,T}(o), p.offsets)) -# The offset (du, dv) of subray k in the cell (iu, iv) of view `view` (numbered as in the full -# geometry, so that subsets draw the same random positions). +""" + point_offset(points, T, iu, iv, view, k) -> (du, dv) + +The offset of ray `k` from the center of the detector cell `(iu, iv)` in view `view`, in element type +`T` and in cell units (`±1/2` are the cell edges). `view` is numbered as in the full geometry, so +that subsets draw the same random positions. Called in kernels: no allocations. +""" +function point_offset end + @inline function point_offset(p::RegularPoints, ::Type{T}, iu, iv, view, k) where {T} ju, jv = rem(k - 1, p.su), div(k - 1, p.su) half = T(1 // 2) @@ -140,9 +171,14 @@ end T(o[1]), T(o[2]) end -# two numbers in [0, 1) from the hash of (cell, view, subray, seed); the cell as one number -# (iu + 65536 iv) keeps the four inputs distinct for detectors of up to 65536 cells per row +""" + random_pair(T, seed, iu, iv, view, k) -> (ru, rv) + +Two uniform random numbers in `[0, 1)`, fixed by the arguments: the same on every device and in +every call, and independent for different cells, views, rays and seeds. For random point strategies. +""" @inline function random_pair(::Type{T}, seed, iu, iv, view, k) where {T} + # the cell as one number keeps the four inputs distinct for up to 65536 cells per row cell = (iu % UInt32) + ((iv % UInt32) << 16) hu, hv, _, _ = pcg4d(cell, view % UInt32, k % UInt32, seed) unit_interval(T, hu), unit_interval(T, hv) @@ -153,7 +189,7 @@ end All subrays of a detector cell count equally. Weighting strategy of [`MultiRay`](@ref). """ -struct UniformWeights end +struct UniformWeights <: AbstractWeights end """ GaussianWeights(sigma_u, sigma_v=sigma_u) @@ -162,7 +198,7 @@ Subrays weighted by a Gaussian of the offset from the cell center, with standard `sigma_u`, `sigma_v` in cell units: a detector response that falls off towards the cell edges. Weighting strategy of [`MultiRay`](@ref). """ -struct GaussianWeights{T} +struct GaussianWeights{T} <: AbstractWeights sigma_u::T sigma_v::T function GaussianWeights{T}(sigma_u, sigma_v) where {T} @@ -174,10 +210,23 @@ end GaussianWeights(sigma_u::Real, sigma_v::Real=sigma_u) = GaussianWeights{float(promote_type(typeof(sigma_u), typeof(sigma_v)))}(sigma_u, sigma_v) -weights_eltype(w, ::Type{T}) where {T} = w +""" + weights_eltype(weights, T) + +`weights` with the floating-point numbers it stores in element type `T`. The default returns +`weights` unchanged, for strategies that store none. +""" +weights_eltype(w::AbstractWeights, ::Type{T}) where {T} = w weights_eltype(w::GaussianWeights, ::Type{T}) where {T} = GaussianWeights{T}(w.sigma_u, w.sigma_v) -# the (unnormalized) weight of a subray at offset (du, dv) +""" + subray_weight(weights, du, dv) + +The weight of a ray at offset `(du, dv)` from the cell center, positive and not normalized (the +weights of a cell are normalized in the kernel). Called in kernels: no allocations. +""" +function subray_weight end + @inline subray_weight(::UniformWeights, du, dv) = one(du) @inline subray_weight(w::GaussianWeights, du, dv) = exp(-((du / w.sigma_u)^2 + (dv / w.sigma_v)^2) / 2) @@ -196,9 +245,9 @@ struct InPlane end The detector geometry `geometry` (`ConeBeam`, `ArcConeBeam`, `ParallelBeam`) with several rays per detector cell. The rays of a cell end at the positions given by the point strategy `points` ([`RegularPoints`](@ref), [`JitteredPoints`](@ref), [`RandomPoints`](@ref), -[`GaussianPoints`](@ref), [`Points`](@ref)) and +[`GaussianPoints`](@ref), [`Points`](@ref), or an own [`AbstractPoints`](@ref)) and are averaged with the weights of `weights` ([`UniformWeights`](@ref), -[`GaussianWeights`](@ref)), normalized to sum 1 per cell. A ray-driven projector built on it +[`GaussianWeights`](@ref), or an own [`AbstractWeights`](@ref)), normalized to sum 1 per cell. A ray-driven projector built on it computes y[i] = sum_k w_k * (line integral along subray k of data element i) @@ -230,7 +279,7 @@ struct MultiRay{T,G<:AbstractDetectorGeometry{T},P,W,V,S} <: AbstractProjectionG span::S # FullCell() or InPlane() end -function MultiRay(geometry, points; weights=UniformWeights()) +function MultiRay(geometry, points::AbstractPoints; weights::AbstractWeights=UniformWeights()) g = projection_geometry(geometry) g isa AbstractDetectorGeometry || throw(ArgumentError("MultiRay needs a detector geometry (cells to sample), got $(nameof(typeof(g)))")) @@ -258,7 +307,13 @@ end Adapt.adapt_structure(to, g::MultiRay) = multiray(adapt(to, g.geometry), g.points, g.weights, adapt(to, g.views), g.span) -# the geometry with its subrays in the plane v = 0, for 2D images (plain geometries are unchanged) +""" + in_plane(geometry) + +`geometry` with the rays of a [`MultiRay`](@ref) in the image plane (offsets along `u` only); other +geometries unchanged. The constructor of a projector calls it for 2D images, which have one layer +of voxels along `z`. +""" in_plane(g) = g function in_plane(g::MultiRay) check_in_plane(g.points) @@ -295,7 +350,7 @@ end end # the weight (not normalized) of subray k in the cell -@inline function subray_weight(g::MultiRay, cell, k) +@inline function cell_weight(g::MultiRay, cell, k) du, dv = subray_offset(g, cell, k) subray_weight(g.weights, du, dv) end diff --git a/src/raytracing/actions.jl b/src/raytracing/actions.jl index 50a9731..9af5d9c 100644 --- a/src/raytracing/actions.jl +++ b/src/raytracing/actions.jl @@ -1,14 +1,28 @@ -# What a ray does with the voxels it visits, shared by the ray-driven methods: sum the weighted -# image values (forward projection) or add its weighted value to them (back projection). +""" + Gather(image) + +What a ray does in the forward projection: it sums the weighted values of the voxels it visits. +""" struct Gather{A} image::A end +""" + Scatter(image, value) + +What a ray does in the back projection: it adds its weighted `value` to the voxels it visits. +""" struct Scatter{A,T} image::A value::T end +""" + visit(action, acc, lin, w) -> acc + +A ray meets the voxel with linear index `lin`, with weight `w`: [`CTOperators.Gather`](@ref) adds to +the accumulated value `acc`, [`CTOperators.Scatter`](@ref) writes to the image. Called in kernels. +""" @inline visit(g::Gather, acc, lin, w) = @inbounds muladd(w, g.image[lin], acc) @inline function visit(s::Scatter, acc, lin, w) diff --git a/src/raytracing/subrays.jl b/src/raytracing/subrays.jl index 79d8228..632c24a 100644 --- a/src/raytracing/subrays.jl +++ b/src/raytracing/subrays.jl @@ -1,7 +1,12 @@ -# The rays of data element i for the ray-driven methods: one ray for plain geometries, the weighted -# subrays of a MultiRay. `segment(action, grid, a, b, raytype)` traces one ray from a to b and -# returns the value of the action (Gather: the weighted sum of the image). +""" + trace_element(segment, action, grid, geometry, i) +Traces the rays of data element `i` for a ray-driven method: one ray for plain geometries, the +weighted rays of a cell for a [`MultiRay`](@ref). `segment(action, grid, a, b, raytype)` traces one +ray from `a` to `b`, calling [`CTOperators.visit`](@ref) for every voxel, and returns the value of +the action. The result is the value of data element `i` for a [`CTOperators.Gather`](@ref). A +projector that traces through this function works on every geometry. Called in kernels. +""" @inline trace_element(segment::F, action, grid, geometry, i) where {F} = segment(action, grid, ray(geometry, i)..., raytype(geometry)) @@ -23,7 +28,7 @@ end cell = element_cell(g, i) wsum = zero(T) for k in 1:nsubrays(g) - wsum += subray_weight(g, cell, k) + wsum += cell_weight(g, cell, k) end for k in 1:nsubrays(g) a, b, w = cell_subray(g, cell, k) diff --git a/test/operators/multiray.jl b/test/operators/multiray.jl index bd9833b..a10473b 100644 --- a/test/operators/multiray.jl +++ b/test/operators/multiray.jl @@ -117,3 +117,27 @@ end # subcells along v would collapse onto the plane @test_throws ArgumentError SiddonProjector(image2, MultiRay(fan(:flat), RegularPoints(2, 2))) end + +# Strategies are open: types of other packages only subtype AbstractPoints and AbstractWeights and +# add a method each +@testitem "MultiRay: own strategies" tags = [:operators, :multiray] setup = [MultiRayProblem] begin + struct Cross{K} <: AbstractPoints{K} + reach::Float64 + end + CTOperators.point_offset(p::Cross, ::Type{T}, iu, iv, view, k) where {T} = + isodd(k) ? (T(p.reach) * (k - 2), zero(T)) : (zero(T), T(p.reach) * (k - 3)) + struct Linear <: AbstractWeights end + CTOperators.subray_weight(::Linear, du, dv) = 1 - abs(du) - abs(dv) + + g = cone() + m = MultiRay(g, Cross{4}(0.25); weights=Linear()) + equal = MultiRay(g, Points([-0.25 0.0 0.25 0.0; 0.0 -0.25 0.0 0.25])) + x, y = rand(image.dims...), rand(datasize(g)...) + for P in (SiddonProjector, JosephProjector) + op = P(image, m; eltype=Float64) + # the weights are equal at equal distance from the center: the same as given points + @test isapprox(op * x, P(image, equal; eltype=Float64) * x) + @test isapprox(sum((op * x) .* y), sum(x .* (adjoint(op) * y))) + end + @test_throws MethodError MultiRay(g, [0.0; 0.0;;]) +end From 1914205e5be93c913a91a8e982a98b58ffd9594f Mon Sep 17 00:00:00 2001 From: Niklas Hackelberg Date: Wed, 4 Nov 2026 20:29:21 +0100 Subject: [PATCH 08/11] Docs: custom components Navigation with a further level. "Implement a CT Operator" is now "Write a Projector", under "Custom Components"; "Write a Sampling Strategy" writes a point and a weighting strategy and lets the sampling projector of the first page use them. The point figure shows the subcells. Co-Authored-By: Claude Opus 5.5 --- docs/literate/explanation/masks.jl | 4 +- docs/literate/howto/write_projector.jl | 2 +- docs/literate/howto/write_sampling.jl | 143 +++++++++++++++++++++++++ docs/make.jl | 45 +++++--- docs/src/api/internals.md | 8 +- docs/src/index.md | 2 +- docs/utils/multiray_figures.jl | 13 ++- docs/utils/sampling_figures.jl | 2 +- 8 files changed, 192 insertions(+), 27 deletions(-) create mode 100644 docs/literate/howto/write_sampling.jl diff --git a/docs/literate/explanation/masks.jl b/docs/literate/explanation/masks.jl index 2dd326a..c540257 100644 --- a/docs/literate/explanation/masks.jl +++ b/docs/literate/explanation/masks.jl @@ -60,7 +60,7 @@ nothing #hide # # 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)): +# outside this package (see [Write a Projector](@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 @@ -71,4 +71,4 @@ nothing #hide # 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. +# masked parts. [Write a Projector](@ref) shows how a projector uses the masks itself. diff --git a/docs/literate/howto/write_projector.jl b/docs/literate/howto/write_projector.jl index 7b69d8c..4201c10 100644 --- a/docs/literate/howto/write_projector.jl +++ b/docs/literate/howto/write_projector.jl @@ -1,4 +1,4 @@ -# # Implement a CT Operator +# # Write a Projector # # This guide shows what a new projector has to provide and what it gets in return, by implementing a # simple one: it samples the image at equally spaced points along each ray, taking the value of the diff --git a/docs/literate/howto/write_sampling.jl b/docs/literate/howto/write_sampling.jl new file mode 100644 index 0000000..72e6fda --- /dev/null +++ b/docs/literate/howto/write_sampling.jl @@ -0,0 +1,143 @@ +# # Write a Sampling Strategy +# +# [`MultiRay`](@ref) places several rays in a detector cell according to a point strategy and weights +# them with a weighting strategy. Both are ordinary types, and new ones do not have to be part of +# CTOperators.jl: a script or a package of your own can define them, and they fit in through multiple +# dispatch. This guide writes one of each, and then lets the sampling projector of +# [Write a Projector](@ref) use them. + +using CTOperators, KernelAbstractions, LinearAlgebra +using CTOperators: Gather, Scatter, VoxelGrid, clip, to_index, trace_element, visit +include(joinpath(pkgdir(CTOperators), "docs", "literate", "howto", "write_projector.jl")) #hide +include(joinpath(pkgdir(CTOperators), "docs", "utils", "multiray_figures.jl")) #hide + +# ## A point strategy +# +# A point strategy is a subtype of [`AbstractPoints`](@ref)`{K}`, with `K` the number of rays per +# cell. It implements [`CTOperators.point_offset`](@ref), which returns the offset of ray `k` from +# the cell center in cell units, where ±1/2 are the cell edges. As an example, `K` points on a circle +# around the cell center, rotated by a random angle in every cell: + +struct RingPoints{K,T} <: AbstractPoints{K} + radius::T + seed::UInt32 +end + +RingPoints(K, radius; seed=0) = RingPoints{K,float(typeof(radius))}(radius, UInt32(seed)) + +@inline function CTOperators.point_offset(p::RingPoints{K}, ::Type{T}, iu, iv, view, k) where {K,T} + turn, _ = CTOperators.random_pair(T, p.seed, iu, iv, view, 1) + s, c = sincos(2 * T(pi) * (turn + T(k) / K)) + T(p.radius) * c, T(p.radius) * s +end + +# [`CTOperators.random_pair`](@ref) gives random numbers that are fixed by the seed, the cell, the +# view and a number of your choice, here always 1, so that all rays of a cell share the same rotation. +# They are the same in the forward and the back projection and on every device; [Several Rays per +# Cell](@ref) explains why a generator with state does not work here. +# +# `point_offset` runs inside the kernels, so it must not allocate. Since `RingPoints` stores a +# floating-point number, it also says how to convert it, so that it follows the element type of the +# projector (some GPUs have no Float64): + +CTOperators.points_eltype(p::RingPoints{K}, ::Type{T}) where {K,T} = RingPoints{K,T}(p.radius, p.seed) + +# ## A weighting strategy +# +# A weighting strategy is a subtype of [`AbstractWeights`](@ref) and implements +# [`CTOperators.subray_weight`](@ref), the positive weight of a ray at an offset. The weights of a +# cell are normalized in the kernel. A response that falls off linearly from the cell center: + +struct LinearWeights <: AbstractWeights end + +CTOperators.subray_weight(::LinearWeights, du, dv) = (1 - abs(du)) * (1 - abs(dv)) + +# ## Using them +# +# The new types work like those of the package, with every projector that supports `MultiRay`. The +# points in one cell, with the marker area showing the weights: + +fig = points_figure([("RingPoints(8, 0.3)", RingPoints(8, 0.3; seed=1), UniformWeights()), #hide + ("RingPoints(8, 0.3; seed=2),\nLinearWeights()", RingPoints(8, 0.3; seed=2), LinearWeights())]) #hide +save("ring_points.png", fig) #hide +nothing #hide +# ![Ring points in one cell](ring_points.png) + +g = MultiRay(geometry, RingPoints(8, 0.3; seed=1); weights=LinearWeights()) +op = JosephProjector(image, g) +x, y = rand(Float32, imagesize(op)), rand(Float32, datasize(op)) +dot(op * x, y) / dot(x, adjoint(op) * y) + +# ## A projector with several rays per cell +# +# The sampling projector of [Write a Projector](@ref) asks every ray for its end points with `ray`, so +# it has one ray per data element. A `MultiRay` has no such ray. To work on every geometry, a projector +# moves the code for one ray into a function and traces through +# [`CTOperators.trace_element`](@ref), which calls it once for a plain geometry and once per ray of +# the cell for a `MultiRay`. +# +# The function takes the two end points. [`CTOperators.visit`](@ref) does what the two kernels did +# so far: with a [`CTOperators.Gather`](@ref) it sums the voxels, with a +# [`CTOperators.Scatter`](@ref) it adds to them. The option `samples` travels in a small struct, +# since a closure in a kernel could not be compiled for GPUs: + +struct SamplingSegment{T} + samples::T +end + +@inline function (s::SamplingSegment)(action, grid::VoxelGrid{T}, a, b, rt) where {T} + p = to_index(grid, a) + d = (b - a) ./ grid.deltas + amin, amax = clip(p, d, grid.dims, rt) + n = amin < amax ? ceil(Int, s.samples * norm(d) * (amax - amin)) : 0 + dalpha = n > 0 ? (amax - amin) / n : zero(T) + w = norm(b - a) * dalpha + acc = zero(T) + for j in 1:n + acc = visit(action, acc, nearest_voxel(grid, p + (amin + (j - T(0.5)) * dalpha) * d), w) + end + acc +end + +# The kernels shrink to one call each: + +@kernel function traced_forward!(y, @Const(x), grid, geometry, samples) + i = @index(Global, Linear) + y[i] = trace_element(SamplingSegment(samples), Gather(x), grid, geometry, i) +end + +@kernel function traced_backward!(x, @Const(y), grid, geometry, samples) + i = @index(Global, Linear) + trace_element(SamplingSegment(samples), Scatter(x, y[i]), grid, geometry, i) +end + +function CTOperators.forward!(y, op::SamplingProjector, x) + traced_forward!(op.backend)(y, vec(x), op.grid, op.geometry, op.samples; ndrange=nrays(op.geometry)) + y +end + +function CTOperators.backward!(x, op::SamplingProjector, y) + fill!(x, zero(eltype(x))) + traced_backward!(op.backend)(vec(x), y, op.grid, op.geometry, op.samples; ndrange=nrays(op.geometry)) + x +end + +# The sampling projector now takes every strategy, also the ones above, and its adjoint stays exact: + +op = SamplingProjector(image, MultiRay(geometry, RingPoints(8, 0.3; seed=1); weights=LinearWeights())) +dot(op * x, y) / dot(x, adjoint(op) * y) + +# On plain geometries it computes what it did before. For 2D images, which have one layer of voxels +# along `z`, a projector's constructor should pass the geometry through +# [`CTOperators.in_plane`](@ref), as the projectors of the package do; it keeps the rays in the image +# plane. +# +# ## What fits where +# +# | to change | subtype | implement | +# |:----------|:--------|:----------| +# | where the rays of a cell end | [`AbstractPoints`](@ref) | [`CTOperators.point_offset`](@ref), [`CTOperators.points_eltype`](@ref) if it stores floats | +# | how they are weighted | [`AbstractWeights`](@ref) | [`CTOperators.subray_weight`](@ref), [`CTOperators.weights_eltype`](@ref) if it stores floats | +# | how a ray meets the voxels | [`AbstractProjector`](@ref) | see [Write a Projector](@ref), tracing through [`CTOperators.trace_element`](@ref) | +# +# These functions are declared `public`: they are not exported, but they are part of the interface. diff --git a/docs/make.jl b/docs/make.jl index ebbe9d5..eccbde5 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -38,23 +38,38 @@ makedocs(; "Tomosynthesis" => "generated/tutorials/tomosynthesis.md", ], "How to" => [ - "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", - "Use Several Rays per Cell" => "generated/howto/multiray.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", + "Geometries" => [ + "Set Up a Geometry" => "generated/howto/geometries.md", + "Use Geometries of Sinograms.jl" => "generated/howto/sinograms.md", + ], + "Reconstruction" => [ + "Ordered Subsets and Normalization" => "generated/howto/subsets.md", + "Use Masks" => "generated/howto/masks.md", + "Use Several Rays per Cell" => "generated/howto/multiray.md", + "Differentiate through a Projector" => "generated/howto/differentiation.md", + ], + "Devices" => [ + "Run on a GPU" => "generated/howto/gpu.md", + ], + "Custom Components" => [ + "Write a Projector" => "generated/howto/write_projector.md", + "Write a Sampling Strategy" => "generated/howto/write_sampling.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", - "Masks" => "generated/explanation/masks.md", - "Several Rays per Cell" => "generated/explanation/multiray.md", - "Workgroups and Ray Order on GPUs" => "generated/explanation/launch.md", + "Scans" => [ + "The Parameters of a Scan" => "generated/explanation/scan_parameters.md", + "Coordinates and Conventions" => "generated/explanation/conventions.md", + ], + "Projection Methods" => [ + "Siddon's Method" => "generated/explanation/siddon.md", + "Joseph's Method" => "generated/explanation/joseph.md", + "Several Rays per Cell" => "generated/explanation/multiray.md", + ], + "Operators" => [ + "Masks" => "generated/explanation/masks.md", + "Workgroups and Ray Order on GPUs" => "generated/explanation/launch.md", + ], ], "Reference" => [ "Images" => "api/images.md", diff --git a/docs/src/api/internals.md b/docs/src/api/internals.md index 5534c2b..130c78f 100644 --- a/docs/src/api/internals.md +++ b/docs/src/api/internals.md @@ -1,8 +1,8 @@ # Internals -For developers of operators (see [Implement a CT Operator](@ref)): the functions each operator +For developers of operators (see [Write a Projector](@ref)): the functions each operator implements, the masks a projector stores, the kernel form of the image grid and the clipping of -rays. None of these is exported. +rays. None of these is exported; those that custom components use are declared `public`. ## Operator implementation @@ -33,15 +33,13 @@ Public = false ```@autodocs Modules = [CTOperators] Pages = ["image/image.jl", "image/voxelgrid.jl"] -Public = false ``` ## Ray tracing ```@autodocs Modules = [CTOperators] -Pages = ["raytracing/raytracing.jl", "raytracing/clip.jl", "raytracing/actions.jl"] -Public = false +Pages = ["raytracing/raytracing.jl", "raytracing/clip.jl", "raytracing/actions.jl", "raytracing/subrays.jl"] ``` ## Siddon and Joseph diff --git a/docs/src/index.md b/docs/src/index.md index 592c182..4e658fe 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -55,7 +55,7 @@ y = op * rand(Float32, 128, 128, 128) # projections, size (256, 128, 360) [use several rays per detector cell](@ref "Use Several Rays per Cell"), 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"). + [write your own projector](@ref "Write a Projector") or [sampling strategy](@ref "Write a Sampling Strategy"). - 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), diff --git a/docs/utils/multiray_figures.jl b/docs/utils/multiray_figures.jl index df7f6db..053f95f 100644 --- a/docs/utils/multiray_figures.jl +++ b/docs/utils/multiray_figures.jl @@ -41,14 +41,23 @@ function rays_figure(; n=8, du=2.0, k=4, size=(820, 470)) fig end -# The points of several strategies in one cell (-1/2 to 1/2 in u and v), the marker area -# proportional to the weight. `strategies` holds (title, points, weights) tuples. +# The points of several strategies in one cell (-1/2 to 1/2 in u and v), with the subcells dotted +# and the marker area proportional to the weight. `strategies` holds (title, points, weights) tuples. function points_figure(strategies; size=(220 * length(strategies), 270)) fig = Figure(; size) for (col, (title, points, weights)) in enumerate(strategies) ax = Axis(fig[1, col]; aspect=DataAspect(), title, limits=(-0.9, 0.9, -0.9, 0.9)) hidedecorations!(ax) poly!(ax, Rect2f(-0.5, -0.5, 1, 1); color=(:steelblue, 0.12), strokecolor=:steelblue, strokewidth=1.5) + # the subcells of the strategies that have them + if hasproperty(points, :su) + for j in 1:points.su-1 + lines!(ax, fill(j / points.su - 0.5, 2), [-0.5, 0.5]; color=:steelblue, linestyle=:dot, linewidth=1) + end + for j in 1:points.sv-1 + lines!(ax, [-0.5, 0.5], fill(j / points.sv - 0.5, 2); color=:steelblue, linestyle=:dot, linewidth=1) + end + end K = typeof(points).parameters[1] offs = [point_offset(points, Float64, 3, 2, 1, k) for k in 1:K] w = [subray_weight(weights, o...) for o in offs] diff --git a/docs/utils/sampling_figures.jl b/docs/utils/sampling_figures.jl index d8bb3b5..7562b24 100644 --- a/docs/utils/sampling_figures.jl +++ b/docs/utils/sampling_figures.jl @@ -1,4 +1,4 @@ -# Figure of the guide "Implement a CT operator" (CairoMakie): the sampling projector on a small 2D +# Figure of the guide "Write a Projector" (CairoMakie): the sampling projector on a small 2D # grid. Included invisibly by the page. using CairoMakie From 9834c1011bae120bd255d4e0902ae7bab0a0b398 Mon Sep 17 00:00:00 2001 From: Niklas Hackelberg Date: Wed, 4 Nov 2026 21:20:27 +0100 Subject: [PATCH 09/11] Geometry: one layer of voxels, documented interface Rays stay in the image plane for every image with one layer of voxels along z, also a 3D image with nz = 1, which was treated like a volume. cell_ray and is_masked are documented and public; the docstrings say what custom strategies must provide. Co-Authored-By: Claude Opus 5.5 --- src/CTOperators.jl | 4 +- src/geometry/geometry.jl | 17 ++++--- src/geometry/multiray.jl | 86 +++++++++++++++++++--------------- src/operators/joseph/joseph.jl | 6 +-- src/operators/masks.jl | 7 ++- src/operators/siddon/siddon.jl | 6 +-- src/raytracing/actions.jl | 2 + src/raytracing/subrays.jl | 15 +++--- test/operators/multiray.jl | 5 ++ 9 files changed, 89 insertions(+), 59 deletions(-) diff --git a/src/CTOperators.jl b/src/CTOperators.jl index 3a6125d..13e504b 100644 --- a/src/CTOperators.jl +++ b/src/CTOperators.jl @@ -19,10 +19,10 @@ export datasize, nrays, ray, raytype, nviews, detectorsize, subset, projection_g export AbstractProjector, AdjointCTOperator, SiddonProjector, JosephProjector, LaunchConfig export geometry, image_geometry, imagesize, rowsums, colsums -# not exported, but stable: the interface for types of other packages +# not exported, but part of the interface: for projectors, geometries and strategies of other packages public forward!, backward!, similar_projector, storage_array, device_geometry public point_offset, random_pair, points_eltype, subray_weight, weights_eltype -public trace_element, in_plane, visit, Gather, Scatter, to_index, clip, VoxelGrid +public trace_element, in_plane, one_layer, cell_ray, is_masked, visit, Gather, Scatter, to_index, clip, VoxelGrid public image_mask, ray_mask, handles_raymask, projector_masks, ProjectorMasks """ diff --git a/src/geometry/geometry.jl b/src/geometry/geometry.jl index f4918d2..2f166ed 100644 --- a/src/geometry/geometry.jl +++ b/src/geometry/geometry.jl @@ -28,6 +28,7 @@ Interface (in addition to that of `AbstractProjectionGeometry`): * `nviews(g)`, `detectorsize(g) = (nu, nv)` * `subset(g, views)`: the geometry of a subset of the views +* [`cell_ray`](@ref CTOperators.cell_ray)`(g, cu, cv, k)`, for several rays per cell """ abstract type AbstractDetectorGeometry{T} <: AbstractProjectionGeometry{T} end @@ -120,12 +121,16 @@ end c + (cu - T(nu + 1) / 2) * u + (cv - T(nv + 1) / 2) * v end -# cell_ray(g, cu, cv, k) -> (a, b): the ray of view k of the detector geometry g that ends at -# position (cu, cv) on the detector, as two points like `ray` (cone beams: source and detector -# point; parallel beams: detector point and that point plus the direction). Positions are in cell -# units: integers are cell centers, so `cell_ray(g, iu, iv, k)` is the ray of cell (iu, iv), equal to -# `ray`; half-integers are cell edges, e.g. `cell_ray(g, iu + 1/2, iv, k)` runs through the midpoint -# of the edge between cells iu and iu + 1. Callable inside kernels. +""" + cell_ray(g::AbstractDetectorGeometry, cu, cv, k) -> (a, b) + +The ray of view `k` that ends at position `(cu, cv)` on the detector, as two points like +[`ray`](@ref): for cone beams the source and the detector point, for parallel beams the detector +point and that point plus the direction. Positions are in cell units: integers are cell centers, so +`cell_ray(g, iu, iv, k)` is the ray of cell `(iu, iv)`; half-integers are cell edges, e.g. +`cell_ray(g, iu + 1/2, iv, k)` runs through the middle of the edge between cells `iu` and `iu + 1`. +Callable inside kernels. A [`MultiRay`](@ref) needs it from its geometry. +""" function cell_ray end function check_views(name, ms...) diff --git a/src/geometry/multiray.jl b/src/geometry/multiray.jl index 381903a..7e1a223 100644 --- a/src/geometry/multiray.jl +++ b/src/geometry/multiray.jl @@ -1,7 +1,7 @@ # Several rays per detector cell: a detector geometry wrapped with a point strategy (where in the -# cell the subrays end) and a weighting strategy. Offsets are in cell units: (0, 0) is the cell -# center, +-1/2 the cell edges. The number of subrays K is a type parameter of the point strategy, -# so the loop over the subrays has a length known to the compiler. +# cell the rays end) and a weighting strategy. Offsets are in cell units: (0, 0) is the cell +# center, +-1/2 the cell edges. The number of rays K is a type parameter of the point strategy, +# so the loop over the rays has a length known to the compiler. """ AbstractPoints{K} @@ -25,8 +25,8 @@ abstract type AbstractWeights end """ RegularPoints(su, sv=su) -The centers of `su × sv` subcells of each detector cell: `su` along `u`, `sv` along `v`. Point -strategy of [`MultiRay`](@ref). +The centers of `su × sv` subcells of each detector cell, `su` along `u` and `sv` along `v`: +`su * sv` rays per cell. Point strategy of [`MultiRay`](@ref). """ struct RegularPoints{K} <: AbstractPoints{K} su::Int @@ -38,8 +38,8 @@ end The points of [`RegularPoints`](@ref), each moved to a random position within its subcell ("jittered sampling", R. L. Cook, "Stochastic sampling in computer graphics", ACM Transactions on -Graphics 5(1), 51-72, 1986). The positions are fixed by `seed`, the detector cell, the view and the -subray, so they are the same in the forward and the back projection, on every device and in +Graphics 5(1), 51-72, 1986). The offsets are fixed by `seed` (default 0), the detector cell, the view +and the ray, so they are the same in the forward and the back projection, on every device and in `subset`. Point strategy of [`MultiRay`](@ref). """ struct JitteredPoints{K} <: AbstractPoints{K} @@ -51,8 +51,8 @@ end """ RandomPoints(n; seed=0) -`n` random points in each detector cell, fixed like those of [`JitteredPoints`](@ref). Point -strategy of [`MultiRay`](@ref). +`n` random points, uniformly distributed in each detector cell, fixed like those of +[`JitteredPoints`](@ref). Point strategy of [`MultiRay`](@ref). """ struct RandomPoints{K} <: AbstractPoints{K} seed::UInt32 @@ -64,7 +64,7 @@ end `n` random points per detector cell, normally distributed around the cell center with standard deviations `sigma_u`, `sigma_v` in cell units: a detector response set by where the rays land, with all rays counting equally. Points can land outside the cell, which models blur across neighboring -cells. The positions are fixed like those of [`JitteredPoints`](@ref). Point strategy of +cells. The offsets are fixed like those of [`JitteredPoints`](@ref). Point strategy of [`MultiRay`](@ref). """ struct GaussianPoints{K,T} <: AbstractPoints{K} @@ -77,7 +77,7 @@ end Points(offsets) Given points, the same in every detector cell: column `k` of the `2 × K` matrix `offsets` is the -position of subray `k` in cell units, `(0, 0)` the cell center and `±1/2` the cell edges. Points +offset of ray `k` in cell units, `(0, 0)` the cell center and `±1/2` the cell edges. Points outside the cell are allowed. Point strategy of [`MultiRay`](@ref). """ struct Points{K,T} <: AbstractPoints{K} @@ -133,9 +133,11 @@ points_eltype(p::Points{K}, ::Type{T}) where {K,T} = Points{K,T}(map(o -> SVecto """ point_offset(points, T, iu, iv, view, k) -> (du, dv) -The offset of ray `k` from the center of the detector cell `(iu, iv)` in view `view`, in element type -`T` and in cell units (`±1/2` are the cell edges). `view` is numbered as in the full geometry, so -that subsets draw the same random positions. Called in kernels: no allocations. +The offset of ray `k` (from `1` to `K`) from the center of the detector cell `(iu, iv)` in view +`view`, in element type `T` and in cell units (`±1/2` are the cell edges). `view` is numbered as in +the full geometry, so that subsets draw the same random offsets. Called in kernels: no allocations. +Fields of a strategy must be plain values: they reach the kernels unchanged, so an array field would +stay on the CPU. """ function point_offset end @@ -175,7 +177,8 @@ end random_pair(T, seed, iu, iv, view, k) -> (ru, rv) Two uniform random numbers in `[0, 1)`, fixed by the arguments: the same on every device and in -every call, and independent for different cells, views, rays and seeds. For random point strategies. +every call, and independent for different arguments. `seed` is a `UInt32`, `(iu, iv)` the detector +cell, `view` the view; `k` is any number, usually the ray of the cell. For random point strategies. """ @inline function random_pair(::Type{T}, seed, iu, iv, view, k) where {T} # the cell as one number keeps the four inputs distinct for up to 65536 cells per row @@ -187,7 +190,7 @@ end """ UniformWeights() -All subrays of a detector cell count equally. Weighting strategy of [`MultiRay`](@ref). +All rays of a detector cell count equally. Weighting strategy of [`MultiRay`](@ref). """ struct UniformWeights <: AbstractWeights end @@ -222,16 +225,16 @@ weights_eltype(w::GaussianWeights, ::Type{T}) where {T} = GaussianWeights{T}(w.s """ subray_weight(weights, du, dv) -The weight of a ray at offset `(du, dv)` from the cell center, positive and not normalized (the -weights of a cell are normalized in the kernel). Called in kernels: no allocations. +The weight of a ray at offset `(du, dv)` from the cell center, not normalized: the kernel divides by +the sum of the weights of the cell, which must not be zero. Called in kernels: no allocations. """ function subray_weight end @inline subray_weight(::UniformWeights, du, dv) = one(du) @inline subray_weight(w::GaussianWeights, du, dv) = exp(-((du / w.sigma_u)^2 + (dv / w.sigma_v)^2) / 2) -# Whether the subrays of a cell move along both detector axes or only along u. For 2D images (one -# voxel layer along z) the rays must stay in the image plane, which the projectors ensure with +# Whether the rays of a cell move along both detector axes or only along u. For images with one +# voxel layer along z the rays must stay in the image plane, which the projectors ensure with # `in_plane`. struct FullCell end struct InPlane end @@ -243,24 +246,25 @@ struct InPlane end MultiRay(geometry, points; weights=UniformWeights()) The detector geometry `geometry` (`ConeBeam`, `ArcConeBeam`, `ParallelBeam`) with several rays per -detector cell. The rays of a cell end at the positions given by the point strategy `points` +detector cell. The rays of a cell end at the offsets given by the point strategy `points` ([`RegularPoints`](@ref), [`JitteredPoints`](@ref), [`RandomPoints`](@ref), -[`GaussianPoints`](@ref), [`Points`](@ref), or an own [`AbstractPoints`](@ref)) and +[`GaussianPoints`](@ref), [`Points`](@ref), or a custom [`AbstractPoints`](@ref)) and are averaged with the weights of `weights` ([`UniformWeights`](@ref), -[`GaussianWeights`](@ref), or an own [`AbstractWeights`](@ref)), normalized to sum 1 per cell. A ray-driven projector built on it +[`GaussianWeights`](@ref), or a custom [`AbstractWeights`](@ref)), normalized to sum 1 per cell. A ray-driven projector built on it computes - y[i] = sum_k w_k * (line integral along subray k of data element i) + y[i] = sum_k w_k * (line integral along ray k of data element i) and the exact adjoint. This averages the line integrals over the cell (closer to what a detector -cell measures than one line through its center) and removes the gaps between rays that occur when -detector cells are large against the voxels. The data has the size of `geometry`. +cell measures than one line through its center) and reduces the gaps between rays that occur when +detector cells are large compared to the voxels. The data have the size of `geometry`. -For 2D images (one voxel layer along `z`) the projectors keep the subrays in the image plane: they -are spread along `u` only, and point strategies with several subcells along `v` are rejected. +For images with one layer of voxels along `z` (2D images, or 3D images with `nz = 1`) the projectors +keep the rays in the image plane: they are spread along `u` only, and point strategies with several +subcells along `v` are rejected. -A `MultiRay` has no single ray per data element (no [`ray`](@ref)); projectors that trace only one -ray per cell reject it. `subset` and `adapt` keep the points and weights; random positions stay the +A `MultiRay` has no single ray per data element (no [`ray`](@ref)): a projector that traces only one +ray per cell fails with a `MethodError` when applied, instead of ignoring the other rays. `subset` and `adapt` keep the points and weights; random offsets stay the same for the views of a subset. # Example @@ -311,8 +315,8 @@ Adapt.adapt_structure(to, g::MultiRay) = in_plane(geometry) `geometry` with the rays of a [`MultiRay`](@ref) in the image plane (offsets along `u` only); other -geometries unchanged. The constructor of a projector calls it for 2D images, which have one layer -of voxels along `z`. +geometries unchanged. The constructor of a projector calls it for images with one layer of voxels +along `z`, see [`CTOperators.one_layer`](@ref). """ in_plane(g) = g function in_plane(g::MultiRay) @@ -321,8 +325,16 @@ function in_plane(g::MultiRay) end # subcells along v would all fall onto v = 0 and trace the same rays several times +""" + one_layer(image::ImageGeom) + +Whether `image` has one layer of voxels along `z`: a 2D image, or a 3D image with `nz = 1`. Rays that +leave the center of the layer see too little of it, so the rays of a cell stay in its plane. +""" +one_layer(image::ImageGeom) = length(image.dims) == 2 || image.dims[3] == 1 + check_in_plane(p::Union{RegularPoints,JitteredPoints}) = - p.sv == 1 || throw(ArgumentError("a 2D image has one row of rays per view: use $(nameof(typeof(p)))(su, 1)")) + p.sv == 1 || throw(ArgumentError("an image with one layer of voxels needs one row of rays: use $(nameof(typeof(p)))(su, 1)")) check_in_plane(p) = nothing function Base.show(io::IO, ::MIME"text/plain", g::MultiRay) @@ -343,19 +355,19 @@ end iu, iv, view, full_view(g.views, view) end -# the offset of subray k in the cell +# the offset of ray k in the cell @inline function subray_offset(g::MultiRay{T}, cell, k) where {T} iu, iv, _, number = cell restrict(g.span, point_offset(g.points, T, iu, iv, number, k)...) end -# the weight (not normalized) of subray k in the cell +# the weight (not normalized) of ray k in the cell @inline function cell_weight(g::MultiRay, cell, k) du, dv = subray_offset(g, cell, k) subray_weight(g.weights, du, dv) end -# subray k in the cell +# ray k in the cell @inline function cell_subray(g::MultiRay{T}, cell, k) where {T} iu, iv, view, _ = cell du, dv = subray_offset(g, cell, k) @@ -366,7 +378,7 @@ end """ subray(g::MultiRay, i, k) -> (a, b, w) -Subray `k` of data element `i` as two points like [`ray`](@ref), and its weight (not normalized). +Ray `k` of data element `i` as two points like [`ray`](@ref), and its weight (not normalized). Callable inside kernels. """ @inline subray(g::MultiRay, i, k) = cell_subray(g, element_cell(g, i), k) diff --git a/src/operators/joseph/joseph.jl b/src/operators/joseph/joseph.jl index 417a4ab..bd81ec2 100644 --- a/src/operators/joseph/joseph.jl +++ b/src/operators/joseph/joseph.jl @@ -67,10 +67,10 @@ 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; for 2D images the -# subrays of a MultiRay stay in the image plane +# from a geometry already in the operator's element type and on its device; for images of one layer +# the rays of a MultiRay stay in the image plane function joseph_projector(image::ImageGeom, geometry::AbstractProjectionGeometry{T}, launch, raymask) where {T} - g = length(image.dims) == 2 ? in_plane(geometry) : geometry + g = one_layer(image) ? in_plane(geometry) : geometry backend = KA.get_backend(storage_array(g)) JosephProjector(image, VoxelGrid{T}(image), g, backend, launch, ProjectorMasks(image, g, raymask)) end diff --git a/src/operators/masks.jl b/src/operators/masks.jl index 9ea6080..02e9281 100644 --- a/src/operators/masks.jl +++ b/src/operators/masks.jl @@ -85,6 +85,11 @@ function masked_backward!(x, op, y) x end -# is ray i masked out? Inside kernels; without a mask the check compiles away +""" + is_masked(raymask, i) + +Whether ray `i` is masked, with `raymask` from [`CTOperators.ray_mask`](@ref). Called in kernels; +without a mask (`nothing`) the check compiles away. +""" @inline is_masked(::Nothing, i) = false @inline is_masked(mask, i) = @inbounds !mask[i] diff --git a/src/operators/siddon/siddon.jl b/src/operators/siddon/siddon.jl index a46560c..e8d81f2 100644 --- a/src/operators/siddon/siddon.jl +++ b/src/operators/siddon/siddon.jl @@ -57,10 +57,10 @@ 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; for 2D images the -# subrays of a MultiRay stay in the image plane +# from a geometry already in the operator's element type and on its device; for images of one layer +# the rays of a MultiRay stay in the image plane function siddon_projector(image::ImageGeom, geometry::AbstractProjectionGeometry{T}, launch, raymask) where {T} - g = length(image.dims) == 2 ? in_plane(geometry) : geometry + g = one_layer(image) ? in_plane(geometry) : geometry backend = KA.get_backend(storage_array(g)) SiddonProjector(image, VoxelGrid{T}(image), g, backend, launch, ProjectorMasks(image, g, raymask)) end diff --git a/src/raytracing/actions.jl b/src/raytracing/actions.jl index 9af5d9c..1d95aa0 100644 --- a/src/raytracing/actions.jl +++ b/src/raytracing/actions.jl @@ -2,6 +2,7 @@ Gather(image) What a ray does in the forward projection: it sums the weighted values of the voxels it visits. +`image` is indexed linearly. """ struct Gather{A} image::A @@ -11,6 +12,7 @@ end Scatter(image, value) What a ray does in the back projection: it adds its weighted `value` to the voxels it visits. +`image` is indexed linearly. """ struct Scatter{A,T} image::A diff --git a/src/raytracing/subrays.jl b/src/raytracing/subrays.jl index 632c24a..84a9f1a 100644 --- a/src/raytracing/subrays.jl +++ b/src/raytracing/subrays.jl @@ -1,16 +1,17 @@ """ trace_element(segment, action, grid, geometry, i) -Traces the rays of data element `i` for a ray-driven method: one ray for plain geometries, the -weighted rays of a cell for a [`MultiRay`](@ref). `segment(action, grid, a, b, raytype)` traces one -ray from `a` to `b`, calling [`CTOperators.visit`](@ref) for every voxel, and returns the value of -the action. The result is the value of data element `i` for a [`CTOperators.Gather`](@ref). A +Traces the rays of data element `i` for a ray-driven method: one ray for plain geometries, the rays +of a cell with normalized weights for a [`MultiRay`](@ref). `segment(action, grid, a, b, raytype)` +traces one ray from `a` to `b`, calling [`CTOperators.visit`](@ref) for every voxel, and returns the +accumulated value of the action. The result is the value of data element `i` for a +[`CTOperators.Gather`](@ref). A projector that traces through this function works on every geometry. Called in kernels. """ @inline trace_element(segment::F, action, grid, geometry, i) where {F} = segment(action, grid, ray(geometry, i)..., raytype(geometry)) -# forward: the weighted mean of the subrays' line integrals +# forward: the weighted mean of the line integrals of the rays @inline function trace_element(segment::F, action::Gather, grid::VoxelGrid{T}, g::MultiRay, i) where {F,T} cell = element_cell(g, i) acc, wsum = zero(T), zero(T) @@ -22,8 +23,8 @@ projector that traces through this function works on every geometry. Called in k acc / wsum end -# back: every subray adds its share w / sum(w) of the value; the weights are summed first, so that -# the shares are known before any subray is traced +# back: every ray adds its share w / sum(w) of the value; the weights are summed first, so that the +# shares are known before any ray is traced @inline function trace_element(segment::F, action::Scatter, grid::VoxelGrid{T}, g::MultiRay, i) where {F,T} cell = element_cell(g, i) wsum = zero(T) diff --git a/test/operators/multiray.jl b/test/operators/multiray.jl index a10473b..a72e86f 100644 --- a/test/operators/multiray.jl +++ b/test/operators/multiray.jl @@ -116,6 +116,11 @@ end end # subcells along v would collapse onto the plane @test_throws ArgumentError SiddonProjector(image2, MultiRay(fan(:flat), RegularPoints(2, 2))) + # a 3D image of one layer is treated like a 2D image + layer = ImageGeom(; dims=(24, 20, 1), deltas=(1.0, 1.0, 1.0)) + @test geometry(SiddonProjector(layer, MultiRay(fan(:flat), pts))).span isa InPlane + @test isapprox(JosephProjector(layer, MultiRay(fan(:flat), pts); eltype=Float64) * reshape(x, 24, 20, 1), + JosephProjector(image2, MultiRay(fan(:flat), pts); eltype=Float64) * x) end # Strategies are open: types of other packages only subtype AbstractPoints and AbstractWeights and From 65697aa79260f4f6a194b8a881e4f63ceadaaa96 Mon Sep 17 00:00:00 2001 From: Niklas Hackelberg Date: Wed, 4 Nov 2026 21:20:27 +0100 Subject: [PATCH 10/11] Docs: comparing projectors, review of the multi-ray pages An explanation page with the back projection of one view and of all views, for Siddon and Joseph with one and with several rays per cell and for detector cells smaller than, equal to and larger than the pixels. Corrections from a cold review: which strategy causes which error, metrics for the coverage, what the figures encode, and the documented conditions of custom strategies. Co-Authored-By: Claude Opus 5.5 --- docs/literate/explanation/comparing.jl | 85 ++++++++++++++++++ docs/literate/explanation/multiray.jl | 40 +++++---- docs/literate/howto/multiray.jl | 117 +++++++++++++++---------- docs/literate/howto/write_projector.jl | 11 ++- docs/literate/howto/write_sampling.jl | 75 +++++++++------- docs/make.jl | 3 + docs/src/api/geometries.md | 3 +- docs/src/api/internals.md | 9 +- docs/src/api/projectors.md | 4 +- docs/src/index.md | 3 +- docs/utils/comparison_figures.jl | 59 +++++++++++++ docs/utils/multiray_figures.jl | 5 +- 12 files changed, 306 insertions(+), 108 deletions(-) create mode 100644 docs/literate/explanation/comparing.jl create mode 100644 docs/utils/comparison_figures.jl diff --git a/docs/literate/explanation/comparing.jl b/docs/literate/explanation/comparing.jl new file mode 100644 index 0000000..c4068de --- /dev/null +++ b/docs/literate/explanation/comparing.jl @@ -0,0 +1,85 @@ +# # Comparing Projectors +# +# How evenly a projector covers the image depends mainly on the ratio of detector cell size to voxel +# size and on the number of rays per cell, less on the projector. This page compares +# [Siddon's Method](@ref) and [Joseph's Method](@ref) with one and with several rays per cell (see +# [Use Several Rays per Cell](@ref)) for cells smaller than, equal to and larger than the pixels of a +# 2D image, with two checks from Miao (2013): the back projection of one view and of all views of +# ones. It compares coverage, not accuracy on real data or speed. + +using CTOperators +include(joinpath(pkgdir(CTOperators), "docs", "utils", "comparison_figures.jl")) #hide + +models = [("Siddon", SiddonProjector, identity), ("Joseph", JosephProjector, identity), + ("Siddon, 4 rays", SiddonProjector, g -> MultiRay(g, RegularPoints(4, 1))), + ("Joseph, 4 rays", JosephProjector, g -> MultiRay(g, RegularPoints(4, 1)))] +ratios = (0.5, 1.0, 2.0); + +# ## One view +# +# The back projection of a single parallel-beam view of ones, at an angle of 0.3 radians, shows how +# the rays of a view cover the pixels: every pixel gets the summed weights of the rays through it. +# +# This coverage should be the same for every pixel. If rays covered the detector continuously, then +# adding up the lengths of all rays inside a pixel would slice the pixel into thin parallel strips +# and sum them: the result is the area of the pixel, at any angle. One ray per cell samples this +# integral with the spacing `du` of the cells, so the sum times `du` should be 1, and any other value +# is an error of the sampling. The figure shows this scaled coverage for the 16 × 16 central pixels, +# from 0 (red) over even coverage 1 (white) to 2 (blue); rows without a number of rays have one ray +# per cell. The angle 0.3 avoids rays parallel to the pixel grid, which hide the effect for cells up +# to the size of a pixel: + +fig = view_figure(models, ratios) #hide +save("comparison_view.png", fig) #hide +nothing #hide +# ![Back projection of one view](comparison_view.png) + +# With cells as large as the pixels or smaller, every method covers the pixels nearly evenly. With +# cells of 2 pixels, one ray per cell misses pixels: Siddon's method leaves many pixels with little or +# no weight, and Joseph's interpolation softens the pattern without removing it. Four rays per cell +# fill the gaps: their spacing is then that of cells of 0.5 pixels. +# +# ## All views +# +# The back projection of all 180 views of ones (the "pixel support", `colsums`), scaled in the same +# way and divided by the number of views, shows what remains after the views are added up. Below each +# panel is the variation, the standard deviation relative to the mean, inside the circle of radius 24 +# pixels that is drawn; the color window is ±20 %: + +fig = support_figure(models, ratios) #hide +save("comparison_support.png", fig) #hide +nothing #hide +# ![Back projection of all views](comparison_support.png) + +# The pattern of the single views does not average out. With cells of 2 pixels, the coverage varies +# by about 10 %, in a cross and a checker pattern around the center; four rays per cell reduce that +# to a fraction of a percent. An iterative reconstruction sees this pattern as part of the model, and +# it shows up as artifacts in the images, see [Use Several Rays per Cell](@ref). +# +# ## How many rays +# +# The variation for cells of 2 pixels and regular points, for several numbers of rays per cell: + +coverage(P, k) = 2 .* colsums(P(comparison_image(64), MultiRay(comparison_scan(64, 2.0, range(0, pi, 181)[1:180]), RegularPoints(k, 1)); eltype=Float64)) ./ 180 +[(rays=k, Siddon=round(variation(coverage(SiddonProjector, k), 24); digits=2), + Joseph=round(variation(coverage(JosephProjector, k), 24); digits=2)) for k in (1, 2, 3, 4, 8)] + +# More rays help, but not steadily: the number matters relative to the pixels, and an odd number of +# rays can be worse than a smaller even one. A spacing of the rays at the object of at most half a +# voxel covers the image well. +# +# ## What to take from it +# +# * **The ratio matters most.** With cells no larger than the voxels, both methods cover the image +# well with one ray per cell, and Joseph's method is the more even of the two. +# * **For cells larger than the voxels,** one ray per cell is not enough for either method; several +# rays per cell are, at the cost of about that many times the computing time (see the +# [benchmarks](https://github.com/JuliaImageRecon/CTOperators.jl/blob/main/benchmark/README.md)). +# Voxels no smaller than the cells avoid the problem as well, at the price of a coarser image. +# * **In a cone beam,** the size to compare with the voxels is the cell size at the object: the cell +# size divided by the magnification. +# +# ## References +# +# - C. Miao, "Comparative studies of different system models for iterative CT image +# reconstruction", master's thesis, Wake Forest University, 2013. diff --git a/docs/literate/explanation/multiray.jl b/docs/literate/explanation/multiray.jl index cc492e4..6282309 100644 --- a/docs/literate/explanation/multiray.jl +++ b/docs/literate/explanation/multiray.jl @@ -13,7 +13,7 @@ include(joinpath(pkgdir(CTOperators), "docs", "utils", "multiray_figures.jl")) # # detector cell as a point: one line from the source to the center of the cell. When the cells are # large compared to the voxels, the lines miss voxels between them, and the back projection covers # the voxels unevenly (De Man & Basu 2004). A real cell integrates over its area. `MultiRay(geometry, points; weights)` replaces the one line by `K` rays -# per cell, ending at the positions given by `points`, and averages their line integrals: +# per cell, ending at the offsets given by `points`, and averages their line integrals: # # ```math # y_i = \sum_{k=1}^K \frac{w_{ik}}{\sum_l w_{il}} \int_{\text{ray } (i, k)} x @@ -36,8 +36,10 @@ include(joinpath(pkgdir(CTOperators), "docs", "utils", "multiray_figures.jl")) # # detector a line through that point. So the rays of a cell sample its actual shape, also on curved # detectors. # -# The point strategies give the offsets of the rays from the cell center, between -1/2 and 1/2; -# the weighting strategies give a weight for every offset. Marker areas show the weights: +# The point strategies give the offsets of the rays from the cell center, where ±1/2 are the cell +# edges; some go beyond them. The weighting strategies give a weight for every offset. In the figure +# the blue square is one cell, the dotted lines its subcells, and the orange points the rays, with +# the marker area relative to the largest weight in the panel: fig = points_figure([("RegularPoints(3, 3)", RegularPoints(3, 3), UniformWeights()), #hide ("JitteredPoints(3, 3)", JitteredPoints(3, 3; seed=1), UniformWeights()), #hide @@ -49,23 +51,25 @@ nothing #hide # ![Point strategies in one cell](multiray_points.png) # * **Regular points** sit at the centers of `su × sv` subcells: the midpoint rule for the integral -# over the cell. For data that are cell averages they are the most accurate with few points. +# over the cell, a sum at equal distances. For data that average over the cell they are the most +# accurate with few points. # * **Jittered points** move each regular point to a random position within its subcell ("jittered # sampling", Cook 1986, also called stratified sampling). They keep the even spread, but neighboring -# cells no longer share the same pattern, so regular structures of the rays do not line up with -# the voxel grid. -# * **Random points** are spread anywhere in the cell. With few points their mean can differ -# noticeably from the cell average. +# cells do not share the same pattern, so the rays do not line up with the voxel grid in a regular +# way. The price is noise in the model. +# * **Random points** are spread anywhere in the cell. With few points their mean can differ noticeably +# from the cell average: with 4 rays per cell the example of [Use Several Rays per Cell](@ref) is +# worse than with one. # * **Gaussian points** are normally distributed around the cell center, and all rays count equally: # the distribution of the positions is the weighting. Points may land outside the cell, which # models blur across neighboring cells. -# * **Given points** (`Points`) are the same in every cell, also outside it. +# * **Given points** (`Points`) are the same in every cell and may lie outside it. # * **Gaussian weights** describe a similar response with fixed positions: the points stay where # they are, and the rays near the edges count less. # -# For 2D images, which have one layer of voxels along `z`, the projectors keep all rays in the image -# plane: the offsets along `v` are set to 0, and strategies with several subcells along `v` are -# rejected. +# For images with one layer of voxels along `z` (2D images, or 3D images with `nz = 1`) the projectors +# keep all rays in the image plane: the offsets along `v` are set to 0, and strategies with several +# subcells along `v` are rejected. A ray away from the center of the layer would see too little of it. # # ## Random positions without a random number generator # @@ -79,9 +83,10 @@ nothing #hide # give the offsets. # # The view is numbered as in the full scan: `subset(op, views)` keeps these numbers, so the rays of -# a subset are those of the same views of the full operator, as ordered-subset methods require. A +# a subset are those of the same views of the full operator, as +# [ordered-subset methods](@ref "Ordered Subsets and Normalization") require. A # different seed gives a different, again fixed, set of positions: `MultiRay(geometry, -# JitteredPoints(4; seed=2))`. +# JitteredPoints(4, 1; seed=2))`. # # ## Implementation # @@ -91,12 +96,13 @@ nothing #hide # * The kernels are those of the plain projectors: each thread handles one data element. For a # `MultiRay` it traces the `K` rays one after the other and adds them up (forward), or first sums # the weights and then lets every ray add its share (back). No memory per ray is needed. -# * The time is about `K` times that of one ray; computing the positions and weights is negligible -# (see the [benchmarks](https://github.com/JuliaImageRecon/CTOperators.jl/blob/main/benchmark/README.md)). +# * The time is about `K` times that of one ray; jittered offsets and Gaussian weights added nothing +# measurable on a GPU (see the +# [benchmarks](https://github.com/JuliaImageRecon/CTOperators.jl/blob/main/benchmark/README.md#Several-rays-per-cell-(GPU))). # * The ray mask applies per data element (a masked element skips all its rays), the clipping to the # bounding box of the image mask per ray. # * A `MultiRay` has no single ray per data element. A projector that traces only one ray per cell -# does not accept it, instead of silently ignoring the extra rays. +# fails with a `MethodError` for `ray` when applied, instead of silently ignoring the other rays. # # ## References # diff --git a/docs/literate/howto/multiray.jl b/docs/literate/howto/multiray.jl index 8848859..a35ee6a 100644 --- a/docs/literate/howto/multiray.jl +++ b/docs/literate/howto/multiray.jl @@ -4,32 +4,35 @@ # ray per cell, through its center. When the cells are large compared to the voxels, these rays # miss voxels between them, and the reconstruction shows patterns that come from the sampling, not # from the object. This guide shows the problem on a small example and how `MultiRay` traces several -# rays per cell instead. [Several Rays per Cell](@ref) explains how the rays are placed and weighted. +# rays per cell instead. [Several Rays per Cell](@ref) explains how the rays are placed and weighted, +# and [Comparing Projectors](@ref) when it matters. using CTOperators, ImagePhantoms, RegularizedLeastSquares, LinearOperatorCollection, CairoMakie using ImageGeoms: axis include(joinpath(pkgdir(CTOperators), "docs", "utils", "tutorial_figures.jl")) #hide include(joinpath(pkgdir(CTOperators), "docs", "utils", "multiray_figures.jl")) #hide +include(joinpath(pkgdir(CTOperators), "docs", "utils", "comparison_figures.jl")) #hide # ## The example # -# A 2D Shepp-Logan phantom of 128 × 128 pixels and a parallel-beam scan with 180 views. The detector -# cells are 2 pixels wide: +# A 2D Shepp-Logan phantom of 128 × 128 pixels and a parallel-beam scan with 180 views (projection +# directions). The detector cells are 2 pixels wide; `u` runs along a detector row, `v` across the +# rows. In a cone beam, the size to compare with the voxels is the cell size at the object, the cell +# size divided by the magnification `dsd / dso`. n = 128 image = ImageGeom(; dims=(n, n), deltas=(1.0, 1.0)) x_true = Float32.(phantom(axis(image, 1), axis(image, 2), shepp_logan(SheppLoganToft(); fovs=(n, n)), 2)) geometry = parallel_beam(; angles=range(0, pi, 181)[1:180], nu=96, nv=1, du=2.0); -# The data should come from cells that average over their width, as a real detector does. They are -# simulated on a 4 times finer image with 16 rays per cell, so that the data do not come from any of -# the models compared below: +# The data average the line integrals over the width of each cell. They are simulated on a 4 times +# finer image with 16 rays per cell, so that they do not come from any of the models compared below: fine = ImageGeom(; dims=(4n, 4n), deltas=(0.25, 0.25)) x_fine = Float32.(phantom(axis(fine, 1), axis(fine, 2), shepp_logan(SheppLoganToft(); fovs=(n, n)), 2)) y = JosephProjector(fine, MultiRay(geometry, RegularPoints(16, 1))) * x_fine; -# The reconstructions use 40 CGNR iterations: +# The reconstructions use 40 iterations of CGNR (conjugate gradients on the normal equations): function reconstruct(op) solver = createLinearSolver(CGNR, op; AHA=normalOperator(op), iterations=40) @@ -38,27 +41,28 @@ end; # ## One ray per cell # -# So far each data element was one ray, through the center of its cell. Positions on the detector -# are counted in cells: cell 2 has its center at position 2 and its edges at 1.5 and 2.5. With cells -# twice as wide as the pixels, the rays of a view pass through every other column of pixels (left): +# Without `MultiRay`, each data element (one number of the projection data: one cell in one view) is +# one ray, through the center of its cell. Positions on the detector are counted in cells: cell 2 has +# its center at position 2 and its edges at 1.5 and 2.5. With cells twice as wide as the pixels, the +# rays of a view miss every other line of pixels (left; gray: pixels, orange: rays, blue: cells): fig = rays_figure() #hide save("multiray_rays.png", fig) #hide nothing #hide # ![Rays through the cell centers and four rays per cell](multiray_rays.png) -# Which pixels the rays miss changes from view to view, but not evenly. The column sums of the -# projector (`colsums`, the summed weights of all rays in every pixel) show it: with even coverage -# they would be smooth, here they vary by about 10 % in a fine pattern. +# From view to view the rays shift and miss other pixels, but the misses do not cancel. The column +# sums of the system matrix show it: `colsums` is the back projection of data that are all ones, in +# every pixel the summed weights of all rays, and with even coverage it would be flat. op = SiddonProjector(image, geometry) x_center = reconstruct(op); # ## Several rays per cell # -# `MultiRay` wraps the geometry and traces several rays per cell, here 4 at the centers of 4 -# subcells along `u` (right in the figure above). The projection of a cell is the mean of its rays' -# line integrals; the data keep their size. +# `MultiRay` wraps the geometry and traces several rays per cell, here at the centers of 4 equal +# parts of the cell along `u` (right in the figure above). The value of a cell is the mean of the +# line integrals of its rays; the data keep their size. op_multi = SiddonProjector(image, MultiRay(geometry, RegularPoints(4, 1))) x_multi = reconstruct(op_multi); @@ -67,16 +71,20 @@ x_multi = reconstruct(op_multi); c, c_multi = colsums(op), colsums(op_multi) #hide m = sum(c[49:80, 49:80]) / 32^2 #hide -fig = images_figure([("column sums, 1 ray per cell", c), ("column sums, 4 rays per cell", c_multi)]; #hide +label(name, img) = "$name\nvariation $(round(variation(img, 60); digits=2)) %" #hide +fig = images_figure([(label("1 ray per cell", c), c), (label("4 rays per cell", c_multi), c_multi)]; #hide colorrange=(0.8m, 1.2m)) #hide save("multiray_colsums.png", fig) #hide nothing #hide # ![Column sums with one and with four rays per cell](multiray_colsums.png) -# Both images are shown in the window of ±20 % around the mean. With 4 rays per cell the column sums -# vary by less than 1 % inside the field of view. The -# reconstructions, in the color window 0 to 0.4 with the relative error -# ``\|x - x_\text{true}\| / \|x_\text{true}\|``: +# The titles give the variation, the standard deviation relative to the mean inside a circle of radius +# 60 pixels; the window is ±20 % around the mean. With one ray per cell the column sums show rings and +# a cross around the center, with peaks of about ±30 %. With 4 rays per cell the variation drops by a +# factor of about 30, and the peaks to about ±1 %. +# +# The reconstructions, in the color window 0 to 0.4, which shows the soft tissue and clips the skull, +# with the relative error ``\|x - x_\text{true}\| / \|x_\text{true}\|``: fig = images_figure([("phantom", x_true), ("1 ray per cell", x_center), ("4 rays per cell", x_multi)]; #hide colorrange=(0, 0.4), truth=x_true) #hide @@ -84,48 +92,63 @@ save("multiray_reconstructions.png", fig) #hide nothing #hide # ![Reconstructions with one and with four rays per cell](multiray_reconstructions.png) -# The model with several rays per cell is closer to what the cells measure, and the error drops from -# about 0.33 to 0.27. `JosephProjector` takes a `MultiRay` the same way. In 2D the rays stay in the -# image plane, so the points are spread along `u` only (`RegularPoints(4, 1)`); for a 3D image, -# `RegularPoints(4)` places 4 × 4 rays per cell. +# The model with several rays per cell is closer to what the cells measure: the grainy pattern inside +# the head and between the ellipses goes away, and the error drops from 0.33 to 0.27. The streaks +# outside the head remain, and more rays do not lower the error further (16 rays give 0.27): what +# remains comes from the width of the cells and the number of iterations. +# +# `JosephProjector` takes a `MultiRay` the same way. The rays of a 2D image stay in the image plane, +# so the points are spread along `u` only (`RegularPoints(4, 1)`; `RegularPoints(4)` would be an +# error). For a 3D image, `RegularPoints(4)` places 4 × 4 = 16 rays per cell, and the time grows +# accordingly. # -# ## Other placements and weights +# ## Other point strategies # -# Besides the regular points, a cell can be sampled at random positions, fixed by a seed: +# Besides `RegularPoints`, the rays of a cell can end at random offsets, fixed by a seed (default 0): # -# * `JitteredPoints(4, 1; seed=1)`: each of the regular points moved to a random position within its -# subcell; +# * `JitteredPoints(4, 1; seed=1)`: each regular point moved to a random position within its part of +# the cell; # * `RandomPoints(4; seed=1)`: 4 random points anywhere in the cell; -# * `GaussianPoints(4, 0.3; seed=1)`: 4 random points, normally distributed around the cell center, -# which can also land in neighboring cells (blur); -# * `Points(offsets)`: given positions, the same in every cell; +# * `GaussianPoints(4, 0.3; seed=1)`: 4 random points, normally distributed around the cell center +# with a standard deviation of 0.3 cells, which can also land in neighboring cells (blur). # -# and the rays can be weighted, e.g. by a detector response that falls off towards the cell edges: +# `Points(offsets)` takes given offsets, a `2 × K` matrix in cell units, the same in every cell. How +# the strategies compare on the data above: -g = MultiRay(geometry, JitteredPoints(4, 1; seed=1); weights=GaussianWeights(0.4)) -op_jittered = SiddonProjector(image, g) -x_jittered = reconstruct(op_jittered); +op_jittered = SiddonProjector(image, MultiRay(geometry, JitteredPoints(4, 1; seed=1))) +op_random = SiddonProjector(image, MultiRay(geometry, RandomPoints(4; seed=1))) +x_jittered, x_random = reconstruct(op_jittered), reconstruct(op_random); #- -fig = images_figure([("4 regular rays", x_multi), ("4 jittered rays,\nGaussian weights", x_jittered)]; #hide +fig = images_figure([("4 regular rays", x_multi), ("4 jittered rays", x_jittered), ("4 random rays", x_random)]; #hide colorrange=(0, 0.4), truth=x_true) #hide save("multiray_jittered.png", fig) #hide nothing #hide -# ![Regular and jittered rays](multiray_jittered.png) +# ![Regular, jittered and random rays](multiray_jittered.png) + +# Regular points fit best here: for data that average over the cell, they sample it at equal +# distances. Random offsets add noise to the model, visible as grain: 4 jittered rays lose most of the +# gain, and 4 random rays are worse than one ray per cell. Random strategies need many more rays per +# cell, so use `RegularPoints` unless the detector calls for something else. +# +# ## Weights +# +# The rays of a cell can also count differently, e.g. for a detector response that falls off towards +# the cell edges. The weights are normalized per cell: + +op_weighted = SiddonProjector(image, MultiRay(geometry, RegularPoints(4, 1); weights=GaussianWeights(0.4))); -# With data that are exact cell averages, as here, the regular points fit best: they are the -# midpoint rule for the average over the cell. The jittered points with Gaussian weights describe a -# different detector, one whose response falls off towards the cell edges, and the error is larger. -# Random positions break up the regular pattern of the rays; which placement and weights fit a -# scanner depends on its detector. +# On the data above this changes little (error 0.27), since these data give all parts of a cell the +# same weight. # # ## Cost and devices # -# * The time grows with the number of rays per cell: 4 rays take about 4 times as long; the -# placement and the weights cost nothing measurable (see the -# [benchmarks](https://github.com/JuliaImageRecon/CTOperators.jl/blob/main/benchmark/README.md)). +# * The time grows with the number of rays per cell: 4 rays take about 4 times as long. Jittered +# points and Gaussian weights added nothing measurable on a GPU (see the +# [benchmarks](https://github.com/JuliaImageRecon/CTOperators.jl/blob/main/benchmark/README.md#Several-rays-per-cell-(GPU))). # No extra memory is needed. # * `MultiRay` works on every device (`storage=CuArray` etc.), and the random positions are the same # on all of them, in the forward and the back projection, and for `subset(op, views)`. -# * Masks work as before: a masked data element skips all its rays. +# * Masks work as on plain geometries: a masked data element skips all its rays. +# * Only detector geometries have cells: a `RayList` cannot be wrapped. diff --git a/docs/literate/howto/write_projector.jl b/docs/literate/howto/write_projector.jl index 4201c10..956ae64 100644 --- a/docs/literate/howto/write_projector.jl +++ b/docs/literate/howto/write_projector.jl @@ -3,7 +3,9 @@ # This guide shows what a new projector has to provide and what it gets in return, by implementing a # simple one: it samples the image at equally spaced points along each ray, taking the value of the # nearest voxel. The method is too crude for practice, but its kernels are short, so the parts every -# projector needs stand out. +# projector needs stand out. A projector does not have to live in CTOperators.jl: it fits in through +# multiple dispatch from a script or another package. [Write a Sampling Strategy](@ref) continues +# with this projector. using CTOperators, KernelAbstractions, Atomix, LinearAlgebra, ImagePhantoms using ImageGeoms: axis @@ -45,8 +47,9 @@ nothing #hide # LinearOperatorCollection.jl) and automatic differentiation (the ChainRulesCore rule); a projector # also with `subset` for ordered subsets and `adapt` to move it to another device. # -# The sampling projector is an `AbstractProjector`. It works with every geometry, since it only asks -# each ray for its end points. +# The sampling projector is an `AbstractProjector`. It works with every geometry that has one ray per +# data element, since it only asks each ray for its end points; [Write a Sampling Strategy](@ref) +# extends it to several rays per cell. # ## The struct # @@ -179,7 +182,7 @@ size(op), datasize(op) x, y = rand(Float32, imagesize(op)), rand(Float32, datasize(op)) dot(op * x, y) / dot(x, adjoint(op) * y) -# Subsets, normalizations and moving to another device come for free: +# Subsets, normalizations and moving to another device follow from the interface: datasize(subset(op, 1:10)), size(colsums(op)) diff --git a/docs/literate/howto/write_sampling.jl b/docs/literate/howto/write_sampling.jl index 72e6fda..989e93d 100644 --- a/docs/literate/howto/write_sampling.jl +++ b/docs/literate/howto/write_sampling.jl @@ -3,8 +3,9 @@ # [`MultiRay`](@ref) places several rays in a detector cell according to a point strategy and weights # them with a weighting strategy. Both are ordinary types, and new ones do not have to be part of # CTOperators.jl: a script or a package of your own can define them, and they fit in through multiple -# dispatch. This guide writes one of each, and then lets the sampling projector of -# [Write a Projector](@ref) use them. +# dispatch. This guide writes one of each, and then lets `SamplingProjector` from +# [Write a Projector](@ref) use them. It continues that page: to run this one, run the code there +# first; `image`, `geometry`, `SamplingProjector` and `nearest_voxel` come from there. using CTOperators, KernelAbstractions, LinearAlgebra using CTOperators: Gather, Scatter, VoxelGrid, clip, to_index, trace_element, visit @@ -14,9 +15,11 @@ include(joinpath(pkgdir(CTOperators), "docs", "utils", "multiray_figures.jl")) # # ## A point strategy # # A point strategy is a subtype of [`AbstractPoints`](@ref)`{K}`, with `K` the number of rays per -# cell. It implements [`CTOperators.point_offset`](@ref), which returns the offset of ray `k` from -# the cell center in cell units, where ±1/2 are the cell edges. As an example, `K` points on a circle -# around the cell center, rotated by a random angle in every cell: +# cell. It implements [`CTOperators.point_offset`](@ref). This function gets the element type `T` of +# the projector, the cell `(iu, iv)`, the `view` (numbered as in the full scan) and the ray `k` from 1 +# to `K`, and returns the offset `(du, dv)` from the cell center in cell units, where ±1/2 are the +# cell edges. As an example, `K` points on a circle around the cell center, rotated by a random angle +# in every cell: struct RingPoints{K,T} <: AbstractPoints{K} radius::T @@ -32,33 +35,37 @@ RingPoints(K, radius; seed=0) = RingPoints{K,float(typeof(radius))}(radius, UInt end # [`CTOperators.random_pair`](@ref) gives random numbers that are fixed by the seed, the cell, the -# view and a number of your choice, here always 1, so that all rays of a cell share the same rotation. +# view and a number, usually the ray; here always 1, so that all rays of a cell share one rotation. # They are the same in the forward and the back projection and on every device; [Several Rays per # Cell](@ref) explains why a generator with state does not work here. # -# `point_offset` runs inside the kernels, so it must not allocate. Since `RingPoints` stores a -# floating-point number, it also says how to convert it, so that it follows the element type of the -# projector (some GPUs have no Float64): +# `point_offset` runs inside the kernels, so it must not allocate, and the fields of the type must be +# plain values, not arrays. Since `RingPoints` stores a floating-point number, it also says how to +# convert it, so that it follows the element type of the projector (some GPUs have no Float64): CTOperators.points_eltype(p::RingPoints{K}, ::Type{T}) where {K,T} = RingPoints{K,T}(p.radius, p.seed) # ## A weighting strategy # # A weighting strategy is a subtype of [`AbstractWeights`](@ref) and implements -# [`CTOperators.subray_weight`](@ref), the positive weight of a ray at an offset. The weights of a -# cell are normalized in the kernel. A response that falls off linearly from the cell center: +# [`CTOperators.subray_weight`](@ref), the weight of a ray at an offset. The kernel divides by the sum +# of the weights of a cell, which must not be zero. A response that falls off linearly from the cell +# center: struct LinearWeights <: AbstractWeights end -CTOperators.subray_weight(::LinearWeights, du, dv) = (1 - abs(du)) * (1 - abs(dv)) +CTOperators.subray_weight(::LinearWeights, du, dv) = max(1 - abs(du), 0) * max(1 - abs(dv), 0) # ## Using them # -# The new types work like those of the package, with every projector that supports `MultiRay`. The -# points in one cell, with the marker area showing the weights: - -fig = points_figure([("RingPoints(8, 0.3)", RingPoints(8, 0.3; seed=1), UniformWeights()), #hide - ("RingPoints(8, 0.3; seed=2),\nLinearWeights()", RingPoints(8, 0.3; seed=2), LinearWeights())]) #hide +# The new types work like those of the package, with every projector that supports `MultiRay`. In +# one cell, the blue square, with the marker area relative to the largest weight in the panel: the +# ring with two seeds, which rotate it, as other cells do; and `LinearWeights` on regular points, +# where the corners count less. + +fig = points_figure([("RingPoints(8, 0.3; seed=1)", RingPoints(8, 0.3; seed=1), UniformWeights()), #hide + ("RingPoints(8, 0.3; seed=2)", RingPoints(8, 0.3; seed=2), UniformWeights()), #hide + ("RegularPoints(3, 3),\nLinearWeights()", RegularPoints(3, 3), LinearWeights())]) #hide save("ring_points.png", fig) #hide nothing #hide # ![Ring points in one cell](ring_points.png) @@ -70,16 +77,16 @@ dot(op * x, y) / dot(x, adjoint(op) * y) # ## A projector with several rays per cell # -# The sampling projector of [Write a Projector](@ref) asks every ray for its end points with `ray`, so -# it has one ray per data element. A `MultiRay` has no such ray. To work on every geometry, a projector -# moves the code for one ray into a function and traces through -# [`CTOperators.trace_element`](@ref), which calls it once for a plain geometry and once per ray of -# the cell for a `MultiRay`. +# `SamplingProjector` asks every ray for its end points with `ray`, so it has one ray per data +# element. A `MultiRay` has no such ray. To work on every geometry, a projector moves the code for one +# ray into a function and traces through [`CTOperators.trace_element`](@ref), which calls it once for +# a plain geometry and once per ray of the cell for a `MultiRay`. # -# The function takes the two end points. [`CTOperators.visit`](@ref) does what the two kernels did -# so far: with a [`CTOperators.Gather`](@ref) it sums the voxels, with a -# [`CTOperators.Scatter`](@ref) it adds to them. The option `samples` travels in a small struct, -# since a closure in a kernel could not be compiled for GPUs: +# The function takes an action, the grid, the two end points and the ray type `rt` (see +# [`raytype`](@ref)). [`CTOperators.visit`](@ref) does what the two kernels of +# [Write a Projector](@ref) do: with a [`CTOperators.Gather`](@ref) it sums the voxels, with a +# [`CTOperators.Scatter`](@ref) it adds to them. The option `samples` travels in a small callable +# struct with one concrete field, which `trace_element` passes on: struct SamplingSegment{T} samples::T @@ -122,15 +129,21 @@ function CTOperators.backward!(x, op::SamplingProjector, y) x end -# The sampling projector now takes every strategy, also the ones above, and its adjoint stays exact: +# `SamplingProjector` takes every strategy, also the ones above, and its adjoint stays exact: op = SamplingProjector(image, MultiRay(geometry, RingPoints(8, 0.3; seed=1); weights=LinearWeights())) dot(op * x, y) / dot(x, adjoint(op) * y) -# On plain geometries it computes what it did before. For 2D images, which have one layer of voxels -# along `z`, a projector's constructor should pass the geometry through -# [`CTOperators.in_plane`](@ref), as the projectors of the package do; it keeps the rays in the image -# plane. +# On plain geometries it computes the same as the kernels of [Write a Projector](@ref), up to rounding. +# For images with one layer of voxels along `z`, the constructor of a projector keeps the rays in the +# image plane, as those of the package do: +# +# ```julia +# g = device_geometry(geometry, T, storage) +# g = CTOperators.one_layer(image) ? CTOperators.in_plane(g) : g +# ``` +# +# A custom strategy gets offsets along `v` of 0 there, whatever it returns. # # ## What fits where # diff --git a/docs/make.jl b/docs/make.jl index eccbde5..e84e7fe 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -16,6 +16,8 @@ DocMeta.setdocmeta!(CTOperators, :DocTestSetup, :(using CTOperators); recursive= makedocs(; modules=[CTOperators], sitename="CTOperators.jl", + pagesonly=true, + checkdocs=:exports, authors="Niklas Hackelberg and contributors", repo=Remotes.GitHub("JuliaImageRecon", "CTOperators.jl"), format=Documenter.HTML(; @@ -65,6 +67,7 @@ makedocs(; "Siddon's Method" => "generated/explanation/siddon.md", "Joseph's Method" => "generated/explanation/joseph.md", "Several Rays per Cell" => "generated/explanation/multiray.md", + "Comparing Projectors" => "generated/explanation/comparing.md", ], "Operators" => [ "Masks" => "generated/explanation/masks.md", diff --git a/docs/src/api/geometries.md b/docs/src/api/geometries.md index 2b0865c..377ea69 100644 --- a/docs/src/api/geometries.md +++ b/docs/src/api/geometries.md @@ -23,7 +23,8 @@ Private = false ## Several rays per cell A detector geometry wrapped with a point strategy (where in each cell the rays end) and a -weighting strategy. +weighting strategy, together with the functions that custom strategies implement (see +[Write a Sampling Strategy](@ref)). ```@autodocs Modules = [CTOperators] diff --git a/docs/src/api/internals.md b/docs/src/api/internals.md index 130c78f..657bae4 100644 --- a/docs/src/api/internals.md +++ b/docs/src/api/internals.md @@ -1,8 +1,10 @@ # Internals -For developers of operators (see [Write a Projector](@ref)): the functions each operator -implements, the masks a projector stores, the kernel form of the image grid and the clipping of -rays. None of these is exported; those that custom components use are declared `public`. +For developers of operators and strategies (see [Write a Projector](@ref) and +[Write a Sampling Strategy](@ref)): the kernel form of the image grid, the clipping and tracing of +rays, the masks a projector stores, and private functions. The public functions they implement are +on [Operators](operators.md) and [Projection Geometries](geometries.md). None of these is exported; +those that custom components use are declared `public`. ## Operator implementation @@ -17,6 +19,7 @@ Public = false ```@docs CTOperators.ProjectorMasks CTOperators.projector_masks +CTOperators.is_masked ``` ## Geometries diff --git a/docs/src/api/projectors.md b/docs/src/api/projectors.md index d1e12d7..68c0b90 100644 --- a/docs/src/api/projectors.md +++ b/docs/src/api/projectors.md @@ -1,7 +1,7 @@ # Projectors -The projectors of the package. Both trace one ray per detector cell through the voxel grid and -have an exact adjoint; [Siddon's Method](@ref) and [Joseph's Method](@ref) explain how they +The projectors of the package. Both trace the rays of a geometry through the voxel grid, one per +detector cell or several with [`MultiRay`](@ref), and have an exact adjoint; [Siddon's Method](@ref) and [Joseph's Method](@ref) explain how they weight the voxels. ## Siddon diff --git a/docs/src/index.md b/docs/src/index.md index 4e658fe..7b64027 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -59,7 +59,8 @@ y = op * rand(Float32, 128, 128, 128) # projections, size (256, 128, 360) - 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), - [several rays per cell](@ref "Several Rays per Cell"), and + [several rays per cell](@ref "Several Rays per Cell"), how they + [compare](@ref "Comparing Projectors"), and [workgroups and ray order on GPUs](@ref "Workgroups and Ray Order on GPUs"). - The reference lists all types and functions by domain: [Images](api/images.md), [Projection Geometries](api/geometries.md), [Operators](api/operators.md), [Projectors](api/projectors.md) and, for implementers, [Internals](api/internals.md). diff --git a/docs/utils/comparison_figures.jl b/docs/utils/comparison_figures.jl new file mode 100644 index 0000000..2c64ce5 --- /dev/null +++ b/docs/utils/comparison_figures.jl @@ -0,0 +1,59 @@ +# Figures of the explanation page "Comparing Projectors" (CairoMakie): the same checks for several +# projectors and ratios of detector cell size to pixel size. Included invisibly by the page. + +using CairoMakie + +# a 2D image of unit pixels and a parallel-beam scan with cells of width du (in pixels) that covers it +comparison_image(n) = ImageGeom(; dims=(n, n), deltas=(1.0, 1.0)) +comparison_scan(n, du, angles) = parallel_beam(; angles, nu=ceil(Int, 1.5n / du), nv=1, du) + +# the variation of an image inside the circle of radius r around its center, in percent +function variation(img, r) + n = size(img, 1) + inside = [img[i, j] for i in 1:n, j in 1:n if hypot(i - (n + 1) / 2, j - (n + 1) / 2) <= r] + m = sum(inside) / length(inside) + 100 * sqrt(sum((inside .- m) .^ 2) / length(inside)) / m +end + +# `models` holds (name, projector, geometry -> geometry) tuples, `ratios` the cell sizes in pixels. +# One panel per model and ratio; `panel!(ax, op, du)` draws it. +function comparison_grid(panel!, models, ratios, n, angles, window; size=(250 * length(ratios) + 140, 250 * length(models))) + fig = Figure(; size) + Colorbar(fig[1:length(models), length(ratios)+1]; colormap=:RdBu, limits=window, label="coverage") + for (i, (name, P, wrap)) in enumerate(models), (j, du) in enumerate(ratios) + ax = Axis(fig[i, j]; title=i == 1 ? "cells of $du pixels" : "", ylabel=j == 1 ? name : "") + op = P(comparison_image(n), wrap(comparison_scan(n, du, angles)); eltype=Float64) + panel!(ax, op, du) + end + fig +end + +function image_panel!(ax, img, window; crop=nothing) + ax.aspect = DataAspect() + ax.xticksvisible = ax.yticksvisible = false + ax.xticklabelsvisible = ax.yticklabelsvisible = false + shown = isnothing(crop) ? img : img[crop, crop] + heatmap!(ax, shown; colormap=:RdBu, colorrange=window) +end + +# the back projection of one view of ones, scaled so that even coverage is 1 +function view_figure(models, ratios; n=64, angle=0.3, crop=25:40) + window = (0, 2) + comparison_grid(models, ratios, n, [angle], window) do ax, op, du + image_panel!(ax, du .* colsums(op), window; crop) + end +end + +# the back projection of all views of ones (the "pixel support"), scaled so that even coverage is 1, +# with its variation below the panel +function support_figure(models, ratios; n=64, angles=range(0, pi, 181)[1:180]) + window = (0.8, 1.2) + comparison_grid(models, ratios, n, angles, window) do ax, op, du + c = du .* colsums(op) ./ length(angles) + image_panel!(ax, c, window) + circle = 24 .* (cos.(range(0, 2pi, 100)), sin.(range(0, 2pi, 100))) + lines!(ax, (n + 1) / 2 .+ circle[1], (n + 1) / 2 .+ circle[2]; color=:black, linewidth=0.7) + v = variation(c, 24) + ax.xlabel = v < 0.05 ? "variation < 0.05 %" : "variation $(round(v; digits=1)) %" + end +end diff --git a/docs/utils/multiray_figures.jl b/docs/utils/multiray_figures.jl index 053f95f..058e322 100644 --- a/docs/utils/multiray_figures.jl +++ b/docs/utils/multiray_figures.jl @@ -20,8 +20,9 @@ function rays_figure(; n=8, du=2.0, k=4, size=(820, 470)) lines!(ax, [w, w], [0, n]; color=:gray75, linewidth=0.8) lines!(ax, [0, n], [w, w]; color=:gray75, linewidth=0.8) end - # the detector: cells c = 1, ..., ncells with centers at x = (c - (ncells + 1)/2) du + n/2 - x(c) = (c - (ncells + 1) / 2) * du + n / 2 + # the detector: cells c = 1, ..., ncells; the rays through their centers run in the middle of + # pixels, so that they do not hide the lines of the grid + x(c) = (c - (ncells + 1) / 2) * du + n / 2 + 0.5 for c in 1:ncells poly!(ax, Rect2f(x(c - 0.5), -1.2, du, 0.5); color=isodd(c) ? (:steelblue, 0.25) : (:steelblue, 0.45), strokecolor=:steelblue, strokewidth=0.8) From 3bc13dea6bc17ed5a2dc0e65ee0bfff16da7d577 Mon Sep 17 00:00:00 2001 From: Niklas Hackelberg Date: Wed, 4 Nov 2026 21:55:02 +0100 Subject: [PATCH 11/11] Docs: state what was measured A pass over all pages against one rule: a cause, a general rule or a recommendation needs a paper, otherwise the text says what was measured on what setup. Wrong statements corrected, among them the lines a ray mask leaves out, the short fan-beam scan and the claims about Siddon's continuity; unsupported causes dropped; numbers given with their setup; papers cited for what they say. Co-Authored-By: Claude Opus 5.5 --- benchmark/README.md | 18 +++++-------- docs/literate/explanation/comparing.jl | 33 ++++++++++++------------ docs/literate/explanation/joseph.jl | 13 +++++----- docs/literate/explanation/launch.jl | 8 +++--- docs/literate/explanation/masks.jl | 10 +++---- docs/literate/explanation/multiray.jl | 32 ++++++++++------------- docs/literate/explanation/siddon.jl | 15 ++++++----- docs/literate/howto/gpu.jl | 4 +-- docs/literate/howto/masks.jl | 17 ++++++------ docs/literate/howto/multiray.jl | 30 +++++++++++---------- docs/literate/howto/sinograms.jl | 6 ++--- docs/literate/tutorials/cone_beam.jl | 6 ++--- docs/literate/tutorials/fan_beam.jl | 6 ++--- docs/literate/tutorials/helical.jl | 7 +++-- docs/literate/tutorials/parallel_beam.jl | 7 ++--- docs/literate/tutorials/sparse_views.jl | 10 +++---- src/geometry/multiray.jl | 5 ++-- src/operators/joseph/joseph.jl | 4 +-- 18 files changed, 113 insertions(+), 118 deletions(-) diff --git a/benchmark/README.md b/benchmark/README.md index 9eab0fe..e81ecbb 100644 --- a/benchmark/README.md +++ b/benchmark/README.md @@ -84,9 +84,9 @@ cone128: 128³ volume, cone beam with a 128² detector and 180 views (2.95 M ray volume, parallel beam with a 128² detector and 180 views. Phantom error: relative error of the projected voxelized phantom against its analytic line integrals. -CTOperators' Joseph gives the same phantom error as parallelproj's (the same method) and runs at the -same speed; Siddon is about twice as fast as Joseph on the CPU (one memory read per visited voxel -against four per sample). +CTOperators' Joseph gives the same phantom error as parallelproj's (the same method). It runs at the +same speed for cone128 and 1.2 to 1.4× faster for parallel128. On the CPU, Siddon is 1.4 to 1.7× +faster than Joseph forward and 2.0 to 2.5× back. ### GPU @@ -122,21 +122,17 @@ Observations: for maximal occupancy) make the forward projection 1.2–1.9× faster even in data order; tiles of 8 × 8 detector cells per workgroup instead of consecutive cells of a row give Joseph another 1.4–1.8×, Siddon little. Together at 256³: Joseph about 2×, Siddon 1.2–1.6× faster, at 128³ up - to 1.25×. The back projection, dominated by its atomic additions, gains nothing (1–6 % slower, - close to the noise); it runs with the same configuration to keep one code path. On the CPU tiles + to 1.25×. The back projection gains nothing (1–6 % slower, close to the noise); it runs with the same configuration to keep one code path. On the CPU tiles change the forward projection by less than 4 % (8², 16², 32², 24 threads), so the CPU keeps data order and the backend's workgroups. A count of the cache lines per load instruction of a - wavefront (32 lanes) gives about 17 in both orders: the gain of tiles is not better coalescing - within one load; finding the cause needs hardware counters. + wavefront (32 lanes) gives about 17 in both orders. The cause of the gain was not measured. - Ray setup is negligible: computing and clipping all rays of cone256 takes 0.8 ms, and a precomputed ray list (540 MB) gives the same times. - From 128³ to 256³ the work grows 16× (8× the rays, each twice as long), our GPU times 17–21× - forward (26–32× before) and 17–22× back. A 128³ Float32 image (8 MB) fits into the - GPU's 80 MB last-level cache, a 256³ image (64 MB) hardly. + forward (26–32× before) and 17–22× back. - RadonKA, slice by slice, is faster than our parallel-beam forward projection at 256³ (91 ms 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. +- Joseph's back projection takes 4–5× its forward projection, Siddon's about 2.3×. ### Masks (GPU) diff --git a/docs/literate/explanation/comparing.jl b/docs/literate/explanation/comparing.jl index c4068de..11fbd02 100644 --- a/docs/literate/explanation/comparing.jl +++ b/docs/literate/explanation/comparing.jl @@ -1,7 +1,7 @@ # # Comparing Projectors # -# How evenly a projector covers the image depends mainly on the ratio of detector cell size to voxel -# size and on the number of rays per cell, less on the projector. This page compares +# How evenly a ray-driven projector covers the image depends on the ratio of detector cell size to +# voxel size (Miao 2013). This page compares # [Siddon's Method](@ref) and [Joseph's Method](@ref) with one and with several rays per cell (see # [Use Several Rays per Cell](@ref)) for cells smaller than, equal to and larger than the pixels of a # 2D image, with two checks from Miao (2013): the back projection of one view and of all views of @@ -51,10 +51,10 @@ save("comparison_support.png", fig) #hide nothing #hide # ![Back projection of all views](comparison_support.png) -# The pattern of the single views does not average out. With cells of 2 pixels, the coverage varies -# by about 10 %, in a cross and a checker pattern around the center; four rays per cell reduce that -# to a fraction of a percent. An iterative reconstruction sees this pattern as part of the model, and -# it shows up as artifacts in the images, see [Use Several Rays per Cell](@ref). +# With cells of 2 pixels, the coverage varies +# by 11 % for Siddon's method and 8.5 % for Joseph's, in a cross and a checker pattern around the center; four rays per cell reduce that +# to a fraction of a percent. [Use Several Rays per Cell](@ref) shows a reconstruction with this kind +# of pattern. # # ## How many rays # @@ -64,20 +64,19 @@ coverage(P, k) = 2 .* colsums(P(comparison_image(64), MultiRay(comparison_scan(6 [(rays=k, Siddon=round(variation(coverage(SiddonProjector, k), 24); digits=2), Joseph=round(variation(coverage(JosephProjector, k), 24); digits=2)) for k in (1, 2, 3, 4, 8)] -# More rays help, but not steadily: the number matters relative to the pixels, and an odd number of -# rays can be worse than a smaller even one. A spacing of the rays at the object of at most half a -# voxel covers the image well. +# On this example the variation does not fall steadily with the number of rays: for Siddon's method, 3 +# rays give more than 2. # -# ## What to take from it +# ## Summary of the example # -# * **The ratio matters most.** With cells no larger than the voxels, both methods cover the image -# well with one ray per cell, and Joseph's method is the more even of the two. -# * **For cells larger than the voxels,** one ray per cell is not enough for either method; several -# rays per cell are, at the cost of about that many times the computing time (see the +# * With cells no larger than the pixels, the variation with one ray per cell is about 1 % or less +# for both methods. +# * With cells of 2 pixels, it is 11 % (Siddon) and 8.5 % (Joseph) with one ray per cell and below +# 0.25 % with four, at +# about four times the computing time (see the # [benchmarks](https://github.com/JuliaImageRecon/CTOperators.jl/blob/main/benchmark/README.md)). -# Voxels no smaller than the cells avoid the problem as well, at the price of a coarser image. -# * **In a cone beam,** the size to compare with the voxels is the cell size at the object: the cell -# size divided by the magnification. +# * In a cone beam, the size to compare with the voxels is the cell size at the object: the cell size +# divided by the magnification. # # ## References # diff --git a/docs/literate/explanation/joseph.jl b/docs/literate/explanation/joseph.jl index 633a599..0849ba7 100644 --- a/docs/literate/explanation/joseph.jl +++ b/docs/literate/explanation/joseph.jl @@ -60,9 +60,9 @@ round.(weights; digits=2) # ## Properties # -# - **Smooth.** The projection changes continuously when a ray moves; Siddon's jumps when the ray -# crosses into the next voxel. Rays that lie exactly between two voxels, ambiguous for Siddon, -# take the mean of both. +# - **Smooth.** The projection changes continuously when a ray moves. Siddon's does too, except for +# rays parallel to the voxel planes, where it jumps as the ray crosses a plane. Rays that lie exactly +# between two voxels, ambiguous for Siddon, take the mean of both. # - **Exact for linear images.** If the voxel values are samples of a linear function of the # position, the interpolation reproduces it, and the sum per layer integrates it exactly, for rays # that cross the volume from face to face and stay within the outermost voxel centers across the @@ -92,9 +92,8 @@ save("coverage_comparison.png", fig) #hide nothing #hide # ![Coverage of the pixels by the rays of one view, Siddon and Joseph](coverage_comparison.png) -# For reconstructions this means that Joseph's method is less prone to checkerboard patterns, but at -# angles near 45° it needs as many rays per voxel as Siddon's. Projectors that integrate over the -# detector cells (footprint methods) avoid the problem altogether. +# De Man & Basu (2004) report that a distance-driven projector, which integrates over the detector +# cells, removes the artifacts of ray-driven back projection in their examples. # # ## Implementation # @@ -113,3 +112,5 @@ nothing #hide # # - P. M. Joseph, "An improved algorithm for reprojecting rays through pixel images", IEEE Trans. # Med. Imaging 1(3), 192–196, 1982. +# - B. De Man, S. Basu, "Distance-driven projection and backprojection in three dimensions", Phys. +# Med. Biol. 49(11), 2463–2475, 2004. diff --git a/docs/literate/explanation/launch.jl b/docs/literate/explanation/launch.jl index dd1a8a4..6cd3988 100644 --- a/docs/literate/explanation/launch.jl +++ b/docs/literate/explanation/launch.jl @@ -46,9 +46,8 @@ nothing #hide # [benchmarks](https://github.com/JuliaImageRecon/CTOperators.jl/blob/main/benchmark/README.md) # on an AMD Radeon RX 7900 XT with 256³ images, the forward projection became faster with workgroups # of 64 threads instead of the size KernelAbstractions.jl chooses, for both methods. Tile order made -# it faster again for Joseph's method, but changed little for Siddon's. The back projection, whose -# time goes mostly into the atomic additions, became neither faster nor slower. On the CPU, neither -# choice made a noticeable difference. +# it faster again for Joseph's method, but changed little for Siddon's. The back projection became +# neither faster nor slower. On the CPU, tiles changed the forward projection by less than 4 %. # # ## Configuration # @@ -58,8 +57,7 @@ nothing #hide # - on the CPU: the workgroup size that KernelAbstractions.jl chooses, and data order. # # A projector moved with `adapt` switches to the defaults of its new device. The GPU defaults were -# chosen on the GPU above; on other GPUs they are a starting point. Tuning is worth it when -# projections dominate the run time of your reconstruction. +# chosen on the GPU above; on other GPUs they are a starting point. # # To choose the settings yourself, pass a [`LaunchConfig`](@ref) with the keyword `launch`. It # applies to the forward and the back projection, on any device, and projectors made from the diff --git a/docs/literate/explanation/masks.jl b/docs/literate/explanation/masks.jl index c540257..f6abe4a 100644 --- a/docs/literate/explanation/masks.jl +++ b/docs/literate/explanation/masks.jl @@ -21,8 +21,7 @@ include(joinpath(pkgdir(CTOperators), "docs", "utils", "mask_figures.jl")) #hide # 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. +# matrix, 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 @@ -35,7 +34,8 @@ include(joinpath(pkgdir(CTOperators), "docs", "utils", "mask_figures.jl")) #hide # # * **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. +# the rays masked, the back projection got 17 to 20 % faster and the forward projection 4 to 14 % +# (256³ volume, AMD Radeon RX 7900 XT). # * **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). @@ -47,8 +47,8 @@ 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 +# For a small region of interest this saves most of the work: a centered mask of 64³ voxels in a 256³ +# volume, on the same GPU, made the projections 30 to 43 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. # diff --git a/docs/literate/explanation/multiray.jl b/docs/literate/explanation/multiray.jl index 6282309..df9e02f 100644 --- a/docs/literate/explanation/multiray.jl +++ b/docs/literate/explanation/multiray.jl @@ -10,9 +10,9 @@ include(joinpath(pkgdir(CTOperators), "docs", "utils", "multiray_figures.jl")) # # ## What is computed # # A ray-driven projector such as [Siddon's Method](@ref) or [Joseph's Method](@ref) treats a -# detector cell as a point: one line from the source to the center of the cell. When the cells are -# large compared to the voxels, the lines miss voxels between them, and the back projection covers -# the voxels unevenly (De Man & Basu 2004). A real cell integrates over its area. `MultiRay(geometry, points; weights)` replaces the one line by `K` rays +# detector cell as a point: one ray from the source to the center of the cell. When the cells are +# large compared to the voxels, the back projection of such rays has gaps (Miao 2013, p. 29). A real +# cell integrates over its area. `MultiRay(geometry, points; weights)` replaces the one ray by `K` rays # per cell, ending at the offsets given by `points`, and averages their line integrals: # # ```math @@ -20,8 +20,8 @@ include(joinpath(pkgdir(CTOperators), "docs", "utils", "multiray_figures.jl")) # # ``` # # The back projection is the exact adjoint: every ray adds its share ``w_{ik} / \sum_l w_{il}`` of -# ``y_i``. The rays of a cell approximate the integral over the cell, in the classification of -# system models by Miao (2013) the "area integral" model, by sampling. +# ``y_i``. Miao (2013) calls the integral over the area of a ray path the "area integral" model; the +# rays of a cell sample such an integral. # # The weights average **line integrals**, which keeps the operator linear. A detector averages # intensities, ``\exp(-\int x)``, which for a cell that sees an edge of the object is not the same @@ -51,15 +51,12 @@ nothing #hide # ![Point strategies in one cell](multiray_points.png) # * **Regular points** sit at the centers of `su × sv` subcells: the midpoint rule for the integral -# over the cell, a sum at equal distances. For data that average over the cell they are the most -# accurate with few points. +# over the cell, a sum at equal distances. # * **Jittered points** move each regular point to a random position within its subcell ("jittered -# sampling", Cook 1986, also called stratified sampling). They keep the even spread, but neighboring -# cells do not share the same pattern, so the rays do not line up with the voxel grid in a regular -# way. The price is noise in the model. -# * **Random points** are spread anywhere in the cell. With few points their mean can differ noticeably -# from the cell average: with 4 rays per cell the example of [Use Several Rays per Cell](@ref) is -# worse than with one. +# sampling", Cook 1986, also called stratified sampling). They keep the even spread, and neighboring +# cells do not share the same pattern. +# * **Random points** are spread anywhere in the cell, so their mean differs from the average over the +# cell by a random amount. [Use Several Rays per Cell](@ref) shows an example with 4 rays per cell. # * **Gaussian points** are normally distributed around the cell center, and all rays count equally: # the distribution of the positions is the weighting. Points may land outside the cell, which # models blur across neighboring cells. @@ -69,7 +66,8 @@ nothing #hide # # For images with one layer of voxels along `z` (2D images, or 3D images with `nz = 1`) the projectors # keep all rays in the image plane: the offsets along `v` are set to 0, and strategies with several -# subcells along `v` are rejected. A ray away from the center of the layer would see too little of it. +# subcells along `v` are rejected. Joseph's method interpolates across the layer, so a ray away from +# its center is weighted down, and a tilted ray may leave the layer. # # ## Random positions without a random number generator # @@ -91,8 +89,8 @@ nothing #hide # ## Implementation # # * The number of rays per cell `K` is a type parameter of the point strategy. The loop over the -# rays of a cell has a length known to the compiler, and with `K = 1` the projector is as fast as -# on the plain geometry. +# rays of a cell has a length known to the compiler. With `K = 1` the projector was as fast as on +# the plain geometry in our benchmark (128³ on a GPU, see below). # * The kernels are those of the plain projectors: each thread handles one data element. For a # `MultiRay` it traces the `K` rays one after the other and adds them up (forward), or first sums # the weights and then lets every ray add its share (back). No memory per ray is needed. @@ -112,5 +110,3 @@ nothing #hide # Techniques 9(3), 21–38, 2020. # - C. Miao, "Comparative studies of different system models for iterative CT image # reconstruction", master's thesis, Wake Forest University, 2013. -# - B. De Man, S. Basu, "Distance-driven projection and backprojection in three dimensions", Phys. -# Med. Biol. 49(11), 2463–2475, 2004. diff --git a/docs/literate/explanation/siddon.jl b/docs/literate/explanation/siddon.jl index 61c5021..ba91398 100644 --- a/docs/literate/explanation/siddon.jl +++ b/docs/literate/explanation/siddon.jl @@ -32,7 +32,7 @@ bonito_page("siddon") #hide # al. (1998) avoid both: starting from the voxel where the ray enters, they take the nearest next # crossing of the three axes, which says which axis the ray crosses next, and step that axis's voxel # index by one. Every step costs a few additions and comparisons, and a ray through a volume of n³ -# voxels takes about 3n steps. +# voxels takes at most about 3n steps. # # In the figure a ray crosses a 2D image. Orange circles mark where the ray crosses the lines `x = # const` between pixels, red diamonds the lines `y = const`; the pixels are colored by the length of @@ -66,18 +66,17 @@ round.(weights; digits=2) # Whether a voxel is hit by many or by few rays then depends on how the rays fall onto the voxel # grid. The figure shows, for one parallel-beam view at 45°, the summed lengths of all rays in every # pixel (scaled by the cell size, so that even coverage is 1). With one ray per pixel the coverage -# varies between about 0.83 and 1.34, in a regular pattern; more rays per pixel even it out. +# varies between about 0.83 and 1.38, in a regular pattern; more rays per pixel even it out here. fig = coverage_figure(16) #hide save("coverage.png", fig) #hide nothing #hide # ![Coverage of the pixels by the rays of one view](coverage.png) -# In an iterative reconstruction this pattern can appear as a checkerboard in the image (De Man & -# Basu 2004 describe the effect). Remedies are more rays than voxels (finer detector sampling, or -# several rays per cell, see [Use Several Rays per Cell](@ref)), -# [Joseph's Method](@ref), whose interpolation evens out the coverage at most angles, or projectors -# that integrate over the detector cells, such as distance-driven or separable-footprint projectors. +# Ray-driven back projection tends to Moiré patterns, and an iterative reconstruction can show such +# artifacts (De Man & Basu 2004). More rays than voxels reduce them (finer detector sampling, or several +# rays per cell, see [Use Several Rays per Cell](@ref)), as do projectors that integrate over the +# detector cells, such as distance-driven ones (De Man & Basu 2004; Miao 2013). # # ## Implementation # @@ -101,3 +100,5 @@ nothing #hide # 6(1), 89–94, 1998. # - B. De Man, S. Basu, "Distance-driven projection and backprojection in three dimensions", Phys. # Med. Biol. 49(11), 2463–2475, 2004. +# - C. Miao, "Comparative studies of different system models for iterative CT image +# reconstruction", master's thesis, Wake Forest University, 2013. diff --git a/docs/literate/howto/gpu.jl b/docs/literate/howto/gpu.jl index d539dc8..fc89016 100644 --- a/docs/literate/howto/gpu.jl +++ b/docs/literate/howto/gpu.jl @@ -73,7 +73,7 @@ # # ## Precision # -# `Float32` is the default and the right choice on GPUs; `Float64` (`eltype=Float64`) on the CPU or +# `Float32` is the default; `Float64` (`eltype=Float64`) on the CPU or # on GPUs with fast double precision gives reference results. The element types of the data must # match the projector's. # @@ -86,6 +86,6 @@ # detector cells) were chosen on an AMD Radeon RX 7900 XT. See # [Workgroups and Ray Order on GPUs](@ref). # - The projector holds no workspace: its memory is the geometry (a few numbers per view for -# detector geometries, two points per ray for ray lists). +# detector geometries, two points per ray for ray lists) and the masks. # - Creating a projector is cheap; for ordered subsets, create the projectors of all subsets once # (see [Ordered Subsets and Normalization](@ref)). diff --git a/docs/literate/howto/masks.jl b/docs/literate/howto/masks.jl index 3e5bddc..5de847e 100644 --- a/docs/literate/howto/masks.jl +++ b/docs/literate/howto/masks.jl @@ -28,16 +28,16 @@ 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): +# mask `MaskCircle()` of ImageGeoms.jl restricts the image to the circle inscribed in it (a cylinder in +# 3D), which here lies inside the circle all views cover: 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: +# The reconstruction is zero outside the circle. The phantom lies inside it, so the data are the same +# as without the mask, but fewer voxels are unknown: after 30 iterations the error is 0.062 with the +# mask and 0.075 without it. fig = images_figure([("image mask", Float32.(image.mask), (0, 1)), ("reconstruction", x)]; #hide colorrange=(0, 0.4)) #hide @@ -79,15 +79,16 @@ nothing #hide # 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. +# ray mask the rings almost disappear (error 0.08; 1.15 without the ray mask, 0.075 with intact data): +# the lines of the defective cells are not measured, but their wrong values no longer enter the fit. # 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 +# center. The rays then only traverse this box: a centered mask of 64³ voxels in a 256³ volume, on an +# AMD Radeon RX 7900 XT, made the projections 30 to 43 times faster in our # [benchmarks](https://github.com/JuliaImageRecon/CTOperators.jl/blob/main/benchmark/README.md). roi = falses(n, n) diff --git a/docs/literate/howto/multiray.jl b/docs/literate/howto/multiray.jl index a35ee6a..4db4ae8 100644 --- a/docs/literate/howto/multiray.jl +++ b/docs/literate/howto/multiray.jl @@ -1,10 +1,9 @@ # # Use Several Rays per Cell # # A detector cell measures the radiation that falls onto its whole area. The projectors trace one -# ray per cell, through its center. When the cells are large compared to the voxels, these rays -# miss voxels between them, and the reconstruction shows patterns that come from the sampling, not -# from the object. This guide shows the problem on a small example and how `MultiRay` traces several -# rays per cell instead. [Several Rays per Cell](@ref) explains how the rays are placed and weighted, +# ray per cell, through its center. When the cells are large compared to the voxels, the back +# projection of such rays has gaps (Miao 2013, p. 29). This guide shows the effect on a small example +# and how `MultiRay` traces several rays per cell instead. [Several Rays per Cell](@ref) explains how the rays are placed and weighted, # and [Comparing Projectors](@ref) when it matters. using CTOperators, ImagePhantoms, RegularizedLeastSquares, LinearOperatorCollection, CairoMakie @@ -92,10 +91,10 @@ save("multiray_reconstructions.png", fig) #hide nothing #hide # ![Reconstructions with one and with four rays per cell](multiray_reconstructions.png) -# The model with several rays per cell is closer to what the cells measure: the grainy pattern inside -# the head and between the ellipses goes away, and the error drops from 0.33 to 0.27. The streaks -# outside the head remain, and more rays do not lower the error further (16 rays give 0.27): what -# remains comes from the width of the cells and the number of iterations. +# With 4 rays per cell the grainy pattern inside the head and between the ellipses goes away, and the +# error drops from 0.33 to 0.27. The streaks outside the head remain. On this example, 16 rays give +# 0.27 as well, and the gap depends on when the iterations stop: with 10 iterations the errors are +# 0.35 and 0.28, with 200 they are 0.33 and 0.32. # # `JosephProjector` takes a `MultiRay` the same way. The rays of a 2D image stay in the image plane, # so the points are spread along `u` only (`RegularPoints(4, 1)`; `RegularPoints(4)` would be an @@ -127,10 +126,9 @@ save("multiray_jittered.png", fig) #hide nothing #hide # ![Regular, jittered and random rays](multiray_jittered.png) -# Regular points fit best here: for data that average over the cell, they sample it at equal -# distances. Random offsets add noise to the model, visible as grain: 4 jittered rays lose most of the -# gain, and 4 random rays are worse than one ray per cell. Random strategies need many more rays per -# cell, so use `RegularPoints` unless the detector calls for something else. +# On these data, whose cells average evenly over their width, regular points give the smallest error +# of the three, 4 jittered rays lie close to one ray per cell, and 4 random rays above it. This is one +# example; other detectors and data can behave differently. # # ## Weights # @@ -139,8 +137,7 @@ nothing #hide op_weighted = SiddonProjector(image, MultiRay(geometry, RegularPoints(4, 1); weights=GaussianWeights(0.4))); -# On the data above this changes little (error 0.27), since these data give all parts of a cell the -# same weight. +# On the data above the error is 0.27, as with uniform weights. # # ## Cost and devices # @@ -152,3 +149,8 @@ op_weighted = SiddonProjector(image, MultiRay(geometry, RegularPoints(4, 1); wei # on all of them, in the forward and the back projection, and for `subset(op, views)`. # * Masks work as on plain geometries: a masked data element skips all its rays. # * Only detector geometries have cells: a `RayList` cannot be wrapped. +# +# ## References +# +# - C. Miao, "Comparative studies of different system models for iterative CT image +# reconstruction", master's thesis, Wake Forest University, 2013. diff --git a/docs/literate/howto/sinograms.jl b/docs/literate/howto/sinograms.jl index 8d1b2fb..bf686f0 100644 --- a/docs/literate/howto/sinograms.jl +++ b/docs/literate/howto/sinograms.jl @@ -80,9 +80,9 @@ solver = createLinearSolver(CGNR, op; AHA=normalOperator(op), iterations=20) x_cg = reshape(solve!(solver, vec(y)), size(x_true)) x_fdk = Float32.(fdk(plan_fbp(rg_helix, image), y)); -# Axial (top) and sagittal (bottom, `z` from bottom to top) slices: CGNR recovers the whole volume -# (with some ringing at the skull, where the voxel model cannot match the finer data); FDK sees only -# a band of slices, smeared by the motion of the source. +# Axial (top) and sagittal (bottom, `z` from bottom to top) slices: CGNR recovers the whole volume, +# with ringing visible in the figure; FDK sees only a band of slices, smeared by the motion of the +# source. fig = Figure(size=(900, 600)) for (j, (name, img)) in enumerate((("phantom", x_true), ("CGNR, helix", x_cg), ("FDK, circle assumed", x_fdk))) diff --git a/docs/literate/tutorials/cone_beam.jl b/docs/literate/tutorials/cone_beam.jl index 29e1abb..6d4fff7 100644 --- a/docs/literate/tutorials/cone_beam.jl +++ b/docs/literate/tutorials/cone_beam.jl @@ -71,7 +71,7 @@ nothing #hide # ![Phantom and reconstruction of the circular scan](cone_circle.png) # The middle disc is reconstructed well. Away from the middle the discs lose their edges and blur -# into each other along `z`, although every voxel was in every view. The reason is what a circular +# into each other along `z`, although every disc was in every view. The reason is what a circular # orbit measures: the data determines the object only in the plane of the orbit. Above and below, # no ray runs parallel to the discs, and an iterative method cannot recover what no measurement # contains (Tuy's condition: every plane that cuts the object has to meet the path of the source; @@ -97,5 +97,5 @@ nothing #hide # ![Large and small cone angle](cone_angles.png) # With nearly parallel rays every slice is measured almost as in a fan-beam scan, and the discs come -# back. Scanners with a large cone instead move the source along `z` as well: -# [Helical Cone-Beam CT](@ref) continues with the same phantom. +# back. A helical scan moves the source along `z` as well: [Helical Cone-Beam CT](@ref) continues with +# the same phantom. diff --git a/docs/literate/tutorials/fan_beam.jl b/docs/literate/tutorials/fan_beam.jl index a19614a..6f5a23a 100644 --- a/docs/literate/tutorials/fan_beam.jl +++ b/docs/literate/tutorials/fan_beam.jl @@ -86,9 +86,9 @@ save("fan_reconstructions.png", fig) #hide nothing #hide # ![Reconstructions of the three scans](fan_reconstructions.png) -# The short scan is as good as the full one: it measured every line, some twice. Half a turn misses -# lines near the edges of the fan, and the reconstruction shows streaks and wrong values on the side -# that was seen least. +# The short scan measured every line, some twice; its error is 0.022, against 0.013 for the full +# scan. Half a turn misses lines near the edges of the fan: its error is 0.139, and the reconstruction +# shows streaks and wrong values on the side that was seen least. # # Filtered back projection (FBP, available in # [Sinograms.jl](https://github.com/JuliaImageRecon/Sinograms.jl); this package provides the diff --git a/docs/literate/tutorials/helical.jl b/docs/literate/tutorials/helical.jl index 278edb3..b991c36 100644 --- a/docs/literate/tutorials/helical.jl +++ b/docs/literate/tutorials/helical.jl @@ -80,10 +80,9 @@ nothing #hide # ![Three pitches](helical_pitch.png) # A voxel near the axis stays in the cone while the source rises by the detector height: for one -# turn at a pitch of 96 mm, for half a turn at 192 mm. Half a turn is less than the half turn plus -# fan angle that [Fan Beam with a Curved Detector](@ref) found necessary, and the reconstruction -# shows shading. Clinical scanners quote the pitch relative to the detector height and use values -# between about 0.5 and 1.5. +# turn at a pitch of 96 mm, for half a turn at 192 mm. At 192 mm the error is 0.064, against 0.026 at +# 48 mm and 0.033 at 96 mm, and the axial slice shows shading. Clinical scanners quote the pitch +# relative to the detector height. # # An object moving along the rotation axis during a circular scan, on a table or a conveyor belt, # sees the same helix: only the relative motion counts. [Getting Started](@ref) reconstructs such a diff --git a/docs/literate/tutorials/parallel_beam.jl b/docs/literate/tutorials/parallel_beam.jl index 71d3ded..83f7359 100644 --- a/docs/literate/tutorials/parallel_beam.jl +++ b/docs/literate/tutorials/parallel_beam.jl @@ -50,9 +50,10 @@ size(y) # Stacked view by view, the detector rows form the sinogram. Every point of the object traces a sine # curve in it: the position on the detector at which the point appears moves back and forth as the -# scanner turns; the farther the point is from the center, the larger the swing. The brightest -# parts, at the outer edges of the sinogram, come from rays that graze the skull of the phantom and -# run a long way inside it. +# scanner turns; the farther the point is from the center, the larger the swing. The bright lines +# at the outer edges of the sinogram come from rays that graze the skull of the phantom and run a +# long way inside it. The largest values, in the middle of the first and last views, come from rays +# along the long axis of the head. fig = sinogram_figure(y; angles) #hide save("parallel_sinogram.png", fig) #hide diff --git a/docs/literate/tutorials/sparse_views.jl b/docs/literate/tutorials/sparse_views.jl index 8e224ca..b82b27a 100644 --- a/docs/literate/tutorials/sparse_views.jl +++ b/docs/literate/tutorials/sparse_views.jl @@ -3,8 +3,8 @@ # Fewer views mean less dose and faster scans, and some scanners cannot turn all the way round. This # tutorial takes views away from a 2D parallel-beam scan like the one of [Parallel Beam in 2D](@ref) # (with a finer detector) in two ways, few views over the full half turn and many views over a -# limited arc, and shows what a regularized reconstruction can recover: much in the first case, -# little in the second. +# limited arc, and shows what a regularized reconstruction recovers on this example: much in the +# first case, little in the second. using CTOperators, ImagePhantoms, RegularizedLeastSquares, LinearOperatorCollection, CairoMakie using ImageGeoms: axis @@ -89,8 +89,8 @@ nothing #hide # With 20 views TV removes most of the streaks and comes close to the phantom, though small details # are blurred: the data and the prior together nearly determine it. Over the limited angle TV -# sharpens the measured edges but cannot restore those the scan never saw. Missing directions are a -# harder loss than missing views. +# sharpens the measured edges but cannot restore those the scan never saw. On this example TV lowers +# the error from 0.40 to 0.055 for 20 views and from 0.40 to 0.33 for the limited arc. # -# The phantom suits TV unusually well: real anatomy has texture and smooth transitions, and the +# The phantom is piecewise constant; this result does not carry over to images with texture. The # weight `λ` trades noise and streaks against loss of detail. diff --git a/src/geometry/multiray.jl b/src/geometry/multiray.jl index 7e1a223..a64d488 100644 --- a/src/geometry/multiray.jl +++ b/src/geometry/multiray.jl @@ -328,8 +328,9 @@ end """ one_layer(image::ImageGeom) -Whether `image` has one layer of voxels along `z`: a 2D image, or a 3D image with `nz = 1`. Rays that -leave the center of the layer see too little of it, so the rays of a cell stay in its plane. +Whether `image` has one layer of voxels along `z`: a 2D image, or a 3D image with `nz = 1`. The rays +of a cell stay in its plane: Joseph's method weights rays away from the center of the layer down, and +a tilted ray may leave the layer. """ one_layer(image::ImageGeom) = length(image.dims) == 2 || image.dims[3] == 1 diff --git a/src/operators/joseph/joseph.jl b/src/operators/joseph/joseph.jl index bd81ec2..f45a45b 100644 --- a/src/operators/joseph/joseph.jl +++ b/src/operators/joseph/joseph.jl @@ -7,8 +7,8 @@ line integral is the sum over the layers of the interpolated value times the len the layer. Compared with [`SiddonProjector`](@ref): * the image model is smooth, so the projection changes continuously when a ray moves; -* the rays of a view cover the voxels more evenly, except near 45°, where the two methods give - similar weights; +* for one parallel view with one ray per voxel, the rays cover the voxels more evenly, except near + 45°, where the two methods give similar weights; * the projection is exact for images that are linear functions of the position rather than for voxel images.