Reconstructing a volume from its X-ray projections

The registration problem gave us a volume and asked where it was. Reconstruction takes away the volume. We keep the calibrated source and detector positions, supply several radiographs, and estimate the attenuation field that could have produced them. The forward calculation is familiar, but the unknown now has one coefficient per voxel, which rather changes the storage bill.

Chapter 4 defined the sampled field and Chapter 5 its discrete derivatives. We will use those operators directly, first with a primary monochromatic count model. Chapter 8’s spectral model then lets us ask a harder question: which material mixtures are consistent with energy-resolved measurements?

We begin with a nonnegative voxel field, a count likelihood and a quadratic spatial penalty. These choices give a complete reconstruction loop whose gradient and constraints we can derive explicitly. We then replace each attenuation coefficient with material fractions and reuse the same projection machinery.

12.1 Specify the unknown volume and measured projections

Let $\boldsymbol\mu\in\mathbb{R}^{N}$ contain nonnegative attenuation coefficients in $\mathrm{mm}^{-1}$. They are the samples of the trilinearly interpolated field from Chapter 4. The grid shape is $(n_z,n_y,n_x)$, x varies fastest in memory, and its origin is the first sample centre. Spacing is specified along the grid’s x, y and z axes. An orientation matrix places those axes in object coordinates, and a fixed object-to-world transform places the volume in the calibrated acquisition.

For view $k$, write the predicted optical depths as

$A_k$ includes ray intersections with the volume support, the chosen midpoint quadrature and interpolation weights. Its entries have units of millimetres. $\mathbf n_{0,k}$ contains the calibrated expected open-beam counts, and $\boldsymbol\lambda_k$ is the predicted mean count image. The observed integer counts are $\mathbf y_k$. Conditioning on the open-beam calibration treats its uncertainty as negligible. A noisy flat-field estimate would require another part of the statistical model.

This model assumes independent photon counts, primary transmission and a monochromatic beam or an effective-energy approximation justified for the acquisition. An energy-integrating detector’s processed pixel values do not become Poisson counts because we rename the file. Scatter, electronic noise and spectral hardening can all produce structured discrepancies that the optimiser will otherwise try to explain with incorrect attenuation.

The inputs are calibrated geometry, observed counts, open-beam means and an initial attenuation field in inverse millimetres. A uniform field is a useful starting point because it leaves the measurements and regulariser to determine internal structure.

Figure 12.1 separates known geometry from unknown coefficients. A coefficient whose interpolation support misses every measured ray has a zero column in the combined projector. No choice of optimiser can make those measurements constrain it.

12.2 Derive the discrete volume gradient

For the fixed acquisition, $A_k$ is linear in the voxel coefficients. The exponential that follows it is nonlinear. A detector-space objective $\ell_k$ supplies a cotangent $\mathbf s_k=\partial\ell_k/\partial\boldsymbol\lambda_k$, giving

We use ordinary Euclidean inner products on the stored coefficient and detector vectors. Consequently, $A_k^{\mathsf T}$ is the matrix transpose of the implemented discrete map. A continuous backprojection formula, or the transpose of a different interpolation rule, need not give this derivative. If we instead defined volume and detector inner products through mass matrices $M_v$ and $M_d$, their weighted adjoint would be $M_v^{-1}A_k^{\mathsf T}M_d$. Those weights cannot be added twice, once in the objective and again as an unexplained voxel-volume factor.

Each reverse ray revisits its quadrature samples and scatters the seed through the same interpolation weights. Support clipping and the constant extension between the outer sample centres and half-cell faces remain the forward operator’s rules. Siddon’s radiological-path calculation is a useful reference for a different, cellwise-constant discretisation, but using its traversal does not make it the exact transpose of our sampled trilinear field. [1]

