## 1.

## Introduction

Cone-beam computed tomography (CBCT) finds application in many clinical interventions to provide up-to-date imaging (often in conjunction with hardware/instrumentation) for improved localization and assessment of treatment delivery.^{1}^{–}^{5} Intraoperative CBCT imaging platforms are often based on C-arm devices with flat-panel detectors that are capable of radiography and fluoroscopy in addition to the three-dimensional (3-D) CBCT. Advanced versions of these systems are motorized on several axes to permit nonisocentric source–detector orbits, to expand the field of view (FOV) in 3-D image reconstructions, and to quickly reposition between specific radiographic/fluoroscopic views. Such systems include floor- and ceiling-mounted C-arms as well as robotic C-arms^{6} with many degrees of motion freedom. In addition to C-arm systems, some CBCT mammography^{7} and body imaging systems^{8} are also capable of complex source–detector trajectories.

The additional flexibility provided by these systems permits more general orbits beyond the traditional circular and helical source–detector trajectories that have been the norm for CT for decades. To date, these flexible orbits have mainly been used to address the FOV and the sampling issues in interventional CBCT. For example, noncircular trajectories have been used to provide extended axial^{9} and elliptical^{10} FOVs and to improve 3-D sampling and data completeness^{11}^{–}^{14} to avoid cone-beam artifacts that arise from traditional circular cone-beam orbits.

Tilted circular orbits are commonly used for their ability to positively impact image quality and/or localization. For example, tilting the CT gantry relative to the patient’s longitudinal axis improves the image quality adjacent to the skull base,^{15} reduces eye lens dose,^{16} improves localization in CT-guided biopsies,^{17}^{,}^{18} and reduces metal artifacts associated with prostheses.^{19} These examples suggest that modifications of the orbit beyond simple tilts may also provide clinical advantages. However, selection of an optimal trajectory presents many challenges. For example, the above-mentioned simple tilt examples depend on the patient anatomy and/or the interventional procedure, e.g., aligning the gantry along the canthomeatal line, to the axis of the biopsy needle, etc. Thus, optimal trajectories are likely patient- and imaging-task-dependent. Moreover, data acquired from noncircular trajectories can be difficult to reconstruct, as the sampling conditions for traditional filtered-backprojection methods no longer apply.

There is a growing trend in the use of task-based measures in performance assessment.^{20}^{–}^{22} Such measures have also been used in prospective task-driven optimization of system design,^{23}^{–}^{25} regularization of model-based reconstruction,^{26} and CT data acquisition parameters like dual-energy imaging,^{27} tube current modulation,^{28} and fluence-field modulation.^{29} Thus, we expect that task-based performance models similarly provide a basis for optimizing the source–detector orbit in CBCT.

Interventional imaging presents an ideal opportunity to customize orbits to the patient and diagnostic task for a number of reasons: (1) most patients in interventional imaging have previous diagnostic imaging studies, meaning that a detailed representation of the patient anatomy is available; (2) additional information, including surgical plans, the location and sizes of implants or tools, and particular anatomical targets, is known prior to the intervention; and (3) interventional imaging tasks tend to be well defined, including the particular volume of interest and specific image features that need to be identified or localized.

In this work, we present a mathematical model that leverages prior knowledge, including anatomical dependencies, to predict the imaging task performance for different trajectories. This predictor is integrated into an optimization framework that seeks the scan trajectory that maximizes performance. This framework is applied in simulation studies illustrating behavior and performance in simple digital phantoms. This work is built on previous reports on task-driven orbits,^{30}^{–}^{32} with a comprehensive introduction to the framework as well as an improved implementation and application to optimization for multiple tasks varying in location and/or spatial-frequency content.

## 2.

## Methods for Task-Driven Trajectory Optimization

## 2.1.

### Overview and Proposed Imaging Workflow

CBCT image quality can be highly dependent on patient size, anatomical site, and the presence of interventional hardware in the FOV. Even within a scan for a single patient, the data fidelity can vary widely with orders of magnitude differences in the noise for different measurements. This work seeks to develop a new imaging workflow in which CBCT scans are driven by the particular patient anatomy, selecting projections that maximize data fidelity for a specific imaging task. In general, this requires some knowledge of patient anatomy for prospective trajectory design. Such information is often available in interventional imaging but is typically not used directly by the imaging device.

Figure 1 illustrates a proposed imaging workflow that leverages preoperative imaging and planning data and contrasts the proposed methodology with a conventional workflow. Conventionally, a patient’s diagnosis is obtained via CT (or other modality) to define and plan the particular interventional approach. In many procedures involving implantation of exogenous devices, this includes specification of the particular hardware required for the procedure (e.g., device size, location for deployment, and model number). Minimally invasive treatment increasingly relies on intraoperative imaging or immediate postoperative imaging for assessment of the surgical product and identification of possible complications. Intraoperative imaging has the advantage of detecting complications while the patient is still in the operating room with better opportunity to revise if necessary. Image quality for intraoperative and postoperative assessment is often challenged by the surgical tools and implants delivered during the procedure—often metallic and/or high density. Most importantly, image quality often suffers most in the vicinity of the implant, which, unfortunately, is the part of the image where complications are most likely to be found. In the conventional workflow, the imaging system ignores the wealth of knowledge about the patient anatomy, planned hardware delivery, and specifics of the imaging task. In contrast, a task-driven imaging workflow leverages this information to improve performance for pertinent imaging tasks associated with the particular intervention.

