2.1.

Four-Dimensional Digital Subtraction Angiography

4-D DSA is made possible by combining the 3-D DSA volume and the 2-D rotational projections containing temporal information using the 4-D DSA algorithm further discussed in Fig. 2 and described by Eqs. (1) and (2).

The projection sampling of a C-Arm CT system is shown in Eq. (1), where $T$ is the sampling period. The constraining image, $i_C$, is first generated by applying a threshold to the 3-D DSA volume. A projection containing temporal information, $p_t$, is convolved using the operator “*” with a 2-D square blurring kernel, $k_B$, of uniform intensity. The forward projection of the constraining image $p_C$ is also blurred by the kernel $k_B$. The convolution of projections acts to increase SNR (Ref. 7) of the projection image. These projection data are then divided element by element and backprojected, $f_{\mathrm{BP}}$, into 3-D space at the view angle $\theta$ to create a weighting volume $i_W$. The $f_{\mathrm{BP}}$ operator backprojects a projection value at detector coordinates $u$ and $v$ to a 3-D coordinate $x$, $y$, and $z$ in 3-D space given the view angle $\theta$. The backprojection operation does not perform any filtering or weighting of the projection data. This work implemented the p-matrix voxel driven backprojection algorithm as described in Scherl et al.,¹¹ Wiesent et al.,¹² and Galigekere et al.¹³ The Hadamard product “$°$” of the weighting volume, $i_W$, and the constraining image, $i_C$, generates raw 4-D time frame, $i_{4\mathrm{D}}$, data, which can be further processed. The postprocessing includes corrections for vessel overlap in a projection, which causes a reconstruction artifact, higher than actual values during overlap periods, in the TACs of 4-D DSA.

The 2-D projections record the contrast enhanced dynamics of the vascular network under study following the single injection of contrast near the start of the acquisition. The 2-D projections are used to generate the 3-D DSA volume through conventional techniques.^12,14 The 3-D DSA volume is then sparsified using a simple thresholding approach to generate a 3-D volume for use as a volumetric constraining image, $i_C$. Constraining volume generation using 3-D DSA has the effect of removing background noise, thus allowing the viewing of the vascular network. The process is similar to that of window and leveling techniques currently used to view current 3-D reconstruction volumes.

The angle $\theta (t)$ represents the scan acquisition trajectory view angle at the acquisition time frame $t$. There is a direct relationship between the projection angle $\theta$ and the time index $t$; therefore, the terms will be used interchangeably. For forward projections, $p_C[\theta (t)]$ of the constraining image are taken at the same angle as the temporal projection $p_t(t)$. The constraining image, $i_C$, is a static image with no temporal variation, where $p_C(t)$ has no meaning other than to simplify the equations and shall be equivalent to $p_C[\theta (t)]$.

The equations, descriptive geometry, and reconstruction algorithm for 3-D DSA are discussed in numerous sources, notably Refs. 12 and 14. The 3-D DSA reconstruction process involves cosine weighting, Parker-weighting,¹⁵ and convolution of the projections with a one-dimensional (1-D) filtering kernel followed by a backprojection into 3-D space. Convolution occurs only along the dimension perpendicular to the axis of rotation and is not a 2-D convolution. The 3-D DSA and 4-D DSA reconstructions used the same 1-D Hamming filter as the 4-D DSA constraining image is derived from the 3-D DSA.

The convolved temporal and constraining projections have a consistent attenuation where vessels occur. However, due to vessel overlap in a given projection, attenuation is increased in overlap areas, which creates artifacts in the per voxel TAC of 4-D DSA when reconstructed. To mitigate this effect, the convolved temporal projections are then normalized by dividing the constraining projections element by element to aid in minimizing artifacts due to vessel overlap to create the normalized 2-D temporal time frame projection $p_{4\mathrm{D}}(t)$ as shown in Eq. (3). Normalization does not completely eliminate overlap as the constraining image is a time average integral through the vessels.

