# Neural blind deconvolution for simultaneous partial volume correction and
super-sampling of PSMA PET images

Caleb Sample1,2, Arman Rahmim1,3,4, Carlos Uribe3,4,5, François Bénard3,4,6,
Jonn Wu7,8, Haley Clark1,2,8 1Department of Physics and Astronomy, Faculty of
Science, University of British Columbia, Vancouver, BC, CA
2Department of Medical Physics, BC Cancer, Surrey, BC, CA
3Department of Radiology, Faculty of Medicine, University of British Columbia,
Vancouver, BC, CA
4Department of Integrative Oncology, BC Cancer Research Institute, Vancouver,
CA
5Department of Functional Imaging, BC Cancer, Vancouver, BC, CA
6Department of Molecular Oncology, BC Cancer, Vancouver, BC, CA
7Department of Radiation Oncology, BC Cancer, Vancouver, BC, CA
8Department of Surgery, Faculty of Medicine, University of British Columbia,
Vancouver, BC, CA

###### Abstract

Objective: We aimed to simultaneously mitigate partial volume effects (PVEs)
in prostate-specific membrane antigen (PSMA) positron emission tomography
(PET) images while performing super-sampling. Approach: Blind deconvolution is
a method of estimating the hypothetical ”deblurred” image along with the blur
kernel (related to the point spread function) simultaneously. Traditional
maximum a posteriori blind deconvolution methods require stringent assumptions
and suffer from convergence to a trivial solution. A promising method of
modelling the deblurred image and kernel with independent neural networks,
called ”neural blind deconvolution” was demonstrated on 2D natural images in
2020. In this work, we adapt neural blind deconvolution for PVE correction of
PSMA PET images, along with simultaneous super-sampling. We compare this
methodology with several interpolation methods, using blind image quality
metrics, and test the model’s ability to predict kernels by re-running the
model after applying artificial ”pseudo-kernels” to deblurred images. Main
Results: Our results demonstrate improvements in image quality over other
interpolation methods in terms of blind image quality metrics and visual
assessment. Predicted kernels were similar between patients, and the model
accurately predicted several artificially-applied pseudo-kernels,
Significance: The intrinsically low spatial resolution of PSMA PET leads to
PVEs which negatively impact uptake quantification in small regions. The
proposed method can be used to mitigate this issue, and can be
straightforwardly adapted for other medical imaging modalities.

Keywords: PSMA PET, PVE, neural blind deconvolution, salivary glands,
deblurring

##  1 Introduction:

Prostate-specific membrane antigen (PSMA) positron emission tomography (PET)
is an imaging modality typically used for the detection of prostate cancer
[1]. For these purposes, one side effect is the accumulation of PSMA ligands
in salivary and lacrimal gland tissue, resulting in unwanted dose to these
regions [2, 3, 4]. However, this can be exploited for PSMA PET to identify
intra-salivary gland parenchyma, which is an area of active research in the
context of radiotherapy treatment planning [5, 6, 7].

Quantifying uptake statistics in the relatively small size of salivary glands
is difficult, due to partial volume effects (PVEs) associated with the
intrinsically low resolution of PET imaging [8], which is typically over 5 mm
[9]. PVEs in PET images primarily correspond to spill-over between image
voxels due to the point-spread function of the scanner[10]. Patient motion,
and poisson-distributed noise due to low photon counts, further exacerbates
blur in PET images. While official SNMMI guidelines recommend using the
maximum standard uptake value ($SUV_{max}$) in single voxels for uptake
metrics [11], maximum uptake values are affected by ”spill-in” from
neighbouring voxels and would be expected to differ if PVEs could be
mitigated.

There have been attempts to mitigate partial volume effects in PET images
using analytical methods [10, 12, 13, 14, 15], with most traditional
approaches involving estimation of the point spread function. These methods
are not applicable in the case of collaborative studies involving multiple
scanners, or where the point spread function is unknown.

The ill-posed problem of simultaneously estimating the theoretical ”deblurred”
image, $\vect{x}$, along with the spread-function or blur kernel, $\vect{k}$,
from the original image $\vect{y}$, is referred to as blind deconvolution.

$\vect{y}=\vect{x}\ast\vect{k}$ (1)

Traditional maximum a posteriori (MAP)-based methods require estimations for
the prior distributions of the kernel and deblurred image and special
optimization techniques to avoid convergence towards a trivial solution.

In 2020, Ren et al. [16] developed ”neural blind deconvolution” for estimating
the deblurred image and blur kernel from 2D natural images using two neural
networks which are optimized simultaneously, as opposed to traditional MAP-
based methods which generally employ alternating optimization to avoid a
trivial solution[16]. Their method does not involve a separate training set,
but learns to predict each case independently. Neural blind deconvolution out-
performed other traditional methods based on prediction’s peak signal-to-noise
ratio, SSIM, and error ratio.

In this work, we build off of the development of 2D neural blind deconvolution
[16, 17], implementing changes to the network architecture and optimization
procedures for suitability with 3D PSMA PET medical images. We also
incorporate super-sampling to enhance spatial resolution.

##  2 Methods:

###  2.1 2.1 Overview of Neural Blind Deconvolution

The method of neural blind deconvolution as first demonstrated by Ren et al.
[16] estimates a ”deblurred” version, $\vect{x}$, of an actual image,
$\vect{y}$, by training neural networks on a case-specific basis.
Specifically, two neural networks, $G_{x}$ and $G_{k}$ are trained to predict
$\vect{x}$, and a blur kernel, $\vect{k}$, whose convolution together yields a
close estimate of the original image, $\vect{y}$. The use of neural networks
for image prediction was motivated by the work of Ulyanov et al. on deep image
priors [18]

