Estimation of High Framerate Digital Subtraction Angiography Sequences at Low Radiation Dose

2.1. Phase Decomposition using Independent Component Analysis

In order to obtain an estimate of the contrast progression per image, we start by decomposing the DSA sequence into arterial, capillary and venous phases as illustrated in Figure 2. Because blood flow propagation behaves differently during these three phases, this decomposition permits us to adapt the interpolation mechanism to each phase. We perform the decomposition using Independent Component Analysis (ICA) [4]. ICA is a statistical method that extracts subcomponents of a 1D-signal under the assumption that these subcomponents are independent. The DSA images composing a sequence are first stacked and vectorized to produce a unique 1D-signal. Using ICA on this signal (with three classes) will generate three separate signals corresponding to the arterial, capillary and venous phases. These signals are transformed back to 2D to obtain three distincts images which encode signal information, as well as noise and outliers. To clearly separate meaningful information from the noise and outliers, we use image-base binary thresholding on the histogram distribution. The three thresholded images are used as binary masks to estimate time-density curves (TDCs) [4]. Using the TDCs, we can estimate contrast agent volume that runs through the vascular network between two consecutive images. This amount is normalized between 0 and 1 w.r.t the phase’s peak volume. Assuming a sequence ${\left\{I_i\right\}}_{i=1}^N$ of N DSA images, we can now define the function $v_{(i,j)}=p\left({\left\{I_i\right\}}_{i=1}^N,i,j\right)$ that estimates for any consecutive images I_i and I_j, the contrast volume v_i,j, with j ∈(1, N) and j > i. This volume will be used to control the contribution of the input images to the intermediate image.

2.2. Training and Optimization

Given the training set ${\mathcal{T}}^e={\left\{\left(I_j,L_j\right)\right\}}_j$ composed of DSA images I_j and their corresponding binary labels L_j, with $j\in \left(1,N_{{\mathcal{T}}^e}\right)$, we first start by training a region-of-interest extractor e(I; θ_r) with θ_r being the learned parameters for network e. Because a binary segmentation may discard vessels with small diameters and information about the vessels/background boundaries, we use the extractor to generate a per-pixel entropy map M that will associate with each pixel of an image I the probability of that pixel being part of the vessels or the background (See Figure 3). The higher the probability the richer the information around the pixel. We used a patch-based method [1] that optimizes a mean-square error loss function over the parameters θ_r of the network e.

Then given the training set ${\mathcal{T}}^g={\left\{\left(X_j,Y_j\right)\right\}}_j$, we train a generator network g(X; M; θ_g) over the parameters θ_g to interpolate new images. X_j is composed of a pair of images I_k and I_k+2 and the contrast agent volume v_(k,k+2) with $k\in \left(1,N_{{\mathcal{T}}^g}\right)$. Y_j corresponds to the output I_k+1 that consists of a skipped image that represents the true intermediate image. Finally, $\mathcal{M}={\{{\mathbf{\text{M}}}_k\}}_{k=1}^{N_{{\mathcal{T}}^g}}$ is the set of entropy maps. Note that M is not considered an input to train g but will be used in the optimization loss function.

Given an intermediate image I_k and our predicted intermediate image ${\widehat{I}}_k$, the optimization loss function over the parameters θ_g of the network g is as follows:

where α and β are meta-parameters to control the interaction of the loss function components. The reconstruction loss models how good the generation of the image is. We opted for a Charbonnier loss function and confirmed the observations in [17] that it performs better than an ℓ₁ loss function. We penalize this loss using the entropy maps $\mathcal{M}$ to enforce the network to focus on rich vessel information. Moreover, using the reconstruction error alone might produce blur in the predictions. We thus use a perceptual loss [14] to make interpolated images sharper and preserve their details, where ϕ denote the conv4_3 features of an ImageNet pretrained VGG16 model [19].

2.3. Network Details

We adopt a U-Net architecture [18] for both e and g networks. The network e architecture is straightforward and built upon a ResNet34 pretrained model [13] and a Softmax final activation function to generate the entropy map. The network g consists of 6 hierarchies in the encoder, composed of two convolutional and one ReLU layers. Each hierarchy except the last one is followed by an average pooling layer with a stride of 2 to decrease the spatial dimension. There are 6 hierarchies in the decoder part. A bilinear upsampling layer is used at the beginning of each hierarchy, doubling the spatial dimension, followed by two convolutional and ReLU layers. The input consists of a stacked pair of images, while the contrast volume scalar is pointwise added as a 2D features array at the bottom of the network where features are the most dense and spatial dimensions are the least.