What is ultimately desired out of the division process is a normalized projection in the set of real numbers, $p_{4\mathrm{D}}\in {\mathfrak{R}}^n$ and $p_{4\mathrm{D}}\in [\mathrm{0,1}]$. An ideal normalized projection allows the activation or deactivation of the vessel network in the constraining image when backprojected to create the weighting image. In practice, this is challenging to obtain as the constraining volume is a time average of the reconstructed temporal projections. Thus, the reprojection of the constraining image is a projection of a time averaged volume. Normalization allows the weighting of the image by a factor between zero (off) and one (on). A method found to work is described in Eq. (4). The denominator is increased by a factor that is five percent of the max value of the projection data in the denominator.

While normalization of the projections can aid in decreasing the vessel overlap artifact, algorithms have been created to further reduce the artifact. The single projection (SP) algorithm described in Eq. (2) is the basic algorithm of 4-D DSA. It became clear early on in the investigation of 4-D DSA that this algorithm would require further refinement to account for vessel overlap in the projections. 4-D DSA separation angle (SA) and minimum within a search angle (MSA) algorithms have been developed to further process the 4-D data generated by Eq. (2) temporally on a voxel-by-voxel basis.

The separation angle algorithm is described in Eq. (5) and takes two projections separated by an angle. The square root is used as an approximation to retain the actual values as two projections and thus two weighting volumes are used. A typical value for $\tau$ was one that would yield an angular separation in the acquisition trajectory of 30 deg. The SA is selected to provide a physical projection angle separation while still maintaining acceptable temporal resolution. Ideally, the angle selected would be 90 deg, providing ideal spatial separation; however, selecting this value is at the cost of temporal resolution.

The minimum search angle algorithm as described by Eq. (6) aims to decrease vessel overlap artifacts by searching for the minimum value of the temporal data for a given voxel over a search angle to mitigate the contrast curve intensity spikes that can occur when vessels overlap in a projection. The value $t_{\min }$ is the minimum temporal index over a time range $t^{\prime }\in [t_i\dots t_{i+\tau }]$, where the argument minimum of $i_{4\mathrm{D}}(t^{\prime })$ is at a minimum on a per voxel TAC basis. The case as $t_i$ approaches the last temporal index $t_{i\max }$, the temporal delta $\tau$ is allowed to approach zero, $\tau \to 0$, as the number of remaining time frames decrease. For a given voxel TAC, the time index selected for the time frame is the time frame with the lowest TAC value. This time index thus selects the lowest TAC value over the search range. Thus, $t_{\min }(t_i)$ represents a 3-D array of time indices for each voxel TAC where the value in each voxel TAC is a minimal value over the temporal search range.

While in an ideal mathematical sense MSA is equivalent to SP and SA algorithms when a perfect static object is used, when used in real-world data, this is not the case due to imperfections in scan geometry and data acquisition. However, limiting spatial resolution results for the other algorithms were found to be identical to the results of MSA. The effects of the algorithms were not the focus of this research. Recently, various interpolation algorithms have been developed and investigated to interpolate over areas of vessel overlap in the TAC data. The impact of these interpolation algorithms on resolution were not investigated in this work.

2.2.

Frequency Response

A simplified form of the equation for 4-D DSA is shown in Eq. (7). The division by the reprojection of the constraining volume has been omitted for the current purposes of discussion; the theta notation and time variable have been dropped for the sake of simplicity.