$G_{x}(\theta_{x})=\vect{x}\,\,\,\,,\,\,\,\,G_{k}(\theta_{k})=\vect{k}\,\,\,\,,\,\,\,\,\vect{y}=\vect{x}\ast\vect{k}$
(2)

where $\theta_{x}$ and $\theta_{k}$ represent trainable model parameters of
$G_{x}$ and $G_{k}$, respectively. The following implicit constraints exist on
network outputs which are satisfied automatically due to the networks’
architectures.

$0\leq G_{x}(\theta_{x})\leq 1\,\,,\,\,G_{k}(\theta_{k})\geq
0\,\,,\,\,\sum_{i}G_{k}(\theta_{k})_{i}=1$ (3)

Previously, [16] [17] such networks were trained using a loss function
including a mean squared error (MSE) fidelity term along with regularization
terms, $R(\vect{x},\vect{k})$, where

$\vect{x},\vect{k}=\arg\min_{\vect{x},\vect{k}}||\vect{x}\ast\vect{k}-\vect{y}||_{2}^{2}+R(\vect{x},\vect{k})$
(4)

The mathematical formulation of our problem differs slightly as we have
adapted $G_{x}$ to yield an output image, $\vect{x}$ which is twice as large
as $\vect{y}$ in each dimension to accomplish super-sampling. Our $\vect{x}$
is then linearly down-sampled to the original image size,
$\vect{x_{\downarrow}}$, before being convolved with the kernel. Furthermore,
we incorporate a component of mean absolute error into our loss function’s
fidelity term.

$\vect{x},\vect{k}=\arg\min_{\vect{x},\vect{k}}||\vect{x_{\downarrow}}\ast\vect{k}-\vect{y}||_{2}^{2}+||\vect{x_{\downarrow}}\ast\vect{k}-\vect{y}||_{1}+R(\vect{x},\vect{k})$
(5)

This can be written in terms of the model parameters $\theta_{x}$ and
$\theta_{k}$ as

$\theta_{x},\theta_{k}=\arg\min_{\theta_{x},\theta_{k}}||G_{x}(\theta_{x})_{\downarrow}\ast
G_{k}(\theta_{k})-\vect{y}||_{2}^{2}+||G_{x}(\theta_{x})_{\downarrow}\ast
G_{k}(\theta_{k})-\vect{y}||_{1}+R(G_{x}(\theta_{x}),G_{k}(\theta_{k}))$ (6)

which can then be optimized via backpropagation of loss gradients.

###  2.2 2.2 Model Architecture

Neural blind deconvolution was originally demonstrated on 2D images using a
”U-Net”[19]-style, symmetric, multiscale, autoencoder network for $G_{x}$ and
a shallow fully-connected network for $G_{k}$. We built off of this
methodology while implementing some notable changes to optimize model
performance on PSMA PET images. A schematic diagram of the network design is
shown in Figure 1.

[width=1.1]network_architecture.png

Figure 1: The blind deconvolution architecture used for deblurring PSMA PET
images is illustrated. $G_{x}$ is an asymmetric, convolutional auto-encoder
network for predicting the deblurred image, $\vect{x}$, and $G_{k}$ is a
fully-connected network for predicing the blur kernel, $\vect{k}$. The
convolution of $\vect{x}$ and $\vect{k}$ is trained to match the original
image by using back-propagated model gradients from the loss function for
iterative optimization.

####  2.2.1 2.2.1 $\mathbf{G_{x}}$ Architecture

All operations involving $G_{x}$ were adapted for use with 3D images. The
network includes only 3 layers in the encoder-decoder series (as opposed to 5
in Ren et al.’s implementation) [16]. A decrease in network layers is
justified by the low resolution of PET images, which requires a smaller
receptive field than needed for high-resolution images to localize edge
features [20]. An additional layer was added to the end of the decoder
section, including an up-sampling operation, and double convolution, to output
a final $\vect{x}$ twice as large in each dimension as the original input.
This design allows $\vect{x}$ to be simultaneously super-sampled and deblurred
within the blind deconvolution process.

Rather than iteratively increasing channels towards deeper regions of the
encoder and decoder, we found that reversing the design such that the channel
count decreases towards inward regions of the network produced favourable
results. For each double convolution, the first convolution decreases the
image size via a stride of 2. We included a high number of skip connection
channels (64), as it was found to improve prediction accuracy.

Furthermore, we scaled the final sigmoid layer used in previous versions of
neural blind deconvolution [16, 17], as a regular sigmoid constrains the
$\vect{x}$ to have the same maximum value as $\vect{y}$. This is not
necessarily true, nor expected for the case of partial volume effect
correction. We therefore scaled the final layer by $3/2$.

####  2.2.2 2.2.2 $\mathbf{G_{k}}$ Architecture

$G_{k}$ is a single layer perceptron, whose architecture Ren et al. [16] found
to outperform the auto-encoder style design of $G_{x}$ for the case of
predicting the kernel. Model input is a 1-dimensional vector of length 3375,
which is fed through a fully connected layer with 16,875 nodes, to an output
of length 3375 which is passed through a SoftMax operation before being
reshaped into a $15\times 15\times 15$ 3D kernel image.

###  2.3 2.2 Loss Function and Optimization

