2021

[1]Angelina Senchukova

1]School of Engineering Sciences, Lappeenranta–Lahti University of Technology
LUT, Yliopistonkatu 34, Lappeenranta, FI-53850, South Karelia, Finland

2]Finnos Oy, Tukkikatu 5, Lappeenranta, FI-53900, South Karelia, Finland

# Geometry Parameter Estimation for Sparse X-ray Log Imaging

<EMAIL_ADDRESS>Jarkko Suuronen<EMAIL_ADDRESS>Jere
Heikkinen<EMAIL_ADDRESS>Lassi Roininen<EMAIL_ADDRESS>[  [

###### Abstract

We consider geometry parameter estimation in industrial sawmill fan-beam X-ray
tomography. In such industrial settings, scanners do not always allow
identification of the location of the source-detector pair, which creates the
issue of unknown geometry. This work considers an approach for geometry
estimation based on the calibration object. We parametrise the geometry using
a set of 5 parameters. To estimate the geometry parameters, we calculate the
maximum cross-correlation between a known-sized calibration object image and
its filtered backprojection reconstruction and use differential evolution as
an optimiser. The approach allows estimating geometry parameters from full-
angle measurements as well as from sparse measurements. We show numerically
that different sets of parameters can be used for artefact-free
reconstruction. We deploy Bayesian inversion with first-order isotropic Cauchy
difference priors for reconstruction of synthetic and real sawmill data with a
very low number of measurements.

###### keywords:

geometry calibration, fan-beam computed tomography, Bayesian inverse problems,
differential evolution

##  1 Introduction

A particular case of sparse X-ray computed tomography (CT) problems is found
in the sawmill industry, where X-ray scanners are employed for non-invasive
detection of log features and defects such as knots, rotten parts, and foreign
objects. Furthermore, X-ray scanners allow timber tracking from the raw
material to the resulting product (wood fingerprinting) Zol2019 ; Flod2008 .
Non-destructive feature detection and timber tracking provide information for
selection of optimal sawing strategies and efficient process control.

CT image reconstruction can be described mathematically as an inverse problem
Tar2004 ; Kaip2005 ; Stu2010 . The objective of two-dimensional X-ray CT is to
reconstruct the cross-sectional image of a scanned object from a collection of
one-dimensional line projections. The projections are obtained by measuring
the intensity loss of the beam of X-rays that penetrates the object from
different view angles. The projections are collected to a sinogram.

In the sawmill industry, highly limited (sparse) tomographic data is used for
log imaging. The sparse tomography is known to be a severely ill-posed
problem. Sparse data lead to severe artefacts when the image is reconstructed
by conventional methods like filtered backprojection Nat1986 .

To solve ill-posed inverse problems, a Bayesian approach is often applied
Stu2010 . The Bayesian framework naturally allows for incorporating a prior
distribution which encodes the possible internal structure of the object.
First-order isotropic Cauchy difference priors are shown to promote sharp
edges, which enables detection of sharp features and material interfaces with
a very low number of measurements in the sawmill applications Suur2020 .

The CT scanner geometry describes the location and orientation of the source-
detector pair for each projection of the CT scan with respect to the reference
system of the scanned object. Geometric inaccuracies are known to cause errors
and severe artefacts in reconstructions LiGull1994 ; Huili1998 . In practical
industrial settings in sawmills, the structure of an X-ray scanning device
does not always allow for the direct measurement of the important geometry
parameters such as the distance between the source and detector or the
detector tilt. Moreover, geometric errors can arise due to factors such as
vibration during motion of either the scanner or scanned object Jereb2012 .
There is, therefore, a need to address the issue of unknown geometry
parameters.

In this work, we estimate the fan-beam geometry defined by 5 parameters: the
first scanning angle, the distance between the centre of rotation and the
detector, the source shift, the detector shift, and the detector tilt. In the
particular application to sawmill log imaging discussed in the paper, a single
source-detector pair is fixed while the scanned object is rotated to obtain
its projections at different angles.

Figure 1: A scheme of the X-ray measurement geometry. The point source $S$ and
the flat-line detector $D$ are placed around the log. The log is rotated
around the centre of rotation COR to obtain different projections. The
distance between the source and COR is $r_{S}$ and the distance between COR
and the detector is $r_{D}$. In the geometry, there are possible shifts of
source $h_{S}$ and detector $h_{D}$ in tangent directions, and the detector
tilt $\alpha_{D}$

###  1.1 Related Work

Existing studies on geometry calibration can be roughly divided into two
groups: 1) methods using calibration phantoms, and 2) self-calibration
(marker-free) methods. Calibration phantom-based methods, which are the more
common approach, use one or more calibration objects (calibration phantoms)
incorporating a known arrangement of radiopaque markers such as steel ball
bearings. An extensive review of such approaches is given in Ferr2015 . An
overview of marker-free methods for estimating the crucial geometry parameter,
the axis of rotation, is given in Zem2023 .

A pioneering study Gull1987  on imaging system geometry estimation utilises
single photon emission computed tomography. In the study, the Levenberg-
Marquardt optimisation algorithm is used to calibrate a fan-beam geometry from
projection measurements of reference objects consisting of point markers. The
same principles are also appropriate for calibrating X-ray CT systems.

A calibration method for cone-beam geometry that employs a calibration phantom
to estimate the geometry parameters is proposed in Cho2005 . The method uses a
transparent cylindrical acrylic tube with 24 steel ball bearings as a
calibration phantom, but the approach can be adapted also to other phantom
configurations.

A self-calibration method is introduced in Ouadah2016 , where calibration is
done in the cone-beam setting without calibration objects. The method performs
calibration by registering two-dimensional projection data to a previously
acquired three-dimensional image of the scanned object.

Joint estimation of uncertain view angles determining the geometry of the
forward CT model and the attenuation coefficient function of a scanned object
is discussed in Uribe2021 ; Hansen2021 . In Ped2023 , the geometry estimation,
in particular, the center-of-rotation offset, is performed jointly with
reconstruction based on a Bayesian approach that also provides information on
the uncertainty of geometry parameters.

In Gen2021 , data-driven estimation of unknown fan-beam geometry is done by
employing a neural network that learns the unknown forward operator from
training data consisting of sinogram-image pairs. Another learned approach to
tackle parametric uncertainties in the measurement operator is introduced in
Xie2020 . The approach allows for joint automatic geometry parameter
calibration along with reconstruction of the unknown CT image.

###  1.2 Our Contributions

We propose a new approach for identifying unknown geometry parameters in fan-
beam X-ray log tomography by using a calibration object of known size that is
scanned jointly with a log.