To implement the reverse pass, activate volume derivatives and hold pose fixed. Clear the voxel gradient once, then add each view’s contribution with projection_vjp. Recomputing interpolation weights avoids storing either a system matrix or a ray-by-sample tape. The listing composes this transpose operation with the count-space derivative:

    def evaluate(self, volume: Any, volume_wp: Any, *, gradient: bool) -> float:
        """Evaluate every supplied view before an update; add the prior once."""
        if gradient:
            self.gradient.zero_()
        loss = 0.0
        for view in self.views:
            project_optical_depth(
                volume_wp,
                self.pose,
                workspace=view["projection"],
                out_L=view["depth"],
                stream=self.stream,
            )
            transmit(
                view["depth"],
                view["beam"],
                workspace=view["transmission"],
                out_counts=view["prediction"],
                stream=self.stream,
            )
            evaluate_objective(
                view["prediction"],
                view["observed"],
                workspace=view["objective"],
                out_loss=view["loss"],
                out_seed=view["image_seed"] if gradient else None,
                stream=self.stream,
            )
            if gradient:
                transmission_vjp(
                    view["depth"],
                    view["beam"],
                    seed_counts=view["image_seed"],
                    workspace=view["transmission"],
                    out_grad_L=view["depth_seed"],
                    stream=self.stream,
                )
                projection_vjp(
                    volume_wp,
                    self.pose,
                    adj_L=view["depth_seed"],
                    workspace=view["projection"],
                    out_mu=self.gradient_wp,
                    accumulate=True,
                    stream=self.stream,
                )
            loss += float(view["loss"].numpy()[0])
        loss += self.regularise(volume, gradient=gradient)
        if not math.isfinite(loss):
            raise NumericalError("The objective is non-finite; inspect counts and model range.")
        if gradient and not bool(self.torch.isfinite(self.gradient).all().item()):
            raise NumericalError("The volume gradient overflowed; discard this evaluation.")
        return loss

Figure 12.2 follows these quantities through both directions. Before accepting numerical results, test $\langle A\mathbf u,\mathbf v\rangle=\langle\mathbf u,A^{\mathsf T}\mathbf v\rangle$ with independently chosen vectors and check the composed objective along volume perturbations. A reverse pass that merely returns finite numbers has met a rather modest ambition.

12.3 Choose the likelihood and regularisation

For independent Poisson observations, the code minimises the half-deviance, which differs from the negative log-likelihood by an observation-only constant:

For $y=0$, the contribution is $\lambda$, and we never evaluate $\log 0$ as an intermediate observation term. For positive $y$, the mean must be positive. The count derivative is $1-y/\lambda$, so the algebraically composed optical-depth derivative simplifies to $y-\lambda$. Thus

Under these fixed-geometry primary-count assumptions, the likelihood is convex in attenuation. That does not guarantee strict convexity or uniqueness. A null space survives wherever different coefficient fields have identical measured projections. Spectral mixing and joint geometry estimation will not generally retain this convenient structure.

The factored chain evaluates $1-y/\lambda$ before multiplying by $-\lambda$. At very small predicted means, that intermediate can overflow even when the simplified depth derivative $y-\lambda$ is representable. Check the numerical range when composing separate operators; where needed, evaluate the simplified expression directly or retain higher-precision intermediates.

Any voxel change in the likelihood’s null space is invisible to the measurements. A spatial penalty states which possible fields we prefer. Here we favour smaller differences between neighbouring coefficients, using quadratic neighbouring-voxel regularisation:

$\mathcal E_a$ contains each neighbouring pair along axis $a$ once, $h_a$ is that axis’s spacing and $\Delta V=h_xh_yh_z$. Edges beyond the grid are omitted. This adds no boundary penalty against exterior zero, even though the projector’s field is zero outside its support. $R$ has units of $\mathrm{mm}^{-1}$, so $\beta$ has units of millimetres when the count objective is dimensionless. Replacing millimetres with centimetres therefore changes the numerical regularisation coefficient as well as the grid metadata.

For one edge, the regulariser adds $\beta\Delta V(\mu_i-\mu_j)/h_a^2$ to voxel $i$ and its negative to voxel $j$. The implementation uses FP64 differences and loss reductions, with FP32 updates to the stored volume gradient.

    def regularise(self, volume: Any, *, gradient: bool) -> float:
        """Quadratic physical-space differences; omit edges beyond the grid."""
        self.volume64.copy_(volume)
        penalty = 0.0
        for low, high, grad_low, grad_high, scratch, coefficient in self.edges:
            difference, derivative = scratch
            self.torch.sub(high, low, out=difference)
            if gradient:
                # Scale in FP64 before the FP32 accumulation boundary.
                difference.mul_(coefficient)
                derivative.copy_(difference)
                grad_low.sub_(derivative)
                grad_high.add_(derivative)
                # Restore unscaled differences for the objective.
                self.torch.sub(high, low, out=difference)
            difference.square_()
            self.torch.sum(difference, dim=(0, 1, 2), out=self.scalar)
            penalty += 0.5 * coefficient * float(self.scalar.item())
        return penalty

