#  Image denoising and model-independent parameterization for improving IVIM
MRI

Caleb Sample1,2, Jonn Wu3,4, Haley Clark1,2,4 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 Radiation Oncology, BC Cancer, Vancouver, BC, CA
4Department of Surgery, Faculty of Medicine, University of British Columbia,
Vancouver, BC, CA

###### Abstract

Variability of IVIM parameters throughout the literature is a long-standing
issue, and perfusion-related parameters are difficult to interpret. We
demonstrate for improving the analysis of intravoxel incoherent motion imaging
(IVIM) magnetic resonance (MR) images, using image denoising and a
quantitative approach that does not require imposing specific exponential
models. IVIM images were acquired for 13 head-and-neck patients prior to
radiotherapy. Of these, 5 patients also had post-radiotherapy scans acquired.
Image quality was improved prior to parameter fitting via denoising. For this,
we employed neural blind deconvolution, a method of undertaking the ill-posed
mathematical problem of blind deconvolution using neural networks. The signal
decay curve was then quantified in terms of area under the curve ($AUC$)
parameters. Denoised images were assessed in terms of blind image quality
metrics, and correlations between their derived parameters in parotid glands
with radiotherapy dose levels. We assessed the method’s ability to recover
artificial pseudokernels which had been applied to denoised images. $AUC$
parameters were compared with the apparent diffusion coefficient ($ADC$),
biexponential, and triexponential model parameters, in terms of their
correlations with dose, and their relative contributions to the total variance
of the dataset, obtained through singular value decomposition. Image denoising
resulted in improved blind image quality metrics, and higher correlations
between IVIM parameters and dose. $AUC$ parameters were more correlated with
dose than traditional IVIM parameters, and captured the highest proportion of
the dataset’s variance. V This method of describing the signal decay curve
with model-independent parameters like the $AUC$, and preprocessing images
with denoising techniques, shows potential for improving reproducibility and
functional utility of IVIM imaging.

##  1 Introduction

Intravoxel Incoherent Motion (IVIM) Magnetic Resonance Imaging (MRI) is a
diffusion imaging technique which quantifies the translational motion of
molecules within imaging voxels [1]. This translational motion is due to the
random diffusion of water molecules within tissue, as well as perfusion
effects. These translational motions are inferred by analyzing signal decay in
each voxel as a function of applied diffusion gradient strengths (b-values).
Le Bihan first introduced a method of quantifying these in vivo motions using
Stejskal and Tanner’s [2] diffusion gradient method in 1986 [3]. This paper
introduced the concept of an ”apparent diffusion coefficient” (ADC). The
pseudodiffusion coefficient quantifies the combined effects of real diffusion
due to random molecular motion and the randomly oriented microcirculation of
blood in capillary beds. By acquiring signals using multiple diffusion
b-values, the ADC can be determined using the equation,

$S(b)/S(0)=e^{-ADC\cdot b}$ (1)

which is analagous to Stejkal and Tanner’s in vitro diffusion equation for a
still liquid, $S(b)/S(0)=e^{-D\cdot b}$ [2], where D is the self-diffusion
coefficient of the liquid. The ADC is a simple metric for describing the
signal decay curve in vivo; however, it does not model real diffusion, and is
only a method of approximating the biological mechanisms responsible for the
observed signal decay seen with increasing diffusion b-values. An important
difference is that Stejskal and Tanner’s diffusion equation comes from
directly solving the Bloch equations [4] for the magnetization vector,
$\vec{M_{0}}$, as a function of time, while Le Bihan’s application of this
equation to in vivo motion is an approximation based on this derivation. The
most useful way to interpret the $ADC$ is to think of it as simply a parameter
for approximating the in vivo signal versus b-value curve as an exponential
decay curve.

It is theorised that contributions of diffusion and perfusion motion to the
signal decay can be disentangled by fitting appropriate models to the signal
decay curve in each voxel, as originally posed by Le Bihan in 1988 [5]. This
operates by modelling the collective perfusion of blood through
microcapillaries as a ”pseudodiffusion” process, with a pseudodiffusion
coefficient, $D^{*}$. $D^{*}$ is typically on the order of 10 times larger
than $D$ [1], and therefore, the signal decay is largely attributed to
perfusion at low b-values, and diffusion at high b-values [6]. A common
approach for separating diffusion and pseudodiffusion in practice is to apply
a biexponential signal decay model to the data,

$S(b)/S(0)=f\,e^{-b\cdot D^{*}}+(1-f)\,e^{-b\cdot D}$ (2)

where $f$ denotes the fraction of the signal decay attributed to
pseudodiffusion within a voxel.

Reproducibility of the pseudodiffusion coefficient remains a long-standing
issue while D and ADC have been shown to exhibit higher stability [7, 8, 9,
10, 11, 12]. $f$ and $D^{*}$ have been shown to have inferior reproducibility
to $D$; however, $f$ has been shown to be empirically useful for predicting
age in nervous tissue [13, 14].

In recent years, the literature has seen a shift towards adoption of a
triexponential model [15, 16, 17, 18, 19, 20] of the form

$S(b)/S(0)=f_{1}\,e^{-b\cdot D_{2}^{*}}+f_{2}\,e^{-b\cdot
D_{2}^{*}}+(1-f_{1}-f_{2})\,e^{-b\cdot D}$ (3)

Model-fit uncertainty can be reduced by introducing an extra exponential term
[16, 17, 18, 19, 20]; however, it remains to be seen whether these parameters
are reproducible or have practical physiological interpretations.

It is uncertain whether attempts to separate diffusion from perfusion are
effective, as it has been reported that fitting a mono-exponential function,
$S(b)/S(0)=(1-f)\,e^{-b\cdot D}$ is more reliable than a bi-exponential fit
for differentiating pathological grades of esophageal squamous cell carcinoma
(ESCC) [12]. Heightened measurement error of signal at low b-values [21, 22]
leads to sub-optimal conditions for estimating perfusion effects, and
estimation of $ADC$ in low b-value regions has shown poor reproducibility
[23]. Noise-levels have also been shown to greatly impact parameter estimates
[24].

Poor reproducibility of perfusion-related parameters could be partially
explained by over-simplification of physiological processes inherent in simple
exponential models, as it has been shown that the optimal choice of model is
dependent on tissue type [25]. It is without doubt that, in reality, each
imaging voxel containing biological tissue will encompass a complex
arrangement of micro-structures, all having various perfusion fractions and
molecular motion. Kuai et al [15] used simulated data with 2-5 perfusion
components to show that the variance of $f$ and $D^{*}$ tend to increase as
the number of perfusion components increases, or the difference between
pseudodiffusion components increases.

Reproducibility of IVIM parameters is also impacted by a lack of
standardization for voxel sizes and b-value distributions, which several
studies have attempted to optimize [26, 27, 28, 29, 30, 31]. Variation in
voxel size is certain to impact multi-exponential model fits. For example,
suppose there exists IVIM images with a voxel size of $2$ x $2$ x $2$ $mm^{3}$
and it is fitted with a bi-exponential model. That image is then immediately
re-acquired with voxels of size $1$ x $1$ x $1$ $mm^{3}$ and again modeled
with a bi-exponential fit. By comparing the two images, it is clear that the
original image had been modelling octo-exponential components with a bi-
exponential function. Averaging the parameters from the smaller voxels over
the larger ones will, in general, yield different parameters than those of
larger voxels exactly paired. Voxel-by-voxel analysis of IVIM images is
difficult due to uncertainty, and parameters are therefore averaged over
regions-of-interest (ROIs). This issue is exacerbated by partial volume
effects (PVEs) from various tissue types and bleed-in from neighbouring
voxels.

The goal of applying various signal-decay models to IVIM data is to describe
the signal decay curve, which in turn, quantifies microscopic fluid motion
within tissue. Bi- and tri-exponential parameters are particularly enticing as
they appear to have a clear physiological interpretation. However, there is no
consensus that $f$ is correlated with blood vessel density [24, 32, 33], and
the ability of parameters to be interpreted physiologically is generally
unimportant for practical applications. There is furthermore no consensus on
the optimal IVIM model for describing the signal decay curve, with many
alternative versions being recommended [34, 35, 36, 37, 38].

IVIM parameters tend to be used to construct models for predicting various
clinical end-points, without being directly interpreted. While model
complexity has risen since the conception of IVIM MRI, the overall goal
remains simple: to best describe the signal decay curve. We hypothesize that
describing the signal decay curve without imposing any specific models onto
the data will lead to higher reproducibility, and better modelling of
outcomes.

The purpose of this study is to introduce two methods for improving IVIM MRI
efficacy and reproducibility. First, to improve image quality and reduce the
effect of image aberrations on parameter estimates, we apply a denoising
algorithm to each patient’s diffusion images prior to model-fitting. Second,
we introduce several model-independent parameters for quantifying the signal
decay curve of IVIM MRI images, and compare them to bi-exponential, tri-
exponential, and ADC parameters.

##  2 Methods

###  2.1 2.1 Dataset

This study was approved by an institutional review board, and written,
informed consent was obtained from all patients. The inclusion criteria for
recruitment were: (1) a diagnosis of nasopharynx, base-of-tongue, or tonsil
cancer to be treated with external beam radiotherapy; (2) expected to receive
dose in parotid glands during treatment; (3) scheduled to have MR images
acquired prior to radiotherapy and at 3 months following the last treatment
day. IVIM MR images encompassing the parotid glands were acquired during
routine clinical appointments for 12 head-and-neck cancer patients prior to
radiotherapy (Age: 35-78Y, average: 60.4 ; sex: 11M, 1F). Of these, 5 had
follow-up images acquired at 3 months post-RT (Age: 49-78Y, average: 59.3 ;
sex: 5M, 1F). 4 of these 5 patients received weekly chemotherapy concurrent
with radiotherapy (three patients: weekly 40 $mg/m^{2}$ cisplatin, one
patient: weekly 20-40 $mg/m^{2}$ per week and 240 mg carboplatin).

Images were acquired on a 1.5 T Magnetom Sola scanner (Siemens Healthineers)
with an echo planar imaging sequence (EPI). Diffusion weighted images were
acquired with a voxel size of 1.6 x 1.6 x (3.4-3.9) mm3. This slight inter-
patient variation in slice thickness was imposed to encompass parotid glands
within 32 slices while maximizing the signal-to-noise ratio (SNR). The
repetition time (TR) and echo time (TE) were 4900 ms and 105 ms, with a flip
angle of 90∘. Images were acquired for 16 b-values (0, 20, 30, 40, 50, 60, 70,
80, 90, 100, 120, 150, 250, 400, 800, 1000). These values were chosen to
follow recommendations for effectively separating perfusion parameters [39,
40], typically with a partition margin around 50 $s/mm^{2}$ [28]. Signals were
acquired in 6 directions, spanning the 3 principal Cartesian axes and their
vectors. Two signal averages were computed in each direction. Signals from
various directions were combined by taking the geometric mean [41].

Parotid glands were manually contoured by a graduate physics student on
clinical T1 images. These images were acquired with a turbo spin-echo (TSE)
sequence (TR: 564 ms, TE: 20 ms, flip angle: 148∘). T1 images were acquired
with a voxel size of 0.6 x 0.6 x 3 mm.

DICOM radiotherapy dose and planning structure set files used for external
beam radiotherapy treatment planning were exported from Varian Eclipse (Varian
Medical Systems, Inc.) for computing dose statistics inside parotid glands.
Whole-mean parotid gland dose levels were calculated within computed
tomography (CT) contours, manually defined by a single, senior clinical
radiation oncologist.

###  2.2 2.2 Denoising

