# A direct solution to the interpolative inverse non-uniform fast fourier
transform (iNFFT) problem, for spectral analyses of non-equidistant time-
series data

Michael Sorochan Armstrong, José Carlos Pérez-Girón, José Camacho,Regino
Zamora Michael Sorochan Armstrong and José Camacho are part of the
Computational Data Science (CoDaS) Lab in the Department of Signal Theory,
Telematics, and Communications at the University of Granada.
E-mail: mdarmstr (at) go.ugr.es José Carlos Pérez-Girón and Regino Zamora are
part of the Inter-university Institute for Research on the Earth System in
Andalucia through the University of Granada.

###### Abstract

A simple least-squares optimisation enables the determination of the spectrum
for irregularly sampled data that is readily reconstructed using an adjoint
transformation of the Non-Uniform Fast Fourier Transform (NFFT). This is an
improvement upon previously reported iterative methods for such problems, and
is competitive in terms of time complexity with more recently proposed direct
NFFT inversions when considering comparable matrix pre-computation steps. The
software is highly portable, and available as a convenient Python package
using standard libraries. Given its mathematical simplicity however, it can be
conveniently implemented on any platform.

###### Index Terms:

Sensors, Digital processing, Non-uniformity correction, Remote monitor

##  1 Introduction

Remote sensor data is streamed from regions that are difficult to access for
conventional measurements. Due to interruptions in either network
connectivity, power or other environmental conditions there may be significant
periods of time during which data acquisition is interrupted. This may not
pose problems for ”binned” time-series analysis of the data, if the
acquisition rate of the data is relatively high compared with the expected
duration of a typical outage.

Binning data does partially resolve the spatio-temporal lag of phenomena
measured across different points for the same event which is assumed to take
place within a pre-defined period of time. However regardless of whether or
not missing data is present, this does not make use of the full acquisition
rate of the instrument, which negates some benefits of remote measurement.
Instead, it is possible to perform an analysis of the spectrum or frequency
domain information of the remote measurements and reconstruct the original
data based on parallel analysis of similar remote stations.

Spectral analysis of irregularly sampled data is an active area of research,
and is often referenced broadly as a least-squares analysis, owing to a need
for numerical estimation. This type of analysis commonly finds applications
wherever there is an interruption in measurement - as is commonly seen in
climate applications [1], and other environmental monitoring [2].

###  1.1 Fourier transforms of irregularly sampled data

Non-Uniform Discrete Fast Fourier Transfroms (NFFTs) are a family of
techniques used to analyse irregularly-sampled time series data, which depart
from the traditional assumptions of the Discrete Fourier Transform, (DFTs) and
the related Fast Transforms (FFTs) [3] which assume that the time domain data
are of consistent density. In the case of a time series signal
$f_{j}\in\mathbf{R}$ where $j/M\in[-\frac{1}{2},\frac{1}{2})$, the DFT is
calculated by evaluating the trigonometric polynomial (i.e. frequency)
components, $h_{k}\in\mathbf{C}$ where $k=-N/2,...,N/2-1$:

$h_{k}=\sum_{j=-M/2}^{M/2-1}f_{j}e^{-2\pi\mathbf{i}k\frac{j}{M}}$ (1)

and its inversion evaluates $f_{j}$ to yield the Fourier coefficients:

$f_{j}=\frac{1}{M}\sum_{k=-N/2}^{N/2-1}f_{k}e^{2\pi\mathbf{i}\frac{j}{M}k}.$
(2)

When the number of Fourier coefficients is less than the number of
observations, (i.e. $N<M$), the equality in Equation 2 holds only up to an
approximation. The convention is sometimes reversed [4, 5] by the signage of
the exponential term and the normalisation coefficient. Equations 1 and 2 can
be written as a linear transformation $f_{k}=Ff_{j}$, where from Equations 1
and 2:

$\displaystyle F=\left(e^{-2\pi\mathbf{i}k\frac{j}{M}}\right)_{j,k\in M,N}$
(3) $\displaystyle
F^{-1}=\frac{1}{M}\left(e^{2\pi\mathbf{i}\frac{j}{M}k}\right)_{j,k\in M,N}.$
(4)

A full spectrum of Fourier coefficients in $N$ can be calculated up to $M=N$
where the maximum frequency of $N/2$ is the Nyquist frequency; beyond which
the sampling rate is not sufficient to warrant calculation of further
coefficients. A full Fourier representation is theoretically sufficient to
reconstruct the original data, $f_{j}$, in its entirety; however there are
some notable exceptions in cases of discontinuous data that are not readily
represented by a series of trigonometric polynomials [6]. Since these matrices
are circulant, the operations: $h_{k}=Ff_{j}$ and $f_{j}=F^{-1}h_{k}$ can be
performed in $\mathcal{O}(MlogM)$ time complexity using a “full” Fast Discrete
Fourier Transform (here abreviated as FFT), which can divide the operation
into smaller operations assuming constant spectral density via the $j/M$ term.
Even in cases where the input data is continuous, but a fewer number of
Fourier coefficients are calculated, a “truncated” Fourier transform performs
the operation in $\mathcal{O}(NlogM)$ time [7] [8].

###  1.2 Non-Uniform Fast Fourier Transforms

It is possible to adapt the linear transformations in Equations 3 for
irregularly sampled data in the time domain:
$f(x_{j}),x_{j}/M\in[-\frac{1}{2},\frac{1}{2})$ or irregularly sampled data in
the frequency domain: $h(x_{k}),x_{k}/N\in[-\frac{1}{2},\frac{1}{2})$ as a
Non-uniform Discrete Fourier Transform (NDFT). However for an inconsistent
sampling rate, the assumptions of the FFT break down and these transformations
would otherwise scale relative to a naïeve matrix multiplication in
$\mathcal{O}(N^{3})$ time via the following transformation matrices:

$A_{I}=\left(e^{-2\pi\mathbf{i}k\frac{x_{j}}{M}}\right)_{j,k\in M,N}$ (5)
$A_{II}=\left(e^{-2\pi\mathbf{i}x_{k}\frac{j}{M}}\right)_{j,k\in M,N}$ (6)
$A_{III}=\left(e^{-2\pi\mathbf{i}x_{k}\frac{x_{j}}{M}}\right)_{j,k\in M,N}.$
(7)

Equations 5, 6, and 7 refer to non-uniform transforms of types: $I$
(irregularly-sampled time-series data; equispaced frequency-domain data), $II$
(equispaced time-series data, discontinuous frequency-domain data), and $III$
(discontinuous time and frequency-domain data) which are equivalent to the
transform:

$Af(x_{j})=\hat{\eta}(x_{k})$ (8)

which evaluates $f(x_{j})$ at the points $x_{j}$ to yield the Fourier
coefficients $\hat{\eta}(x_{k})$ evaluated at the frequencies $x_{k}$. $x_{j}$
is non-equispaced, and $k$ is equispaced in the case of a type-I transform,
and for a type-II transform $x_{k}$ and $j$ follow the opposite nomenclature.
Equispaced measurements in time or frequency are abbreviated as $f_{j}$ or
$h_{k}$ respectively. The use of Non-Uniform Discrete Fast Fourier Transforms
(NUDFTs, or NFFTs) complete the operation in
$\mathcal{O}(NlogN+|log(1/\epsilon)|M)$ time complexity, by approximating the
discontinuities as continuous via a window function; although from this it
follows that Equation 8 is not exact. We define a forward NDFT as one that
transforms a time-series signal into the frequency domain111The direction of
this transformation utilises opposite signage to the operations described in
the literature by Potts et al [9, 10]; the definition proposed herein is more
consistent with a more common intuition of the Fourier transform. The
corresponding adjoint transformation is defined by a Hermitian of the forward:

$A^{H}h_{k}=\hat{\phi}(x_{j}).$ (9)

The normalisation of a forward, or adjoint transform is dependent on the
inverse of its transformation matrix as $AA^{H}$ (and so cannot be normalized
by $\frac{1}{N}$ \- see Appendix A) and is not calculated as part of a
standard transform. And so the Hermitian of the forward transformation matrix
(the adjoint transform) is not also equal to its inverse, as would be the case
for equidistant time domain data.

###  1.3 Inverse NFFTs

A fundamental property of the DFT is its self-inversion property; that is,
$F^{H}=F^{-1}$, or $FF^{H}=MI_{N}$. This identity implies orthogonality. But
in the case of NDFTs, $AA^{H}\neq MI_{N}$, since $A$ is not directly
invertible via a constant scaling factor of its Hermitian (Appendix A).

Standard | Inverse
---|---
| Forward NFFT
---
$Af(x_{j})=\hat{\eta}_{k}$
| Inverse NFFT
---
$\hat{f}(x_{j})=A^{-1}h_{k}$
| Adjoint NFFT
---
$A^{H}h_{k}=\hat{\phi}(x_{j})$
| Inverse adjoint NFFT
---
$\hat{h}_{k}=\left(A^{H}\right)^{-1}f(x_{j})$
TABLE I: Type I NFFTs: Standard, Inverse, Adjoint, and Inverse Adjoint. Note that in all cases, the time domain is discontinuous and the frequency domain is continuous. Standard | Inverse
---|---
| Forward NFFT
---
$Af_{j}=\hat{\eta}(x_{k})$
| Inverse NFFT
---
$\hat{f}_{j}=A^{-1}h(x_{k})$
| Adjoint NFFT
---
$A^{H}h(x_{k})=\hat{\phi}_{j}$
| Inverse adjoint NFFT
---
$\hat{h}(x_{k})=\left(A^{H}\right)^{-1}f_{j}$
TABLE II: Type II NFFTs: Standard, Inverse, Adjoint, and Inverse Adjoint. Note
that in all cases, the frequency domain is discontinuous and the time domain
is continuous.

Inverting $A$ directly via a pseudo-inverse is computationally inefficient,
relative to the standard set by an inverse Fourier transform in the
equidistant case. Rather, it is more convenient to find a series of Fourier
coefficients such that the original time series data is recoverable via the
adjoint transform. This can be done in the least squares sense by minimising:

$\operatorname*{arg\,min}_{\hat{h}_{k}}||f(x_{j})-A^{H}\hat{h}_{k}||_{W}$ (10)

such that an invertible $\hat{h}_{k}$ is sought such that:

$A^{H}\hat{h}_{k}=f(x_{j})$ (11)

The norm: $||g||_{W}$ relates the inner product of the function $g$ multiplied
by the weight function $W$ as an $M\times M$ diagonal matrix as:
$||g||_{W}=(g*W)^{H}(g*W)$. The weights function in this instance is implied
to operate in the time domain, which from the convolution of:
$A(f(x_{j})*W)=\hat{h}_{k}\hat{W}$.

A number of proposals to solve this optimisation problem, and others have been
put forward, but the method varies depending on the condition of the inverse
problem, especially for direct solutions [9]. Tables I and II are provided to
clear any potential confusion, for transformations of types I & II. For the
purposes of this work and to explicate the far more common scenario where
irregularly sampled time-series data are provided, the focus will be
exclusively on type-I transforms, although it would be relatively
straightforward to apply the same thinking to type-II transforms.

###  1.4 Optimal interpolation problem

In instances where a function that approximates an irregularly sampled dataset
is sought, and such a function is used to interpolate between the gaps in the
observed data, the smallest possible magnitude of of $N$ Fourier coefficients
is calculated via by a weight function, $\hat{w}_{k}$ that also satisfies
Equation 10:

$\operatorname*{arg\,min}_{\hat{f}_{k}|A^{H}\hat{f}_{k}=f(x_{j})}{\sum_{k\in
I_{N}}\frac{|\hat{f}_{k}|^{2}}{\hat{w}_{k}}}$ (12)