The notation shown shall use lowercase to represent equations in image space with uppercase variables representing the Fourier transformed variables in frequency space. The temporal projection data, $p_t$, is first convolved with the blurring kernel, $k_B$. This operation is a convolution in image space, but is a multiplication in frequency space as described by Eqs. (8), (9), and (10). This is essentially a low-pass filter operation in which frequency content in the projection $P_t$ is limited to the bandwidth of the filter kernel $K_B$. Basic properties of MTF and frequency analysis will be used to discuss the effects of 4-D DSA on the limiting spatial resolution. The equation ¹⁶ using the Fourier transform, $\mathcal{F}$, for finding the MTF given the point spread function (PSF) data is shown in Eq. (11). The solution for the MTF of the system is shown in Eqs. (12) and (13). For the purpose of discussion, the limiting spatial resolution of two data sets in object space, $a$ and $b$, having been transformed into frequency space, $A$ and $B$, shall be defined here as $f_A$ and $f_B$, respectively. The maximum limiting spatial frequency $f_M$ resulting from product of the two frequency responses is defined in Eq. (14). The limiting spatial resolution in MTF analysis is either the point where the MTF reaches the first zero when using a pure mathematical approach or in analysis of real-world data, the first point on the MTF curve that reaches 10% of the value at zero spatial frequency starting from zero spatial frequency. The maximum limiting spatial frequency, $f_M$, of a convolution operation in image space is the minimum of the spatial frequency of the data being convolved, sets $a$ and $b$, as the multiplication in frequency space is dependent on the first zero crossing of either of the two data sets as described by Eq. (14). The limiting spatial frequency for the resulting 4-D DSA temporal projection, $p_{4\mathrm{D}}$, to be backprojected is therefore described by Eq. (15). The next step is to backproject the temporal projection, $p_{4\mathrm{D}}$, to the object plane. The scaling of the image detected by the detector at the image plane is scaled into the object plane and is accomplished through the backprojection operation. The magnification, $m$, represents the geometric scaling factor. The backprojection operation is not a simple 1-D operator, but acts to project the data along direct rays between the image plane and the source. However, for the purposes of discussion and to aid in clarity of understanding the backprojection operation effects on 4-D DSA, the backprojection operation shall be represented as a scaling of image data to the object plane. This analysis is akin to taking a 1-D centered radial sample of the transverse plane such that at the angle that is acquired, the central slice object data axis is exactly parallel with the detector as shown in Fig. 3.

The effect of this is shown as scaling the projection of an object of width $a_I$ from the image plane $I$ to the object plane $O$, as shown in Fig. 3, resulting in an object width of $a_O$. This results in a frequency scaling in frequency space and is described by Eq. (16) and application to 4-D DSA in Eq. (17).

The last operation to be performed is that of multiplying the constraining image with the backprojection of the temporal projection. A multiplication in image space is a convolution in frequency space. The wire object in the constraining image and the weighting volume image were modeled with a pair of Jinc functions¹⁷ to aid in determining the effects of the multiplication in image space. Various sizes were used in the model and resulted in the maximum spatial frequency tracking with object that had the highest spatial frequency or the smaller of the two objects in image space. Therefore, it was determined that maximum spatial frequency of the multiplication of the constraining image with the weighting volume is the maximum of the two spatial frequencies of the two objects when convolved in frequency space as described by Eq. (18) and with application to 4-D DSA shown in Eq. (19). The resulting limiting resolution of 4-D DSA, $f_{4\mathrm{D}}$, is the maximum of the limiting resolution of the constraining image, $i_C$, and the minimum of the limiting resolution of the blurring kernel and the temporal projection data, $p_t$, scaled to the image plane by the scale factor, $m$, as shown in Eq. (20). The 4-D DSA process includes a number of nonlinear operations such as thresholding and division of the temporal projections by reprojections of the constraining volume. Although 4-D DSA is a nonlinear system, it is treated as linear at an operating point. The standard approach¹⁸ of scanning a fine highly attenuating wire centered in and perpendicular to the transverse plane was performed. The effects of the 4-D DSA reconstruction on volumetric spatial resolution were investigated using both an ISPH and the scan of a physical 76 micron (um) diameter tungsten wire centered in a cylindrical supporting structure. ISPH was modeled after the PPH as cylinder centered in the transverse plane with the axis parallel to the axis of rotation extending throughout the entire volume of interest. The physical phantom is shown in Fig. 4. The constraining image inherits the reconstruction parameters from the 3-D DSA as it is a 3-D DSA where a threshold has been applied. The spatial resolution of a 4-D temporal volume can be obtained by the same means as for the 3-D DSA. The 3-D DSA is compared with the 4-D temporal frames (temporally enhanced 3-D volumes) of 4-D DSA reconstructions. This was done on a volumetric basis using a single transverse slice of a physical or digital phantom of a small wire. Resulting PSF data were radially averaged and Fourier analysis performed to generate averaged MTF data for each 3-D DSA volume and 4-D DSA temporal slice. The limiting spatial resolution was automatically determined by finding the spatial resolution when 10% of the magnitude at zero spatial frequency was reached. Numerical simulations and reconstruction of the real-world phantom were performed using a combination of MATLAB^® (The Mathworks Inc., Natick, Massachusetts) and CUDA (NVIDIA Corp. Santa Clara, California) based software. The spatial resolution at the center of the central slice of the C-Arm CT biplane system (Artis zee, Siemens Healthcare, Forchheim, Germany) used is determined by Eq. (21) and is based on the Nyquist sampling criteria using the detector spacing ($du$), imaging geometry source to image distance (SID), and source to object distance (SOD) to determine the minimal detectable distance achievable in the transverse plane. This equation does not account for blurring effects due to focal spot, geometric instabilities of the C-Arm, or projection filtering. The minimal resolvable distance calculation is shown in Eq. (22). The scan geometry was determined using a standard 4-D DSA acquisition procedure without zoom of the detector/C-Arm with SID set to the maximum of 1200 mm and SOD held constant at 750 mm resulting in the magnification factor $m$, and collimation set to maximum field of view with $2\times 2$ binning yielding a detector resolution of $1240\times 960$ detector elements with isotropic detector element size of 0.3080 mm ($2480\times 1920$ native nonbinned resolution with 0.154 mm isotropic detector elements). The minimal resolvable distance is calculated with Eq. (22), which was used to ensure the reconstruction grid was smaller than this dimension. The selection of the maximum allowable wire size was found to be 86 micron using the approach of Nickoloff.¹⁹ The reconstruction grid (slice) was 512 by 512 voxels with a 0.0376 mm isotropic voxel size. The voxel size was made considerably smaller than the minimum resolvable object size, $d_{\min }$, to ensure proper reconstruction and 4-D DSA reconstruction artifacts could be properly investigated. Off axis and off central slice resolution was not investigated in this research and is a topic for future investigation and research. The physical wire phantom was built using a Scientific Instruments (New Jersey) 0.076 mm tungsten wire part number W91 surrounded in an air-filled thin-walled plastic cylinder, 47 mm O.D., as a supporting structure as shown in Fig. 3. The phantom was then scanned on C-Arm CT biplane system. The scan procedure was 8.2 s in duration providing 248 projections, at a frame rate of 30 fps, at 60.4 kVp and 168 mAs.