To mitigate the influence of partial volume and other blurring effects on IVIM
parameter estimates, images were denoised using neural blind deconvolution.
Blind deconvolution is the ill-posed mathematical problem of estimating
hypothetical denoised images, $\vect{x}$, and their associated blur kernel,
$\vect{k}$, which convolve to yield the original image,
$\vect{y}=\vect{x}\ast\vect{k}$. This method is particularly useful in
situations where spatial resolution and precision of signal localization is
low, hence many studies have attempted to employ it for positron emission
tomography (PET) [42, 43, 44, 45, 46]. Neural blind deconvolution is the
method of solving this problem using neural networks which are simultaneously
optimized to predict $\vect{x}$ and $\vect{k}$. It was originally implemented
by Ren et al. [47] in 2020 and applied to images from a 2-dimensional natural
image database, and then significantly improved by Kotera et al. [48] in 2021.
Neural blind deconvolution is an unsupervised learning methodology which
trains and predicts on a case-by-case basis, without the requirement of a
separate training set with ground truth images.

Neural blind deconvolution was adapted for 3-dimensional prostate specific
membrane antigen (PSMA) positron emission tomography (PET) in 2023, while
incorporating simultaneous super-sampling into the methodology [49]. It was
shown to improve blind image quality metrics, and strengthen correlations
between PSMA PET uptake and sub-regional importance estimates in the parotid
gland for predicting post-radiotherapy xerostomia (subjective dry mouth) [50].

Here we further build off of this methodology for suitability with IVIM MRI.
First of all, we adapt the network to handle 4-dimensional input images (b, x,
y, z), composed of 3-dimensional images for all b-values stacked. Images were
normalized as $\vect{x}\rightarrow\frac{\vect{x}-\mu^{body}}{\sigma_{body}}$,
where $\mu_{body}$ and $\sigma_{body}$ are the mean and standard deviation of
signal in body voxels. The network architecture was maintained to predict a
single kernel, such that a single blur kernel was predicted using all b-value
images for each patient.

As the traditional formulation of the blind deconvolution problem,
$\vect{y}=\vect{x}\ast\vect{k}$, accounts for bleeding of voxel signals, it
does not account for background noise. As MRI signals contain rician noise
[51], we modified the mathematical formulation of the problem to include a
noise term:

$\vect{y}=\vect{x}\ast\vect{k}+\vect{\delta}$ (4)

This required the introduction of an additional symmetric convolutional auto-
encoder network, $G_{\delta}$, for predicting the noise term $\vect{\delta}$.
The networks used are thus,

$G_{x}(\theta_{x})=\vect{x}\,\,\,\,,\,\,\,\,G_{\delta}(\theta_{\delta})=\vect{\delta}\,\,\,\,,\,\,\,\,G_{k}(\theta_{k})=\vect{k}\,\,\,\,,\,\,\,\,\vect{y}=\vect{x}\ast\vect{k}+\vect{\delta}.$
(5)

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

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

The optimization problem can be written in terms of the model parameters
$\theta_{x}$, $\theta_{\delta}$ and $\theta_{k}$ as

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

The full architecture is summarized visually in Figure 1.

[width=0.8]network_architecture.png

Figure 1: The blind deconvolution architecture used for denoising IVIM images
is illustrated. $G_{x}$ and $G_{\delta}$ are symmetric, convolutional auto-
encoder networks for predicting the denoised image, $\vect{x}$, and its
additive noise contributions, $\vect{\delta}$. $G_{k}$ is a fully-connected
network for predicting the blur kernel, $\vect{k}$. The fidelity loss term
used in optimization compares the original image with
$\vect{x}\ast\vect{k}+\vect{\delta}$ by using back-propagated model gradients
from the loss function for iterative optimization.

Optimization was performed using a single NVIDIA GeForce GTX 1060 GPU and a
2.8GHz Intel CoreTM i5-8400 CPU. As neural blind deconvolution is a
computationally expensive procedure, further exacerbated by stacking multiple
b-value images, images were cropped to extend at least 5 voxels outside of
parotid gland borders in all 3 cartesian directions. This was performed
separately for the left and right parotid glands, such that neural blind
deconvolution was performed separately for each gland. Decreasing the image
size in this manner was necessary to limit GPU memory usage.

A multi-scale optimization procedure was implemented, as recommended for
neural blind deconvolution by Kotera et al. [48] and recently implemented for
PET denoising [49, 50]. This involved pre-training networks using images that
were first down-sampled by a factor of 2, then up-sampling final outputs by
$\sqrt{2}$ before performing a second round of pre-training. After these two
pre-training rounds, the outputs, $G_{x}$, $G_{k}$, and $G_{\delta}$, are used
as initial inputs for the regular training procedure. An 11 x 11 x 11 kernel
was used, which was down-sampled to 5 x 5 x 5, and then 7 x 7 x 7 during pre-
training stages. The algorithm for updating network weights to predict
denoised diffusion images is summarized in Table 1

justification=raggedright  Diffusion Denoising Algorithm Input: Original 4D
(b, z, y, x) image, $\vect{y}$ Output: Deblurred 4D (b, z, y, x)image,
$\vect{x}$, 4D (b, z, y, x) noise image, $\vect{\delta}$, and 3D blur kernel,
$\vect{k}$ Pre-training 1\. Downsample spatial dimensions of $y$ to 1/2
resolution 2\. Initialize $\vect{z_{x}}$, $\vect{z_{\delta}}$ from uniform
distribution to match size of $\vect{y}$ 3\. Initialize $\vect{z_{k}}$ as size
$5\times 5\times 5$ Gaussian kernel with standard deviation of one voxel 4\.
for` i = 1 to 1000:` 5\.  $\vect{x}=G^{i}_{x}(\vect{z_{x}})$ 6\.
$\vect{k}=G^{i}_{k}(\vect{z_{k}})$ 7\.
$\vect{\delta}=G^{i}_{\delta}(\vect{z_{\delta}})$ 8\.  Compute loss and back-
propagate gradients 9\.  Update $G^{i}_{x}$, $G^{i}_{\delta}$ and $G^{i}_{k}$
using the ADAM optimizer [52] 10\. $\vect{x}=G^{1000}_{x}(\vect{z_{x}})$,
$\vect{k}=G^{1000}_{k}(\vect{z_{k}})$,
$\vect{\delta}=G^{1000}_{\delta}(\vect{z_{\delta}})$ 11\. Upsample $\vect{k}$
by 7/5 and upsample spatial dimensions of $\vect{x}$ and $\vect{\delta}$ by a
factor of $\sqrt{2}$ 12\. Downsample spatial dimensions of $\vect{y}$ by
$\sqrt{2}$ to match new convolution size 13\. $\vect{z_{x}}=\vect{x}$,
$\vect{z_{\delta}}=\vect{\delta}$, $\vect{z_{k}}=\vect{k}$ 14\. Repeat steps 4
through 10. 15\. Upsample $\vect{k}$ by 11/7 and upsample spatial dimensions
of $\vect{x}$ and $\vect{\delta}$ by a factor of $\sqrt{2}$ 16\.
$\vect{z_{x}}=\vect{x}$, $\vect{z_{\delta}}=\vect{\delta}$,
$\vect{z_{k}}=\vect{k}$ Main Training 17\. for` i = 1 to 5000:` 18\.
$\vect{x}=G^{i}_{x}(\vect{z_{x}})$ 19\.  $\vect{k}=G^{i}_{k}(\vect{z_{k}})$
20\.  $\vect{\delta}=G^{i}_{\delta}(\vect{z_{\delta}})$ 21\.  Compute loss and
back-propagate gradients 22\.  Update $G^{i}_{x}$, $G^{i}_{\delta}$ and
$G^{i}_{k}$ using the ADAM optimizer [52] 23\.
$\vect{x}=G^{5000}_{x}(\vect{z_{x}})$, $\vect{k}=G^{5000}_{k}(\vect{z_{k}})$,
$\vect{\delta}=G^{1000}_{\delta}(\vect{z_{\delta}})$

Table 1: The optimization algorithm for updating network weights to predict
deblurred PSMA PET images. This algorithm builds off of Ren et al.’s proposed
joint optimization algorithm [47] and implements modifications suggested by
Kotera et al. [48].

The regular training procedure consisted of 5000 iterations, as in previous
implementations. This was split into two stages, with two different loss
functions. These stages and their corresponding loss terms are summarized in
Table 2. It is necessary to include regularization terms in the loss function
to avoid convergence to a trivial solution ($\vect{k}$ becomes a unit
impulse). Therefore, the mean squared error (MSE) of kernel values above 0.7
were penalized. We only penalized values above this cutoff to avoid a trivial
solution while still allowing small contributions from neighbouring voxels to
go unpenalized. We furthermore penalized the MSE of noise voxels above
$\mu_{\delta}+\sigma_{\delta}$, where $\mu_{\delta}$ and $\sigma_{\delta}$ are
the mean and standard deviation of the normalized original image voxels
located outside of the body. A total variation (TV) loss term for $\vect{x}$
was included after the 3500th iteration, as this has been shown to improve
results when employed in later stages of optimization [48]. Furthermore, our
fidelity loss function switched from MSE to the structural similarity index
metric (SSIM) after the 3500th iteration, as recommended [48] for improving
final image quality.

justification=raggedright  Iteration$<3500$ Iteration $\geq 3500$ Fidelity
Term MSE SSIM Regularization Terms 1\. Kernel 1\. Kernel 2\. Noise 2\. Noise
3\. TV

Table 2: Optimization consisted of 5000 iterations, split into 3 stages, each
having a different loss function. Up to the 3500th iteration, the fidelity
loss function was mean squared error (MSE). The structural similarity index
metric (SSIM) was used in remaining iterations. In all iterations, a kernel
and noise regularization term were included. The MSE of kernel values were
penalized to avoid convergence to the trivial solution. The MSE of noise
values above $\mu_{\delta}+\sigma_{\delta}$, where $\mu_{\delta}$ and
$\sigma_{\delta}$ are the mean and standard deviation of the normalized
original image voxels located outside of the body. Finally, the total
variation (TV) between denoised image voxels was penalized after the 3500th
iteration, as it has been shown to improve results when employed in later
stages of optimization [48]

###  2.3 2.3 Evaluating Denoised Images

As neural blind deconvolution is an unsupervised training methodology, ground
truth images were not used for direct evaluation of denoised images. We
compared image quality between denoised and original IVIM images in terms of
two blind image quality metrics. The Blind/Referenceless Image Spatial Quality
Evaluator (BRISQUE) [53] and Contrastive Language Image Pre-training (CLIP)
[54] score were evaluated. BRISQUE ranks image quality on a scale between 0
and 100 based on a series of features derived from various signal intensity
and distribution statistics. These features are used to predict the deviation
of the input image from a natural, undistorted image. CLIP uses a
language/vision neural network trained using millions of natural image/caption
pairs to assess either image quality or caption quality. We furthermore
assessed the quality of IVIM parameter maps predicted using original and
denoised images in terms of BRISQUE and CLIP metrics.

We furthermore assessed the similarity between predicted kernels amongst
patients. As all patient images were acquired on the same scanner, predicted
blur kernels were expected to be similar. As the blur kernel contains
components due to both scanner uncertainty, as well as patient motion, we
anticipated predicted blur kernels to be similar, yet nonidentical, between
patients. We compared blur kernels quantitatively by computing the inner
product between normalized kernels.

To ensure the model’s ability for predicting accurate blur kernels, we
employed the same strategy used for evaluating neural blind deconvolution for
PSMA PET images [49]. This involves generating artificial kernels, convolving
them with previously denoised images, then re-starting the neural blind
deconvolution process. The objective is then to test how accurately
pseudokernels can be predicted, measured using the inner product between
generated and predicted pseudokernels. Four types of pseudokernels were
generated: a regular Gaussian with a standard deviation of 1 voxel in each
direction, and three oblong Gaussians, each having a standard deviation
truncated by a factor of three in one of the three Cartesian directions.
Pseudokernels were all normalized to a unit sum.

