2.1

Data acquisition

Our photon counting micro-CT system uses a SANTIS 1604 detector (DECTRIS Ltd.). The detector is constructed with a 1 mm thick CdTe sensor, 150-μm isotropic pixels (1035x257), and four independent energy thresholds, here set to 20, 25, 34, and 50 keV to bracket the K-edge of iodine. The projection data set was acquired with a G297 X-ray tube (Varian Medical Systems; fs = 0.3/0.8 mm; tungsten rotating anode; filtration: 0.1 mm Cu; 80 kVp, 2.5 mA, 90 ms exposure / projection). The source-to-detector and source-to-object distances were 821 mm and 674 mm, respectively. To minimize ring artifacts in our reconstructions, we scanned using a helical trajectory with 3 rotations (45 seconds/rotation) and 1.25 cm of total translation during scanning, acquiring a total of 1500 projections per threshold. This acquisition protocol was repeated 4 times, consecutively, yielding a total of 6000 projections over 540 seconds. The absorbed radiation dose associated with imaging was ~570 mGy. To image perfusion, we injected 0.5 ml of ISOVUE-370 (Bracco Diagnostic Inc.) over 10 secs via a tail vein catheter using a computer-controlled injector (Nexus 3000; Chemyx Inc.), starting 22.5 seconds after beginning the scan.

The animal scan conducted for this work followed protocols approved by the Duke University Institutional Animal Care and Use Committee. The animal model was a syngeneic transplant oral cavity squamous cell carcinoma model (MOC2)⁶ initiated in a C57BL/6J mouse by implantation of 30,000 cells into the buccal mucosa (median time to tumor onset: 2.2 weeks). The model exhibits metastasis to the cervical lymph nodes, similar to the natural history of human HNSCC. Imaging was performed when caliper measurements of the primary tumor fell in the range of 50-75 mm³.

2.2

Reconstruction

To handle the channels of photon-counting perfusion data (3D volumes with n_x voxels, n_t time points, and n_e energy thresholds), we organize the reconstructions into a 3D array:

where, e.g., x_t=1,e=2 denotes a column vector of the reconstruction at perfusion time point 1 and energy threshold 2. Each energy threshold slice corresponds to a 2D matrix. Following from this notation (array, A; matrix, A; vector, a; scalar, a; index, a), algebraic reconstruction minimizes the reconstruction error relative to log-transformed projection data, y, under an L2 norm penalty applied separately to each channel:

R is the system-specific projection matrix, which is universal for all channels when working with co-registered photon counting projection data. The weight array Q is constructed from weight matrices, Q_t,e, the product of diagonal matrices indicating the target energy threshold and perfusion phase.

We use multi-channel regularization to improve reconstruction performance. Specifically, we perform iterative reconstruction using the split Bregman optimization framework,¹⁰ and we employ an extended version of our RSKR regularizer to enforce gradient sparsity and data consistency between channels.⁷ Our iterative reconstruction is governed by the following objective function:

X_nxnt×nedenotes the unfolding of the array X into a matrix with n_xn_t rows and n_e columns (X_nxne×nt: n_xn_e rows, n_t columns). These unfoldings and their corresponding regularization terms enable penalization of the global spectral rank (||.||_*, nuclear norm; λ,₁), the spectrally joint intensity gradient sparsity (BTV: bilateral total variation⁷; λ_S), and the patch-wise temporal rank (λ_*2; n_p total patches; P_p : patch extraction operator for patch p). Within the split Bregman framework, it is feasible to split each regularization term into a separate regularization sub-step (each with its own set of auxiliary variables, D and V) to incrementally solve this objective function. Practically, however, the computational cost of the data fidelity updates (Fig. 1, step 6; here, evaluated on 48 volumes) and the often slow rate of convergence of shrinkage methods (singular value thresholding to reduce the nuclear norm¹¹) make this approach intractable.

Fig. 1.

The split Bregman method applied for iterative reconstruction of photon-counting CT perfusion data. Starting from the left column, an Initialization procedure estimates the reconstructed data, X, and the regularization parameters, G (here α = 0.007, h₀ = 1.5, γ = 0.5, ν = 0.7). Outer Bregman iterations approximately solve the Objective function through regularization (step 4, D), residual (5, V), and data fidelity (6, X) updates. More details on the inner Bregman iterations for regularization are provided in the text (steps 4a-4i). Six iterations of the Bi-CGSTAB(2) solver (algorithm 3.1 in ¹⁴) are used to update X at each time point and energy during Initialization (step 1) and for each data fidelity update (step 6). Additional notation: μ_water,e, the expected attenuation of water at energy threshold e; ∇, compute high-pass only component of the 3D Haar wavelet transform.; i, singular value, vector index (n_e total); x_:,e + v_:,e time indices assumed equivalent for noise estimation.

