# Pixel-wise Color Constancy via Smoothness Techniques in Multi-Illuminant
Scenes

###### Abstract

Most scenes are illuminated by several light sources, where the traditional
assumption of uniform illumination is invalid. This issue is ignored in most
color constancy methods, primarily due to the complex spatial impact of
multiple light sources on the image. Moreover, most existing multi-illuminant
methods fail to preserve the smooth change of illumination, which stems from
spatial dependencies in natural images. Motivated by this, we propose a novel
multi-illuminant color constancy method, by learning pixel-wise illumination
maps caused by multiple light sources. The proposed method enforces smoothness
within neighboring pixels, by regularizing the training with the total
variation loss. Moreover, a bilateral filter is provisioned further to enhance
the natural appearance of the estimated images, while preserving the edges.
Additionally, we propose a label-smoothing technique that enables the model to
generalize well despite the uncertainties in ground truth. Quantitative and
qualitative experiments demonstrate that the proposed method outperforms the
state-of-the-art.

Index Terms—  Color Constancy, Multi-Illuminant Estimation, Total Variation,
Bilateral Filter, Label-smoothing

##  1 Introduction

Color constancy focuses on eliminating the effects of illumination on the
observed colors in captured scenes [1]. The human visual system is inherently
developed to perceive accurate surface colors to some degree, despite varying
illuminations causing different color casts [2]. Inspired by this ability,
digital cameras aim to counterbalance these color casts in images as if they
are captured under canonical illuminant [3]. This process is known as white-
balancing, ensuring color constancy.

In a digital camera, the individual pixel $(x,y)$ of a captured image consists
of three elements, spectral power distribution (SPD) of the light source
$\phi(x,y,\lambda)$, spectral reflectance of the surface $S(x,y,\lambda)$, and
camera sensor spectral sensitivity $C(\lambda)$, all are through the visible
spectrum, $w$, based on the wavelength of the light source, $\lambda$ [4]. The
image formulation can be expressed as follows:

$\displaystyle\textbf{I}_{\text{RGB}}(x,y)=\int_{w}\Phi(x,y,\lambda)S(x,y,\lambda)C(\lambda)\,d\lambda,$
(1)

where $\textbf{I}_{\text{RGB}}(x,y)$ is the RGB image taken under an unknown
illuminant. If we assume that the illumination is uniform and single across
the entire image, von Kries transformation [4] holds for white-balancing of
the image. This transformation is expressed as:

$\displaystyle\mathbf{I}_{\text{RGB}}^{\text{a}}$
$\displaystyle=\mathbf{I}^{\text{a}}\cdot\mathbf{I}_{\text{RGB}}^{\text{C}},\hskip
10.00002pt\mathbf{I}^{\text{a}}=\begin{bmatrix}e_{1}&0&0\\\ 0&e_{2}&0\\\
0&0&e_{3}\end{bmatrix},$ (2)

where $\mathbf{I}_{\text{RGB}}^{\text{a}}$ is the image captured under unknown
illuminant $a$, $\mathbf{I}_{\text{RGB}}^{\text{C}}$ is the image under
canonical illumination, and $\mathbf{I}^{\text{a}}$ is the white-balance gain
matrix [5]. However, uniform and single illumination assumptions are
frequently violated in real-world applications, including indoor areas where
multiple light sources are present [6, 7, 8, 9]. Therefore, the color casts
under the exposition of multiple light sources cannot be corrected with only
one von Kries transformation, but a pixel-wise varying transformation [10].
For instance, in the presence of two illuminants, the white-balance gain
matrix is modified as follows:

$\displaystyle\mathbf{I}^{\text{ab}}(x,y)=\begin{bmatrix}e_{1}(x,y)&0&0\\\
0&e_{2}(x,y)&0\\\ 0&0&e_{3}(x,y)\end{bmatrix},$ (3)

where $\mathbf{I}^{\text{ab}}(x,y)$ is the pixel-wise white-balance gain
matrix, referred to as the illumination map in this work, under two
illuminants, specifically denoted as $a$ and $b$ light sources. In the context
of multi-illuminant color constancy, the goal is to map an input image to a
white-balanced image, achieved by estimating the illumination map affected by
multiple sources [10].

There has been many research performed on color constancy in past years [11,
12, 13, 14]. With the progress in deep learning methodologies, several
learning-based illumination estimation methods have been suggested [9, 10, 15,
16, 17, 18, 19, 8]. One of the techniques is to divide input images into
patches and estimate the illuminants within each patch using Convolutional
Neural Networks (CNNs) [8]. While this method stands as the pioneer in
employing CNNs for addressing color constancy under multiple illuminants, like
most patch-wise models, it does not capture the spatial dependencies among
patches and fails to preserve the smoothness of illumination changes in non-
overlapping patches. A recent work in [10] presented a Generative Adversarial
Network (GAN)-based approach by introducing an additional angular loss term to
estimate the illumination map. While this is an innovative solution, the
estimated images frequently introduce noticeable spatial artifacts and
compromise naturalistic images. The method in [19] constructed a multi-
illuminant dataset and proposed U-Net [20] for per-pixel white-balancing.
Despite the dataset being extensive for multi-illuminant color constancy
research, the absence of object reflectivity leads to inaccuracies in the
ground-truth model. Another work in [9] focused on the detection of achromatic
pixels and surface color similarity to estimate local illuminants and predict
the similarity between object colors. However, the overall performance of this
approach does not surpass the work in [19].

In this paper, we tackle the challenge of multi-illuminant color constancy by
ensuring pixel-wise smoothness with Total Variation [21]. This choice is
motivated by the observation that illuminants change smoothly across the image
and remain consistent between nearby pixels [22]; which is missed in most
existing works. Furthermore, we utilize a bilateral filter [23] to retain
edges and maintain pixel-wise consistency, thereby improving the overall
visual realism of the resulting images. To train and evaluate our multi-
illuminant estimation model, we employ the LSMI dataset [19]. To mitigate the
possible inaccuracies in the ground truth illumination model caused by the
lack of surface reflectivity, we propose using a label-smoothing technique
[24].

In summary, the main contributions of this paper are as follows:

  * •

We introduce a pixel-level learning-based illumination estimation model to
address the multi-illuminant color constancy problem.

  * •

The method we propose ensures that illuminant changes remain consistent among
neighboring pixel locations, achieved through the application of total
variation loss during the training of the estimation model.

  * •

After the training, we apply bilateral filtering to the estimated image,
delivering denoising across the entire image while enhancing photorealism.

  * •

To address possible noises in the ground truth illumination model, we propose
a label-smoothing approach by introducing Gaussian noise injected into the
ground truths.

##  2 Related Work

Many earlier methods in color constancy assume that there is one global
illuminant in the scene, and the uniformity of colors is preserved. Previous
studies [25, 26] are based on the average color of the image, gray, which is
used to estimate the illuminant. White Patch [27] and Max-RGB [28] address the
earlier statistical-based approaches to the color constancy problem. While
these methods do not require a dataset, they do not deliver visually
satisfying results if multiple illuminations are present. Recent statistical
approaches include convex functions [29], pixel-level probabilistic model
[30], and additional mathematical terms by improving Gray-world and Max-RGB
assumptions [31, 32] to estimate the illuminant of an image. In contrast,
learning-based approaches have emerged as leading solutions to color constancy
problems recently. Those studies [33, 34] use CNNs to estimate the illuminant.
However, these methods depend on the assumption of uniform color constancy.

Uniformity of illumination is repeatedly violated in real-world settings [6,
9, 7]. Therefore, non-uniform illuminant estimation, referred to as multi-
illuminant estimation, is examined in recent studies [35, 36, 7, 15, 8, 10,
37, 19, 38, 9]. These studies are limited due to the lack of available multi-
illuminant datasets. For instance, in [10], the author proposes the AngularGAN
approach to estimate multi-illuminant colors in images whose ground truths are
generated by tinting the illumination maps. In [19], they propose the Large
Scale Multi-Illuminant (LSMI) Dataset and a learning-based algorithm. The
dataset consists of real-life images illuminated by multiple sources, along
with pixel-wise ground truth illumination maps. However, their learning-based
method overlooks the absence of object reflectivity, resulting in possible
noisy ground truth modeling. Another study in [38] proposes an N-white
balancing to determine white points of illuminants in different segments of
the image. A recent work [9] suggests a multi-task learning method. They
employ auxiliary tasks to detect achromatic pixels and to predict surface
color similarity. Most of these methods evaluate their models with the LSMI
dataset. However, they neglect to consider pixel-wise consistent illumination
variation.

##  3 Proposed Method

###  3.1 Method Overview

The proposed method models the problem by predicting the pixel-wise
illuminants for each input image, through supervised learning. To ensure
illuminant smoothness and enhance visual quality over neighboring pixels, we
integrate a total variation loss during training and apply bilateral filtering
after training. Label smoothing is employed for increasing robustness against
noisy ground truth distribution. See Figure 1 for an overview of the pipeline.

Fig. 1: Overview of the proposed method: The U-Net multi-illuminant estimation
model processes input images and smoothed ground truth images within the log-
chrominance color space. During training, the model learns an illumination
mapping between input and ground truth images by minimizing Total Variation
Loss and L2 loss. The Bilateral Filter is then applied to the estimated
images. The final output is the white-balanced image.

###  3.2 Ground Truth Illumination Model

Under the assumption that the scene is lit by two illuminants, the ground
truth illumination model is defined as follows [19]:

$\displaystyle\textbf{I}{}_{\text{RGB}}^{\text{ab}}(x,y)=[\alpha(x,y)\cdot
l_{a}+[1-\alpha(x,y)]\cdot
l_{b}]\cdot\textbf{I}_{\text{RGB}}^{\text{ab,gt}}(x,y)$ (4)
$\displaystyle\textbf{I}{}_{\text{RGB}}^{\text{ab}}(x,y)=\textbf{{I}}^{\text{ab}}(x,y)\cdot\textbf{I}_{\text{RGB}}^{\text{ab,gt}}(x,y),$
(5)

