Several rays per detector cell (MultiRay)

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)
 
@@ -157,3 +153,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.

docs/literate/explanation/comparing.jl

@@ -0,0 +1,84 @@
+# # Comparing Projectors
+#
+# 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
+# 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)
+
+# 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
+#
+# 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)]
+
+# On this example the variation does not fall steadily with the number of rays: for Siddon's method, 3
+# rays give more than 2.
+#
+# ## Summary of the example
+#
+# * 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)).
+# * 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.

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.

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

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.
 #
@@ -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.