The optimal interpolation problem has been solved earlier using the Conjugate
Gradient Method that minimises the successive error at each iteration, rather
than the residuals relative to the original data [10].

This work demonstrates the quadratic form of Equation 12, and offers a
convenient method for solving the interpolative iNFFT without the need to
iterate. In addition, it avoids problems with dimensionality by inverting the
$AA^{H}$ matrix, which sacrifices time complexity for simplicity; this is an
acceptable trade-off for small values of $N$ and competitive with faster
algorithms for the direct inversion of non-uniform fast Fourier transforms
which require similar $\mathcal{O}(N^{2})$ matrix pre-computations[9].

###  1.5 Admissible weight functions, $\hat{W}$

Weight functions are used to represent the weight matrix $\hat{W}$ along the
diagonal as a high-pass filter. These filters prevent the Fourier coefficients
from fitting a discontinuity in the observed data. Admissible weight functions
are defined by Kunis and Potts [10] as any function $g\in V$ that satisfies:

$\displaystyle|g^{\beta-1}|_{V}=\int_{\mathbf{R}}|dg^{\beta-1}(z)|$ (13)
$\displaystyle=sup\sum_{j=0}^{n-1}|g^{\beta-1}(z_{j+1})-g^{\beta-1}(z_{j})|<\infty$
(14)

i.e. that for $|z|<\frac{1}{2}$ the maximum bounded set of the absolute
difference between $g^{\beta-1}(z_{j+1})$ and $g^{\beta-1}(z_{j})$ does not
diverge when the $\beta^{th}$-fold derivative is taken. The authors further
proposed the practical requirement that all admissible weight functions be
normalised by their $L1$ norm to 1. Two possible kernels that satisfy this
condition are the modified Fejer and Sobolev kernels, where the Fejer kernel
is a kernel of order $\beta=2$ of polynomial order $N$:

$B_{2,N}(x)=\frac{2(1+e^{-2\pi\mathbf{i}x})}{N^{2}}\left(\frac{\sin(\frac{N}{2}\pi
x}{\sin(\pi x)}\right)^{2}$ (15)

Defining the norm of the Kernel function as the convolution of $g*f_{j}$ via
the adjoint transformation of $g$:

$||g||_{1,N}=\sum_{k=-N/2}^{k=N/2}g\left(\frac{k}{N}\right)e^{2\pi\mathbf{i}kx}$
(16)

can be used to define the Sobolev kernel, first as an intermediate and
generalisable weight function:

$g_{\alpha,\beta,\gamma}(z):=c_{\alpha,\beta,\gamma}\frac{\left(\frac{1}{4}-z^{2}\right)^{2}}{\gamma+|z|^{2\alpha}}$
(17)

which can be used recursively to solve for the constant:
$c_{\alpha,\beta,\gamma}$ by taking the inverse of the sum of the right-most
term. The Sobolev kernel is defined as:

$S_{\alpha,\beta,\gamma,N}:=\frac{1+e^{-2\pi\mathbf{i}x}}{2||g_{\alpha,\beta,\gamma}||_{1,N}}\sum_{k=-N/2}^{N/2}g\left(\frac{k}{N}\right)e^{2\pi\mathbf{i}kx}$
(18)

for some values of $\alpha,\beta,\gamma$. Selection of $\alpha$ and $\beta$
are related to the dimensionality of the transform and conditions for a
suitable weight function, but a useful value for $\gamma$ is $10^{-2}$ for
1-dimensional transforms as suggested previously by Kunis and Potts [10] (See
also: Figure 1).

Figure 1: Various values for $\gamma$ and their impact on the Sobolev kernel
function. For lower values of $\gamma$, the filter weights higher frequencies
more strongly.

##  2 Optimization of the optimal interpolation problem

###  2.1 With respect to the predicted Fourier coefficients,
$\mathbf{\hat{f}}$

Let $||\hat{f}||_{\hat{W}^{-1}}$ be the weighted norm of the predicted Fourier
coefficients, $\hat{f}$ and let $||A^{H}\hat{f}-y||_{F}^{2}$ be the
constraint. Let $y=f(x_{j})$, as the provided discontinuous data in the time
domain. The optimisation problem can be described as:

$\operatorname*{arg\,min}_{\hat{f}}\left(||\hat{f}||_{\hat{W}^{-1}}+||A^{H}\hat{f}-y||_{F}^{2}\right)$
(19)

which through their respective norms is equivalent to:

$\operatorname*{arg\,min}_{\hat{f}}\left(\hat{f}^{H}\hat{W}^{-1}\hat{f}+\left(A^{H}\hat{f}-y\right)^{H}\left(A^{H}\hat{f}-y\right)\right)$
(20)

which unpacked as a polynomial expression:

$\hat{f}^{H}\hat{W}^{-1}\hat{f}+\hat{f}^{H}AA^{H}\hat{f}-\hat{f}^{H}Ay-y^{H}A^{H}\hat{f}+y^{H}y.$
(21)

To find the optimal $\hat{f}$ to minimise Equation 20, let the function $\phi$
represent the cost function for the following derivatives (assuming
$\phi=tr\phi(\hat{f})$):

$\displaystyle\frac{\partial\phi}{\partial\hat{f}}\hat{f}^{H}\hat{W}^{-1}\hat{f}=\hat{W}^{-1}\hat{f}+\hat{f}^{H}\hat{W}^{-1}$
(22) $\displaystyle=2\hat{W^{-1}}\hat{f}$ (23)
$\displaystyle\frac{\partial\phi}{\partial\hat{f}}\hat{f}^{H}A^{H}A\hat{f}=2AA^{H}\hat{f}$
(24) $\displaystyle\frac{\partial\phi}{\partial\hat{f}}\hat{f}^{H}Ay=Ay$ (25)
$\displaystyle\frac{\partial\phi}{\partial\hat{f}}y^{H}A^{H}\hat{f}=Ay$ (26)
$\displaystyle\frac{\partial\phi}{\partial\hat{f}}y^{H}y=0$ (27)

Taking $\partial\phi/\partial\hat{f}=0$:

$\displaystyle 2\hat{W^{-1}}\hat{f}+2AA^{H}\hat{f}-2Ay=0$ (28)
$\displaystyle\hat{W^{-1}}\hat{f}+AA^{H}\hat{f}=Ay$ (29)
$\displaystyle\hat{f}=\frac{Ay}{\hat{W^{-1}}+AA^{H}}$ (30)

where in Equation 31 the identity $AA^{H}=MI_{N}$ holds for equispaced
observations as demonstrated in Appendix A. Assuming this identity, even as an
approximation in cases were irregular sampling is rare - a solution for
$\hat{f}$ can be found by bringing the $\hat{W}^{-1}$ out of the inverse term.

####  2.1.1 Optimal interpolation problem for equispaced nodes

Factoring out the constant $M$:

$\hat{f}=\frac{Ay}{M\left(\frac{1}{M}\hat{W^{-1}}+I\right)}$ (31)

it is possible to use one of the Searle identities[11]:
$(I+X^{-1})^{-1}=X(X+I)^{-1}$ so that Equation 31 becomes the simplified
expression for the gradient as a function of the weights:

$\nabla\hat{f}_{\hat{W}}=\frac{Ay\hat{W}}{M\hat{W}+I_{N}}$ (32)

since $\frac{1}{M}\hat{W}^{-1}=M\hat{W}$ for a diagonal matrix.

####  2.1.2 Optimal interpolation problem for non-equispaced nodes

When $f_{j}\in[-\frac{1}{2},\frac{1}{2})$ is not equispaced, the identity:
$AA^{H}=MI_{N}$ does not hold. The most straightforward way of solving for
$\hat{f}$ in Equation 28 is via a numerical approximation of $AA^{H}$; however
warranting further consideration is the fact that it exists as part of a total
matrix inverse as: $(\hat{W}^{-1}+AA^{H})^{-1}$. A direct inverse of this
expression is not numerically stable, since
$\lim_{N/2\rightarrow|N|}\hat{W}_{n}=0$,
$\lim_{N/2\rightarrow|N|}\hat{W}_{n}^{-1}=\infty$. We are then motivated to
find a more stable numerical representation of Equation 28 via the Kailath
identity [12, 11]:

$\left(X+YZ\right)^{-1}=X^{-1}-X^{-1}Y\left(I+ZX^{-1}Y\right)^{-1}ZX^{-1}$
(33)

in which we can substitute $Y,Z$ in Equation 33 via the upper and lower
triangular matrices of $AA^{H}$ via the $LU$ decomposition:

$\displaystyle AA^{H}=LU$ (34)
$\displaystyle\left(\hat{W}^{-1}+LU\right)^{-1}=\hat{W}-\hat{W}L\left(I+U\hat{W}L\right)^{-1}\hat{W}$
(35)

substituting Equation 34 into the denominator of Equation 28, the function
that minimises the cost function for the optimal interpolation problem as a
function of $\hat{f}$ in Equation 19:

$\displaystyle\hat{f}=Ay\left(\hat{W}-\hat{W}L\left(I+U\hat{W}L\right)^{-1}\hat{W}\right)$
(36)

where it is clear that in Equation 36 we avoid the problem of directly
inverting the kernel function along the diagonal of $\hat{W}$, as a substitute
for Equation 28.

###  2.2 Second derivative test - equispaced case

From 28, the second derivatives with respect to $\hat{f}$ are as follows:

$\displaystyle\frac{\partial^{2}\phi}{\partial\hat{f}^{2}}\hat{f}^{H}\hat{W}^{-1}\hat{f}=\hat{W}^{-1}$
(37)
$\displaystyle\frac{\partial^{2}\phi}{\partial\hat{f}^{2}}\hat{f}^{H}AA^{H}\hat{f}=MI_{N}$
(38)

since the inverse of any kernel function
$\lim_{N/2\rightarrow|N|}\hat{W}_{n}=0$,
$\lim_{N/2\rightarrow|N|}\hat{W}_{n}^{-1}=\infty$, we are motivated to reduce
the form:

$\displaystyle\nabla^{2}\hat{f}=0$ (39) $\displaystyle-\hat{W}^{-1}=MI_{N}$
(40) $\displaystyle M\hat{W}+I_{N}=\nabla^{2}\hat{f}$ (41)

$M\hat{W}+I_{N}$ is a constant, diagonal matrix that is positive semi-definite
and indicates that the step direction minimises Equation 20. For the non-
equispaced case, the positive definite-ness of the solution has not been
proven, but will nonetheless be assumed.

##  3 Materials and Methods

####  3.0.1 Remote Sensor Data

Soil temperature data were collected in the Sierra Nevada National Park in
Andalucía, Spain, as part of the Smart EcoMountains program through the
Lifewatch ERIC Program, from solar-powered remote multi-parametric micro-
stations located within different environments (shrub, rock, and bare-earth)
and altitudes (from 1500m to 2300m) from 29-07-2020 until 04-05-2023.
Measurements were continuously collected at a nominal 10 minute interval with
micro-sensors (ibuttons, thermochron-type; Maxim integrated) connected to a
network that transmits the collected information when connectivity is
available. All the collected data were visually inspected for extreme values
outside the normal ranges of the variable in the time series due to sensor
failures, identifying a single sensor reporting erroneous values that were
reported as NaNs.

####  3.0.2 Computational information

All calculations were run with a B550 AORUS ELITE AX V2 motherboard with 128
GiB DDR4 RAM, an AMD Ryzen 9 5950x 32 core processor with an NVIDIA GeForce
RTX 3080 graphics processor, and a 2 TB NVMe Samsung SSD 980. Software was
installed on Ubuntu 22.04.3 LTS 64-bit operating system “Jammy Jellyfish”.

Graphical inspections and data cleaning were conducted with the tidyverse
2.0.0 package [13] implemented in the R 4.3.1 software environment [14].

Calculations were performed using Python 3.10. Data was parsed and cleaned
from the original .csv values using Pandas 2.0.3. Various numerical operations
were performed using the vectorized functions available in NumPy 1.25 [15],
with additional statistical functionality from SciPy 1.11 [16]. Graphics were
created using Matplotlib 3.7.1 [17].

##  4 Results and Discussion

Since it could not be assumed that $AA^{H}=MI_{N}$ due to the obvious gaps in
sensor data, the solution to the inverse problem via Equation 36 was used for
these data. In order to validate the solution using previously unconsidered
data, measurements were randomly removed from a largely continuous dataset.
The interpolative iNFFT was trained on the remaining data, before the data
that was previously removed was used as a test set to evaluate the agreement
of the interpolated function with the observational data. In the second test,
a contiguous block of data was removed and a similar procedure was followed.
An example the interpolation performance over a missing contiguous block can
be seen in Figure 2, where a large number values failed a QC check over a
period of several months.

All training data were mean-centred prior to analysis; and the means were
added to the reconstructed data prior to recording the residuals. The results
summarise cross validation routines for randomly sampled data in Table III,
and randomly sampled contiguous blocks in Table III at 3 levels of missing
data (10, 20, and 30%) performed for a total of 7 times to get the average
error, and the uncertainty, as one standard deviation associated with the
error measurement. Both relative error for the predicted set, and the
correlation coefficient for the reconstructed data on the full data set,
versus the partially sampled dataset were also performed to intuit the error
due to the constraint of the reconstruction (via the weighted kernel
approximation) versus the sparsity of the observed data.

Figure 2: Results of the interpolative iNFFT algorithm using the Sobolev kernel with $\gamma=1e-2$, and $N=1024$. The mean absolute fractional error: $\frac{1}{M}\sum_{j=0}^{M}|(y_{pred}-y_{obs})/y_{obs}|$ is $6.76\times 10^{-2}$ for the iNFFT reconstruction, versus $7.48\times 10^{-2}$ for the equivalent weighted truncated inverse FFT (iFFT) reconstruction. Note the stronger tendency for the iFFT to regress towards the mean for larger gaps. | $r^{iNFFT}_{|rand}$ | $r^{iFFT}_{|rand}$
---|---|---
10 | 1.00E+00 $\pm$ 3E-07 ** | 9.99E-01 $\pm$ 4E-05
20 | 1.00E+00 $\pm$ 9E-07 ** | 9.99E-01 $\pm$ 6E-05
30 | 1.00E+00 $\pm$ 2E-06 ** | 9.99E-01 $\pm$ 8E-05
TABLE III: Comparison of the iNFFT, and iFFT reconstruction errors based on correlation between the trigonometric polynomial trained on all available data, versus the randomly sampled data at 10, 20, and 30% external sampling. * and ** denote a permutational, one-sided t-test indicating $p<0.05$ and $p<0.01$, respectively ($N=20,\mathcal{\pi}$ (permutations) $=10,000$) | $r^{iNFFT}_{|blck}$ | $r^{iFFT}_{|blck}$
---|---|---
10 | 9.91E-01 $\pm$ 6E-03 ** | 9.92E-01 $\pm$ 4E-02
20 | 9.61E-01 $\pm$ 3E-02 ** | 8.72E-01 $\pm$ 5E-02
30 | 8.83E0-1 $\pm$ 7E-02 ** | 7.97E-01 $\pm$ 4E-02
TABLE IV: Comparison of the iNFFT, and iFFT reconstruction errors based on correlation between the trigonometric polynomial trained on all available data, versus the contiguous block sampled data at 10, 20, and 30% external sampling. * and ** denote a permutational, one-sided t-test indicating $p<0.05$ and $p<0.01$, respectively ($N=20,\pi=10,000$) | $(1-Err)^{iNFFT}_{|rand}$ | $(1-Err)^{iFFT}_{|rand}$
---|---|---
10 | 9.86E-01 $\pm$ 3E-04 | 9.86E-01 $\pm$ 3E-02
20 | 9.73E-01 $\pm$ 5E-04 | 9.73E-01 $\pm$ 6E-02
30 | 9.59E-01 $\pm$ 6E-04 | 9.59E-01 $\pm$ 7E-02
TABLE V: Comparison of the iNFFT, and iFFT reconstruction errors based on external prediction for randomly selected data at 10, 20, and 30% external sampling. * and ** denote a permutational, one-sided t-test indicating $p<0.05$ and $p<0.01$, respectively ($N=20,\pi=10,000$) | $(1-Err)^{iNFFT}_{|blck}$ | $(1-Err)^{iFFT}_{|blck}$
---|---|---
10 | 9.91E-01 $\pm$ 6E-03 * | 9.28E-01 $\pm$ 5E-02
20 | 9.61E-01 $\pm$ 3E-02 * | 8.67E-01 $\pm$ 7E-02
30 | 8.83E-01 $\pm$ 7E-02 | 8.12E-01 $\pm$ 4E-02
TABLE VI: Comparison of the iNFFT, and iFFT reconstruction errors based on
external prediction for contiguous blocks of selected data at 10, 20, and 30%
external sampling. * and ** denote a permutational, one-sided t-test
indicating $p<0.05$ and $p<0.01$, respectively ($N=20,\pi=10,000$)

The relative error for the randomly sampled data is on the same order of
magnitude as reported in the original study for the iterative solution to the
optimal interpolation problem. However, it is clear that a more critical
metric for determining the accuracy of the interpolation is by removing
contiguous blocks of data, and comparing the reconstruction against the true
values for those data. Despite this, and considering that the interpolation
was “smoothed” using a rather aggressive low-pass filter, the results are
reasonable and corroborate that the calculated spectrum of the data is
relatively consistent with what data is observed, which is consistent with
what we would expect a true inverse transform to return.

The results also suggest that the interpolative algorithm performs
significantly better relative to the equivalent truncated inverse Fast Fourier
Transform (iFFT) using an equivalent number of coefficients and weight
function in virtually all cases. The differences are most differentiated
relative to correlation with the “true” underlying trigonometric polynomial,
as shown in Tables III and Tables IV; only in the case of randomly sampled
data is the difference inconsequential (Table V). Despite poorer evidence of
significance, on average the iNFFT algorithm in the most challenging use cases
for up to 30% missing data.

##  5 Conclusion

A convenient mathematical expression to solve for the spectrum of irregularly
sampled time-domain data has been demonstrated in a way that satisfies the
conditions for an inverse non-uniform Fourier transform. This methodology has
been derived from cost function proposed in the original study, and via its
idealised quadratic form strongly suggests that its minimum is unique and
convex. Using a randomly sample time series from data with few missing values,
the observed relative error is on the same order of magnitude as the most
current solution to this problem [10]. Upon cross-validation with contiguous
blocks, we report a relative error that is much higher, and suggest this is a
more critical, but informative lens through which to view interpolation
problems.

The proposed formula will find use in a variety of signal processing
applications, where the use of an FFT is impractical given its assumption of a
consistent time-domain density. The software package used in this work is
available online: https://github.com/CoDaSLab/intrp_infft_1d

##  Appendix A Proof of hermitian self-adjoint product identity for
equidistant time-domain measurements

Consider the case of an NFFT matrix and its hermitian, where:

$\displaystyle A=\left(e^{-2\pi\mathbf{i}k_{f}\frac{j}{M}}\right)_{j\in
M,k_{f}\in N}$ (42) $\displaystyle
A^{H}=\left(e^{2\pi\mathbf{i}k_{a}\frac{j}{M}}\right)_{j\in M,k_{a}\in N}$
(43)

and $j=-\frac{M}{2},...,\frac{M}{2}-1$, $k_{f}$ (representing the frequency-
wise columns of the forward transform) $=k_{a}$ (representing the same columns
for the adjoint transform) $=-\frac{N}{2},...,\frac{N}{2}-1$. Multiplying the
$N\times M$ matrix $A$ by its adjoint yields the $N\times N$ matrix:

$AA^{H}=\sum_{j=-M/2}^{M/2-1}e^{-2\pi\mathbf{i}jk_{f}}e^{2\pi\mathbf{i}k_{a}j}$
(44)

Equation 44 simplifies to:

$AA^{H}=\sum_{j=-M/2}^{M/2-1}e^{2\pi\mathbf{i}j(k_{a}-k_{f})}.$ (45)

It is clear from Equation 45 that $e^{2\pi\mathbf{i}j(0)}=1$ where
$k_{f}=k_{a}$. This indicates that the product of $AA^{H}$ along the diagonal
is $M$ for the summation of $j\in M$ in 45. To prove that the off-diagonal
matrices are zero, consider $\theta=2\pi(k_{a}-k_{f})$ as a complex number,
$z=\cos{\theta}+\mathbf{i}\sin{\theta}$ following Euler’s identity:

$\sum_{j=-M/2}^{M/2-1}e^{\mathbf{i}\theta j}=z$ (46)

Multiplying both sides of Equation 46 by $e^{\mathbf{i}\theta}$ yields the
following expression:

$ze^{\mathbf{i}\theta}=\sum_{j=-M/2}^{M/2-1}e^{\mathbf{i}\theta(j+1)}\\\ $
(47)

since $\sum_{j=0}^{M-1}e^{\mathbf{i}\theta(j+1)}$ is a $j^{th}$ integer
multiple of $\theta$,
$\sum_{j=0}^{M-1}e^{\mathbf{i}\theta(j+1)}=\sum_{j=0}^{M-1}e^{\mathbf{i}\theta
j}$. It follows that:

$ze^{\mathbf{i}\theta j}=z.$ (48)

Evaluating Equation 48, either $\theta=0$ in the case of Equation 45, or for
instances in which $k_{f}\neq k_{a}$, $z=0$. This proves that the off-diagonal
entries of $A^{H}A$ are zero, proving the identity for instances where the
$j^{th}$ time domain observations are equidistant:

$AA^{H}=MI_{N}$ (49)

##  Appendix B $AA^{H}\neq MI_{N}$ in the non-equidistant case

Consider again the equispaced example, where $j\in[-\frac{M}{2},\frac{M}{2})$.
We seek to prove that for some complex number, $z$ evaluated as the summation
of $\sum_{j=-M/2}^{M/2-1}e^{\mathbf{i}\theta j}=z$, $z=1$ when $\theta=0$ and
$z=1$ when $\theta\neq 0$. We prove the same identity, but in this case using
an evaluation of the finite geometric series, where:

$\sum_{i=n}^{n-1}r^{k}=\frac{1-r^{n}}{1-r}.$ (50)

since it is obvious that when $\theta=0$, the only possible evaluation for $z$
is $1$, consider only when $\theta\neq 0$, which we can write as:

$z=e^{-\mathbf{i}\theta\frac{M}{2}}+\sum_{j=-M/2+1}^{M/2-1}e^{\mathbf{i}\theta
j}$ (51)

to evaluate for $z$ for some constant integer $\theta$, consider the finite
geometric series:

$\sum_{j=-M/2}^{M/2-1}e^{\mathbf{i}\theta
j}=\frac{1-e^{\mathbf{i}\theta(M-1)}}{1-e}=0.$ (52)

Because $\sum_{-\frac{M}{2}+1}^{\frac{M}{2}-1}$ does not evaluate to a
convenient integer, evaluating it terms of known quantities:

$\sum_{j=-M/2}^{M/2-1}e^{\mathbf{i}\theta
j}=e^{-\mathbf{i}\theta\frac{M}{2}}+\sum_{j=-M/2+1}^{M/2-1}e^{\mathbf{i}\theta
j}$ (53)

since $e^{-\mathbf{i}\theta\frac{M}{2}}=-1$ and the LHS of Equation 53 was
evaluated in Equation 52, solving for
$\sum_{j=-M/2+1}^{M/2-1}e^{\mathbf{i}\theta j}$ yields $1$. Therefore, for
some complex number $z$, when $\theta\neq 0$, $z=0$.

Using this alternative formula, it is easy to prove that in the non-
equidistant case where $x_{j}\in[-\frac{M}{2},\frac{M}{2}-1)$ the geometric
series is not guaranteed, or even likely to sum to $M-1$. And so from Equation
52:

$\forall(\theta\neq 0),\exists\left(z:z=\sum_{j\in x_{j}}e^{\mathbf{i}\theta
j}\neq 0\right)$ (54)

## Acknowledgments

The proof for the equidistant case was adapted from a proof by David J. Fleet
and Allan D. Jepson from the University of Toronto [18].

## References

  * [1] C. Haley, “Missing-data nonparametric coherency estimation,” _IEEE Signal Processing Letters_ , vol. 28, pp. 1704–1708, 2021.
  * [2] S. Dilmaghani, I. C. Henry, P. Soonthornnonda, E. R. Christensen, and R. C. Henry, “Harmonic analysis of environmental time series with missing data or irregular sample spacing,” _Environmental science & technology_, vol. 41, no. 20, pp. 7030–7038, 2007.
  * [3] J. W. Cooley and J. W. Tukey, “An algorithm for the machine calculation of complex fourier series,” _Mathematics of computation_ , vol. 19, no. 90, pp. 297–301, 1965.
  * [4] F. Rasheed, P. Peng, R. Alhajj, and J. Rokne, “Fourier transform based spatial outlier mining,” in _Intelligent Data Engineering and Automated Learning-IDEAL 2009: 10th International Conference, Burgos, Spain, September 23-26, 2009. Proceedings 10_. Springer, 2009, pp. 317–324.
  * [5] C. Van Loan, _Computational frameworks for the fast Fourier transform_. SIAM, 1992.
  * [6] H. Carslaw, “A historical note on gibbs’ phenomenon in fourier’s series and integrals,” _Bulletin of the American Mathematical Society_ , 1925.
  * [7] J. Markel, “Fft pruning,” _IEEE transactions on Audio and Electroacoustics_ , vol. 19, no. 4, pp. 305–311, 1971.
  * [8] J. C. Bowman and Z. Ghoggali, “The partial fast fourier transform,” _Journal of Scientific Computing_ , vol. 76, pp. 1578–1593, 2018.
  * [9] M. Kircheis and D. Potts, “Direct inversion of the nonequispaced fast fourier transform,” _Linear Algebra and its Applications_ , vol. 575, pp. 106–140, 2019.
  * [10] S. Kunis and D. Potts, “Stability results for scattered data interpolation by trigonometric polynomials,” _SIAM Journal on Scientific Computing_ , vol. 29, no. 4, pp. 1403–1419, 2007.
  * [11] K. B. Petersen, M. S. Pedersen _et al._ , “The matrix cookbook,” _Technical University of Denmark_ , vol. 7, no. 15, p. 510, 2008.
  * [12] C. M. Bishop, _Neural networks for pattern recognition_. Oxford university press, 1995.
  * [13] H. Wickham, M. Averick, J. Bryan, W. Chang, L. D. McGowan, R. François, G. Grolemund, A. Hayes, L. Henry, J. Hester, M. Kuhn, T. L. Pedersen, E. Miller, S. M. Bache, K. Müller, J. Ooms, D. Robinson, D. P. Seidel, V. Spinu, K. Takahashi, D. Vaughan, C. Wilke, K. Woo, and H. Yutani, “Welcome to the tidyverse,” _Journal of Open Source Software_ , vol. 4, no. 43, p. 1686, 2019.
  * [14] R Core Team, _R: A Language and Environment for Statistical Computing_ , R Foundation for Statistical Computing, Vienna, Austria, 2023. [Online]. Available: https://www.R-project.org/
  * [15] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, “Array programming with NumPy,” _Nature_ , vol. 585, no. 7825, pp. 357–362, Sep. 2020. [Online]. Available: https://doi.org/10.1038/s41586-020-2649-2
  * [16] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,” _Nature Methods_ , vol. 17, pp. 261–272, 2020.
  * [17] J. D. Hunter, “Matplotlib: A 2d graphics environment,” _Computing in Science & Engineering_, vol. 9, no. 3, pp. 90–95, 2007.
  * [18] “Notes on fourier analysis,” http://www.cs.toronto.edu/~jepson/csc320/notes/fourier.pdf, accessed: 2023-09-01.


# Segmentation method for cerebral blood vessels from MRA using hysteresis

Georgia Kenyon1,3, Stephan Lau1,2, Michael A. Chappell3 and Mark Jenkinson1,2
1Australian Institute for Machine Learning, School of Computer and
Mathematical Sciences, The University of Adelaide, Australia,
<EMAIL_ADDRESS>Australian Health and Medical Research
Institute, Adelaide, Australia3Sir Peter Mansfield Imaging Centre, School of
Medicine, University of Nottingham, UK

###### Abstract

Segmentation of cerebral blood vessels from Magnetic Resonance Imaging (MRI)
is an open problem that could be solved with deep learning (DL). However,
annotated data for training is often scarce. Due to the absence of open-source
tools, we aim to develop a classical segmentation method that generates vessel
ground truth from Magnetic Resonance Angiography for DL training of
segmentation across a variety of modalities. The method combines size-specific
Hessian filters, hysteresis thresholding and connected component correction.
The optimal choice of processing steps was evaluated with a blinded scoring by
a clinician using 24 3D images. The results show that all method steps are
necessary to produce the highest (14.2/15) vessel segmentation quality score.
Omitting the connected component correction caused the largest quality loss.
The method, which is available on GitHub, can be used to train DL models for
vessel segmentation.

Clinical relevance— We provide and validate an open-access method to
efficiently segment vessels from brain MRA images for improved diagnostics and
intervention planning.

Keywords— Magnetic Resonance Angiography, segmentation, brain, blood vessel,
Hessian filter, hysteresis thresholding, connected component correction

##  I INTRODUCTION

Cerebrovascular diseases such as aneurysms and stenosis affect blood flow and
vessels in the brain. This may lead to a lack of sufficient blood flow, and as
a consequence cause stroke. Optimal treatment methods require imaging
techniques such as 3D time-of-flight (TOF) magnetic resonance angiography
(MRA) prior to intervention. However, visualization of blood vessels from raw
data - without appropriate pre-processing can be difficult, due to noise,
uneven contrast and limited signal from smaller vessels. This requires high-
level, automated pre-processing and segmentation methods to be developed.

Segmentation can be a challenging task, particularly for images of lower
quality. Deep learning provides a solution to this, as models can learn more
abstract patterns, such as vessels presented as small voxel intensity changes,
from data directly. However, DL requires large amounts of labelled data for
training, which is often infeasible. Manual annotation of vessels in 3D images
is a very tedious, time-consuming job that requires a high degree of
expertise- with labeling taking 60–80 min per MRA image [1]. Due to this,
developing (semi)-automatic approaches to provide vascular ground truth to
supervised approaches has become an important task.

Some existing works utilize DL for the segmentation of MRA, applying sparse
label techniques [2], synthetic or limited training sets [1] or unsupervised
learning [3], but these have not been made openly available. Existing vessel
data sets are limited or unsuitable for semantic segmentation tasks. For
example, Brava [4], an open-source vessel reconstruction data set, contains
reconstructions, but without the paired MR images necessary for DL training.
The freely available ’IXI’ data set [5], offers a large subset of MRA, along
with additional paired MR acquisitions, however these are unprocessed and do
not contain voxel-wise segmentations.

Multiple works have been published that use a combination of traditional
segmentation methods on MRA. They apply pre-processing filters that highlight
’vesselness’- tubular and high contrast features, use thresholding techniques
and apply voxel correction [6, 7]. However, there are no open-source methods
easily available for reuse. In 2016, Lu et al. [7] utilized the Hessian
filter, to enhance the vessels in MRA, prior to segmentation. This filter uses
edge-based enhancement for tubular structures, by taking advantage of the
strong gradients at the edges of vessels in comparison to the background. In
2020, Avadiappan et al. [8] employed an alternative method to implement the
concept of ’vesselness’, using an adaptive Frangi filter to preserve the radii
of vessels. The adaptive filter was implemented to quantify vessel radii, as
it was summarized that existing methods underestimated the segmentation for
vessels of larger radii and overestimated thin vessels. Multiple methods apply
thresholding techniques to extract vessels, applying Gaussian Mixture models
over the intensity distributions to separate anatomic structures [7, 9, 6].
However, these single threshold methods can be limited when extracting fine
detail. Thresholding with hysteresis, developed for the Canny Edge detector
[10], applies a high and low threshold, to preserve intensity information.
Chang et al. [11] applied a modified version of the Canny Edge Detector to
detect small retinal blood vessels, generating a dynamic hysteresis threshold
value based on local neighborhood of a pixel, improving the detection of
smaller edges. In 2021, Ooi et al. [12] applied the edge detector to segment
retinal images, after hand-selecting the vessels, to highlight finer details
and draw out vessel edges. Voxel correction methods, such as erosion, are also
applied to segmentation outputs to remove noise after thresholding [7].
However, small vessel misclassification as noise can limit these methods.

This paper aims to develop a method, shown in Figure 1, to segment cerebral
blood vessels from MRA images. The method provides a solution to data
labelling required to train deep learning models for segmentation tasks using
non-MRA modalities (e.g. T2-weighted images) by utilizing paired data sets
such as IXI. The method includes a novel combination of processing elements
that exploit the benefits of the Hessian Matrix to enhance vessels regardless
of vessel radii, and segment using hysteresis thresholds. The method also
utilizes the benefits of connected component correction to remove noise and
artifacts, whilst maintaining small vessel integrity.

##  II MATERIALS & METHODS

Figure 1: Proposed method structure (left) with axial view of cerebral blood
vessels (red) segmented with the proposed method from an example MRA (right)
shown with one MRA slice (grey) for anatomical reference.

###  II-A Imaging Data & Pre-processing

3D TOF-MRA was extracted from the IXI database [5]. Only 3T and 1.5T images
acquired in Hammersmith Hospital (IXI-HH) or Guy’s Hospital (IXI-Guys), with a
resolution of 0.47 x 0.47 x 0.8 $mm^{3}$, were used. The Brain Extraction Tool
(BET) in FSL [13] was used to skull-strip the TOF-MRA, to exclude non-vascular
skull voxels of similar intensity. Bias field correction was then applied,
using FSL’s FAST.

###  II-B Vessel-Enhancement Filter

A 3D line-enhancement filter based on the eigenvalues of the Hessian Matrix
was used to enhance the tubular structure of vessels in the image [14]. An
anisotropic Gaussian filter was applied during Hessian, to reduce noise, and
was varied using the algorithm’s standard deviation parameter $\sigma$. The
$\sigma$ values were chosen based on the image resolution, to optimize the
extraction of the vessels specific to their radii [15], since low $\sigma$
values would highlight small, thin structures, whilst larger values would
highlight the larger structures. Due to the variation of the vessel structures
within the MRA, two values for sigma were applied, producing two image outputs
per MRA input. The value $\sigma_{l}=0.47$ was chosen to conserve fine
vascular detail of the small vessels, whilst $\sigma_{u}=0.94$ was chosen to
promote the tubular vessels with larger diameters, reaching two voxels, such
as in the Circle of Willis [14].

###  II-C Hysteresis Thresholding

Both vessel enhancement outputs were thresholded at a fraction of the 99.9th
percentile intensity of the image, respectively. The intensity value at the
99.9th percentile was used instead of the maximum to prevent high intensity
noise and outliers impacting this value. The decision to use fractions of the
intensity was made so that the threshold applied uniformly across the data
set. Importance was placed on improving the quality of the segmentation for
larger vessels from the $\sigma_{u}$ filtered images, whilst the extraction of
finer vessels was the focus for thresholding the $\sigma_{l}$ filtered images.
The absolute value of the threshold is obtained from the relative threshold
${th_{p}}$ by applying it to the image’s 99.9th percentile intensity value.

Hysteresis thresholding was applied to improve the detection of smaller edges
from the vessels. The method took two threshold values: a low-threshold value
(LTV) and a high-threshold value (HTV). A voxel was included if its intensity
was above the HTV, or between the LTV and HTV and connected to a voxel above
the HTV.

To preserve the fine detail of the smaller vessels for the $\sigma_{l}$
outputs, $\overline{th_{p}}$ was applied as the HTV (e.g. HTV = 67), whilst
the LTV was chosen at 10% below $\overline{th_{p}}$ (e.g. LTV = 57). The lower
range was chosen to promote inclusion of extra voxels of lower intensities if
they were connected to an established edge. The hysteresis threshold range for
the data produced using the $\sigma_{u}$ filter, utilized the
$\overline{th_{p}}$ as the low-threshold value for hysteresis (e.g. LTV = 39),
and a HTV value ten percent above this (e.g. HTV = 49). This was chosen to
remove unwanted noise and extract the vessels with larger intensity
variations. All vessel data was binarized after hysteresis, with background
and brain voxels set to zero, and vessels set to one. The two binary outputs
were then merged (with a union operation) to form one image.

###  II-D Connected Component Selection

A voxel correction algorithm was implemented, to reduce the noise and non-
vessel structures left behind after hysteresis thresholding. This was achieved
using 3D connected component correction, which exploits the continuous
structure of the blood vessels to remove small, disconnected voxel clusters of
high intensity left behind by thresholding, suspected to be noise. Any
clusters of size less than $10mm^{3}$ were removed.

###  II-E Clinical Expert Scoring

To determine an optimal set of hysteresis thresholds, outputs over a range of
values were presented to a clinical expert, who then chose the option with the
highest anatomical accuracy.

An ablation experiment was also performed, removing the four steps of the
method one at a time, and the 24 images from the IXI data set were reprocessed
by each alternative method. 3D segmentation outputs were displayed to a
trained clinician using FSLeyes [13] and overlaid over the original MRA to
enhance detail visibility. The clinician was instructed to score the output
images of the alternative methods using metrics determined from literature
[16], as shown in Table I.

TABLE I: Quality Scores for Ablation Experiment Metrics |  Description: score out of 3
---|---
Connectivity |  Poor level of connectivity (1), mostly connected (2), fully connected (3).
Segmentation |  Segmented 50% or less (1), 60-80% (2), more than 80% of the vessels (3).
Fine Detail |  Limited/no small vessel segmentation (1), small gaps in small vessels (2), lots of fine/small vessel detail (3).
False Positives |  Lots of (1), limited (2), minimal false positive vessel segmentations (3).
Noise |  Lots of noise (1), some/acceptable levels of noise (2), minimal Noise (3).

##  III RESULTS

###  III-A Hysteresis Threshold Range

The clinician was tasked to select the optimal threshold range for two sets of
Hessian filtered images ($\sigma_{l}$ and $\sigma_{u}$) out of five options.
Figure 2 is an example of the thresholding study used, with the output
segmentation for five threshold ranges provided to the clinician to score
alongside the original MRA.

Figure 2: Example of hysteresis thresholds for validation of $\sigma_{u}$
filtered images: a) Original MRA, b) High threshold (H)=69, Low threshold
(L)=59, c) H=59, L=49, d) H=49, L=39, e) H=39, L=29, f) H=29, L=19