Original and denoised images were further compared in terms of their practical
ability for model-building. Beginning with the hypothesis that IVIM parameters
are related to histological features within imaging voxels, we tested the
ability of radiotherapy dose levels to predict changes in IVIM parameters
following treatment. This was assessed in two different ways. First of all,
Spearman’s rank correlation coefficient ($r_{s}$) was compared between the
whole-mean dose and relative changes in mean IVIM parameters within each
parotid gland ($\frac{P_{f}-P_{i}}{{P_{i}}})$, where $P_{i}$ and $P_{f}$ are
the parameters before and after radiotherapy).

The dose response of parameters were further analysed using a multiple
regression analysis. Post-radiotherapy parameters were used as dependent
variables while pre-radiotherapy parameters and mean dose levels served as two
independent variables. The best-fit slope with respect to dose was then
compared between parameters. For comparison’s sake, slopes were calculated
with respect to parameters that had been normalized to statistical Z-values,
with statistics derived using pre-radiotherapy voxels within each patient’s
parotid glands. Images were normalized to statistical Z-value rather than by
the maximum image voxel, as it was found that maximum values tended to
substantially exceed average values within parotid glands, leading to
compression of voxel values.

###  2.4 2.4 Extracting Exponential Model Parameters

To serve as a basis of comparison for newly derived, model-independent IVIM
parameters, we extracted ADC, bi-exponential, and tri-exponential parameter
maps for all patients. Voxel-wise parameter fitting was performed using
IVIM3-NET models [55], trained for each model type using all acquired original
and denoised IVIM images. Voxels were first filtered to exclude voxels located
outside the body from training. This resulted in approximately 3.6 million
training voxels. Voxels for each b-value were normalized to signals acquired
with $b=0$. IVIM3-NET hyper-parameters were tuned according to previously
derived optimal values [55]. Models consist of fully-connected networks with 2
hidden layers, which take normalized signal versus b-value arrays as input,
and predict the appropriate number of parameters for each model. Each layer
consists of a linear function, an exponential linear unit (ELU) activation
function [56], and batch normalization. Sigmoid functions were used to
constrain parameter predictions within boundaries chosen to over-encompass
anticipated ranges (Table3). Network weights were then optimized via back-
propagation using a batch size of 256 and an ADAM optimizer [52].

To mediate the effect of image noise at $S(b=0)$ affecting all normalized
measurements acquired at other b-values, IVIM3-NET models were also used to
simultaneously predict $S(b=0)$ along with other model parameters. This was
shown to improve fit results in a previous study [55]. Hence, the loss term is
of the form,

$L=\left(\frac{S(b)}{\hat{S_{0}}}-\hat{S}(b)\right)^{2}$ (8)

where $S(b)$ is the actual signal in a given voxel, and $\hat{S_{0}}$ and
$\hat{S}(b)$ are model predictions. For example, a bi-exponential model
outputs $\hat{S_{0}}$, $f$, $D^{*}$ and $D$, which yield
$\hat{S}(b)=fe^{-b\cdot D^{*}}+(1-f)e^{-b\cdot D}$

Dark_Grayrgb0.5,0.5,0.5 Grayrgb0.7,0.7,0.7 whiterrgb0.8,0.8,0.8 whitergb1,1,1
a¿Grayc b¿whiterc d¿whitec

justification=raggedright  Dark_GrayModel Parameter Lower Bound Upper Bound
GrayBi-exponential $D$ (mm2/s) $7\times 10^{-5}$ $4\times 10^{-3}$ Gray
$D^{*}$ (mm2/s) $2\times 10^{-3}$ $0.15$ Gray $f$ $0$ $0.5$ whiteTri-
exponential $D$ (mm2/s) $7\times 10^{-5}$ $4\times 10^{-3}$ white
$D^{*}_{1}$(mm2/s) $3\times 10^{-3}$ $2\times 10^{-2}$ white
$D^{*}_{2}$(mm2/s) $2\times 10^{-2}$ 0.7 white $f_{1}$ 0 0.5 white $f_{2}$ 0
0.5 GrayADC $ADC$(mm2/s) $7\times 10^{-5}$ $5\times 10^{-3}$

Table 3: Parameter bounds set for IVIM3-NET predictions of exponential model
parameters. Bounds were set to over-encompass anticipated parameter ranges, in
order to compensate for diminishing gradients at sigmoid asymptotes [57]

###  2.5 2.5 Model-independent IVIM parameters

To describe the signal decay curve without imposing specific mathematical
models onto the data, we define the area under the curve ($AUC$) parameter.
The $AUC$ captures the general fall-off of signal with increasing b-value, and
can be equivalently described as the integral of $S(b)/S(0)$. The anticipated
relationship between $AUC$ and exponential parameters can be found by
integrating over exponential functions between $b=0$ to $b=1000$. For the ADC
model,

$AUC=\int^{1000}_{0}\frac{S_{ADC}(b)}{S_{ADC}(0)}db=\frac{1}{ADC}\left[1-e^{-ADC\cdot
1000}\right].$ (9)

$AUC$ and ADC are thus inversely proportional. For the biexponential model,

$AUC=\int^{1000}_{0}\frac{S_{bi}(b)}{S_{bi}(0)}db=\frac{f}{D^{*}}\left[1-e^{-1000\cdot
D^{*}}\right]+\frac{1-f}{D}\left[1-e^{-1000\cdot D}\right]$ (10)

. $AUC$ decreases with increasing D, and is roughly independent of $D^{*}$
assuming $D^{*}\approx 10D$, and f small). For the triexponential model,

$AUC=\int^{1000}_{0}\frac{S_{tri}(b)}{S_{tri}(0)}db=\frac{f_{1}}{D_{1}^{*}}\left[1-e^{-1000\cdot
D_{1}^{*}}\right]+\frac{f_{2}}{D_{2}^{*}}\left[1-e^{-1000\cdot
D_{2}^{*}}\right]+\frac{1-f_{1}-f_{2}}{D}\left[1-e^{-1000\cdot D}\right].$
(11)

$AUC$ is roughly independent of $D_{1}^{*}$ and $D_{2}^{*}$ (assuming
$D_{i}^{*}\approx 10D$ and $f_{i}$ small, $i\in\\{1,2\\}$).

The $AUC$ was calculated for each voxel using a middle Riemann sum. The AUC
was furthermore calculated in 3 subsets of b-values, defining a low-, middle-,
and high-range $AUC$ ($AUC_{l}$, $AUC_{m}$, $AUC_{h}$, respectively). The
ranges were defined as:

  * •

Low: $b\leq 120$

  * •

Medium: $120\leq b\leq 400$

  * •

High: $400\leq b\leq 1000$

The method of calculating $AUC$, $AUC_{l}$, $AUC_{m}$, and $AUC_{h}$ using
middle Riemann sums is shown graphically in Figure 2.

[width=0.9]auc_graph.png

Figure 2: $AUC$ parameters were calculated using middle Riemann sums. On the
left, the $AUC$ within a single voxel is calculated graphically. In a separate
voxel on the right, $AUC_{L}$, $AUC_{M}$, and $AUC_{H}$ are calculated
(yellow, rose, blue), which sum to the regular $AUC$.

###  2.6 2.6 Comparing Parameters

Whole-mean parameters were computed within all parotid glands. To compare
functional utility of $AUC$ parameters with exponential model parameters, we
employed several strategies. First of all, we included $AUC$ parameters in the
previously mentioned dose-response comparison of IVIM parameters in parotid
glands. This allowed the practical utility of parameters for use in dose
response studies to be analysed, based on the hypothesis that parameter
statistics are related to tissue composition and characteristics at a
histological level.

We further tested the relative importance of IVIM parameters in capturing the
variance of the dataset by evaluating contributions to the primary principal
component obtained through singular value decomposition. The primary principal
component is the optimal linear combination of all input features for
describing the maximum variation in the data. This was computed by first
creating a 2-dimensional array, where each row represents a single voxel in an
image and each column represents a given feature. Voxels within the body in
all cropped images were included as rows in this array. All biexponential,
triexponential, $ADC$, and $AUC$ features were included as columns. This
yielded an $m\times 13$ matrix for which the singular value decomposition was
computed. The squared projection of the primary principal component on each of
the 13 parameters was then calculated and plotted.

##  3 Results

Neural blind deconvolution for denoising diffusion MR images appeared to
effectively suppress random noise while enhancing fine detail within images.
This improvement is supported quantitatively by an increase in blind image
quality metrics, BRISQUE and CLIP for denoised images (Table 4).

justification=raggedright  Denoised Original Significance BRISQUE $36.4\pm
7.8$ $14.6\pm 6.9$ $p<0.001$ CLIP $0.37\pm 0.03$ $0.35\pm 0.03$ $p<0.001$

Table 4: Blind image quality metrics, CLIP and BRISQUE are compared for
original ($\vect{y}$) and denoised ($\vect{x}$) diffusion MR images. Results
were averaged over all slices and b-values. The significance of the difference
was tested using a paired t-test.

Predicted blur kernels were similar among all patients, having a mean inter-
patient inner product of $0.88\pm 0.10$. Furthermore, each patient had two
predicted kernels (for diffusion images cropped over the left and right
parotid, separately) whose mean inner product was $0.97\pm 0.02$. Projections
of the mean predicted kernel are shown in Fig 4. Kernel voxels were mainly
confined to within axial planes, with the highest variance in the anterior-
posterior direction. The general shape of pseudokernels was predicted using
neural blind deconvolution; however, it can be seen that predicted kernel
voxels were mostly confined to within one voxel of the centre value (Fig 5).

Axial slices of denoised and original diffusion images corresponding to
several different b-values are contrasted in Fig 3. Image denoising resulted
in significantly ($p<0.001$) more monotonically decreasing signal versus
b-value curves, as shown in Fig 6. This was assessed by first taking
Spearman’s rank correlation coefficient, $r_{s}$, between signal and b-values,
averaged over all voxels within body tissue. A paired t-test was then used to
assess differences between correlations in denoised and original images.
$r_{s}$ was ${-}0.80\pm 0.31$ in denoised images and $-0.60\pm 0.29$ in
original images.

$ADC$, biexponential, triexponential, and all $AUC$ parameter statistics in
parotid glands are summarized in Table 5.

Grayrgb0.5,0.5,0.5 whiterrgb0.8,0.8,0.8 whitergb1,1,1 a¿Grayc b¿whiterc
d¿whitec