where $l_{a}$ and $l_{b}$ are chromaticities of $a$ and $b$ illuminants,
$\alpha(x,y)$ is the illumination weight of the illuminant $a$ at pixel
$(x,y)$, $\textbf{I}{}_{\text{RGB}}^{\text{ab}}(x,y)$ is the input RGB image,
$\textbf{I}_{\text{RGB}}^{\text{ab,gt}}(x,y)$ is the ground truth RGB image,
and $\textbf{I}^{\text{ab}}(x,y)$ is the pixel-wise ground truth illumination
map. Thus, $\alpha$ defines $\textbf{I}^{\text{ab}}$ that does the white
balancing operation to the input image.

###  3.3 Color Space Transformation

The proposed algorithm takes a linear-raw RGB input image and the
corresponding ground truth image. Then it converts them to log-chrominance
space. This conversion includes a logarithm that corresponds to the
chromaticity of illumination in terms of camera-sensor response. The
conversion results in a simpler color space than the RGB triplet although the
absolute scale ambiguity is not preserved [39]. The log-chrominance space is
defined as follows:

$\displaystyle
u(x,y)=\log\left(\frac{I_{\text{R}}(x,y)+\epsilon}{I_{\text{G}}(x,y)+\epsilon}\right)$
(6) $\displaystyle
v(x,y)=\log\left(\frac{I_{\text{B}}(x,y)+\epsilon}{I_{\text{G}}(x,y)+\epsilon}\right),$
(7)

where $\epsilon$ is an infinitesimal number to avoid division by zero error
and $I_{R}(x,y),I_{G}(x,y),I_{B}(x,y)$ are spatial pixels of each R, G, and B
channels, respectively. Note that the estimation model utilizes the log-
chrominance color space, uv, for the input, the ground truth, and the
predicted output. Also, the estimated illumination map is transformed back to
the RGB color space.

###  3.4 Estimation Model

The estimation model, based on the U-Net model [20] and detailed in [19],
serves as the baseline model in this work. It comprises 16 convolutional
steps, with 8 dedicated to downsampling and 8 for upsampling. The network
takes the input and ground truth images in log-chrominance space and shrinks
them to 256x256 image size. Following this, the network proceeds to estimate
the illuminants, producing the output image in the same color space.

###  3.5 Label Smoothing

In the context of constructing the ground truth illumination map (referenced
in Section 3.2), the alpha parameter, $\alpha$, plays a crucial role. However,
the ground truths in LSMI dataset [19] are extracted only considering the
reflected light and not the incoming light, thus the surface reflectivity is
not accounted for. Therefore, the actual ground truth is considered to be a
noisy version of the ground truth in [19]. To address this, we propose
applying label-smoothing [24], which not only prevents the network from
becoming over-confident to the noisy ground truth but also improves the
generalization of estimating illumination maps.

Label smoothing is introduced by adding a random Gaussian noise to the
$\alpha$ in every training batch. The distribution has zero mean and standard
deviation equal to one-tenth of the magnitude of the corresponding
$\alpha_{raw}(x,y)$. That is:

$\displaystyle\sigma_{n}(x,y)=\frac{\alpha_{raw}(x,y)}{w_{n}}\quad w_{n}=10$
(8) $\displaystyle\textit{X}(x,y)\sim\mathcal{N}(0,\sigma_{n}^{2}(x,y))$ (9)
$\displaystyle\alpha_{smooth}(x,y)=\alpha_{raw}(x,y)+\textit{X}(x,y),$ (10)

where $w_{n}$ is the smoothing constant, $\sigma_{n}$ is the standard
deviation, $\textit{X}(x,y)$ is the random Gaussian noise, $\alpha_{raw}(x,y)$
is the raw illumination weight, and $\alpha_{smooth}(x,y)$ is the smoothed
illumination weight, all evaluated at the $(x,y)$ pixel.

Fig. 2: Distributions of $\alpha_{raw}(x,y)$ and $\alpha_{smooth}$ parameters

Figure 2 illustrates the distributions of raw alpha and smoothed alpha values
during a single training iteration. The visualization demonstrates that label
smoothing transforms the $\alpha$ distribution to closely resemble a Gaussian
distribution, indicating a smoother distribution of hard labels.

###  3.6 Total Variation Loss

We anticipate that the illumination map, derived from multiple sources, should
demonstrate piece-wise smoothness due to the gradual variation of illuminants
across the image [21]. To accomplish this, the Total Variation (TV) Loss is
proposed. This loss imposes an extra penalty to the overall loss function,
ensuring smoothness across illuminant changes. TV Loss is defined as:

$\displaystyle L_{\text{TV}}=\sum_{x,y\in I}|\nabla_{X}f(x,y)|$
$\displaystyle+|\nabla_{Y}f(x,y)|$ (11) $\displaystyle=\sum_{i,j\in
I}|\hat{\textbf{I}}_{\text{RGB}}^{(i,j)}-\hat{\textbf{I}}_{\text{RGB}}^{(i+1,j)}|+\sum_{i,j\in
I}|\hat{\textbf{I}}_{\text{RGB}}^{(i,j)}-\hat{\textbf{I}}_{\text{RGB}}^{(i,j+1)}|,$
(12)

where $\hat{\textbf{I}}{}_{\text{RGB}}^{(i,j)}$ is the estimated RGB value at
$(i,j)$. The overall loss function consisting of a linear combination of
L2-loss and TV Loss becomes:

$\displaystyle L=L_{2}+\lambda_{\text{TV}}\cdot L_{\text{TV}},$ (13)

where $\lambda_{\text{TV}}$ is the weight of the TV Loss. Overall, the TV Loss
serves to regularize the estimation of the illumination map, accounting for
spatial dependencies.

###  3.7 Bilateral Filtering

To achieve visually realistic white-balanced images, we propose using a
bilateral filter, known for its spatial and range filtering properties [40].
These properties play a pivotal role in maintaining pixel-wise consistency and
preserving edges across the image. The formulation of the bilateral filter is
as follows:

$\displaystyle
BF(\sigma_{s},\sigma_{r})=\frac{1}{w}\sum_{\mathbf{p},\mathbf{q}\in
S}G_{\sigma_{s}}(\|\mathbf{p}-\mathbf{q}\|)\cdot
G_{\sigma_{r}}(|I_{\mathbf{p}}-I_{\mathbf{q}}|)\cdot I_{\mathbf{q}}.$ (14)

In this equation, $I_{\mathbf{p}}$ and $I_{\mathbf{q}}$ represent the
intensity at pixel $\mathbf{p}$ and $\mathbf{q}$, and $w$ is a normalization
factor to ensure the weights sum up to 1. The parameters $G_{\sigma_{s}}$ and
$G_{\sigma_{r}}$ are Gaussian functions assuring smoothing and edge-preserving
with standard deviations $\sigma_{s}$ and $\sigma_{r}$, respectively.
Bilateral filtering in the context of color constancy is given as:

$\displaystyle\hat{\textbf{I}}^{f}=BF(\sigma_{s},\sigma_{r})*\hat{\textbf{I}}.$
(15)

Let $\hat{\textbf{I}}$ represent the output RGB image of the estimation model,
$BF(\sigma_{s},\sigma_{r})$ the bilateral filter with parameters $\sigma_{s}$
and $\sigma_{r}$ set as 75, and ($*$) denotes a convolution with a 9-pixel
diameter. These parameter choices optimize noise filtering and edge
preservation.

##  4 Experiments and Discussion

The proposed method is trained and assessed on the Large Scale Multi-
Illuminant (LSMI) dataset [19], which comprises real-world 7,486 images
acquired by three distinct cameras under 1-3 varied illuminants, containing
approximately 3000 scenes. Our evaluation focuses exclusively on images
captured by the Samsung Galaxy Note 20 Ultra, resulting in a subset of 2360
images illuminated by only two illuminant sources. The resulting dataset is
partitioned into training, validation, and test sets with a ratio of
0.75:0.2:0.05, respectively. Note that color checkers in images are masked
with black color.

In the training process, the estimation model is trained for 2000 epochs with
a batch size of 32. Two different learning rates, namely $0.0001$ and
$0.0005$, are determined based on the weight of the TV Loss. These two
settings give rise to two distinct proposed models, referred to as Pixel-wise
Color Constancy (PWCC), PWCC_v1 and PWCC_v2. Those parameters and various
settings can be investigated in table 1. Moreover, the learning rates exhibit
a decay factor of 800, beginning from the 800th epoch. The training is
performed with 4 NVIDIA GeForce RTX 2080 Ti GPUs.

Table 1: Different settings for the Estimation Model of our method on Two-Illuminant images of LSMI dataset [19]. Method | Filtering | Label Smooth | $\lambda_{\text{TV}}$ | Learning Rate
---|---|---|---|---
PWCC_v1 | Yes | No | 2e-4 | 5e-4
PWCC_v2 | Yes | Yes | 2e-3 | 1e-4
Fig. 3: Visualization results on different two-illuminant images of LSMI
Dataset [19]. From left to right, input image, baseline (LSMI-U-Net [19]),
PWCC_v1, PWCC_v2 predictions, and the ground truth image, respectively.

###  4.1 Evaluation Metrics

The proposed approach is evaluated using a total of 112 two illuminant test
images including corresponding ground truth images. Evaluation is
quantitatively done by calculating recovery angular errors between ground
truth and predicted illumination maps, as it is an intensity-independent error
measure [41] and is commonly used in evaluating the performances of color
constancy algorithms. For qualitative comparison, predicted images of our best
two models along with the input and ground truth images are compared with the
prediction of the baseline method, that is, LSMI-U-Net described in [19].
Quantitative Evaluation We report the mean, median, worst 25%, and best 25% of
angular recovery error, $e_{recovery}$, between the ground truth and the
estimated illumination maps. $e_{recovery}$ is defined as follows:

$\displaystyle
e_{\text{recovery}}(\textbf{I}^{\text{ab}},\hat{\textbf{I}}{}^{\text{ab}})=\cos^{-1}\frac{\textbf{I}^{\text{ab}}\cdot\hat{\textbf{I}}{}^{\text{ab}}}{\|\textbf{I}^{\text{ab}}\|\>\|\hat{\textbf{I}}{}^{\text{ab}}\|},$
(16)

where $\textbf{I}^{\text{ab}}$ and $\hat{\textbf{I}}{}^{\text{ab}}$ are the
ground truth and the estimated illumination maps for an image, and ($\cdot$)
is the dot product.
$e_{\text{recovery}}(\textbf{I}^{\text{ab}},\hat{\textbf{I}}{}^{\text{ab}})$
is calculated for each image and used for the above-mentioned reports.