The results from the scoring of the hysteresis thresholds are shown in Table
II and Table III. The spread of clinician preference was one threshold step
wider for $\sigma_{u}$ in Table III, than for $\sigma_{l}$ in Table II.

TABLE II: Percentage of Cases Selected by the Clinician as the Best Hysteresis Threshold Option for ($\sigma_{l}$) Lower Threshold % | Upper Threshold % | Clinician Preference
---|---|---
77 | 87 | 0%
67 | 77 | 17%
57 | 67 | 83%
47 | 57 | 0%
37 | 47 | 0%
TABLE III: Percentage of Cases Selected by the Clinician as the Best Hysteresis Yhreshold Option for ($\sigma_{u}$) Lower Threshold % | Upper Threshold % | Clinician Preference
---|---|---
59 | 69 | 0%
49 | 59 | 9%
39 | 49 | 58%
29 | 39 | 33%
19 | 29 | 0%

###  III-B Ablation Study

The results of the ablation study are shown in Table IV, with the full method
achieving the highest quality score. The output without connected component
selection scored lowest but had a lower score variance than the remaining
alternatives.

TABLE IV: Average Clinician Quality Score for Alternative Methods Alternative Method | Quality Score, out of 15 | Score s.d
---|---|---
Full Method | 14.2 | 0.9
Removal of $\sigma_{u}$ | 12.4 | 1.1
Removal of $\sigma_{l}$ | 12.8 | 0.9
No Hysteresis Thresholding | 14.1 | 1.1
No Connected Component | 11.0 | 0.2