The main optimization procedure consisted of 5000 iterations broken into 3
stages, stage 1: $0\leq$ step $<300$, stage 2: $300\leq$ step $<2000$, stage
3: $2000\leq$ step $<3500$, stage 4: step $\geq 3500$. 5000 iterations was
chosen as it had been used in previous neural blind deonvolution methods, and
furthermore because we found this number of iterations to be sufficient for
the training loss to level out. We implemented a multi-scale pre-training
procedure similar to the method recommended by Kotera et al. [17]. This
involves initially training $G_{x}$ and $G_{k}$ to predict $\vect{x}$ and
$\vect{k}$ using $2\times$ downsampled PSMA PET images, which are then
upsampled by $\sqrt{2}$ and used as input for a second round of pre-training.
The output is then upsampled back to the original size to serve as inputs for
the regular training stage. We modify the procedure further such that after
every up-sampling during pre-training, the models, $G_{x}$ and $G_{k}$ are
trained to predict their own inputs for 100 iterations, which was sufficient
to allow training in the next scale to begin by predicting up-sampled versions
of the predictions in the previous layer. The size of $\vect{k}$ varies with
the size of $\vect{x}$, such that three kernel sizes are used, $7\times
7\times 7$, $10\times 10\times 10$, and the original size $15\times 15\times
15$. For the down-sampled pre-training of $\vect{x}$ and $\vect{k}$, 1000
iterations were used in each stage. This number of iterations was chosen to
allow the pre-training image to approximately converge to its final estimate.
The initial input vector for $G_{k}$ was a Gaussian Kernel with a standard
deviation of two voxels, and the initial input vector for $G_{x}$ was a
pseudorandom-generated array, sampled from a uniform distribution. The loss
function components for each stage are summarized in Table 1. The parameters
of $G_{x}$ and $G_{k}$ were updated using the ADAM optimizer [21].

Table 1: Loss Function components used for optimization of $G_{x}(\theta_{x})$ and $G_{k}(\theta_{k})$. To guide training during the first 300 optimization steps, the Structural Similarity Index Measure (SSIM) is minimized between the deblurred image $\vect{x_{\downarrow}}$ and the registered CT’s texture map. Conditions are invoked on the kernel estimate until the final 1500 steps of optimization. TV loss is introduced after the first 2000 steps, to avoid convergence towards a trivial solution. $0\leq$ step $<300$ |  $300\leq$ step $<2000$ |  $2000\leq$ step $<3500$ |  step $\geq 3500$
---|---|---|---
MSE + MAE | MSE + MAE | MSE + MAE | SSIM
CT Texture SSIM | Kernel | Kernel | TV
Kernel |  | TV |

It was previously reported by [17] that adopting a total variation (TV) loss
term biases outputs towards a trivial solution during the initial stages of
training, but improves image quality as predictions near the final, correct
solution. We experienced similar findings, and therefore implemented a TV loss
term only after the first 2000 iterations. After the first 3500 iterations,
the fidelity term of the loss function transitions from a combination of MSE
and MAE to the structural similarity index metric, which was also found to
boost performance during final stages [17].

During the first 300 iterations of training, we employ registered computed
tomography (CT) images to provide anatomical guidance. This is implemented by
including a SSIM loss term involving $\vect{x}$ and a CT texture map of the
grey level run length matrix (GLRLM) with a long run emphasis, computed with
pyradiomics [22]. We opted to use the texture map as opposed to the regular CT
image as CT texture maps have been demonstrated to have a closer correlation
with PSMA PET uptake [23] than standard images. CT guidance is used to steer
outputs towards the correct solution during the the first 300 steps of
training before being dropped, to prevent competition with the fidelity term
as outputs move closer towards final estimates. A previous study using
traditional MAP-based blind deconvolution of PET images demonstrated the
utility of magnetic resonance images for guiding optimization [24, 25, 26].

A cumbersome challenge in performing a blind deconvolution is avoiding a
trivial solution, $\vect{x}\ast\vect{k}=\vect{x}\ast\delta=\vect{y}$ [27, 28,
29]. It is therefore necessary to impose regularization on kernel predictions.
Previously, [16, 17] applied a MSE constraint to the kernel to avoid
convergence to a delta function. We also applied a MSE term; however, we only
penalized kernel values larger than 0.7, to avoid a trivial solution while not
penalizing small- to mid-range values in the kernel.

###  2.4 2.3 Evaluation Metrics

As our model was developed on true diagnostic PSMA PET images, ground truth
deblurred images are not available, and we therefore adopted a variety of
strategies and metrics for evaluating the quality of predicted deblurred
images and the model’s ability to predict accurate blur kernels.

We compare two quantitative blind image quality metrics between $\vect{x}$ and
$\vect{y}$. First of all, we compare images via the Blind/Referenceless Image
Spatial Quality Evaluator (BRISQUE) score [30], which varies between 0 and
100, indicating the lowest and highest possible quality, respectively.
Secondly, we compare images using the Contrastive Language-Image Pre-training
(CLIP) metric [31] which varies between 0 and 1. BRISQUE takes an input image
and computes a set of features using the distribution of intensity and
relationships between neighbouring voxels. It then predicts its image quality
from its predicted deviation from a natural, undistorted image. CLIP is a
multi-modal (language and vision) neural network trained with millions of
image/text pairs. The model can be used for predicting both image quality and
image caption quality. The clip network predicts image quality by ranking the
applied captions ”good photo” and ”bad photo”. CLIP was compared with numerous
other quality metrics and was found to perform second to only BRISQUE.

To validate the model’s ability to accurately recover the blur kernel, rather
than arbitrary or semi-random shapes imposed by the loss
function/optimization, we generated four types of ”pseudo-kernels” and
convolved them with predicted deblurred images, which were then fed back into
the blind deconvolution model. We could then test the models ability to
predict a kernel by comparing the newly predicted kernels with pseudo-kernels.
Ideally, the model would predict the pseudo-kernel as the blur kernel, and
thus they would be exactly the same. The generated pseudo-kernel shapes
included a regular gaussian, and 3 gaussians skewed in either the x, y, and z
direction. The similarity between pre-constructed and predicted pseudo-kernels
was measured by taking the inner product of normalized kernel images.

The theoretical blur kernel of PET images has a component resulting from the
point spread function of the scanner, as well as other limitations and patient
motion. Since all PET images were acquired on the same scanner, with the same
image collection protocl, it is expected that predicted blur kernels should
have high inter-patient similarity, Kernels were compared between patients
quantitatively by taking the inner product of normalized kernel images.