Table 2: Recovery angular errors on Multi-Illuminant images of LSMI dataset [19] Method | Mean | Median | W.25% | B.25%
---|---|---|---|---
Gray World [26] | 11.3 | 8.8 | 20.74 | 4.93
White Patch [27] | 12.8 | 14.3 | 23.49 | 5.6
LSMI U-Net [19] | 2.31 | 1.91 | 4.24 | 1.01
PWCC_v1 | 2.08 | 1.74 | 3.8 | 0.9
PWCC_v2 | 2 | 1.7 | 3.8 | 0.86

Table 2 reports the recovery angular errors obtained from various methods. We
compare our proposed methods with Gray-World [26], White Patch [27], and the
reproduced results of the state-of-the-art multi-illuminant method LSMI U-Net
[19]. It is observed that our proposed methods, PWCC, outperform the competing
methods. Specifically, PWCC_v2 achieves an impressive mean error of 2, which
is a 13% improvement over the best-competing method. A notable observation is
that not only the proposed method improve mean performance, but it also
significantly improves the results of the worst 25%, which are often the
corner cases.

Qualitative Evaluation Figure 3 showcases three visual examples, from the
proposed methods and the LSMI-U-Net as the baseline. It can be observed that
our method consistently outperforms the baseline, demonstrating better white
balancing across various input images. For instance, in the top row image, the
baseline cannot fully remove the greenish color cast, while the proposed shows
a much better result. Moreover, for the third-row image, the baseline
overcompensates for the reddish color cast, which also removes the natural
color on the wall. However, in some test data instances, such as the second
row, the two methods show competing performances.

##  5 Conclusion

In this paper, a multi-illuminant pixel-wise color constancy approach is
proposed. We highlighted the importance of smooth illumination changes in
natural images and proposed total variation loss for regularization during
training. We proposed to filter the output white-balanced images of the
estimation model with a Bilateral Filter, to guarantee the pixel-wise
consistency across the image and create a visually natural appearance.
Moreover, the usage of the label smoothing technique ensures that the model
generalizes better on the noisy ground truth data. Experimental results
demonstrate that the proposed method outperforms the current state-of-the-art
method. Future work includes proposing an estimation model accounting for
surface reflectivity to enhance the accuracy of ground truth modeling for the
multi-illuminant dataset used in this paper.

##  6 ACKNOWLEDGEMENT

This project has received funding from the European Union’s Horizon 2020
research and innovation program under the Marie Skłodowska-Curie grant
agreement No [101022466], and the NSF-Business Finland Center for Big Learning
(CBL), Advanced Machine Learning for Industrial Applications (AMaLIA) under
Grant 97/31/2023.

## References

  * [1] G. D. Finlayson, S. D. Hordley, and P. M. Hubel,  “Color by correlation: A simple, unifying framework for color constancy,”  Pattern Analysis and Machine Intelligence, 2001.
  * [2] J.T. Nikkanen,  Computational color constancy in mobile imaging,  Ph.D. thesis, Tampere University, 2013.
  * [3] M. Ebner,  Color constancy,  John Wiley & Sons, 2007.
  * [4] J. von Kries,  “Influence of adaptation on the effects produced by luminous stimuli,”  Source of Color Science, 1970.
  * [5] O. Yilmaz,  “Spectral ray tracing for generation of spatial color constancy training data,”  M.S. thesis, 2022.
  * [6] R. O. Dror, A. S. Willsky, and E. H. Adelson,  “Statistical characterization of real-world illumination,”  Journal of Vision, 2004.
  * [7] S. Beigpour, C. Riess, J. Van De Weijer, and E. Angelopoulou,  “Multi-illuminant estimation with conditional random fields,”  IEEE Transactions on Image Processing, 2013.
  * [8] S. Bianco, C. Cusano, and R. Schettini,  “Single and multiple illuminant estimation using convolutional neural networks,”  IEEE Transactions on Image Processing, 2017.
  * [9] Shuwei L., Jikai W., Michael S. B., and Robby T. T.,  “Mimt: Multi-illuminant color constancy via multi-task local surface and light color learning,” 2023.
  * [10] O. Sidorov,  “Conditional gans for multi-illuminant color constancy: Revolution or yet another approach?,”  in Proceedings of the IEEE/CVF workshops, 2019.
  * [11] E. H Land,  “The retinex theory of color vision,”  Scientific american, 1977.
  * [12] S. D. Hordley,  “Scene illuminant estimation: past, present, and future,”  Color Research & Application, 2006.
  * [13] A. Gijsenij, T. Gevers, and J. Van De Weijer,  “Computational color constancy: Survey and experiments,”  IEEE Transactions on Image Processing, 2011.
  * [14] Y. Fan, H. Ren, and H. Huang,  “Fusion temporal color constancy,”  in IEEE ICIP. IEEE, 2022.
  * [15] Y. Hu, B. Wang, and S. Lin,  “Fc4: Fully convolutional color constancy with confidence-weighted pooling,”  in CVPR, 2017.
  * [16] F. Laakom, J. Raitoharju, A. Iosifidis, J. Nikkanen, and M. Gabbouj,  “Color constancy convolutional autoencoder,”  in IEEE Symposium Series on Computational Intelligence, 2019.
  * [17] F. Laakom, N. Passalis, J. Raitoharju, J. Nikkanen, A. Tefas, A. Iosifidis, and M. Gabbouj,  “Bag of color features for color constancy,”  IEEE Transactions on Image Processing, 2020.
  * [18] T. Akazawa, Y. Kinoshita, and H. Kiya,  “Multi-color balance for color constancy,”  in IEEE ICIP, 2021.
  * [19] D. Kim, J. Kim, S. Nam, D. Lee, Y. Lee, N. Kang, H.-E. Lee, B. Yoo, J.-J. Han, and S. J. Kim,  “Large scale multi-illuminant (lsmi) dataset for developing white balance algorithm under mixed illumination,”  in IEEE/CVF ICCV, 2021.
  * [20] O. Ronneberger, P. Fischer, and T. Brox,  “U-net: Convolutional networks for biomedical image segmentation,”  in Medical Image Computing and Computer-Assisted Intervention. Springer, 2015.
  * [21] M. Javanmardi, M. Sajjadi, T. Liu, and T. Tasdizen,  “Unsupervised total variation loss for semi-supervised deep learning of semantic segmentation,”  arXiv:1605.01368, 2016.
  * [22] K. Barnard, G. Finlayson, and B. Funt,  “Colour constancy for scenes with varying illumination,”  in ECCV. Springer, 1996.
  * [23] A. Choudhury and G. Medioni,  “Color constancy using denoising methods and cepstral analysis,”  in ICIP. IEEE, 2009.
  * [24] R. Müller, S. Kornblith, and G. E. Hinton,  “When does label smoothing help?,”  Advances in Neural Information Processing Systems, 2019.
  * [25] G. Buchsbaum,  “A spatial processor model for object colour perception,”  Journal of the Franklin Institute, 1980.
  * [26] J. Cepeda-Negrete and R. E. Sanchez-Yanez,  “Gray-world assumption on perceptual color spaces,”  in Image and Video Technology: Pacific-Rim Symposium, 2013.
  * [27] A. Rizzi, C. Gatta, and D. Marini,  “Color correction between gray world and white patch,”  in Human Vision and Electronic Imaging. SPIE, 2002.
  * [28] E. H. Land and J. J. McCann,  “Lightness and retinex theory,”  Josa, 1971.
  * [29] Z. Abedini and M. Jamzad,  “Single and multiple illuminant estimation using convex functions,”  Computer Vision and Image Understanding, 2023.
  * [30] F. Laakom, J. Raitoharju, A. Iosifidis, U. Tuna, J. Nikkanen, and M. Gabbouj,  “Probabilistic color constancy,”  in ICIP. IEEE, 2020.
  * [31] G. D. Finlayson,  “Corrected-moment illuminant estimation,”  in IEEE ICCV, 2013.
  * [32] O. Ulucan, D. Ulucan, and M. Ebner,  “Color constancy beyond standard illuminants,”  in ICIP. IEEE, 2022.
  * [33] S. Bianco and C. Cusano,  “Quasi-unsupervised color constancy,”  in IEEE/CVF CVPR, 2019.
  * [34] Y.-C. Lo, C.-C. Chang, H.-C. Chiu, Y.-H. Huang, C.-P. Chen, Y.-L. Chang, and K. Jou,  “Clcc: Contrastive learning for color constancy,”  in IEEE/CVF CVPR, 2021.
  * [35] A. Gijsenij, R. Lu, and T. Gevers,  “Color constancy for multiple light sources,”  IEEE Transactions on Image Processing, 2011.
  * [36] H. R. V. Joze and M. S. Drew,  “Exemplar-based color constancy and multiple illumination,”  Pattern Analysis and Machine Intelligence, 2013.
  * [37] I. Domislović, D. Vrsnak, M. Subašić, and S. Lončarić,  “Outdoor daytime multi-illuminant color constancy,”  in International Symposium on Image and Signal Processing and Analysis. IEEE, 2021.
  * [38] T. Akazawa, Y. Kinoshita, S. Shiota, and H. Kiya,  “N-white balancing: White balancing for multiple illuminants including non-uniform illumination,”  IEEE Access, 2022.
  * [39] J. T. Barron,  “Convolutional color constancy,”  in IEEE ICCV, 2015.
  * [40] C. Tomasi and R. Manduchi,  “Bilateral filtering for gray and color images,”  in ICCV. IEEE, 1998.
  * [41] S. D. Hordley and G. D. Finlayson,  “Re-evaluating colour constancy algorithms,”  in ICPR. IEEE, 2004.


We use a similar dataset acquired from <en.boardgamearena.com> as in [Hu et al., 2021]. The dataset consists of 240,954 2-player Hanabi games. We randomly sample 1,000 games to create a validation set and another 4,000 games for the test set. The training set contains the remaining 235,954 games with an average score of 15.88. Each game records the AOH $\tau_i, i\in \{1, 2\}$ for both players. The IL policy $\pi_\theta$ is parameterized by a neural network $\theta$ and is trained to minimize the cross-entropy loss

$$ \mathcal{L(\theta)} = - \mathbb{E}_{\tau_i \sim D}\sum_{t=0}^{T} \pi_\theta(a_i^t|\tau_i^t) $$