##  IV DISCUSSION

A hysteresis range of 57% to 67% of the 99.9th percentile was determined to be
the optimal choice for the $\sigma_{l}$ Hessian filtered image, as the
clinician concluded that it was the optimal threshold range for the majority
(83%) of the images displayed. A range of 39% to 49% of the 99.9th percentile
was chosen for the $\sigma_{u}$ filtered images, as it was considered optimal
for the majority (58%) of the images.

The ablation study confirmed that all steps of the method, in the specific
sequence established, were required to optimize the segmentation output. The
method produced outputs such as the slices presented in Figure 3 \- and the
full method achieved the highest quality score, as shown in Table IV. Figure 3
c & d highlighted the necessity for two Hessian filters to be applied to MRA
prior to thresholding. The different $\sigma$ values for larger and smaller
vessels optimized the thresholding in different ways. The green rings
highlighted complete loss of small vessel detail when a filter of $\sigma_{l}$
was omitted, and under-prediction of larger vessels when the $\sigma_{u}$
filter output was removed. The study corroborated the literature [8] that
noted the importance of using different sigma values for Hessian filters for
vessels of different radii. Removal of hysteresis thresholding was shown to
have the smallest effect on vessel segmentation, with a score of 14.1 compared
to 14.2 for the full method. However, the maximum score was only achieved with
the full method, and examples presented to the clinician highlighted the loss
of vessel edges when the hysteresis thresholding step was omitted, confirming
a similar result in retinal imaging [11].

Morphological operations (e.g. erosion) were not used due to the thin
structure of the vessels being eroded too aggressively. The connected
component correction solved this issue by removing clusters of noise without
affecting the structural integrity of the remaining vessels. The ablation
study showed that component correction was the method’s most important step,
with the most severe quality loss compared to removing any of the other steps
(score went from 14.2 to 11.0). The green ring over the MRA segmentation f) in
Figure 3 highlights the advantage of the step, as lots of noise remains on
removal of it. This confirms the importance of removing scattered errors, as
noted by Lu et al. [7].

Future work will use more participants, from different centers and data sets
and fine-tune the parameters. However, although the quality score of 14.2 out
of 15 leaves some room for improvement, this may only be achievable with deep
learning. Our method provides sufficient ground truth for the training of a
deep learning model on both MRA and non-MRA images, and future DL models can
surpass this performance by using these segmentations as weak or noisy labels,
and generalize these across even larger data sets.

Figure 3: Example slice of alternative method outputs: a) Original MRA, b)
Full method c) No $\sigma_{l}$ filter, d) No $\sigma_{u}$ filter, e) No
hysteresis thresholding, f) No component correction

##  V CONCLUSION

We presented and evaluated a novel method for the automated segmentation of
cerebral blood vessels from MRA images. The method and IXI segmentations are
available on GitHub, and can be used as weak/noisy labels for training deep
learning models to segment MRA and non-MRA images.

## ACKNOWLEDGMENT

We would like to thank Minyan Zeng from the Medical School at The University
of Adelaide for the blinded labelling.

## References

  * [1] M. Livne, J. Rieger, O. U. Aydin, A. A. Taha, E. M. Akay, T. Kossen, J. Sobesky, J. D. Kelleher, K. Hildebrand, D. Frey _et al._ , “A u-net deep learning framework for high performance vessel segmentation in patients with cerebrovascular disease,” _Front. Neurosci._ , vol. 13, p. 97, 2019\.
  * [2] B. Zhang, S. Liu, S. Zhou, J. Yang, C. Wang, N. Li, Z. Wu, and J. Xia, “Cerebrovascular segmentation from tof-mra using model-and data-driven method via sparse labels,” _Neurocomputing_ , vol. 380, pp. 162–179, 2020\.
  * [3] S. Fan, Y. Bian, H. Chen, Y. Kang, Q. Yang, and T. Tan, “Unsupervised cerebrovascular segmentation of tof-mra images based on deep neural network and hidden markov random field model,” _Front. Neuroinform._ , vol. 13, p. 77, 2020.
  * [4] S. N. Wright, P. Kochunov, F. Mut, M. Bergamino, K. M. Brown, J. C. Mazziotta, A. W. Toga, J. R. Cebral, and G. A. Ascoli, “Digital reconstruction and morphometric analysis of human brain arterial vasculature from magnetic resonance angiography,” _Neuroimage_ , vol. 82, pp. 170–181, 2013.
  * [5] “IXI dataset,” available online, accessed 11 November 2022. [Online]. Available: https://brain-development.org/ixi-dataset/
  * [6] N. Li, S. Zhou, Z. Wu, B. Zhang, and G. Zhao, “Statistical modeling and knowledge-based segmentation of cerebral artery based on TOF-MRA and MR-T1,” _Comput. Methods Programs Biomed._ , vol. 186, p. 105110, 2020\.
  * [7] P. Lu, J. Xia, Z. Li, J. Xiong, J. Yang, S. Zhou, L. Wang, M. Chen, and C. Wang, “A vessel segmentation method for multi-modality angiographic images based on multi-scale filtering and statistical models,” _Biomed. Eng. Online_ , vol. 15, no. 1, pp. 1–18, 2016.
  * [8] S. Avadiappan, S. Payabvash, M. A. Morrison, A. Jakary, C. P. Hess, and J. M. Lupo, “A fully automated method for segmenting arteries and quantifying vessel radii on magnetic resonance angiography images of varying projection thickness,” _Front. Neurosci._ , vol. 14, p. 537, 2020.
  * [9] L. Wen, X. Wang, Z. Wu, M. Zhou, and J. S. Jin, “A novel statistical cerebrovascular segmentation algorithm with particle swarm optimization,” _Neurocomput._ , vol. 148, pp. 569–577, 2015.
  * [10] J. F. Canny, “Finding edges and lines in images,” Massachusetts Institute of Technology Cambridge Artificial Intelligence Lab, Tech. Rep., 1983.
  * [11] S. H. Chang, L. Gong, M. Li, X. Hu, and J. Yan, “Small retinal vessel extraction using modified canny edge detection,” in _2008 Int. Conf. on Audio, Language and Image Processing_. IEEE, 2008, pp. 1255–1259.
  * [12] A. Z. H. Ooi, Z. Embong, A. I. Abd Hamid, R. Zainon, S. L. Wang, T. F. Ng, R. A. Hamzah, S. S. Teoh, and H. Ibrahim, “Interactive blood vessel segmentation from retinal fundus image based on canny edge detector,” _Sensors_ , vol. 21, no. 19, p. 6380, 2021.
  * [13] M. Jenkinson, C. F. Beckmann, T. E. Behrens, M. W. Woolrich, and S. M. Smith, “FSL,” _NeuroImage_ , vol. 62, no. 2, pp. 782–790, 2012, 20 YEARS OF fMRI. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1053811911010603
  * [14] Y. Sato, S. Nakajima, H. Atsumi, T. Koller, G. Gerig, S. Yoshida, and R. Kikinis, “3D multi-scale line filter for segmentation and visualization of curvilinear structures in medical images,” in _CVRMed-MRCAS’97_. Springer, 1997, pp. 213–222.
  * [15] A. Deshpande, N. Jamilpour, B. Jiang, P. Michel, A. Eskandari, C. Kidwell, M. Wintermark, and K. Laksari, “Automatic segmentation, feature extraction and comparison of healthy and stroke cerebral vasculature,” _NeuroImage Clin._ , vol. 30, p. 102573, 2021.
  * [16] S. Moccia, E. De Momi, S. El Hadji, and L. S. Mattos, “Blood vessel segmentation algorithms — review of methods, datasets and evaluation metrics,” _Comput. Methods Programs Biomed._ , vol. 158, pp. 71–91, 2018\.