Compared to other existing works on geometry calibration, the difference in
this approach is that: 1) it employs an easy-to-produce phantom for
calibration of X-ray CT geometry with a fixed source-detector pair and
rotating scanned object; 2) the approach itself is based on calculating the
maximum cross-correlation between the ground truth image of the calibration
object and its filtered backprojection (FBP) reconstruction; and 3) the
approach allows estimating geometry parameters even for a limited number of
projection angles.

Additionally, we demonstrate a number of possible artefacts that may appear in
the reconstructions in the eventuality of misspecified geometry parameters
when using a fan-beam configuration with a flat line detector. In extreme
cases, misspecifications can lead to full unidentifiability, and in less
severe cases, the reconstructions can contain duplicated features, halos,
shifts, and blurring. We aim to demonstrate these features and present a way
to overcome the issue with suitable parameter choices.

Finally, we consider the sparse-angle tomography and employ maximum a
posteriori estimates with a recently developed class of priors, first-order
isotropic Cauchy difference priors Suur2022 . We demonstrate the amenability,
for our specific use case, of this method over conventional methods (filtered
backprojection, Tikhonov regularisation) when the number of projection data is
extremely limited.

###  1.3 Outline of the Paper

The rest of this paper is organised as follows: In Section 2, we review the
tomography reconstruction model and methods. The new geometry parameter
estimation algorithm is presented in Section 3. Synthetic and real data CT
examples are demonstrated in Section 4. Section 5 concludes the study.

##  2 Tomography reconstruction

###  2.1 Fan-beam geometry setting

The schematic plot of the X-ray scanner geometry used in the numerical
experiments is shown in Figure 1. An X-ray source, placed around the log,
emits fan-shaped X-ray beams that are recorded by the corresponding flat line
detector.

The geometry of the fan-beam model is determined by seven parameters in total.
There are two known parameters:

  * •

$r_{S}=859.46$ (mm) — distance between the X-ray source and centre of rotation
(COR), and

  * •

$n_{D}=768$ (px) — number of detector elements, the pixel size is $2$mm,

and the other five parameters are treated as unknown:

  * •

$\alpha_{0}$ (rad) — the first scanning angle; the projection angles are
defined w.r.t. the first scanning angle, and are assumed to be uniformly
distributed between the first and last angles,

  * •

$r_{D}$ (mm) — distance between the COR and detector,

  * •

$h_{S}$ (mm)— the source shift (in the tangent direction),

  * •

$h_{D}$ (mm)— the detector shift (in the tangent direction), and

  * •

$\alpha_{D}$ — defines the detector tilt. In the default configuration, the
initial source-to-detector vector is $[0,1]$, and the initial detector axis is
perpendicular to it (i.e. $[1,\alpha_{D}=0]$). However, if the detector is
tilted then component $\alpha_{D}$ should be adjusted accordingly.

We further denote those five unknown parameters as a 5-dimensional vector

$\boldsymbol{\theta}=[\alpha_{0},r_{D},h_{S},h_{D},\alpha_{D}]^{T}.$

###  2.2 Discrete CT imaging model

Using the Beer-Lambert law Nat1986 , the intensity loss of X-ray beams passing
through a scanned object can be represented as

$\frac{I_{\varphi,\tau}}{I_{0}}=\exp\left(-\int_{l_{\varphi,\tau}}x(s)ds\right),$
(1)

where $x$ is the attenuation coefficient describing the internal structure of
the object, $I_{0}$ and $I_{\varphi,\tau}$ are the initial and final
intensities, respectively, and $l_{\varphi,\tau}$ denotes the line
corresponding to the acquisition angle $\varphi\in[0,2\pi)$ and detector
location $\tau\in\mathbb{R}$.

Taking logarithms in Equation (1) we obtain the following

$-\ln\left(\frac{I_{\varphi,\tau}}{I_{0}}\right)=\int_{l_{\varphi,\tau}}x(s)ds,$
(2)

where the left-hand side $y:=-\ln(I_{\varphi,\tau}/I_{0})$ is measured data
known as projection. The collection of projections is commonly referred to as
a sinogram. The right-hand side of (2) is known as the Radon transform

$\mathcal{R}(x):=\int_{l_{\varphi,\tau}}x(s)ds,$

that maps the unknown attenuation coefficient function $x$ into the set of its
line integrals Rad1986 . In other words, the known measured data is
represented by the Radon transform, which corresponds to the forward CT model.

In practice, only a finite number of X-ray projections is available and the
sinogram measurements are recorded by the pixel detectors. If the number of
acquisition angles is $K$, the number of detector pixels is $M$, and a scanned
object is discretised into a $N\times N$ grid, we obtain the discrete model
for the X-ray CT imaging

$\mathbf{y}=\mathbf{A}_{\boldsymbol{\theta}}\mathbf{x}+\mathbf{e},$ (3)

where $\mathbf{x}\in\mathbb{R}^{N^{2}\times 1}$ is the unknown vector of
attenuation coefficients,
$\mathbf{A}_{\boldsymbol{\theta}}\in\mathbb{R}^{KM\times N^{2}}$ is a so-
called system matrix (forward matrix) describing intersection lengths for each
of $KM$ lines. The system matrix is assumed constant throughout the
reconstruction and depends on the scanning geometry parameters
$\boldsymbol{\theta}$ that are not known exactly. Vector
$\mathbf{y}\in\mathbb{R}^{KM\times 1}$ specifies the measurement data
(projections) taken during the scanning process, and $\mathbf{e}$ denotes
measurement noise.

The corresponding inverse problem aims to reconstruct $\mathbf{x}$ from the
available measured sinogram data $\mathbf{y}$. This can be done either in
Bayesian estimation or a classical deterministic framework Kaip2005 . In the
paper, we use and compare three different reconstruction methods (filtered
backprojection, Tikhonov regularisation and MAP estimate) discussed in Section
2.3.

###  2.3 Reconstruction methods

The traditional method to reconstruct the cross-sectional image of an object
from the measurement data is filtered backprojection (FBP). If $\mathbf{y}$
denotes the measurement data and $\boldsymbol{\theta}$ is the geometry
parameter vector, then the discretised form of the FBP formula is defined as

$\mathbf{x}_{\mathrm{FBP}}=\frac{1}{2}\mathbf{B}_{\boldsymbol{\theta}}(\boldsymbol{\phi}\bar{\ast}\mathbf{y})=:\mathrm{FBP}(\boldsymbol{\theta}\mathchar
24635\relax\;\mathbf{y}),$ (4)

where $\boldsymbol{\phi}$ is a convolutional filter in the projection domain,
$\mathbf{B}_{\boldsymbol{\theta}}$ is discrete backprojection, and
$\bar{\ast}$ denotes the discrete convolution.