An overview of the task-driven optimization is illustrated in Fig. 2. The approach combines an anatomical model of the patient (as well as planning information), which is important for predicting the fidelity of the projection data, with an imaging system model, which accounts for the entire imaging chain, including any parameters to be optimized (e.g., source–detector trajectory). These models incorporate prediction of imaging performance, including local spatial resolution and noise in the reconstructed image volume. These measures of imaging performance may then be used to compute task performance using a particular observer model. With the ability to model the end-to-end system from data collection to observer performance, various parameters may be tuned in an iterative process to find the optimal source–detector trajectory that maximizes the imaging performance. The modeling and predictive framework in Fig. 2 are detailed in the following sections.

## 2.2.

### Task-Based Performance Prediction and Optimization

Various mathematical observer models have been used for performance prediction. The authors have previously demonstrated basic agreement between human observer performance and a non-prewhitening observer model over a broad range of imaging conditions^{33} and have therefore elected to use this model to evaluate imaging tasks. Alternative metrics could be considered in future work, e.g., a prewhitening observer model to examine fundamental signal and noise content or a channelized Hoteling observer model to potentially capture aspects more closely related to a human observer. A metric based on large-area transfer characteristics (e.g., contrast-to-noise ratio) captures only the low-frequency performance and would likely miss aspects related to spatial resolution and frequency response. The non-prewhitening model chosen here derives detectability in terms of the local spatial resolution and noise properties of the reconstructed image as well as a task function that specifies the spatial frequencies of interest. The detectability index for the non-prewhitening observer model may be written as

## Eq. (1)

$${d}_{j}^{\prime 2}(\mathrm{\Omega})=\frac{{\{\iiint {[{\mathrm{MTF}}_{j}(\mathrm{\Omega})\xb7{H}_{\text{Task}(j)}]}^{2}\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{d}{f}_{x}\text{\hspace{0.17em}}\mathrm{d}{f}_{y}\text{\hspace{0.17em}}\mathrm{d}{f}_{z}\}}^{2}}{\iiint {\mathrm{NPS}}_{j}(\mathrm{\Omega})\xb7{[{\mathrm{MTF}}_{j}(\mathrm{\Omega})\xb7{H}_{\text{Task}(j)}]}^{2}\mathrm{d}{f}_{x}\text{\hspace{0.17em}}\mathrm{d}{f}_{y}\text{\hspace{0.17em}}\mathrm{d}{f}_{z}},$$The detectability index provides an objective for the design of an optimal source–detector orbit. The most straightforward objective that seeks to optimize detectability index for a single location in the imaging volume may be written as

## Eq. (2)

$$\hat{\mathrm{\Omega}}={\mathrm{arg}\mathrm{max}}_{\mathrm{\Omega}\in {\mathrm{\Omega}}_{\text{feasible}}}{d}_{j}^{\prime 2}(\mathrm{\Omega}),$$Although single-location optimization may be appropriate for some imaging tasks, this objective does not consider performance at any other location in the image, thus leading to solutions that may be highly optimized to a single point and sacrifice image quality everywhere else in the image. As such, we also consider multilocation objectives. There are many possible choices for a multilocation objective; however, the principal concern is how to weigh the relative importance of performance at different locations in the FOV. We explore the following three choices in this work:

*Maximum mean detectability* (*maxi-mean*)—in which the average detectability index over an ensemble of locations within a specified region of interest (ROI), ${j}_{\mathrm{ROI}}$, is computed and maximized. Mathematically, this objective is written as

## Eq. (3)

$$\hat{\mathrm{\Omega}}={\mathrm{argmax}}_{\mathrm{\Omega}\in {\mathrm{\Omega}}_{\text{feasible}}}\text{\hspace{0.17em}}{\text{mean}}_{j\in {j}_{\mathrm{ROI}}}\{{d}_{j}^{\prime 2}(\mathrm{\Omega})\}.$$This objective treats all performance gains equally throughout the ROI, including potential solutions where detectability is decreased in one region for a larger gain in another region.

*Maximum median detectability* (*maxi-median*)—in which the median detectability index over the regional ensemble is maximized. Mathematically, it can be written as

## Eq. (4)

$$\hat{\mathrm{\Omega}}={\mathrm{argmax}}_{\mathrm{\Omega}\in {\mathrm{\Omega}}_{\text{feasible}}}\text{\hspace{0.17em}}{\text{median}}_{j\in {j}_{\mathrm{ROI}}}\{{d}_{j}^{\prime 2}(\mathrm{\Omega})\}.$$This is similar to the maxi-mean objective, except that it permits larger outliers. For example, if overall detectability can be increased at the cost of a significant decrease at a few locations, this would be deemed an advantage.

*Maximum minimum detectability* (*maxi-min*)—in which one seeks to achieve the highest minimum detectability over an ensemble of locations. Mathematically, it can be written as

## Eq. (5)

$$\hat{\mathrm{\Omega}}={\mathrm{arg}\mathrm{max}}_{\mathrm{\Omega}\in {\mathrm{\Omega}}_{\text{feasible}}}\text{\hspace{0.17em}}{\text{minimum}}_{j\in {j}_{\mathrm{ROI}}}\{{d}_{j}^{\prime 2}(\mathrm{\Omega})\}.$$This objective indicates that it is the location of minimum performance that drives the design, and that imaging performance cannot be sacrificed in one location (in the designated region of interest) for an improvement in another location.