For the parameters of the loss function in Eq. 1, we empirically chose, using a validation set, α = 1 and β = 0.01. We optimize the loss using gradient descent [15] for 80 epochs. using mini-batches of size 16 and with a learning rate of 0.001.

2.4. Final Composition

To increase the framerate of the DSA sequence ${\left\{I_i\right\}}_{i=1}^N$ we infer an intermediate image for each pair of successive images following:

where $\left({\widehat{\theta }}_r,{\widehat{\theta }}_e\right)$ are the network parameters found during the training. Using the entropy to enforce learning features on the regions with rich vessel information has the drawback of mis-interpolating the background and creates a visually disturbing effect. Thus, we add a blending step to build the final image ${\widehat{I}}_{(i,i+1)}^f$ that linearly interpolates the pixels with low entropy as follows:

where ⊙ is an element-wise multiplication function and x ∈ (1, w × h) represents a pixel of the entropy map of size w × h. The parameter δ that controls the interpolation is set empirically to δ = 0.6 while η, a threshold parameter to discard regions with high entropy is set to η = 0.1.

Dataset:

Our dataset is composed of 32 DSA sequences for a total of 3216 DSA images (after geometric transformations augmentation). These images were randomly split and 25% of the images were kept for validation. All images were acquired on human subjects using a standard clinical protocol with a biplane General Electric imaging system. We used both frontal and lateral image sequences. The sequences were acquired at 1-3 fps and captured the full cycle of contrast inflow and washout after injection. The dataset ${\mathcal{T}}^e$ is composed of the the totality of images with their corresponding labeled images while the dataset ${\mathcal{T}}^g$ is composed of pair of images with a leave-one-image-out strategy, skipping every other image to be used as output for the model.

Ablation study on the impact of different components design:

We first perform an ablation study to analyze the contribution of each component of our approach. We test the impact of having a phase-constrained model, the impact of using a Charbonnier loss instead of a mean-square loss (MSE) and finally the impact of using the entropy maps. To this end, we train five variants of our model: woP-MSE which is a traditional U-Net wihtout phase constraints and optimized using an MSE loss function, P-MSE and P-eMSE which are phase-constrained models using MSE loss function with and without the entropy maps respectively and finally, P-Ch and P-eCh (full model) which are phase-contrained models using Charbonnier loss function with and without the entropy maps respectively. To quantify the accuracy of our method we use the Peak Signal-to-Noise Ratio (PSNR) as a measure of corruption and noise and the interpolation error (IE), which is defined as the root-mean-squared difference between the ground-truth image and the interpolated image, as a measure of accuracy.

We can observe from Table 1 that removing the phase knowledge from the model harms performance, particularly the interpolation error, while using the Charbonnier loss slightly improves the results w.r.t MSE loss. We also verify that adding the entropy maps improves the predictions, which validate our hypothesis to enforce a learning on regions with rich vessel information.

Comparison with state-of-the-art methods:

We then compare our approach with state-of-the-art methods including neural and non-neural approaches. In addition to our approach (Our), we include a simple bilinear interpolation method (Lin), an optical-flow based interpolation (OF) [12] and a CNN-based method (SConv) [16] for slow motion video generation. We conduct the comparison on four 7.5 fps DSA sequences that are not used during the training. For each triplet of images we leave the middle image out to serve as ground truth. We report the interpolation error in Figure 4 that shows that our model achieves the best performance on all sequences, on all images. We can notice that our method is only slightly impacted by the contrast agent progression over time, as opposed to the other methods. The performance of our model validates the generalization ability of our approach. Furthermore, in addition to the quantitative measurements, Figure 5 shows the visual differences between our and others methods and highlights its efficiency.

Iterative interpolation:

Finally, using our solution we successfully produce high frame rate sequences of up to 17 fps from 3-fps DSA sequences by iteratively interpolating intermediate images from previously predicted images. Figure 6 shows examples from our dataset including arterial, capillary and venous phases with and without the presence of AVMs.