justification=raggedright  Denoised Images Original Images whiteParameter Pre-
RT Post-RT p-value Pre-RT Post-RT p-value p-value (Deblur versus Orig) ADC
($\times 10^{-3}mm^{2}/s$) $1.23\pm 0.37$ $1.98\pm 0.20$ $<0.001$ $1.27\pm
0.36$ $1.96\pm 20$ $<0.001$ $>0.05$ D (biexponential, $\times
10^{-4}mm^{2}/s$) $4.8\pm 1.1$ $6.6\pm 1.2$ $<0.01$ $4.7\pm 1.2$ $6.5\pm 1.2$
$<0.01$ $<0.001$ $D^{*}$ (biexponential, $\times 10^{-2}mm^{2}/s$) $1.6\pm
0.8$ $1.1\pm 0.2$ $>0.05$ $2.2\pm 0.9$ $1.4\pm 0.4$ $<0.01$ $<0.001$ $f$
(biexponential) $0.25\pm 0.04$ $0.32\pm 0.03$ $<0.001$ $0.26\pm 0.03$ $0.32\pm
0.03$ $<0.001$ $<0.001$ $D$ (Triexponential, $\times 10^{-4}mm^{2}/s$) $4.6\pm
1.1$ $6.3\pm 1.1$ $<0.001$ $4.5\pm 1.2$ $6.3\pm 1.1$ $<0.001$ $<0.001$
$D_{1}^{*}$ (Triexponential, $\times 10^{-3}mm^{2}/s$) $6.2\pm 0.9$ $5.2\pm
0.6$ $>0.05$ $6.9\pm 1.0$ $6.1\pm 0.7$ $<0.01$ $<0.001$ $f_{1}$
(triexponential) $0.20\pm 0.04$ $0.27\pm 0.03$ $<0.001$ $0.20\pm 0.03$
$0.26\pm 0.03$ $<0.001$ $p<0.001$ $D_{2}^{*}$ (Triexponential, $\times
10^{-1}mm^{2}/s$) $0.70\pm 0.30$ $0.55\pm 0.3$ $>0.05$ $1.3\pm 0.3$ $0.9\pm
0.3$ $>0.05$ $p<0.001$ $f_{2}$ (Triexponential) $0.04\pm 0.01$ $0.03\pm 0.01$
$>0.05$ $0.06\pm 0.01$ $0.042\pm 0.001$ $<0.01$ $p<0.001$ $AUC$ ($s/mm^{2}$)
$652\pm 41$ $587\pm 32$ $<0.001$ $642\pm 39$ $580\pm 32$ $<0.001$ $p<0.001$
$AUC_{L}$ ($s/mm^{2}$) $108\pm 2.4$ $108\pm 2.4$ $>0.05$ $105\pm 1.6$ $107\pm
1.7$ $>0.05$ $p<0.001$ $AUC_{M}$ ($s/mm^{2}$) $205\pm 8.7$ $191\pm 9.0$
$<0.01$ $201\pm 7.7$ $189\pm 7.6$ $<0.01$ $p<0.001$ $AUC_{H}$ ($s/mm^{2}$)
$443\pm 39.4$ $382\pm 28$ $<0.01$ $437\pm 37$ $378\pm 29$ $<0.01$ $p<0.001$

Table 5: All IVIM parameter averages and standard deviations inside parotid
glands are listed. Averages were calculated over parotid glands from pre-
radiotherapy (pre-RT) and post-RT scans separately, along with the
significance of their difference. Results are calculated using both denoised
and original images. Lastly, the final column lists the significance of
differences between denoised and original values before radiotherapy.
Significance values were determined using a paired t-test.

Image denoising was found to strengthen correlations between whole-mean
parotid dose levels and changes in IVIM parameters (Fig 9). Changes in $AUC$,
$AUC_{M}$, and $AUC_{H}$ and $ADC$ had the strongest correlations with mean
dose levels. Correlations between dose and $AUC$ parameters were particular
strong in denoised images. Regression slopes for changes in parameters versus
dose are also shown in Fig 9. Changes in $ADC$ had the highest slope with
respect to dose in both denoised and original images. $AUC$, $AUC_{M}$,
$AUC_{H}$, $D^{*}_{2}$, $f$, and $D_{biexp}$ had similarly high dose slopes,
while $D^{*}_{biexp}$, $D_{triexp}$, $D^{*}_{1}$, $f_{1}$, $f_{2}$, and
$AUC_{L}$ had lower slopes. Fig 9 also shows the relative contribution of each
parameter to the first principal component obtained via the SVD analysis.
$AUC$, $AUC_{H}$ and $ADC$ captured the highest proportion of variance in the
data, with $f_{1}$ also capturing a relatively high proportion of the
variance. Pairwise correlations between all parameter combinations are
displayed in Fig 9. $AUC$, $AUC_{M}$, and $AUC_{H}$ were highly correlated
with one another, and also with $ADC$, $D$ (biexponential and triexponential),
$f$, and $f_{1}$, while $AUC_{L}$ was most correlated with $f_{2}$ and
$D^{*}_{1}$.

$AUC$ parameter maps derived using denoised and original diffusion images are
shown in Fig 8. Parameter maps of biexponential parameters corresponding to
the same slice are shown in Fig 7.

[height=1.2]deblurred_and_original_v_b.png

Figure 3: Axial slices through parotid glands of denoised (left) and original
(right) diffusion images at 4 different b-values are shown. Denoiding appears
to effectively suppress noise while enhancing fine detail within images. Here
in particular, the posterior region of the parotid gland becomes discernible
only upon denoising.

[width=1]kernel.png

Figure 4: The average blur kernel predicted for patients is projected onto the
three primary imaging planes. The kernel was mostly confined to axial planes,
which was expected due to elongation of voxels in the superior-inferior
direction (slice thickness ${>}2\times$ axial slice pixel spacing). The kernel
had a greater spread in the anterior-posterior direction than the left-right
direction.

[width=1]pseudokernels.png

Figure 5: The neural blind deconvolution denoising process was tested by
convolving fixed pseudokernels with previously denoised diffusion images
before restarting the denoising process. Four pseudokernels, including a
regular Gaussian, as well as a Gaussian elongated in each of the 3 primary
image axes, were applied.

[width=1.1]signal_w_decay.png

Figure 6: Axial slices through the parotid gland of denoised and original
diffusion images with $b=1000$ are shown, along with the signal decay curve of
each in a single voxel (indicated in the denoised image). Signal values were
normalized to statistical Z-values within the body. Denoising appears to
suppress random noise while illuminating fine detail within images. The loss
function used during image denoising did not impose a decrease in signal value
with increasing b-value, yet signal values in denoised images appear to
naturally decrease more monotonically and more smoothly than in original
images.

[width=1.1]biexp_image_slices.png

Figure 7: Axial slices through the parotid gland of biexponential model
parameter maps are shown, derived using denoised (top) and original (bottom)
diffusion MR images.

[width=1.1]auc_image_slices.png

Figure 8: Axial slices through the parotid gland of biexponential parameter
maps are shown, derived using denoised (top) and original (bottom) diffusion
MR images. Axial $b=0$, $AUC$, $AUC_{L}$, $AUC_{M}$, and $AUC_{H}$ image
slices derived from denoised (top) and original (bottom) are shown.

[width=0.9]svd.png

Figure 9: Pairwise Spearman’s rank correlation coefficients, $r_{s}$, between
dose and relative changes in biexponential, triexponential, ADC, and AUC
parameters derived from IVIM images are displayed (top) for denoised and
original images (left and right). Best-fit regression slopes of parameters
versus dose are shown underneath in green. The singular value decomposition
(SVD) of the data was computed, and the relative variance captured by each
parameter in the first principal component was plotted (bottom).

##  4 Discussion

There is high variability of derived IVIM MRI parameters in the literature.
The proposed model-independent parameter, the $AUC$, proved to be an effective
parameter for characterizing the signal decay versus b-value curve. In terms
of stand-alone model-building potential, $AUC$ parameters may be more
effective than traditional IVIM parameters for predicting various clinical
outcomes, based on their high correlations with radiotherapy dose levels and
high contributions to the primary principal component of the total dataset.
Changes in the $AUC$, $AUC_{M}$ and $AUC_{H}$ following radiotherapy were more
correlated with dose than diffusion coefficients derived from exponential
models. Furthermore, changes in the $ADC$ following radiotherapy were more
correlated with dose than diffusion coefficients derived from biexponential
and triexponential models. Perfusion fractions for biexponential and
triexponential models, $f$ and $f_{1}$, were more correlated with dose than
their corresponding model’s diffusion coefficient, $D$.

Neural blind deconvolution for denoising IVIM MRI prior to parameter fitting
proved to be effective. Denoising enhanced fine detail in images while
suppressing noise, resulting in higher blind image quality metrics, BRISQUE
and CLIP [53, 54]. denoising resulted in higher correlations between dose and
IVIM parameter values, especially with $AUC$, $AUC_{M}$, and $AUC_{H}$. This
suggests that denoising may be beneficial for future studies seeking to use
IVIM imaging to assess radiotherapy-induced tissue damage. $AUC$ values
significantly decreased following radiotherapy. $AUC_{L}$ was found to be the
least effective $AUC$ variant, possibly due to the low amount of dataset
variability it accounted for, and its lower correlation with dose compared
with other $AUC$ variants. This is similar to the limited functional utility
of pseudodiffusion coefficients, which also describe the behaviour of the
signal decay curve in the low b-value regime.

In both denoised and original images, the newly proposed Area Under the Curve
($AUC$) of the normalized signal decay curve, and the $AUC$ in the high
b-value region ($AUC_{H}$) captured the largest proportion of the variance.
The $ADC$ captured a higher proportion of the variance than the diffusion
coefficients of biexponential and triexponential models. Perfusion fractions,
$f$ (biexponential) and $f_{1}$ (triexponential) captured a higher proportion
of the variance than diffusion or pseudodiffusion coefficients.
Pseudodiffusion coefficients and $AUC_{L}$ constituted a particularly low
proportion of the data’s variance. These variables are less correlated with
other variables, as seen in the correlations grid.

The results of this study suggest that within our dataset, $AUC$ parameters
and the $ADC$ are the most effective IVIM parameters for modelling dose-
related outcomes. However, this study included a small number of patients
having post-radiotherapy data, rendering our dose response quantities prone to
uncertainty. Further validation is needed to assess the efficacy of $AUC$
parameters and neural blind deconvolution for improving the practical utility
of IVIM imaging. It appears that denoising images could improve parameter
reproducibility, but more work is required for this to be confirmed.

Predicted blur kernels were highly similar between patients, and kernels
predicted for images encompassing the left and right parotid gland of each
patient were nearly identical (mean inner product: $0.97\pm 0.02$). Predicted
kernels were mainly confined to axial planes, which is expected, due to the
axial elongation of the voxel geometry used in this study. To avoid
convergence towards a trivial solution, where the kernel becomes a delta
function, a regularization term penalizing the MSE of kernel values above 0.7
was included. Previous neural blind deconvolution studies have simply applied
MSE to all kernel voxels [48, 47], but we only penalized voxels with values
above a given cutoff in order to allow small contributions from neighbouring
voxels to go entirely unpenalized.

This method successfully approximated the general shape of all four
pseudokernel types, albeit with truncation at the edges. Predictions were
mostly confined to within one voxel of the center value. This indicates that
denoising is effective for mitigating cross-talk between neighbouring voxels,
but does not account for non-neighbouring voxel interactions. Similar
truncation is reported for PSMA PET images [49]; however, voxels appear
slightly more centralized in this study. This may be partially explained by
the smaller kernel size employed in this study. As we were anticipating a
relatively small deblur kernel for IVIM MRI images (due to generally sharper
MR images than PET), while also considering memory constraints, we decreased
the kernel size to 11x11x11 from the 15x15x15 size used for PSMA PET [49, 50].
Overall, this denoising process appears to be effective for correcting cross-
talk between neighbouring voxels, but more distant effects are not accounted
for.

Parameters obtained from denoised and original images had similar population
means and standard deviations, but paired t-tests between values found all
parameters to be significantly different, except for the $ADC$. $AUC$ values
tended to be slightly higher in denoised images than original images, although
differences did not exceed one standard deviation. In all cases, differences
in parameters following radiotherapy followed the same trend. Mean ADC values
agreed with those previously reported in the literature [58, 59, 60]. $D$ and
$D^{*}$ were lower than previously reported values [6, 61, 60], both before
and after radiotherapy. However, $D$ and $D^{*}$ were within a standard
deviation of the mean values reported by Kimura et al. [60]. Furthermore,
diffusion coefficients increased after radiotherapy with similar proportions
to those found previously [6]. $f$ within parotid glands was lower than values
reported by Zhou et al. [6] and Beckert et al. [61] but smaller than values
reported by Kimura et al. [60].