Both the single-location and multilocation objectives are investigated below. To estimate the solution to any of these objective functions, we used the covariance matrix adaptation evolution strategy (CMA-ES) algorithm.^{34} In CMA-ES, a population sample is randomly drawn according to a multivariate normal distribution at each iteration. The best solutions in the population are used to estimate the local covariance matrix of the objective function in an adaptive manner. The mean, covariance matrix, step size, and evolutionary paths are then updated to generate the next population with the goal of maximizing the number of successful samples in each successive population. This is repeated until convergence, defined as the iteration beyond which changes in function evaluation are negligible. This stochastic approach is attractive, as only function evaluations are required and populations of solutions are used, which helps to find global optima in nonlinear and nonconvex objectives. Another attractive feature of CMA-ES is its improved robustness with increased function evaluations (e.g., increasing population size), which can be tuned to help avoid local optima.

## 2.3.

### Imaging System and Reconstruction Model

To use the performance objectives defined in the previous section, one must be able to predict imaging properties (viz., local spatial resolution and noise). This necessitates an imaging system model that includes all of the pertinent imaging physics and acquisition parameters as well as the reconstruction process. Starting with data acquisition, we adopt the following forward model for a CBCT system for which mean measurements are modeled as

## Eq. (6)

$${\overline{y}}_{i}={g}_{i}\text{\hspace{0.17em}}\mathrm{exp}\text{\hspace{0.17em}\hspace{0.17em}}(-{[\mathbf{A}(\mathrm{\Omega})\mu ]}_{i}),$$Traditional analytic reconstruction methods are challenged by orbits that deviate from standard designs (e.g., circular and helical). However, model-based reconstruction methods are straightforward to apply on unusual and even incomplete data orbits—providing the “best” possible estimates, given the data that was collected. Thus, we adopt a statistically motivated model-based reconstruction for this work. Presuming a Poisson noise model for the measurements, we may write the following penalized-likelihood estimator:

where $L(y;\mu )=\sum _{i=1}^{N}{y}_{i}\text{\hspace{0.17em}}\mathrm{log}({g}_{i}\text{\hspace{0.17em}}\mathrm{exp}(-{[\mathbf{A}(\mathrm{\Omega})\mu ]}_{i}))-{g}_{i}\text{\hspace{0.17em}}\mathrm{exp}(-{[\mathbf{A}(\mathrm{\Omega})\mu ]}_{i})$ and denotes the log-likelihood function, which is a function of the pre-log-transformed measurements $y$, and $R(\mu )=\beta {\mu}^{T}\mathbf{R}\mu $ represents a roughness penalty used in controlling the noise in the reconstruction with a strength parameter $\beta $. In this work, we focus on a quadratic penalty that is specified by the matrix $\mathbf{R}$, a constant matrix that defines how voxels are combined and penalized such that $R(\mu )=\frac{1}{2}\sum _{j}\sum _{k}{w}_{j,k}{({\mu}_{j}-{\mu}_{k})}^{2}$, where ${w}_{j,k}=1$ for the six nearest neighbors in 3-D space and 0 otherwise. Although there are many more sophisticated regularization schemes, this particular choice of roughness penalty is well suited to previously developed imaging performance predictors.^{35}

Another advantage of the penalized-likelihood framework is that arbitrary source–detector trajectories may be reconstructed without modifying the underlying algorithm used in solving Eq. (2). In contrast, many direct reconstruction approaches will not implicitly handle noncircular or nonhelical trajectories without substantial modification or rederivation of the algorithm. In this work, we have solved Eq. (7) iteratively using the ordered subsets, separable quadratic surrogate approach discussed in Ref. 36. At each iteration, all voxels are updated simultaneously, requiring one forward-projection and one backprojection to compute the likelihood gradient. The penalty gradient and curvature are computed directly from the image. In experiments described below, the penalized-likelihood estimator is run to convergence using a specific number of iterations and a zero image has been used for initialization.

## 2.4.

### Parameterization of the Source–Detector Trajectory

As described above, the system geometry associated with a particular source–detector trajectory was parameterized by the vector $\mathrm{\Omega}$. Parameterization of the trajectory can take many possible forms that depend on the particular capabilities of the imaging system. In this work, we considered two parameterizations of the orbit and concentrated on orbits that sample a sphere around a common center of rotation whose x-ray source positions are specified by the coordinate pair $(\theta ,\varphi )$ with rotation angle $\theta $ and tilt angle $\varphi $. Other geometrical parameters, such as source–detector distance and translations, remained fixed in the current studies. We chose orbits that are continuous functions of the rotation angle $\theta $ such that the gantry tilt $\varphi $ is defined by the rotation angle.

A simple parameterization involves *periodic basis functions* using constant, sine, and cosine terms such that

## Eq. (8)