# Free fermionic webs of heterotic T–folds

Alon E. Faraggi<EMAIL_ADDRESS>Dept. of Mathematical Sciences,
University of Liverpool, Liverpool L69 7ZL, UK
and CERN, Theoretical Physics Department, CH–-1211 Geneva 23, Switzerland
Stefan Groot Nibbelink<EMAIL_ADDRESS>School of Engineering and
Applied Sciences and Research Centre Innovations in Care, Rotterdam University
of Applied Sciences, Postbus 25035, 3001 HA Rotterdam, the Netherlands
Benjamin Percival<EMAIL_ADDRESS>Dept. of Mathematical Sciences,
University of Liverpool, Liverpool L69 7ZL, UK

###### Abstract

Moduli stabilisation is key to obtaining phenomenologically viable string
models. Non–geometric compactifications, like T–duality orbifolds (T–folds),
are capable of freezing many moduli. However, in this Letter we emphasise that
T–folds, admitting free fermionic descriptions, can be associated with a large
number of different T–folds with varying number of moduli, since the fermion
pairings for bosonisation are far from unique. Consequently, in one
description a fermionic construction might appear to be asymmetric, and hence
non–geometric, while in another it admits a symmetric orbifold description. We
introduce the notion of intrinsically asymmetric T–folds for fermionic
constructions that do not admit any symmetric orbifold description after
bosonisation. Finally, we argue that fermion symmetries induce mappings in the
bosonised description that extend the T–duality group.

††preprint: LTH-1343††preprint: CERN-TH-2023-115

##  I  Introduction

String theory realises a unification of gravity, gauge interactions and their
charged matter via the properties of Conformal Field Theories (CFTs) residing
on its two dimensional (2D) worldsheet. Heterotic strings on toroidal
orbifolds [1, 2] led to some of the most realistic string–derived models to
date [3, 4, 5]. However, orbifolds and other geometrical backgrounds result in
free moduli (such as the metric, B–field or Wilson lines) on which detailed
physics, like gauge and Yukawa couplings, depend.

Strings on tori and their orbifolds admit exact quantisation. This was
instrumental in the discovery of T–dualities [6], like the famous
$R\rightarrow 1/R$ duality, which sets the effective minimum of the radius $R$
equal to the string scale. Investigations of string backgrounds had a profound
impact on mathematics as mirror symmetry showed, which was argued to be a form
of T–duality [7].

Modding out T–duality symmetries may lead to exotic non–geometric backgrounds
[8, 9], dubbed T–folds. Hence, the landscape of string vacua may be much
vaster than suggested by geometrical compactifications alone. Even though
non–geometric constructions have been studied far less than their geometric
counterparts, they may be vital for phenomenological string explorations, as
they are capable of freezing many moduli.

Such T–folds may have different actions on their left– and right–moving
bosonic coordinate fields, and are thus referred to as asymmetric orbifolds
[10, 11]. If only order two symmetries are modded out, an alternative
fermionic description may be obtained by bosonisation, a CFT equivalence of
chiral bosons and fermions in 2D [12]. This led to a detailed dictionary
between these two formulations explicated for symmetric
$\mathbbm{Z}_{2}\\!\times\\!\mathbbm{Z}_{2}$ orbifolds [13]. Asymmetric
boundary conditions in the fermionic formalism have profound phenomenological
consequences, such as doublet–triplet splitting mechanism [14, 15], Yukawa
coupling selection [16], and moduli fixing [17].

Although a similar dictionary for asymmetric orbifolds is not this letter’s
aim, heterotic bosonisation ambiguities suggest identifications of seemingly
unrelated T–folds. This sheds new light on non–geometric moduli stabilisation.
Fermionic symmetries parameterising these ambiguities, suggest an extension of
the T–duality group.

##  II  Order two bosonic T–fold models

The bosonic formulation of the heterotic string [18] describes $d$–dimensional
Minkowski space by coordinate fields $x=(x^{\mu=2\ldots d-1})$ in the
light–cone gauge. The internal coordinate fields
$X=(X_{\text{R}}|X_{\text{L}})$, with right– and left–chiral parts,
$X_{\text{R}}=(X_{\text{R}}^{i=1\ldots D})$ and
$X_{\text{L}}=(X_{\text{L}}^{a=1\ldots D+16})$, $D=10-d$, are subject to torus
periodicities

$\displaystyle X\sim X+2\pi\,N\,,\qquad N\in\mathbbm{Z}^{2D+16}\,.$ (1)

The worldsheet supersymmetry current,

$\displaystyle T_{\text{F}}(z)=i\,\psi_{\mu}\partial
x^{\mu}+i\,\chi^{i}\partial X_{\text{R}}^{i}\,,$ (2)

involves the real holomorphic superpartners $\psi=(\psi^{\mu})$ and
$\chi=(\chi^{i})$ of $x$ and $X_{\text{R}}$, respectively. Here,
$(\bar{\partial})\partial$ denotes (anti–)holomorphic worldsheet derivative
and repeated indices are summed over.

An order two generator, defining the orbifold action

$\displaystyle X\sim e^{2\pi i\,v}\,X-2\pi\,V\,,$ (3)

with
$v=(v_{\text{R}}|v_{\text{L}}),\,V=(V_{\text{R}}|V_{\text{L}})\in\frac{1}{2}\,\mathbbm{Z}^{2D+16}$,
is called a shift, a twist, or a roto–translation, if $V\not\equiv 0\equiv v$,
$v\not\equiv 0\equiv V$, or $v,V\not\equiv 0$, respectively. ($\equiv$ means
equal up to integral vectors.)

An orbifold is called symmetric if there is a basis such that the left– and
right–twist parts are equivalent according to

$\displaystyle v_{\text{L}}\equiv(v_{\text{R}},0^{16})\,,$ (4)

for all its generators simultaneously [19, 20, 21], and asymmetric if no such
basis exists. The addition of $0^{16}$ is essential, as the vectors
$v_{\text{L}}$ and $v_{\text{R}}$ have unequal lengths.

##  III  Real free fermionic models

In the free fermionic formulation [22, 23, 24] the internal degrees of freedom
are described by real holomorphic fermions $f=(y,w)$ with $y=(y^{i})$ and
$w=(w^{i})$ and real anti–holomorphic fermions
${\bar{f}}=({\bar{f}}^{u=1\ldots 2D+32})$. Worldsheet supersymmetry is
realised non–linearly

$\displaystyle T_{\text{F}}(z)=i\psi_{\mu}\partial
x^{\mu}+i\chi^{i}y^{i}w^{i}\,.$ (5)

A fermionic model is defined by a set of basis vectors with entries $0$ or $1$
for real fermions. Each basis vector
${\boldsymbol{\beta}}=(\beta|{\bar{\beta}})$ defines boundary conditions

$\displaystyle f\sim-e^{\pi i\,\beta}\,f~{},\qquad{\bar{f}}\sim-e^{\pi
i\,{\bar{\beta}}}\,{\bar{f}}\,.$ (6)

##  IV  Bosonisations

###  IV.1 Holomorphic bosonisation

Assuming that the fermions $\chi$ are identical in the supercurrents (2) and
(5) and they generate the same worldsheet supersymmetry in the bosonic and
fermionic descriptions, bosonisation uniquely relates the currents

$\displaystyle J^{i}=\ :\\!(\lambda^{i})^{*}\lambda^{i}\\!:\ =\
:\\!y^{i}w^{i}\\!:\ \cong i\,\partial X_{\text{R}}^{i}\,,$ (7a) and complex
fermions $\displaystyle\lambda^{i}=\frac{1}{\sqrt{2}}(y^{i}+i\,w^{i})\cong\
:\\!e^{i\,X_{\text{R}}^{i}}\\!:$ (7b)

to normal ordered exponentials of chiral bosons. Here $\cong$ emphasises that
these expressions are not identities but rather that both sides have identical
operator product expansions in either formulation.

The bosonisation formulae relates the boundary conditions in both
descriptions. The torus periodicities (1) reflect the $2\pi$ ambiguities of
$X_{\text{R}}$ in the complex exponentials (7b). Comparing the orbifold
conditions (3) of the right–moving bosons $X_{\text{R}}$ with boundary
conditions (6) of the holomorphic fermions $y$ and $w$ in (7) leads to the
following identifications:

$\displaystyle
v_{\text{R}}=\mbox{$\frac{1}{2}$}\beta(w)-\mbox{$\frac{1}{2}$}\beta(y)\,,\qquad
V_{\text{R}}=\mbox{$\frac{1}{2}$}(1^{D})-\mbox{$\frac{1}{2}$}\beta(y)\,.$ (8)

###  IV.2 Anti–holomorphic bosonisation

Contrary to the holomorphic side, the pairing of the anti–holomorphic fermions
is arbitrary. Associating odd and even fermion labels to the real and
imaginary parts of complex fermions results in an anti–holomorphic
bosonisation procedure given by:

$\displaystyle\overline{J}^{a}=\
:\\!({\bar{\lambda}}^{a})^{*}{\bar{\lambda}}^{a}\\!:\ =\
:\\!\bar{f}^{2a-1}\bar{f}^{2a}\\!:\ \cong
i\,\bar{\partial}X_{\text{L}}^{a}~{},$ (9a) with
$\displaystyle{\bar{\lambda}}^{a}=\frac{1}{\sqrt{2}}({\bar{f}}^{2a-1}+i{\bar{f}}^{2a})\cong\
:\\!e^{iX_{\text{L}}^{a}}\\!:~{},$ (9b)

for $a=1,\ldots D+16$.

Then by similar arguments as above, the torus periodicities (1) for
$X_{\text{L}}$ follow. And splitting
${\bar{\beta}}=({\bar{\beta}}_{\text{o}},{\bar{\beta}}_{\text{e}})$ in two
$(D+16)$–dimensional vectors,
${\bar{\beta}}_{\text{o}}=({\bar{\beta}}^{1,3\ldots 2D+31})$ and
${\bar{\beta}}_{\text{e}}=({\bar{\beta}}^{2,4\ldots 2D+32})$ leads to the
identifications:

$\displaystyle
v_{\text{L}}=\mbox{$\frac{1}{2}$}{\bar{\beta}}_{\text{e}}-\mbox{$\frac{1}{2}$}{\bar{\beta}}_{\text{o}}\,,\qquad
V_{\text{L}}=\mbox{$\frac{1}{2}$}(1^{D+16})-\mbox{$\frac{1}{2}$}{\bar{\beta}}_{\text{o}}\,.$
(10)