It is not uncommon for IVIM parameters reported in the literature to be in
conflict [27], which was the motivation for this work. For example, the
significant variability of reported IVIM parameters in the liver has created
difficulties establishing normal, baseline perfusion parameters [62]. The low
SNR associated with echo-planar imaging is likely a significant contributor to
these issues. Disregarding time constraints, a spin echo sequence may lead to
much greater reliability in parameter estimates. Variations in SNR between
scanners used in different studies has been reported to affect parameter
estimates [62, 25]. Variability of b-value distributions [63] and voxel sizes
also contribute to variability among parameter estimates in the literature.

The intrinsically low SNR of echo planar imaging sequences used for IVIM
imaging makes denoising an important endeavour. Incorporating a noise term
into the neural blind deconvolution methodology allowed for corrections of
non-negative, scanner related noise to be corrected for. Noise suppression
after denoising is clearly visible in voxels located outside of body regions
(Fig 3). A notable result is the increased monotonicity of the signal decay
curve found after image denoising (Fig 6). It is expected that increasing
diffusion gradient strength should result in strictly decreasing signal
acquisition; however, this expectation was not enforced in the loss function
used for optimization. This perceivable suppression of noise in the decay
curve seems to validate the denoising process.

Characterizing the signal decay curve non-parametrically can be advantageous,
as it has the potential to improve the utility and reproducibility of IVIM
parameters without presuming the underlying mathematical model. While it is
more desirable in theory to have physically interpretable parameters, it has
not proven to be practical for IVIM imaging. The AUC and its sub-quantities
are only one method of describing the signal decay curve. While this study
affirms their utility, there remains room for further exploration and
validation of other model-independent parameters. Future studies can continue
to improve upon this methodology using the code available online [64].

##  5 Conclusion

Characterizing the signal decay curve using model-independent metrics, such as
the the $AUC$ parameters, and denoising images prior to parameter estimation,
could improve reproducibility and the practical utility of IVIM imaging for
developing outcome-predictive models. IVIM parameters have shown substantial
variability throughout the literature. A shift in the methodology for
characterizing the signal decay curve offers a potential solution to mitigate
this inconsistency.

##  6 Acknowledgements

We thank the entire MRI staff, and radiation oncologists at BC Cancer
Vancouver for helping with the data collection. We thank Mrs. Rosalie Segal
and the late Mr. Joseph Segal for funding this research. We acknowledge an in-
kind software loan from Siemens Healthcare Limited.

## References

## References

  * [1] D. Le Bihan, “What can we see with ivim mri?,” NeuroImage, vol. 187, pp. 56–67, 2019.  Physiological and Quantitative MRI.
  * [2] E. O. Stejskal and J. E. Tanner, “Spin diffusion measurements: Spin echoes in the presence of a time-dependent field gradient,” The Journal of Chemical Physics, vol. 42, pp. 288–292, Jan. 1965.
  * [3] D. L. Bihan, E. Breton, D. Lallemand, P. Grenier, E. Cabanis, and M. Laval-Jeantet, “MR imaging of intravoxel incoherent motions: application to diffusion and perfusion in neurologic disorders.,” Radiology, vol. 161, pp. 401–407, Nov. 1986.
  * [4] F. Bloch, “Nuclear induction,” Phys. Rev., vol. 70, pp. 460–474, Oct 1946\.
  * [5] D. Le Bihan, E. Breton, D. Lallemand, M.-L. Aubin, J. Vignaud, and M. Laval-Jeantet, “Separation of diffusion and perfusion in intravoxel incoherent motion mr imaging.,” Radiology, vol. 168, no. 2, pp. 497–505, 1988.
  * [6] N. Zhou, C. Chu, X. Dou, M. Li, S. Liu, L. Zhu, B. Liu, T. Guo, W. Chen, J. He, J. Yan, Z. Zhou, X. Yang, and T. Liu, “Early evaluation of irradiated parotid glands with intravoxel incoherent motion MR imaging: correlation with dynamic contrast-enhanced MR imaging,” BMC Cancer, vol. 16, Nov. 2016.
  * [7] D. M. Koh, E. Scurr, D. J. Collins, A. Pirgon, B. Kanber, N. Karanjia, G. Brown, M. O. Leach, and J. E. Husband, “Colorectal hepatic metastases: quantitative measurements using single-shot echo-planar diffusion-weighted MR imaging,” European Radiology, vol. 16, pp. 1898–1905, May 2006.
  * [8] M. A. Troelstra, A.-M. V. Dijk, J. J. Witjes, A. L. Mak, D. Zwirs, J. H. Runge, J. Verheij, U. H. Beuers, M. Nieuwdorp, A. G. Holleboom, A. J. Nederveen, and O. J. Gurney-Champion, “Self-supervised neural network improves tri-exponential intravoxel incoherent motion model fitting compared to least-squares fitting in non-alcoholic fatty liver disease,” Frontiers in Physiology, vol. 13, Sept. 2022.
  * [9] Y.-H. Nai, X. Wang, J. Gan, C. P. L. Lian, R. F. Kirwan, F. S. L. Tan, and D. J. Hausenloy, “Effects of fitting methods, high b-values and image quality on diffusion and perfusion quantification and reproducibility in the calf,” Computers in Biology and Medicine, vol. 157, p. 106746, May 2023\.
  * [10] E. Scalco, G. Rizzo, and A. Mastropietro, “The quantification of IntraVoxel incoherent motion – MRI maps cannot preserve texture information: An evaluation based on simulated and in-vivo images,” Computers in Biology and Medicine, vol. 154, p. 106495, Mar. 2023.
  * [11] O. J. Gurney-Champion, R. Klaassen, M. Froeling, S. Barbieri, J. Stoker, M. R. W. Engelbrecht, J. W. Wilmink, M. G. Besselink, A. Bel, H. W. M. van Laarhoven, and A. J. Nederveen, “Comparison of six fit algorithms for the intra-voxel incoherent motion model of diffusion-weighted magnetic resonance imaging data of pancreatic cancer patients,” PLOS ONE, vol. 13, p. e0194590, Apr. 2018.
  * [12] N. Liu, X. Yang, L. Lei, K. Pan, Q. Liu, and X. Huang, “Intravoxel incoherent motion model in differentiating the pathological grades of esophageal carcinoma: Comparison of mono-exponential and bi-exponential fit model,” Frontiers in Oncology, vol. 11, Apr. 2021.
  * [13] C. Vieni, B. Ades-Aron, B. Conti, E. E. Sigmund, P. Riviello, T. M. Shepherd, Y. W. Lui, D. S. Novikov, and E. Fieremans, “Effect of intravoxel incoherent motion on diffusion parameters in normal brain,” NeuroImage, vol. 204, p. 116228, Jan. 2020.
  * [14] S. Yamada, S. Hiratsuka, T. Otani, S. Ii, S. Wada, M. Oshima, K. Nozaki, and Y. Watanabe, “Usefulness of intravoxel incoherent motion MRI for visualizing slow cerebrospinal fluid motion,” Fluids and Barriers of the CNS, vol. 20, Mar. 2023.
  * [15] Z.-X. Kuai, W.-Y. Liu, and Y.-M. Zhu, “Effect of multiple perfusion components on pseudo-diffusion coefficient in intravoxel incoherent motion imaging,” Physics in Medicine & Biology, vol. 62, pp. 8197–8209, Oct. 2017.
  * [16] J.-P. Cercueil, J.-M. Petit, S. Nougaret, P. Soyer, A. Fohlen, M.-A. Pierredon-Foulongne, V. Schembri, E. Delhom, S. Schmidt, A. Denys, S. Aho, and B. Guiu, “Intravoxel incoherent motion diffusion-weighted imaging in the liver: comparison of mono-, bi- and tri-exponential modelling at 3.0-t,” European Radiology, vol. 25, pp. 1541–1550, Dec. 2014.
  * [17] M. C. Wurnig, M. Germann, and A. Boss, “Is there evidence for more than two diffusion components in abdominal organs? – a magnetic resonance imaging study in healthy volunteers,” NMR in Biomedicine, vol. 31, Nov. 2017.
  * [18] O. Chevallier, N. Zhou, J.-P. Cercueil, J. He, R. Loffroy, and Y. X. J. Wáng, “Comparison of tri-exponential decay versus bi-exponential decay and full fitting versus segmented fitting for modeling liver intravoxel incoherent motion diffusion MRI,” NMR in Biomedicine, vol. 32, July 2019.
  * [19] A. J. Riexinger, J. Martin, S. Rauh, A. Wetscherek, M. Pistel, T. A. Kuder, A. M. Nagel, M. Uder, B. Hensel, L. Müller, and F. B. Laun, “On the field strength dependence of bi- and triexponential intravoxel incoherent motion (IVIM) parameters in the liver,” Journal of Magnetic Resonance Imaging, vol. 50, pp. 1883–1892, Apr. 2019.
  * [20] A. Riexinger, J. Martin, A. Wetscherek, T. A. Kuder, M. Uder, B. Hensel, and F. B. Laun, “An optimized b-value distribution for triexponential intravoxel incoherent motion (IVIM) in the liver,” Magnetic Resonance in Medicine, vol. 85, pp. 2095–2108, Nov. 2020.
  * [21] D. L. Bihan, R. Turner, C. T. W. Moonen, and J. Pekar, “Imaging of diffusion and microcirculation with gradient sensitization: Design, strategy, and significance,” Journal of Magnetic Resonance Imaging, vol. 1, pp. 7–28, Jan. 1991.
  * [22] J. Pekar, C. T. W. Moonen, and P. C. M. van Zijl, “On the precision of diffusion/perfusion imaging by gradient sensitization,” Magnetic Resonance in Medicine, vol. 23, pp. 122–129, Jan. 1992.
  * [23] D.-M. Koh, M. Blackledge, D. J. Collins, A. R. Padhani, T. Wallace, B. Wilton, N. J. Taylor, J. J. Stirling, R. Sinha, P. Walicke, M. O. Leach, I. Judson, and P. Nathan, “Reproducibility and changes in the apparent diffusion coefficients of solid tumours treated with combretastatin a4 phosphate and bevacizumab in a two-centre phase i clinical trial,” European Radiology, vol. 19, pp. 2728–2738, June 2009.
  * [24] M. Iima, K. Yano, M. Kataoka, M. Umehana, K. Murata, S. Kanao, K. Togashi, and D. L. Bihan, “Quantitative non-gaussian diffusion and intravoxel incoherent motion magnetic resonance imaging,” Investigative Radiology, vol. 50, pp. 205–211, Apr. 2015.
  * [25] Y.-P. Liao, S. ichi Urayama, T. Isa, and H. Fukuyama, “Optimal model mapping for intravoxel incoherent motion MRI,” Frontiers in Human Neuroscience, vol. 15, Feb. 2021.
  * [26] J. A. U. Perucho, H. C. C. Chang, V. Vardhanabhuti, M. Wang, A. S. Becker, M. C. Wurnig, and E. Y. P. Lee, “B-value optimization in the estimation of intravoxel incoherent motion parameters in patients with cervical cancer,” Korean Journal of Radiology, vol. 21, no. 2, p. 218, 2020.
  * [27] A. Lemke, B. Stieltjes, L. R. Schad, and F. B. Laun, “Toward an optimal distribution of b values for intravoxel incoherent motion imaging,” Magnetic Resonance Imaging, vol. 29, pp. 766–776, July 2011.
  * [28] J. Raju, U. A. C., and A. John, “A novel approach for b-value optimization in intravoxel incoherent motion imaging using metaheuristic algorithm,” Expert Systems with Applications, vol. 168, p. 114270, Apr. 2021.
  * [29] C. Paganelli, M. A. Zampini, L. Morelli, G. Buizza, G. Fontana, L. Anemoni, S. Imparato, G. Riva, A. Iannalfi, E. Orlandi, and G. Baroni, “Optimizing b-values schemes for diffusion MRI of the brain with segmented intravoxel incoherent motion (IVIM) model,” Journal of Applied Clinical Medical Physics, vol. 24, Apr. 2023.
  * [30] G. Zhu, J. J. Heit, B. W. Martin, D. G. Marcellus, C. Federau, and M. Wintermark, “Optimized combination of b‑values for IVIM perfusion imaging in acute ischemic stroke patients,” Clinical Neuroradiology, vol. 30, pp. 535–544, Aug. 2019.
  * [31] A. Riexinger, J. Martin, A. Wetscherek, T. A. Kuder, M. Uder, B. Hensel, and F. B. Laun, “An optimized b-value distribution for triexponential intravoxel incoherent motion (IVIM) in the liver,” Magnetic Resonance in Medicine, vol. 85, pp. 2095–2108, Nov. 2020.
  * [32] H.-J. Lee, S. Y. Rha, Y. E. Chung, H. S. Shim, Y. J. Kim, J. Hur, Y. J. Hong, and B. W. Choi, “Tumor perfusion-related parameter of diffusion-weighted magnetic resonance imaging: Correlation with histological microvessel density,” Magnetic Resonance in Medicine, vol. 71, pp. 1554–1558, June 2013.
  * [33] S. Bisdas, C. Braun, M. Skardelly, J. Schittenhelm, T. H. Teo, C. H. Thng, U. Klose, and T. S. Koh, “Correlative assessment of tumor microcirculation using contrast-enhanced perfusion MRI and intravoxel incoherent motion diffusion-weighted MRI: is there a link between them?,” NMR in Biomedicine, vol. 27, pp. 1184–1191, Aug. 2014.
  * [34] R. P. Kennan, J.-H. Gao, J. Zhong, and J. C. Gore, “A general model of microcirculatory blood flow effects in gradient sensitized MRI,” Medical Physics, vol. 21, pp. 539–545, Apr. 1994.
  * [35] A. Wetscherek, B. Stieltjes, and F. B. Laun, “Flow-compensated intravoxel incoherent motion diffusion imaging,” Magnetic Resonance in Medicine, vol. 74, pp. 410–419, Aug. 2014.
  * [36] R. M. Henkelman, J. J. Neil, and Q.-S. Xiang, “A quantitative interpretation of IVIM measurements of vascular perfusion in the rat brain,” Magnetic Resonance in Medicine, vol. 32, pp. 464–469, Oct. 1994.
  * [37] T. Q. Duong and S.-G. Kim, “In vivo MR measurements of regional arterial and venous blood volume fractions in intact rat brain,” Magnetic Resonance in Medicine, vol. 43, pp. 393–402, Mar. 2000.
  * [38] G. Fournet, J.-R. Li, A. M. Cerjanic, B. P. Sutton, L. Ciobanu, and D. L. Bihan, “A two-pool model to describe the IVIM cerebral perfusion,” Journal of Cerebral Blood Flow & Metabolism, vol. 37, pp. 2987–3000, Jan. 2016\.
  * [39] A. D. Cohen, M. C. Schieke, M. D. Hohenwalter, and K. M. Schmainda, “The effect of low b-values on the intravoxel incoherent motion derived pseudodiffusion parameter in liver,” Magnetic Resonance in Medicine, vol. 73, pp. 306–311, Jan. 2014.
  * [40] C. Ye, D. Xu, Y. Qin, L. Wang, R. Wang, W. Li, Z. Kuai, and Y. Zhu, “Estimation of intravoxel incoherent motion parameters using low b-values,” PLOS ONE, vol. 14, p. e0211911, Feb. 2019.
  * [41] P. Mukherjee, J. Berman, S. Chung, C. Hess, and R. Henry, “Diffusion tensor MR imaging and fiber tractography: Theoretic underpinnings,” American Journal of Neuroradiology, vol. 29, pp. 632–641, Mar. 2008.
  * [42] 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.
  * [43] J. A. Lee, X. Geets, V. Grégoire, and A. Bol, “Edge-preserving filtering of images with low photon counts,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 30, no. 6, pp. 1014–1027, 2008.
  * [44] 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.
  * [45] 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.
  * [46] S. Guérit, L. Jacques, B. Macq, and J. A. Lee, “Post-reconstruction deconvolution of pet images by total generalized variation regularization,” 2015\.
  * [47] 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.
  * [48] 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.
  * [49] 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.
  * [50] C. Sample, A. Rahmim, F. Bénard, J. Wu, and H. Clark, “Psma pet as a predictive tool for sub-regional importance estimates in the parotid gland,” 2023\.
  * [51] H. Gudbjartsson and S. Patz, “The rician distribution of noisy mri data,” Magnetic Resonance in Medicine, vol. 34, pp. 910–914, Dec. 1995.
  * [52] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” 2014.
  * [53] 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.
  * [54] J. Wang, K. C. K. Chan, and C. C. Loy, “Exploring clip for assessing the look and feel of images,” arXiv, 2022.
  * [55] M. P. T. Kaandorp, S. Barbieri, R. Klaassen, H. W. M. van Laarhoven, H. Crezee, P. T. While, A. J. Nederveen, and O. J. Gurney-Champion, “Improved unsupervised physics-informed deep learning for intravoxel incoherent motion modeling and evaluation in pancreatic cancer patients,” Magnetic Resonance in Medicine, vol. 86, pp. 2250–2265, June 2021.
  * [56] D.-A. Clevert, T. Unterthiner, and S. Hochreiter, “Fast and accurate deep network learning by exponential linear units (elus),” 2015.
  * [57] S. Barbieri, O. J. Gurney-Champion, R. Klaassen, and H. C. Thoeny, “Deep learning how to fit an intravoxel incoherent motion model to diffusion-weighted MRI,” Magnetic Resonance in Medicine, vol. 83, pp. 312–321, Aug. 2019.
  * [58] M. Bruvo and F. Mahmood, “Apparent diffusion coefficient measurement of the parotid gland parenchyma,” Quantitative Imaging in Medicine and Surgery, vol. 11, pp. 3812–3829, Aug. 2021.
  * [59] A. Nada, A. Youssef, A. Basmy, A. Amin, and A. Shokry, “Diffusion-weighted imaging of the parotid gland: Can the apparent diffusion coefficient discriminate between normal and abnormal parotid gland tissues?,” Erciyes Tıp Dergisi/Erciyes Medical Journal, vol. 39, pp. 125–130, Sept. 2017\.
  * [60] M. Kimura, H. Yabuuchi, H. Narita, Y. Kurihara, K. Hisada, N. Sakai, K. Nagatomo, R. Mikayama, M. Masaki, and H. Kimura, “Intravoxel incoherent motion diffusion-weighted imaging (IVIM-DWI) of the major salivary glands: an assessment of the optimal number and combination of b-values,” Polish Journal of Radiology, vol. 87, no. 1, pp. 246–256, 2022.
  * [61] A. S. Becker, A. Manoliu, M. C. Wurnig, and A. Boss, “Intravoxel incoherent motion imaging measurement of perfusion changes in the parotid gland provoked by gustatory stimulation: A pilot study,” Journal of Magnetic Resonance Imaging, vol. 45, pp. 570–578, July 2016.
  * [62] A. Cieszanowski, K. Pasicz, J. Podgórska, E. Fabiszewska, W. Skrzyński, J. Jasieniak, A. Anysz-Grodzicka, I. Grabska, J. Pałucki, M. Naduk-Ostrowska, B. Jagielska, and P. Kukołowicz, “Reproducibility of intravoxel incoherent motion of liver on a 3.0t scanner: free-breathing and respiratory-triggered sequences acquired with different numbers of excitations,” Polish Journal of Radiology, vol. 83, pp. 437–445, 2018.
  * [63] H. Dyvorne, G. Jajamovich, S. Kakite, B. Kuehn, and B. Taouli, “Intravoxel incoherent motion diffusion imaging of the liver: Optimal b-value subsampling and impact on parameter precision and reproducibility,” European Journal of Radiology, vol. 83, pp. 2109–2113, Dec. 2014.
  * [64] Sample, C. (2023). `neural_blind_deconv_PSMA`. GitHub. https://github.com/samplecm/neural_blind_deconv_PSMA.