Several frequency domain filters are commonly used in FBP, including the most
basic Ram-Lak filter Kak2001  and its windowed modifications such as Shepp-
Logan, Cosine, Hamming and Hann filters. The filter function significantly
affects the robustness of the reconstruction and should be carefully selected
based on the particular reconstruction task.

In sparse CT where the number of projections is extremely limited, the FBP
method leads to severe artefacts in reconstructions due to ill-posedness of
the problem Nat1986 . Different regularisation techniques have been used to
solve ill-posed problems. One of the most common regularisation techniques,
Tikhonov regularisation, is defined as a solution of the following
minimisation problem

$\mathbf{x}_{\mathrm{Tikh}}=\operatorname*{arg\,min}_{\mathbf{x}}\left\\{\|\mathbf{A}_{\boldsymbol{\theta}}\mathbf{x}-\mathbf{y}\|^{2}+\alpha\|\mathbf{x}\|^{2}\right\\}.$
(5)

The regularisation parameter $\alpha>0$ is chosen depending on the noise
strength to find a trade-off between the solution stability and its ability to
explain the measurements.

In the Bayesian formulation of inverse problem (3), the observed data
$\mathbf{y}$, the unknown parameter $\mathbf{x}$ and noise $\mathbf{e}$ are
treated as random variables. The noise $\mathbf{e}\sim
N(\mathbf{0},\mathbf{\Sigma})$ is assumed to be a zero-mean Gaussian
multivariate with covariance matrix $\mathbf{\Sigma}$.

According to Bayes’ theorem, the posterior distribution of unknown
$\mathbf{x}$ can be formulated as

$\pi_{\text{post}}(\mathbf{x})=\pi(\mathbf{x}|\mathbf{y})=\frac{\pi_{\text{pr}}(\mathbf{x})\pi(\mathbf{y}|\mathbf{x})}{\pi(\mathbf{y})}.$

The prior distribution $\pi_{\text{pr}}(\mathbf{x})$ incorporates a priori
information about the unknown $\mathbf{x}$ before the measurements were
obtained, and $\pi(\mathbf{y}|\mathbf{x})$ is the data-likelihood that
contains information about the forward operator
$\mathbf{A}_{\boldsymbol{\theta}}$ and assumptions about the statistical
distribution of noise $\mathbf{e}$ Kaip2005 ; Stu2010 .

The features on the boundaries between different regions are of particular
interest in reconstructions. In sawmill applications, sharp-interface
reconstructions are needed for measuring the size and shape of tree knots.
Various edge-preserving regularisation techniques have been proposed to avoid
smoothing of edges in CT Yu2002 ; Xu2019 ; Hama2013 .

As an edge-preserving method, we employ the maximum a posteriori (MAP)
estimate of the posterior distribution

$\mathbf{x}_{\text{MAP}}=\operatorname*{arg\,max}_{\mathbf{x}\in\mathbb{X}}\pi_{\text{post}}(\mathbf{x}),$
(6)

with the first-order isotropic Cauchy difference prior Suur2020 ; Suur2022 ;
Markkanen19  and the following data likelihood term

$\pi(\mathbf{y}|\mathbf{x})\propto\exp\left(-\frac{1}{2}\|\mathbf{A}_{\boldsymbol{\theta}}\mathbf{x}-\mathbf{y}\|^{2}\right).$

Unlike the standard Tikhonov regularisation (5), which assumes no special
structure for $\mathbf{x}$ except that it should be close to zero, the first-
order Cauchy difference prior incorporates prior information on how the first-
order derivatives of $\mathbf{x}$ will be distributed. Employing Cauchy
distribution for the distribution of the differences effectively allows the
existence of discontinuities in $\mathbf{x}$ by favouring features that are
piecewise close to constant.

If $\mathbf{x}$ is assumed to be discretised in a two-dimensional uniform grid
indexed by $(i,j)$ with size $N\times N$, then the isotropic Cauchy difference
prior distribution can be expressed as follows Suur2022

$\begin{split}\pi_{\textrm{pr}}(\mathbf{x})&\propto\prod_{i=1}^{N-1}\prod_{j=1}^{N-1}\big{(}\beta^{2}+(x_{i+1,j}-x_{i,j})^{2}+\\\
&(x_{i,j+1}-x_{i,j})^{2}\big{)}^{-3/2},\end{split}$ (7)

where $\beta^{2}$ is a parameter that controls the strength of the prior in
terms of how close to zero the first-order differences should be.

Figure 2: Images of the L-shaped calibration phantom (left) and the
calibration phantom with a hole (right). All dimensions are given in mm

##  3 Geometry parameter estimation

For geometry calibration, we use two different calibration phantoms of known
size and shape: an L-shaped calibration phantom and a calibration phantom with
a hole (see Figure 2).

###  3.1 Geometry estimation model

The parameters of the partially unknown geometry of the measurement equipment
are estimated using an objective function that evaluates the correlation of a
reference image to an FBP reconstruction with given parameters. More
precisely, we compare both the reference image and its mirrored image to the
reconstruction, because the ambiguity of the direction of the detector plates
and the rotation direction of the object inside the scanning device may render
the reconstructions mirrored as well. The ambiguity of the handedness of the
reconstructions can be taken into account in the objective function by
calculating the maximum of correlations of both handednesses.

We seek to find a vector of unknown geometry parameters $\boldsymbol{\theta}$
based on the calibration phantom image $\mathbf{X}_{\mathrm{ref}}$ of $N\times
N$ pixels and its mirrored image $\widehat{\mathbf{X}}_{\mathrm{ref}}$. To
separate the reference image from its flattened vector form
$\mathbf{x}_{\mathrm{ref}}$, we denote the image ($N\times N$ pixels) with the
capitalised $\mathbf{X}_{\mathrm{ref}}$.

Let
$\mathbf{f}_{\boldsymbol{\theta}}=\mathrm{FBP}(\boldsymbol{\theta}\mathchar
24635\relax\;\mathbf{y})$ be the corresponding FBP reconstruction that depends
on the input X-ray (sinogram) data $\mathbf{y}$ and geometry parameters
$\boldsymbol{\theta}$. Then the objective function $J(\boldsymbol{\theta})$ is
defined as the negative maximum of the correlation between
$\mathbf{X}_{\mathrm{ref}}$ and $\mathbf{f}_{\boldsymbol{\theta}}$ and the
correlation between $\widehat{\mathbf{X}}_{\mathrm{ref}}$ and
$\mathbf{f}_{\boldsymbol{\theta}}$

$J(\boldsymbol{\theta})=-\max\left\\{\mathrm{Tr}\left({\mathbf{X}}_{\mathrm{ref}}^{T}\mathbf{f}_{\boldsymbol{\theta}}\right),\mathrm{Tr}\left(\widehat{\mathbf{X}}_{\mathrm{ref}}^{T}\mathbf{f}_{\boldsymbol{\theta}}\right)\right\\},$
(8)