$$\varphi (\theta )=\sum _{i=1}^{K}{\mathrm{\Omega}}_{i}{b}_{i}(\theta )\phantom{\rule[-0.0ex]{1em}{0.0ex}}{b}_{1}(\theta )=1,\phantom{\rule[-0.0ex]{1em}{0.0ex}}{b}_{2}(\theta )=\mathrm{sin}\text{\hspace{0.17em}}\theta ,\phantom{\rule[-0.0ex]{1em}{0.0ex}}{b}_{3}(\theta )=\mathrm{cos}\text{\hspace{0.17em}}\theta ,\phantom{\rule[-0.0ex]{1em}{0.0ex}}{b}_{4}(\theta )=\mathrm{sin}\text{\hspace{0.17em}}2\theta ,\dots ,$$A second parameterization of the source–detector trajectory uses $B$-*spline basis functions*, where the individual parameters ${\mathrm{\Omega}}_{i}$ define a limited set of knot locations. Each knot is fixed to a single rotation angle with equal spacing throughout the orbit and is allowed to vary in tilt angle. This parameterization may more easily admit nonperiodic designs while maintaining relatively low dimensionality. The orbit may be similarly constrained to those feasible with a given C-arm gantry and to avoid table collision. Using either the periodic or B-spline basis functions as presented here not only reduces the dimensionality of the parameterization but also imposes the additional constraint that the trajectory be smoothly changing, which is beneficial for practical implementation of a task-driven orbit from a mechanical point of view.

## 2.5.

### Imaging Performance Prediction and Anatomical Modeling

For prospective design of source–detector trajectories, the imaging performance must be estimated for various orbits. Although exhaustive simulation of projection data, reconstruction, and assessment is possible, it is more practical to estimate the imaging properties of the reconstruction directly. Previous work derived closed-form approximations for local spatial resolution^{35} and noise^{37} in penalized-likelihood reconstructions of the type in Eq. (7). Specifically, the local impulse response, ${l}_{j}$, and local covariance, ${c}_{j}$, may be approximated as

## Eq. (9)

$${l}_{j}\approx {[\mathbf{A}{(\mathrm{\Omega})}^{T}\mathbf{D}\{\overline{y}\}\mathbf{A}(\mathrm{\Omega})+\beta \mathbf{R}]}^{-1}\mathbf{A}{(\mathrm{\Omega})}^{T}\mathbf{D}\{\overline{y}\}\mathbf{A}(\mathrm{\Omega}){e}_{j},$$## Eq. (10)

$${c}_{j}\approx {[\mathbf{A}{(\mathrm{\Omega})}^{T}\mathbf{D}\{\overline{y}\}\mathbf{A}(\mathrm{\Omega})+\beta \mathbf{R}]}^{-1}[\mathbf{A}{(\mathrm{\Omega})}^{T}\mathbf{D}\{\overline{y}\}\mathbf{A}(\mathrm{\Omega})]{[\mathbf{A}{(\mathrm{\Omega})}^{T}\mathbf{D}\{\overline{y}\}\mathbf{A}(\mathrm{\Omega})+\beta \mathbf{R}]}^{-1}{e}_{j},$$^{38}extended these predictors to nonideal detectors and validated the predictions in flat-panel CBCT reconstructions.

Thus, we may predict local noise and resolution properties prospectively, given knowledge of ${\mu}_{\text{prior}}$ and a system model through which $y$ may be simulated. As previously discussed, preoperative CT (or CBCT acquired earlier in the procedure) can provide an anatomical model for prediction, recognizing that the model may be mismatched to the intraoperative data for a number of reasons. First is registration of the previous image to the current measurements. For many anatomical sites (e.g., in neurointerventions), a rigid registration may be sufficient, allowing the designed trajectory to be transformed into the intraoperative patient coordinates using, for example, 3-D to two-dimensional (2-D) registration, as shown in Ref. 39. In this work, we will presume that an accurate registration is achieved and will focus on the subsequent improvements in imaging performance gained from a task-driven design of the orbit.

Other mismatches in patient anatomy include the delivery of hardware and changes in anatomy that might be found in the intraoperative data. Again, the workflow presented in Fig. 1 offers a means to model implanted hardware. As the preoperative scan is often used for planning, e.g., to determine the size and location of an implant, it is relatively straightforward to include these attenuation changes in a modified anatomical model. Modeling of significant attenuation changes like metal implants is important as these changes have a significant impact on the statistics of the data (i.e., the diagonal weighting $\mathbf{D}\{\overline{y}\}$ in the predictors). In contrast, soft-tissue differences like hemorrhage (an important complication one would like to detect) have a relatively small effect on noise in the projection data. This suggests that those more subtle changes do not need to be modeled explicitly for the proposed trajectory design.

## 2.6.

### Approximate Predictors and Practical Implementation

Although the predictors in Eqs. (9) and (10) describe the basic imaging performance metrics required for evaluation of detectability index, these expressions contain large matrix inverses that challenge efficient computation of the local spatial resolution and noise for such a large optimization space. Previous work used local Fourier approximation^{40}^{,}^{41} to yield approximate forms for the local ${\mathrm{MTF}}_{j}$ and local ${\mathrm{NPS}}_{j}$ as

## Eq. (11)

$${\mathrm{MTF}}_{j}=\mathcal{F}\{{l}_{j}\}\approx \frac{\mathcal{F}\{\mathbf{A}{(\mathrm{\Omega})}^{T}\mathbf{D}\{\overline{y}\}\mathbf{A}(\mathrm{\Omega}){e}_{j}\}}{\mathcal{F}\{\mathbf{A}{(\mathrm{\Omega})}^{T}\mathbf{D}\{\overline{y}\}\mathbf{A}(\mathrm{\Omega}){e}_{j}+\beta \mathbf{R}{e}_{j}\}},$$## Eq. (12)