Rhombic dodecahedral robots. When bees build honeycombs (A) they fashion three
rhombuses at the bottom of each hexagonal chamber (B), minimizing wax and
tessellating two layers back-to-back. Completing this pattern in 3D (C) and
closing each cell yields rhombic dodecahedra (D), a hitherto nonexistent shape
in robotics. We built dozens of dodecahedral cells with internal power and
actuation (E) that attach to form self-motile superstructures (F).

# A non-cubic space-filling modular robot

Author Names Omitted for Anonymous Review. Paper-ID [add your ID here]  Tyler
Hummer and Sam Kriegman  Northwestern University

###### Abstract

Space-filling building blocks of diverse shape permeate nature at all levels
of organization, from atoms to honeycombs, and have proven useful in
artificial systems, from molecular containers to clay bricks. But, despite the
wide variety of space-filling polyhedra known to mathematics, only the cube
has been explored in robotics. Thus, here we roboticize a non-cubic space-
filling shape: the rhombic dodecahedron. This geometry offers an appealing
alternative to cubes as it greatly simplifies rotational motion of one cell
about the edge of another, and increases the number of neighbors each cell can
communicate with and hold on to. To better understand the challenges and
opportunities of these and other space-filling machines, we manufactured 48
rhombic dodecahedral cells and used them to build various superstructures. We
report locomotive ability of some of the structures we built, and discuss the
dis/advantages of the different designs we tested. We also introduce a
strategy for genderless passive docking of cells that generalizes to any
polyhedra with radially symmetrical faces. Future work will allow the cells to
freely roll/rotate about one another so that they may realize the full
potential of their unique shape.

Figure 1: Exploded view of internal components within an active cell. Active
cells include: a hexagonal PCB (A) as a structural foundation and to create
all electrical connections; an analog ON/OFF switch (B) for simple mode
control; an eccentric rotating mass (ERM) motor (C) to provide mechanical
energy to the morphological system; a 100 mAh LiPo battery cell (D) with high
discharge capabilities; surface mounted light emitting diodes (LEDs; E) to
indicate active cells’ current states; a 3D printed PCB mount (F) to hold
components in place and mount outer shells; M2 heat inserts (G) to mount the
outer shells; two outer shell halves (H) with 48 face-embedded magnets; and M2
screws (I) to complete cell assembly.  Figure 2: The potential for self
reconfiguration in future work. (A and B:) The two kinds of rolling motions
for cubes. (C:) 2D projection of rhombic dodecahedron rotation (hexagon).

##  I Introduction

When Johannes Kepler peered into the hexagonal cells of honeycombs he saw
three rhombi at their base (Fig. A non-cubic space-filling modular robotA,B).
Capping a cell with three more rhombi, Kepler sketched a rhombic dodecahedron.
But the shape was already on his mind. Studying how snowflakes could be built
from microscopic spheres of ice, Kepler made the following conjecture [9]: no
packing of congruent spheres has density greater than that of the face-
centered cubic packing (Fig. A non-cubic space-filling modular robotC), in
which each sphere touches 12 others, and which fills just over 74% of
space.111Kepler’s conjecture was formally proven four centuries later, in 2017
[6]. When uniformly compressed, these spheres will deform into space-filling
polyhedra with twelve faces: rhombic dodecahedra (Fig. A non-cubic space-
filling modular robotD).