Quadratic smoothing penalises sharp genuine boundaries alongside noise. Total variation or a learned prior would change that bias and require a different solver or derivative. Here, nonnegativity is enforced by projection after an update. Unlike an exponential parameterisation, projection permits exact zero attenuation without driving a latent variable towards minus infinity.

Sum the likelihood over the included pixels. Treat missing and saturated measurements through masking or censoring, while retaining zero counts as valid observations. Choose $\beta$ to balance data agreement against spatial smoothing; inspect the resulting boundaries as well as the scalar loss.

12.4 Plan the GPU reconstruction iteration

Keep the volume, trial field, gradient and update displacement on the GPU throughout the solve. The same applies to observations and per-view derivative buffers. In the accompanying implementation, Torch and Warp share array storage and one CUDA stream. This lets us compose their operations without copying whole fields between them.

An iteration evaluates every view at the accepted volume, adds regularisation once, and chooses a projected step. With $\mathbf g=\nabla\Phi(\boldsymbol\mu)$,

The line search accepts strict decrease satisfying $\Phi(\boldsymbol\mu_\alpha)\leq\Phi(\boldsymbol\mu)+c\,\mathbf g^{\mathsf T}\mathbf d_\alpha$, for $0<c<1$. It reduces $\alpha$ geometrically when this condition fails. The gradient has units of millimetres, so $\alpha$ has units of $\mathrm{mm}^{-2}$. The stopping quantity is $\|\mathbf d_{\alpha_0}\|_\infty/\alpha_0$ at a fixed reference step, measured in millimetres. A small objective change alone does not establish constrained stationarity.

    def solve(
        self, callback: Callable[[int, dict[str, Any]], None] | None = None
    ) -> dict[str, Any]:
        """Solve with optional observation of accepted updates only.

        The callback receives the stable one-based accepted iteration and a
        detached history dictionary after ``mu`` has been updated. It may read
        or export the accepted fields but must not mutate solver state. Callback
        exceptions propagate with the last accepted volume intact.

        An optional fixed mapping step separates the stationarity diagnostic
        from the initial Armijo trial. Omitting it preserves the original
        diagnostic at ``initial_step_mm_inverse_squared``.
        """
        policy = self.policy
        mapping_step = (
            policy.initial_step_mm_inverse_squared
            if policy.mapping_step_mm_inverse_squared is None
            else policy.mapping_step_mm_inverse_squared
        )
        history: list[dict[str, Any]] = []
        reason = "iteration_budget"
        for iteration in range(policy.iterations):
            loss = self.evaluate(self.mu, self.mu_wp, gradient=True)
            step = policy.initial_step_mm_inverse_squared
            mapping, _ = self.displacement(mapping_step)
            if not math.isfinite(mapping):
                raise NumericalError("The trial update overflowed; reduce the initial step.")
            if mapping <= policy.gradient_mapping_tolerance_mm:
                reason = "projected_gradient_tolerance"
                break
            accepted = False
            for attempt in range(policy.maximum_backtracks):
                if step == 0:
                    break
                _, slope = self.displacement(step)
                if not math.isfinite(slope) or slope >= 0:
                    step *= policy.backtracking_factor
                    continue
                # Trial value evaluation leaves the accepted gradient unchanged.
                # Domain/range errors abort with the last accepted volume intact.
                trial_loss = self.evaluate(self.trial, self.trial_wp, gradient=False)
                if trial_loss < loss and trial_loss <= loss + policy.armijo * slope:
                    self.mu.copy_(self.trial)
                    history.append(
                        {
                            "iteration": iteration + 1,
                            "loss_before": loss,
                            "loss_after": trial_loss,
                            "step_mm_inverse_squared": step,
                            "backtracks": attempt,
                            "mapping_mm_before": mapping,
                        }
                    )
                    accepted = True
                    if callback is not None:
                        callback(iteration + 1, dict(history[-1]))
                    break
                step *= policy.backtracking_factor
            if not accepted:
                reason = "line_search_failed"
                break
        # Trial buffers may describe a rejected volume; refresh accepted outputs.
        final_loss = self.evaluate(self.mu, self.mu_wp, gradient=True)
        final_mapping, _ = self.displacement(mapping_step)
        if final_mapping <= policy.gradient_mapping_tolerance_mm:
            reason = "projected_gradient_tolerance"
        return {
            "termination": reason,
            "accepted_steps": len(history),
            "final_objective": final_loss,
            "final_gradient_mapping_mm": final_mapping,
            "history": history,
        }