with stochastic gradient descent. The training set $D$ is dynamically augmented with color shuffling [Hu et al., 2020] where a random color permutation that changes both observation and action space is applied to each trajectory $\tau_i$ sampled from the dataset before feeding it to the network. Hu et al., 2021 shows that this data augmentation method greatly reduces overfitting and leads to better policies. The network $\theta$ uses the Public-LSTM structure in [Hu et al., 2021], which eliminates the need to re-unroll LSTM on sampled trajectories from the beginning of the game as they share the same public observations as the real trajectory.

To construct the approximate belief model $\hat{\mathcal{B}}$, we train a neural network $\phi$ to predict each player's own hand, which is the only hidden information in Hanabi.
The network predicts each card in hand from oldest to newest auto-regressively. We denote the hidden cards as $\{h_i^{j}\}$ with $h^1_i$ being the oldest and $h^m_i$ being the newest, $m \leq 5$. Then the cross-entropy loss for $\phi$ in an $N$-player setting becomes

\begin{align*}
\mathcal{L(\phi)} = -\mathbb{E}_{\tau \sim \pi_\theta}\left[ \frac{1}{N}\sum_{i=1}^{N} \sum_{t=1}^{T} \sum_{j=1}^{m} \log P_{\phi}(h^j_i | \tau_i^t, h^1_i, \ldots, h_i^{j-1}) \right].
\end{align*}

In Hanabi, it is sufficient to reconstruct $\tau^t$ given $\tau^t_i$ and $\{h^i_j\}$.
The model is trained on infinite stream of data generated by $\pi_\theta$ through self-play to avoid overfitting. We train two variants of the belief model, one with sampled action $\textbf{a} \sim \pi_\theta$ while the other with greedy action $\textbf{a} = \argmax \pi_\theta$. They are used by piKL-SPARTA and piKL-SPARTA-G(reedy) respectively.

We first run piKL-SPARTA on the test set to compare its ability to predict human moves against the IL policy. For each $\tau_i^t$ in the test set, we sample $\frac{K}{|A_i|}$ hands from the belief model $P_{\phi}$ where $K$ is the total number of searches for this step and ${|A_i|}$ is the number of legal actions of the search player. In all our experiments we set $K=10,000$. In practice, we sample $\frac{2K}{|A_i|}$ hands and take the top $\frac{K}{|A_i|}$ samples not contradicting the public knowledge. If the belief model fails to produce any samples that comply with the public knowledge, we revert back to the blueprint. We then compute the expected value for each search action by unrolling the IL policy until the end of the game for each sampled trajectory and compare $a_{\text{pred}} = \argmax_a \pi_{\theta}(a | \tau_i^t) \cdot \exp [\frac{Q_{\pi}(\tau_i^t, a)}{\lambda}]$ against the human move $a^t_i$ from the data. We experiment with different $\lambda$ to study its effect on prediction accuracy.

We also evaluate piKL-SPARTA in self-play to compare its performance under different $\lambda$ against the IL policy. We run piKL-SPARTA and the IL policy on 4,000 games with different seeds for the deck. At each step, the IL policy acts following $a_t \sim \pi_\theta(a | \tau^t_i)$ while piKL-SPARTA acts following Eq. (<ref>) again using $K=10,000$ rollouts.

The $\lambda$ in these experiments are significantly higher than those in Diplomacy or implicitly in MCTS in Chess and Go[Via the relationship $\lambda \approx c_{\text{puct}} \sqrt{N}$ where $N$ is the number of MCTS iterations] because $\lambda$ represents the scale of utility difference that offsets a particular KL penalty, and the range of the utilities $Q_{\pi}(\tau_i^t, a)$ in Hanabi is $[-25, 25]$ whereas the range in other games are either $[0, 1]$ or $[-1, 1]$.

We experiment with both 1p piKL-SPARTA, where one player uses piKL-SPARTA and the other follows the blueprint, and 2p piKL-SPARTA, where both players run the same single-agent version of piKL-SPARTA independently. The latter is not theoretically sound and 2p SPARTA has in past papers produced worse performance than 1p SPARTA [Lerer et al., 2020] because 1p SPARTA assumes that partners are playing according to the blueprint policy while in actuality the partners are playing according to a SPARTA policy. Despite the lack of theoretical soundness, we are interested in whether 2p piKL-SPARTA can obtain empirical improvement over the 1p version, since piKL-SPARTA regularizes towards the blueprint and therefore the mismatch between assuming the partner follows the blueprint and the policy they actually play would likely be less severe.

Lastly, in Dec-POMDPs, it is common to select actions greedily in self-play instead of sampling from a mixed policy, since in fully cooperative settings there is no need to avoid being deterministic or predictable to an adversary. Therefore, we also experimented with piKL-SPARTA-G(reedy) where every agent, including the baseline IL policy, play according to the $\argmax$ of their action distributions at every step, and the belief model for piKL-SPARTA-G is trained on trajectories produced by a greedy IL policy.

§.§ Results

Top-1 test accuracy and self-play score of IL policy and piKL-SPARTA in Hanabi. Blue dot is the IL policy and green dots are piKL-SPARTA with different $\lambda$. The self-play score is evaluated with sampling based 2p piKL-SPARTA and sampling based SL policy. Additional evaluations with more $\lambda$ and more algorithm variants are presented in Table <ref> and Table <ref>. Error bar is 1 standard error.

Subset of Test Set     $\lambda=$0     $\lambda=$0.5     $\lambda=$1    $\lambda=$2     IL Policy

All games       25.30% $\pm$ 0.16%       56.04% $\pm$ 0.24%      60.22% $\pm$ 0.21%     62.44% $\pm$ 0.19%     63.63% $\pm$ 0.18%

Games w/score $\geq$ 10     27.01% $\pm$ 0.18%     58.86%  $\pm$ 0.25%     62.62% $\pm$ 0.22%     64.51% $\pm$ 0.20%     65.29% $\pm$ 0.19%

Games w/score $\geq$ 20     28.87% $\pm$ 0.19%     61.86% $\pm$ 0.25%     65.21% $\pm$ 0.23%     66.81% $\pm$ 0.21%     67.39% $\pm$ 0.19%

Subset of Test Set     $\lambda=$5     $\lambda=$10     $\lambda=$20     $\lambda=$50     IL Policy

All games       63.50% $\pm$ 0.18%     63.68% $\pm$ 0.18%     63.71% $\pm$ 0.18%     63.68% $\pm$ 0.18%     63.63% $\pm$ 0.18%

Games w/score $\geq$ 10     65.33% $\pm$ 0.19%     65.43% $\pm$ 0.19%     65.42% $\pm$ 0.19%     65.35% $\pm$ 0.19%     65.29% $\pm$ 0.19%

Games w/score $\geq$ 20     67.48% $\pm$ 0.19%     67.54% $\pm$ 0.19%     67.53% $\pm$ 0.19%     67.45% $\pm$ 0.19%     67.39% $\pm$ 0.19%

Human prediction accuracy of unregularized SPARTA ($\lambda = 0$) and of piKL-SPARTA with different $\lambda$ and the IL policy on test set. Each row represents their accuracy on a subset filtered by the final score of the games. pikl-SPARTA achieves similar prediction accuracy as IL for most $\lambda$ and is far more accurate than unregularized SPARTA.

$\lambda=$0     $\lambda=$0.5     $\lambda=$1     $\lambda=$2     $\lambda=$5     $\lambda=$10     IL Policy

1p piKL-SPARTA     20.56 $\pm$ 0.05     21.16 $\pm$ 0.06     20.02 $\pm$ 0.09     17.71 $\pm$ 0.12     13.53 $\pm$ 0.16     11.04 $\pm$ 0.17     2*8.81 $\pm$ 0.15

2p piKL-SPARTA     20.23 $\pm$ 0.04     23.11 $\pm$ 0.03     22.60 $\pm$ 0.03     21.63 $\pm$ 0.05     18.65 $\pm$ 0.11     15.23 $\pm$ 0.15

1p piKL-SPARTA-G     22.78 $\pm$ 0.03     23.07 $\pm$ 0.03     22.76 $\pm$ 0.03     22.41 $\pm$ 0.04     21.91 $\pm$ 0.05     21.45 $\pm$ 0.06     2*19.72 $\pm$ 0.10

2p piKL-SPARTA-G     19.98 $\pm$ 0.04     23.72 $\pm$ 0.02     23.39 $\pm$ 0.03     22.93 $\pm$ 0.03     22.36 $\pm$ 0.03     21.98 $\pm$ 0.04

Performance of unregularized SPARTA ($\lambda = 0$) and piKL-SPARTA under different $\lambda$ and IL policy evaluated on 4,000 self-play games, reported $\pm$ one standard error. In the top two rows, both piKL-SPARTA and IL policy sample actions according to their action distribution respectively. In the bottom two rows, both algorithms take greedy actions; “-G" is short for “-Greedy". All numbers shown in each table section use the same learned belief model for fair comparison. At the high lambdas that maintain or improve human accuracy, piKL-SPARTA significantly improves playing strength over IL, while at lower lambda values piKL-SPARTA outscores unregularized SPARTA. For reference, unregularized greedy 1p SPARTA with exact beliefs rather than learned beliefs gets 23.49 $\pm$ 0.02.

Table <ref> summarizes the results for human move prediction. On the full test set, piKL-SPARTA with $\lambda = 10$ and $\lambda = 20$ outperform the IL policy slightly, although not statistically confidently, and vastly outperform unregularized SPARTA. We also investigate the prediction accuracy on games with score $\geq 10$ or $\geq 20$ that more predominantly come from more experienced human players. The improvement of piKL-SPARTA over the IL policy may be larger on these games, although this result is noisy and at best should be considered only mildly suggestive since filtering on final score also introduces other major confounding factors as well.

In Table <ref>, we show the self-play performance piKL-SPARTA and IL policy. The top two rows show the results of the sampling version while the bottom two rows show those of the greedy version. The conclusions are consistent across both cases. Both 1p and 2p variants of piKL-SPARTA outscore IL for all $\lambda$ tested. Together with Table <ref> we find with $\lambda = 10$, piKL-SPARTA maintains IL prediction accuracy and outscores IL greatly in self-play, an overall improvement without any tradeoff. For smaller $\lambda = 2$ or $\lambda = 5$, piKL-SPARTA further outscores IL while losing some prediction accuracy on human moves, but remains vastly more accurate than unregularized SPARTA ($\lambda = 0$).