###  2.5 2.4 Data Set

This study was approved by an institutional review board. Full-body
[18F]DCFPyL PSMA PET/CT images were de-identified for 30 previous prostate
cancer patients (Mean Age 68, Age Range 45-81; mean weight: 90kg, weight range
52kg-128kg). Patients had been scanned, two hours following intravenous
injection, from the thighs to the top of the skull on a GE Discovery MI (DMI)
scanner. PET images were reconstructed using VPFXS (OSEM with time-of-flight
and point spread function corrections) (pixel spacing: 2.73/3.16mm, slice
thickness: 2.8/3.02mm). For appropriate comparison of voxel values between
patients, patient images were all re-sampled to a slice thickness of 2.8 and
pixel spacing of 2.73 using linear interpolation. Helical CT scans were
acquired on the same scanner (kVP: 120, pixel spacing: 0.98mm, slice
thickness: 3.75mm).

Images were cropped to within 6 slices below the bottom of the submandibular
glands and 6 slices above the top of the parotid glands. Cropping was employed
to avoid exceeding time constraints and the 6 GB memory of the NVIDIA GeForce
GTX 1060 GPU used for training. Registered CT images were used for delineating
parotid and submandibular glands for calculating uptake statistics. Limbus AI
[32] was used for preliminary auto-segmentation of the glands on CT images,
which were then manually refined by a single senior Radiation Oncologist, Jonn
Wu. Images were scaled so that voxels ranged between [0,1] prior to
deblurring, and afterwards rescaled to standard uptake values normalized by
lean body mass ($SUV_{lbm}$).

##  3 Results:

For the two tested blind image quality metrics, BRISQUE and CLIP, deblurred
images out-ranked original images upsampled using nearest-neighbours, linear,
quadratic, and cubic interpolation. These results are summarized in Table 2.

justification=raggedright  Original Polynomial Interpolation Proposed Method
$\vect{y}$ $\vect{y_{\uparrow_{1}}}$ $\vect{y_{\uparrow_{2}}}$
$\vect{y_{\uparrow_{3}}}$ $\vect{x}$ CLIP $0.19\pm 0.03$ $0.33\pm 0.05$
$0.32\pm 0.05$ $0.32\pm 0.05$ $0.39\pm 0.04$ BRISQUE $75.8\pm 8.1$ $89.6\pm
5.4$ $85.0\pm 6.3$ $84.8\pm 6.32$ $90.5\pm 7.5$

Table 2: Blind image quality metrics, CLIP and BRISQUE are compared for
original ($\vect{y}$) and deblurred ($\vect{x}$) PSMA PET images. For
appropriate comparison, $y$ was up-sampled by 2 in each dimension to match the
size of $\vect{x}$. In the first column, $\vect{y}$ is up-sampled using
nearest neighbours interpolation to remain visually identical. We also
compared the quality of $\vect{y}$ when up-sampled using linear, quadratic and
cubic interpolation
($\vect{y_{\uparrow_{1}}}$,$\vect{y_{\uparrow_{2}}}$,$\vect{y_{\uparrow_{3}}}$).

[width=1.1]x_k_conv_y.png

Figure 2: From left to right, for an individual patient, a deblurred image
slice, and axial projection of the predicted kernel are shown, along with the
convolution of the deblurred image and kernel, which is trained to match the
original image, on the right. The inferior aspect of the tubarial glands are
visible between the parotid glands.

An example of $\vect{x}$, $\vect{k}$, $\vect{x}\ast\vect{k}$, and ${y}$ are
shown for a single patient in Figure 2 and visual comparisons of original and
deblurred PET image slices including parotid glands and submandibular glands
are shown for 3 different patients in Figure 3

[width=1.1]par_sm_comparisons_axial.png

Figure 3: Axial slices intersecting parotid glands (left) and submandibular
glands (right) for 3 different patients are shown on original and deblurred
images.

Predicted blur kernels displayed low inter-patient variability, having a mean
inner product between different patients of 0.73. A heat-map of all inter-
patient inner products is found in Figure 4.

[width=1]kernel_heatmap.png

Figure 4: The inner product of recovered blur kernels predicted for all 30
patients is represented by a heat map. Blur kernels for PET images are caused
by both scanner limitations as well as patient motion, so kernels are expected
to be similar, but not identical. Diagonal elements are all equal to unity by
definition. The mean inner product, excluding the diagonal terms, is 0.73.

Projections of the mean predicted kernel onto the three standard cartesian
planes are found in Figure 5. Predicted kernels were found to be mostly
symmetric, with small deviations in skewness found from patient-to-patient.
The distribution of the kernel in various patient directions is summarized in
Table 3, with maximum variation found in the anterior-posterior direction.
Kernels were found to display a centred peak, which rapidly fell off to less
than 0.01 when two voxels from the center in any direction.

The mean absolute difference of the total activity found within all voxels of
original and deblurred images varied by less than $0.05\pm 0.07$%. The
predicted blurred image, $\vect{x}\ast\vect{k}$, and the original image,
$\vect{y}$, had a MAE and MSE of $0.013\pm 0.004$ $SUV_{lbm}$ and $0.003\pm
0.001$ $SUV_{lbm}$, respectively.

[width=1]kernel.png

Figure 5: The mean predicted kernel is projected onto the 3 standard patient
planes. The predicted blur kernel was found to be mostly symmetric.

justification=raggedright  Anterior Posterior Superior Inferior Left Right
$0.52\pm 0.07$ $0.48\pm 0.07$ $0.49\pm 0.05$ $0.51\pm 0.05$ $0.51\pm 0.05$
$0.49\pm 0.05$