$${\mathrm{NPS}}_{j}=\mathcal{F}\{{c}_{j}\}\approx \frac{\mathcal{F}\{\mathbf{A}{(\mathrm{\Omega})}^{T}\mathbf{D}\{\overline{y}\}\mathbf{A}(\mathrm{\Omega}){e}_{j}\}}{{|\mathcal{F}\{\mathbf{A}{(\mathrm{\Omega})}^{T}\mathbf{D}\{\overline{y}\}\mathbf{A}(\mathrm{\Omega}){e}_{j}+\beta \mathbf{R}{e}_{j}\}|}^{2}},$$## Eq. (13)

$$\mathcal{F}\{\mathbf{A}{(\mathrm{\Omega})}^{T}\mathbf{D}\{w\}\mathbf{A}(\mathrm{\Omega}){e}_{j}\}={\mathbf{L}}_{j}w,$$*precomputation approach*becomes tractable without additional modification.

For a large number of potential views, memory limitations may challenge the precomputation approach. In the spirit of Refs. 43 and 44, one can also recognize that we may replace ${\mathbf{L}}_{j}$ with an analytic form. That is, each column of ${\mathbf{L}}_{j}$ is the Fourier transform of the backprojection of a single projection view consisting of a single point projection. The backprojection of a point is a line in 3-D space through the original location $j$ connecting both the source and detector. The Fourier transform of this line is a Fourier plane centered at the origin at the same angle as the line. However, for a discrete system model with finite-sized detector elements, voxels, etc., this Fourier plane will not be infinitely thin. Whereas Schmitt derived closed-form analytic expressions that specify the form and profile through this plane, we opted to use a plane whose profile is Gaussian with a width specified by fitting an empirical calculation of the Fourier-domain projection-backprojection. This results in a very fast analytic form for ${\mathbf{L}}_{j}$ that may be computed, given the location $j$ and the coordinates of the x-ray source. This *on-the-fly computation approach* is another efficient alternative for calculating Eqs. (11) and (12). The difference between the precomputation and on-the-fly approaches is largely computational (without significant difference in noise and resolution estimates) and represents a classic computing trade-off between storage and speed. As both approaches are potentially useful, both are presented in studies below.

## 3.

## Experimental Methods

Three experiments were conducted to investigate various aspects of the proposed task-driven trajectory design process, summarized in Fig. 4. Each experiment is detailed in the following subsections.

Two simulation configurations were used in these studies. Studies that investigated basic trajectory design behavior used a compact system geometry to better illustrate dependencies on location and task, specifically 700-mm source-to-detector distance and 350-mm source-to-axis distance. A flat-panel detector with $560\times 1000\text{\hspace{0.17em}\hspace{0.17em}}\text{pixels}$ at 1-mm pitch was simulated. This geometry is referred to as the *compact geometry*. A second system geometry emulated a $C$-*arm geometry* with 1200-mm source-to-detector distance, 800-mm source-to-axis distance, and a $960\times 1240\text{\hspace{0.17em}\hspace{0.17em}}\text{pixel}$ detector at 0.308-mm pitch.

## 3.1.

### Location Dependence: A Sphere in a Cylinder

The first experiment used a simple object to illustrate location dependence. Specifically, a 20-cm diameter cylinder was simulated with 1-mm voxels and an attenuation of $0.05\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{mm}}^{-1}$ [Fig. 4(a)]. A relatively high attenuation was used to exaggerate the location-dependent effects. Two 3-mm spheres with $0.03\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{mm}}^{-1}$ contrast relative to background were added to the cylinder, centered in the axial slice, at the central slice and at 15 cm below the central slice. This study used the compact geometry, and bare-beam fluence was set at ${10}^{5}\text{\hspace{0.17em}\hspace{0.17em}}\text{photons}$ per pixel.

Trajectory design was conducted using the single-location objective in Eq. (2) using nine periodic basis functions given in Eq. (8), constrained to tilts in the range of $\text{\hspace{0.17em}\hspace{0.17em}}\varphi =-50\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}$ to 50 deg, and with rotation angles $\theta =1\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}$ to 360 deg with 1 deg increments. The task function in the objective corresponded to the 3-mm spherical stimulus [Fig. 4(d)] and optimization was performed for each stimulus location individually. For computation of detectability, the on-the-fly method was applied with statistical weights ($w$, calculated as described in Sec. 2.6) sampled over 110 equally spaced rotation angles and 51 equally spaced tilt angles. To estimate the solution of Eq. (2), the CMA-ES optimization was applied using a population size of 40 without restarts due to the simplicity of the search space.

Penalized-likelihood reconstruction using the optimized trajectory found for each sphere location (central slice and 15 cm below the central slice) was performed using dynamically relaxed ordered subsets with the number of subsets changing every five iterations through the sequence {54, 24, 12, 6, 4, 2, 1} for a total of 50 iterations. This schedule was chosen to accelerate the convergence of the simple object. Quadratic regularization with a regularization strength of $\beta ={10}^{6}$ was applied using 1-mm isotropic voxels on a $240\times 240\times 500$ grid. This choice of regularization strength was empirically selected based on visual assessment of noise and resolution in the images.

## 3.2.

### Task Dependence: Line Pairs in a Cylinder

The second study used the same 20-cm cylinder as the first experiment—in this case involving a line-pair stimulus placed in the center of the central slice. Specifically, the line pairs were a $20\times 20\times 20\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{mm}}^{3}$ cube with $2\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{lp}/\mathrm{mm}$ and a contrast of $0.015\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{mm}}^{-1}$ relative to background. The cube was rotated 20 deg around the x-axis to angulate the features relative to the coordinate axes [Fig. 4(b)].