where $\mathrm{Tr}$ means the trace of a matrix, and the superscript $T$ means
transpose. In other words, the objective function evaluates the Frobenius
inner product of the FBP reconstruction with given geometry parameters and two
possible calibration object candidates. It returns the negated maximum of the
values. In this study, we do not employ the normalised correlation in the
method, because the unnormalised correlation proved effective and robust. For
the second, usage of the normalised correlation requires addressing the
exceptional case when the standard deviation of the reconstruction is zero.
This might be the case when the differential evolution proposes a set of
geometry parameters that results in an empty reconstruction in the domain of
the reference object, although such a case is very unlikely.

We also experimented with an objective function based on sinogram data
directly rather than FBP, namely

$\tilde{J}(\boldsymbol{\theta})=-\max\left\\{\mathrm{}\left(\mathbf{A}_{\boldsymbol{\theta}}{\mathbf{x}}_{\mathrm{ref}}\right)^{T}\mathbf{y},\mathrm{}\left(\mathbf{A}_{\boldsymbol{\theta}}\widehat{\mathbf{x}}_{\mathrm{ref}}\right)^{T}\mathbf{y}\right\\},$
(9)

where $\mathbf{A}_{\boldsymbol{\theta}}$ is the forward CT operator matrix
that depends on the geometry parameters $\boldsymbol{\theta}$. However, the
objective function (9) demonstrated worse performance in comparison with the
objective function (8). More specifically, FBP reconstructions using the set
of found geometry parameters contained severe artifacts and were often
misaligned with the reference object. We argue this is because the
optimisation schemes cannot propose and find good sets of geometry parameters
to align the reference and measured sinograms when the correlation is
performed in the projection space, where the differences between the perturbed
and unperturbed geometry parameters appear relatively subtle. Instead, the
effect of perturbation of the geometry parameters is typically more
substantial in the reconstruction space, which helps the optimiser to find the
potential global optimum of the cross-correlation, which, in turn, ensures the
proper alignment of the geometry.

From a numerical viewpoint, simultaneous evaluation of two calibration object
scenarios in the objective function (8) may pose additional challenges to the
optimisation algorithms as the likelihood of getting trapped to a local
minimum might increase. Evaluating the two calibration object cases separately
is plausible. However, even an objective function that evaluates the
correlation with one calibration object at times is non-differentiable, so we
use the objective function presented in Equation (8) for the unified
optimisation workflow.

###  3.2 Differential Evolution

The aim of geometry parameter estimation task is to find parameters
$\boldsymbol{\theta}^{*}$ that minimise the objective function
$J(\boldsymbol{\theta})$

$\boldsymbol{\theta}^{*}=\operatorname*{arg\,min}_{\boldsymbol{\theta}}J(\boldsymbol{\theta}).$
(10)

To find a solution to the optimisation problem (10), we use the differential
evolution (DE) method StornPrice2020 . The method belongs to evolutionary
algorithms that apply natural mechanisms of evolution theory such as mutation
and selection to evolve a solution to an optimisation problem.

Let $S\subset\mathbb{R}^{D}$ be the search space of the optimisation problem
to be solved. Then the DE algorithm keeps a population of $N_{\mathrm{pop}}$
individuals ($D$-dimensional vectors)

$\mathbf{x}_{i}=[x_{i,1},x_{i,2},\ldots,x_{i,D}]^{T},\quad
i=1,\ldots,N_{\mathrm{pop}},$

each of these individuals stands for a possible solution of the problem. At
every pass through the population, the algorithm mutates each candidate
solution by other candidate solutions to create a mutant candidate.

Different mutation strategies are used in DE Red2007 ; Leon2014 . In this
work, we adhere to the DE/best/1/bin strategy and use its implementation in
the SciPy Python library SciPy2020  (the pseudo-code is listed in Algorithm
1). In this strategy, the best member of the population,
$\mathbf{x}_{\mathrm{best}}$, is mutated (the Mutation stage) by the
difference of two other randomly selected members $\mathbf{x}_{k_{1}}$ and
$\mathbf{x}_{k_{2}}$, ($k_{1}\neq k_{2}\neq\mathrm{best}$) to form a mutant
vector $\mathbf{m}_{i}$ as follows

$\mathbf{m}_{i}=\mathbf{x}_{\mathrm{best}}+\mu(\mathbf{x}_{k_{1}}-\mathbf{x}_{k_{2}}),$

where $\mu\in[0,2]$ is the mutation constant that controls the amplification
of the differential variation $(\mathbf{x}_{k_{1}}-\mathbf{x}_{k_{2}})$.

Figure 3: A digital log phantom representing the cross-sectional image of a
tree log (image resolution $256\times 256$ pixels)

Next, at the Crossover stage, the trial vector is formed as follows

$u_{i,j}=\begin{cases}m_{i,j},&\mathrm{if}\;r_{i,j}<P_{\mathrm{cross}},\\\
x_{i,j},&\mathrm{otherwise},\end{cases}$

for $j=1,\ldots,D-1$, where $r_{i,j}\sim U(0,1)$ is a random number drawn from
the uniform distribution on the interval $[0,1]$, and
$P_{\mathrm{cross}}\in[0,1]$ is a crossover (recombination) probability
StornPrice2020 . In the numerical experiments, we use
$P_{\mathrm{cross}}=0.7$. The final entity $u_{i,D}$ is always calculated as

$u_{i,D}=m_{i,D}.$

Eventually, at the Selection stage, the fitness of the trial candidate is
assessed: if the trial candidate is better than the original candidate then it
takes its place.

Figure 4: Example artefacts in reconstructions when geometry parameters are
incorrectly specified: (a) the source and detector shifts, detector tilt are
misspecified; (b) the detector shift and tilt are misspecified; (c) the
detector tilt is misspecified; (d) the detector and source shifts are
misspecified; (e) the detector radius is slightly perturbed; (f) all the
parameters are specified correctly. FBP reconstruction (f) with true
parameters is used as a reference for computing $\varepsilon_{\mathrm{rel}}$
and SSIM Algorithm 1 Differential Evolution (DE)

1:Input parameters: $D$ – problem dimensionality; $N_{\mathrm{pop}}$ –
population size; $\mu$ – mutation constant; $P_{\mathrm{cross}}$ – crossover
probability; $f$ – objective function under consideration;

2:Initialisation:  Create population with $N_{\mathrm{pop}}$ individuals and
initialise them with random positions in the search space;

3:while the termination condition is not satisfied do

4: for $i=1..N_{\mathrm{pop}}$ do

5: Generation: Take the best individual

6: $\mathbf{x}_{\mathrm{best}}$ in the current population;

7: Choose 2 random integers

