# Vertex Fitting In Low-Material Budget Pixel Detectors

Andrea Loreti  Department of Physics, University of Liverpool, The Oliver
Lodge Laboratory, Liverpool L69 7ZE, United Kingdom

###### Abstract

This paper provides a detailed description of a vertex fitting algorithm
suitable for precision measurements in low-energy particle physics
experiments. An accurate reconstruction of low-momentum trajectories can be
accomplished by reducing the material budget of the detector to a few per mill
of the radiation length. This limits the multiple scattering undertaken by
particles inside the detector and improves the vertex fitting accuracy.
However, for sufficiently light detection systems, additional sources of
errors, such as the intrinsic spatial resolution of the sensors, must be
considered in the reconstruction of the vertex parameters. The algorithm
developed in this work provides a complete treatment of multiple scattering
and spatial resolution in the context of vertex fitting for light pixel
detectors. In addition to this, a study of the vertex reconstruction in the
low-material budget pixel detector of the Mu3e experiment is presented.

##  I Introduction

With the increase of the instantaneous and integrated beam luminosities, the
requirements of particle physics experiments for precise tracking and
vertexing detectors, with high radiation tolerance, have become more stringent
e.g., SNOEYS2023168678 ; Hartmut2018 ; CARNESECCHI2019608 ; MOSER201685 . In
this regard, silicon pixel sensors can provide high granularity, low material
budget structures and the radiation-hardness that most experiments need, e.g.,
MOSER201685 ; SPANNAGEL2019612 ; Abelev_2014 ; ARNDT2021165679 . It is
important, however, that precise detection systems are developed in
conjunction with equally performing analysis methods for the reconstruction of
particle trajectories and decay vertices, e.g. FRUHWIRTH1987444 ;
BILLOIR1992139 ; Waltenberger2007 ; RevModPhys.82.1419 . To this aim, several
fitting algorithms have been designed and optimized over the years to deal
with hit selection, pattern recognition, errors calculations and high track
multiplicity (see for instance Mankel_2004 ; RevModPhys.82.1419  and
references therein). The practical implementation of these methods must be
tailored around the actual detector and magnetic field configuration of the
experiments. This makes tracking and vertexing a topic which is in continuous
evolution adapting itself to the new challenges set by upcoming experiments.

This study addresses the problem of vertex fitting in the low-material budget
pixel detector of the Mu3e experiment ARNDT2021165679 . As explained in
section II, the detector design has been optimized to minimize Multiple
Coulomb Scattering (MS) effects on particle trajectories and signal
kinematics. However, for light detectors such as Mu3e, the intrinsic pixel
resolution becomes another limiting factor in the vertex reconstruction which
cannot be ignored. Under these circumstances, the vertex fitting should
account for MS and pixel resolution in the error calculations as well as for
any energy loss in the detector which may hamper tracks momentum
reconstruction. The present algorithm deals with all these sources of errors
and it is illustrated in section III whilst in section IV, a comparative study
is made among different inclusive scenarios that encompass: (A) MS only; (B)
MS and pixel resolution together; (C) all sources of errors (MS, pixel
resolution and energy losses) are included.

##  II The low-material budget Mu3e pixel detector

The Mu3e experiment aims to find or exclude the rare Charged Lepton Flavour
(CLF) violating muon decay:

$\mu^{+}\rightarrow e^{+}e^{-}e^{+}$ (1)

at Branching Ratios (BR) $>10^{-16}$ ARNDT2021165679 . This threshold is $4$
orders of magnitude smaller than previous experimental upper limits
($10^{-12}$) BELLGARDT19881  and $38$ orders of magnitude larger than
theoretical Standard Model (SM) calculations (BR $=10^{-54}$), e.g.,
MARCIANO1977303 ; RevModPhys.73.151 . However, new theoretical models predict
the existence of extra degrees of freedom beyond the SM which may bring CLF
violation within the reach of near future experiments such as Mu3e, e.g.,
KAKIZAKI2003210 ; DEGOUVEA201375 . Consequentially, an observation of
$\mu^{+}\rightarrow e^{+}e^{-}e^{+}$ at single event sensitivities aimed by
the Mu3e experiment will imply scenarios of new physics.

The process in (1) yields a relatively simple decay topology with the 3 final
state leptons produced at the same vertex of interaction and momentum vectors,
$\vec{p}$, determined by the energy and momentum conservation for decaying
muons at rest. The main background processes in Mu3e measurements are the muon
internal conversion $\mu^{+}\rightarrow
e^{+}e^{-}e^{+}+\nu_{e}+\bar{\nu}_{\mu}$ (BR $\approx 10^{-5}$) and the
combination of one electron and 2 positrons from independent sources, e.g.,
one Bhabha electron plus two Michel positrons $\mu^{+}\rightarrow
2\times(e^{+}+\nu_{e}+\bar{\nu}_{\mu})+e^{-}$ BELLGARDT19881 .

Figure 1: Scheme of the central barrel of the Mu3e pixel detector, side and
x-y views ARNDT2021165679 .

The aimed single event sensitivity of $2\cdot 10^{-15}$, during phase I of the
experiment, can be achieved with an energy-momentum resolution of $\lesssim 1$
MeV/c and by using precise vertexing and timing systems ARNDT2021165679 .

The energy spectrum of the decay particles in the Mu3e experiment extends up
to $m_{\mu}/2$, where $m_{\mu}$ is the muon mass. In this low-energy region,
MS poses a serious challenge to the reconstruction of particle trajectories
and signal kinematics. To minimize MS, Mu3e uses a low-material budget pixel
detector ($0.1\%$ of the radiation length, Xo, per layer). This is made of
high-voltage monolithic active pixel sensors PERIC2007876  that can be thinned
down to $50\text{\,}\mathrm{\SIUnitSymbolMicro m}$ or $0.05\%$Xo. The rest of
the material budget of the detector is used in the flex-tape that provides
mechanical support and the electrical routing to the sensors.

Figure 1 shows the schematic of the foreseen Mu3e tracker central station
which is important for vertex fitting and track reconstruction ARNDT2021165679
. Two recurl stations, one up-stream and one down-stream, will also be part of
the final detector design. These increase the angular acceptance of the
experiment and allow to measure long-armed trajectories to achieve improved
momentum resolution.

The layers have cylindrical symmetry and are concentrically placed around the
target, a hollow double cone made of Mylar 100 mm in length and with a base
radius of 19 mm. The target is placed in a solenoid magnetic field of 1 T with
the base at a minimum distance of $\approx 4$ mm from the innermost layer of
the pixel tracker. Particle trajectories bend inwards following helical
trajectories around the field lines possibly making multiple re-curls. Each
layer of the pixel detector is sectioned in sub-elements called ladders. A
ladder is a series of chips mounted on the same flex-tape. There are 8, 10, 24
and 28 ladders for layer 1, 2, 3 and 4, respectively. For instance, the
innermost layer of the tracker, crucial for vertex fitting, is made of 8
ladders each one tilted by a 45∘ angle with respect to the neighbours. This
configuration forms a 8-sided surface which extends for $\sim 12$ cm or 6
chips length, see figure 1.

The intrinsic detector spatial resolution is set by the pixel sensitive area,
$80\times 80$$\text{\,}\mathrm{\SIUnitSymbolMicro m}$2. Pixel resolution
becomes more important for high-momentum trajectories for which MS is lower.
However, for low-material budget detectors such as Mu3e, this effect cannot be
ignored at any momentum and it must be treated simultaneously with MS.

##  III Vertex fitting in the Mu3e detector

Vertex reconstruction can be accomplished in two steps: vertex finding and
vertex fitting, see e.g., RevModPhys.82.1419 . The former consists in grouping
trajectories that have been most likely produced in the same decay process.
The latter involves finding the most likely vertex coordinates
$(x,y,z)_{\text{v}}$ and the initial momentum vectors of all clustered tracks.
In Mu3e, vertex finding is accomplished by considering all possible
combinations of two positive and one negative tracks in the detector within
time frames of 64 ns. For the vertex fitting, a least-squares optimization
algorithm has been developed based on the method illustrated in BILLOIR1992139
.

###  III.1 Track parameters and uncertainties

Trajectories are defined by 6 parameters ($[x,y],z,\phi,\lambda,k)$. These are
the coordinates of one point along the track, the angles $\phi$ and $\lambda$
defining the direction of the tangent vector to the trajectory and the factor
$k=(p/q)^{-1}$ where $q$ is the charge of the particle and $p$ is the
magnitude of the momentum vector. The following relationships among the
momentum components in the global Cartesian frame and the angles $\phi$ and
$\lambda$ hold true:

$\left\\{\begin{aligned} p_{x}&=p\,\text{cos}(\lambda)\,\text{cos}(\phi)\,,\\\
p_{y}&=p\,\text{cos}(\lambda)\,\text{sin}(\phi)\,,\\\
p_{z}&=p\,\text{sin}(\lambda)\,,\\\
R_{\perp}&=\frac{p\,\text{cos}(\lambda)}{q\,B}\,,\\\ \phi&\in[0,2\pi]\,,\\\
\lambda&\in[-\pi/2,\pi/2]\,,\end{aligned}\right.$ (2)

where $B=B_{z}$ is the homogeneous magnetic field directed along the beam-line
$\hat{z}$ and $R_{\perp}$ is the transverse radius of the trajectory.

The input measurements to fit are the values of the track parameters at a
given Reference Surface (RS) along with their covariance matrix ($\Xi$):
($[x,y],z,\phi,\lambda,k,\Xi)_{\text{meas}}$. In this study, the RS is the
innermost layer of the pixel detector which is the closest one to the expected
real vertex position. The coordinates $x,y$ are not independent and can be
given in a single expression by using a local reference frame with center in
the middle of the pixel area $(\bar{x},\bar{y},\bar{z})$ and base vectors
$(\hat{u},\hat{z^{\prime}})$. The vector $\hat{z^{\prime}}$ is parallel to the
global coordinate $\hat{z}$ while $\hat{u}$ is perpendicular to it and
parallel to the pixel surface, see figure 2. The equations that link the local
coordinates to the global ones are:

$\begin{array}[]{lll}u=\left(x-\bar{x}\right)\text{cos}(\gamma)+\left(y-\bar{y}\right)\text{sin}(\gamma)\,,\\\
z^{\prime}=z-\bar{z}\,,\\\ \phi^{\prime}=\phi\,,\\\
\lambda^{\prime}=\lambda\\\ k^{\prime}=k\end{array}$ (3)

where the angle $\gamma$ is the orientation of the pixel with respect to the
$x$ axis of the global reference frame. Following the equations in 3, the
track parameters at the RS become $(u,z,\phi,\lambda,k)$ 111The coordinate
$z^{\prime}$ can be replaced by $z$ given that the constant $\bar{z}$ in eq. 3
does not contribute to the calculations of the derivatives and residuals
carried out in section III.2.

Figure 2: Sketch of the RS in the Mu3e vertex fitting. The axes of the global
(local) reference frame are drawn in black (green); the $\hat{z^{\prime}}$
axis perpendicular to the $x-y$ plane is not shown. The parameters
(${x,y},z)_{\text{meas}}$ are placed in the geometrical centre of the pixel
surface $\square$. The fit parameters ($u,z,\phi,\lambda,k)_{\text{fit}}$ are
obtained by forward propagating vertex parameters
($x,y,z,\phi,\lambda,k)_{\text{v}}$ via a function
$h=h(\textbf{v}=(x,y,z)_{\text{v}},\textbf{t}(\phi_{\text{v}},\lambda_{\text{v}}),k_{\text{v}}))$.