Keep the accepted field and gradient intact while evaluating line-search trials. Once the solve stops, recompute projections at the accepted field so the displayed images cannot come from the last rejected trial.

Allocate the working fields before entering the optimisation loop. The main storage grows with the voxel count and retained detector images, rather than with the number of iterations. The host needs only scalar losses and line-search decisions; transferring a whole volume after every trial would add work without helping the optimisation.

When all observations fit on the GPU, retain them for repeated use. Otherwise process bounded groups, accumulating the complete gradient before updating the volume. Updating after every group defines a different optimisation schedule, with its own step-size and convergence behaviour.

12.5 Reconstruct with fewer or restricted views

Sparse angular sampling leaves gaps between otherwise broadly distributed views. Limited-angle acquisition restricts the directions themselves. Both can produce an apparently stable minimiser, but they constrain different structures poorly. Adding more views inside a narrow angular interval does not supply the missing directions outside it.

DiffVox also reconstructs a voxel grid through differentiable rendering. Its formulation combines an image-domain L1 loss with total variation, and uses Adam to optimise the grid. [37] We instead use a Poisson likelihood, a quadratic penalty and a projected line search. The grid representation is shared; the statistical model and update rule are explicit choices. Total variation favours piecewise-constant structure, whereas our quadratic penalty produces a smoother objective and the neighbour-difference gradient derived above.

With fewer views, each full-gradient iteration costs less. To compare acquisitions, hold the total photon budget fixed and examine reconstruction error together with numerical progress. This separates the effect of missing directions from the cost of reaching a solution.

12.6 Extend reconstruction to spectral material decomposition

A single attenuation field cannot explain how the same object changes appearance across distinct spectra. Represent attenuation instead as a combination of energy-dependent basis curves. Alvarez and Macovski developed the energy-selective basis formulation. The inverse problem concerns basis coefficients, whose physical interpretation depends on the chosen model. [11]

For a mixture of reference materials, use volume fractions $f_m(\mathbf x)$, with $f_m\geq0$ and $\sum_m f_m\leq1$. An equivalent-material basis is another choice: its coefficients scale reference attenuation curves and may require only nonnegativity. Choose the domain to match that physical interpretation. With linear reference attenuation $a_m(E)$ in $\mathrm{mm}^{-1}$,

The material paths $b_{kpm}$ are in millimetres. For this count-channel model, $s_{kce}$ is the expected incident photon population integrated over energy bin $e$, multiplied by the calibrated probability of detection in channel $c$. Energy integration or quadrature weights are already absorbed into that population, so multiplying by a bin width again would count it twice. An energy-weighted signal requires different response units and generally a different noise model. Unassigned fraction means zero contribution from the represented attenuators. It is not an implicit, calibrated air material.

Convert mass attenuation in $\mathrm{cm}^2/\mathrm g$ to linear attenuation by multiplying by the reference density in $\mathrm g/\mathrm{cm}^3$ and dividing by ten for $\mathrm{mm}^{-1}$. Curves already expressed in linear units include that density. With these linear reference curves, the reconstructed coefficients are dimensionless fractions or equivalent-material amounts. A mass-density basis instead pairs density fields with mass attenuation curves and compatible length units.

The per-ray material sensitivity is

Locally stable first-order separation of unrestricted material paths requires sufficient rank in this sensitivity matrix. Having as many channels as materials is only a counting condition: nearly proportional columns still make their separation unstable. Inspect the noise-weighted singular values after declaring coefficient scales, across representative material thicknesses. Spectral hardening changes those sensitivities with the object itself. Independent Poisson channels give Fisher weights $1/\lambda_c$, while detector effects that correlate channels require the appropriate likelihood instead. Alvarez’s analysis of noninvertibility gives a reason to examine the map beyond a favourable local derivative. [12]