The phantom was placed in the field of view for the mask run and removed from the field of view for the fill run. This was done to allow the same processing pipeline to be used for current 4-D DSA processing while also setting the automatic exposure control to the object during the mask run. The mask and fill projections were then interchanged to allow correct subtraction and sign of the value in the projections. Normalized projections using Eq. (4) were used in this analysis.

The steps to generate MTFs are as follows. Select a backprojection filtering kernel and perform a reconstruction at the correct zoom factor to satisfy the Nyquist criteria and reconstruction grid requirements. Acquire a central transverse slice of the temporal frame of the reconstructed object and repeat the process for each 4-D volume. Sum all transverse slice time frames extracted from each 4-D time frame volume to make a composite PSF. Find the center of the composite PSF. Preprocess the transverse slices of the temporal time frame if necessary. Radially average the temporal slice in the spatial domain using the center coordinates found using the composite time frame. Generate the measured modulation transfer function (MTF) from the radially averaged PSF by taking the fast Fourier transform (FFT). Then take the absolute value of the result of the FFT. The system MTF is then determined by dividing the measured MTF by an estimate of the MTF of the object (Jinc). The maximum spatial resolution is determined from the resulting MTF plot by finding the first zero or value on the frequency axis where a ten percent of the zero frequency value is reached.

Generation of a composite frame stated above was performed to aid in finding an accurate center of the PSF. Composite frame generation was performed to avoid projection angle dependent modulation in some reconstructed frames. The modulation was induced by the backprojection of blurred temporal projection data where blurring as a result of the blurring kernel, of size zero or two, is less than that induced blurring due to geometric instability. Preprocessing of the transverse slice as stated above can include cropping the image from 512 to 64 voxels on center to avoid any ripples in the frequency response curve near zero spatial frequency due to the standard deviation increasing with increasing $N$ as described in Ref. 17.