The covariance matrix can be written as the sum of two terms, i.e.,
$\Xi_{\text{meas}}=\Xi_{\text{track}}+\Xi_{\text{phys}}$. The term
$\Xi_{\text{track}}$ accounts for the uncertainties accumulated during track
fitting while $\Xi_{\text{phys}}$ accounts for MS and pixel resolution at the
RS. Pixel resolution contributes to the smearing of the hit position by
$\sigma_{l}=l/\sqrt{12}$, where $l$ is the length of the pixel side. The MS
changes the direction of the track upon crossing the RS. The tilt can be
approximated by a Gaussian distribution with zero mean and standard deviation
$\theta_{\text{MS}}$ Yao_2006 :

$\theta_{\text{MS}}=13.6\,\left[\text{MeV/c}\right]\,\frac{q}{p}\sqrt{\frac{x}{\beta^{2}\text{X}_{o}}}\left(1+0.038\,\text{ln}\left(\frac{d}{\text{X}_{o}}\right)\right)\,.$
(4)

In the previous equation $\beta$ is the relativistic velocity and $d$ is the
distance travelled by the particle in the RS. From eq. 4, the changes in
$\phi$ and $\lambda$ due to MS, $\sigma_{\phi}$ and $\sigma_{\lambda}$, can be
obtained by projecting $\theta_{\text{MS}}$ onto the transverse and
longitudinal planes, respectively, see e.g., VALENTAN2009728 :

$\sigma_{\lambda}=\theta_{\text{MS}}\,,\\\ $ (5)

and

$\sigma_{\phi}=\theta_{\text{MS}}/\text{cos}(\lambda)\,.$ (6)

In addition to MS and pixel resolution, an error $\sigma_{k}$ on the track
curvature (or $k$) is introduced in $\Xi_{\text{meas}}$ to account for the
energy lost by electrons and positrons in the detector, see e.g., Mankel_2004
. Although energy losses in a thin silicon layer like the RS are negligible
RevModPhys.60.663 , they are nevertheless included in the present vertex
fitting algorithm to provide a complete treatment of the errors. Without loss
of generality, the covariance matrix in this study is written as
$\Xi_{\text{meas}}$ =
diag$\left\\{\sigma_{u}^{2},\sigma_{z}^{2},\sigma_{\phi}^{2},\sigma_{\lambda}^{2},\sigma_{k}^{2}\right\\}$
with $\sigma_{u}=\sigma_{z}=\sigma_{l}$.

###  III.2 Least-squares algorithm

In the vertex fitting, a map that links the track vertex parameters
($x,y,z,\phi^{i},\lambda^{i},k^{i})_{\text{v}}$ to the parameters at the RS
($u^{i},z^{i},\phi^{i},\lambda^{i},k^{i})_{\text{fit}}$ is needed, as shown in
figure 2. In what follows, $i=1,2,3$ for the ($e^{-},e^{+},e^{+}$) in a Mu3e
decay. Starting from the analytical expression of a helical trajectory of a
charged particle in a magnetic field, this map can be written as:
$h_{j}=h(\textbf{v},\textbf{t}^{i},k^{i}_{\text{v}})_{j}$, where
$h_{j}=(u,z,\phi,\lambda,k)_{\text{fit}}$ for $j=1,2,...,5$,
$\textbf{v}=(x,y,z)_{\text{v}}$ and $\textbf{t}=(\phi,\lambda)_{\text{v}}$,
see section VII for the analytic expressions of $h_{j}$. In first
approximation, this function can be linearized near some initial guessed
vertex parameters ($\textbf{v}_{o},\textbf{t}_{o},k_{o}$):

$h(\textbf{v}_{o}+\delta\textbf{v},\textbf{t}_{i,o}+\delta\textbf{t}_{i},k_{i,o}+\delta
k_{i})_{j}\simeq
h(\textbf{v}_{o},\textbf{t}_{i,o},k_{i,o})_{j}+D_{i}\,\delta\textbf{v}+E_{i}\,\delta\textbf{t}_{i}\,+F_{i}\,\delta
k_{i},$ (7)

where, for each track $i$, the matrices $D_{i}$, $E_{i}$ and $F_{i}$ have
dimensions $5\times 3$, $5\times 2$ and $5\times 1$, respectively, and are
calculated as follows:

$\displaystyle\begin{array}[]{lll}D^{j}_{n}=\frac{\partial
h_{j}}{\partial\textbf{v}_{n}}\hskip
14.22636pt\textbf{v}_{n}=x_{\text{v}},y_{\text{v}},z_{\text{v}}\,\hskip
14.22636pt\text{for }n=1,2,3,\\\ E^{j}_{m}=\frac{\partial
h_{j}}{\partial\textbf{t}_{m}}\hskip
14.22636pt\textbf{t}_{m}=\phi_{\text{v}},\lambda_{\text{v}}\,\hskip
14.22636pt\text{for }m=1,2,\\\ F^{j}=\frac{\partial h_{j}}{\partial
k_{\text{v}}}\end{array}$ (11)

From the definitions in eq. 7, a quadratic cost function, $\chi^{2}$, is
defined:

$\chi^{2}=\sum_{i}\left(\delta
q_{i}-D_{i}\delta\textbf{v}-E_{i}\delta\textbf{t}_{i}-F_{i}\delta
k_{i}\right)^{T}W_{i}\left(\delta
q_{i}-D_{i}\delta\textbf{v}-E_{i}\delta\textbf{t}_{i}-F_{i}\delta
k_{i}\right).$ (12)

where $W=\Xi^{-1}$ is the weight matrix (see for instance BILLOIR1985115 ) and
$\delta q$ is the residual at the RS:

$\displaystyle\delta q^{i}$
$\displaystyle=(u^{i},z^{i},\phi^{i},\lambda^{i},k^{i})_{\text{meas}}-h(\textbf{v}_{o},\textbf{t}^{i}_{o},k_{o}^{i})\,.$
(13)

Equation 12 is the total normalized error, due to the approximation in eq. 7,
expressed as a function of $\delta\textbf{v}$, $\delta\textbf{t}$ and $\delta
k$. The accuracy of such an approximation can be seen in figure 3(a) where the
RMS of the deviations $(x,y,z)_{\text{meas}}-(x,y,z)_{\text{fit}}$ at the RS
are plotted versus fit iteration number.

The corrections to the initial guessed vertex parameters can be found by
minimizing the cost function, i.e., by solving the following system of
equations:

$\left\\{\begin{aligned}
&\sum_{i}A_{i}^{T}\delta\textbf{v}+\sum_{i}B^{T}_{i}\delta\textbf{t}_{i}+\sum_{i}H^{T}_{i}\delta
k_{i}=\sum_{i}T_{i},\\\
&B_{i}\delta\textbf{v}+C_{i}\delta\textbf{t}_{i}+G_{i}\delta k_{i}=U_{i},\\\
&H_{i}\delta\textbf{v}+G_{i}^{T}\delta\textbf{t}_{i}+L_{i}\delta
k_{i}=Z_{i}.\end{aligned}\right.$ (14)

where

$\displaystyle A_{i}^{T}=D_{i}^{T}W_{i}D_{i}\,,\hskip
14.22636ptB_{i}=E_{i}^{T}W_{i}D_{i}\,\hskip
14.22636ptC_{i}=E_{i}^{T}W_{i}E_{i}\,,$ (15) $\displaystyle
G_{i}=E_{i}^{T}W_{i}F_{i}\,,\hskip 14.22636ptH_{i}=F^{T}W_{i}D_{i}\,,\hskip
14.22636ptL_{i}=F_{i}^{T}W_{i}F_{i}\,,$ $\displaystyle
T_{i}=D_{i}^{T}W_{i}\delta q_{i}\,,\hskip 14.22636ptU_{i}=E_{i}^{T}W_{i}\delta
q_{i}\,\hskip 14.22636ptZ_{i}=F_{i}^{T}W_{i}\delta q_{i}.$

The solutions of system described in eq. 14 are:

$\displaystyle\delta\textbf{v}$
$\displaystyle=\left[\left(\sum_{i}A^{T}-\sum_{i}B^{T}_{i}C^{-1}_{i}B_{i}\right)+\left(\sum_{i}B^{T}_{i}C^{-1}_{i}G_{i}-\sum
H_{i}^{T}\right)N_{1}N_{3}\right]^{-1}$ (16)
$\displaystyle\times\left(\sum_{i}T_{i}-\sum_{i}B_{i}^{T}C^{-1}_{i}U_{i}+\left(\sum_{i}B^{T}_{i}C^{-1}_{i}G_{i}-\sum
H_{i}^{T}\right)N_{1}N_{2}\right),$ $\displaystyle\delta k_{i}$
$\displaystyle=N_{1}\left(N_{2}-N_{3}\delta\textbf{v}\right),$
$\displaystyle\delta\textbf{t}_{i}$
$\displaystyle=C_{i}^{-1}\left(U_{i}-B_{i}\delta\textbf{v}-G_{i}\delta
k_{i}\right),$

where:

$\displaystyle N_{1}$
$\displaystyle=\left(L_{i}-G^{T}_{i}C^{-1}_{i}G_{i}\right)^{-1}\,,$ (17)
$\displaystyle N_{2}$
$\displaystyle=\left(Z_{i}-G^{T}_{i}C^{-1}_{i}U_{i}\right)\,$ $\displaystyle
N_{3}$ $\displaystyle=\left(H_{i}-G^{T}_{i}C^{-1}_{i}B_{i}\right)\,.$

From eq. 16 the covariance matrices for the track parameters at the vertex can
be calculated:

$\displaystyle\text{Cov}(\textbf{v})$
$\displaystyle=\left[\left(\sum_{i}A^{T}-\sum_{i}B^{T}_{i}C^{-1}_{i}B_{i}\right)+\left(\sum_{i}B^{T}_{i}C^{-1}_{i}G_{i}-\sum
H_{i}^{T}\right)N_{1}N_{3}\right]^{-1},$ (18) $\displaystyle\text{Cov}(k_{i})$
$\displaystyle=L_{i}^{-1}+\left(N_{1}N_{3}\right)\text{Cov}(\textbf{v})\,\left(N_{1}N_{3}\right)^{T},$
$\displaystyle\text{Cov}(\textbf{t}_{i})$
$\displaystyle=C_{i}^{-1}+\left(C_{i}^{-1}B_{i}\right)\text{Cov}(\textbf{v})\left(C_{i}^{-1}B_{i}\right)^{T}+\left(C_{i}^{-1}G_{i}\right)\text{Cov}(k)\left(C_{i}^{-1}G_{i}\right)^{T}.$

###  III.3 Algorithm testing

Geant4 based Monte-Carlo (MC) simulations of Mu3e decays have been carried out
for testing the vertex fitting described in section III.2. The procedure
followed by the test was:

  1. 1.

Hit parameters at the RS, $(u_{i},z_{i},\phi_{i},\lambda_{i},k_{i})_{MC}$,
were obtained from MC trajectories.

  2. 2.

$(u_{i},z_{i})_{MC}$ were smeared by using pixel resolution, i.e., by adding a
random offset derived from a uniform distribution with mean and standard
deviation $(0,\sigma_{l})$. The $(\phi_{i},\lambda_{i})$ angles were tilted
according to the errors in equations 5 and 6 for values of the MS pooled from
a Gaussian distribution with $\hat{\mu}=0$ and $\sigma=\theta_{\text{MS}}$.
The error on $k_{i}$ was obtained by drawing from a Gaussian distribution with
mean zero and a standard deviation $\sigma_{k}$.

  3. 3.

From point 1 and 2, $\Xi_{\text{meas}}$ was derived together with the weight
matrix $W$.

  4. 4.

The initial parameters $\textbf{v}_{o}=x_{o},y_{o},z_{o}$ were obtained as the
average coordinates of the tracks intersection points. Since two tracks can
have up to two intersections, the one that is met first when back propagating
the track parameters from the RS to the target is retained in the calculation
of the average value. If two tracks do not intersect, the point of closest
approach is used. The vectors $\textbf{t}_{o}$ were extracted at the point of
closest approach of the trajectories to $\textbf{v}_{o}$ whilst $k_{o}$ was
directly obtained from the track reconstruction carried out before the vertex
fitting.

  5. 5.

From equations 7 and 11, the residuals at the RS were calculated and thus the
corrections $\delta\textbf{v}$, $\delta k_{i}$ and $\delta\textbf{t}_{i}$ in
eq. 16. Only a few iterations were required to minimize the $\chi^{2}$ in eq.
12, as it can be seen in figure 3(b).

(a)  (b)

Figure 3: Typical RMS of the residuals $(x,y,z)_{meas}-(x,y,z)_{fit}$ at the
RS (a) and the $\chi^{2}$ minimization (b) as a function of the fit iteration
number.

(a)  (b)

(c)  (d)

(e)  (f)

Figure 4: Pull distributions (a-f) of the track parameters at the vertex:
($x,y,z,\phi,\lambda,k)_{\text{v}}$ for test trajectories obtained by
following the steps described in the text.

A precise determination of the fit errors depends on the correct
characterization of $\Xi_{\text{meas}}$ and its propagation, e.g.,
WOLIN1993493 Lund_2009 . In the present study, the covariance matrix of the
vertex parameters was calculated by propagating $\Xi_{\text{meas}}$ from the
RS to the vertex point by using equations 18. A key test of the algorithm
developed in this section consists in plotting the pull distributions of the
track parameters. The pull of a variable $X$ with expected value $\mu_{X}$ and
standard error $\sigma_{X}$ is:

$P=\frac{\left(X-\mu_{X}\right)}{\sigma_{X}}.$ (19)

If the error and the residual in eq. 19 are well characterized, the pulls are
normally distributed.

Figure 4 shows the normal distributions for the pulls of all the track
parameters at the vertex $(x,y,z,\phi,\lambda,k)_{\text{v}}$ as calculated by
the present vertex fitting algorithm.

##  IV Comaprative study of error sources

In section III.1, the explicit form of the covariance matrix was given by
including the contribution of MS, pixel resolution and energy losses at RS:
$\Xi_{\text{meas}}=\text{diag}\left\\{\sigma_{u}^{2},\sigma_{z}^{2},\sigma_{\phi}^{2},\sigma_{\lambda}^{2},\sigma_{k}^{2}\right\\}$.
In this section, the relative weights of these errors on the determinations of
the fit vertex parameters and their uncertainties are discussed. Three
different inclusive scenarios have been considered: (A) MS only; (B) MS and
pixel resolution are both included; (C) all sources of errors (MS, pixel
resolution and energy losses) are considered. These three scenarios are
summarized in table 1.

It must be noted that the kernel of $\Xi_{\text{meas}}$ grows larger going
from scenario (A) to (C) along with the dimensionality of the problem, see for
example (BILLOIR1992139, ). For instance, if the energy loss and pixel
resolution at the RS are ignored, MS remains the only source of uncertainty
against which all measurements are fixed but
$(\phi_{i},\lambda_{i})_{\text{meas}}$, i.e.,
$\Xi_{\text{meas}}=\text{diag}\left\\{\sigma_{\phi}^{2},\sigma_{\lambda}^{2}\right\\}$.
These 6 angles can be fitted with 3 vertex variables $(x,y,z)_{\text{v}}$ thus
simplifying eq. 16, see e.g., schenk2013vertex .

Table 1: Scenrarios A,B,C in the comparative study of the vertex fitting. Scenario | $\Xi_{\text{meas}}$ | measurements | fit parameters | errors
---|---|---|---|---
A | $\text{diag}\left\\{\sigma_{\phi}^{2},\sigma_{\lambda}^{2}\right\\}$ | $(\phi,\lambda)_{\text{meas}}$ | $(x,y,z)_{\text{v}}$ | MS
B | $\text{diag}\left\\{\sigma_{u}^{2},\sigma_{z}^{2},\sigma_{\phi}^{2},\sigma_{\lambda}^{2}\right\\}$ | $(u,z,\phi,\lambda)_{\text{meas}}$ | $(x,y,z,\phi,\lambda)_{\text{v}}$ | MS, pixel
C | $\text{diag}\left\\{\sigma_{u}^{2},\sigma_{z}^{2},\sigma_{\phi}^{2},\sigma_{\lambda}^{2},\sigma_{k}^{2}\right\\}$ | $(u,z,\phi,\lambda,k)_{\text{meas}}$ | $(x,y,z,\phi,\lambda,k)_{\text{v}}$ |  MS, pixel, $\Delta E$

Panels (a-e) in figure 5 show the the deviations between MC vertex parameters
of simulated Mu3e decays and those obtained from vertex fitting in scenarios
(A) and (B), respectively. From figure 5(a-c), it can be seen that the fit
accuracy for the determinations of $(x,y,z)_{\text{v}}$ does not improve when
the pixel resolution is included in covariance matrix. However, a significant
improvement is found on the determinations of
$(\phi_{i},\lambda_{i})_{\text{v}}$, as shown in figure 5(d,e). The results
for scenario (C) are statistically the same as for scenario (B). The former
being the only case in which the fit attempts to optimize the track parameter
$k$, see table 1. As expected, having neglected the energy losses at the RS,
track curvatures do not vary significantly throughout the vertex fitting, as
it can be seen in figure 5(f).

The improved fit accuracy of case (C) with respect to scenario (A) for the
coordinates of the momentum vector $p_{x}$, $p_{y}$ and $p_{z}$ is shown in
figure 6(a-c). This improvement reflects also onto the determination of the
average total momentum which is $\sim$10$\%$ smaller in scenario (C) than the
corresponding average in scenario (A), and thus closer to real MC value (in
the hypothesis of decaying muons at rest), see figure 6(d).

The invariant mass of simulated Mu3e decays, calculated from the fit vertex
parameters in scenarios (A) and (C), is shown in figure 7. As expected, no
significant difference is seen in the two fit scenarios. In fact, the
magnitude of the invariant mass is dominated by the muon rest mass over which
the fit accuracy has little leverage.

(a)  (b)

(c)  (d)

(e)  (f)

Figure 5: Deviations between the fit and MC track parameters at the vertex.
The legends show the mean and standard deviation obtained from a Gaussian fit
in scenarios (A) [full gray] and (B or C) [empty red], respectively.

(a)  (b)

(c)  (d)

Figure 6: Deviations between the fit and MC momentum coordinates $p_{x}$,
$p_{y}$ and $p_{z}$ panels (a-c) and total momenutm $P$ in panel (d) for
scenario (A) [full gray] and (C) [empty red], respectively. In (a,b,c), the
legends show the mean and standard deviation obtained from a Gaussian fit
whilst the legend in (d) shows the average and stadard deviation of the
distributions. Figure 7: Rreconstructed invariant mass for vertices with
$\chi^{2}<15$ and total momentum $P<4$ MeV/c, scenario (A) [full gray] and (C)
[empty red], respectively.

##  V Conclusion

In this paper, a simple least-squares method has been described which can be
applied to reconstruct decay vertices in experiments equipped with pixel
detectors. The relative weights of 1) MS, 2) pixel resolution and 3) energy
losses to the final reconstruction accuracy has been investigated in the case
study of the Mu3e low-material budget pixel detector. The exhaustive errors
treatment of the present study goes beyond the MS-only approximation showing a
significant improvement of the fit accuracy when the intrinsic pixel
resolution is accounted for. This should encourage a rigorous treatment of the
pixel resolution in the development of future reconstruction algorithms
concerning precise particle physics measurements at low-energy.