8: $k_{1}$, $k_{2}\in(1,N_{\mathrm{pop}})$ such that

9: $k_{1}\neq k_{2}\neq\mathrm{best}$;

10: Mutation: Form a mutant vector

11:
$\mathbf{m}_{i}\leftarrow\mathbf{x}_{\mathrm{best}}+\mu(\mathbf{x}_{k_{1}}-\mathbf{x}_{k_{2}})$;

12: Crossover:

13: for $j=1..D-1$ do

14: Generate $r_{i,j}\sim U(0,1)$;

15: if $r_{i,j}<P_{\mathrm{cross}}$ then

16: $u_{i,j}\leftarrow m_{i,j}$;

17: else

18: $u_{i,j}\leftarrow x_{\mathrm{best},j}$;

19: end if

20: end for

21: $u_{i,D}\leftarrow m_{i,D}$;

22: Selection:

23: if $f(\mathbf{u}_{i})\leq f(\mathbf{x}_{\mathrm{best}})$ then

24: $\mathbf{x}_{\mathrm{best}}\leftarrow\mathbf{u}_{i}$;

25: end if

26: end for

27:end while

The choice of the DE method for our problem of geometry parameter optimisation
stems from the fact that DE is able to cope with nonlinear and non-
differentiable objective functions and requires only a few control variables.
It is also a global optimisation algorithm, which is important for our
geometry parameter search. There are several local optima in the parameter
space when evaluating the correlation between a reference object and the
parametrised FBP reconstructions. Consequently, local optimisation algorithms
that do not require the objective function to be differentiable, such as the
Nelder-Mead method Nelder1965 , are not suitable for this task. It must be
recalled that the objective function likely does not even have a unique
optimum: there might be multiple parameter vectors that give almost the same
cross-correlation with the reference object when used in the FBP
reconstruction function. However, the uniqueness of the geometry parameters of
the measurement equipment is not interesting as such, since we are only
interested in the quality of reconstructions.

In addition, we claim that the nature of the objective function (8) favors DE
over simulated annealing (SA), which is another common global optimisation
algorithm SA1989 ; Horst2013 . Differential evolution utilises a set of
particles that mimic the natural selection, evolution and mutation of the best
candidates for the global maximum of the objective function. When some of the
geometry parameter candidates have been mutated closer to a possible optimum,
the evolution steps help to find the precise maximum in a more robust manner
than SA.

##  4 Results and discussion

In the numerical experiments, we used the Operator Discretization Library
(ODL) JonasAdler2017 . ODL is a Python library developed at KTH Royal
Institute of Technology, Stockholm. The library focuses on inverse problems
including CT reconstruction.

Figure 5: Full-angle FBP-reconstructions of log phantom (image resolution $256\times 256$ pixels): true geometry parameters (framed subplot) and different parameter vectors found by the calibration method with the L-shaped phantom. FBP reconstruction with true parameters is used as a reference for computing $\varepsilon_{\mathrm{rel}}$ and SSIM | Geometry parameters |  |
---|---|---|---
| $\alpha_{0}$ | $r_{D}$ | $h_{S}$ | $h_{D}$ | $\alpha_{D}$ | $\varepsilon_{\mathrm{rel}}$ | SSIM
| (rad) | (mm) | (mm) | (mm) |  |  |
True parameters | $2.55$ | $715$ | $320$ | $44$ | $0.28$ |  |
| $2.45$ | $743$ | $425$ | $135$ | $0.40$ | $0.073$ | $0.829$
| $2.73$ | $678$ | $157$ | $-89$ | $0.11$ | $0.064$ | $0.847$
L-shaped phantom | $3.00$ | $621$ | $-81$ | $-268$ | $-0.16$ | $0.072$ | $0.796$
| $2.60$ | $670$ | $266$ | $-8$ | $0.16$ | $0.083$ | $0.785$
| $2.89$ | $641$ | $21$ | $-192$ | $-0.08$ | $0.083$ | $0.764$
| $2.79$ | $641$ | $-134$ | $-89$ | $-0.02$ | $0.081$ | $0.693$
| $3.16$ | $593$ | $-223$ | $-370$ | $-0.34$ | $0.065$ | $0.834$
Phantom with a hole | $2.90$ | $638$ | $8$ | $-204$ | $-0.04$ | $0.122$ | $0.709$
| $2.57$ | $699$ | $305$ | $28$ | $0.24$ | $0.077$ | $0.789$
| $2.58$ | $716$ | $291$ | $22$ | $0.27$ | $0.045$ | $0.928$
Table 1:  Relative error $\varepsilon_{\mathrm{rel}}$ and SSIM for full-angle
FBP reconstructions of log phantom obtained with different parametrisations.
The parametrisations were found using the L-shaped calibration phantom and the
calibration phantom with a hole. True parameters are provided as a reference.
Geometry parameter vectors for the L-shaped phantom correspond to those in
Figure 5

We simulated the scanning geometry, computed forward projections, Tikhonov and
FBP reconstructions using ODL. To compute MAP estimates we used the low-memory
variant of the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm implemented
in Julia.

To assess the quality of reconstruction, we report relative reconstruction
error defined as

$\varepsilon_{\mathrm{rel}}=\frac{\|\hat{\mathbf{x}}-\mathbf{x}_{\mathrm{true}}\|_{2}}{\|\mathbf{x}_{\mathrm{true}}\|_{2}},$

where image reconstruction $\hat{\mathbf{x}}$ and reference image
$\mathbf{x}_{\mathrm{true}}$ are in the vector form. As the other image
quality metric, we employ the structural similarity index measure (SSIM) since
it is shown to mimics well the perceived quality of an image by the human
visual system WangBovik2004 ; WangBovik2009 . In numerical experiments, the
full-angle FBP reconstruction is used as a reference.

###  4.1 Synthetic data analysis

Here, we illustrate the proposed method for geometry calibration with an
application to a digital log phantom (Figure 3). The phantom consists of
concentric rings modelling tree growth rings, two ellipses representing knots
inside a tree trunk, and one small ellipse with higher density representing a
foreign object, for example, a metallic bullet.

Figure 6: Pairwise plot of geometry parameters
$\boldsymbol{\theta}=[\alpha_{0},r_{D},h_{S},h_{D},\alpha_{D}]^{T}$ and the
objective function $J(\boldsymbol{\theta})$. Optimal geometry parameter values
(in blue) were obtained by running the calibration method with full-angle (360
projections) measurements of the L-shaped calibration phantom 100 times. True
values of $\boldsymbol{\theta}$ and $J(\boldsymbol{\theta})$ are marked as red
dots Figure 7: Pairwise plot of geometry parameters
$\boldsymbol{\theta}=[\alpha_{0},r_{D},h_{S},h_{D},\alpha_{D}]^{T}$ and the
objective function $J(\boldsymbol{\theta})$. Optimal geometry parameter values
(in blue) were obtained by running the calibration method with sparse (20
projections) measurements of the L-shaped calibration phantom 100 times. True
values of $\boldsymbol{\theta}$ and $J(\boldsymbol{\theta})$ are marked as red
dots Figure 8: Full-angle reconstructions of log phantom. Left: geometry
parameters are estimated using full-angle measurements (360 projection
angles); right: geometry parameters are estimated using sparse measurements
(20 projection angles)