Table 1:  Fermionic symmetry induced bosonic coordinate field transformations. (Only non–inert fields are given.)  Fermionic symmetry | Action on left–moving bosons
---|---
$(2a\\!{\,\mbox{-}}\\!1\,2b\\!{\,\mbox{-}}\\!1)(2a\,2b)$ | $X_{\text{L}}^{a}\leftrightarrow X_{\text{L}}^{b}$
$(2a)$ | $X_{\text{L}}^{a}\rightarrow-X_{\text{L}}^{a}$
$(2a\\!{\,\mbox{-}}\\!1)$ | $X_{\text{L}}^{a}\rightarrow\pi-X_{\text{L}}^{a}$
$(2a\\!{\,\mbox{-}}\\!1\,2a)$ | $X_{\text{L}}^{a}\rightarrow\mbox{$\frac{1}{2}$}\pi-X_{\text{L}}^{a}$
Table 2:  Fermionic symmetry induced bosonic boundary condition mappings. (Only non–inert entries of the vectors $v_{\text{L}}$ and $V_{\text{L}}$ modulo integral vectors are given.)  Fermionic symmetry | Action on twist and shift entries
---|---
$(2a\\!{\,\mbox{-}}\\!1\,2b\\!{\,\mbox{-}}\\!1)(2a\,2b)$ | $v_{\text{L}}^{a}\leftrightarrow v_{\text{L}}^{b}\,,~{}V_{\text{L}}^{a}\leftrightarrow V_{\text{L}}^{b}$
$(2a)$ | $v_{\text{L}}^{a}\rightarrow-v_{\text{L}}^{a}+2\,V_{\text{L}}^{a}\,,~{}V_{\text{L}}^{a}\rightarrow V_{\text{L}}^{a}$
$(2a\\!{\,\mbox{-}}\\!1\,2a)$ | $v_{\text{L}}^{a}\rightarrow-v_{\text{L}}^{a}\,,~{}V_{\text{L}}^{a}\rightarrow V_{\text{L}}^{a}-v_{\text{L}}^{a}$
(2a 2b) | $\begin{array}[]{l}v_{\text{L}}^{a}\rightarrow v^{b}_{\text{L}}+V^{a}_{\text{L}}-V^{b}_{\text{L}}\,,~{}V^{a}_{\text{L}}\rightarrow V^{a}_{\text{L}}\,,\\\\[-2.15277pt] v_{\text{L}}^{b}\rightarrow v^{a}_{\text{L}}+V^{b}_{\text{L}}-V^{a}_{\text{L}}\,,~{}V^{b}_{\text{L}}\rightarrow V^{b}_{\text{L}}\end{array}$
(2a - 1 2b -1) | $\begin{array}[]{l}v_{\text{L}}^{a}\rightarrow v^{a}_{\text{L}}+V^{a}_{\text{L}}-V^{b}_{\text{L}}\,,~{}V^{a}_{\text{L}}\rightarrow V^{b}_{\text{L}}\,,\\\\[-2.15277pt] v_{\text{L}}^{b}\rightarrow v^{b}_{\text{L}}+V^{b}_{\text{L}}-V^{a}_{\text{L}}\,,~{}V^{b}_{\text{L}}\rightarrow V^{a}_{\text{L}}\end{array}$
$(2a\\!{\,\mbox{-}}\\!1)$ | $v_{\text{L}}^{a}\rightarrow v_{\text{L}}^{a}-2\,V_{\text{L}}^{a}\,,~{}V_{\text{L}}^{a}\rightarrow-V_{\text{L}}^{a}$

##  V  Extension of the T–duality group

###  V.1 Fermionic inversions and permutations

On the anti–holomorphic side, the fermionic symmetries contain inversions and
permutations: $(u)$ denotes the fermion inversion
${\bar{f}}^{u}\rightarrow-{\bar{f}}^{u}$. The permutation $(u_{1}\ \cdots\
u_{p})$ acts as
${\bar{f}}^{u_{1}}\rightarrow{\bar{f}}^{u_{2}}\cdots\rightarrow{\bar{f}}^{u_{p}}\rightarrow{\bar{f}}^{u_{1}}$
leaving the remaining fermions inert. The permutation group contains elements
which consist of multiple factors like this provided their entries are all
distinct. It is generated by permutations of two elements $(uv)$. The induced
fermionic symmetry actions within the bosonic formulation can be identified
using the bosonisation (9).

###  V.2 Induced bosonic coordinate transformations

The fermionic symmetries, that leave these fermion bosonisation pairs intact,
realise mappings of the bosonic coordinate fields $X_{\text{L}}$ to
themselves. Their generators and their realisations on the bosonic coordinates
are listed in Table 1. The bosonic transformations above the middle line of
this table are part of the T–duality group, while those below involve
translations as well.

###  V.3 Induced mappings of bosonic boundary conditions

Other fermionic symmetries break up fermion bosonisation pairs and hence
correspond to mappings between different coordinate fields between which no
obvious coordinate transformation exists. However, all fermionic symmetries,
generated by inversions and permutations, map the boundary conditions of one
orbifold theory to another. The mappings induced by the generators of the
fermionic symmetries are collected in Table 2. The transformations induced by
permutations $(2a\,2b)$ and $(2a\\!{\,\mbox{-}}\\!1\,2b\\!{\,\mbox{-}}\\!1)$
combined (in whatever order) leads to the boundary condition mapping
associated with $(2a\\!{\,\mbox{-}}\\!1\,2b\\!{\,\mbox{-}}\\!1)(2a\,2b)$ as
the group property would suggest. Since some actions can be interpreted as
T–duality transformations, while others cannot, this hints at an extension of
the T–duality group.

The Table 2 mappings $(2a\\!{\,\mbox{-}}\\!1\,2a)$,
$(2a\\!{\,\mbox{-}}\\!1\,2b\\!{\,\mbox{-}}\\!1)$ and $(2a\,2b)$ are of special
significance: they mix the twist and shift vector entries. The action of
$(2a\\!{\,\mbox{-}}\\!1\,2a)$ recalls that the shift part of a
roto–translation in directions, where the twist acts non–trivially, can be
removed via the associated coordinate transformation in Table 1. The actions
$(2a\\!{\,\mbox{-}}\\!1\,2b\\!{\,\mbox{-}}\\!1)$ and $(2a\,2b)$ imply that a
pure shift boundary condition can be turned into a roto–translation. By
combining these mappings, a web of equivalent (mostly asymmetric) orbifold
theories emerge.

Since all these T–folds are just different bosonic representations of the same
fermionic theory, their physical properties are identical, even though they
may not look alike. For example, their modular invariance conditions may seem
to disagree, as the number of non–zero entries in the twist vectors under
mappings, like $(2a\,2b)$, change. However, since only the part of the shift
of (3), on which the twist acts trivially, takes part in the modular
invariance condition [21], their consistency conditions are numerically
identical.

###  V.4 A free fermionic T–fold web

Table 3:  Fermionic basis vectors ${\boldsymbol{\beta}}$ with the twists $v$ and shifts $V$ corresponding to $\mathbf{1}\\!\\!-\\!\\!{\boldsymbol{\beta}}$ obtained via (8) and (10).  Fermionic basis vectors | Twist and | shift vectors
---|---|---
${\boldsymbol{\beta}}$ | $v\\!\\!=\\!\\!(v_{\text{R}}|v_{\text{L}})$ | $V\\!\\!=\\!\\!(V_{\text{R}}|V_{\text{L}})$
$\mathbf{1}\\!\\!=\\!\\!\\{\psi^{1\ldots 4}\chi^{1\dots 4}y^{1\dots 4}w^{1\dots 4}|\bar{f}^{1\ldots 40}\\}$ | $(0^{4}|0^{20})$ | $\mbox{$\frac{1}{2}$}(1^{4}|1^{20})$
$\mathbf{S}\\!\\!=\\!\\!\\{\psi^{1\ldots 4}\chi^{1\dots 4}\\}$ | $(0^{4}|0^{20})$ | $(0^{4}|0^{20})$
${\boldsymbol{\xi}}\\!\\!=\\!\\!\\{\bar{f}^{9\dots 40}\\}$ | $(0^{4}|0^{20})$ | $\mbox{$\frac{1}{2}$}(0^{4}|0^{4}1^{16})$
$\mathbf{b}\\!\\!=\\!\\!\\{\chi^{1\ldots 4}w^{1,\ldots 4}|\bar{f}^{1\ldots 4}\bar{f}^{9\ldots 12}\\}$ | $\mbox{$\frac{1}{2}$}(1^{4}|0^{20})$ | $\mbox{$\frac{1}{2}$}(0^{4}|1^{2}0^{2}1^{2}0^{14})$

The basis vectors, ${\boldsymbol{\beta}}$, for a simple illustrative 6D
fermionic model are given in Table 3 together with associated twist and shift
vectors using the odd–even pairings (9). Within this bosonisation the model is
understood as an asymmetric orbifold. The interpretation may change by
applying fermionic symmetries.

The permutations
$(2\,6)^{p_{1}}(4\,8)^{p_{2}}(10\,14)^{p_{1}}(12\,16)^{p_{4}}$ with
$p_{i}=0,1$, map the twist vector
$v_{\text{L}}(\mathbf{1}\\!-\\!\mathbf{b})=\mbox{$\frac{1}{2}$}(0^{20})\rightarrow$

$\displaystyle|p_{1}p_{2}\,p_{3}p_{4})=\mbox{$\frac{1}{2}$}(p_{1}p_{2}p_{1}p_{2}\,p_{3}p_{4}p_{3}p_{4}\,0^{12})\,,$
(11)

while the other twists and shifts remain the same, since $(2a\,2b)$ leave
$V_{\text{L}}$ inert (see Table 2). When these permutations are successively
switched on, the T–fold web, given in Figure 1, is obtained.

For the cases with two non–zero $p_{i}$, (11) implies that
$v_{\text{L}}=(v_{\text{R}},0^{16})$ possibly up to a change of basis. Thus
the resulting bosonic models are interpreted as symmetric orbifolds. In
particular, the model obtained after the fermionic permutation $(2\,6)(4\,8)$
is conventionally considered as the bosonic representation of this fermionic
model in which ${\boldsymbol{\xi}}$ just separates out the $SO(32)$ gauge
group, while for all the other cases with two non–zero $p_{i}$,
${\boldsymbol{\xi}}$ acts as an asymmetric Wilson line.

Table 4 provides an overview of all inequivalent T–fold models associated with
this fermionic model. It indicates in how many directions $\mathbf{b}$ and
${\boldsymbol{\xi}}$ act as left–moving twists. Apart from the sixteen models
depicted in Figure 1 (of which nine are inequivalent), ${\boldsymbol{\xi}}$
has an asymmetric twist action as can inferred from this table. The total
number of inequivalent T–fold models associated with the fermionic basis
vectors given in Table 3 is $213$. This number rapidly increases for fermionic
models defined with more basis vectors. For example, for the fermionic model
in which ${\boldsymbol{\xi}}$ is split into ${\boldsymbol{\xi}}_{1}$ and
${\boldsymbol{\xi}}_{2}$ the number of inequivalent bosonisations becomes
$11\,273$ and for the NAHE set [25, 26, 27] is $85\,735$.

Figure 1:  Web of T–folds associated with the fermionic model given in Table 3 in which only $\mathbf{b}$ acts as a twist (11).  Table 4:  The number of $\mathbf{b}$ and ${\boldsymbol{\xi}}$ twists for the inequivalent T–folds associated with the fermionic model in Table 3.  $\boldsymbol{~{}\mathbf{b}~{}}$ $\boldsymbol{{\boldsymbol{\xi}}}$ | 0 | 2 | 4 | 6 | 8 | Sum
---|---|---|---|---|---|---
0 | 1 | 2 | 3 | 2 | 1 | 9
2 | 2 | 11 | 18 | 12 | 3 | 46
4 | 3 | 18 | 32 | 19 | 6 | 78
6 | 2 | 12 | 19 | 18 | 7 | 58
8 | 1 | 3 | 6 | 7 | 5 | 22
Sum | 9 | 46 | 78 | 58 | 22 | 213