Projection-domain decomposition first estimates material paths on corresponding rays, then reconstructs each field. Its convenience depends on genuinely corresponding acquisitions, and the intermediate estimates carry correlated, nonuniform uncertainty. Joint material-volume reconstruction instead composes the material projector and spectral adjoint directly, with the chosen field constraints and spatial priors. The fraction example uses the simplex constraint. It can combine different view geometries, but remains dependent on their combined information. Priors may couple otherwise ambiguous coefficients: that is additional information assumed by the reconstruction.

Changing the unknown also changes the regulariser’s units. For dimensionless material coefficients, the physical-gradient construction in §12.3 has units of millimetres, so its multiplier has units $\mathrm{mm}^{-1}$. For attenuation coefficients, that construction has units $\mathrm{mm}^{-1}$ and its multiplier has units of millimetres. Copying the monochromatic regularisation number into the material problem would quietly change what is being penalised.

The two material sensitivity columns can differ greatly in scale and be nearly proportional. A Euclidean step limited by the strongly constrained combination can then make slow progress along the weak combination. The material solver therefore permits a symmetric positive definite matrix $H$ in its projection. With a positive spatial factor $D_i$, the metric at voxel $i$ is $H_i=D_iH$. It changes which feasible update we propose; the objective remains the summed Poisson half-deviance and the spatial penalty.

For the fraction problem, let $\mathcal C=\{\mathbf u\in\mathbb R^2:\mathbf u\geq0,\ \mathbf1^{\mathsf T}\mathbf u\leq1\}$ and $\mathbf g_i=\nabla_{\mathbf f_i}\Phi$. The ordinary metric trial is

where $\Pi^{H_i}$ denotes projection in the $H_i$ norm. For a fixed gradient and step, larger metric curvature reduces the unconstrained move in that direction; the projection enforces feasibility. The line search then decides whether to accept this trial.

Spatial scaling matters because rays through different material thicknesses carry different count sensitivity. At the accepted field, let $\mathbf j_{kpc}$ be the two-component material-path derivative above. The expected Poisson Fisher block is $B_{kp}=\sum_c\mathbf j_{kpc}\mathbf j_{kpc}^{\mathsf T}/\lambda_{kpc}$. Choose $\kappa_{kp}$ as the largest eigenvalue of $H^{-1/2}B_{kp}H^{-1/2}$. The nonnegative interpolation weights give the spatial row bound

For one ray with weights $a_i\geq0$, weighted Cauchy–Schwarz gives $\|\sum_i a_i\mathbf u_i\|_H^2\leq(\sum_i a_i)\sum_i a_i\|\mathbf u_i\|_H^2$. Combining this with $B\preceq\kappa H$ and summing rays produces the backprojected row bound. For the quadratic regulariser, let $d_i$ be the sum of incident edge weights $\beta\Delta V/h_a^2$ at voxel $i$. Adding $2d_i/\lambda_{\min}(H)$ to the raw scale bounds that regulariser’s curvature as well.

Apply a positive floor to the spatial scale and divide by its mean to obtain $\mathbf D$. This leaves the mean step size to $\alpha$ while retaining relative scaling across voxels. The curvature bound above applies to the raw scale. Because the Fisher matrix approximates local curvature, the line search still decides whether a proposed update decreases the actual objective.

The metric sets relative update sizes; the Barzilai–Borwein rule supplies an overall length from the change in gradient over the preceding step. The original two-point construction approximates a secant equation using one previous iterate and gradient. [42] Applied in fixed metric coordinates, it gives $\alpha=(\sum_iD_i\mathbf s_i^{\mathsf T}H\mathbf s_i)/(\mathbf s^{\mathsf T}\mathbf y)$, where $\mathbf s$ and $\mathbf y$ are consecutive field and gradient differences. We use this as a proposal inside the constrained line search. If the denominator is nonpositive or nonfinite, fall back to ordinary step growth.

The optional inertial proposal also extrapolates from the preceding accepted field before metric projection. Extrapolation can continue useful progress, but it can overshoot. The line search therefore tests the true objective with strict Armijo decrease, using the gradient at the current accepted field. If an inertial proposal fails, it retries the same step length without extrapolation.

When refreshing the metric, reset the two-point and inertial histories: differences measured under the previous scaling no longer describe the same step geometry.

The material loop therefore consists of a forward and reverse sweep, a scaled feasible proposal, and a line search. Retain the fields and proposal histories on the GPU, returning only scalar reductions for step selection. Refresh the spatial metric between stretches of iterations, since computing it after every update would add another substantial projection calculation.