Trajectory design was performed using the single-location objective, nine periodic bases, and the same angular sampling, constraints, geometry, and bare-beam fluence as in the first experiment. Two task functions were examined: the 3-mm sphere detection task from the first experiment and a line-pair discrimination task. The discrimination task was defined as the difference between the true rotated line-pair stimulus and a stimulus consisting of an identically rotated solid cube of the same dimensions and attenuation, i.e., a binary hypothesis of line pairs “present” or “absent” [Fig. 4(e)]. As in the first experiment, on-the-fly detectability computations were adopted and the same CMA-ES parameters were used. Penalized-likelihood ordered-subset reconstructions were computed as in the first experiment, using the {54, 24, 12, 6, 4, 2, 1} schedule, 50 iterations to reach convergence, $\beta ={10}^{6}$, and 1-mm isotropic voxels on a $240\times 240\times 500$ grid.

## 3.3.

### Multiple Locations: Elliptical Cylinder with Needle

The third experiment examined task-driven imaging in a situation for which the precise task location was unknown. We therefore chose several potential task locations and solved for a single trajectory that maximized a multilocation objective. A digital elliptical cylinder phantom was modeled with major and minor axis diameters of 25 and 17.5 cm, respectively, and with a height of 25 cm. The elliptical cylinder had a 12-mm outer shell with attenuation coefficient of $0.04\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{mm}}^{-1}$ and the core was filled with low-contrast spheres of diameter 80 mm and attenuation coefficients ranging from 0.0175 to $0.0225\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{mm}}^{-1}$. A high-contrast cylinder (attenuation coefficient $0.2\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{mm}}^{-1}$) with a diameter of 5 mm and a length of 75 mm was added at a 10-deg angle in the central coronal plane to simulate a needle entering the body [shown in Fig. 4(c)]. Nine 8-mm spheres with $0.08\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{mm}}^{-1}$ contrast were added along the length and tip of the cylinder and their nine locations were used for multilocation optimization. We used the C-arm geometry with bare-beam fluence lowered to ${10}^{4}\text{\hspace{0.17em}\hspace{0.17em}}\text{photons}$ per pixel to mimic scanning at lower dosage.

Trajectory design was first performed for each task location using the single-location objective in Eq. (2), followed by the three multilocation objectives in Eqs. (3)–(5). The B-spline basis functions were used with eight equally spaced knots, and the trajectories were constrained to tilt angles in the range of $\varphi =-30\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}$ to 30 deg and rotation angles of $\theta =1\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}$ to 360 deg. The task function at all nine locations consisted of midfrequency content [Fig. 4(f)] to discriminate the presence of a single object or two separate objects (needle and sphere), as described in the ICRU Report 54.^{45} The precomputation approach was used with $w$ sampled over 72 equally spaced rotation angles and 13 equally spaced tilt angles to compute detectability. The CMA-ES optimization algorithm was used to estimate $\hat{\mathrm{\Omega}}$ using a population size of 200 with six restarts and random initialization to ensure that the optimal detectability values were achieved due to the more complicated search space compared to the previous two experiments.

Quadratic penalized-likelihood reconstructions for the nominal circular orbit, nine independent single-location optimizations, and three variations of multilocation optimization were performed using 10 ordered subsets over 200 iterations. The total number of iterations was increased and the number of subsets was decreased from previous experiments to achieve stability in convergence of the structurally more complicated object. Regularization strength $\beta $ was set at ${10}^{5.5}$ to decrease noise in the highly attenuating object, and 0.5-mm isotropic voxels were used in a $512\times 512\times 512$ grid.

To illustrate the convergence of the CMA-ES algorithm, the optimization for a single sphere was performed as described above, using a circular orbit for initialization as opposed to random initialization. This allowed comparison of the detectability index and the reconstructed images associated with the initial circular orbit, an intermediate (suboptimal) orbit, and the resulting optimal orbit.

## 4.

## Results

## 4.1.

### Location Dependence: A Sphere in a Cylinder

Results from the first experiment are summarized in Fig. 5. Task-driven source trajectories are shown for the cylindrical object and two optimizations: maximum detectability of the spherical stimulus at locations #1 [Fig. 5(b)] and #2 [Fig. 5(c)]. Source orbits are shown for each case in magenta as well as a sampled sphere of all possible source locations (blue dots). Note that the designed orbits suggest that a simple equatorial source trajectory around each stimulus is optimal for detection of that stimulus. Reconstructions using the two orbits confirm that conspicuity of each sphere is maximized for such an equatorial orbit, with reduced detectability for the stimulus in the location that is not optimized. That is, the stimulus at location #2 is more difficult to detect when using the orbit designed to maximize detectability at location #1, and vice versa.

This simple, idealized simulation offers basic intuition to the optimization by considering the fluence through each stimulus location for all potential views. This fluence is equal to the statistical weighting—$w$ in Eq. (13)—and we see that the orbit (dashed magenta line) in each case maximizes fluence through the target location, thereby maximizing data fidelity by selecting projection views with the shortest path length through the object. This experiment illustrates the importance of location in task-driven designs and shows that using a single-location objective can maximize task performance at a given location but may do so at the cost of decreased performance at other locations.

## 4.2.