Since Kepler’s discovery of the shape in the early 17th century, rhombic
dodecahedra have been found in crystals of viruses [5], proteins [4],
carbohydrates [25], metals [17] and minerals [3]. These self-assembling
crystalline structures have been of keen interest in several fields of
engineering due to their exceptional optical, catalytic, electronic and
magnetic properties [24, 8], which suggests that space filling in general, and
the dodecahedral morphology in particular, may also be useful for self-
assembling modular robots [22].

Like crystals, modular robots can “grow” to form multicellular bodies by
ceaselessly combining and separating unit cells. Such robots—if space-
filling—could form gapless containers, pipes, mirrors, shields, submersibles,
spaceships and extraterrestrial habitats on the fly. And by virtue of broader
surface contacts between neighbors, space-filling robots could build much more
stable and energetically favorable superstructures than non-space-filling
robots [14].

Many unit shapes have been explored in robotics: triangles [2], hexagons [18],
spheres [23], cubes [19], cuboids [1], cylinders [13] semi-cylinders [20] and
cylindroids [7]. But only two of these shapes were built to attach in 3D space
(spheres and cubes), and only one was capable of filling it: cubes.

The cube is a parallelohedron [21]: a congruent polyhedron that can be
translated to completely fill space, without rotations. In other words, space
can be partitioned into a regular lattice by cubes. In addition to the cube,
which has square faces, there are just four other parallelohedra: the rhombic
dodecahedron (rhombic faces), the hexagonal prism (hexagonal and cubic faces),
the elongated dodecahedron (rhombic and octagonal faces), and the truncated
octahedron (cubic, rhombic, and hexagonal faces). The cube and the rhombic
dodecahedron are the simplest parallelohedra. Their isometric faces (all
squares or all rhombi) allow two cells to connect on all sides, which
simplifies the assembly and reconfiguration of multicellular bodies.

There is an infinitude of other polyhedra that, alone or in combination with
additional polyhedra, via translations, rotations or reflections, along
periodic or aperiodic tilings, can fill space without any gaps. Out of all of
these possible shapes, only the cube has been taken seriously in robotics.
Indeed, the vast majority of modular robots built to date have been either
cubes or cuboids [15].

Rhombic dodecahedra have an advantage kinematically over cubes that was first
noted by Yim et al. [26]: Cubes require two kinds of inherently difficult
rolling/rotational motions (90$\degree$ and 180$\degree$; Fig. 2A and B) to
move about each other, and can incur considerable sliding friction against any
cubes perpendicular to the axis of rotation. Rhombic dodecahedra, in contrast,
require only a single simple rotational movement (120°; Fig. 2C), and
eliminate the need for any part of the cells to slide against each other,
while preserving the good property of cubes: Once two opposing faces of two
cells are aligned and connected, all rotations about one of the four shared
edges will result in another pair of aligned faces.

After neatly illustrating this, Yim et al. [26] went on to show how to
coordinate the self-assembly of virtual rhombic dodecahedra toward target
shapes, but, until now, no such robots have been created. Thus here we design
and deploy the first rhombic dodecahedra robots. While they cannot yet
assemble themselves, their physicality helps us envisage a future in which
they do.

##  II Methods

In this section we detail the design and construction of the robot’s cells.

###  II-A Active and passive cells.

We 3D printed two kinds of rhombic dodecahedra: passive cells (purple) and
active cells (white). Both active and passive cells consist of a hollow
plastic dodecahedral shell (Fig. 1H) with four circular indentations on each
face. Inside each indentation a neodymium magnet was inserted and glued in
place (Fig. 3A). Passive cells were printed as a single piece of PLA plastic.

Active cells contain a hand-soldered hexagonal printed-circuit-board (PCB), an
eccentric rotating mass (ERM) vibration motor, and a 3.7V 100mAh rechargeable
LiPo battery mounted to the outer shell, which is printed in two parts and
screwed together. An exploded view of an active module is provided in Fig. 1.
The analog PCB has an ON/OFF toggle switch and LEDs that indicate whether the
module is in the ON or OFF state. When switched ON, the motor rotates at 19000
RPM consuming 195 mA, which allows for a total of $\sim$25 min of continuous
operation (approximately half a million revolutions) before depleting its
rechargeable battery. Each active module costs less than $\$$10 USD to build,
with the most expensive component being the magnets.

###  II-B Docking.

Figure 3: Genderless passive docking was achieved by arranging and embedding
axially-poled neodymium magnets into the outer shell of the cells (A). The
arrangement of north (N) and south (S) poled magnets takes advantage of the
symmetry line that exists along the long axis of the rhombic faces. (B:) A jig
was designed and 3D printed to simplify cell assembly. The jig consists of a
base (B1) to hold and supply magnets, a slider base (B2) with identical magnet
configurations as shown in A that automatically loads itself when slid through
the base, and a slider top (B3) that can be removed to allow all four magnets
in a cell face to be embedded simultaneously.

Genderless passive magnetic docking allows any two cells to dock at any
aligned face, without rotations (Fig. 3A), in a fashion similar to that of
biological cells [10, 11].

Initial prototypes incorporated diametrically-poled neodymium magnets that
could freely rotate into the necessary orientation to enable cell-cell
docking. However, the scale and shape of the cells caused the orientation of
magnets to lock pole-to-pole within the cell, which prevented the magnets from
rotating and thus inhibited the cell from docking with another. To rectify
this, four axially-poled neodymium magnets were embedded on each of the cell’s
12 faces (48 total magnets per cell).

There are many ways to arrange these magnets to ensure a stable connection.
For example, magnets could be placed in the four corners of each face.
However, this would not always allow for two geometrically aligned cells to
connect magnetically. While a new cell could always be rotated to connect to
the end of a 1D chain, in 2D and 3D structures a valid rotation will not
always be possible, rendering connections unstable (one or more pairs of
magnets are repulsive) or impossible (no pairs of magnets are attractive).

To ensure genderless docking, pairs of magnets along the two opposing faces of
two aligned cells must preserve their unlike poles under 180$\degree$
rotations. The magnet configuration we chose respects this constraint (Fig.
3A). With this magnet configuration, cells can connect within the entire face-
centered cubic lattice-configuration that is characteristic of rhombic
dodecahedra (Fig. A non-cubic space-filling modular robotC). This provides
orientation-invariant and configuration-invariant attachment of cells, or
groups of cells (Fig. 4), and is not limited to rhombic dodecahedra. Indeed,
this strategy generalizes to any polyhedra with faces that possess at least
two degrees of rotational symmetry. Examples of such shapes include the
regular icosahedron, which has 20 equilateral triangular faces with three
degrees of rotational symmetry on each face; the gyrobifastigium, which has
four equilateral triangular faces and four square faces (four degrees of
rotational symmetry); and the truncated octahedron, which has eight hexagonal
faces (six degrees of rotational symmetry) and six square faces.

Despite the simplicity that passive magnetic docking enables, complications
still arose. The first was that, due to the magnet’s smooth surface finish,
glue could not properly adhere to their surface, and as a result, embedded
magnets were prone to falling out. This problem was ameliorated by sanding
down the surface finish of the magnets to increase their roughness and surface
friction prior to mounting. This greatly reduced the number of magnets that
dislodged from the face during reconfiguration.

Figure 4: Robotic crystals formed by various numbers (n) of space-filling unit
cells.

Another challenge was that embedding magnets manually, one by one, is a time
consuming and labor intensive task that required upwards of 30 min per cell.
We addressed this problem by designing a jig to hold the magnets and guide
them into place (Fig. 3B). The jig consists of only three parts: a base that
holds and supplies the magnets, a slider that has identical magnet
configurations as the cells, and a removable slider top. The slider is moved
through the slot in the base allowing the magnets to reload themselves into
the desired location. The slider top can then be removed and the entire face
of a cell can be mounted simultaneously. With these refinements, one cell can
be embedded with magnets in 10 min, and is now primarily limited by the curing
time of the glue used for mounting.

Figure 5:  Locomotion. Behavior was tracked for four different designs (A-D).
Across six independent trials, design’s center of mass was tracked over an
evaluation period of 1 min. The initial orientation for each trial—the top-
down view used for motion tracking—as well as a side view of each morphology,
is shown for each respective design. The origin of each trial is marked with a
red triangle indicating the initial position and orientation of the robot. The
robot’s position over the trial period, along with the final position and
orientation, are marked with grayscale lines and corresponding triangles.
Active cells (white) were rotated into six different random orientations, one
for each trial, to understand the influence of motor orientation on behavior.
In smaller body plans (e.g. A), the orientation of active cells were found to
have greater influence on behavior (e.g. the direction of rotation) when
compared to larger bodies (e.g. C and D). Potential correlations between the
type of surface contact for a given design were also identified. For example,
point contacts seem to produce a high rate of rotation (A); edge contacts
produced large, sweeping rotations (B and C); and face contact produced mostly
translational motion (D).

##  III Results

Four unique morphologies comprising active and passive cells were configured
manually and tested for their locomotive abilities (Fig. 5). Six trials were
recorded, each lasting 1 min. At the beginning of each trial, the body was
placed at a marked central location with the same initial orientation. To
isolate the influence of motor orientation from that of body shape and motor
position, the active cells within each design were rotated into a different
random orientation for each trial, changing the orientation and directionality
of the ERM motor within them. Although an exhaustive test of motor
orientations would be ideal, given that each active cell has 12 sides, each of
which can be attached in two ways, testing each motor orientation of each
active cell was not feasible under the given time constraints.

Motion of the bodies was tracked using Blender222www.blender.org and extracted
for analysis.333github.com/Amudtogal/blenderMotionExport Tracking boxes within
Blender were placed along all of the upwards facing cells of the body, and
then joined to track the center of mass. One advantage of the chosen magnet
arrangement is that it offers a plethora of markers for motion tracking
without needing to add sticker tags. The mirror finish of the magnets can
cause issues with light reflections at certain orientation, but the effects of
this can be minimized with a sufficient number of tracking markers.

With these four unique morphologies, four distinct behaviors emerged. Design A
(Fig. 5A) was the smallest morphology we tested, comprising just three cells:
two passive and one active. It was extremely motile, traveling a mean distance
of 126 cm $\pm$ 34 SD over the six trials; however, the motion was largely
rotational rather than translational, resulting in a mean net displacement of
only 6 cm $\pm$ 4 SD. The direction of rotation, clockwise or
counterclockwise, was greatly influenced by the orientation of the motor
within the passive cell. Out of the six trials, 66.67% random rotations of the
single active cell caused Design A to turn clockwise, and in the other 33.33%
it turned counterclockwise.

Design B (Fig. 5B) is also quite small (three passive, one active) and is also
very sensitive to motor orientation. Exactly half of the random rotations
applied to its single active cell caused Design B to turn clockwise, and in
the other half of the trials it turned counterclockwise.

TABLE I: Morphological and behavioral characterization of the four designs we tracked. | Design A | Design B | Design C | Design D
---|---|---|---|---
No. of passive cells | 2 | 3 | 5 | 7
No. of active cells | 1 | 1 | 2 | 3
Ratio of passive to active | 2 to 1 | 3 to 1 | 2.5 to 1 | 2.33 to 1
Body length (cm) | 9.5 | 11.5 | 15 | 21
Body weight (g) | 77 | 98 | 175 | 252
Type of surface contacts | Point | Edge | Edge | Face
Avg. distance traveled (cm) | 126 +/- 34 SD | 76 +/- 14 SD | 106 +/- 22 SD | 79 +/- 26 SD
Avg. net displacement (cm) | 6 +/- 4 SD | 22 +/- 21 SD | 9 +/- 5 SD | 16 +/- 10 SD
No. of trials | 6 | 6 | 6 | 6