Stationarity is assessed separately: at the fixed diagnostic step $\alpha_0$, compute $\|\mathbf f-\Pi_{\mathcal C}(\mathbf f-\alpha_0\nabla\Phi)\|_\infty/\alpha_0$, with ordinary Euclidean projection onto the per-voxel simplex $\mathcal C$. Neither $H$ nor $\mathbf D$ enters this diagnostic. Otherwise, making the update metric sufficiently large could manufacture a small reported step while leaving the same unresolved gradient. That would be an admirably cheap reconstruction, provided nobody expected it to be finished.

Figure 12.3 follows both material fields from the uniform start and compares them with the assigned reference. Inspect the materials together: a weakly constrained combination can exchange one component for another while leaving a small image residual.

For recorded mean signals, the effective spectral weights can also vary over detector pixels. Writing $K_{kpce}$ for the effective spectral kernel in stored-signal units gives

$K$ includes the integrated energy-bin populations and spatial open-beam factors. Restrict the sum to entries selected by the fixed validity mask $m_{kpc}$ before computing residuals. Choosing $w$ as inverse squared predicted air signal normalises the residual scale; inverse noise-covariance weighting would require a statistical model for the measured channels. Cumulative channels share photons and are therefore correlated.

12.7 Validate the volume independently

Assess the reconstruction in two domains. Reserved views test how well it predicts unseen projections; a known reference field tests material or attenuation recovery in physical units. The worked example displays both, since a good projection fit can still leave uncertainty inside the volume.

While building the projector, compare its forward values and directional derivatives with an independent integral on simple mathematical fields. These checks test the discretisation before it becomes part of a large inverse problem. Then use a different quadrature or independently implemented projector when constructing a numerical reference.

For acquired data, separate plausible contributors to error through controlled studies: count noise with fixed calibration, perturbed geometry with otherwise fixed inputs, and forward-model mismatch such as unmodelled hardening. Report attenuation bias and spatial resolution with task-relevant quantities, not only an image similarity score. A smoothed volume can look agreeably quiet while erasing the structure the reconstruction was supposed to recover.

Repeated acquisitions or independent noise replicates show how much the recovered field varies. A local Hessian gives a complementary approximation around one fitted solution, conditional on the chosen measurement model and prior.

12.8 Complete the reconstruction study

The scalar algorithm is now complete: render all views, pull the count residuals back to the volume, add the spatial gradient, and accept a projected step. For the material version, replace the scalar transmission with the spectral mixture and project each voxel onto its material domain. The ray geometry and interpolation code remain the same.

The material worked case uses a CT-derived water/bone assignment. The anatomy is averaged to a 16³ field, which defines both the simulated object and reconstruction grid. Its fractions are assigned model values rather than measured patient composition.

Use 48 fitting views on three full-azimuth orbits at tilts of −15°, 0° and 15°, with a 48 × 48 detector covering the box. Reserve twelve views at different yaw angles for prediction. This keeps the example small enough to inspect both the field and the optimisation history.

Start both noise replicates from uniform water. Prepare the spatial metric there, then refresh it after 100 accepted updates so its weighting can reflect the evolving field.

The material solver supplies the joint reconstruction step for specified acquisitions. Finish both fits before comparing either with the known reference field.

Both fits start from uniform water and meet the projected-gradient stopping criterion defined in §12.6. In the first fit, the objective is already within 1% of its final value after 100 updates, while the projected-gradient measure is still far above its stopping threshold. A flat loss curve is therefore insufficient reason to stop.

Compare the completed fields with their assigned reference and render the reserved views. Using generating mean images for that comparison measures prediction error without adding another observation’s count noise.

Armijo requires the objective to decrease at accepted updates; the projected-gradient measure need not decrease at every update. Use the viewer to compare materials, replicates and slices. Its fraction and signed-error windows remain fixed.

The recovered fields agree with the coarse assignment and predict the reserved views. Noise and the spatial penalty leave a nonzero objective. To explore the regulariser’s effect, vary its strength and watch the trade-off between smooth regions and sharp boundaries.

Throughout this chapter, someone else chose the acquisition. Chapter 13 makes that choice part of the problem: which next view, and which allocation of exposure, would improve a specified registration or reconstruction task?