Through qualitative analysis of the games played by both methods, we find that IL makes many mistakes due to both sampling low probability actions and the fact that the training set contains bad moves. piKL-SPARTA, on the other hand, avoids many of these mistakes while still otherwise following the same strategies and humanlike signaling conventions. It is particularly good at preventing catastrophic failures where the agents lose all life tokens and points since the $Q$ values for good and bad actions differ much more widely in those cases.

The 2p piKL-SPARTA variants, despite being theoretically unsound, result in a further score improvement. Meanwhile, we notice that 2p versions of the plain SPARTA ($\lambda = 0$) underperform their 1p variants, which is consistent with the observations from [Lerer et al., 2020] despite us using an approximate learned belief model instead of exact belief. This suggests that regularization towards the blueprint IL policy is successful in keeping  the trajectories similar enough to those from the blueprint that the unsound assumption that the partner plays the blueprint does not cause major problems, while still allowing both players to deviate enough to correct major blunders that they would otherwise make.

Lastly, we notice that 1p piKL-SPARTA with $\lambda=0.5$ outperforms the unregularized version in self-play. The reason could be that the quality of the samples from the learned belief model worsen as the trajectories it sees at test time becomes more off-distribution for smaller $\lambda$. For reference, we also run the original SPARTA which does not have this problem as it computes beliefs analytically. The original SPARTA gets 23.49 $\pm$ 0.02, which is better than 1p piKL-SPARTA with $\lambda=0.5$ but worse than 2p piKL-SPARTA with the same $\lambda$. The computational cost of 1p SPARTA with exact beliefs and 2p piKL-SPARTA with approximate learned beliefs is roughly the same because the exact beliefs computation accounts for 50% of the entire computation. Therefore, 2p piKL-SPARTA could be a preferable option to improve performance in Dec-POMDPs via self play without incurring the huge cost of joint SPARTA or RL-search [Fickinger et al., 2021].


]PSMA PET as a predictive tool for sub-regional importance estimates in the parotid gland

, , , ,

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 Molecular Oncology, BC Cancer, Vancouver, BC, CA

6Department of Radiation Oncology, BC Cancer, Vancouver, BC, CA

7Department of Surgery, Faculty of Medicine, University of British Columbia, Vancouver, BC, CA

Objective: Xerostomia (subjective sensation of oral dryness) and radiation-induced salivary gland dysfunction remain a common side effect for head-and-neck radiotherapy patients, and attempts have been made to quantify the variation of the dose response within parotid glands. Here, we aim to compare several models of parotid gland regional importance with prostate specific membrane antigen (PSMA) positron emission tomography (PET), which has high concentrations of uptake in salivary glands and has been previously suggested to relate to gland functionality. Furthermore, we develop a predictive model of Clark et al.'s relative importance using radiomic features, and demonstrate a methodology for predicting patient-specific importance deviations from the population. Approach: Intra-parotid uptake was compared with four regional importance models using [18F]DCFPyL PSMA PET images. The correlation of uptake and importance was ascertained when numerous non-overlapping sub-regions were defined, while a paired t-test was used when binary regions were defined. Radiomic PSMA PET/CT features for Clark et al.'s sub-regions were used to develop a predictive model of population importance using a double cross-validation methodology. We demonstrate a method for supplementing population importance estimates using individual patient features.
Main Results: Clark et al.'s relative importance regions were significantly ($p < 0.02$) anti-correlated with PSMA PET uptake. Van Luijk et al.'s critical regions had significantly lower ($p < 0.01$) uptake than in non-critical regions. Kernel Ridge Regression with principal component analysis feature selection performed best over test sets (Mean Absolute Error = 0.08), with gray level co-occurrence matrix (GLCM) features being particularly important. Deblurring PSMA PET images with neural blind deconvolution strengthened correlations and improved model performance. Significance: This study suggests that regions of relatively low PSMA PET concentration in parotid glands may exhibit relatively high dose-sensitivity. We've demonstrated the ability of PSMA PET radiomic features for predicting relative importance within the parotid glands.

§ INTRODUCTION:

Intensity Modulated Radiotherapy (IMRT) allows for the creation of treatment plans with high dose conformity in cancerous regions while minimizing dose to healthy tissue [1]. However, high dose levels in cancerous tissue inevitably tail off into healthy tissue, so treatment planners and oncologists must prioritize which healthy regions to spare. Treatment planning in the head-and-neck region is particularly challenging, as there are many organs in close proximity which often abut or overlap with tumour volumes [2]. Dose levels in the salivary glands are of particular concern, as xerostomia (subjective sensation of oral dryness) remains a common side effect for head-and-neck cancer patients [3]. Dose to the largest salivary glands, the parotids, is the greatest risk factor for post-treatment xerostomia [4].

The current standard of care is to minimize the whole-mean dose to parotid glands [5], which were previously considered to have a uniform dose response [6]. However, there have been numerous attempts in recent years to quantify the relative importance of various parotid gland sub-regions for predicting post-treatment complications [7, 8, 9, 10].

Prostate specific membrane antigen (PSMA) positron emission tomography (PET) has high ligand accumulation in the parotid glands [11, 12, 13], and has been suggested to relate to whole-gland functionality [14, 15, 16]. Furthermore, uptake within parotid glands has been found to be non-uniform, with high uptake regions tending towards lateral, posterior, and superior regions [17]. We hypothesize that intra-parotid gland uptake variability of PSMA PET is predictive of functional importance.

The purpose of this study is two-fold. First, we use a data set of 30 PSMA PET images to compare intra-parotid PSMA PET uptake trends with several regional importance estimates from the literature [8, 7, 9, 10]. Second, we develop a population-level model of Clark et al.'s [8] regional importance using radiomic features from  PSMA PET and Computed Tomography (CT). We also demonstrate how such a model can be used to predict a patient's deviation from population-derived importance estimates, creating a single metric with population-derived, and patient-specific components.

§ METHODS:

§.§ 2.1 Dataset

This study was approved by an institutional review board. The data set included de-identified [18F]DCFPyL PSMA PET/CT images for 30 previous prostate cancer patients (Mean Age 68, Age Range 45-81; mean weight: 90 kg, weight range 52 - 128 kg). Scans were acquired 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.16 mm, slice thickness: 2.8 - 3.02 mm).  Helical CT scans were acquired on the same scanner (kVP: 120, pixel spacing: 0.98 mm, slice thickness: 3.75 mm). Images were scaled to standard uptake values normalized by lean body mass ($SUV_{lbm}$). Registered CT images were used for delineating parotid and submandibular glands. Limbus AI [18] was used for preliminary auto-segmentation of the glands, which were then manually refined by a single senior radiation oncologist, Jonn Wu.

§.§ 2.2 Correction of Partial Volume Effects

One weakness of PET as an imaging modality is its intrinsically low spatial resolution [19]. The burden of partial volume effects is less pronounced when analyzing large geometric regions with homogeneous uptake, but cannot be ignored when attempting to compare heterogeneous uptake in small regions-of-interest (ROIs), such as sub-regions of the parotid glands. Recently, a method has been developed for simultaneous deblurring and super-sampling of PSMA PET images using neural blind deconvolution [20], which we employ in this work for pre-processing PSMA PET images. This model has been shown to illuminate fine uptake trends within small regions of parotid glands. We performed all calculations with both "enhanced" images, and unmodified original images.

§.§ 2.3 Comparison of PSMA PET with Parotid Gland Importance Models

PSMA PET uptake trends were compared with four models of intra-parotid gland importance found in the literature, detailed in the following sub-sections. Uptake metrics included the mean, median, and maximum, calculated in ROIs defined according to each specific model. As importance models were all population-level estimates, uptake metrics were averaged over the 60 parotid glands from the 30 patients.

§.§.§ Clark et al.'s Model

Clark et al.'s model [8] estimates the relative importance of 18 equal-volume sub-regions of the contralateral parotid gland for predicting salivary dysfunction following radiotherapy. Stimulated saliva measurements were collected for 332 patients before and at one year after radiotherapy. The relative differences were predicted with conditional inference trees using radiotherapeutic dose levels in parotid gland sub-regions. Parotid glands were sub-segmented using nested planar segmentation (planes: 2 axial, 1 coronal, 2 sagittal). Regions of high relative importance tended towards caudal-anterior regions, as shown in Fig <ref>.

3D rendering of voxels within a parotid gland corresponding to Clark et al's relative importance sub-regions are shown from two angles.

We sub-segmented parotid glands according to the same regimen, and voxels within each of the 18 sub-regions were used to calculate uptake statistics. To test whether intra-parotid PSMA PET uptake is related to Clark et al.'s importance estimates, Spearman's rank correlation coefficient, $r_s$, was computed between uptake in the 18 sub-regions along with their corresponding importance estimates.

§.§.§ Han et al.'s model

Han et al. [7] assess the relative importance of 9 parotid gland sub-regions for predicting injury ($\ge$ grade 2 xerostomia at 6 months post-radiotherapy) and recovery ($\ge$ grade 2 xerostomia at 6 months post-radiotherapy, followed by $<$ grade 2 xerostomia at 18 months post-radiotherapy). Sub-regions were defined by first applying a 3 mm margin to whole parotid glands, then dividing glands into three radial sectors (anterior, medial, posterior), then further dividing these sectors along the inferior-superior axis into 3 equal-length regions. Voxels within a parotid gland corresponding to these importance regions are shown in Fig <ref>.

3D rendering of voxels within a parotid gland corresponding to Han et al.'s relative importance sub-regions are shown from two angles.

Han et al. [7] determined the relative importance of 9 dose-volume statistics in 10% volume increments from D10 (Minimum dose to 10% volume) to D90 (Minimum dose to 90% volume) in each sub-region. For our purposes, the mean importance computed over all dose statistics was used as a single relative importance estimate for each sub-region. Spearman's rank correlation coefficient, $r_s$, was calculated between uptake and relative importance for predicting injury, and recovery.

§.§.§ Van Luijk et al.'s Model