##  VI Acknowledgements

I am grateful to the STFC grant for supporting this work. I wish to thank the
members of the Mu3e Software and Analysis group for providing the simulation
and track reconstruction software behind this study. I also want to thank Joel
Goldstein, Niklaus Berger and Gavin Hesketh for all the detailed and useful
discussions about this work. I also thank Naik Paras and Andre Schoning for
their careful reading and useful comments.

##  VII Appendix 1: forward propagation of track parameters

The map $h(\textbf{v},\textbf{t},k)$ propagates a trajectory from the vertex
to the hit at the RS. In this section, its analytical expression is derived by
starting from the propagation of a track in the transverse plane and then
along the beam direction.

#### Propagation in the transverse plane

In this study, trajectories are helices with symmetry axis $\hat{z}$ and
transverse radius $R_{\perp}$ which sign is given by the charge $q$ of the
particle, see eq. 2. We write $R_{\perp}=R$ cos($\lambda$) such as:

$R:=\left\\{\begin{aligned} &+\frac{p}{|q|B}\hskip 8.5359pt\text{if q
}>0\text{ c.c.w rotation}\,\\\ &-\frac{p}{|q|B}\hskip 8.5359pt\text{if q
}<0\text{ c.w rotation}\,.\end{aligned}\right.$ (20)

In the x-y plane, a helix is a circumference with center in $(x_{c},y_{c})$
and radius $R_{\perp}$. It is not difficult to prove that:

$\displaystyle x_{c}=x_{\text{v}}-R_{\perp}\text{sin}(\phi_{\text{v}})\,,$
(21) $\displaystyle
y_{c}=y_{\text{v}}+R_{\perp}\text{cos}(\phi_{\text{v}})\,,$

Figure 8: Sketch of a track with negative radius, as defined in eq. 20,
projected on the x-y plane. The $\phi_{\text{v}}$ values are obtained by
subtracting $\pi/2$ from the phase angle of the radial vector.

The transport equations in the transverse plane is obtained by calculating the
coordinates $(x_{q},y_{q})$ of the intersection point between the track
originating from $(x_{\text{v}},y_{\text{v}})$ and the detector ladder. In the
x-y plane, the ladder profile is a line characterised by the parameters
$y_{o}$ and $m$, see figure 9:

$\left\\{\begin{aligned}
&(y_{q}-y_{c}(\textbf{v},\textbf{t},k))^{2}+(x_{q}-x_{c}(\textbf{v},\textbf{t},k))^{2}=R_{\perp}^{2}\,,\\\
&y_{q}=m\,x_{q}+y_{o}\,,\\\ &m=\text{tg}(\gamma)\,.\end{aligned}\right.$ (22)

Figure 9: A sketch representing the intersection between a trajectory and the
detector RS in the x-y plane.

In the previous equation, $\gamma$ is the angle of the detector ladder with
respect to the global $\hat{x}$ axis. The parameter $y_{o}$ can be calculated
by using $y_{\text{meas}}=tg(\gamma)x_{\text{meas}}+y_{o}$. In conclusion, the
solutions of the system of equations 22 can be written as:

$\begin{array}[]{l}x_{q}=\left(\begin{array}[]{c}-\frac{\mathrm{y_{o}}-\frac{\mathrm{y_{o}}+m\,\mathrm{x_{v}}+\sigma_{1}+m^{2}\,\mathrm{y_{v}}+\sigma_{2}-\sigma_{3}}{m^{2}+1}}{m}\\\
-\frac{\mathrm{y_{o}}-\frac{\mathrm{y_{o}}+m\,\mathrm{x_{v}}-\sigma_{1}+m^{2}\,\mathrm{y_{v}}+\sigma_{2}-\sigma_{3}}{m^{2}+1}}{m}\end{array}\right)\\\
\end{array}$ (23)

and

$y_{q}=\left(\begin{array}[]{c}\frac{\mathrm{y_{o}}+m\,\mathrm{x_{v}}+\sigma_{1}+m^{2}\,\mathrm{y_{v}}+\sigma_{2}-\sigma_{3}}{m^{2}+1}\\\
\frac{\mathrm{y_{o}}+m\,\mathrm{x_{v}}-\sigma_{1}+m^{2}\,\mathrm{y_{v}}+\sigma_{2}-\sigma_{3}}{m^{2}+1}\end{array}\right)\,,\\\
\mathrm{}\\\ $ (24)

where

$\begin{array}[]{l}\mathrm{}\\\
\;\;\sigma_{1}=m\,\left[-R^{2}\,m^{2}\,{\mathrm{cos}\left(\mathrm{\lambda_{v}}\right)}^{2}\,{\mathrm{sin}\left(\mathrm{\phi_{v}}\right)}^{2}+R^{2}\,m^{2}\,{\mathrm{cos}\left(\mathrm{\lambda_{v}}\right)}^{2}\right.\\\
\left.-2\,R^{2}\,m\,{\mathrm{cos}\left(\mathrm{\lambda_{v}}\right)}^{2}\,\mathrm{cos}\left(\mathrm{\phi_{v}}\right)\,\mathrm{sin}\left(\mathrm{\phi_{v}}\right)-R^{2}\,{\mathrm{cos}\left(\mathrm{\lambda_{v}}\right)}^{2}\,{\mathrm{cos}\left(\mathrm{\phi_{v}}\right)}^{2}\right.\\\
\left.+R^{2}\,{\mathrm{cos}\left(\mathrm{\lambda_{v}}\right)}^{2}+2\,R\,m^{2}\,\mathrm{x_{v}}\,\mathrm{cos}\left(\mathrm{\lambda_{v}}\right)\,\mathrm{sin}\left(\mathrm{\phi_{v}}\right)+2\,R\,m\,\mathrm{x_{v}}\,\mathrm{cos}\left(\mathrm{\lambda_{v}}\right)\,\mathrm{cos}\left(\mathrm{\phi_{v}}\right)\right.\\\
\left.+2\,R\,m\,\mathrm{y_{o}}\,\mathrm{cos}\left(\mathrm{\lambda_{v}}\right)\,\mathrm{sin}\left(\mathrm{\phi_{v}}\right)-2\,R\,m\,\mathrm{y_{v}}\,\mathrm{cos}\left(\mathrm{\lambda_{v}}\right)\,\mathrm{sin}\left(\mathrm{\phi_{v}}\right)\right.\\\
\left.+2\,R\,\mathrm{y_{o}}\,\mathrm{cos}\left(\mathrm{\lambda_{v}}\right)\,\mathrm{cos}\left(\mathrm{\phi_{v}}\right)-2\,R\,\mathrm{y_{v}}\,\mathrm{cos}\left(\mathrm{\lambda_{v}}\right)\,\mathrm{cos}\left(\mathrm{\phi_{v}}\right)-m^{2}\,{\mathrm{x_{v}}}^{2}\right.\\\
\left.-2\,m\,\mathrm{x_{v}}\,\mathrm{y_{o}}+2\,m\,\mathrm{x_{v}}\,\mathrm{y_{v}}-{\mathrm{y_{o}}}^{2}+2\,\mathrm{y_{o}}\,\mathrm{y_{v}}-{\mathrm{y_{v}}}^{2}\right]^{1/2}\\\
\mathrm{}\\\
\;\;\sigma_{2}=R\,m^{2}\,\mathrm{cos}\left(\mathrm{\lambda_{v}}\right)\,\mathrm{cos}\left(\mathrm{\phi_{v}}\right)\\\
\mathrm{}\\\
\;\;\sigma_{3}=R\,m\,\mathrm{cos}\left(\mathrm{\lambda_{v}}\right)\,\mathrm{sin}\left(\mathrm{\phi_{v}}\right)\,.\end{array}$
(25)

A choice between the two solutions in equations 23 and 24 can be made by
accounting for the track direction of motion and the vertex position relative
to the hit.

For what is concerned with $\phi_{\text{v}}$, the following expression can be
written, see figure 8:

$\phi_{q}=\text{atan}\left(\frac{y_{q}-y_{c}(x_{\text{v}},y_{\text{v}},\phi_{\text{v}},\lambda_{\text{v}},k_{\text{v}})}{x_{q}-x_{c}(x_{\text{v}},y_{\text{v}},\phi_{\text{v}},\lambda_{\text{v}},k_{\text{v}})}\right)+\text{sign}(R)\pi/2\,.$
(26)

#### Propagation along the beam axis

The propagation of helical trajectories along the $\hat{z}$ axis is
characterized by the following equations VALENTAN2009728 :

$\left\\{\begin{aligned} &\lambda_{q}=\lambda_{\text{v}}\,,\\\
&z_{q}=z_{\text{v}}+R_{\perp}\,\text{tan}(\lambda_{\text{v}})\left(\phi_{q}-\phi_{\text{v}}\right)\,.\end{aligned}\right.$
(27)

## References

## References

  * [1] W. Snoeys.  Monolithic CMOS sensors for high energy physics — challenges and perspectives.  NIM-A, 1056:168678, 2023.
  * [2] H. F. W. Sadrozinski, A. Seiden, and N. Cartiglia.  4D tracking with ultra-fast silicon detectors.  Rep. Prog. Phys., 81:026101, 2018.
  * [3] F. Carnesecchi et al.  Development of ultra fast silicon detector for 4D tracking.  NIM-A, 936:608–611, 2019.  Frontier Detectors for Frontier Physics: 14th Pisa Meeting on Advanced Detectors.
  * [4] H. G. Moser.  The Belle II DEPFET pixel detector.  NIM-A, 831:85–87, 2016.  Proceedings of the 10th International “Hiroshima” Symposium on the Development and Application of Semiconductor Tracking Detectors.
  * [5] S. Spannagel.  Technologies for future vertex and tracking detectors at CLIC.  NIM-A, 936:612–615, 2019.  Frontier Detectors for Frontier Physics: 14th Pisa Meeting on Advanced Detectors.
  * [6] B. Abelev and The ALICE Collaboration.  Technical design report for the upgrade of the ALICE inner tracking system.  Journal of Physics G: Nuclear and Particle Physics, 41(8):087002, 2014.
  * [7] K. Arndt et al.  Technical design of the phase I Mu3e experiment.  NIM-A, 1014:165679, 2021.
  * [8] R. Frühwirth.  Application of Kalman filtering to track and vertex fitting.  NIM-A, 262(2):444–450, 1987.
  * [9] P. Billoir and S. Qian.  Fast vertex fitting with a local parametrization of tracks.  NIM-A, 311(1):139–150, 1992.
  * [10] W. Waltenberger, R. Frühwirth, and P. Vanlaer.  Adaptive vertex fitting.  Journal of Physics G: Nuclear and Particle Physics,, 34:N343, 2007\.
  * [11] A. Strandlie and R. Frühwirth.  Track and vertex reconstruction: from classical to adaptive methods.  Rev. Mod. Phys., 82:1419–1458, 2010.
  * [12] R. Mankel.  Pattern recognition and event reconstruction in particle physics experiments.  Reports on Progress in Physics, 67(4):553, 2004.
  * [13] U. Bellgardt et al.  Search for the decay ${\mu}^{+}$ ${\rightarrow}$ ${e}^{+}$${e}^{+}$${e}^{-}$.  Nuclear Physics B, 299(1):1–6, 1988.
  * [14] W.J. Marciano and A.I. Sanda.  Exotic decays of the muon and heavy leptons in gauge theories.  Physics Letters B, 67(3):303–305, 1977.
  * [15] Y. Kuno and O. Yasuhiro.  Muon decay and physics beyond the standard model.  Rev. Mod. Phys., 73:151–202, 2001.
  * [16] M. Kakizaki, Y. Ogura, and F. Shima.  Lepton flavor violation in the triplet Higgs model.  Physics Letters B, 566(3):210–216, 2003.
  * [17] A. de Gouvêa and P. Vogel.  Lepton flavor and number conservation, and physics beyond the standard model.  Progress in Particle and Nuclear Physics, 71:75–92, 2013.  Fundamental Symmetries in the Era of the LHC.
  * [18] I. Perić.  A novel monolithic pixelated particle detector implemented in high-voltage CMOS technology.  NIM-A, 582(3):876–885, 2007.  VERTEX 2006.
  * [19] W. M. Yao et al.  Review of particle physics.  Journal of Physics G: Nuclear and Particle Physics, 33(1):1, 2006\.
  * [20] M. Valentan, M. Regler, and R. Frühwirth.  Generalization of the Gluckstern formulas II: multiple scattering and non-zero dip angles.  NIM-A, 606(3):728–742, 2009.
  * [21] H. Bichsel.  Straggling in thin silicon detectors.  Rev. Mod. Phys., 60:663–699, 1988.
  * [22] P. Billoir, R. Frühwirth, and M. Regler.  Track element merging strategy and vertex fitting in complex modular detectors.  NIM-A, 241(1):115–131, 1985.
  * [23] E. J. Wolin and L. L. Ho.  Covariance matrices for track fitting with the Kalman filter.  NIM-A, 329(3):493–500, 1993.
  * [24] E. Lund et al.  Transport of covariance matrices in the inhomogeneous magnetic field of the atlas experiment by the application of a semi-analytical method.  Journal of Instrumentation, 4(04):P04016, 2009.
  * [25] Sebastian Schenk.  A vertex fit for low momentum particles in a solenoidal magnetic field with multiple scattering.  Master’s Thesis, Heidelberg University, 2013.


A Deep Learning Generative Model Approach
for Image Synthesis of Plant Leaves

Alessandro Benfenati 1, Davide Bolzi2, Paola Causin 2*, Roberto Oberti 3

1 Dept. of Environmental Science and Policy, Università degli Studi di Milano,
Milano, Italy

2 Dept. of Mathematics, Università degli Studi di Milano, Milano, Italy

3 Dept. of Agricultural and Environmental Sciences - Production, Landscape,
Agroenergy, Università degli Studi di Milano, Milano, Italy

*<EMAIL_ADDRESS>
## Abstract

### Objectives

We generate via advanced Deep Learning (DL) techniques artificial leaf images
in an automatized way. Our aim is to dispose of a source of training samples
in artificial intelligence applications for modern crop management in
agriculture, with focus on disease recognition on plant leaves. Such
applications require large amounts of data and, while leaf images are not
truly scarce, image collection and annotation remains a very time–consuming
process. Data scarcity can be addressed by augmentation techniques consisting
in simple transformations of samples belonging to a small dataset, but the
richness of the augmented data is limited: this motivates the search for
alternative approaches.

### Methods

Pursuing an approach based on DL generative models, we propose a Leaf-to-Leaf
Translation (L2L) procedure structured in two steps: firstly, a residual
variational autoencoder architecture is used to generate synthetic leaf
skeletons (leaf profile and veins) starting from companions binarized
skeletons of real leaf images. In a second step, we perform the process of
translation via a Pix2pix framework, which uses conditional generator
adversarial networks to reproduce the colorization of leaf blades, preserving
the shape and the venation pattern.

### Results

The L2L procedure generates synthetic images of leaves with a realistic
appearance. We address the performance measurement both in a qualitative and a
quantitative way; for this latter evaluation, we employ a DL anomaly detection
strategy which quantifies the degree of anomaly of synthetic leaves with
respect to real samples.

### Conclusions

Generative DL approaches have the potential to be a new paradigm to provide
low-cost meaningful synthetic samples for computer-aided applications. The
present L2L approach represents a step towards this goal, being able to
generate synthetic samples with a relevant qualitative and quantitative
resemblance to real leaves.

## Author summary

In this work we present an end-to-end workflow incorporating state-of-the-art
Deep Learning strategies based on generative methods to produce realistic
synthetic images of leaves. At the best of our knowledge, this is the first
attempt of such an approach to this problem. Our focus application is the
training of neural networks for modern crop management systems in agriculture,
but we believe that many other computer–aided applications may benefit from
it. We take inspiration from previous works carried out on eye retina image
synthesis, an application domain which shares some similarities with the
present problem (a venation pattern over a colorized “fundus”). Our approach
relies on the successive use of autoencoders and generative adversarial
architectures, able to generate leaf images both in the Red-Green-Blue
channels as well as in the Near-Infra-Red. The generated leaf images have a
realistic appearance even if they sometimes suffer from small inconsistencies,
especially discolored patches. A quantitative evaluation via an anomaly
detection algorithm shows that in average a synthetic sample is classified as
such only in 25% of the cases.

## Introduction

The ability to generate meaningful synthetic images of leaves is highly
desirable for many computer-aided applications. At the best of our knowledge,
attempts at generating synthetic images of leaves have been made mostly in the
field of computer graphics and were aimed at creating the illusion of
realistic landscapes covered with plants, trees or meadows. These efforts were
mainly based on procedural mathematical models describing the venation
structure and the color/texture of the leaf. A specific type of formal
grammar, called L–grammar, was developed to generate instructions to draw a
leaf. Several instances of the profile of leaves of a certain species were
created upon random variations of parameters of the directives of a certain
L–grammar [1]. Biologically-motivated models were also proposed. A main point
of these approaches is the representation of the interaction between auxin
sources distributed over the leaf blade and the formation of vein patterns
[2]. Some attempts were also carried out using finite elements to build
mechanistic models of the leaf blade, tuned on its structural parameters [3].
After generating the leaf shape and venation pattern - regardless of the
approach- texture and colors were rendered by a color palette prescribed by
the user or generated according to a pseudo–random algorithm. A color model
based on convolution sums of divisor functions was proposed in [4], while a
shading model based on the PROSPECT model for light transmission in leaves
[5], was proposed in [6]. “Virtual rice” leaves were created in [7] based on a
RGB-SPAD model.

In this work we aim at introducing a radically different approach by
generating artificial images of leaves by automatized techniques based on Deep
Learning (DL) techniques. Our focus is mainly to enrich dataset of leaf images
for neural networks training, even if we deem that the present approach may be
of interest also in a wide range of other fields, starting from computer
graphics. The motivation underlying this work is that DL methods require a
large amount of data – often of the order of hundreds of thousands of images –
to avoid overfitting phenomena. Data augmentation is a common remedy, usually
consisting in simple transformations such as random rotations, translations or
deformations of the original images. However, the richness of the augmented
dataset is limited and more sophisticated approaches for synthesizing
additional training data have a greater potential to improve the training
process. In this respect, DL generative models represent attractive methods to
produce synthetic images (with corresponding labels) using the information
from a limited set of real, unlabeled images of the same domain. This idea is
not new - in absolute - but it has been used mainly in the field of medicine,
where data may be extremely scarce and difficult to obtain (see, e.g., the
recent review [8]). In the present context, while scarcity of leaf images may
be not a real issue, what is more relevant is to avoid the huge mole of work
required to collect, examine and annotate images. This is especially true when
image segmentation should be performed, which is a pixel-wise problem: the
acquisition of annotated segmentation masks is exceedingly costly and time
consuming, as a human expert annotator has to label every pixel manually. For
our model we take inspiration from [9] (and reference therein), where the
authors synthesised eye retina images. The fundus of the eye shares indeed
several characteristics with our problem: a fine network of hierarchically
organized blood vessels (leaf veins) superposed to a colored background (leaf
blade). In addition, in our problem the leaf blade is also characterized by a
specific silhouette that must be represented as well. We propose a Leaf-to-
Leaf Translation (L2L) approach to obtain synthetic colorized leaf blades
organized in two steps: first we use a residual variational autoencoder
architecture to generate fake leaf skeletons starting from binarized companion
skeletons of real leaf images. In a second step we perform the process of
translation via a Pix2pix framework, which uses conditional generator
adversarial networks (cGANs) to reproduce the specific color distribution of
the leaf blade, preserving leaf shape and venation pattern. We carry out both
qualitative and quantitative evaluations of the degree of realism of the
synthetic samples of leaves. Specifically, a DL-based anomaly detection
strategy is used to evaluate the distance (“anomaly”) between synthetic and
real samples. The results show a good degree of realism, that is a low anomaly
score, and indicate that with the present approach one can significantly
enrich a small dataset and improve the training performance of DL
architectures.

## Materials and methods

### Dataset

Grapevine leaves were imaged via a QSi640 ws-Multispectral camera (Atik
Cameras, UK) equipped with a Kodak 4.2 Mp micro-lens image sensor and 8
spectral selection filters operating in the bands 430 to 740 nm. For the
purpose of this experiment, leaves were imaged singularly on a dark
background, under controlled diffuse illumination conditions. Images were
acquired in the single spectral channels 430 nm (blue, B), 530 nm (green, G),
685 nm (red, R) and 740 nm (near–infrared, NIR). These channels are typically
considered when dealing with the task of recognition of plant diseases in a
multispectral analysis approach [10, 11]. A set of RGB images of the same
leaves in standard CIE color space were also acquired for reference. Camera
parameters were set and image collection was performed via an in–house
developed acquisition software written in MATLAB. Reflectance calibration was
carried out by including in each image 3 reflectance references (Spectralon R
= 0.02, R = 0.50 and R = 0.99; Labsphere, USA). We obtained photos of 80
leaves with a resolution of $2048\times 2048$ pixels and 8 bit for each
channel. Preprocessing operations were performed on each image: removal of hot
pixels, normalization along each channel according to the reference probes,
creation of a companion binarized skeleton image. For this latter procedure,
the NIR channel was used, since it presents an a high contrast between the
leaf and background. The skeleton comprises the profile of the leaf and the
vein pattern. Images and companion skeletons were resized at $256\times 256$
resolution. Fig 1 shows the original images in the RGB and RGNIR spaces, the
normalized NIR channel and the corresponding companion skeleton. Before using
the generative algorithms, we performed standard data augmentation by randomly
flipping each image horizontally and vertically, rotating by an angle randomly
chosen in $[-\pi/4,\pi/4]$ and finally zooming with a random amount in the
range $[-20\%,+20\%]$. The dataset was thus increased in this way from 80 to
240 samples.

Fig 1: Sample of grapevine leaf from the dataset. A: RGB image; B: RGNIR
image; C: normalized and cropped NIR image; D: companion skeleton. In the
skeleton binarized image, the white color identifies the leaf profile and
veins, the black color identifies other parts of the leaf and the background.

### Generative methods for L2L translation

The authors of [9, 12] generated artificial patterns of blood vessels along
with corresponding eye fundus images using a common strategy which divides the
problem of the image generation into two sub–problems, each one addressed by a
tailored DL architecture: first they generate the blood vessel tree, then they
color the eye fundus. We adopt this very approach, first generating the leaf
profile and veins and then coloring the leaf blade. Also in our experience
this approach has turned out to be more effective than generating the
synthetic image altogether.

#### Skeleton Generation

According to the above considerations, the generation of a realistic leaf
skeleton is the first step towards the final goal of our work. For this task,
we use a convolutory autoencoder architecture, that is, a network trained to
reconstruct its input. An autoencoder (AE) is composed of two submodels: 1) an
encoder $Q$ that maps the training dataset to a latent (hidden) representation
$z$; 2) a decoder $P$ that maps $z$ to an output that aims to be a plausible
replica of the input. We have experimented that simple autoencoders cannot
generate realistic skeletons. For this reason, we use a more sophisticated
architecture, called Residual Variational Autoencoder (ResVAE, see Fig 2).

Fig 2: Illustration of the ResVAE framework (training).

This learning framework has already been successfully applied to image
recognition, object detection, and image super-resolution (see, e.g., [13]).
In the data generation framework, AEs learn the projection of the initial data
into a _latent subspace_ , and then a sample of this subspace is randomly
extracted to build up a new instance of the initial data. Instead of learning
such projection, VAEs learn the probability distribution of the latent
variables given the input $x$. As a matter of fact, a variational autoencoder
can be defined as an autoencoder whose training is regularized to avoid
overfitting and ensure that the latent space has good properties that enable
the generative process. To achieve this goal, instead of encoding an input as
a single point, VAEs encode it as a (Gaussian) distribution over the latent
space, where $p(z|x)$ represents the probability of the latent variable $z$
given the input $x$. The decoding part consists in sampling a variable from
$p(z|x)$ and then providing a reconstruction $\widehat{x}$ of the initial data
$x$. We associate to this framework the following loss function

$\mathcal{L}_{\text{VAE}}(x,\hat{x})=\mathcal{L}_{L_{2}}(x,\hat{x})+\beta\mathcal{L}_{KL}\left(p(z|x),\mathcal{N}(0,1)\right),$
(1)

where the first term $\mathcal{L}_{L_{2}}=||x-\widehat{x}||^{2}$ is the
$L_{2}$ norm of the reconstruction loss, and the second term
$\mathcal{L}_{KL}=KL[N(\mu_{x},\sigma_{x}),N(0,1)]$ is the Kullback–Leibler
(KL) divergence [14, 15] The KL divergence enhances sparsity in neurons
activation to improve the quality of the latent features keeping the
corresponding distribution close to the Gaussian distribution
$\mathcal{N}(0,1)$. The tunable regularization hyperparameter $\beta$ is used
to weigh the two contributions[16]. With respect to VAEs, ResVAEs additionally
employ residual blocks and connection skips. The idea beyond residual blocks
is the following [17]: normal layers try to directly learn an underlying
mapping, say $h(x)$, while residual ones approximate a residual function
$r(x)=h(x)-x$. Once the learning is complete, $r(x)$ is added to the input to
retrieve the mapping: $h(x)=r(x)+x$. In our architecture, residual blocks are
concatenated to the decoder to increase the capacity of model [13]. The
connection skips allow to back–propagate the gradients more efficiently giving
the bottleneck more access to the simpler features extracted earlier in the
encoder. The resulting ResVAE compresses $256\times 256$ leaf skeleton images
to a low dimension latent vector of size 32 and then it reconstructs it to
$256\times 256$ images. We refer to S1 Appendix for specifications of the
present ResVAE architecture.

#### Translation to colorized leaf image

We consider the colorization of the leaf out of an existing skeleton as an
image-to- image translation problem, which implies to learn a mapping from the
binary vessel map into another representation. Similarly to what observed in
[9] for retinal image generation, many leaf images can share a similar binary
skeleton network due to variations in color, texture, illumination. For this
reason, learning the mapping is an ill-posed problem and some uncertainty is
present. We learn the mapping via a Pix2pix net (also known as conditional
GAN, (cGAN)), an unsupervised generative model which represents a variation of
a standard GAN. As such it includes two deep neural networks, a generator $G$
and discriminator $D$. The generator aims to capture the data distribution,
while the discriminator estimates the probability that a sample actually came
from the training data rather than from the generator. In order to learn a
generative distribution over the data $x$, the generator builds a mapping
$G(z;\theta_{G})$ from a prior noise distribution $p_{z}$ to the image data
space, $\theta_{G}$ being the generator parameters. The discriminator outputs
the probability that $x$ came from the real data distribution $p_{data}(x)$
rather from the generated one. We denote by $D(x;\theta_{D})$ the
discriminator function, $\theta_{D}$ being the discriminator parameters. In
standard GANs, the optimal mappings $G^{*}$ is obtained as the equilibrium
point of the min–max game:

$(G^{*},D^{*})=\arg\displaystyle\min_{G}\max_{D}\mathcal{L}_{GAN}(D,G),$

where we have defined the objective function

$\mathcal{L}_{GAN}(D,G):=\mathbb{E}_{x\sim p_{data}(x)}[\log
D(x;\theta_{D})]+\mathbb{E}_{z\sim p_{z}(z)}[\log(1-D(G(z;\theta_{G})))].$ (2)

In the conditional framework, an extra variable $y$ is added as a further
source of information on $G$, which combines the noise prior $p_{z}(z)$ and
$y$. The objective function thus becomes

$\mathcal{L}_{cGAN}(D,G)=\mathbb{E}_{x\sim p_{data}(x)}[\log
D(x;\theta_{D})]+\mathbb{E}_{z\sim p_{z}(z)}[\log(1-D(G(z|y;\theta_{G})))].$
(3)

Previous approaches have found it beneficial to mix the GAN objective with a
more traditional loss, such as $L_{2}$ distance [18]. The discriminator’s job
remains unchanged, but the generator is bound not only to fool the
discriminator but also to stay near the ground truth output in an $L_{2}$
sense. In this work we rather explore the use of the $L_{1}$ distance rather
than $L_{2}$ as $L_{1}$ promotes sparsity and at the same time it encourages
less blurring [19]:

$\mathcal{L}_{L_{1}}(G)=\mathbb{E}_{x,y,z}[||y-G(z|y;\theta_{G})||].$ (4)

The final objective is thus

$(G^{*},D^{*})=\arg\min_{G}\max_{D}\mathcal{L}_{cGAN}(D,G)+\lambda\mathcal{L}_{L_{1}}(G)$
(5)

where $\lambda$ is a regularization hyperparameter. In our implementation the
extra information corresponds to the leaf skeletons which condition $G$ in the
image generation task to preserve leaf shape and venation pattern. The
discriminator is provided with skeleton plus generated image pairs and must
determine whether the generated image is a plausible (feature preserving)
translation. Fig 3 shows the training process of the cGAN. We refer to S2
Appendix for specifications of the Pix2pix architecture we adopted.

Fig 3:  Illustration of the Pix2Pix framework (training).

#### L2L workflow: from random samples to leaf images

Upon training of the ResVAE and Pix2pix architectures, we dispose of an end-
to-end procedure for the generation of synthetic leaves. The procedure, which
is completely unsupervised, can be summarized as follows (see also Fig 4):

  1. 1.

Load weights of the trained ResVAE decoder and Pix2pix generator.

  2. 2.

Draw a random vector from a normal distribution whose parameters are chosen
according to the ResVAE latent space representation (note that its size equals
the dimension of the latent space used in the ResVAE, 32 in the present case).

  3. 3.

Input the random vector in the trained ResVAE decoder and generate a leaf
skeleton

  4. 4.

Input the leaf skeleton into the trained generator of the Pix2Pix net to
translate it into a fully colorized leaf.

Fig 4:  L2L workflow illustration. A random input vector is drawn from the
ResVAE latent space representation and is input into the trained ResVAE
decoder. This latter outputs a synthetic leaf skeleton, which in turn is fed
into the trained generator of the Pix2Pix and translated into a corresponding
colorized leaf.

## Results and Discussion

The proposed technique can be employed to generate as many synthetic leaf
images as the user requires. The model has been implemented with Keras 111Code
to reproduce our experiments will be made available upon publication of this
work.. Upon generation of the synthetic images, their quality is assessed
performing both experimental qualitative (visual) and quantitative evaluations
as follows.

### Visual qualitative evaluation

Consistency test. Beforehand, we have evaluated the consistency of the
methodology by verifying that the net has learned to translate a leaf sample
comprised in the training set into itself. Fig 5 shows an example of this
test. The generated leaf is very similar to the real one, except for some vein
discoloration and a small blurring effect effect, which is a well–known
product of AEs employed in image generation [20].

Fig 5:  Consistency test. The companion binarized leaf skeleton of a real leaf
is passed through the generator of the Pix2Pix net to check whether the
synthetic colorized leaf blade is similar to the original one. A: companion
skeleton of a real leaf; B: synthetic colorized blade generated; C: original
leaf.

Translation from unseen real companion skeleton. Having ensured that the model
has learned to translate on the training data, we verify that it is able to
produce reliable synthetic images using skeletons obtained from leaves that
are not part of the training dataset. Fig 6 shows an instance of colorized
leaf obtained from this test.

Fig 6:  Translation from unseen real companion skeleton. A binarized leaf
skeleton companion of a real leaf not belonging to the training set is passed
through the generator of the Pix2Pix net to check. A: companion skeleton; B:
synthetic colorized blade.

Full L2L translation Fig 7 shows several instances of synthetic colorized
leaves obtained starting from different random latent vectors. Note that the
generated leaf images differ in terms of their global appearance, that is the
model generalizes and does not trivially memorizes the examples. As a note,
one should observe that some discolored parts may be appear. Moreover,
sometimes the skeletons show small artifacts consisting in not–connected
pixels positioned outside the leaf boundary (not appearing in Fig 7). This
latter issue will be addressed via a refinement algorithm explained below.

Fig 7:  Full L2L translation results. Examples of synthetic colorized leaves
along with the corresponding synthetic companion skeletons.

L2L-RGNIR translation As mentioned above, applications in crop management
require to have at disposal images also in the NIR channel. To do this, we use
the L2L generation procedure as for the RGB channels starting from RGNIR
images as Pix2Pix targets. Since the same leaf skeletons are used, it is not
necessary to re-train the ResVAE if this procedure has been already carried
out for the RGB case. Fig 8 shows some results of this model.

Fig 8:  L2L-RGNIR translation results. Examples of synthetic leaves colorized
in the RGNIR channels along with the corresponding synthetic companion
skeletons.

Refinement algorithm. We have already discussed the fact that synthetically
generated images may sometimes present artifacts (leaf regions that appear
detached from the leaf blade). Obviously this is not realistic and we need to
remove such artifacts. The refinement algorithm is implemented at present in a
procedural way and it is based on the idea of finding the contours of all the
objects and removing all objects laying outside the leaf contour. Note that
this procedure must pay attention to leave internal holes intact, because in
nature such holes are the result of the superposition of leaf lobes or due to
several abiotic/biotic conditions. Fig 9 shows the first leaf in Fig 7 which
presents artifacts (panel A, zoomed area including the artifact in panel B)
and its cleaned counterpart (panel C).

Fig 9:  Refinement algorithm. The generative procedure sometimes produces
artifacts, that is leaf regions that appear outside the leaf blade. These
artifacts are corrected by procedurally finding the contours of all the
objects in the image and removing the objects outside the leaf contour. A:
first leaf in Fig 7 presenting artifacts; B: inset showing the magnified
artifacts; C: cleaned leaf.

### Quantitative Quality Evaluation

In order to assess quantitatively the quality of the generated leaves, we
employ a DL–based anomaly detection strategy. This approach is discussed in
detail in [21], here we briefly recall the main points. The strategy consists
in training an AE to compress real leaf images in a latent subspace and then
reconstruct the images using the latent representation (see Skeleton
Generation section for the same concept). Once the network is trained in this
way, we feed it with a synthetic image generated by our procedure. The AE
encodes it in the latent space and tries to recover the original image
according to its training rules. Since the net has been trained to be the
identity operator for real images, if the artificial images are substantially
different, an anomalous reconstruction is obtained. Fig 10 provides a visual
schematization of this approach. The figure also details the score system used
to detect the anomaly.

Fig 10: AE for anomaly detection. The AE is trained with images of real leaves
to be the identity operator of the input. A synthetic leaf with a low level of
similarity is recognized as an anomaly if fed into the trained AE and its
anomaly score $s_{x}$ is high.

The degree of anomaly is quantified via the ROC curve and its area, the AUC
index [22]. For this latter, we found AUC=0.25, which means that for a random
synthetic image fed into the AE, there is a 25% of possibility to classify it
as an anomaly, that is to be synthetic instead or real. While this result does
not indicate a perfect reconstruction of the real leaves, it shows that the
synthetic leaves are a reasonably accurate surrogate of real leaves and can be
used for a first massive training at a very low cost. A successive refinement
can then be applied using a limited number of real leaves upon transfer
learning.

Fig 11: Quantification of anomaly via ROC curve and AUC index. A point on the
ROC curve represents - for a certain threshold on the anomaly score - the
false positive rate (genuinely real images) vs the true positive rate
(genuinely synthetic images). The value AUC=0.25 means that a synthetic image
is (mis–)classified as synthetic in the 25% of cases. The dotted line
represents the result one would obtain by tossing a coin to decide whether an
image is artificial or real.

## Conclusion

Goal of this work was to explore advanced DL generative methods to produce
realistic images of leaves to be used in computer–aided applications The main
focus was on the generation of artificial samples of leaves to be used to
train DL networks for modern crop management systems in precision agriculture.
Disposing of synthetic samples which have a reasonable resemblance to real
samples alleviates the burden of manually collecting and annotating hundreds
of data. The Pix2pix net performs good translations from the leaf skeletons
generated by the ResVAE, except for some discolored parts, both for the
colorization of RGB and RGNIR images. Also, the leaves generated by ResVAE
have sometimes pixels positioned outside the boundary which, if not corrected,
can cause artifacts in the synthetic leaves. An easy procedure has been
proposed as well to correct these artifacts. We believe that the generative
approach can significantly contribute to automatize the process of building a
low-cost training set for DL applications. Several computer–aided applications
may also benefit of such a strategy, where many samples are required, possibly
with different degree of accuracy in the representation.

## Author contribution

Conceptualization: Alessandro Benfenati, Paola Causin
Dataset: Alessandro Benfenati, Davide Bolzi, Paola Causin, Roberto Oberti
Methodology: Alessandro Benfenati, Davide Bolzi, Paola Causin
Implementation: Davide Bolzi
Analysis: Alessandro Benfenati, Davide Bolzi, Paola Causin, Roberto Oberti
Writing: Alessandro Benfenati, Paola Causin

## Supporting information

##### S1 Appendix

Implementation and training of the ResVAE neural network. The architecture,
inspired by the one described in [23], is shown in Fig 12.

Fig 12:  ResVAE architecture. Building blocks of the encoder and decoder
components of the ResVAE. The convolutional filters have kernels of size
$4\times 4$. The Residual block is formed by 5 convolutional layers of 16
filters each with kernel of size $4\times 4$ and stride equal to 1, followed
by a Batch Normalization layer and LeakyReLU activation function.

The training is performed via a stochastic gradient descent strategy, with
gradients computed by standard back–propagation; we use the Adam optimizer
with learning rate $\eta=0.001$ and we train the model for 2000 epochs with a
batch size of 64. After a hyper–parameter search, $\beta$ in the loss function
(1) was set to 75.

##### S2 Appendix

Implementation and training of the Pix2pix neural network.

The Pix2pix is a GAN architecture designed for image-to-image translation,
originally presented in [19] and comprising a generator and a discriminator.
The discriminator is deep neural network that performs image classification.
It takes both the source image (leaf skeleton) and the target image (colorized
leaf) as input and predicts the likelihood of whether the target image is real
or a fake translation of the source image. We use a PatchGAN model which tries
to establish whether each $N\times N$ (local) patch in the image is real or
fake. We run this discriminator convolutionally across the image, averaging
all responses to provide the ultimate output of the discriminator. The
generator is an encoder-decoder model using a U-Net architecture with feature-
map concatenation between two corresponding blocks of the encoder/decoder. The
encoder and decoder of the generator are comprised of standardized blocks of
convolution, batch normalization, dropout, and activation layers. We proceed
as suggested in [19]: the generator is updated via a weighted sum of both the
adversarial loss and the $L_{1}$ loss, where the parameter $\lambda$ in the
loss function eq5 is set to 100 in order to encourage the generator to produce
plausible translations of the input image, and not just plausible images in
the target domain. We initialize the generator/discriminator weights with a
normal distribution of zero mean and standard deviation $\sigma=0.002$; we use
the Adam optimizer with a learning rate $\eta=0.0002$ and we train the
generator/discriminator paired model for 12000 training steps, using a batch
size of 1. Fig 13 shows the generator and discriminator architectures.

Fig 13:  Pix2pix architecture. Building blocks of the generator and
discriminator components.

## Acknowledgments

We acknowledge support from the SEED PRECISION project (PRecision crop
protection: deep learnIng and data fuSION), funded by Università degli Studi
di Milano. AB and PC are part of the GNCS group of INDAM (Istituto Nazionale
di Alta Matematica ”Francesco Severi”).

## References

  * 1. Peyrat A, Terraz O, Merillou S, Galin E.  Generating vast varieties of realistic leaves with parametric 2Gmap L-systems.  Vis Comput. 2008;24(7):807–816.
  * 2. Runions A, Fuhrer M, Lane B, Federl P, Rolland-Lagan AG, Prusinkiewicz P.  Modeling and visualization of leaf venation patterns.  In: ACM SIGGRAPH 2005 Papers; 2005. p. 702–711.
  * 3. Samee SB.  Modeling and Simulation of Tree Leaves Using Image-Based Finite Element Analysis.  PhD Thesis, University of Cincinnati; 2012.
  * 4. Kim D, Kim J.  Procedural modeling and visualization of multiple leaves.  Multimed Sys. 2017;23(4):435–449.
  * 5. Féret JB, Gitelson A, Noble S, Jacquemoud S.  PROSPECT-D: towards modeling leaf optical properties through a complete lifecycle.  Remote Sens Environ. 2017;193:204–215.
  * 6. Miao T, Zhao C, Guo X, Lu S.  A framework for plant leaf modeling and shading.  Math Comput Model. 2013;58(3-4):710–718.
  * 7. Yi Wl, He Hj, Wang Lp, Yang Hy.  Modeling and simulation of leaf color based on virtual rice.  DEStech Trans Mater Sci Eng. 2016;(mmme).
  * 8. Taghanaki SA, Abhishek K, Cohen JP, Cohen-Adad J, Hamarneh G.  Deep semantic segmentation of natural and medical images: a review.  Artif Intell Rev. 2021;54(1):137–178.
  * 9. Costa P, Galdran A, Meyer MI, Niemeijer M, Abràmoff M, Mendonça AM, et al.  End-to-end adversarial retinal image synthesis.  IEEE Trans Med Imaging. 2017;37(3):781–791.
  * 10. Oberti R, Marchi M, Tirelli P, Calcante A, Iriti M, Borghese AN.  Automatic detection of powdery mildew on grapevine leaves by image analysis: Optimal view-angle range to increase the sensitivity.  Comput Electron Agric. 2014;104:1–8.
  * 11. Mahlein AK, Kuska MT, Behmann J, Polder G, Walter A.  Hyperspectral sensors and imaging technologies in phytopathology: state of the art.  Annu Rev Phytopathol. 2018;56:535–558.
  * 12. Sengupta S, Athwale A, Gulati T, Zelek J, Lakshminarayanan V.  FunSyn-Net: enhanced residual variational auto-encoder and image-to-image translation network for fundus image synthesis.  In: Medical Imaging 2020: Image Processing. vol. 11313. International Society for Optics and Photonics; 2020. p. 113132M.
  * 13. Cai L, Gao H, Ji S.  Multi-stage variational auto-encoders for coarse-to-fine image generation.  In: Proceedings of the 2019 SIAM International Conference on Data Mining. SIAM; 2019. p. 630–638.
  * 14. Asperti A, Trentin M.  Balancing Reconstruction Error and Kullback-Leibler Divergence in Variational Autoencoders.  IEEE Access. 2020;8:199440–199448.
  * 15. Benfenati A, Ruggiero V.  Image regularization for Poisson data.  Journal of Physics: Conference Series. 2015 nov;657:012011.  Available from: https://doi.org/10.1088/1742-6596/657/1/012011.
  * 16. Higgins I, Matthey L, Pal A, Burgess C, Glorot X, Botvinick M, et al.  $\beta$-VAE: Learning basic visual concepts with a constrained variational framework.  In: ICLR 2017 Conference Proceedings; 2017. .
  * 17. He K, Zhang X, Ren S, Sun J.  Deep residual learning for image recognition.  In: Proceedings of the IEEE conference on computer vision and pattern recognition; 2016. p. 770–778.
  * 18. Kurach K, Lučić M, Zhai X, Michalski M, Gelly S.  A large-scale study on regularization and normalization in GANs.  In: International Conference on Machine Learning. PMLR; 2019. p. 3581–3590.
  * 19. Isola P, Zhu JY, Zhou T, Efros AA.  Image-to-image translation with conditional adversarial networks.  In: Proceedings of the IEEE conference on computer vision and pattern recognition; 2017. p. 1125–1134.
  * 20. Huang H, Li Z, He R, Sun Z, Tan T.  Introvae: Introspective variational autoencoders for photographic image synthesis.  arXiv preprint arXiv:180706358. 2018;.
  * 21. Benfenati A, Causin P, Oberti R, Stefanello G.  Unsupervised feature–oriented deep learning techniques for powdery mildew recognition based on multispectral imaging.  submitted;.
  * 22. Bradley AP.  The use of the area under the ROC curve in the evaluation of machine learning algorithms.  Pattern Recognit. 1997;30(7):1145–1159.
  * 23. Chollet F. Variational AutoEncoder; 2020.  https://keras.io/examples/generative/vae/.


$h_{\pi\pi,0}=$ | $-0.09\pm 0.60$ | $\begin{bmatrix}[r]1&-0.5&-0.5&0.6\;\;\;\;\;&-0.2&0.3&-0.4&0.2&{\bf-1}&-0.5\\\ &1&-0.4&0.4\;\;\;\;\;&0&-0.6&{\bf 1}&-0.2&0.5&{\bf 1}\\\ &&1&{\bf-1}\;\;\;\;\;&0.1&0.3&-0.5&0.1&0.5&-0.5\\\ &&&1\;\;\;\;\;&-0.1&-0.3&0.5&-0.2&-0.6&0.4\\\\[6.88889pt] &&&&1&-0.2&0&-0.2&0.1&0\\\ &&&&&1&-0.6&0.6&-0.2&-0.6\\\ &&&&&&1&-0.3&0.4&{\bf 1}\\\ &&&&&&&1&-0.1&-0.2\\\ &&&&&&&&1&0.5\\\ &&&&&&&&&1\end{bmatrix}$
---|---|---
$h_{\pi\pi,1}=$ | $0.001\pm 0.400$
$h_{\overline{K}K,0}=$ | $0.8\pm 1.6$
$h_{\overline{K}K,1}=$ | $-0.28\pm 0.50$
$m=$ | $0.1338\,(5)\cdot a_{t}^{-1}$
$g_{\pi\pi}=$ | $0.441\,(9)$
$g_{K\overline{K}}=$ | $0.17\,(30)$
$\gamma_{\pi\pi,\pi\pi}=$ | $(2.9\pm 0.9)\cdot a_{t}^{2}$
$\gamma_{\pi\pi,K\overline{K}}=$ | $-(2.4\pm 5.0)\cdot a_{t}^{2}$
$\gamma_{K\overline{K},K\overline{K}}=$ | $-(2.2\pm 4.0)\cdot a_{t}^{2}$
$\chi^{2}/N_{\text{dof}}=\frac{126.9}{32-4}=4.53$ ,

(51)

where we present also the parameters of the scattering amplitude to illustrate
the correlation between the functions $\mathcal{F}_{a}$ and $\mathcal{M}$
181818The quoted $\chi^{2}$ describes only the variation of the smooth
function parameters, $h$. The slight variations of the correlation matrix in
Eq. 51 with respect to what is reported in Eq. 16 are due to the fact that
only a subset of 348 configurations out of the 400 available to calculate the
spectrum are used for the extraction of the form factors.. Note that some of
the smooth function parameters, $h$, are maximally correlated or
anticorrelated with parameters in the scattering amplitude.

Even though the smooth functions are individually consistent with zero, when
weighted by the finite-volume correction factors or the scattering amplitude,
the resulting values are not compatible with zero. For this to occur, it is
necessary, although not sufficient, that the functions $\mathcal{F}_{a}$ have
a significant correlation with the scattering amplitude $\mathcal{M}$, which
can be seen in Eq. 51.

##  Appendix H Spacelike form factor of the pion and renormalization constant

The pion form-factor in the _spacelike_ region was extracted from three-point
correlation functions, $\langle 0|\Omega_{\pi}(\Delta
t)\,\mathcal{J}(t)\,\Omega^{\dagger}_{\pi}(0)|0\rangle$, computed with a
single value of $\Delta t=32\,a_{t}$. Details of the computational approach
are presented in Ref. Radhakrishnan  _et al._ (2022) and Ref. Shultz  _et al._
(2015). In order to cover a wide kinematic region, correlators were computed
giving access to matrix elements,
$\matrixelement{\pi(\mathbf{p}_{1})}{\mathcal{J}^{i}_{\rho,\text{lat}}}{\pi(\mathbf{p}_{2})}$,
for combinations of pion momenta up to $|\mathbf{p}_{i}|^{2}\leq
6\,\left(\tfrac{2\pi}{L}\right)^{2}$, and current momentum insertion up to
$|\mathbf{p}_{1}-\mathbf{p}_{2}|^{2}\leq 4\,\left(\tfrac{2\pi}{L}\right)^{2}$.

In a previous analysis of some of these correlation functions in Ref.
Radhakrishnan  _et al._ (2022), the leading time dependence was removed by
forming the combination,

$\widetilde{C}_{\text{3pt}}(\Delta t,t)=\frac{\langle 0|\Omega_{\pi}(\Delta
t)\,\mathcal{J}(t)\,\Omega^{\dagger}_{\pi}(0)|0\rangle}{e^{-E_{\mathbf{p}_{1}}(\Delta
t-t)}e^{-E_{\mathbf{p}_{2}}t}}\,,$ (52)

where $E_{\mathbf{p}}$ corresponds to the energy of a single-pion state of
momentum $\mathbf{p}$, and where the normalization of the optimized operators
follows the same convention as in the main text.

It can be the case that the timeslice-to-timeslice data correlation for this
quantity is considerable, resulting in fits with reasonable values of
$\chi^{2}$ which undershoot the data. One such case is presented in panel (a)
of Fig. 24.

An alternative approach is to form a ratio using optimized two-point
correlation functions,

$R_{\text{3pt}}(\Delta
t,t)=4E_{\mathbf{p}_{1}}E_{\mathbf{p}_{2}}\frac{\expectationvalue{\Omega_{\pi}(\Delta
t)\,\mathcal{J}(t)\,\Omega^{\dagger}_{\pi}(0)}{0}}{\expectationvalue{\Omega_{\pi}(\Delta
t\\!-\\!t)\,\Omega^{\dagger}_{\pi}(0)}{0}\expectationvalue{\Omega_{\pi}(t)\,\Omega^{\dagger}_{\pi}(0)}{0}}\,,$
(53)

which will have the same constant contribution, but differing excited-state
contributions. This combination proves to have much smaller timeslice-to-
timeslice data correlation, and fits follow more closely the data points. This
is illustrated in panel (b) of Fig. 24. Fits to a constant, and a constant
with an excited state exponential at source or sink or both are carried out
for a range of time windows, and the results averaged using the AIC as in the
two-point function case discussed previously. The columns on the right
describe the time window of the fit $[t_{\text{min}},$ $t_{\text{max}}]$, and
the number of exponentials at the source, $n_{\text{src}}$, and sink,
$n_{\text{snk}}$.

(a)  (b)

Figure 24:  Fits to three-point correlation functions with
$\mathbf{p}_{1}=\mathbf{p}_{2}=\tfrac{2\pi}{L}[110]$ (averaged over rotations
and directions of current insertion) and fixed $\Delta t=32\,a_{t}$ using
either (a) Eq. 52 or (b) Eq. 53. Fitted constant value in this case
corresponds to $1/Z_{V}^{\ell}$, the vector current renormalization constant.
Variation of fit window is shown in the right columns, along with the model-
averaged result.

The difference with respect to the previous method using
$\widetilde{C}_{\text{3pt}}(\Delta t,t)$ is modest, but is the origin of any
differences in the current analysis with that in Ref. Radhakrishnan  _et al._
(2022), such as for the light-quark vector current renormalization factor as
shown in Fig. 25.

Figure 25: Determination of light-quark vector current renormalization factor
$Z_{V}^{l}$ by correlated fitting of extractions from six values of pion
momentum.