Table 3: The fraction of the kernel contained within a plane oriented in each
Cartesian direction and intersecting the center of the kernel image, averaged
over all patients. For example, the average blur kernel has a sum of 0.52
anterior to the center point.

Pseudo-kernels applied to deblurred PET images were predicted well, with mean
inner products of predicted and generated pseudo-kernels above 0.90 for each
of the 4 types of pseudo-kernel applied (regular Gaussian, and skewed Gaussian
along 3 patient axes). These results are displayed in Figure 6

[width=1.1]pseudokernel_plots.png

Figure 6: To validate the model’s ability to predict accurate blur kernels, 4
variously skewed pseudo-kernels were applied to deblurred images before re-
running the blind deconvolution. Generated pseudo-kernels ($k^{*}$, left) and
their corresponding predicted kernel ($k$, right) are shown in separate rows
for each of the 4 pseudo-kernel shapes, in axial, coronal, and sagittal image
slices. The dot-product of normalized kernels, $k^{*}\cdot k$, was used to
assess kernel similarity.

Deblurred images appeared to be visually sharper than original images. Since
$\vect{x}$ is double the size of $\vect{y}$, we up-sampled $\vect{y}$ using
nearest-neighbours, linear, quadratic, and cubic interpolation for
comparison’s sake. Side by side images are found in Figure 7.

[width=1.1]interp_comparison.png

Figure 7: Predicted deblurred images, $\vect{x}$ are twice the original
resolution, and we therefore tested whether perceived improvements in image
quality were only due to the size difference by up-sampling the original image
using nearest neighbours, linear, quadratic, and cubic interpolation. The
deblurred image (right) appeared to have the best quality.

Deblurred and original images of parotid glands are shown in axial, coronal
and sagittal plane slices in Figure 8.

Cohort uptake statistics within CT-defined salivary glands, calculated using
both original and deblurred images, are listed in Table 4. Maximum uptake
values in both parotid and submandibular glands were significantly higher in
deblurred images ($p<0.01$), and mean uptake values were insignificantly
higher.

justification=raggedright  Original Deblurred Parotid Submandibular Parotid
Submandibular $\overline{SUV_{max}}$ $16.2\pm 4.6$ $15.9\pm 3.9$ $21.0\pm 6.2$
$20.2\pm 5.2$ $\overline{SUV_{mean}}$ $6.7\pm 1.8$ $7.1\pm 1.9$ $7.1\pm 1.8$
$7.8\pm 2.2$

Table 4: PSMA PET uptake statistics in CT-defined salivary glands are
summarized using both original and deblurred images.

Uptake plots are shown over corresponding lines through parotid and
submandibular glands in Figures 9 and 10. The gradient of uptake from 0.5
$SUV_{lbm}$ to the full width at half-maximum (FWHM) is approximately twice as
steep for deblurred images than their original counter-parts. The deblurred
images display internal heterogeneity within parotid and submandibular glands
that appears unresolved in original images.

[width=1.]multi_angle_comparison.png

Figure 8: PSMA PET image Slices through parotid glands are shown in the axial
(top), coronal (middle) and sagittal (bottom) planes, in both original (left)
and deblurred (right) images.

[width=1]par_fwhm.png

Figure 9: An axial slice through the original (left) and deblurred (right)
PSMA PET images, intersecting the parotid glands is shown, along with a red
line whose corresponding uptake plot is shown below. The deblurred images are
found to have smaller partial volume effects, as demonstrated by a more rapid
ascent to the full width at half maximum, as shown on lateral edges. In this
case, the rise to FWHM was twice as steep in the deblurred images.

[width=1]sm_fwhm.png

Figure 10: An axial slice through the original (left) and deblurred (right)
PSMA PET images, intersecting the submandibular glands is shown, along with a
red line whose corresponding uptake plot is shown below. The deblurred images
are found to have smaller partial volume effects, as demonstrated by a more
rapid ascent to the full width at half maximum, as shown on lateral edges. In
this case, the rise to FWHM was more than twice as steep in the deblurred
images.

##  4 Discussion:

Based on blind image quality metrics and visual appearance, neural blind
deconvolution helped mitigate PVEs and improve overall image quality. Quality
improvements were a demonstrated result of the blind deconvolution process,
and not only up-sampling, as shown by comparing results with original images
up-sampled using various interpolation orders. The intrinsically low spatial
resolution, and consequent PVEs [33], render PVE correction a necessary task
when quantifying uptake in small image regions. Furthermore, PVE correction of
PET images has been shown to increase the statistical significance of
correlations between uptake bio-markers and various clinical endpoints and
prognostic indices [34, 35, 36, 37, 38].

Predicted blur kernels were consistent with anticipated sizes and shapes,
being small (generally no more than 3 voxels across in any direction with
values larger than 0.01) and mostly symmetric. By observing uptake fall-off at
body borders in original images, it appears that PVEs cause small spill-out to
distances larger than 1 voxel, but neural blind deconvolution cannot detect
such small blur components. For model validation, it was essential that
predicted kernels shared similarities between patients. The low but non-zero
inter-patient variability of predicted kernels was as expected, which can be
considered to be an unknown combination of the point spread function, scanner
limitations, patient motion and other artifacts.

An advantage of neural blind deconvolution over traditional blind
deconvolution methods [39] for mitigating PVEs in PET images is that it does
not require prior assumptions of the PSF to be imposed. A review of various
blind deconvolution algorithms (excluding neural blind deconvolution) applied
to natural images [40] found that the shift-invariant assumption for the blur
kernel is often incorrect, which possibly explains a portion of the variance
seen in predicted kernels. However, due to the static nature of the PET
scanner, it is likely that this variance is less pronounced here than in cases
of natural image acquisition using hand-held devices.