Van Luijk et al. [9] used stimulated saliva measurements and radiotherapeutic dose levels to locate "critical" regions within parotid glands which are most predictive of salivary outcome at one year post radiotherapy. The study did not specify a well-defined location of this critical region over the population, however, it is stated to be in close proximity of the Stensen's duct, adjacent to the dorsal side of the mandible. For our purposes, we approximated the critical region by applying a 9 mm margin to the mandible, which was intersected with the top half of the parotid gland (Fig <ref>). The 9 mm margin was found to consistently intersect a region of the parotid gland  approximately corresponding Van Luiijk et al.'s [9] critical regions.

The approximate location of Van Luijk et al.'s critical region of the parotid gland used for computing uptake statistics is shown

Uptake statistics were compared within expanded critical and non-critical regions using a paired t-test.

§.§.§ Buettner et al.'s Model

Buettner et al. [10] evaluated the predictive ability of various dose "moments" in a regression model for post-treatment xerostomia in 63 head-and-neck cancer patients, treated with either IMRT or conventional radiotherapy. Important variables included mean dose to the superficial lobe, skewness of dose in the cranial-caudal direction within the deep and superficial lobe, and relative concentration of dose in the caudal-medial region of the deep lobe. While the parotid glands were only segmented into superficial and deep lobes, dose moments calculated within these regions evaluated the spatial variance of the dose response.

§.§ 2.4 Development of a predictive Model for parotid gland relative importance using PSMA PET and CT

To demonstrate the predictive ability of PSMA PET for parotid gland regional functionality, we develop a model for predicting Clark et al.'s [8] relative importance estimates using radiomic features extracted from both PSMA PET and CT images. We then demonstrate how such a model can be used for predicting patient-specific perturbations away from population-level importance estimates.

§.§.§ Feature Extraction

The standardized pyradiomics library [21] was used for computing radiomic features of PSMA PET and CT within all 18 of Clark's sub-regions. The full set of Gray Level Co-occurrence Matrix (GLCM), Gray Level Run Length Matrix (GLRLM), Gray Level Size Zone Matrix (GLSZM), Gray Level Dependence Matrix (GLDM), and first order features, computed with original squared, square root, and wavelet image types, were computed for a total of 1060 features. For gray level discretization, a fixed bin width was chosen over a fixed bin count, as a fixed bin width has been shown to have better reproducibility [22], especially when chosen to yield a bin count between 16-128 [23]. We therefore set the bin count to 0.2 for original, 0.1 for square root, and 1 for square. For the wavelet features, a fixed bin count of 100 was used, due to uncertainty in the expected range.

Radiomic features were calculated for individual patients, and averaged over each parotid gland for all 18 patients. This yielded a population-level design matrix of shape (36, 1060) prior to feature selection. Features were calculated for both enhanced and original PSMA PET images.

§.§.§ Double Cross Validation

The small size of our data set made it inappropriate to define a single test set for final performance evaluation, and we therefore used double cross validation, or sometimes called nested cross validation, where our test set was rotated through 9 outside folds, each having its own inner cross validation loop for tuning the feature selection algorithm, model, and hyper-parameters (Fig <ref>).

For model testing and validation, a double cross validation scheme was employed, where the outside test set is rotated through a 9-fold cross validation loop, each including its own 8-fold inner cross validation loop for parameter tuning.

§.§.§ Feature Selection and Models

To avoid overfitting, the large number of features extracted must be pruned using feature selection methods prior to model training. For this purpose, we include 3 feature selection algorithms within the cross validation loops, including a linear combination filter, a pairwise correlation filter, and principle component analysis as used by Delzell et al. [24] to predict lung cancer using radiomic features.  The linear combination (lincom) filter uses QR decomposition to iteratively remove features which are linear combinations of others. The pairwise correlation filter tests the correlation between features and removes those who are correlated above a specified cutoff. Principal component analysis changes the basis of the feature space to capture a large portion of the variance using a smaller number of feature vectors. For more information, refer to the work by Delzell et al. [24]

Five different regression model types were included within cross validation. This included a linear regression, support vector machine, random forest, conditional inference tree, and kernel ridge model. Model performance is highly dependent on a variety of hyper-parameters which can be tuned for the different models and feature selection algorithms. A set of different hyper-parameters for each model were iterated over within cross validation. All models, feature selection algorithms, and their corresponding hyper-parameters tested are listed in Table <ref> Models were scored according to the mean absolute error (MAE).

Models and feature selection (F.S) algorithms, along with their corresponding hyper-parameters tested in cross-validation, are shown.

Dark_GrayModel   Hyper-parameters

GraySupport Vector Machine   $\epsilon = [0.01,0.05,0.1]$

Gray   Kernel = [Linear, radial basis function, sigmoid, poly]

Gray   Degree = [2,3]

Gray   $\gamma = $ [scale, auto]

Gray   coef0 = $[-0.3, -0.2, -0.1, 0]$

Random Forest   number estimators $= [3,5,7,10]$

max depth $= [3,5,8, None]$

Criterion = [Absolute Error, Squared Error]

GrayConditional Inference Trees   Max Depth $= [5,10,15,20,25,None]$

Gray   Criterion = [Absolute Error, Squared Error]

Kernel Ridge   $\alpha = [0.1, 0.5, 1, 5, 10]$

Kernel = [Linear, Radial Basis Function, Sigmoid]

GrayLinear Regression   N/A

Dark_GrayF.S. Algorithm   Hyper-parameters

Pairwise Correlation Filter   feature count = [1,2,3,4,5,6], cutoff = [0.85, 0.88, 0.9, 0.92]

GrayPCA   feature count = [1,2,3,4,5,6,8,10,15,20,30]

Linear Combination Filter   correlation cutoff = [0.05, 0.1, 0.2, 0.3]

§.§.§ Error Analysis

For estimating the uncertainty of model predictions, we employ the methodology described by Cawley et al. [25] for kernel ridge models. This involves computing the leave-one-out absolute error for each sample in the data set, and then training a second kernel ridge model for predicting the absolute error of predictions based on the same input features and using this as an estimate of prediction variance.

§.§.§ Demonstrating a method for predicting patient-specific importance perturbations

Finally, we demonstrate how patient-specific deviations from population-level importance estimates can be obtained, to create a single, combined importance estimate including both population-level and patient-specific components. This is obtained by first computing and processing a patient's radiomic features according to the feature selection algorithm employed by the final model. As the model has been trained to understand the relationship between specific features and relative importance estimates, inputting patient specific features into the model for all 18 parotid gland sub-regions provides an estimate of relative importance for said patient. A single combined importance estimate is created by taking the population level sub-region importance estimates, $I^P_j$, $j \in \mathbb{Z}, 1 \leq j \leq 18$ and the difference in patient specific and population estimates, $\Delta_j$, $j \in \mathbb{Z}, 1 \leq j \leq 18$ and computing

\begin{equation}
I = \begin{cases}
I^P_j & \Delta_j < 0 \\
\frac{2I^P_j}{1+e^{-2\Delta_j}} &\Delta_j > 0 \\
\end{cases}
\end{equation}

This defines a minimum importance estimate using the population-level estimates, and increases estimates in regions of high patient-specific importance, levelling off as the patient-specific estimate approaches about 3x the maximum population-level importance estimate. Using this approach, patient-specific predictions can be used to supplement or perturb population-level estimates in regions predicted to be of high radiotherapeutic importance. Combined importance estimates for sub-regions are never lower than population-level estimates, to avoid potential negative impacts associated with under-estimating importance.

§ RESULTS:

§.§ 3.1 Comparison of PSMA PET with Importance Models

Overall, uptake of PSMA PET was found to be inversely proportional with sub-region importance estimates from the literature. These trends appeared stronger when enhanced images were used for uptake calculations.

Clark et al.'s [8] importance predictions in the 18 equal volume regions were significantly anti-correlated with mean and median uptake (Table <ref>). A scatter plot of importance vs mean uptake is shown in Fig <ref>.

Spearman's rank correlation coefficients for PSMA PET uptake with Clark's relative importance estimates. Correlations are calculated for mean, median, and maximum uptake, normalized by lean body mass. Results are calculated using enhanced (de-blurred and super-sampled) and original PSMA PET images.

Mean   Median   Maximum

Enhanced   $r_s = -0.56, p = 0.015$   $r_s = -0.55, p = 0.016$   $r_s = -0.25, p = 0.31$

Original   $r_s = -0.50, p = 0.03$   $r_s = -0.51, p = 0.03$   $r_s = -0.30, p = 0.22$

Clark's relative importance vs mean PSMA PET uptake in 18 equal-volume parotid gland sub-regions, averaged over 30 patients. Relative importance was found to have a significant ($p=0.015$) anti-correlation with regional PSMA PET uptake. Calculations were performed with de-blurred PSMA PET images. A best fit line is shown in red.

Han et al's [7] model predictions for relative importance in 9 regions of unequal volume were not significantly correlated with uptake metrics. There did however exist weak trends of anti-correlation of uptake with importance for predicting injury, and direct correlation for predicting recovery <ref>.
Uptake levels were found to be approximately two times higher in Van-Luijk et al.'s [9] non-critical region than in the critical region ($p < 0.01$) (Table <ref>).

Spearman's rank correlation coefficients for PSMA PET uptake with Han's relative importance estimates for predicting both injury, and recovery. Correlations are calculated for mean, median, and maximum uptake, normalized by lean body mass, using enhanced (de-blurred and super-sampled) and original PSMA PET images.

2c$r_s, p$ (Mean & Imp)   2c$r_s, p$ (Median & Imp)   2c$r_s, p$ (Maximum & Imp)

Injury   Recovery   Injury   Recovery   Injury   Recovery

Enhanced   -0.27, 0.49   0.39, 0.29   -0.07, 0.86   0.12, 0.76   -0.56, 0.11   0.49, 0.18

Original   -0.41, 0.26   0.42, 0.26   -0.19, 0.61   0.18, 0.64   -0.56, 0.11   0.49, 0.18

Comparisons for mean, median, and maximum uptake, normalized by lean body mass, in Van Luijk et al.'s [9] critical and non-critical parotid gland sub-regions. Calculations are shown using enhanced (de-blurred and super-sampled) and original PSMA PET images.

2cMean   2cMedian   2cMaximum

Critical   Non-Critical   Critical   Non-Critical   Critical   Non-Critical

Enhanced   $3.1 \pm 2.1$   $7.2 \pm 2.1$   $2.7 \pm 2.1$   $7.2 \pm 2.2$   $9.4 \pm 5.1$   $20.0 \pm 6.0$