##  VI  Moduli

The unfixed Narain moduli ($m_{ij}=g_{ij}+b_{ij}$ with metric $g_{ij}$,
B–field $b_{ij}$ and Wilson lines $m_{ix=1\ldots 16}$) of a T–fold correspond
to the operators,

$\displaystyle m_{ia}\,\partial
X_{\text{R}}^{i}\,\bar{\partial}X_{\text{L}}^{a}\,,$ (12)

left inert by (3). Symmetric orbifolds always leave at least the diagonal
metric moduli $m_{ii}=g_{ii}$ free, asymmetric orbifolds may fix all moduli.

This would suggest that the number of frozen moduli may vary dramatically
depending on which bosonic description of a given fermionic model is used.
There is no paradox here either: the unfixed scalar deformations of the
fermionic model can be identified by the Thirring interactions,

$\displaystyle m_{iuv}\,y^{i}w^{i}\,{\bar{f}}^{u}\\!{\bar{f}}^{v}\,,$ (13)

left inert by (6). Thus, the total number of massless untwisted scalars is
bosonisation independent, and therefore identical in any bosonic realisation.
Which of them are interpreted as free Narain deformations, however, does
depend on the choice of bosonisation, as $X_{\text{L}}$ in (12) does.

##  VII  Intrinsically asymmetric T–folds

The previous section showed that whether a real fermionic model should be
considered as a symmetric or asymmetric model is very much bosonisation
dependent. A free fermionic model is called intrinsically asymmetric if for
any bosonisation it corresponds to an asymmetric orbifold. An intrinsically
asymmetric T–fold is a bosonic model associated with an intrinsically
asymmetric fermionic model.

In light of the observation below (12), a fermionic model that admits a
symmetric interpretation has at least inert Thirring interactions (13) with
different $u\neq v$ for each $i$. If not, the fermionic model is intrinsically
asymmetric and hence in any bosonic realisation all Narain moduli are frozen.
This is, in particular, the case when no Thirring interactions (13) are
invariant under (6). An example of such a model is given in ref. [28].

Simple examples of intrinsically asymmetric free fermionic models can be
obtained by taking basis vectors that act as purely holomorphic twists. (For
example, consider the twist basis vector $\mathbf{b}=\big{\\{}\chi^{1\ldots
4},y^{1\ldots 4}\big{\\}}$ in 6D or $\mathbf{b}_{1}=\big{\\{}\chi^{1\ldots
4},y^{1\ldots 4}\big{\\}}$ and $\mathbf{b}_{2}=\big{\\{}\chi^{3\ldots
6},y^{3\ldots 6}\big{\\}}$ in 4D.) As there are no invariant Thirring
interactions (13) possible, the corresponding T–fold models are necessarily
intrinsically asymmetric.

##  VIII  Discussion

This letter focused on heterotic T–folds that admit fermionic descriptions.
Even though the key observation that bosonisation in a fermionic CFT is not
unique is not new, its striking consequences seem not to have been appreciated
so far: a single free fermionic model can be associated with a large number of
seemingly unrelated bosonic theories. Some may admit a symmetric orbifold
interpretation while most others are asymmetric, but in many different ways.

In light of this, studies of non–geometric constructions, and T–folds in
particular, may need to be revised, since seemingly different non–geometries
may, in fact, be equivalent. In particular, in the bosonic orbifold literature
it would be inconceivable that symmetric and asymmetric orbifolds can be
identified. Moreover, the number of frozen moduli turns out to be a
bosonisation dependent quantity; only the total number of massless untwisted
scalars is identical in any description. Only for an intrinsically asymmetric
T–fold all Narain moduli are fixed in any bosonic description. In addition,
the induced bosonic actions of fermionic symmetries hints at an extension of
the T–duality group of toroidal and $\mathbbm{Z}_{2}$ orbifold
compactifications.

The findings presented here were derived at free fermionic points. However,
the induced transformations of the bosonic boundary conditions by the
fermionic symmetries may be considered without referring to the fermionic
description. Hence, it is an interesting question whether the suggested
extension of the T–duality group discussed above is a general duality symmetry
of string theory or exists at free fermionic points only.

Our analysis is partially motivated by quasi-realistic model building using
the free fermionic formulation, see e.g. [28], to give rise to some central
features of the Standard Model and its supersymmetric extensions, such as the
existence of three generations charged under the Standard Model gauge group
with potentially viable Yukawa couplings to Higgs doublets. While this paper
focused on the moduli of the internal manifolds, there exist free fermionic
models in which the moduli space is further restricted [29], which shows the
need for a deeper understand the moduli space in these quasi–realistic example
which our analysis may provide. Moreover, the methods adopted in the
supersymmetric cases considered here can also be utilised in
non–supersymmetric string constructions as well as in tachyon free models that
are obtained from compactifications of tachyonic ten dimensional string vacua
[30]. In this respect, the understanding of the correspondence between the
fermionic and bosonic representations of string vacua is essential to obtain a
more profound understanding of the string dynamics at the Planck scale.

##  IX acknowledgments

AEF would like to thank the CERN theory division for hospitality and support.
SGN would like to thank the organizers of the workshop Gauged Linear Sigma
Models @ 30 for their kind invitation and the warm hospitality of the Simons
Center at Stony Brook University during this event where part of this work was
done. The work of BP is supported by EPSRC grant EP/W522399/1.

## References

  * Dixon  _et al._ [1985] L. J. Dixon, J. A. Harvey, C. Vafa, and E. Witten, Strings on orbifolds, Nucl. Phys. B261, 678 (1985).
  * Dixon  _et al._ [1986] L. J. Dixon, J. A. Harvey, C. Vafa, and E. Witten, Strings on orbifolds. 2, Nucl. Phys. B274, 285 (1986).
  * Faraggi  _et al._ [1990] A. E. Faraggi, D. V. Nanopoulos, and K.-j. Yuan, A Standard like model in the 4D free fermionic string formulation, Nucl. Phys. B335, 347 (1990).
  * Lebedev  _et al._ [2007] O. Lebedev _et al._ , A mini-landscape of exact MSSM spectra in heterotic orbifolds, Phys. Lett. B645, 88 (2007), arXiv:hep-th/0611095 .
  * Blaszczyk  _et al._ [2010] M. Blaszczyk _et al._ , A Z2xZ2 standard model, Phys. Lett. B683, 340 (2010), arXiv:0911.4905 [hep-th] .
  * Duff [1990] M. J. Duff, Duality Rotations in String Theory, Nucl. Phys. B 335, 610 (1990).
  * Strominger  _et al._ [1996] A. Strominger, S.-T. Yau, and E. Zaslow, Mirror symmetry is T duality, Nucl. Phys. B 479, 243 (1996), arXiv:hep-th/9606040 .
  * Dabholkar and Hull [2006] A. Dabholkar and C. Hull, Generalised T-duality and non-geometric backgrounds, JHEP 05, 009, arXiv:hep-th/0512005 .
  * Grana  _et al._ [2009] M. Grana, R. Minasian, M. Petrini, and D. Waldram, T-duality, Generalized Geometry and Non-Geometric Backgrounds, JHEP 04, 075, arXiv:0807.4527 [hep-th] .
  * Narain  _et al._ [1987] K. S. Narain, M. H. Sarmadi, and C. Vafa, Asymmetric orbifolds, Nucl. Phys. B288, 551 (1987).
  * Ibáñez  _et al._ [1988a] L. E. Ibáñez, J. Mas, H. P. Nilles, and F. Quevedo, Heterotic strings in symmetric and asymmetric orbifold backgrounds, Nucl. Phys. B301, 157 (1988a).
  * Di Francesco  _et al._ [1997] P. Di Francesco, P. Mathieu, and D. Senechal,  _Conformal field theory_ (Springer, New York, USA, 1997) 890 p.
  * Athanasopoulos  _et al._ [2016] P. Athanasopoulos, A. E. Faraggi, S. Groot Nibbelink, and V. M. Mehta, Heterotic free fermionic and symmetric toroidal orbifold models, JHEP 04, 038, arXiv:1602.03082 [hep-th] .
  * Faraggi [1994] A. E. Faraggi, Proton stability in superstring derived models, Nucl. Phys. B 428, 111 (1994), arXiv:hep-ph/9403312 .
  * Faraggi [2001] A. E. Faraggi, Doublet–triplet splitting in realistic heterotic string derived models, Physics Letters B 520, 337 (2001).
  * Faraggi [1993] A. E. Faraggi, Yukawa couplings in superstring-derived standard-like models, Phys. Rev. D 47, 5021 (1993).
  * Faraggi [2005] A. E. Faraggi, Moduli fixing in realistic string vacua, Nuclear Physics B 728, 83 (2005).
  * Gross  _et al._ [1985] D. J. Gross, J. A. Harvey, E. J. Martinec, and R. Rohm, Heterotic string theory. 1. The free heterotic string, Nucl. Phys. B256, 253 (1985).
  * Ibáñez  _et al._ [1988b] L. E. Ibáñez, J. Mas, H.-P. Nilles, and F. Quevedo, Heterotic strings in symmetric and asymmetric orbifold backgrounds, Nucl. Phys. B301, 157 (1988b).
  * Groot Nibbelink and Vaudrevange [2017] S. Groot Nibbelink and P. K. S. Vaudrevange, T-duality orbifolds of heterotic Narain compactifications, JHEP 04, 030, arXiv:1703.05323 [hep-th] .
  * Groot Nibbelink [2021] S. Groot Nibbelink, A worldsheet perspective on heterotic T-duality orbifolds, JHEP 04, 190, arXiv:2012.02778 [hep-th] .
  * Antoniadis  _et al._ [1987] I. Antoniadis, C. Bachas, and C. Kounnas, Four-dimensional superstrings, Nucl. Phys. B289, 87 (1987).
  * Kawai  _et al._ [1987] H. Kawai, D. C. Lewellen, and S.-H. Henry Tye, Construction of fermionic string models in four dimensions, Nucl. Phys. B288, 1 (1987).
  * Antoniadis and Bachas [1988] I. Antoniadis and C. Bachas, 4d fermionic superstrings with arbitrary twists, Nucl. Phys. B298, 586 (1988).
  * Ferrara  _et al._ [1987] S. Ferrara, L. Girardello, C. Kounnas, and M. Porrati, The Effective Interactions of Chiral Families in Four-dimensional Superstrings, Phys. Lett. B 194, 358 (1987).
  * Antoniadis  _et al._ [1988] I. Antoniadis, J. R. Ellis, J. S. Hagelin, and D. V. Nanopoulos, GUT Model Building with Fermionic Four-Dimensional Strings, Phys. Lett. B 205, 459 (1988).
  * Faraggi and Nanopoulos [1993] A. E. Faraggi and D. V. Nanopoulos, Naturalness of three generations in free fermionic $Z_{2}^{n}\otimes Z_{4}$ string models, Phys. Rev. D48, 3288 (1993).
  * Faraggi [1992] A. E. Faraggi, A new standard-like model in the four dimensional free fermionic string formulation, Physics Letters B 278, 131 (1992).
  * Cleaver  _et al._ [2008] G. B. Cleaver, A. E. Faraggi, E. Manno, and C. Timirgaziu, Quasi-realistic heterotic-string models with vanishing one-loop cosmological constant and perturbatively broken supersymmetry?, Phys. Rev. D 78, 046009 (2008).
  * Faraggi  _et al._ [2020] A. Faraggi, V. Matyas, and B. Percival, Stable three generation standard-like model from a tachyonic ten dimensional heterotic-string vacuum, The European Physical Journal C 80, 10.1140/epjc/s10052-020-7894-x (2020).