The model’s capacity to predict pseudo-kernels convolved with previously
deblurred images validated its ability to predict various kernel shapes rather
than simply Gaussian shapes. All generated pseudo-kernel shapes were predicted
with high accuracy, but it is notable that the highest accuracy (mean inner
product = 0.98) was found for pseudo-kernel’s that were stretched onto
primarily one axial slice, and all kernel’s stretched onto primarily one axis
were predicted with a higher accuracy than the symmetric Gaussian kernel. This
may be due to the slightly larger slice spacing than pixel spacing (ratio of
1.025), as it is unclear what mechanism in the model architecture would result
in increased predictive power of axially flattened kernels.

While in traditional U-Net as well as previous neural blind deconvolution
studies, the channel count increases towards deeper regions of the encoder and
decoder [19, 16, 17], we discovered that larger channel counts towards outside
regions of the network, decreasing towards the bottom of the encoder and
decoder, performed most accurately. Furthermore, including a large number of
skip connection channels (64 each) greatly improved performance. This is
likely due to increased high-level feature extraction at the top of the
network allowing for ultimately better image reconstruction, along with
bottle-necking in deeper regions leading to improved learning of important
features [41]. Our modelling power was limited by time and GPU resources, and
we were unable to test whether scaling all channel counts further could
improve results. The model code is available for download at
https://github.com/samplecm/neural_blind_deconv_PSMA.

Image quality appeared higher within the high uptake regions of salivary
glands as well as other low uptake regions. Before adopting the combination of
MAE and MSE used in the loss function, we used only MSE as in previous
literature[16, 17], but found that this lent itself poorly to the highly
skewed distribution of PSMA PET uptake. Adding in MAE allowed appropriate
penalization of discrepancy in lower-uptake regions, leading to higher overall
quality. We found that decreasing the number of autoencoder layers in $G_{x}$
from 5 to 4 or 3 had no perceived affect on the model’s predictive power, and
opened up GPU memory for more channels to be added to the remaining layers. We
experimented with various numbers of hidden layers and channel counts in
$G_{k}$ and found that a single layer with 5 times the number of channels in
the flattened kernel image worked best.

PET imaging suffers from intrinsically low resolution, resulting in apparent
partial volume effects. We have demonstrated neural blind deconvolution to be
a viable post-reconstruction method for mitigating these effects. As opposed
to traditional maximum a posteriori methods for estimating deblurred images,
neural blind deconvolution allows for simultaneous estimation of $\vect{x}$
and $\vect{k}$ without convergence towards a trivial solution. We have built
off of previously demonstrated 2D natural image architectures [16, 17] to
create a suitable architecture for 3D PET images. Inclusion of the CT texture
map for steering optimization during early stages helped guide training
towards a centered kernel, as opposed to an off-centered kernel and deblurred
image, which when convolved, predict a centered image for comparison with the
original image. Use of CT images for guiding blind deconvolution relies on
accurate image registration. The GLRLM with a long-run emphasis was used
rather than the original CT im age, due to its lower contrast which matches
PSMA PET better, and due to our previous findings which demonstrated a strong
correlation between PSMA PET uptake and the CT GLRLM [23]. Initially, a joint
entropy loss function was used rather than SSIM, however, this relies on
voxel-binning which causes problems during gradient calculation for back-
propagation. Regularization of the kernel was necessary for avoiding a trivial
solution. However, once a MSE penalty was applied to kernel values greater
than 0.7, the model converged to non-trivial solutions for all patients.
Maximum kernel values were generally less than 0.4, so our kernel constraint
did not arbitrarily constrain the maximum value found in predicted kernels to
a set value.

PSMA PET uptake is heavily biased towards the salivary glands than other
regions in the head and neck, as seen in figures. A previous study has
demonstrated heterogeneous uptake of PSMA-PET in parotid glands [23], and
mitigating partial volume effects using neural blind deconvolution appears to
make this effect more pronounced. The tubarial glands [42] appeared with
greater definition in deblurred images, as two distinct regions of high
uptake. Maximum uptake regions were better localized in deblurred images,
resulting in increased overall maximum values, as well as slight increases in
whole-gland means. These results are too be expected, due to mitigated spill-
out and spill-in of voxel values around maximum values, and spill-out near
gland edges. The total activity in original and deblurred images remained
constant, which was necessary in the context of PET imaging. Originally, we
included a separate total activity constraint in the loss function, but found
this redundant since our fidelity term accomplished this goal.

Super-sampling during neural blind deconvolution is not limited to double
scaling, and the demonstrated model architecture can easily be amended for
accommodating additional resizing. Furthermore, the methodology could be
adapted for kernels to be first convolved with deblurred images before down-
sampling to compute the fidelity loss. Our methodology was limited by time and
computation resources.

##  5 Conclusion:

In this work, we have adapted neural blind deconvolution for simultaneous PVE
correction and supersampling of 3D PSMA PET images. PVE correction allows for
fine detail recovery within the imaging field, such as uptake patterns within
salivary glands. Leveraging the power of deep learning without the need for a
large training data set or prior probabilistic assumptions, makes neural blind
deconvolution a powerful PVE correction method. Our results demonstrate
improvements in quality metrics of deblurred images over other commonly used
super-sampling techniques. The model code is available online for further
studies: https://github.com/samplecm/neural_blind_deconv_PSMA.

## References