Original   $3.2 \pm 2.0$   $6.3 \pm 1.7$   $2.9 \pm 2.1$   $6.2 \pm 1.8$   $7.6 \pm 4.0$   $15.6 \pm 4.5$

Uptake statistics in parotid gland sub-regions corresponding to Buettner et al.'s important regions [10] are shown in Table <ref>. As shown previously, [17], uptake in parotid glands appears skewed towards lateral, superior regions. Buettner et al. found the superficial and the deep lobe to be important, with importance concentrating slightly towards caudal regions. The caudal-medial sub-region of the deep lobe, which was predicted to be of high importance, was found to have significantly lower (p < 0.01) uptake than in the caudal-lateral sub-region. Caudal regions, in general, tend to have lower uptake than superior regions of the glands.

Buettner et al. [10] found that dose to the superficial lobe, relative concentration of dose in the caudal/cranial region of both superficial and deep lobes, and relative concentration in the caudal-medial region of the deep lobe, were predictive of post-treatment xerostomia for head-and-neck radiotherapy patients. It is unclear whether xerostomia is directly or inversely proportional to these metrics, so we simply report differences in PSMA uptake within corresponding regions. Correlations are calculated for mean, median, and maximum uptake, normalized by lean body mass, using enhanced (de-blurred and super-sampled) and original PSMA PET images.

white      Sup   Deep   Sup Cranial   Sup Caudal   Deep Cranial   Deep Caudal   Deep Caudal-Medial   Deep Caudal-Lateral

Enhanced   Mean   $7.4 \pm 2.2$   $5.7 \pm 2.0$   $8.4 \pm 2.7$   $7.2 \pm 2.1$   $6.2 \pm 2.0$   $5.1 \pm 2.2$   $4.2 \pm 2.3$   $5.9 \pm 2.6$

Medial   $7.6 \pm 2.4$   $5.3 \pm 2.3$   $8.8 \pm 2.8$   $7.0 \pm 2.2$   $5.8 \pm 2.3$   $4.6 \pm 2.6$   $5.7 \pm 2.6$   $5.7 \pm 3.0$

Maximum   $19.8 \pm 6.0$   $16.3 \pm 4.7$   $18.8 \pm 5.6$   $18.3 \pm 5.8$   $15.7 \pm 4.6$   $13.8 \pm 4.6$   $14.8 \pm 4.7$   $15.2 \pm 4.2$

Original   Mean   $6.5 \pm 1.9$     $5.2 \pm 1.7$     $8.0 \pm 2.5$     $6.3 \pm 1.8$     $5.6 \pm 1.9$     $4.5 \pm 1.9$     $3.5 \pm 1.9$   $5.4 \pm 2.3$

Median   $6.5 \pm 1.9$     $5.1 \pm 1.9$     $8.4 \pm 2.7$     $6.3 \pm 1.9$     $5.6 \pm 2.1$     $4.2 \pm 2.1$     $3.1 \pm 2.0$     $5.4 \pm 2.7$

Maximum   $10.9 \pm 5.6$   $8.9 \pm 4.5$   $14.8 \pm 4.0$     $14.1 \pm 4.0$     $11.9 \pm 3.2$   $10.5 \pm 3.5$   $8.1 \pm 3.6$   $10.2 \pm 3.6$

§.§ 3.2 Model Performance for Predicting Parotid Gland Relative Importance

For each of the 9 test-sets of the outer cross validation loop, the M.A.E, along with best model and feature selection algorithm, as determined via inner cross validation, are shown in Table <ref> The average M.A.E was 0.08 using enhanced images, and 0.15 with original images. Overall, the best performing model and feature selection algorithm was the kernel ridge regressor with principal component analysis using 20 features. A performance comparison between all models and feature selection algorithms is shown in Fig <ref>. The most important features for principal component analysis (determined by projecting principal components scaled by their singular values onto original feature axes), are shown in Table <ref>. The overall-best hyper-parameters were found to be a polynomial kernel of degree 2, with $\alpha = 0.1$, and coef0$=1$.

Mean absolute error for each test set of the outer cross validation is shown, along with the best performing model and feature selection algorithm determined during the inner cross validation. Results are shown for models created with both enhanced and original PSMA PET images.

white   3cEnhanced Images   3cOriginal Images

Fold   M.A.E   Model   F.S Algorithm   M.A.E   Model   F.S Algorithm

1   0.08   K.R   P.C.A   0.31   C.I.T   P.C.A

2   0.10   K.R   P.C.A   0.29   R.F   P.W.C

3   0.06   K.R   P.C.A   0.08   K.R   P.C.A

4   0.04   K.R   P.C.A   0.11   K.R   P.C.A

5   0.06   C.I.T   P.C.A   0.24   R.F   P.C.A

6   0.09   K.R   P.C.A   0.03   C.I.T   P.C.A

7   0.07   K.R   P.C.A   0.02   C.I.T   P.C.A

8   0.13   C.I.T   P.C.A   0.04   K.R   P.C.A

9   0.08   K.R   P.C.A   0.20   C.I.T   P.C.A

whiteAverage   0.08         0.15

Relative importance of radiomic PSMA PET / CT features determined via principal components analysis for modelling of parotid gland sub-region importance.

Feature Name   Modality   Relative Importance

1.   Original GLCM - Inverse Difference   PSMA PET   1.0

2.   Square Root GLCM - Run Variance   PSMA PET   0.98

3.   Original GLCM - Inverse Difference Moment   PSMA PET   0.97

4.   Square Root GLCM Inverse Difference   PSMA PET   0.97

5.   Original GLCM Inverse Variance   PSMA PET   0.96

6.   Square Root GLCM Inverse Difference Moment   PSMA PET   0.95

7.   Square Root First Order RMS   CT   0.95

8.   Square GLSZM GLNUN   PSMA PET   0.94

9.   Original GLRLM Long Run Emphasis   PSMA PET   0.94

10.   Original GLDM Large Dependence   PSMA PET   0.93

11.   Original GLRLM Run Variance   PSMA PET   0.91

12.   Wavelet HHH GLCM Joint Average   CT   0.91

13.   Wavelet GLCM Sum Average   CT   0.91

14.   Wavelet HHH GLSZM High Gray Level Zone Emphasis   CT   0.91

15.   Square Root GLCM Inverse Variance   CT   0.90

16.   Wavelet HHH GLDM High Gray Level Emphasis   CT   0.90

17.   Wavelet HHH GLRLM High Gray Level Run Emphasis   CT   0.90

18.   Wavelet LLL GLCM Dependence Entropy   CT   0.90

19.   Wavelet HHH GLCM Autocorrelation   CT   0.89

20.   Square Root First Order Entropy   CT   0.88

To compare performance across models and feature selection algorithms, the mean absolute error (MAE) of the top-performing model and feature selection algorithm  from each test fold was computed and averaged over all folds. Overall, kernel ridge regression and principal component analysis using 20 features demonstrated the best performance. Error bars correspond to the standard deviation of prediction accuracy over various hyper-parameters.

All model predictions for importance in the 9 test sets are collected and plotted against Clark's estimates in Fig <ref>. Prediction error was estimated by a separate kernel ridge model which was trained to predict model estimation error, as previously described.

Model-predicted relative importance estimates for all sub-regions are plotted along with Clark et al.'s (in this case, ground truth) importance estimates. Predictions are shown with their associated model predicted errors. Predictions were made using the best performing model and feature selection algorithm found with nested cross validation - a kernel ridge regressor model and principal component analysis feature selection using 20 features.

§.§ Estimating perturbations of patient importance from the population

The best performing model, feature selection algorithm, and hyper-parameters as determined via cross validation, were then used to train a final population-level model using the entire population-level data set. Importance estimates for parotid gland sub-regions of individual patients could then be obtained by inputting a given patient's radiomic features into the model. Examples of how individual predictions deviate from population-level estimates, along with prediction errors estimated with the kernel ridge error model, are demonstrated for six different patients in Fig <ref>

Importance estimates obtained for six individual patients using the population-level model, are shown. Population level estimates for the 18 sub-regions are shown as purple squares, with the patient specific estimates as gold circles. The population-level model captures the relationship between important radiomic features and importance estimates, and can be used for estimating approximate shifts in importance estimates for an individual patient's parotid gland sub-regions. Error estimates were obtained with the kernel ridge error model.

Fig <ref> illustrates an example of how an individual patient's parotid gland radiomic features can be used to supplement the population importance estimate, using the equation described in the methods.

We demonstrate how patient-specific parotid sub-regional importance estimates (left), can be used to supplement population-level importance estimates (middle) using the formula described in the methods, such that final estimates are never lower than population-level estimates, but further increased in regions where patient-specific estimates are high. Only positive perturbations to regional importance were made, to avoid negatively impacting patients in the case where importance estimates are to be used for designing dose constraints.

§ DISCUSSION:

The results of this work demonstrate that intra-parotid PSMA PET uptake may be inversely related to regional importance. Clark et al.'s [8] and Han et al.'s [7] relative importance estimates were practically advantageous in that they define numerous non-overlapping sub-regions of parotid glands where correlations between relative importance and PSMA PET uptake could be assessed. Both models predict higher importance towards medial and caudal-middle regions of the gland, while Clark et al.'s predicts much higher importance in the anterior half of the gland. Clark et al.'s model was significantly ($p < 0.02$) anti-correlated with PSMA PET uptake, while Han et al.'s model was also anti-correlated, but not significantly. It should be noted, however, that regions of high importance predicted by both models correspond to regions of lower than average uptake [17]. Han et al.'s sub-segmentation method yields sub-regions of unequal volume, which also creates problems when comparing uptake statistics. Due to the shape of the parotid glands, sub-regions within center regions will be larger than those at the top and especially the bottom, when sub-segmenting according to equal superior-inferior length. Note that both models (Clark et al. and Han et al.) were developed using salivary measurements. The anti-correlation relationship of importance with PSMA PET uptake is unexpected and suggests an underlying connection that could relate to salivary gland functionality.