We simulated full-angle tomography measurements of the log phantom and images
of two calibration phantoms using ODL with true geometry parameter values
$[\alpha_{0},r_{D},h_{S},h_{D},\alpha_{D}]^{T}=[2.55,715,320,44,0.28]^{T}$. To
avoid the inverse crime Kaip2005 , we created the synthetic measurement data
on a grid with $N=1013$, then interpolated to a coarser grid with $N=256$ and
added relative error noise of $2\%$.

Figure 9: Reconstructions of the log phantom obtained by different
reconstruction methods with various numbers of projection angles (image
resolution $256\times 256$ pixels). The geometry parameters were obtained by
the calibration method using the full-angle measurements of the L-shaped
calibration phantom

Figure 4 shows the effect of misspecified geometry parameters on the quality
of full-angle FBP reconstructions of log phantom. We tested various
combinations of geometry parameters (zero value means that the corresponding
geometry parameter is misspecified, otherwise — specified correctly):

  1. (a)

$[\alpha_{0},r_{D},0,0,0]^{T}$

  2. (b)

$[\alpha_{0},r_{D},h_{S},0,0]^{T}$

  3. (c)

$[\alpha_{0},r_{D},h_{S},h_{D},0]^{T}$

  4. (d)

$[\alpha_{0},r_{D},0,0,\alpha_{D}]^{T}$

  5. (e)

$[\alpha_{0},r_{D}+\varepsilon,h_{S},h_{D},\alpha_{D}]^{T}$, where
$\varepsilon$ is a small positive number (pertubation);

  6. (f)

$[\alpha_{0},r_{D},h_{S},h_{D},\alpha_{D}]^{T}$.

In the figure, reconstructions corresponding to cases (b), (c) and (e) contain
severe artefacts such as duplicated features in the log phantom and halos,
whereas reconstructions with misspecified parameter vectors (a) and (d) are
completely unidentifiable.

To tackle the issue of unknown geometry parameters, we ran the parameter
search according to the method described in Section 3. In the parameter
search, we aimed to find optimal values for five unknown geometry parameters,
that is, $\boldsymbol{\theta}=[\alpha_{0},r_{D},h_{S},h_{D},\alpha_{D}]^{T}$,
based on simulated X-ray data of two calibration phantoms with true parameter
values that we set to
$[\alpha_{0},r_{D},h_{S},h_{D},\alpha_{D}]^{T}=[2.55,715,320,44,0.28]^{T}$.

Figure 5 shows the full-angle FPB reconstructions with different geometry
parameters found using the proposed calibration method. The geometry parameter
search was done using the simulated full-angle measurements of the L-shaped
calibration phantom, the reconstruction with the true geometry parameters is
provided in the figure as a reference. It can be noted that all the different
parametrisations provide satisfactory quality in the FBP reconstructions
although there is no unique optimal parameter set.

We also ran the geometry parameter search using the other calibration phantom
— phantom with a hole. The obtained results were similar to those of L-shaped
calibration phantom in Figure 5, so we omit the plot. Instead, in Table 1, we
report relative reconstruction error $\varepsilon_{\mathrm{rel}}$ and SSIM for
FBP reconstructions with different geometry parameters found by the geometry
parameter search using two calibration phantoms. It can be concluded that both
calibration phantoms perform well in the geometry parameter search algorithm
and the obtained estimates for the geometry parameters result in the
reconstructions of good quality.

Figure 10: Two-dimensional reconstructions of the log slice obtained by
different reconstruction methods with various numbers of projection angles
(image resolution $256\times 256$ pixels). The geometry parameters were
obtained by the calibration method using the full-angle measurements of the
L-shaped calibration phantom

To provide a view on the uncertainty of the found parameters, we ran the
optimisation algorithm with L-shaped phantom 100 times and created the
pairwise plot of optimal parameters together with the true parameter values
(see Figure 6). It can be noted that the parameters are highly correlated and
their distribution is very narrow. When using the calibration phantom with a
hole we observed similar behaviour of the parameters so we omitted the
corresponding plot.

In addition, we tested the performance of the calibration method when the
sparse tomography measurements of the calibration phantoms were used. In
Figure 7, we show the results of geometry parameter search using sparse
measurements (20 projection angles) of L-shaped phantom. The results are
similar to those in experiment in Figure 6, where the full-angle measurements
of the calibration phantom were used. In Figure 8, two full-angle FBP
reconstructions of log phantom are shown for comparison: geometry parameters
for the left reconstruction were obtained using the full-angle measurements of
L-shaped phantom while for the right reconstruction, sparse measurements of
the phantom were used. Red contour lines highlighting the phantom features are
used to show that there is a small angular misalignment between two
reconstructions. Despite the misalignment both reconstructions are of good
quality, therefore, the proposed calibration method works in the sparse-angle
scenario as well.

Although the scanning geometry is estimated correctly using the calibration
method, severe artefacts still appear in the reconstructions if classical
reconstruction approaches such as FBP and Tikhonov regularisation are used and
the number of projection angles is very low. To tackle the issue of artefacts,
we employed Bayesian inversion with edge-preserving Cauchy priors from Section
2.3. As shown in Figure 9, in the sparse case (20 projection angles), the FBP
reconstruction contains severe streaking artefacts and the Tikhonov
regularised reconstruction is extremely blurry, while the quality of the MAP
estimate with a Cauchy prior is still good enough and allows for detection of
knots and other features in the log phantom.

For MAP reconstructions in the figure, the relative error
$\varepsilon_{\mathrm{rel}}$ grows representing the expected quality
deterioration as the number of projection angles decreases. However, SSIM
rises suggesting the quality improvement, which contradicts the visual quality
assessment and $\varepsilon_{\mathrm{rel}}$. This phenomena can be explained
by the reduced ability of SSIM to represent the visual quality for CT
reconstructions PamNou2015  and the images containing salt-and-pepper
(impulse) noise LinChen2022 .

###  4.2 Real data analysis

In this section, the proposed method for geometry calibration is applied to
real X-ray data of a log. The X-ray measurements were obtained with an X-ray
log scanning system provided by Finnos Oy. A tree log was scanned jointly with
the wooden calibration phantoms (Figure 2) attached to the saw cut of the one
side of the log. The calibration phantom plane was, therefore, perpendicular
to the log axis of rotation.