### Task Dependence: Line Pairs in a Cylinder

The results for a task with strong directional dependence (a line-pair stimulus within a cylinder) are summarized in Fig. 6. Single-location task-driven designs are performed for two tasks: the sphere detection task from Sec. 3.1 [Fig. 6(b)] and a discrimination task corresponding to the frequency content and angulation of the line-pair stimulus [Fig. 6(c)]. For the task function corresponding to the line-pair stimulus, the task-driven orbit (orbit A) is tilted to match the angulation of the line-pair stimulus. Example reconstructions illustrate that the task-driven design outperforms the nominal orbit (orbit S), thus showing increased noise between the line pairs for the trajectory optimized for the spherical task [Fig. 6(b)] and improved discrimination of line pairs for the orbit optimized according to the line-pair task [Fig. 6(c)].

## 4.3.

### Multiple Locations: Elliptical Cylinder with Needle

Figure 7 summarizes the task-driven trajectories resulting from optimization with respect to a single-location objective (in which ${d}^{\prime}$ is optimized for each location independently), and for the three variations of a multilocation objective [maxi-mean, maxi-median, and maxi-min shown in Eqs. (3)–(5), in which ${d}^{\prime}$ is optimized over all locations simultaneously]. The fluence through the stimulus location for all potential views is shown for each of the nine task locations with the optimal orbits shown in magenta. The orbits resulting from a single-location objective are shown in the left-most column, and the single orbit resulting from a multilocation objective is shown in the three right columns for the maxi-mean, maxi-median, and maxi-min objectives for all nine locations. The single-location objective shows a different orbit for each location, and these can be considered the best possible trajectory for each location. Note that at some of the locations the tilt angle differs for the starting and ending views (i.e., locations 4, 5, and 7). The 2-D B-spline basis function representing the source trajectory is not constrained to be equal at the start and end vertices, allowing a discontinuity. In these cases, an orbit with a discontinuity produces a higher ${d}^{\prime}$ value at the task location. We may theorize that allowing discontinuities at any point in the orbit (not only at start and end vertices) might increase ${d}^{\prime}$; however, doing so loses the advantage of a low-dimensional parametrization and increases the difficulty of implementation on a robotic C-arm.

When comparing the three multilocation objectives (in which a single trajectory is generated for all stimulus locations), each objective function yields a different result, as expected. The trajectory generated when using a maxi-median objective is dictated by the median ${d}^{\prime}$ value, which in this case corresponds to location 3 for the optimal trajectory. Note that the orbit is not the same as that for location 3 achieved when using a single-location objective. This is because further changes to the orbit would cause all ${d}^{\prime}$ values to shift, which would generate a new median value that does not maximize the objective function. The maxi-min objective is similarly driven by a single location, i.e., location 1. Location 1 is at the tip of the needle and is subsequently obscured by the angulated needle to some extent at all possible views. This results in the lowest fluence, compared to all other locations for all possible views. As the ${d}^{\prime}$ value at location 1 is significantly smaller than the values at locations 2 to 9 due to the lower fluence, it remains the minimum ${d}^{\prime}$ value for all trajectories examined during optimization. This results in an equivalent optimization as when using a single-location objective at location 1, and the resulting orbit is almost identical to that case. When using a maxi-mean objective, all locations jointly drive the orbit, and no single location skews the result.

Figure 8(a) shows detectability index for the nine task locations corresponding to a circular trajectory, using a single-location objective, and using the three multilocation objective functions. The single-location objective is seen to provide the highest ${d}^{\prime}$ value and represents the upper limit of improvement in ${d}^{\prime}$. The respective advantage of each multilocation objective function is evident, with maxi-mean having the highest mean ${d}^{\prime}$ value at 3.1, maxi-median having the highest median ${d}^{\prime}$ value at 3.7, and maxi-min having the highest minimum ${d}^{\prime}$ value at 1.3. Figure 8(b) shows the percentage increase (or decrease) in ${d}^{\prime}$ compared to the location-matched ${d}^{\prime}$ value for a circular orbit. For a single-location objective, all ${d}^{\prime}$ values increase, ranging from 5.5% at location 8 to 121.1% at location 1. In contrast, for all three multilocation objectives there are locations that exhibit a decrease in detectability. For maxi-mean, the change in ${d}^{\prime}$ ranges from $-8.4\%$ at location 4 to $+37.7\%$ at location 1. For maxi-median, the range is from $-46.1\%$ at location 4 to $+45.4\%$ at location 9. For maxi-min, the range is from $-35.4\%$ at location 3 to $+120.9\%$ at location 1.

The reconstructed image in the region around each stimulus is shown for each orbit in Fig. 8(c), as well as for a nominal circular orbit for comparison (top row). Either the axial, coronal, or sagittal plane is chosen at each location for visual representation of improvements, although the reconstructed plane could be arbitrary. Improvements are particularly apparent for locations 2, 7, and 9.

Figure 9 demonstrates convergence of the CMA-ES algorithm for the single-location objective performed for the task defined at location 7. Images reconstructed for a circular orbit, suboptimal intermediate orbit, and optimal task-driven orbit are shown in Figs. 9(a)–9(c). The ${d}^{\prime}$ values for each are 2.85, 3.35, and 4.08, respectively. Visualization of the spherical stimulus alongside the highly attenuating needle is progressively improved as the algorithm converges on the optimal orbit, consistent with the increase in ${d}^{\prime}$. A plot of detectability index versus iteration is shown in Fig. 9(d), which shows a fairly smooth convergence for the CMA-ES optimizer when initialized with the circular orbit. The ${d}^{\prime}$ values corresponding to the three orbits (A–B–C) in Figs. 9(a)–9(c) are indicated in the plot.