We were able to approximate the location of Van Luijk et al.'s [9] critical regions and compute PSMA PET uptake statistics within and outside said regions. Uptake statistics within critical regions were significantly lower than non-critical regions, supporting the anti-correlation trend between importance and uptake seen with Clark et al. and Han et al.'s sub-regions. Sparing dose in Van Luijk' et al.'s [9] critical regions of parotid glands for radiotherapy patients was recently shown to insignificantly impact patient outcomes [26]. However, dose to critical regions was more predictive of salivary dysfunction than whole-gland dose. Comparing uptake with importance in regions defined by Buettner et al.'s [10] analysis also pointed towards an anti-correlative relation of importance with PSMA PET uptake.

Simultaneous deblurring and super-sampling of PSMA PET images prior to uptake calculations led to stronger correlations between uptake and importance estimates, and better model performance for predicting regional importance. Better performance using enhanced images was expected, as partial volume effects cause fine detail wash-out in small regions of PET images.

The PSMA ligand binds to the PSMA epitope of acinar and ductal cells in salivary glands [14], so uptake has been suggested to be directly proportional to functional importance [14, 15, 16]. In addition, irradiation of parotid glands has been shown to decrease PSMA PET uptake [27]. It is likely that whole-gland PSMA PET uptake is directly proportional to whole-gland importance, but the results of this work demonstrate that intra-parotid functional importance may tend towards relatively low-uptake regions. It is unclear which physiological mechanisms within the gland could result in this relationship.

The relationship between PSMA PET and relative importance appears non-linear (Fig <ref>) and can be captured using radiomic features and non-linear modelling. Model development for predicting relative importance with PSMA PET and CT radiomic features was successful, yielding a relatively low absolute error for test predictions, considering the small size of the data set. Multilayer perceptrons were original included in the cross validation but were found to perform poorly due to over-fitting. Kernel ridge regression with a quadratic kernel out-performed all other models tested. This suggests that the relationship between PSMA PET uptake and regional importance is not simply linear. This is further supported by the most important features determined by principal component analysis (Table <ref>). In particular, radiomic features of squared PSMA PET uptake were highly predictive of importance. The top 6 most important features were all forms of the GLCM of PSMA PET images. Based on this finding, we recommend using them in future predictive models of functional importance for salivary glands.

In this work we demonstrated how PSMA PET and CT can be used together to predict a patient's deviation from population-level estimates of parotid gland regional importance. The purpose was to present a hypothetical method of extracting patient-specific parotid gland importance estimates for tailoring patient dose constraints for radiotherapy [28]. PSMA PET is not acquired as standard-of-care for head-and-neck patients, and is too likely too costly to add to the standard-of-care, and therefore these images would not be practically available in most clinical situations. However, we believe further PSMA PET studies could be critical in shedding light on importance trends and variability within the glands. The method presented here is not exclusive to PSMA PET models, and could be used to supplement other population-based models.

Our data set was small, comprised of only 60 parotid glands from 30 patients. This necessitated the double cross validation methodology used for model development, where the test set was rotated through, along with inner validation sets for each, to determine model parameters. Outer test sets had no influence on model development and were used to independently test the models predictive accuracy of all sub-regions. Clark et al.'s [8] importance regions are biased towards low numbers, with only a few regions of high importance. We believe double cross-validation is warranted in this case to help mitigate biased error estimates for different importance values.

§ CONCLUSION

In this work, we have compared four models of parotid gland sub-regional importance from the literature, with regional PSMA PET uptake. In general, an inverse proportionality between importance and PSMA PET uptake was observed. We demonstrated the ability of PSMA PET radiomic features to predict regional importance by training a predictive model of Clark et al.'s importance estimates (MAE = 0.08). Lastly, we demonstrated a methodology for supplementing population-level importance estimates using patient-specific radiomic features.

§ ACKNOWLEDGEMENTS

This work was supported by the Canadian Institutes of Health Research (CIHR) Project Grant.

§ REFERENCES

[1]
A. Taylor and M. E. B. Powell, “Intensity-modulated radiotherapy–what is
it?,” Cancer Imaging, vol. 4, pp. 68–73, Mar. 2004.

[2]
S. C. H. Cheng, V. W. C. Wu, D. L. W. Kwong, and M. T. C. Ying, “Assessment of
post-radiotherapy salivary glands,” Br. J. Radiol., vol. 84,
pp. 393–402, May 2011.

[3]
S. J. Ma, C. I. Rivers, L. M. Serra, and A. K. Singh, “Long-term outcomes of
interventions for radiation-induced xerostomia: A review,” World J.
Clin. Oncol., vol. 10, pp. 1–13, Jan. 2019.

[4]
M. S. Chambers, D. I. Rosenthal, and R. S. Weber, “Radiation-induced
xerostomia,” Head & Neck, vol. 29, no. 1, pp. 58–63, 2007.

[5]
J. O. Deasy, V. Moiseenko, L. Marks, K. S. C. Chao, J. Nam, and A. Eisbruch,
“Radiotherapy dose-volume effects on salivary gland function,” Int. J.
Radiat. Oncol. Biol. Phys., vol. 76, pp. S58–63, Mar. 2010.

[6]
A. Eisbruch, J. A. Ship, H. M. Kim, and R. K. Ten Haken, “Partial irradiation
of the parotid gland,” Semin. Radiat. Oncol., vol. 11, pp. 234–239,
July 2001.

[7]
P. Han, P. Lakshminarayanan, W. Jiang, I. Shpitser, X. Hui, S. H. Lee,
Z. Cheng, Y. Guo, R. H. Taylor, S. A. Siddiqui, M. Bowers, K. Sheikh,
A. Kiess, B. R. Page, J. Lee, H. Quon, and T. R. McNutt, “Dose/Volume
histogram patterns in salivary gland subvolumes influence xerostomia injury
and recovery,” Sci. Rep., vol. 9, p. 3616, Mar. 2019.

[8]
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.

[9]
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.

[10]
F. Buettner, A. Miah, S. Gulliford, E. Hall, K. Harrington, S. Webb,
M. Partridge, and M. Nutting, “Novel approaches to improve the therapeutic
index of head and neck radiotherapy: an analysis of data from the parsport
randomised phase iii trial,” Radiother. Oncol., vol. 103, p. 82–87,

[11]
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.

[12]
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.

[13]
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.

[14]
T. J. W. Klein Nulent, M. H. Valstar, B. de Keizer, S. M. Willems, L. A. Smit,
A. Al-Mamgani, L. E. Smeele, R. J. J. van Es, R. de Bree, and W. V. Vogel,
“Physiologic distribution of PSMA-ligand in salivary glands and seromucous
glands of the head and neck on PET/CT,” Oral Surg. Oral Med. Oral
Pathol. Oral Radiol., vol. 125, pp. 478–486, May 2018.

[15]
Y. Zhao, Y. Xia, H. Liu, Z. Wang, Y. Chen, and W. Zhang, “Potential
applications of 68Ga-PSMA-11 PET/CT in the evaluation of salivary gland
uptake function: Preliminary observations and comparison with 99mTcO4 -
salivary gland scintigraphy,” Contrast Media Mol. Imaging, vol. 2020,
p. 1097516, Jan. 2020.

[16]
V. Mohan, W. V. Vogel, G. D. Valk, J. P. de Boer, M. G. E. H. Lam, and
B. de Keizer, “PSMA PET/CT identifies intrapatient variation in salivary
gland toxicity from iodine-131 therapy,” Mol. Imaging, vol. 19,
p. 1536012120934992, Jan. 2020.

[17]
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.” under journal review, 2023.

[18]
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.

[19]
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.

[20]
S. Caleb, U. Carlos, A. Rahmim, F. Benard, J. Wu, and H. Clark, “Neural blind
deconvolution for simultaneous partial volume effect correction and
super-sampling of psma pet images,” arXiv: arXiv:2309.00590v1, 2023.

[21]
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.

[22]
Leijenaar RTH, Nalbantov G, Carvalho S, van Elmpt WJC, Troost EGC, Boellaard R,
et al. ; The effect of SUV discretization in quantitative FDG-PET Radiomics:
the need for standardized methodology in tumor texture analysis; Sci Rep.

[23]
F. Tixier, C. C. Le Rest, M. Hatt, N. Albarghach, O. Pradier, J.-P. Metges,
L. Corcos, and D. Visvikis, “Intratumor heterogeneity characterized by
textural features on baseline 18F-FDG PET images predicts response to
concomitant radiochemotherapy in esophageal cancer,” J. Nucl. Med.,
vol. 52, pp. 369–378, Mar. 2011.

[24]
D. A. P. Delzell, S. Magnuson, T. Peter, M. Smith, and B. J. Smith, “Machine
learning and feature selection methods for disease classification with
application to lung cancer screening image data,” Frontiers in
Oncology, vol. 9, 2019.

[25]
G. C. Cawley, N. L. C. Talbot, and O. Chapelle, “Estimating predictive
variances with kernel ridge regression,” in Machine Learning
Challenges. Evaluating Predictive Uncertainty, Visual Object Classification,
and Recognising Tectual Entailment (J. Quiñonero-Candela, I. Dagan,
B. Magnini, and F. d'Alché Buc, eds.), (Berlin, Heidelberg), pp. 56–77,
Springer Berlin Heidelberg, 2006.

[26]
R. J. H. M. Steenbakkers, M. I. van Rijn-Dekker, M. A. Stokman, R. G. J.
Kierkels, A. van der Schaaf, J. G. M. van den Hoek, H. P. Bijl, M. C. A.
Kramer, R. P. Coppes, J. A. Langendijk, and P. van Luijk, “Parotid gland
stem cell sparing radiation therapy for patients with head and neck cancer: A
double-blind randomized controlled trial,” Int. J. Radiat. Oncol. Biol.
Phys., vol. 112, pp. 306–316, Feb. 2022.

[27]
M. Hotta, K. Nguyen, P. Thin, W. R. Armstrong, A. Gafita, M. R. Benz,
J. Czernin, A. U. Kishan, N. G. Nickols, and J. Calais, “Kinetics of psma
pet uptake in prostate cancer lesions after radiation therapy: A
single-center retrospective study.,” Journal of Clinical Oncology,
vol. 40, no. 6_suppl, pp. 36–36, 2022.

[28]
C. Sample, J. Wu, S. Thomas, and H. Clark, “Incorporating parotid gland
inhomogeneity into head-and-neck treatment optimization through the use of
artificial base plans,” J Appl Clin Med Phys, vol. 22, pp. 141–149,
Mar 2021.