diff --git a/benchmark/README.md b/benchmark/README.md index fa26950..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) @@ -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. diff --git a/docs/literate/explanation/comparing.jl b/docs/literate/explanation/comparing.jl new file mode 100644 index 0000000..11fbd02 --- /dev/null +++ b/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. 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 2dd326a..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. # @@ -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/explanation/multiray.jl b/docs/literate/explanation/multiray.jl new file mode 100644 index 0000000..df9e02f --- /dev/null +++ b/docs/literate/explanation/multiray.jl @@ -0,0 +1,112 @@ +# # 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 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 +# 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``. 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 +# (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, 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 + ("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 +# ![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. +# * **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, 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. +# * **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 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. 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 +# +# 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](@ref "Ordered Subsets and Normalization") require. A +# different seed gives a different, again fixed, set of positions: `MultiRay(geometry, +# JitteredPoints(4, 1; 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. 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. +# * 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 +# fails with a `MethodError` for `ray` when applied, instead of silently ignoring the other 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. diff --git a/docs/literate/explanation/siddon.jl b/docs/literate/explanation/siddon.jl index 32003b3..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,17 +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), -# [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 # @@ -100,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 new file mode 100644 index 0000000..4db4ae8 --- /dev/null +++ b/docs/literate/howto/multiray.jl @@ -0,0 +1,156 @@ +# # 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, 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 +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 (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 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 iterations of CGNR (conjugate gradients on the normal equations): + +function reconstruct(op) + solver = createLinearSolver(CGNR, op; AHA=normalOperator(op), iterations=40) + reshape(solve!(solver, vec(y)), n, n) +end; + +# ## One ray per cell +# +# 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) + +# 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 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); + +#- + +c, c_multi = colsums(op), colsums(op_multi) #hide +m = sum(c[49:80, 49:80]) / 32^2 #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) + +# 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 +save("multiray_reconstructions.png", fig) #hide +nothing #hide +# ![Reconstructions with one and with four rays per cell](multiray_reconstructions.png) + +# 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 +# error). For a 3D image, `RegularPoints(4)` places 4 × 4 = 16 rays per cell, and the time grows +# accordingly. +# +# ## Other point strategies +# +# Besides `RegularPoints`, the rays of a cell can end at random offsets, fixed by a seed (default 0): +# +# * `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 +# with a standard deviation of 0.3 cells, which can also land in neighboring cells (blur). +# +# `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: + +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", 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, jittered and random rays](multiray_jittered.png) + +# 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 +# +# 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))); + +# On the data above the error is 0.27, as with uniform weights. +# +# ## Cost and devices +# +# * 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 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/howto/implement_operator.jl b/docs/literate/howto/write_projector.jl similarity index 95% rename from docs/literate/howto/implement_operator.jl rename to docs/literate/howto/write_projector.jl index 7b69d8c..956ae64 100644 --- a/docs/literate/howto/implement_operator.jl +++ b/docs/literate/howto/write_projector.jl @@ -1,9 +1,11 @@ -# # 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 # 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 new file mode 100644 index 0000000..989e93d --- /dev/null +++ b/docs/literate/howto/write_sampling.jl @@ -0,0 +1,156 @@ +# # 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 `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 +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). 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 + 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, 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, 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 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) = 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`. 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) + +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 +# +# `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 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 +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 + +# `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 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 +# +# | 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/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/docs/make.jl b/docs/make.jl index 2809b5b..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(; @@ -38,24 +40,46 @@ 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", - "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", - "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", + "Comparing Projectors" => "generated/explanation/comparing.md", + ], + "Operators" => [ + "Masks" => "generated/explanation/masks.md", + "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.md b/docs/src/api.md deleted file mode 100644 index 74b9c1f..0000000 --- a/docs/src/api.md +++ /dev/null @@ -1,77 +0,0 @@ -# API - -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. - -## 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 - -```@autodocs -Modules = [CTOperators] -Pages = ["geometry/geometry.jl"] -Private = false -``` - -### Geometry types - -```@autodocs -Modules = [CTOperators] -Pages = ["geometry/raylist.jl", "geometry/conebeam.jl", "geometry/arcconebeam.jl", "geometry/parallelbeam.jl"] -Private = false -``` - -### Trajectories - -```@autodocs -Modules = [CTOperators] -Pages = ["geometry/trajectories.jl"] -Private = false -``` - -## Operators - -### Generic interface - -```@docs -AbstractCTOperator -``` - -```@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"] -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/geometries.md b/docs/src/api/geometries.md new file mode 100644 index 0000000..377ea69 --- /dev/null +++ b/docs/src/api/geometries.md @@ -0,0 +1,43 @@ +# Projection Geometries + +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. + +## Interface + +```@autodocs +Modules = [CTOperators] +Pages = ["geometry/geometry.jl"] +Private = false +``` + +## Geometry types + +```@autodocs +Modules = [CTOperators] +Pages = ["geometry/raylist.jl", "geometry/conebeam.jl", "geometry/arcconebeam.jl", "geometry/parallelbeam.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, together with the functions that custom strategies implement (see +[Write a Sampling Strategy](@ref)). + +```@autodocs +Modules = [CTOperators] +Pages = ["geometry/multiray.jl"] +Private = false +``` + +## Trajectories + +Constructors for the common scans. + +```@autodocs +Modules = [CTOperators] +Pages = ["geometry/trajectories.jl"] +Private = 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..657bae4 --- /dev/null +++ b/docs/src/api/internals.md @@ -0,0 +1,55 @@ +# Internals + +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 + +```@autodocs +Modules = [CTOperators] +Pages = ["operators/operators.jl", "operators/launch.jl"] +Public = false +``` + +## Masks of a projector + +```@docs +CTOperators.ProjectorMasks +CTOperators.projector_masks +CTOperators.is_masked +``` + +## Geometries + +```@autodocs +Modules = [CTOperators] +Pages = ["geometry/geometry.jl", "geometry/raylist.jl", "geometry/conebeam.jl", "geometry/arcconebeam.jl", + "geometry/parallelbeam.jl", "geometry/trajectories.jl", "geometry/multiray.jl"] +Public = false +``` + +## Image grid + +```@autodocs +Modules = [CTOperators] +Pages = ["image/image.jl", "image/voxelgrid.jl"] +``` + +## Ray tracing + +```@autodocs +Modules = [CTOperators] +Pages = ["raytracing/raytracing.jl", "raytracing/clip.jl", "raytracing/actions.jl", "raytracing/subrays.jl"] +``` + +## 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..68c0b90 --- /dev/null +++ b/docs/src/api/projectors.md @@ -0,0 +1,21 @@ +# Projectors + +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 + +```@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..7b64027 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -52,14 +52,17 @@ 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"). + [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), and + the projection methods, [Siddon's Method](@ref) and [Joseph's Method](@ref), + [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 [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 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 new file mode 100644 index 0000000..058e322 --- /dev/null +++ b/docs/utils/multiray_figures.jl @@ -0,0 +1,68 @@ +# 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; 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) + 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), 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] + scatter!(ax, first.(offs), last.(offs); markersize=14 .* sqrt.(w ./ maximum(w)), color=:darkorange) + end + fig +end 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 diff --git a/src/CTOperators.jl b/src/CTOperators.jl index c692d69..13e504b 100644 --- a/src/CTOperators.jl +++ b/src/CTOperators.jl @@ -12,11 +12,19 @@ export AbstractCTOperator export ImageGeom, MaskCircle, MaskEllipse, MaskAllButEdge export AbstractProjectionGeometry, AbstractDetectorGeometry, RaySegment, RayLine export RayList, ConeBeam, ArcConeBeam, ParallelBeam +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 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, one_layer, cell_ray, is_masked, visit, Gather, Scatter, to_index, clip, VoxelGrid +public image_mask, ray_mask, handles_raymask, projector_masks, ProjectorMasks + """ AbstractCTOperator{T} 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..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 @@ -113,12 +114,25 @@ 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::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...) for m in ms size(m, 1) == 3 || throw(DimensionMismatch("$name: per-view vectors must be 3 × nviews, got $(size(m))")) @@ -161,3 +175,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..a64d488 --- /dev/null +++ b/src/geometry/multiray.jl @@ -0,0 +1,385 @@ +# Several rays per detector cell: a detector geometry wrapped with a point strategy (where in the +# 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} + +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` and `sv` along `v`: +`su * sv` rays per cell. Point strategy of [`MultiRay`](@ref). +""" +struct RegularPoints{K} <: AbstractPoints{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 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} + su::Int + sv::Int + seed::UInt32 +end + +""" + RandomPoints(n; seed=0) + +`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 +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 offsets are fixed like those of [`JitteredPoints`](@ref). Point strategy of +[`MultiRay`](@ref). +""" +struct GaussianPoints{K,T} <: AbstractPoints{K} + sigma_u::T + sigma_v::T + seed::UInt32 +end + +""" + Points(offsets) + +Given points, the same in every detector cell: column `k` of the `2 × K` matrix `offsets` is the +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} + 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 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))")) + 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(::AbstractPoints{K}) where {K} = K + +""" + points_eltype(points, T) + +`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)) + +""" + point_offset(points, T, iu, iv, view, k) -> (du, dv) + +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 + +@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 + +# 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]) +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 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 + 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 rays of a detector cell count equally. Weighting strategy of [`MultiRay`](@ref). +""" +struct UniformWeights <: AbstractWeights 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} <: AbstractWeights + 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(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) + +""" + subray_weight(weights, du, dv) + +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 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 + +@inline restrict(::FullCell, du, dv) = (du, dv) +@inline restrict(::InPlane, du, dv) = (du, zero(dv)) + +""" + 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 offsets given by the point strategy `points` +([`RegularPoints`](@ref), [`JitteredPoints`](@ref), [`RandomPoints`](@ref), +[`GaussianPoints`](@ref), [`Points`](@ref), or a custom [`AbstractPoints`](@ref)) and +are averaged with the weights of `weights` ([`UniformWeights`](@ref), +[`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 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 reduces the gaps between rays that occur when +detector cells are large compared to the voxels. The data have the size of `geometry`. + +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)): 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 +```julia +g = MultiRay(circular_cone_beam(; angles, dso, dsd, nu, nv, du), JitteredPoints(2, 2; seed=1)) +op = JosephProjector(image, g) +``` +""" +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::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)))")) + 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) + +""" + 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 images with one layer of voxels +along `z`, see [`CTOperators.one_layer`](@ref). +""" +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 +""" + one_layer(image::ImageGeom) + +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 + +check_in_plane(p::Union{RegularPoints,JitteredPoints}) = + 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) + 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]) + +# 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 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 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 + +# 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) + 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) + +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/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/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..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. @@ -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 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} - backend = KA.get_backend(storage_array(geometry)) - JosephProjector(image, VoxelGrid{T}(image), geometry, backend, launch, ProjectorMasks(image, geometry, raymask)) + 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 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/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/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..e8d81f2 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 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} - backend = KA.get_backend(storage_array(geometry)) - SiddonProjector(image, VoxelGrid{T}(image), geometry, backend, launch, ProjectorMasks(image, geometry, raymask)) + 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 geometry(op::SiddonProjector) = op.geometry diff --git a/src/raytracing/actions.jl b/src/raytracing/actions.jl index 50a9731..1d95aa0 100644 --- a/src/raytracing/actions.jl +++ b/src/raytracing/actions.jl @@ -1,14 +1,30 @@ -# 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. +`image` is indexed linearly. +""" 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. +`image` is indexed linearly. +""" 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/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..84a9f1a --- /dev/null +++ b/src/raytracing/subrays.jl @@ -0,0 +1,39 @@ +""" + 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 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 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) + for k in 1:nsubrays(g) + 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 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) + for k in 1:nsubrays(g) + wsum += cell_weight(g, cell, k) + end + for k in 1:nsubrays(g) + a, b, w = cell_subray(g, cell, k) + segment(Scatter(action.image, action.value * (w / wsum)), grid, a, b, raytype(g)) + end + zero(T) +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), diff --git a/test/geometry/multiray.jl b/test/geometry/multiray.jl new file mode 100644 index 0000000..bfc1eed --- /dev/null +++ b/test/geometry/multiray.jl @@ -0,0 +1,82 @@ +# 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 + +# 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 + 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..a72e86f --- /dev/null +++ b/test/operators/multiray.jl @@ -0,0 +1,148 @@ +# 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)), (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)) + 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))) + # 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 +# 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