To reduce the number of global Bregman iterations (data fidelity update steps) required to achieve convergence (here, 3), we address all of the regularization terms in a single regularization step (step 4). This one step is itself solved by a series of inner Bregman iterations (3-4), which perform image-domain denoising with our RSKR algorithm (steps 4a-4i). RSKR jointly minimizes the spectral rank (λ_*1 term) and the spectral gradient sparsity (λ_S term) by performing joint bilateral filtration on singular vectors computed along the energy dimension (step 4b). As in our previous implementations of RSKR, we use the median absolute deviation (MAD¹²) to estimate the noise level in each channel. These noise estimates are used to weight the singular value decomposition such that high singular values correspond to singular vectors with proportionately high signal-to-noise ratios (steps 4a, 4b). As in past work, a similar noise estimation strategy is used to calibrate the joint bilateral filter parameters (jBF⁷; step 4e), reducing the number of free parameters which must be determined by the user to achieve robust performance.

New to this work, we have improved our implementation of BF to use local noise estimates (vs. one global noise estimate per energy) to deal with spatially variant noise levels in CT data. Our GPU-based implementation of jBF now computes the MAD of image gradients in 32x32 patches with a stride of 16 voxels in axial planes. Overlapping estimates are averaged to prevent sharp transitions in denoising performance, while the noise level is estimated in axial planes, rather than volumetrically, to improve computational efficiency.

Also new to this work, we perform patch-based singular value thresholding (pSVT⁹) along the time dimension of the singular vectors. Notably, pSVT is applied to the singular values of a second singular value decomposition performed on 3³ × n_t patch matrices extracted from L (step 4f). Given a vector of these singular values, ε, computed from patch matrix p, thresholding is applied to the singular values (indexed by i; evaluated separately per energy) as follows:

These thresholded singular values, ε′, are then used to reverse the SVD for each patch matrix, and the central voxel for each time point is stored as the regularized output.

Previously, we have used kernel smoothing to overcome inconsistent denoising performance at high contrast edges when denoising with jBF.¹³ Now, we apply pSVT (v = 0.7). pSVT similarly improves image consistency at edges, but also improves consistency along the time dimension and may better preserve spatial resolution. pSVT addresses the final cost term of the Objective function (Eq. 3; Fig. 1; λ_*2). Like with jBF, MAD-based noise estimates are used to calibrate the threshold values, τ_i,e, used during each inner Bregman iteration. Specifically, the SVD is performed on zero-mean Gaussian noise patch matrices with noise standard deviations matching those measured globally in the corresponding energy channel. Here, the average singular values from 125 noise realizations were used to set the threshold values per energy.

2.3

Spectral processing and analysis

We compare three sets of reconstruction results: analytical reconstructions of each time point and energy (WFBP algorithm¹⁵), iterative reconstruction without pSVT (Fig. 1, steps 4f and 4g skipped), and iterative reconstruction with pSVT (Fig. 1, all steps). Notably, the reference iterative reconstruction (no pSVT) is similar to spectral reconstruction results we have quantitatively validated in previous work⁷ (with the exception of localized noise estimation now performed during jBF). Following reconstruction, we perform non-negative material decomposition using sensitivity matrix inversion and a sub-space projection with iodine (red), photoelectric (PE, green), and Compton scattering (CS, gray) basis functions.³ We characterize the convergence of our iterative reconstructions (Fig. 1, step 6) by measuring the relative error after the final outer Bregman iteration for each perfusion phase and energy threshold (generically, ||Ax − b||₂/||b||₂; Ax: variable terms; b: constant terms). We compare our reconstruction results through standard deviation measurements, residual images, material decomposition errors in reference solutions, and through modulation transfer function measurements.

2.4

Computation

All reconstruction and denoising operations were performed with our custom multi-channel GPU-based reconstruction toolkit.¹⁶ For the mouse data set, a total of 6000 projections were acquired over 12, 360° rotations. The temporal basis functions (rect functions, Q_t,e) assigned equal weight to every projection in each subset of 500 projections, allowing reconstruction of 12 non-overlapping perfusion phases spanning 45 seconds each. Combined, 12 perfusion time points at 4 four energy thresholds were reconstructed (48 total volumes) with a volume size of 360x360x384 voxels and 125-micron, isotropic voxels. For the analytical reconstruction results, a ramp filter was used. For algebraic reconstruction (Fig. 1, steps 1 and 6) six iterations of the Bi-CGSTAB(2) solver were used.¹⁴ Reconstructions were performed on an Ubuntu Linux workstation (v18.04) with four NVIDIA RTX8000 GPUs, 256 GB of system RAM, and two Intel Xeon W-2295 CPUs. Iterative reconstruction of all 48 volumes took 2.9 hours without pSVT (3 outer, 3 inner Bregman iterations) and 11 hours with pSVT (3 outer, 4 inner Bregman iterations).