## References

  * [1] A. Afshar-Oromieh, H. Hetzheim, C. Kratochwil, M. Benesova, M. Eder, O. Neels, M. Eisenhut, W. Kubler, T. Holland-Letz, F. Giesel, W. Mier, K. Kopka, and U. Haberkorn, “The theranostic psmaligand- 617 in the diagnosis of prostate cancer by pet/ct: Biodistribution in humans, radiation dosimetry, and first evaluation of tumor lesions,” J Nucl Med, vol. 56, no. 11, pp. 1697–1705, 2015.
  * [2] J. Trover, M. Beckett, and G. Wright, “Detection and characterization of the prostate-specific membrane antigen (psma) in tissue extracts and body fluids,” Int J Cancer, vol. 62, pp. 552–558, 1995.
  * [3] R. Israeli, C. Powell, J. Corr, W. Fair, and W. Heston, “Expression of the prostate-specific membrane antigen,” Cancer Res, vol. 54, pp. 1807–1811, 1994.
  * [4] P. Wolf, N. Freudenberg, P. Buhler, K. Alt, W. Shultze-Seemann, U. Wetterauer, and U. Elsasser-Beile, “Three conformational antibodies specific for different psma epitopes are promising diagnostic and therapeutic tools for prostate cancer,” Prostate, vol. 70, pp. 562–569, 2010.
  * [5] H. Clark, S. Thomas, S. Reinsberg, V. Moiseenko, A. Hovan, and J. Wu, “Heterogeneous radiotherapy dose-outcomes response in parotid glands,” Converg. Sci. Phys. Oncol, 2018.  https://doi.org/10.1088/20571739/aac8ea/10.1126/scitranslmed.aac4441.
  * [6] P. van Luijk, S. Pringle, J. O. Deasy, V. V. Moiseenko, H. Faber, A. Hovan, M. Baanstra, H. P. van der Laan, R. G. J. Kierkels, A. van der Schaaf, M. J. Witjes, J. M. Schippers, S. Brandenburg, J. A. Langendijk, J. Wu, and R. P. Coppes, “Sparing the region of the salivary gland containing stem cells preserves saliva production after radiotherapy for head and neck cancer,” Science Translational Medicine, vol. 7, no. 305, pp. 305ra147–305ra147, 2015.  https://doi.org/10.1126/scitranslmed.aac4441.
  * [7] F. Buettner, A. Miah, S. Gulliford, E. Hall, K. Harrington, S. Webb, M. Partridge, and M. Nutting, 2012.
  * [8] H. Marquis, K. Willowson, and D. Bailey, “Partial volume effect in spect & pet imaging and impact on radionuclide dosimetry estimates,” Asia Ocean J Nucl Med Biol., vol. 11, no. 1, pp. 44–54, 2023.
  * [9] T. Pan, S. Einstein, S. Kappadath, K. Grogg, C. Lois Gomez, A. Alessio, W. Hunter, G. El Fakhri, P. Kinahan, and O. Mawlawi, “Performance evaluation of the 5-ring ge discovery mi pet/ct system using the national electrical manufacturers association nu 2-2012 standard,” Med Phys., vol. 46, pp. 3025–3033, Jul 2019.
  * [10] K. Erlandsson, P. Buvat, I andPretorius, B. Thomas, and B. Hutton, “A review of partial volume correction techniques for emission tomography and their applications in neurology, cardiology and oncology,” Phys Med Biol., vol. 57, Nov 2012.
  * [11] W. Fendler, M. Eiber, and M. e. a. Beheshti, “68ga-psma pet/ct: Joint eanm and snmmi procedure guideline for prostate cancer imaging: version 1.0,” Eur J Nucl Med Mol Imaging, vol. 44, p. 1014–1024, 2017.
  * [12] N. M. Alpert, A. Reilhac, T. C. Chio, and I. Selesnick, “Optimization of dynamic measurement of receptor kinetics by wavelet denoising,” NeuroImage, vol. 30, no. 2, pp. 444–451, 2006.
  * [13] F. Šroubek, M. Šorel, J. Boldyš, and J. Šroubek, “Pet image reconstruction using prior information from ct or mri,” in 2009 16th IEEE International Conference on Image Processing (ICIP), pp. 2493–2496, 2009.
  * [14] S. Guérit, L. Jacques, B. Macq, and J. A. Lee, “Post-reconstruction deconvolution of pet images by total generalized variation regularization,” 2015\.
  * [15] T. M. Hansen, K. Mosegaard, S. Holm, F. L. Andersen, B. M. Fischer, and A. E. Hansen, “Probabilistic deconvolution of pet images using informed priors,” Frontiers in Nuclear Medicine, vol. 2, 2023.  http://doi.org/10.3389/fnume.2022.1028928.
  * [16] D. Ren, K. Zhang, Q. Wang, Q. Hu, and W. Zuo, “Neural blind deconvolution using deep priors,” in 2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pp. 3338–3347, 2020.  http://doi.org/10.1109/CVPR42600.2020.00340.
  * [17] J. Kotera, F. Sroubek, and V. Smidl, “Improving neural blind deconvolution,” in 2021 IEEE International Conference on Image Processing (ICIP), pp. 1954–1958, 2021.  http://doi.org/10.1109/ICIP42928.2021.9506502.
  * [18] D. Ulyanov, A. Vedaldi, and V. Lempitsky, “Understanding important features of deep learning models for transmission electron microscopy image segmentation,” International Journal of Computer Vision, 2020.
  * [19] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” 2015.
  * [20] J. Horwath, D. Zakharov, R. Megret, and E. Stach, “Understanding important features of deep learning models for transmission electron microscopy image segmentation,” 12 2019.
  * [21] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” 2017.
  * [22] J. J. M. van Griethuysen, A. Fedorov, C. Parmar, A. Hosny, N. Aucoin, V. Narayan, R. G. H. Beets-Tan, J. C. Fillion-Robin, S. Pieper, and H. J. W. L. Aerts, “Computational radiomics system to decode the radiographic phenotype,” Cancer Research, vol. 77, no. 21, p. e104–e107, 2017.
  * [23] S. Caleb, U. Carlos, A. Rahmim, F. Benard, J. Wu, and H. Clark, “Heterogeneous uptake of psma-ligand inside parotid glands and its correlation with ct texture features.” unpublished, 2023.
  * [24] S. Guérit, A. González, A. Bol, J. A. Lee, and L. Jacques, “Blind deconvolution of pet images using anatomical priors,” 2016.
  * [25] A. G. Gillman, M. H. Quah, E. Papoutsellis, K. Thielemans, P. Bourgeat, C. Delplancke, J. Fripp, E. Pasca, and J. A. Dowling, “Anatomically-guided deconvolution of pet using directional total variation regularization,” medRxiv, 2023.
  * [26] T.-A. Song, F. Yang, S. R. Chowdhury, K. Kim, K. A. Johnson, G. El Fakhri, Q. Li, and J. Dutta, “Pet image deblurring and super-resolution with an mr-based joint entropy prior,” 2019.
  * [27] D. Perrone and P. Favaro, “Total variation blind deconvolution: The devil is in the details,” in 2014 IEEE Conference on Computer Vision and Pattern Recognition, pp. 2909–2916, 2014.  http://doi.org/10.1109/CVPR.2014.372.
  * [28] S. Chaudhuri, R. Velmurugan, and R. Rameshan, Blind Deconvolution Methods: A Review, pp. 37–60.  Cham: Springer International Publishing, 2014.  https://doi.org/10.1007/978-3-319-10485-0_3.
  * [29] J. Yu, Z. Chang, and C. Xiao, “Blur kernel estimation using sparse representation and cross-scale self-similarity,” Multimed Tools Appl 78, vol. 13, p. 18549–18570, 2019.
  * [30] A. Mittal, A. K. Moorthy, and A. C. Bovik, “No-reference image quality assessment in the spatial domain,” IEEE Transactions on Image Processing, vol. 21, no. 12, pp. 4695–4708, 2012.  https://doi.org/10.1109/TIP.2012.2214050.
  * [31] J. Wang, K. C. K. Chan, and C. C. Loy, “Exploring clip for assessing the look and feel of images,” 2022.
  * [32] J. Wong, V. Huang, J. Giambattista, T. Teke, C. Kolbeck, J. Giambattista, and S. Artchian, “Training and validation of deep learning-based auto-segmentation models for lung stereotactic ablative radiotherapy using retrospective radiotherapy planning contours,” Front Oncol, vol. 11, no. 626499, 2021.
  * [33] V. Bettinardi, I. Castiglioni, E. Bernardi, and M. C. Gilardi, “Pet quantification: strategies for partial volume correction,” Clinical and Translational Imaging, vol. 2, no. 3, pp. 199–218, 2014.  Color figures online at http://link.springer.com/article/10.1007/s40336-014-0066-y .
  * [34] F. Gallivanone, C. Canevari, I. Sassi, V. Zuber, A. Marassi, L. Gianolli, M. Picchio, C. Messa, M. C. Gilardi, and I. Castiglioni, “Partial volume corrected 18f-fdg pet mean standardized uptake value correlates with prognostic factors in breast cancer.,” The quarterly journal of nuclear medicine and molecular imaging : official publication of the Italian Association of Nuclear Medicine (AIMN) [and] the International Association of Radiopharmacology (IAR), [and] Section of the Society of…, vol. 58 4, pp. 424–39, 2014.
  * [35] M. Hatt, D. Visvikis, O. Pradier, and C. Rest, “Baseline f-18-fdg pet image-derived parameters for therapy response prediction in oesophageal cancer,” European journal of nuclear medicine and molecular imaging, vol. 38, pp. 1595–606, 05 2011.
  * [36] M. Hatt, A. L. Pogam, D. Visvikis, O. Pradier, and C. C. L. Rest, “Impact of partial-volume effect correction on the predictive and prognostic value of baseline 18f-fdg pet images in esophageal cancer,” Journal of Nuclear Medicine, vol. 53, no. 1, pp. 12–20, 2012.
  * [37] M. Hatt, D. Groheux, A. Martineau, M. Espié, E. Hindié, S. Giacchetti, A. de Roquancourt, D. Visvikis, and C. C.-L. Rest, “Comparison between 18f-fdg pet image–derived indices for early prediction of response to neoadjuvant chemotherapy in breast cancer,” Journal of Nuclear Medicine, vol. 54, no. 3, pp. 341–349, 2013.
  * [38] K. OHTAKA, Y. HIDA, K. KAGA, S. OKAMOTO, T. SHIGA, N. TAMAKI, J. MUTO, R. NAKADA-KUBOTA, S. HIRANO, and Y. MATSUI, “Outcome analysis of 18f-fluorodeoxyglucose positron–emission tomography in patients with lung cancer after partial volume correction,” Anticancer Research, vol. 33, no. 11, pp. 5193–5198, 2013.
  * [39] S. Guérit, A. Gonzalez, A. Bol, J. A. Lee, and L. Jacques, “Blind deconvolution of pet images using anatomical priors,” ArXiv, vol. abs/1608.01896, 2016.
  * [40] A. Levin, Y. Weiss, F. Durand, and W. T. Freeman, “Understanding and evaluating blind deconvolution algorithms,” in 2009 IEEE Conference on Computer Vision and Pattern Recognition, pp. 1964–1971, 2009.
  * [41] I. Goodfellow, Y. Bengio, and A. Courville, Deep Learning.  MIT Press, 2016.  http://www.deeplearningbook.org.
  * [42] M. Valstar, B. de Bakker, R. Steenbakkers, K. de Jong, L. Smit, T. Klein Nulent, R. van Es, I. Hofland, B. de Keizer, B. Jaspers, A. Balm, A. van der Schaaf, J. Langendijk, L. Smeele, and W. Vogel, “The tubarial salivary glands: A potential new organ at risk for radiotherapy,” Radiother Oncol., vol. 154, pp. 292–298, 2021.