## 5.

## Discussion

In this work, we have presented a framework for task-driven trajectory design for advanced CBCT imaging systems that leverages knowledge of both the patient anatomy and imaging task and exercises the motion capabilities of motorized, multi-axis CBCT systems to maximize imaging performance. The task-driven approach provides a strategy to overcome traditional imaging limitations via orbital flexibility, e.g., challenges associated with highly attenuating anatomy or implants that can be mitigated through intelligent data collection, selecting the best projection views to accomplish a particular task.

The importance of both the location and spatial-frequency dependence of the imaging task has been investigated. Each of these elements can contribute significantly to which projection views carry the greatest information for a specific task. In general, the task-based approach balances the data fidelity of a view (i.e., noise) with the signal content (frequency response) that each view provides toward the imaging task. Because the trajectory design for a single stimulus location in the image volume can optimize performance at that location to the detriment of performance elsewhere, a number of multilocation design objectives are examined to obtain optimal performance over regions of interest.

The proposed framework is general with respect to various options for parameterization of the source–detector trajectory, including constraints based on system geometry, degrees of freedom, and collision avoidance—permitting application to robotic C-arms and other CBCT systems with additional flexibility to gantry tilt, etc. This includes short-scan geometries, which are popular in CBCT imaging. A short-scan can be similarly optimized as the above-described full-scans with additional parameters such as start and end rotation angles included in the optimization to seek the views that are best for the object and task. Moreover, it is straightforward to further generalize these methods with additional geometric degrees of freedom (e.g., axial translation and variation in source–detector distance), acquisition parameters (e.g., tube current modulation^{28} and fluence-field modulation^{29}), and reconstruction parameters (e.g., regularization strength, $\beta $^{26}). Such factors certainly carry interdependency in a full multivariate optimization over all parameters—well beyond the scope of the current work, which focuses specifically on optimization of the source–detector orbit (holding other parameters fixed).

We note that Fourier approximation of both the MTF and NPS is one potential source of error in orbital design. Specifically, this presumes that the term ${\mathbf{A}}^{T}\mathbf{DA}{e}_{j}$ is locally shift-invariant. For heterogeneous objects, the diagonal weighting is variable with location. However, we note that the data-dependence of this term falls between projection and backprojection operations, which have the effect of strong smoothing even for highly nonuniform weights. The impact of this smoothing is that the Fourier approximation tends to hold well even for very nonuniform objects; however, there is the potential for breakdown of the approximation when shift variance is pushed to the same scale as the spatial resolution of the system.

Consideration should also be given to challenges of accurate CBCT image reconstruction. Physical factors that are not accounted for include x-ray scatter, beam hardening, image lag, detector glare, and patient motion, which may produce artifacts that confound visualization of the task. Recent work has investigated the optimization of orbits based on scatter and scatter-to-primary ratio for imaging the weight-bearing spine using a robotic CBCT system, indicating that these effects are small and can be included in the overall optimization.^{46} Effects that occur from the sampling pattern itself, such as streakiness and cone-beam artifacts, are included in the model of the local MTF and are therefore represented in the task-driven image in such a way as to improve the detectability of the task (although they may not be eliminated from the image).

Accurate geometric calibration of noncircular orbits is also required for accurate image reconstruction. The motion of the C-arm gantry as it moves through complex task-driven orbits may result in gantry wobble. This can be accounted for using a calibration technique like the “self-calibration,” as described in Ref. 47, that performs the geometric calibration subsequent to the scan using the set of acquired projections.

Because the design objectives are generally nonconvex and require many predictions of image quality, different strategies for efficient evaluation and solution of the design objective are presented. Image quality predictors in the current work focus on quadratically penalized model-based image reconstruction; however, future work aims to extend this methodology to nonquadratic regularization.^{44}^{,}^{48}

Task-driven trajectory design has the potential for application in a number of imaging scenarios. This work has focused on theoretical underpinnings of the framework and has illustrated the method in a series of experiments ranging from simple to complex. Although the phantoms used in these scenarios had somewhat uniform backgrounds, the method is designed to incorporate complex anatomic variations in the optimization via the patient model. In fact, the complexity of the model drives the selection of optimal trajectory. This is demonstrated in the Part II companion paper^{49} titled “Task-driven source–detector trajectories in cone-beam computed tomography: II. Application to neuroradiology,” where we investigate application to imaging scenarios in neuroradiology with complex surrounding anatomy. We expect that the methodology also has potentially broad application in other interventional imaging scenarios (e.g., orthopedic procedures). The above-detailed theoretical foundations suggest a new paradigm for interventional imaging wherein preoperative information is included explicitly within a rigorous definition of the imaging task to prospectively drive customized data acquisition that maximizes performance. The framework is an important first step in realizing advanced CBCT capabilities and more fully leveraging the wealth of information available in interventional imaging scenarios.

## Acknowledgments

This research was supported by the National Institutes of Health Grant Nos. R01-EB-017226 and U01-EB-018758, academic-industry partnership with Siemens Healthineers (AX Division, Forcheim, Germany), and a Johns Hopkins University Catalyst Grant.