As a result of the sequential (“step and shoot”) X-ray scanning, we obtained
projection images capturing the log with attached calibration object from 360
different angles. The width of the projection image is fixed and corresponds
to the number of detector elements $n_{D}$. The height of the projection image
depends on the number of slices that have been scanned (the width of 1 slice
is $1$ px, or equivalently, $2$ mm), that is, the length of the log scanned.
Due to the peculiarities of the measurement system, there is an issue of
missing pixel rows at either the top or bottom of the projection image, and
there is a need to overlay images corresponding to projections of the same log
taken from different angles using an image registration technique Zit2003 . In
the context of our task, we adhere to the following image registration method
consisting of four consecutive steps: feature detection, feature matching,
transform model estimation, image resampling and transformation. The algorithm
was implemented using the OpenCV Python library.

Figure 11: Full-angle reconstructions of log slice. Left: geometry parameters
are estimated using full-angle measurements (360 projection angles); right:
geometry parameters are estimated using sparse measurements (20 projection
angles)

We estimated five unknown geometry parameters (the first scanning angle,
distance between the COR and detector, the source shift, the detector shift,
and the detector tilt) using the reference images of the calibration phantoms
and the corresponding X-ray data obtained from the measurement device. In
Figure 10, we show the cross-sectional reconstructions of one slice of the log
obtained after geometry parameter estimation using the full-angle measurements
of the L-shaped calibration phantom. The figure compares different
reconstruction methods (FBP, Tikhonov regularisation, MAP estimate) for
various numbers of projection angles. In sparse measurements scenario (45 and
20 projection angles), the quality of FBP reconstructions is compromised by
streaking artefacts, while Tikhonov and MAP reconstructions show relatively
better quality.

To show the strength of the calibration method applied to sparse data, we
provide two full-angle reconstructions of the same cross-sectional slice of
the log obtained with different parametrisations (see Figure 11). Geometry
parameters of the left reconstruction were obtained using the full-angle
measurements of the the L-shaped calibration phantom, whereas parameters of
the right reconstruction were estimated from sparse measurements (20
projection angles). Both reconstructions are of good quality and can be used
to locate the tree knots and possible defects in the log.

##  5 Conclusions

In this paper, we have proposed an algorithm for the estimation of unknown
geometry parameters in fan-beam X-ray computed tomography (CT) for log
imaging. The algorithm employs easy-to-produce calibration phantoms that are
scanned jointly with the log and applies differential evolution optimisation
to the objective function represented by the maximum cross-correlation between
the ground truth image of the calibration phantom and its filtered
backprojection (FBP) reconstruction. The algorithm is able to estimate
geometry parameters not only from full-angle measurements of the calibration
phantom but also from sparse measurements.

We demonstrated numerically that there might be multiple parameter vectors
that give almost the same cross-correlation with the reference object when
used in the FBP reconstruction, that is, different sets of optimal parameters
deliver an adequate quality of reconstructions. However, the uniqueness of the
geometry parameters of the measurement equipment does not pose any interest,
since we are only interested in the quality of the reconstructions.

Additionally, we demonstrated the effect of misspecified geometry parameters,
including the first scanning angle, the centre of rotation (COR) to detector
distance, the source and detector shifts, and the detector tilt, on the
quality of CT reconstructions. For reconstruction purposes, we employed
Bayesian inversion (the MAP estimate of the posterior distribution with Cauchy
difference priors) as opposed to classical reconstruction approaches such as
FBP and Tikhonov regularisation in order to improve the reconstruction quality
in the sparse CT setting. Bayesian inversion showed better performance on
sparse data than FBP method that suffered from severe streaking artefacts
compromising the quality of reconstructions.

## Acknowledgments

This work was funded by the Academy of Finland (project numbers 336787 and
334816).

## Data availability

All the codes needed to generate the synthetic data and reconstructions are
available on GitHub:
https://github.com/AngelinaSen/geometry_parameter_estimation. The real log
tomography datasets analysed during the current study are not made public but
are available from the corresponding author on reasonable request.