Design C (Fig. 5C) comprised five passive cells and two active cells. Similar
to Design A, Design C was highly motile, traveling an average of 106 cm $\pm$
22 SD across the six trials. However, contrary to the inconsistent direction
of rotation exhibited by Design A, Design C consistently traveled in a
counterclockwise direction of rotation in all six of the trials. This may
imply that as the number of cells within a given morphology grows, the
orientation of the motor within the active cells has less influence on the
overall behavior. Thus, the behavior could become increasingly a result of the
position of active and passive cells.

This may also be supported by the mostly translational motion of Design D
(Fig. 5C), which exhibited a mean distance of travel of 79 cm $\pm$ 26 SD, and
a mean net displacement of 16 cm $\pm$ 10 SD. With seven passive and three
active cells, Design D was the largest of the four designs we tested. In the
majority of trials (66.67%) Design D moved translationally.

Another characteristic of interest is the type of surface contact that a given
morphology maintains with the ground plane. Due to the geometrical properties
of the rhombic dodecahedron, all cells with a given morphology have identical
orientations. As a result, each design can have only a single type of surface
contact: point, edge, or face. Each type of contact was represented by at
least one of the tested designs.

Design A maintained point contact with the ground, and also exemplified the
largest rate of consistent rotation, regardless of direction. Designs B and C
both maintained edge contact. Their behavior, as with Design A, was
characteristically rotational, but with larger, sweeping paths. Design D, on
the other hand, maintained face contact with the ground plane and displayed
more consistent translational motion. These results suggest that the type of
surface contact not only influences the emergent behavior of a given design,
but that surface contacts could be intentionally modified to achieve desired
behaviors.

##  IV Discussion

Inspired by crystals across length scales, as well as sheer intellectual
curiosity, we set out to build a non-cubic space-filling robot. We zeroed in
on the rhombic dodecahedron rather quickly because it seemed to us to be the
simplest option. We later discovered that we were not the first to dream of
robots with this shape. Three decades ago, Yim et al. [26] clearly illustrated
the issues with self-assembling cubes and how the rhombic dodecahedron
obviates them. Lost to history, no such robots were ever built, nor have they
been discussed since.

Here we demonstrated the first robots composed of cells with rhombic
dodecahedron shape. At the moment these self-motile machines are more likely
to raise debates than carry out useful work. What must a machine do, in
addition to moving around, to be considered a bona fide robot? Future work
could closely follow the cadence of this argument, squeezing sensing,
computation and communication inside and along the shell of each module. But,
so long as the individual cells cannot intentionally detach, roll themselves
over 120$\degree$ and reattach, these robots will not be able to do the very
thing that, at least according to Yim et al. [26], makes them so special.

To get the cell rolling, internal mechanisms already used by self-assembling
cubes, such as a flywheel and brakes [19], could be adapted for six axes of
rotation instead of three. But this will likely require a considerable
investment in terms of time, money and effort. Future cells could also be
built with materials that support active de/magnetization and respond to
external forces such as light and rotating magnetic fields [12] to allow the
cells to passively roll up and down the robot’s body. Existing computationally
efficient topology search algorithms [16] could be adapted for internally
modeling the repercussions of such actions prior to reconfiguration.

That said, we would like to emphasize the value of designing, building and
playing with this first generation of non-rolling dodecahedral cells. The
first hand experience of manipulating physical rhombic dodecahedra—manually
rotating and snapping them together in different arrangements—helped us better
understand the mathematics of crystals and fully grasp the kinematics of
robots with dodecahedral morphology. We also found playing with them to be
quite fun. They provide a satisfying tactile experience due to their symmetry,
obtuse angles which let them roll in your hand, and the building tension of
twelve magnetic faces. This suggests that they may provide an engaging and
interactive educational tool for learning about geometry, crystals and robots.

## Acknowledgements

We thank Elaine Liu for help designing early prototypes, Davin Landry for help
with printing and process development for cell assembly, and Stefan Knapik,
David Matthews, and other members of the Center for Robotics and Biosystems at
Northwestern University for helpful discussions. This research was supported
by a seed grant from the Center for Engineering Sustainability and Resilience
at Northwestern University.

## References

  * Baca et al. [2014] José Baca, SGM Hossain, Prithviraj Dasgupta, Carl A Nelson, and Ayan Dutta.  Modred: Hardware design and reconfiguration planning for a high dexterity modular self-reconfigurable robot for extra-terrestrial exploration.  _Robotics and Autonomous Systems_ , 62(7):1002–1015, 2014.
  * Belke et al. [2023] Christoph H Belke, Kevin Holdcroft, Alexander Sigrist, and Jamie Paik.  Morphological flexibility in robotic systems through physical polygon meshing.  _Nature Machine Intelligence_ , pages 1–7, 2023.
  * Deer et al. [2013] William A Deer, Robert A Howie, and Jack Zussman.  _An Introduction to the Rock-Forming Minerals_.  Mineralogical Society of Great Britain and Ireland, 2013.
  * Ford et al. [1984] GC Ford, Pauline May Harrison, DW Rice, JMA Smith, A Treffry, JL White, and J Yariv.  Ferritin: design and formation of an iron-storage molecule.  _Philosophical Transactions of the Royal Society of London. B, Biological Sciences_ , 304(1121):551–565, 1984.
  * Glaser [1945] Otto Glaser.  Genes, viruses, and orderly aggregation: A reconnaissance.  _American Scientist_ , 33(3):175–188, 1945.
  * Hales et al. [2017] Thomas Hales, Mark Adams, Gertrud Bauer, Tat Dat Dang, John Harrison, Hoang Le Truong, Cezary Kaliszyk, Victor Magron, Sean McLaughlin, Tat Thang Nguyen, et al.  A formal proof of the Kepler conjecture.  In _Forum of mathematics, Pi_ , volume 5, page e2, 2017.
  * Hirose et al. [1996] Shigeo Hirose, Takaya Shirasu, and Edwardo F Fukushima.  Proposal for cooperative robot “gunryu” composed of autonomous segments.  _Robotics and Autonomous Systems_ , 17(1-2):107–118, 1996.
  * Huang et al. [2012] Wan-Chen Huang, Lian-Ming Lyu, Yu-Chen Yang, and Michael H Huang.  Synthesis of Cu2O nanocrystals from cubic to rhombic dodecahedral structures and their comparative photocatalytic activity.  _Journal of the American Chemical Society_ , 134(2):1261–1267, 2012.
  * Kepler [1611] Johannes Kepler.  _Strena seu de niue sexangula_.  1611\.
  * Kriegman et al. [2020] Sam Kriegman, Douglas Blackiston, Michael Levin, and Josh Bongard.  A scalable pipeline for designing reconfigurable organisms.  _Proceedings of the National Academy of Sciences_ , 117(4):1853–1859, 2020.
  * Kriegman et al. [2021] Sam Kriegman, Douglas Blackiston, Michael Levin, and Josh Bongard.  Kinematic self-replication in reconfigurable organisms.  _Proceedings of the National Academy of Sciences_ , 118(49):e2112672118, 2021.
  * Li et al. [2020] Chuang Li, Garrett C Lau, Hang Yuan, Aaveg Aggarwal, Victor Lopez Dominguez, Shuangping Liu, Hiroaki Sai, Liam C Palmer, Nicholas A Sather, Tyler J Pearson, et al.  Fast and programmable locomotion of hydrogel-metal hybrids under light and magnetic fields.  _Science Robotics_ , 5(49):eabb9822, 2020.
  * Li et al. [2019] Shuguang Li, Richa Batra, David Brown, Hyun-Dong Chang, Nikhil Ranganathan, Chuck Hoberman, Daniela Rus, and Hod Lipson.  Particle robotics based on statistical mechanics of loosely coupled components.  _Nature_ , 567(7748):361–365, 2019.
  * Liao et al. [2013] Ching-Wen Liao, Yeh-Sheng Lin, Kaushik Chanda, Yen-Fang Song, and Michael H Huang.  Formation of diverse supercrystals from self-assembly of a variety of polyhedral gold nanocrystals.  _Journal of the American Chemical Society_ , 135(7):2684–2693, 2013.
  * Liu et al. [2016] Jinguo Liu, Xin Zhang, and Guangbo Hao.  Survey on research and development of reconfigurable modular robots.  _Advances in Mechanical Engineering_ , 8(8):1687814016659597, 2016.
  * Matthews et al. [2023] David Matthews, Andrew Spielberg, Daniela Rus, Sam Kriegman, and Josh Bongard.  Efficient automatic design of robots.  _Proceedings of the National Academy of Sciences_ , 120(41):e2305180120, 2023.
  * Niu et al. [2009] Wenxin Niu, Shanliang Zheng, Dawei Wang, Xiaoqing Liu, Haijuan Li, Shuang Han, Jiuan Chen, Zhiyong Tang, and Guobao Xu.  Selective synthesis of single-crystalline rhombic dodecahedral, octahedral, and cubic gold nanocrystals.  _Journal of the American Chemical Society_ , 131(2):697–703, 2009.
  * Pamecha et al. [1996] Amit Pamecha, Chih-Jung Chiang, David Stein, and Gregory Chirikjian.  Design and implementation of metamorphic robots.  In _International Design Engineering Technical Conferences and Computers and Information in Engineering Conference_ , volume 97577, 1996.
  * Romanishin et al. [2015] John W Romanishin, Kyle Gilpin, Sebastian Claici, and Daniela Rus.  3D M-Blocks: Self-reconfiguring robots capable of locomotion via pivoting in three dimensions.  In _2015 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 1925–1932, 2015.
  * Salemi et al. [2006] Behnam Salemi, Mark Moll, and Wei-Min Shen.  Superbot: A deployable, multi-functional, and modular self-reconfigurable robotic system.  In _International Conference on Intelligent Robots and Systems_ , pages 3636–3641, 2006.
  * Senechal and Galiulin [1984] Marjorie Senechal and Ravil Vagizovich Galiulin.  An introduction to the theory of figures: The geometry of ES Fedorov.  _Structural Topology 1984 núm 10_ , 1984.
  * Seo et al. [2019] Jungwon Seo, Jamie Paik, and Mark Yim.  Modular reconfigurable robotics.  _Annual Review of Control, Robotics, and Autonomous Systems_ , 2:63–88, 2019.
  * Swissler and Rubenstein [2020] Petras Swissler and Michael Rubenstein.  FireAnt3D: a 3D self-climbing robot towards non-latticed robotic self-assembly.  In _International Conference on Intelligent Robots and Systems (IROS)_ , pages 3340–3347, 2020.
  * Wu et al. [2011] Hsin-Lun Wu, Huei-Ru Tsai, Yun-Ting Hung, Ka-Un Lao, Ching-Wen Liao, Pei-Ju Chung, Jer-Shing Huang, I-Chia Chen, and Michael H Huang.  A comparative study of gold nanocubes, octahedra, and rhombic dodecahedra as highly sensitive sers substrates.  _Inorganic chemistry_ , 50(17):8106–8111, 2011\.
  * Yang et al. [2017] Shenyu Yang, Yun Yan, Jianbin Huang, Andrei V Petukhov, Loes MJ Kroon-Batenburg, Markus Drechsler, Chengcheng Zhou, Mei Tu, Steve Granick, and Lingxiang Jiang.  Giant capsids from lattice self-assembly of cyclodextrin complexes.  _Nature Communications_ , 8(1):15856, 2017.
  * Yim et al. [1997] Mark Yim, John Lamping, Eric Mao, and J Geoffrey Chase.  Rhombic dodecahedron shape for self-assembling robots.  _SPL TechReport P_ , 9710777:1997, 1997.