## References

  *   * (1) Zolotarev, F., Eerola, T., Lensu, L., Kälviäinen, H., Haario, H., Heikkinen, J., Kauppi, T.: Timber tracing with multimodal encoder-decoder networks. In: International Conference on Computer Analysis of Images and Patterns, pp. 342–353 (2019)
  * (2) Flodin, J., Oja, J., Grönlund, A.: Fingerprint traceability of sawn products using x-ray log scanning and sawn timber surface scanning. Forest products journal 58, 100–105 (2008)
  * (3) Tarantola, A.: Inverse Problem Theory and Methods for Model Parameter Estimation. Society for Industrial and Applied Mathematics, Philadelphia (2004)
  * (4) Kaipio, J., Somersalo, E.: Statistical and Computational Inverse Problems. Springer, Dordrecht (2005)
  * (5) Stuart, A.M.: Inverse problems: A Bayesian perspective. Acta Numerica 19, 451–559 (2010). https://doi.org/10.1017/S0962492910000061
  * (6) Natterer, F.: The Mathematics of Computerized Tomography. Wiley, Chicago (1986)
  * (7) Suuronen, J., Emzir, M., Lasanen, S., Särkkä, S., Roininen, L.: Enhancing industrial X-ray tomography by data-centric statistical methods. Data-Centric Engineering 1 (2020)
  * (8) Li, J., Jaszczak, R., Wang, H., Gullberg, G., Greer, K., Coleman, E.: A cone beam SPECT reconstruction algorithm with a displaced center of rotation. Medical physics 21, 145–52 (1994)
  * (9) Wang, H., Smith., M.F., Stone, C.D., Jaszczak, R.J.: Astigmatic single photon emission computed tomography imaging with a displaced center of rotation. Medical physics 25, 1493–1501 (1998)
  * (10) Dennerlein, F., Jerebko, A.: Geometric jitter compensation in cone-beam CT through registration of directly and indirectly filtered projections. In: Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC), pp. 2892–2895 (2012). https://doi.org/10.1109/NSSMIC.2012.6551660
  * (11) Ferrucci, M., Leach, R.K., Giusca, C., Carmignato, S., Dewulf, W.: Towards geometrical calibration of X-ray computed tomography systems — A review 26(9), 092003 (2015). https://doi.org/10.1088/0957-0233/26/9/092003
  * (12) Zemek, M., Šalplachta, J., Zikmund, T., Omote, K., Takeda, Y., Oberta, P., Kaiser, J.: Automatic marker-free estimation methods for the axis of rotation in sub-micron X-ray computed tomography. Tomography of Materials and Structures 1, 100002 (2023). https://doi.org/10.1016/j.tmater.2022.100002
  * (13) Gullberg, G.T., Tsui, B.M.W., Crawford, C.R., Edgerton, E.R.: Estimation of geometrical parameters for fan beam tomography 32(12), 1581–1594 (1987)
  * (14) Cho, Y., Moseley, D., Siewerdsen, J., Jaffray, D.: Accurate technique for complete geometric calibration of cone-beam computed tomography systems. Medical physics 32, 968–83 (2005). https://doi.org/10.1118/1.1869652
  * (15) Ouadah, S., Stayman, J., Gang, G., Ehtiati, T., Siewerdsen, J.: Self-calibration of cone-beam CT geometry using 3D–2D image registration. Physics in Medicine and Biology 61, 2613–2632 (2016). https://doi.org/10.1088/0031-9155/61/7/2613
  * (16) Uribe, F., Bardsley, J.M., Dong, Y., Hansen, P.C., Riis, N.A.B.: A hybrid Gibbs sampler for edge-preserving tomographic reconstruction with uncertain view angles. SIAM/ASA Journal on Uncertainty Quantification 10, 1293–1320 (2021)
  * (17) Riis, N., Dong, Y., Hansen, P.C.: Computed tomography reconstruction with uncertain view angles by iteratively updated model discrepancy. Journal of Mathematical Imaging and Vision 63, 133–143 (2021). https://doi.org/10.1007/s10851-020-00972-7
  * (18) Pedersen, F.H., Jørgensen, J.S., Andersen, M.S.: A bayesian approach to CT reconstruction with uncertain geometry. Applied Mathematics in Science and Engineering 31(1), 2166041 (2023). https://doi.org/10.1080/27690911.2023.2166041
  * (19) Genzel, M., Macdonald, J., März, M.: AAPM DL-Sparse-View CT Challenge submission report: Designing an iterative network for fanbeam-CT with unknown geometry. arXiv (2021). https://arxiv.org/abs/2106.00280
  * (20) Xie, M., Sun, Y., Liu, J., Wohlberg, B., Kamilov, U.S.: Joint reconstruction and calibration using regularization by denoising. arXiv (2020). https://arxiv.org/abs/2011.13391
  * (21) Suuronen, J., Chada, N., Roininen, L.: Cauchy Markov random field priors for Bayesian inversion. Statistics and Computing 32(2) (2022)
  * (22) Radon, J.: On the determination of functions from their integral values along certain manifolds. IEEE Transactions on Medical Imaging 5(4), 170–176 (1986)
  * (23) Kak, A.C., Slaney, M.: Principles of Computerized Tomographic Imaging. Society of Industrial and Applied Mathematics, Philadelphia (2001)
  * (24) Yu, D., Fessler, J.: Edge-preserving tomographic reconstruction with nonlocal regularization. Medical Imaging, IEEE Transactions on 21, 159–173 (2002). https://doi.org/10.1109/42.993134
  * (25) Xu, J., Zhao, Y., Li, H., Zhang, P.: An image reconstruction model regularized by edge-preserving diffusion and smoothing for limited-angle computed tomography. Inverse Problems 35, 085004 (2019). https://doi.org/10.1088/1361-6420/ab08f9
  * (26) Hämäläinen, K., Kallonen, A., Kolehmainen, V., Lassas, M., Niinimäki, K., Siltanen, S.: Sparse tomography. SIAM Journal on Scientific Computing 35, 644–665 (2013). https://doi.org/10.1137/120876277
  * (27) Markkanen, M., Roininen, L., Huttunen, J., Lasanen, S.: Cauchy difference priors for edge-preserving Bayesian inversion. Journal of Inverse and Ill-posed Problems 27 (2019). https://doi.org/10.1515/jiip-2017-0048
  * (28) Storn, R., Price, K.: Differential evolution — A simple and efficient heuristic for global optimization over continuous spaces. Journal of Global Optimization 11, 341–359 (1997)
  * (29) Reddy, M.J., Kumar, D.N.: Multiobjective differential evolution with application to reservoir system optimization. Journal of Computing in Civil Engineering 21, 136–146 (2007). https://doi.org/10.1061/(ASCE)0887-3801(2007)21:2(136)
  * (30) Leon, M., Xiong, N.: Investigation of mutation strategies in differential evolution for solving global optimization problems, vol. 8467 (2014). https://doi.org/%****␣sn-article.bbl␣Line␣475␣****10.1007/978-3-319-07173-2_32
  * (31) Virtanen, P., Gommers, R., Oliphant, T.E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S.J., Brett, M., Wilson, J., Millman, K.J., Mayorov, N., Nelson, A.R.J., Jones, E., Kern, R., Larson, E., Carey, C.J., Polat, İ., Feng, Y., Moore, E.W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E.A., Harris, C.R., Archibald, A.M., Ribeiro, A.H., Pedregosa, F., van Mulbregt, P., SciPy 1.0 Contributors: SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods 17, 261–272 (2020). https://doi.org/10.1038/s41592-019-0686-2
  * (32) Nelder, J.A., Mead, R.: A simplex method for function minimization. Computer Journal 7, 308–313 (1965)
  * (33) Aarts, E., Korst, J.: Simulated Annealing and Boltzmann Machines: A Stochastic Approach to Combinatorial Optimization and Neural Computing. John Wiley & Sons, Inc., New York (1989)
  * (34) Horst, R., Pardalos, P.M.: Handbook of Global Optimization vol. 2. Springer, Dordrecht (2013)
  * (35) Adler, J., Kohr, H., Öktem, O.: Operator Discretization Library (ODL). https://github.com/odlgroup/odl
  * (36) Wang, Z., Bovik, A.C., Sheikh, H., Simoncelli, E.: Image quality assessment: From error visibility to structural similarity. IEEE Transactions on Image Processing 13, 600–612 (2004). https://doi.org/10.1109/TIP.2003.819861
  * (37) Wang, Z., Bovik, A.C.: Mean squared error: Love it or leave it? A new look at signal fidelity measures. IEEE Signal Processing Magazine 26(1), 98–117 (2009). https://doi.org/10.1109/MSP.2008.930649
  * (38) Pambrun, J.-F., Noumeir, R.: Limitations of the SSIM quality metric in the context of diagnostic imaging. (2015). https://doi.org/10.1109/ICIP.2015.7351345
  * (39) Lin, L., Chen, H., Kuruoglu, E., Zhou, W.: Robust structural similarity index measure for images with non-gaussian distortions. Pattern Recognition Letters 163, 10–16 (2022). https://doi.org/10.1016/j.patrec.2022.09.011
  * (40) Zitová, B., Flusser, J.: Image registration methods: A survey. Image and Vision Computing 21(11), 977–1000 (2003)