# Two conjectures on the Stokes complex in three dimensions on Freudenthal
meshes

Patrick E. Farrell Mathematical Institute, University of Oxford, Oxford,
England<EMAIL_ADDRESS>Lawrence Mitchell NVIDIA, Santa
Clara, CA, USA<EMAIL_ADDRESS>L. Ridgway Scott Department of
Computer Science, University of Chicago, Chicago, USA<EMAIL_ADDRESS>
###### Abstract

In recent years a great deal of attention has been paid to discretizations of
the incompressible Stokes equations that exactly preserve the
incompressibility constraint. These are of substantial interest because these
discretizations are pressure-robust, i.e. the error estimates for the velocity
do not depend on the error in the pressure. Similar considerations arise in
nearly incompressible linear elastic solids. Conforming discretizations with
this property are now well understood in two dimensions, but remain poorly
understood in three dimensions. In this work we state two conjectures on this
subject. The first is that the Scott–Vogelius element pair is inf-sup stable
on uniform meshes for velocity degree $k\geq 4$; the best result available in
the literature is for $k\geq 6$. The second is that there exists a stable
space decomposition of the kernel of the divergence for $k\geq 5$. We present
numerical evidence supporting our conjectures.

##  1 Introduction

We consider two closely related problems for a bounded Lipschitz domain
$\Omega\subset\mathbb{R}^{d}$, $d\in\\{2,3\\}$. The first is the
incompressible Stokes equations: given $f\in L^{2}(\Omega;\mathbb{R}^{d})$,
find the velocity $u:\Omega\to\mathbb{R}^{d}$ and pressure
$p:\Omega\to\mathbb{R}$ such that

(1a) $\displaystyle-\nabla\cdot{\boldsymbol{\varepsilon}}u+\nabla p$
$\displaystyle=f\quad$ $\displaystyle\text{ in }\Omega,$ (1b)
$\displaystyle\nabla\cdot u$ $\displaystyle=0\quad$ $\displaystyle\text{ in
}\Omega,$ (1c) $\displaystyle u$ $\displaystyle=0\quad$ $\displaystyle\text{
on }\partial\Omega,$

where ${\boldsymbol{\varepsilon}}u=1/2\left(\nabla u+(\nabla u)^{T}\right)$ is
the symmetric gradient of $u$. The second is the Navier–Cauchy equation of
linear elasticity: given $f\in L^{2}(\Omega;\mathbb{R}^{d})$ and $\gamma>0$,
find the displacement $u:\Omega\to\mathbb{R}^{d}$ that satisfies

(2a)
$\displaystyle-\nabla\cdot{\boldsymbol{\varepsilon}}u-\gamma\nabla\nabla\cdot
u$ $\displaystyle=f\quad$ $\displaystyle\text{ in }\Omega,$ (2b)
$\displaystyle u$ $\displaystyle=0\quad$ $\displaystyle\text{ on
}\partial\Omega.$

Here $\gamma=2\lambda/\mu$, where $\lambda$ and $\mu$ are the Lamé parameters.
As $\gamma\to\infty$, the material is said to be nearly incompressible. The
term $-\nabla\nabla\cdot u$ in Eq. 2a is connected to the incompressibility
constraint Eq. 1b; it arises in the Stokes momentum equation Eq. 1a when
employing an augmented Lagrangian approach [6] to enforcing the divergence-
zero constraint Eq. 1b.

When discretizing Eq. 1, it is highly desirable to choose spaces $V_{h}$ for
the velocity and $\Pi_{h}$ for the pressure that are divergence-free, i.e.
such that the incompressibility constraint Eq. 1b is satisfied exactly on the
discrete level [27]. Achieving this is difficult; no element pair for exact
enforcement is known that is simultaneously inf-sup stable, low-order,
conforming, has polynomial basis functions, and is effective on general
meshes. On simplicial grids with special mesh structure, it is possible to use
the conforming Scott–Vogelius finite element pair
$[\mathrm{CG}_{k}]^{d}$-$\mathrm{DG}_{k-1}$ [48, 49] for $k\geq d$ (for Alfeld
splits [39, 22, 17, 52]) or $k\geq d-1$ (for Powell–Sabin splits [20, 53, 55,
21]). The approach of Guzmán and Neilan [22] is conforming, works for
arbitrary degree and on general meshes, but requires the use of piecewise
polynomial basis functions on each cell (instead of standard polynomials).
Non-conforming divergence-free discretizations are reviewed in John et al.
[27, §4.4]. Finally, another approach is to consider the use of high-order
discretizations, which are attractive for their advantageous computational
properties on modern architectures [49, 48, 37, 54]. Another alternative is to
modify the right-hand side of the problem with an operator that maps
discretely divergence-free test functions to exactly divergence-free ones
[30].

In this paper we state two conjectures regarding the discretization and
multigrid solution of Eq. 1 and Eq. 2. The first is that the Scott–Vogelius
element pair is inf-sup stable on structured uniform tetrahedral meshes
(Freudenthal meshes [7]) for $k\geq 4$; the best available result is that of
Zhang, who proved that the pair is stable for $k\geq 6$ [54]. We conjecture
this on the basis of numerical calculations of the inf-sup constant for
varying $k$. These rely on a new algorithm that can compute the inf-sup
constant for elements that are divergence-free when an exact characterization
of the pressure space is not known. The second is that on the same meshes the
subspace $Z_{h}\subset V_{h}$ of divergence-free functions admits a local
basis defined on the vertex-centred patches for $k\geq 5$, and that the
associated space decomposition is stable. This is significant because
identifying a local basis for the kernel of the divergence operator is crucial
for multigrid algorithms applied to Eq. 2. We conjecture this on the basis of
observed $\gamma$-robustness of a multigrid solution algorithm for Eq. 2, for
which the local kernel decomposition is essential [43]. The existence of a
local basis is known in three dimensions for Alfeld splits [17] and
Worsey–Farin splits [21], but remains an open question for general meshes, and
in particular for the Freudenthal meshes considered here.

##  2 Inf-sup stability of the Scott–Vogelius element

Figure 1: Unit cube for a structured mesh in three dimensions used for
computational experiments. Given $N\in\mathbb{N}$, a $N\times N\times N$ mesh
of cubes is generated, each of which is subdivided into 6 tetrahedra as shown.

The mixed formulation of Eq. 1 is to find $(u,p)\in V\times\Pi\coloneqq
H^{1}_{0}(\Omega;\mathbb{R}^{d})\times L^{2}_{0}(\Omega)$ such that

(3)
$({\boldsymbol{\varepsilon}}u,{\boldsymbol{\varepsilon}}v)_{L^{2}}-(p,{{\nabla\cdot}\,}v)_{L^{2}}-(q,{{\nabla\cdot}\,}u)_{L^{2}}=\langle
f,v\rangle\quad\text{ for all }(v,q)\in V\times\Pi.$

Define $\langle T,U\rangle=(T,U)_{L^{2}}$ for any scalar-, vector-, or tensor-
valued functions $T,U$.

The inf-sup condition determines whether or not a pair of spaces $V_{h}\subset
V$, $\Pi_{h}\subset\Pi$ in a mixed finite element method provide a compatible
discretization [9], since the symmetric gradient term is coercive on $V$. In
the context of Eq. 3, this condition encapsulates a constraint between the
divergence of the velocity space $V_{h}$ and the pressure space $\Pi_{h}$:
there exists $\beta>0$ such that

(4) $\beta\,\sqrt{\langle q,q\rangle}\leq\sup_{v\in
V_{h}}\frac{\langle{{\nabla\cdot}\,}v,q\rangle}{\sqrt{\langle\nabla v,\nabla
v\rangle}}\quad\forall q\in\Pi_{h}.$

In this section we make the following conjecture:

######  Conjecture 1.

Let $V_{h}$ be constructed with continuous Lagrange elements of degree $k$,
and choose $\Pi_{h}={{\nabla\cdot}\,}V_{h}$. For $k\geq 4$, the inf-sup
condition Eq. 4 holds on structured tetrahedral meshes of the form in Fig. 1,
known as the Freudenthal triangulation [7], with a constant that depends only
on $k$.

We conjecture this on the basis of numerical computations with a new algorithm
for calculating the inf-sup constant, which we now describe.

###  2.1 Computing the inf-sup constant

There are several methods to estimate computationally the inf-sup constant
$\beta$ for various spaces and variational problems [39, 41, 4]. In [18],
approximation of the corresponding Ladyzhenskaya inf-sup constant for the
continuous problem is studied. All of these are based on solving
eigenproblems.

In [41], an automated system ASCoT was introduced for computing inf-sup
constants. It takes as input $V_{h}$ and $\Pi_{h}$. However, when
$\Pi_{h}={{\nabla\cdot}\,}V_{h}$, it is advantageous to exploit this structure
in the calculation of the inf-sup constant. More importantly, the space
${{\nabla\cdot}\,}V_{h}$ may not be known a priori, and so a different
algorithm is needed that does not require $\Pi_{h}$ as input.

A technique was developed by one of the authors in [46] that deals
specifically with the $\Pi_{h}={{\nabla\cdot}\,}V_{h}$ case and casts this as
an eigenvalue problem, as follows. Define the bilinear forms

(5) $a(u,v)=\int_{\Omega}\nabla u:\nabla v\,dx,\qquad
b(v,q)=\int_{\Omega}({{\nabla\cdot}\,}v)q\,dx.$

First define

(6) $Z_{h}=\left\\{v\in V_{h}\;:\;b(v,q)=0\;\forall q\in\Pi_{h}\right\\}.$

Define $\kappa$ by

(7) $\kappa=\min_{0\neq v\in
V_{h},\;v\perp_{a}Z_{h}}\frac{({{\nabla\cdot}\,}v,{{\nabla\cdot}\,}v)_{L^{2}}}{a(v,v)}=\min_{0\neq
v\in
Z_{h}^{\perp}}\frac{({{\nabla\cdot}\,}v,{{\nabla\cdot}\,}v)_{L^{2}}}{a(v,v)},$

where $v\perp_{a}Z_{h}$ means that $a(v,w)=0$ for all $w\in Z_{h}$ and

(8) $Z_{h}^{\perp}=\left\\{v\in V_{h}\;:\;a(v,w)=0\;\forall w\in
Z_{h}\right\\}.$

We recall the following Lemma from [46, Lemma 26.1].

######  Lemma 2.1.

Suppose that $\|v\|_{V}=\sqrt{a(v,v)}$ and $\Pi_{h}={{\nabla\cdot}\,}V_{h}$.
Then

(9)
$\beta=\underset{\vbox{\Let@\restore@math@cr\default@tag\halign{\hfil$\m@th\scriptstyle#$&$\m@th\scriptstyle{}#$\cr&{q}\in{\Pi_{h}}\\\
&{q} \neq 0\crcr}}}{\mathrm{inf\mathrlap{\phantom{p}}}}\
\underset{\vbox{\Let@\restore@math@cr\default@tag\halign{\hfil$\m@th\scriptstyle#$&$\m@th\scriptstyle{}#$\cr&{v}\in{V_{h}}\\\
&{v} \neq 0\crcr}}}{\mathrm{sup}}\
\frac{b(v,q)}{\|v\|_{V}\|q\|_{L^{2}}}=\underset{\vbox{\Let@\restore@math@cr\default@tag\halign{\hfil$\m@th\scriptstyle#$&$\m@th\scriptstyle{}#$\cr&{q}\in{\Pi_{h}}\\\
&{q} \neq 0\crcr}}}{\mathrm{inf\mathrlap{\phantom{p}}}}\
\underset{\vbox{\Let@\restore@math@cr\default@tag\halign{\hfil$\m@th\scriptstyle#$&$\m@th\scriptstyle{}#$\cr&{v}\in{Z_{h}^{\perp}}\\\
&{v} \neq 0\crcr}}}{\mathrm{sup}}\
\frac{b(v,q)}{\|v\|_{V}\|q\|_{L^{2}}}\geq\sqrt{\kappa}\geq{\textstyle{1\over
2}}\beta,$

where $\kappa$ is defined in Eq. 7.

This result holds also for the bilinear form

(10)
$a(u,v)=({\boldsymbol{\varepsilon}}u,{\boldsymbol{\varepsilon}}v)_{L^{2}},$

with possibly different constants.

Computing $\kappa$ is equivalent to finding the smallest eigenvalue $\lambda$
of the following eigenproblem: find $0\neq u_{h}\in Z_{h}^{\perp}$ such that

(11)
$({{\nabla\cdot}\,}u_{h},{{\nabla\cdot}\,}v_{h})=\lambda\,a(u_{h},v_{h})\quad\forall
v_{h}\in Z_{h}^{\perp},$

which is equivalent to the Rayleigh quotient minimization Eq. 7. Note that
$\kappa>0$ since $\kappa=0$ leads to the contradiction
${{\nabla\cdot}\,}u_{h}=0$, that is, $u_{h}\in Z_{h}\cap Z_{h}^{\perp}$. Thus
there are no spurious modes when $\Pi_{h}={{\nabla\cdot}\,}V_{h}$.

We can write Eq. 11 in operator form as

(12) $Bu_{h}=\lambda Au_{h}.$

The equation Eq. 12 is a symmetric generalized eigenvalue problem [50], and
its eigenvalues are all real. There are many algorithms for solving symmetric
generalized eigenvalue problems [19, 42, 50]. However, we do not have an
explicit basis for the space $Z_{h}^{\perp}$, so we are forced to use matrix-
free methods. Here we focus on a simple method related to the power method.

Since $A$ is invertible on all of $V_{h}$, it would be attractive to utilize
an iteration in which we invert $A$, and not $B$. We introduce a shift
$\sigma$:

(13) $(B-\sigma A){u}_{h}=\lambda^{\sigma}A{u}_{h}.$

If $\lambda_{1}\leq\lambda_{2}\leq\cdots\leq\lambda_{N}$ are the eigenvalues
of the symmetric problem Eq. 12, that is, for the operator $A^{-1}B$, and
$\lambda^{\sigma}_{1}\leq\lambda^{\sigma}_{2}\leq\cdots\leq\lambda^{\sigma}_{N}$
are the eigenvalues of the shifted problem Eq. 13, then
$\lambda^{\sigma}_{i}=\lambda_{i}-\sigma$ for all $i=1,\dots,N$. Moreover,
$\lambda^{\sigma}_{i}$ are the eigenvalues for the operator $A^{-1}B-\sigma
I$, and the eigenvectors ${u}_{h}^{i}$ for $\lambda^{\sigma}_{i}$ and
$\lambda_{i}$ are the same, for all $i=1,\dots,N$.

###  2.2 Using the power method

We can solve for certain eigenvalues and eigenvectors ${u}_{h}$ via the power
method [45] for Eq. 13, namely to find ${u}^{n}\in Z_{h}^{\perp}$ such that

(14)
$\begin{split}a({u}^{n},{v}_{h})&=({{\nabla\cdot}\,}\hat{u}^{n-1},{{\nabla\cdot}\,}{v}_{h})-\sigma\,a(\hat{u}^{n-1},{v}_{h})\quad\forall\;{v}_{h}\in
Z_{h}^{\perp}\\\
\lambda^{n+1}&=\frac{({{\nabla\cdot}\,}{u}^{n},{{\nabla\cdot}\,}{u}^{n})}{a({u}^{n},{u}^{n})},\qquad\hat{u}^{n}=(\|{u}^{n}\|_{V})^{-1}{u}^{n}.\end{split}$

The choices of initial iterate and termination criterion are discussed
subsequently.

We now consider how to compute this, despite the fact that the space
$Z_{h}^{\perp}$ is not explicitly known. Suppose that $\hat{u}^{0}\in
Z_{h}^{\perp}$ is given and that we solve for ${u}^{1}\in V_{h}$ via

$a({u}^{1},{v}_{h})=({{\nabla\cdot}\,}\hat{u}^{0},{{\nabla\cdot}\,}{v}_{h})-\sigma\,a(\hat{u}^{0},{v}_{h})\quad\forall\;{v}_{h}\in
V_{h},$

which we can do since $a(\cdot,\cdot)$ is coercive on $V_{h}$. Then for all
${v}_{h}\in Z_{h}$

$a({u}^{1},{v}_{h})=({{\nabla\cdot}\,}\hat{u}^{0},{{\nabla\cdot}\,}{v}_{h})-\sigma\,a(\hat{u}^{0},{v}_{h})=0$

since both terms vanish. Thus ${u}^{1}\in Z_{h}^{\perp}$. Moving on, we can
solve for ${u}^{n}\in V_{h}$ via

(15)
$a({u}^{n},{v}_{h})=({{\nabla\cdot}\,}\hat{u}^{n-1},{{\nabla\cdot}\,}{v}_{h})-\sigma\,a(\hat{u}^{n-1},{v}_{h})\quad\forall\;{v}_{h}\in
V_{h}.$

By induction, $a({u}^{n},{v}_{h})=0$ for all ${v}_{h}\in Z_{h}$, so that
${u}^{n}\in Z_{h}^{\perp}$ for all $n>0$.

To start the process, we can solve for ${u}^{0}\in V_{h}$ via

(16)
$a({u}^{0},{v})=({{\nabla\cdot}\,}{w},{{\nabla\cdot}\,}{v})\quad\forall\;{v}\in
V_{h},$

where we must pick some ${w}$ where ${{\nabla\cdot}\,}{w}\neq 0$. In practice,
we chose

(17) ${w}=(\sin(\omega x),\cos(\omega y)),$

for $\omega\in\mathbb{Z}$, although other initializations worked as well. We
then set

(18) $\hat{u}^{0}=(\|{u}^{0}\|_{V})^{-1}{u}^{0}.$

Note that $a({u}^{0},{v})=0$ for all ${v}\in Z_{h}$, which means that
${u}^{0}\in Z_{h}^{\perp}$. The division in Eq. 18 provides a natural check
that ${u}^{0}\neq{\bf 0}$.

###  2.3 Eigenvalue bounds

The eigenvalues $\lambda_{i}$ are bounded above since

$\|{{\nabla\cdot}\,}{v}_{h}\|_{L^{2}}^{2}\leq
Ca({v}_{h},{v}_{h})\quad\forall\;{v}_{h}\in V.$

This is obvious for the gradient form Eq. 5 with $C\leq d$, where $d$ is the
dimension of $\Omega$. For the ${\boldsymbol{\varepsilon}}$ form Eq. 10, it
follows by the Korn inequality [9, section 11.2]. A good estimate of $C$ can
be obtained computationally since it is the largest eigenvalue of Eq. 11. In
all computational examples here, using the gradient form Eq. 5, $C=1$,
apparently due to the homogeneous Dirichlet boundary conditions in $V$. Using
the ${\boldsymbol{\varepsilon}}$ form Eq. 10, which is equivalent to Eq. 5 by
Korn’s inequality, would change the computational results here by at most a
constant factor. When homogeneous Dirichlet boundary conditions are enforced,
it is known [46, (13.12)] that the forms Eq. 10 and Eq. 5 produce identical
values for divergence-free functions.

Thus we can shift by a constant $\sigma$ independent of $h$ to ensure that
$\lambda^{\sigma}_{i}<0$ for all $i$. The algorithm Eq. 14 is the power method
for $A^{-1}B-\sigma I$, and it will generically converge to the eigenvector
associated with the most negative eigenvalue, provided
$\sigma>{\textstyle{1\over 2}}\lambda_{N}$, and so it will generically
converge to $\lambda_{1}^{\sigma}$. Taking $\sigma=0$ will converge to the
largest eigenvalue $\lambda_{N}$, since $0<\lambda_{1}$. We found that taking
$\sigma=0.6$ gave acceptable convergence.

###  2.4 Effect of round-off error

It is essential to start the iteration Eq. 14 with $\hat{u}^{0}\in
Z_{h}^{\perp}$. To see what goes wrong otherwise, write the algorithm in Eq.
14 as

${u}^{k+1}=c_{k}(A^{-1}B-\sigma I){u}^{k},$

where $c_{k}=1/\|{u}^{k}\|_{V}$. Suppose that ${u}^{0}={v}^{0}+{w}^{0}$ where
${v}^{0}\in Z_{h}^{\perp}$ and ${w}^{0}\in Z_{h}$. Then
${u}^{k}={v}^{k}+{w}^{k}$ where ${v}^{k}\in Z_{h}^{\perp}$ and ${w}^{k}\in
Z_{h}$, and

${v}^{k+1}=c_{k}(A^{-1}B-\sigma I){v}^{k},\qquad{w}^{k+1}=c_{k}(-\sigma
I){w}^{k}.$

The reason is that $A^{-1}B$ maps $Z_{h}^{\perp}$ into itself and $Z_{h}$ to
the zero vector. Define $C_{k}=\prod_{i=0}^{k-1}c_{i}$. Then

${v}^{k}=C_{k}(A^{-1}B-\sigma I)^{k}{v}^{0},\qquad{w}^{k}=C_{k}(-\sigma
I)^{k}{w}^{0}.$

In seeking the smallest eigenvalue $\lambda_{1}$, we will need to take
$\sigma$ as large as at least half of the largest eigenvalue $\lambda_{N}$.
Thus ${w}^{k}$ can become significant even if ${w}^{0}$ is on the order of
round-off error. Once ${w}^{k}$ becomes dominant, the Rayleigh quotient in Eq.
14 defining $\lambda^{k}$ will go to zero. Thus it is necessary to monitor the
projection ${z}^{k}$ of ${u}^{k}$ onto $Z_{h}$ and to project ${u}^{k}$ onto
$Z_{h}^{\perp}$ (by simply subtracting ${z}^{k}$) when the divergence-zero
component is too large. This may need to be done for $k=0$ as well.

The projection ${z}^{k}\in Z_{h}$ satisfies

$a({z}^{k},{v})=a({u}^{k},{v})\;\forall\;{v}\in Z_{h}.$

The projection ${z}^{k}$ can be computed via the iterated penalty method. Some
care is required in using the iterated penalty method to do this, since it
will be slow to converge if the inf-sup constant is small. But this appears to
work well in practice, albeit with a large number of iterations required when
the inf-sup constant is smaller.

###  2.5 Full algorithm

We now summarize the full algorithm. First we compute $\hat{u}^{0}$ via Eq. 16
and Eq. 18, where ${w}$ is given in Eq. 17. Next, we solve for ${u}^{n}$ via
Eq. 15, $n\geq 1$. Then we project ${u}^{n}$ onto $Z_{h}$ using the iterated
penalty method [9, §13.1]. To compute the projection of ${u}^{n}$ onto
$Z_{h}$, we solve for ${z}^{\ell},\;{w}^{\ell}\in V_{h}$ such that

(19)
$\begin{split}&a({z}_{\ell},{v})+\rho({{\nabla\cdot}\,}{z}_{\ell},{{\nabla\cdot}\,}{v})=a({u}^{n},{v})-({{\nabla\cdot}\,}{w}_{\ell},{{\nabla\cdot}\,}{v})\quad\forall\;{v}\in
V_{h}\\\ &{w}_{\ell+1}={w}_{\ell+1}+\rho{z}_{\ell},\end{split}$

where we start with ${w}_{0}=0$. The parameter $\rho$ is the penalty parameter
enforcing the incompressibility constraint; in our computations we set
$\rho=10^{4}$. We terminate the iteration on $\ell$ when

(20) $\|{{\nabla\cdot}\,}{z}_{\ell}\|_{L^{2}(\Omega)}\leq\tau,$

where we picked $\tau=10^{-14}$ in our computations. Then
${z}_{\ell}\approx\Pi_{Z_{h}}{u}^{n}$, and if

$\|\nabla{z}_{\ell}\|_{L^{2}(\Omega)}\geq\zeta\|\nabla{u}^{n}\|_{L^{2}(\Omega)},$

we update ${u}^{n}={u}^{n}-{z}_{\ell}.$ In our computations, we picked
$\zeta=10^{-12}$. Once ${u}^{n}$ has been computed, we define

(21)
$\begin{split}&\lambda^{n+1}=\frac{({{\nabla\cdot}\,}{u}^{n},{{\nabla\cdot}\,}{u}^{n})}{a({u}^{n},{u}^{n})},\qquad\hat{u}^{n}=(\|{u}^{n}\|_{V})^{-1}{u}^{n}.\end{split}$

This iteration is continued while $|\lambda^{n+1}-\lambda^{n}|>\epsilon$,
where $\epsilon$ is a pre-specified tolerance, with $\epsilon=10^{-8}$ in our
computations.

###  2.6 2D tests

Table 1: Mesh restrictions for $V_{h}$ constructed with continuous piecewise polynomials of degree $k$, and $\Pi_{h}={{\nabla\cdot}\,}V_{h}$. Here $d$ is the dimension of $\Omega$. In the lowest-order case $k=1$ in two dimensions, the velocity approximation is optimal order despite the fact that the inf-sup constant deteriorates. $d$ | $k$ | inf-sup | mesh restrictions
---|---|---|---
2 |  1 | no, but | optimal velocity approximation on Malkus splits [33, 40]
2 |  1 | yes | Powell–Sabin splits [13, 20, 53]
2 |  2 | yes | some crossed triangles required [39] or Alfeld splits [22, 17]
2 |  3 | yes | under certain conditions [23]
2 | $\geq 4$ | yes |  $p$-robust when no nearly singular vertices [48, 49, 24, 1]
3 | $\geq 1$ | yes | Worsey–Farin splits [13]
3 | $\geq 3$ | yes | Alfeld splits [52, 22, 17]
3 | $\geq 6$ | yes | only one family of meshes known [54]
Table 2: Computation of inf-sup constants on Malkus splits in two dimensions. The shift parameter in Eq. 13 was set to $\sigma=0.6$. The mesh size $N$ refers to the $N\times N$ mesh of quadrilaterals before subdivision. The parameter $\omega$ determines the initial eigenvector approximation defined in Eq. 17. Iterations were continued until the change in eigenvalue was less than $10^{-7}$. $k$ | $N$ | $\omega$ | iterations | inf-sup $\lambda_{1}$ | mesh type | restarts
---|---|---|---|---|---|---
1 | 5 | 5 | 47 | $4.08\times 10^{-1}$ | Malkus splits | 0
1 | 5 | 10 | 40 | $4.08\times 10^{-1}$ | Malkus splits | 0
1 | 5 | 100 | 51 | $4.08\times 10^{-1}$ | Malkus splits | 1
1 | 10 | 10 | 258 | $1.13\times 10^{-2}$ | Malkus splits | 1
1 | 20 | 10 | 642 | $2.98\times 10^{-3}$ | Malkus splits | 1
2 | 10 | 10 | 193 | $1.49\times 10^{-1}$ | Malkus splits | 13
2 | 20 | 10 | 189 | $1.48\times 10^{-1}$ | Malkus splits | 14
2 | 40 | 10 | 187 | $1.48\times 10^{-1}$ | Malkus splits | 19

Figure 2: Two types of regular meshes in two dimensions. The Type II mesh is also known as the Malkus split. Table 3: Computation of inf-sup constants on Type I meshes (see Fig. 2) in two dimensions. The shift parameter in Eq. 13 was set to $\sigma=0.6$. The mesh size $N$ refers to the $N\times N$ mesh of quadrilaterals before subdivision. The parameter $\omega$ determines the initial eigenvector approximation defined in Eq. 17. Iterations were continued until the change in eigenvalue was less than $10^{-7}$ for $k=3$ and $k=4$, but we had to continue iterations until the change in eigenvalue was less than $10^{-10}$ for $k=2$ to obtain a reliable second digit. $k$ | $N$ | $\omega$ | iterations | inf-sup $\lambda_{1}$ | mesh type | restarts
---|---|---|---|---|---|---
2 | 8 | 10 | 13577 | $1.60\times 10^{-3}$ | Type I | 29
2 | 16 | 10 | 50690 | $4.08\times 10^{-4}$ | Type I | 129
2 | 32 | 10 | 19851 | $1.07\times 10^{-4}$ | Type I | 19397
2 | 64 | 10 | 54650 | $2.75\times 10^{-5}$ | Type I | 54649
3 | 3 | 10 | 790 | $8.46\times 10^{-3}$ | Type I | 11
3 | 5 | 10 | 768 | $3.52\times 10^{-3}$ | Type I | 8
3 | 10 | 10 | 1858 | $9.50\times 10^{-4}$ | Type I | 15
4 | 5 | 10 | 585 | $2.59\times 10^{-2}$ | Type I | 12
4 | 10 | 10 | 789 | $2.60\times 10^{-2}$ | Type I | 17
4 | 20 | 10 | 1265 | $2.60\times 10^{-2}$ | Type I | 28

To verify the algorithm, we compare with known results. We summarize some
known results in two and three dimensions in Table 1. In Table 2, we give
results for Lagrange elements of degree $k=1,2$ in two dimensions on Malkus
splits (Type II meshes in the nomenclature of [28] as shown in Fig. 2), for
various meshsizes. A Malkus split mesh starts with a quadrilateral subdivision
and creates a triangulation by subdividing each quadrilateral by adding the
two diagonals connecting opposite vertices. On Malkus splits, for $k\geq 2$,
the inf-sup constant is bounded and converges relatively rapidly in our tests,
but for $k=1$, the inf-sup constant goes to zero with a rate equal to one over
the meshsize parameter [39]. The column “iterations” lists the number of
iterations of the power method required to achieve a change in $\lambda$ less
than $10^{-8}$. The column “restarts” lists the number of times the projection
${z}_{\ell}$ was subtracted from ${u}^{n}$.

In Table 3, we give results for Lagrange elements of degree $k=2,3,4$ in two
dimensions on Type I meshes, that is, meshes consisting of 45-degree right
triangles (Fig. 2). The results match the theory indicated in Table 1: the
inf-sup constant degenerates for $k=2,3$, but does not for $k=4$. There is
renewed interest in the low order cases due to the emergence of the grad-div
penalized Taylor–Hood method [10, 31, 32]. The approximations from these
methods tend to divergence free functions as the penalty is increased, and the
low-order Taylor–Hood methods are widely used. It has been known that
divergence-free quadratics on Type I meshes have reduced approximation order
(this is equivalent to reduced approximation for $C^{1}$ piecewise cubic
approximation for scalar functions [12]). But so far we know of no estimates
for controlled degeneration of the inf-sup constants on general meshes.

###  2.7 3D tests

Table 4: Computation of inf-sup constants for various polynomial degrees $k$ in three dimensions on the Freudenthal meshes depicted in Fig. 1. The shift parameter in Eq. 13 was set to $\sigma=0.6$. The mesh size $N$ refers to the $N\times N\times N$ mesh of cubes before subdivision. Iterations were continued until the change in eigenvalue was less than $10^{-7}$ for $k=4$ and $k=5$, but we had to continue iterations until the change in eigenvalue was less than $10^{-8}$ for $k=3$.  degree $k$ | $N$ | iterations | inf-sup $\lambda_{1}$ | restarts
---|---|---|---|---
3 | 2 | 1026 | $5.75\times 10^{-4}$ | 4
3 | 3 | 3400 | $4.26\times 10^{-4}$ | 14
3 | 4 | 4445 | $2.93\times 10^{-4}$ | 29
3 | 5 | 6044 | $2.11\times 10^{-4}$ | 2048
3 | 6 | 7732 | $1.59\times 10^{-4}$ | 4863
4 | 2 | 802 | $4.28\times 10^{-3}$ | 7
4 | 3 | 2067 | $4.13\times 10^{-3}$ | 16
4 | 4 | 1961 | $4.22\times 10^{-3}$ | 16
4 | 5 | 1325 | $4.27\times 10^{-3}$ | 12
4 | 6 | 1380 | $4.26\times 10^{-3}$ | 12
5 | 2 | 1004 | $6.46\times 10^{-3}$ | 12
5 | 3 | 1003 | $6.59\times 10^{-3}$ | 13
5 | 4 | 877 | $6.34\times 10^{-3}$ | 12
5 | 5 | 938 | $6.26\times 10^{-3}$ | 13
5 | 6 | 978 | $6.23\times 10^{-3}$ | 13

Mesh restrictions in three dimensions for Scott–Vogelius elements are not
fully understood. Our computations for Freudenthal meshes are summarized in
Table 4. They reveal a surprising fact: on this mesh family, the inf-sup
constant is bounded for degrees $k\geq 4$, whereas Zhang [54] was able to
prove inf-sup stability only for $k\geq 6$. Neilan [36, Proposition 6.5]
extends this to more general meshes that satisfy a special condition, which he
states as a conjecture, but still for $k\geq 6$. The results in Table 4 lead
us to make 1.

##  3 Space decompositions for the kernel of the divergence

Our second conjecture relates to multigrid solvers for Eq. 2. In the nearly
incompressible regime, the equation becomes nearly singular, and standard
multigrid methods break down. A key breakthrough for such problems was made by
Schöberl [43], who devised conditions on the relaxation and prolongation
operators that guarantee that a multigrid method is parameter-robust in the
nearly-singular regime.

###  3.1 Background

The key condition for the relaxation is best stated in terms of space
decompositions and subspace corrections [51]. The multigrid relaxation method
we employ will be induced by a space decomposition

(22) $V_{h}=\sum_{i=1}^{J}V_{i}$

where an equation for an approximation to the error is solved over each
subspace $V_{i}$ and the updates combined additively or multiplicatively. Each
$V_{i}$ is assumed small enough so that direct solvers can be afforded. For
example, if $V_{h}=\mathrm{span}\\{\phi_{1},\phi_{2},\dots,\phi_{N}\\}$, and
each $V_{i}$ is chosen to be $V_{i}=\mathrm{span}\\{\phi_{i}\\}$, then
combining the updates additively would yield the Jacobi relaxation, while
combining multiplicatively would yield the Gauss–Seidel relaxation. In a
domain decomposition approach, each $V_{i}$ could be taken as the functions
supported on one subdomain of a given parallel decomposition of the mesh,
combined with a suitable global coarse space. Another important example is to
define each $V_{i}$ as the functions with support on the patch of cells
surrounding each vertex, the so-called _vertex-star_ space decomposition: this
arises in $k$-robust preconditioners for symmetric and coercive problems [38,
44], as the Arnold–Falk–Winther (AFW) relaxation for $H(\mathrm{div})$ and
$H(\mathrm{curl})$ [3], and in Reynolds-robust preconditioners for the
Navier–Stokes equations [6, 16].

Let us consider an abstract nearly singular problem, following Lee et al.
[29]. For $\varepsilon>0$, consider the finite-dimensional linear variational
problem: find $u\in V_{h}$ such that

(23) $a_{0}(u,v)+\varepsilon^{-1}a_{1}(u,v)=(f,v)\text{ for all }v\in V_{h},$

where $a_{0}$ is symmetric and coercive, and $a_{1}$ is symmetric but only
positive semi-definite. In the context of Eq. 2, $a_{1}(u,v)=(\nabla\cdot
u,\nabla\cdot v)_{L^{2}}$. Define the kernel

(24) $\mathcal{N}=\\{u\in V_{h}:a_{1}(u,v)=0\text{ for all }v\in V_{h}\\}.$

In our context, these are divergence-free functions in the finite element
space. A key condition that must be satisfied by the space decomposition for
parameter-robustness in $\varepsilon$ is [29, Assumption (A1)]:

(25) $\mathcal{N}=\sum_{i=1}^{J}\left(\mathcal{N}\cap V_{i}\right).$

In other words, when challenged with a divergence-free function, it must be
possible to decompose this as the sum of functions with the $i^{\mathrm{th}}$
summand drawn from the divergence-free functions in $V_{i}$***This
decomposition must also be stable, but we shall not elaborate here.. The space
decomposition that induces Jacobi or Gauss–Seidel relaxation does not
generally satisfy Eq. 25, because the intersection $\mathcal{N}\cap
V_{i}=\\{0\\}$ for typical finite element methods where the standard bases
used are not divergence-free.

###  3.2 Space decompositions from de Rham complexes

One way to devise space decompositions that satisfy Eq. 25 is by inspecting
discrete subcomplexes of a suitable underlying Hilbert complex. For
concreteness, consider $\Omega\subset\mathbb{R}^{2}$, with $\Omega$ simply
connected. The Stokes complex [27] is given by

(26)
$\displaystyle\mathbb{R}\xrightarrow{\operatorname{id}}H^{2}(\Omega)\xrightarrow{{\rm
curl}\,}[H^{1}(\Omega)]^{2}\xrightarrow{\operatorname{div}}L^{2}(\Omega)\xrightarrow{\operatorname{null}}0.$

This complex is discretized with a discrete subcomplex

(27) $\mathbb{R}\xrightarrow{\operatorname{id}}\Sigma_{h}\xrightarrow{{\rm
curl}\,}V_{h}\xrightarrow{\operatorname{div}}Q_{h}\xrightarrow{\operatorname{null}}0,$

where $\Sigma_{h}\subset H^{2}(\Omega)$, $V_{h}\subset[H^{1}(\Omega)]^{2}$,
and $Q_{h}\subset L^{2}(\Omega)$. These complexes have the property that the
kernel of an operator ${\rm curl}\,,\operatorname{div}$, or
$\operatorname{null}$ is a subspace of the range of the preceeding operator,
e.g. that $\mathrm{ker}(\operatorname{div})\subset\mathrm{range}({\rm
curl}\,)$ [2]. The complex is called _exact_ if the kernel of an operator is
precisely the range of the preceeding operator. The Stokes complex is exact if
the domain $\Omega$ is simply connected, and this property is inherited by Eq.
27 if constructed appropriately (with bounded cochain projections).

In our context, these complexes are useful because they offer a crisp
characterization of $\mathcal{N}=\mathrm{ker}(\operatorname{div})$. Let
$\Sigma_{h}=\mathrm{span}\\{\psi_{1},\dots,\psi_{M}\\}$ and
$V_{h}=\mathrm{span}\\{\phi_{1},\dots,\phi_{N}\\}$. If $u_{h}\in\mathcal{N}$,
then $u_{h}={\rm curl}\,\Psi_{h}$ for $\Psi_{h}\in\Sigma_{h}$ by exactness of
Eq. 27. Expanding $\Psi_{h}$ in terms of its basis functions

(28) $\Psi_{h}=\sum_{i=1}^{M}c_{i}\psi_{i}$

yields an expression for $u_{h}$

(29) $u_{h}=\sum_{i=1}^{M}c_{i}{\rm curl}\,\psi_{i}.$

Suppose the space decomposition Eq. 22 is chosen as $J=N+M$ with

(30) $V_{i}=\begin{cases}\mathrm{span}({\rm curl}\,\psi_{i})&i=1,\dots,M\\\
\mathrm{span}(\phi_{i-M})&i=M+1,\dots,M+N.\end{cases}$

This decomposition would satisfy Eq. 25, with $V_{i}\cap\mathcal{N}=V_{i}$ for
$i\leq M$ or $V_{i}\cap\mathcal{N}=\\{0\\}$ for $i>M$†††This space
decomposition relates to that proposed by Hiptmair [25]. As written it is not
robust in $k$, but (experimentally) can be made so by applying suitable block
Jacobi methods in $V_{h}$ and $\Sigma_{h}$..

Solving the problem Eq. 23 over $\mathrm{span}\\{{\rm curl}\,\psi_{i}\\}$
reduces to solving a problem in a locally-supported subspace of $\Sigma_{h}$,
as described in Hiptmair [25] and Hiptmair–Xu [26] for the $L^{2}$ de Rham
complex. However, for the Stokes complex in three dimensions, explicit
constructions of $\Sigma_{h}$ are poorly understood. The alternative approach
(which we refer to as Pavarino–Arnold–Falk–Winther, PAFW) is to construct the
space decomposition using only knowledge of the supports of the basis
functions $\psi_{i}$. For example, if $\Sigma_{h}$ exists with the property
that each basis function $\psi_{i}$ is a polynomial of degree $k+1$ supported
on a certain region of the domain $\mathrm{supp}(\psi_{i})$, then we may
choose $\\{V_{j}\\}$ so that each $V_{j}$ captures all polynomials of degree
$k$ on one $\mathrm{supp}(\psi_{i})$. In this way, $\forall
i=1,\dots,M\,\exists\,j\text{ s.t. }{\rm curl}\,\psi_{i}\in V_{j}$, and Eq. 25
will also be satisfied, without needing to explicitly construct $\Sigma_{h}$
or $\\{\psi_{i}\\}_{i=1}^{M}$. Only knowledge (or conjecture) of the supports
is required.

To illustrate the PAFW construction, consider the case of the Scott–Vogelius
element pair with $k\geq 4$ and $d=2$. This element constructs $V_{h}$ with
$[\mathrm{CG}_{k}]^{d}$ and $Q_{h}$ with $\mathrm{DG}_{k-1}$ on a mesh
$\mathcal{T}_{h}$. The corresponding $C^{1}$-conforming finite element for
$\Sigma_{h}$ is the Morgan–Scott element of degree $k+1$ [34], which is not
implemented in general purpose finite element software, preventing the
application of Eq. 30. The Morgan–Scott element employs (among others) degrees
of freedom at the vertices of the mesh, and the basis function associated with
such a degree of freedom will have support over the patch of cells sharing
that vertex. This suggests the following vertex-star space decomposition,
which we now describe. Let ${\nu}_{1},\dots,{\nu}_{J}$ be the vertices of
$\mathcal{T}_{h}$. The subspaces are given by

(31) $V_{j}=\\{v_{h}\in
V_{h}:\mathrm{supp}(v_{h})\subset\mathrm{star}({\nu}_{j})\\},$

where the $\mathrm{star}$ operation of a simplex $p$ returns the union of all
simplices containing $p$ as a subsimplex [35, §2]. In other words, when
applied to a vertex ${\nu}_{j}$, the star is the union of the cells and edges
(and faces in $d=3$) containing ${\nu}_{j}$, as well as ${\nu}_{j}$ itself. In
[34], it is proven that the Morgan–Scott basis functions are each supported in
$\mathrm{star}({\nu}_{j})$ for some $j$, provided $k\geq 4$. Thus Eq. 31 also
satisfies the kernel decomposition Eq. 25. Numerical experiments indicate that
Eq. 31 does indeed yield $\varepsilon$-robust convergence for $k\geq 4$, as
expected from the Morgan–Scott theory, while it does not for $k<4$ [14, §4.2].

There are other examples where the vertex-star space decomposition does not
yield $\varepsilon$-robust convergence. For $k=1$ on Malkus splits, the basis
with smallest possible support has support that is larger than
$\mathrm{star}({\nu}_{j})$, although it is in
$\mathrm{star}(\mathrm{closure}(\mathrm{star}({\nu}_{j})))$ [47], where the
$\mathrm{star}$ of a set of simplices is the union of their stars, and closure
is defined in [14].

###  3.3 Second conjecture

The preceding discussion indicates that one may experimentally investigate
whether $\Sigma_{h}$ permits a basis with support captured by the stars of the
vertices, by applying the multigrid algorithm and observing whether the
convergence is $\varepsilon$-robust or not.

The continuous Stokes complexes in 3D [15] is given by

(32) $\mathbb{R}\xrightarrow{\operatorname{id}}H^{2}(\Omega)\xrightarrow{{\rm
grad}\,}H^{1}({\rm curl}\,,\Omega)\xrightarrow{{\rm
curl}\,}[H^{1}(\Omega)]^{3}\xrightarrow{\operatorname{div}}L^{2}(\Omega)\xrightarrow{\operatorname{null}}0.$

Here

(33) $H^{1}({\rm curl}\,,\Omega)\coloneqq\\{v\in[H^{1}(\Omega)]^{3}:{\rm
curl}\,v\in[H^{1}(\Omega)]^{3}\\}.$

This complex is discretized with a discrete subcomplex

(34) $\mathbb{R}\xrightarrow{\operatorname{id}}S_{h}\xrightarrow{{\rm
grad}\,}\Sigma_{h}\xrightarrow{{\rm
curl}\,}V_{h}\xrightarrow{\operatorname{div}}\Pi_{h}\xrightarrow{\operatorname{null}}0.$

The space $S_{h}$ consists of $C^{1}$ piecewise polynomials of degree $k+2$.
The space $\Sigma_{h}$ consists of continuous vector-valued piecewise
polynomials of degree $k+1$ which have a continuous curl. Note that
$(S_{h})^{3}\subset\Sigma_{h}$.

The potential space $\Sigma_{h}$ is known on special meshes (Alfeld or
Worsey–Farin splits) in three dimensions [17, 21, 8]. We are not aware of any
results that characterize $\Sigma_{h}$ on Freudenthal meshes.

Based on the numerical experiments we will report, we make the following
conjecture:

######  Conjecture 3.1.

Let $V_{h}\subset[H^{1}(\Omega)]^{3}$ be constructed with continuous Lagrange
elements of degree $k$, and let $\Sigma_{h}\subset H^{1}({\rm curl}\,,\Omega)$
be the space preceeding $V_{h}$ in a subcomplex of the three-dimensional
Stokes complex. For $k\geq 5$, there is a local basis for $\Sigma_{h}$,
supported on the stars of vertices, on the Freudenthal meshes depicted in Fig.
1.

###  3.4 Numerical experiments

We provide numerical evidence for 3.1. We solve the problem: find $u\in
V_{h}\subset[{H}^{1}_{0}(\Omega)]^{3}$ such that:

(35) $(\nabla u,\nabla
v)+\gamma({{\nabla\cdot}\,}u,{{\nabla\cdot}\,}v)=(1,v)\quad\forall v\in
V_{h}.$

The problem becomes nearly singular as $\gamma\to\infty$. Here $V_{h}$ is
constructed with continuous Lagrange elements of degree $k$. We employ the
solver denoted in Fig. 3.

Krylov solver: CGTwo-grid V-cycleRelaxation: damped additive vertex-star
iterationCoarse grid: Cholesky factorization Figure 3: Solver diagram for Eq.
35.

The relaxation on each level is one application of damped Richardson iteration
preconditioned by the additive Schwarz method with vertex-star subspaces
chosen as Eq. 31. The damping factor is set to $1/M$ where $M$ is the maximum
number of patches that a given vertex is contained in. Here $M=3$ for $d=2$
and $M=4$ for $d=3$. The coarse grid solve was computed with CHOLMOD [11] via
PETSc [5]. The code was run in serial.

Table 5: Iteration counts for Eq. 35 on a uniform refinement of a Type I $4\times 4$ mesh. The iteration counts are $\gamma$-robust for $k\geq 4$, but not for $k<4$. $k\backslash\gamma$ | $10^{0}$ | $10^{1}$ | $10^{2}$ | $10^{3}$ | $10^{4}$ | $10^{5}$
---|---|---|---|---|---|---
2 | 11 | 17 | 35 | 54 | 87 | 100
3 | 10 | 13 | 23 | 49 | 62 | 83
4 | 9 | 10 | 13 | 13 | 13 | 12
5 | 9 | 10 | 11 | 12 | 11 | 10

###  3.5 Two dimensions

The convergence results for the solver depicted in Figure 3 are shown in Table
5, as a function of $k$ and $\gamma$. The Krylov method was terminated when
the unpreconditioned residual was reduced by 8 orders of magnitude. The coarse
mesh was a Type I mesh [28] of size $4\times 4$. The main observation is that
$k\geq 4$ is required for $\gamma$-robust performance. This matches well with
the existing theory described above: the key property that changes between
these two cases is the existence of the Morgan–Scott element for $\Sigma_{h}$
with a basis supported on vertex-stars for $k\geq 4$.

###  3.6 Three dimensions

In three dimensions we again choose $V_{h}=[\mathrm{CG}_{k}]^{d}$. We employ
an analog of a Type I coarse mesh of size $4\times 4\times 4$, with one
refinement. We show iteration counts for the solution of the problem of Eq. 35
in Table 6. As in two dimensions, the solver was terminated when the
unpreconditioned residual was reduced by 8 orders of magnitude. These results
lead us to 3.1.

Table 6: Iteration counts for Eq. 35 on 2 uniform refinements of a Type I $4\times 4\times 4$ mesh. The iteration counts are $\gamma$-robust for $k\geq 5$, but not for $k<5$. $k\backslash\gamma$ | $10^{0}$ | $10^{1}$ | $10^{2}$ | $10^{3}$ | $10^{4}$ | $10^{5}$
---|---|---|---|---|---|---
2 | 14 | 23 | 52 | 118 | 300 | 849
3 | 11 | 16 | 29 | 66 | 173 | 458
4 | 11 | 13 | 19 | 26 | 54 | 110
5 | 11 | 12 | 16 | 19 | 20 | 19
6 | 10 | 11 | 14 | 15 | 16 | 15
7 | 10 | 11 | 13 | 14 | 14 | 13

Note that our 1 is that the global inf-sup condition holds on this mesh family
for $k\geq 4$. Thus there is a gap for $k=4$: we conjecture that the global
inf-sup condition holds, but that there is not a basis for the potential space
supported on vertex-stars.

##  4 Conclusions

Our computational experiments have suggested two conjectures regarding the
three-dimensional Stokes complex on Freudenthal meshes. We conjecture that the
inf-sup condition holds for a velocity space $V_{h}$ consisting of continuous
piecewise polynomials of degree $k$ and pressure space
${{\nabla\cdot}\,}V_{h}$ for $k\geq 4$, and that there is a local basis for
the associated potential space $\Sigma_{h}$ for $k\geq 5$.

##  5 Acknowledgements

The authors would like to thank F. Wechsung, M. Neilan, and C. Parker for
useful discussions.

## References

  * [1] M. Ainsworth and C. Parker, Unlocking the secrets of locking: finite element analysis in planar linear elasticity, Computer Methods in Applied Mechanics and Engineering, 395 (2022), p. 115034, https://doi.org/10.1016/j.cma.2022.115034.
  * [2] D. N. Arnold, Finite element exterior calculus, Society for Industrial and Applied Mathematics, 2018, https://doi.org/10.1137/1.9781611975543.
  * [3] D. N. Arnold, R. S. Falk, and R. Winther, Multigrid in $H(\mathrm{div})$ and $H(\mathrm{curl})$, Numerische Mathematik, 85 (2000), pp. 197–217, https://doi.org/10.1007/pl00005386.
  * [4] D. N. Arnold and M. E. Rognes, Stability of Lagrange elements for the mixed Laplacian, Calcolo, 46 (2009), pp. 245–260, https://doi.org/10.1007/s10092-009-0009-6, http://dx.doi.org/10.1007/s10092-009-0009-6.
  * [5] S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan, J. Sarich, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, and J. Zhang, PETSc/TAO users manual, Tech. Report ANL-21/39 - Revision 3.17, Argonne National Laboratory, 2022\.
  * [6] M. Benzi and M. A. Olshanskii, An augmented Lagrangian-based approach to the Oseen problem, SIAM Journal on Scientific Computing, 28 (2006), pp. 2095–2113, https://doi.org/10.1137/050646421.
  * [7] J. Bey, Simplicial grid refinement: On Freudenthal’s Algorithm and the optimal number of congruence classes, Numerische Mathematik, 85 (2000), pp. 1–29, https://doi.org/10.1007/s002110050475.
  * [8] D. Boffi, S. Gong, J. Guzmán, and M. Neilan, Convergence of Lagrange finite element methods for Maxwell eigenvalue problem in 3D, 2022, https://arxiv.org/abs/2204.10876.
  * [9] S. C. Brenner and L. R. Scott, The Mathematical Theory of Finite Element Methods, Springer-Verlag, third ed., 2008, https://doi.org/10.1007/978-0-387-75934-0.
  * [10] M. A. Case, V. J. Ervin, A. Linke, and L. G. Rebholz, A connection between Scott–Vogelius and grad-div stabilized Taylor–Hood FE approximations of the Navier–Stokes equations, SIAM Journal on Numerical Analysis, 49 (2011), pp. 1461–1481, https://doi.org/10.1137/100794250.
  * [11] Y. Chen, T. A. Davis, W. W. Hager, and S. Rajamanickam, Algorithm 887: CHOLMOD, supernodal sparse Cholesky factorization and update/downdate, ACM Transactions on Mathematical Software, 35 (2008), pp. 1–14, https://doi.org/10.1145/1391989.1391995.
  * [12] C. de Boor and K. Höllig, Approximation order from bivariate $c^{1}$-cubics: a counterexample, Proceedings of the American Mathematical Society, 87 (1983), pp. 649–655, https://doi.org/10.1090/S0002-9939-1983-0687634-4.
  * [13] M. Fabien, J. Guzmán, M. Neilan, and A. Zytoon, Low-order divergence-free approximations for the Stokes problem on Worsey–Farin and Powell–Sabin splits, Computer Methods in Applied Mechanics and Engineering, 390 (2022), p. 114444, https://doi.org/10.1016/j.cma.2021.114444.
  * [14] P. E. Farrell, M. G. Knepley, L. Mitchell, and F. Wechsung, PCPATCH: Software for the topological construction of multigrid relaxation methods., ACM Transactions on Mathematical Software, 47 (2021), pp. 1–22, https://doi.org/10.1145/3445791.
  * [15] P. E. Farrell, L. Mitchell, L. R. Scott, and F. Wechsung, A Reynolds-robust preconditioner for the Scott–Vogelius discretization of the stationary incompressible Navier–Stokes equations, SMAI Journal of Computational Mathematics, 7 (2021), pp. 75–96, https://doi.org/10.5802/smai-jcm.72.
  * [16] P. E. Farrell, L. Mitchell, and F. Wechsung, An augmented Lagrangian preconditioner for the 3D stationary incompressible Navier–Stokes equations at high Reynolds number, SIAM Journal on Scientific Computing, 41 (2019), pp. A3073–A3096, https://doi.org/10.1137/18M1219370.
  * [17] G. Fu, J. Guzmán, and M. Neilan, Exact smooth piecewise polynomial sequences on Alfeld splits, Mathematics of Computation, 89 (2020), pp. 1059–1091, https://doi.org/10.1090/mcom/3520.
  * [18] D. Gallistl, Rayleigh–Ritz approximation of the inf-sup constant for the divergence, Mathematics of Computation, 88 (2019), pp. 73–89, https://doi.org/10.1090/mcom/3327.
  * [19] G. H. Golub and Q. Ye, An inverse free preconditioned Krylov subspace method for symmetric generalized eigenvalue problems, SIAM Journal on Scientific Computing, 24 (2002), pp. 312–334, https://doi.org/10.1137/S1064827500382579.
  * [20] J. Guzmán, A. Lischke, and M. Neilan, Exact sequences on Powell–Sabin splits, Calcolo, 57 (2020), pp. 1–25, https://doi.org/10.1007/s10092-020-00361-x.
  * [21] J. Guzmán, A. Lischke, and M. Neilan, Exact sequences on Worsey–Farin splits, Mathematics of Computation, 91 (2022), pp. 2571–2608, https://doi.org/10.1090/mcom/3746.
  * [22] J. Guzmán and M. Neilan, Inf-sup stable finite elements on barycentric refinements producing divergence-free approximations in arbitrary dimensions, SIAM Journal on Numerical Analysis, 56 (2018), pp. 2826–2844, https://doi.org/10.1137/17M1153467.
  * [23] J. Guzmán and L. R. Scott, Cubic Lagrange elements satisfying exact incompressibility, SMAI Journal of Computational Mathematics, 4 (2018), pp. 345–374, https://doi.org/10.5802/smai-jcm.38.
  * [24] J. Guzmán and L. R. Scott, The Scott–Vogelius finite elements revisited, Mathematics of Computation, 88 (2019), pp. 515–529, https://doi.org/10.1090/mcom/3346.
  * [25] R. Hiptmair, Multigrid method for Maxwell's equations, SIAM Journal on Numerical Analysis, 36 (1998), pp. 204–225, https://doi.org/10.1137/s0036142997326203.
  * [26] R. Hiptmair and J. Xu, Nodal auxiliary space preconditioning in $H(\mathrm{curl})$ and $H(\mathrm{div})$ spaces, SIAM Journal on Numerical Analysis, 45 (2007), pp. 2483–2509, https://doi.org/10.1137/060660588.
  * [27] V. John, A. Linke, C. Merdon, M. Neilan, and L. G. Rebholz, On the divergence constraint in mixed finite element methods for incompressible flows, SIAM Review, 59 (2017), pp. 492–544, https://doi.org/10.1137/15m1047696.
  * [28] M.-J. Lai and L. L. Schumaker, Spline functions on triangulations, no. 110 in Encyclopedia of Mathematics and its Applications, Cambridge University Press, 2007, https://doi.org/10.1017/CBO9780511721588.
  * [29] Y.-J. Lee, J. Wu, J. Xu, and L. Zikatanov, Robust subspace correction methods for nearly singular systems, Mathematical Models and Methods in Applied Sciences, 17 (2007), pp. 1937–1963, https://doi.org/10.1142/s0218202507002522.
  * [30] A. Linke, A divergence-free velocity reconstruction for incompressible flows, Comptes Rendus Mathematique, 350 (2012), pp. 837–840, https://doi.org/10.1016/j.crma.2012.10.010.
  * [31] A. Linke, L. G. Rebholz, and N. E. Wilson, On the convergence rate of grad-div stabilized Taylor–Hood to Scott–Vogelius solutions for incompressible flow problems, Journal of Mathematical Analysis and Applications, 381 (2011), pp. 612–626, https://doi.org/10.1016/j.jmaa.2011.03.019.
  * [32] S. E. Malick, A connection between grad-div stabilized FE solutions and pointwise divergence-free FE solutions on general meshes, SIAM Undergraduate Research Online, 9 (2016), pp. 367–383, https://doi.org/10.1137/16S015176.
  * [33] D. S. Malkus and E. T. Olsen, Linear crossed triangles for incompressible media, North-Holland Mathematics Studies, 94 (1984), pp. 235–248, https://doi.org/10.1016/S0304-0208(08)72627-6.
  * [34] J. Morgan and L. R. Scott, A nodal basis for ${C}^{1}$ piecewise polynomials of degree $n\geq 5$, Mathematics of Computation, 29 (1975), pp. 736–740, https://doi.org/10.2307/2005284.
  * [35] J. R. Munkres, Elements of Algebraic Topology, CRC Press, 1984.
  * [36] M. Neilan, Discrete and conforming smooth de Rham complexes in three dimensions, Mathematics of Computation, 84 (2015), pp. 2059–2081, https://doi.org/10.1090/S0025-5718-2015-02958-5.
  * [37] M. Neilan, The Stokes complex: A review of exactly divergence-free finite element pairs for incompressible flows, in 75 Years of Mathematics of Computation: Symposium on Celebrating 75 Years of Mathematics of Computation, November 1-3, 2018, the Institute for Computational and Experimental Research in Mathematics (ICERM), Providence, Rhode Island, vol. 754, American Mathematical Soc., 2020, p. 141, https://doi.org/10.1090/conm/754/15142.
  * [38] L. F. Pavarino, Additive Schwarz methods for the $p$-version finite element method, Numerische Mathematik, 66 (1993), pp. 493–515, https://doi.org/10.1007/BF01385709.
  * [39] J. Qin, On the convergence of some low order mixed finite elements for incompressible fluids, PhD thesis, Penn State, 1994.
  * [40] J. Qin and S. Zhang, Stability and approximability of the $\mathcal{P}$1–$\mathcal{P}$0 element for Stokes equations, International Journal for Numerical Methods in Fluids, 54 (2007), pp. 497–515, https://doi.org/10.1002/fld.1407.
  * [41] M. E. Rognes, Automated testing of saddle point stability conditions, in Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book, A. Logg, K. Mardal, and G. Wells, eds., Springer, 2012, pp. 657–671, https://doi.org/10.1007/978-3-642-23099-8_36.
  * [42] A. Sameh and Z. Tong, The trace minimization method for the symmetric generalized eigenvalue problem, Journal of Computational and Applied Mathematics, 123 (2000), pp. 155–175, https://doi.org/10.1016/S0377-0427(00)00391-5.
  * [43] J. Schöberl, Multigrid methods for a parameter dependent problem in primal variables, Numerische Mathematik, 84 (1999), pp. 97–119, https://doi.org/10.1007/s002110050465.
  * [44] J. Schöberl, J. M. Melenk, C. Pechstein, and S. Zaglmayr, Additive Schwarz preconditioning for $p$-version triangular and tetrahedral finite elements, IMA Journal of Numerical Analysis, 28 (2008), pp. 1–24, https://doi.org/10.1093/imanum/drl046.
  * [45] L. R. Scott, Numerical Analysis, Princeton University Press, 2011.
  * [46] L. R. Scott, Introduction to Automated Modeling with FEniCS, Computational Modeling Initiative, 2018.
  * [47] L. R. Scott, ${C}^{1}$ piecewise polynomials satisfying boundary conditions, Research Report UC/CS TR-2019-18, Dept. Comp. Sci., Univ. Chicago, 2019.
  * [48] L. R. Scott and M. Vogelius, Conforming finite element methods for incompressible and nearly incompressible continua, in Large Scale Computations in Fluid Mechanics, B. E. Engquist, et al., eds., vol. 22 (Part 2), Providence: AMS, 1985, pp. 221–244.
  * [49] L. R. Scott and M. Vogelius, Norm estimates for a maximal right inverse of the divergence operator in spaces of piecewise polynomials, $M^{2}AN$ (formerly R.A.I.R.O. Analyse Numérique), 19 (1985), pp. 111–143, http://eudml.org/doc/193439.
  * [50] G. W. Stewart, Matrix Algorithms. Volume II: Eigensystems, SIAM, 2001, https://doi.org/10.1137/1.9780898718058.
  * [51] J. Xu, Iterative methods by space decomposition and subspace correction, SIAM Review, 34 (1992), pp. 581–613, https://doi.org/10.1137/1034116.
  * [52] S. Zhang, A new family of stable mixed finite elements for the 3D Stokes equations, Mathematics of Computation, 74 (2005), pp. 543–554, https://doi.org/10.1090/S0025-5718-04-01711-9.
  * [53] S. Zhang, On the P1 Powell–Sabin divergence-free finite element for the Stokes equations, Journal of Computational Mathematics, 26 (2008), pp. 456–470, https://www.global-sci.org/intro/article_detail/jcm/8636.html.
  * [54] S. Zhang, Divergence-free finite elements on tetrahedral grids for $k\geq 6$, Mathematics of Computation, 80 (2011), pp. 669–695, https://doi.org/10.1090/S0025-5718-2010-02412-3.
  * [55] S. Zhang, Quadratic divergence-free finite elements on Powell–Sabin tetrahedral grids, Calcolo, 48 (2011), pp. 211–244, https://doi.org/10.1007/s10092-010-0035-4.


MIDLMedical Imaging with Deep Learning 2019 MIDL 2019 – Extended Abstract
Track Denis Prokopenko1,2<EMAIL_ADDRESS>
Joël Valentin Stadelmann2<EMAIL_ADDRESS>
Heinrich Schulz3<EMAIL_ADDRESS>
Steffen Renisch3<EMAIL_ADDRESS>
Dmitry V. Dylov1<EMAIL_ADDRESS>
1 Skolkovo Institute of Science and Technology, Moscow, Russian Federation
2 Philips Innovation Labs RUS, Moscow, Russian Federation
3 Philips GmbH Innovative Technologies, Hamburg, Germany

# Synthetic CT Generation from MRI Using Improved DualGAN

###### Abstract

Synthetic CT image generation from MRI scan is necessary to create
radiotherapy plans without the need of co-registered MRI and CT scans. The
chosen baseline adversarial model with cycle consistency permits unpaired
image-to-image translation. Perceptual loss function term and coordinate
convolutional layer were added to improve the quality of translated images.
The proposed architecture was tested on paired MRI-CT dataset, where the
synthetic CTs were compared to corresponding original CT images. The MAE
between the synthetic CT images and the real CT scans is $61$ HU computed
inside of the true CTs body shape.

###### keywords:

Machine Learning, Deep Learning, Radiology, Computed Tomography, Magnetic
Resonance Imaging.

##  1 Introduction

Radiotherapy (RT) requires a personalised preparation of the treatment plan
and, especially, a pre-treatment assessment of the radiation dose. The method
is based on two examinations: magnetic resonance imaging (MRI) and computed
tomography (CT). The MRI helps to locate the tumour and to outline its shape
due to superior soft tissue contrast [Karlsson et al.(2009)Karlsson, Karlsson,
Nyholm, Amies, and Zackrisson]. The CT scan helps to obtain an attenuation map
of a body part from which the radiation plan is derived. Following the MRI and
CT procedures, their relative alignment, and dose calculation, the patient
then undergoes the RT [Coy and Kennelly(1980), Khuntia et al.(2006)Khuntia,
Brown, Li, and Mehta]. However, an additional bias appears due to shifting
errors between body alignments in MRI and CT devices. In the modern approach
of synthetic CT generation, a patient has to participate only in an MRI
procedure, with the advantage of having an error-free registration between the
CT and the MRI volumes [Jonsson et al.(2010)Jonsson, Karlsson, Karlsson, and
Nyholm].

This paper focuses on the unpaired GAN-based approach to the MRI-to-CT image
translation. Our contributions include the upgrade of cycle consistent model
with an additional loss function term: the perceptual loss function term that
enables a comparison of the high-level image representations obtained with the
VGG-16 pre-trained model [Johnson et al.(2016)Johnson, Alahi, and Fei-Fei]. We
also include a coordinate convolutional layer [Liu et al.(2018)Liu, Lehman,
Molino, Such, Frank, Sergeev, and Yosinski], which helps to localise spatially
convolutional properties [Hofmann et al.(2008)Hofmann, Steinke, Scheel,
Charpiat, Farquhar, Aschoff, Brady, Schölkopf, and Pichler]. The proposed
architecture was trained on the brain data with a set of almost overlapping
MRI and CT scans. We present meaningful results of synthetic CT image
generation.

##  2 Methods

The DualGAN architecture [Yi et al.(2017)Yi, Zhang, Tan, and Gong] consists of
two image generators, which form a cycle, and two discriminators, see
img:dualgan. The first generator $G_{MRI\rightarrow CT}$ is trained to
translate MRI images to CT images; the second generator $G_{CT\rightarrow
MRI}$ translates images from CT to MRI domain. The cycle allows comparing the
reconstructed image and original to evaluate the quality of generators without
the need in paired data.

Figure 1: DualGAN architecture.

The perceptual loss function (VGG) builds on the idea of a feature matching,
where high-level representations of two images are compared by mean squared
error. Although extracted features by pre-trained VGG-16 are not optimised for
tomographic images, their use is still relevant because the features extracted
in an identical way for both compared images.

The coordinate convolutional layer (CC) allows taking into account the spatial
information of the images by concatenating two additional $x$ and $y$
coordinates slices with the tensor representation of the image. The coordinate
convolutional layer in application to the MRI-to-CT translation helps to
distinguish black pixels of MRI image, which could represent either a bone or
air.

##  3 Experiments

Three medical datasets were obtained to train and test the considered methods
of unpaired MRI to CT translation. Each set was divided into a train and test
parts in a 7:3 ratio. The volumes were normalised and preprocessed.

The CPTAC Phase 3 dataset includes the MRI T1-weighted images of $7$ patients.
Each 3D volume of a patient contains $22-24$ slices in the axial anatomical
plane comprised between the lower nose and the top of the head region
[CPTAC(2018), Clark et al.(2013)Clark, Vendt, Smith, Freymann, Kirby, Koppel,
Moore, Phillips, Maffitt, Pringle, et al.].

The head-and-neck cancer dataset consists of CT scans of $61$ patients. 3D
volumes include 61 - 94 slices from the shoulders, neck and lower part of a
head up to the eyes [Vallières et al.(2017)Vallières, Kay-Rivest, Perrin,
Liem, Furstoss, Khaouam, Nguyen-Tan, Wang, and Sultanem, Vallières et
al.(2017)Vallières, Kay-Rivest, Perrin, Liem, Furstoss, Aerts, Khaouam,
Nguyen-Tan, Wang, Sultanem, et al.].

The private dataset contained CT images, masks and paired MRI T1-weighted
images of $10$ patients. The volumes consist of $66-137$ slices in the axial
plane, which are comprised between the teeth and the top of the head region.
Train part was used in an unpaired way.

All models were trained separately on CPTAC with Head-and-neck cancer datasets
together and the private dataset. The qualitative results and quantitative
evaluation of the performance were obtained on the test part of the private
dataset.

In this work, four different architectures were considered and compared:
DualGAN, DualGAN with the coordinate convolutional layer (DualGAN, CC),
DualGAN with the perceptual loss function term (DualGAN, VGG) and DualGAN with
perceptual loss function term and coordinate convolutional layer (DualGAN,
VGG, CC). All of them work with three-channel images allowing the application
of VGG network.

Performance evaluation was calculated using one channel image representations.
The quantitative results of different translation configurations can be seen
in tab:compare, while qualitative results are presented in img:test, which
shows by column the real MRI slices, the synthetic CT images, the real paired
CT slices and the difference between synthetic and original CTs.

Table 1: Performance comparison of different MRI to CT translation
configurations.

Configuration | MAE, HU $\downarrow$ | PSNR, dB $\uparrow$ | SSIM $\uparrow$
---|---|---|---
DualGAN | $62.95\pm 1.17$ | $16.82\pm 0.83$ | $0.79\pm 0.02$
DualGAN, CC | $63.27\pm 2.00$ | $16.96\pm 0.71$ | $0.78\pm 0.03$
DualGAN, VGG | $66.52\pm 3.74$ | $16.74\pm 1.08$ | $0.78\pm 0.03$
DualGAN, VGG, CC | $60.83\pm 2.20$ | $17.21\pm 1.00$ | $0.80\pm 0.03$

Figure 2: Tanslation made by DualGAN with perceptual loss function term and
coordinate convolutional layer.

##  4 Conclusions

The presented model achieves MAE of $61$ HU and SSIM of $0.8$, which compares
favourably to scientific literature [Wolterink et al.(2017)Wolterink, Dinkla,
Savenije, Seevinck, van den Berg, and Isgum]. The translation architecture
transforms the initial image retaining the structural information. While the
DualGAN with the perceptual loss function term and the DualGAN with the
coordinate convolutional layer are not themselves superior to the DualGAN
model, their triad combination showed meaningful improvement. Besides, our
model works with the unpaired images, and the visual examination confirmed
that the use of the perceptual loss term and the coordinate convolutional
layer enhances the appearance of the resulting images. The error of
translation increases with the complexity of the inner structures, for
instance, in the nose and teeth regions.

##  5 Acknowledgement

Data used in this publication were generated by the National Cancer Institute
Clinical Proteomic Tumor Analysis Consortium (CPTAC).

## References

  * [Clark et al.(2013)Clark, Vendt, Smith, Freymann, Kirby, Koppel, Moore, Phillips, Maffitt, Pringle, et al.] Kenneth Clark, Bruce Vendt, Kirk Smith, John Freymann, Justin Kirby, Paul Koppel, Stephen Moore, Stanley Phillips, David Maffitt, Michael Pringle, et al.  The cancer imaging archive (tcia): maintaining and operating a public information repository.  _Journal of digital imaging_ , 26(6):1045–1057, 2013.
  * [Coy and Kennelly(1980)] Peter Coy and GM Kennelly.  The role of curative radiotherapy in the treatment of lung cancer.  _Cancer_ , 45(4):698–702, 1980.
  * [CPTAC(2018)] National Cancer Institute Clinical Proteomic Tumor Analysis Consortium CPTAC.  Radiology Data from the Clinical Proteomic Tumor Analysis Consortium Glioblastoma Multiforme [CPTAC-GBM] collection [Data set]. The Cancer Imaging Archive. https://doi.org/10.7937/k9/tcia.2018.3rje41q1, 2018.
  * [Hofmann et al.(2008)Hofmann, Steinke, Scheel, Charpiat, Farquhar, Aschoff, Brady, Schölkopf, and Pichler] Matthias Hofmann, Florian Steinke, Verena Scheel, Guillaume Charpiat, Jason Farquhar, Philip Aschoff, Michael Brady, Bernhard Schölkopf, and Bernd J Pichler.  Mri-based attenuation correction for pet/mri: a novel approach combining pattern recognition and atlas registration.  _Journal of nuclear medicine_ , 49(11):1875–1883, 2008.
  * [Johnson et al.(2016)Johnson, Alahi, and Fei-Fei] Justin Johnson, Alexandre Alahi, and Li Fei-Fei.  Perceptual losses for real-time style transfer and super-resolution.  In _European Conference on Computer Vision_ , pages 694–711. Springer, 2016.
  * [Jonsson et al.(2010)Jonsson, Karlsson, Karlsson, and Nyholm] Joakim H Jonsson, Magnus G Karlsson, Mikael Karlsson, and Tufve Nyholm.  Treatment planning using mri data: an analysis of the dose calculation accuracy for different treatment regions.  _Radiation oncology_ , 5(1):62, 2010.
  * [Karlsson et al.(2009)Karlsson, Karlsson, Nyholm, Amies, and Zackrisson] Mikael Karlsson, Magnus G Karlsson, Tufve Nyholm, Christopher Amies, and Björn Zackrisson.  Dedicated magnetic resonance imaging in the radiotherapy clinic.  _International Journal of Radiation Oncology* Biology* Physics_ , 74(2):644–651, 2009.
  * [Khuntia et al.(2006)Khuntia, Brown, Li, and Mehta] Deepak Khuntia, Paul Brown, Jing Li, and Minesh P Mehta.  Whole-brain radiotherapy in the management of brain metastasis.  _Journal of Clinical Oncology_ , 24(8):1295–1304, 2006.
  * [Liu et al.(2018)Liu, Lehman, Molino, Such, Frank, Sergeev, and Yosinski] Rosanne Liu, Joel Lehman, Piero Molino, Felipe Petroski Such, Eric Frank, Alex Sergeev, and Jason Yosinski.  An intriguing failing of convolutional neural networks and the coordconv solution.  _arXiv preprint arXiv:1807.03247_ , 2018.
  * [Vallières et al.(2017)Vallières, Kay-Rivest, Perrin, Liem, Furstoss, Aerts, Khaouam, Nguyen-Tan, Wang, Sultanem, et al.] Martin Vallières, Emily Kay-Rivest, Léo Jean Perrin, Xavier Liem, Christophe Furstoss, Hugo JWL Aerts, Nader Khaouam, Phuc Felix Nguyen-Tan, Chang-Shu Wang, Khalil Sultanem, et al.  Radiomics strategies for risk assessment of tumour failure in head-and-neck cancer.  _Scientific reports_ , 7(1):10117, 2017.
  * [Vallières et al.(2017)Vallières, Kay-Rivest, Perrin, Liem, Furstoss, Khaouam, Nguyen-Tan, Wang, and Sultanem] Martin Vallières, Emily Kay-Rivest, Léo Jean Perrin, Xavier Liem, Christophe Furstoss, Nader Khaouam, Phuc Félix Nguyen-Tan, Chang-Shu Wang, and Khalil Sultanem.  Data from Head-Neck-PET-CT. The Cancer Imaging Archive. https://doi.org/10.7937/K9/TCIA.2017.8oje5q00, 2017.
  * [Wolterink et al.(2017)Wolterink, Dinkla, Savenije, Seevinck, van den Berg, and Isgum] Jelmer M. Wolterink, Anna M. Dinkla, Mark H. F. Savenije, Peter R. Seevinck, Cornelis A. T. van den Berg, and Ivana Isgum.  Deep MR to CT synthesis using unpaired data.  _CoRR_ , abs/1708.01155, 2017.  URL http://arxiv.org/abs/1708.01155.
  * [Yi et al.(2017)Yi, Zhang, Tan, and Gong] Zili Yi, Hao Zhang, Ping Tan, and Minglun Gong.  Dualgan: Unsupervised dual learning for image-to-image translation.  _CoRR_ , abs/1704.02510, 2017.  URL http://arxiv.org/abs/1704.02510.


# FABLE : Fabric Anomaly Detection Automation Process

Simon Thomine  University of technology of Troyes
AQUILAE
Troyes, France
<EMAIL_ADDRESS>Hichem Snoussi  University of technology of Troyes
Troyes, France
<EMAIL_ADDRESS>Mahmoud Soua  AQUILAE
Troyes, France
<EMAIL_ADDRESS>
###### Abstract

Unsupervised anomaly in industry has been a concerning topic and a stepping
stone for high performance industrial automation process. The vast majority of
industry-oriented methods focus on learning from good samples to detect
anomaly notwithstanding some specific industrial scenario requiring even less
specific training and therefore a generalization for anomaly detection. The
obvious use case is the fabric anomaly detection, where we have to deal with a
really wide range of colors and types of textile and a stoppage of the
production line for training could not be considered. In this paper, we
propose an automation process for industrial fabric texture defect detection
with a specificity-learning process during the domain-generalized anomaly
detection. Combining the ability to generalize and the learning process offer
a fast and precise anomaly detection and segmentation. The main contributions
of this paper are the following: A domain-generalization texture anomaly
detection method achieving the state-of-the-art performances, a fast specific
training on good samples extracted by the proposed method, a self-evaluation
method based on custom defect creation and an automatic detection of already
seen fabric to prevent re-training.

###### Index Terms:

Domain-Generalization, unsupervised, anomaly, unseen, knowledge distillation,
student-teacher, memory banks, fabric, automation.

##  I Introduction

Unsupervised anomaly detection in industry is a vast topic, since there are a
lot of possible applications. In this paper, we focus on fabric anomaly, which
is a concerning topic for industry. The specificity of fabric is the pattern
in the structure and if we manage to understand that pattern we can extract
anomalies. Several methods have been introduced for industry anomaly detection
using MVTEC AD [1] the dataset that gathers textures (carpet, leather, grid,
wood, and tile) and objects (bottle, cable, capsule, hazelnut, metal nut,
pill, screw, toothbrush, transistor and zipper). These methods could achieve
high performance. However, they rely on object/texture specific unsupervised
learning without generalization capacity. Recently, knowledge-distillation
based methods have been introduced for the unsupervised anomaly detection task
[2]. It consists of a student-teacher model focusing on the bottom layers of
the network as they represent the edges, color and shapes information. We used
the same approach to design a domain-generalized texture anomaly detection
method with the ability to detect defects on unseen textures and to select
good samples for a texture-specific unsupervised anomaly detection model. In
fabric industry, many types and colors of fabric are analyzed, and it would be
impossible to rely on a specific training on good samples for each type of
fabric without slowing the industrial process.
Therefore, we propose a complete data processing chain for a robust, fast and
adaptive texture specific anomaly detection and localization. Our method is
based on four main modules: a domain-generalized texture anomaly detector, a
fast texture specific training/inference, an auto-evaluation process of our
specific model and an automatic already-seen fabric detection to avoid
retraining an existing model.
The paper is organized as follows. In section II, we review the related work
especially on MVTEC dataset and present the different approaches proposed in
literature for domain-generalized and classic unsupervised anomaly detection.
In section III, we present an enhanced domain-generalized texture defect
detection method. In section IV, we present the specific learning method, the
auto-evaluation process and the already seen texture recognition. Section V is
dedicated to the analysis of the results. Section VI concludes the paper.

##  II Related works

As our proposed methods address two specific tasks, we first present the state
of the art on domain-generalized texture anomaly detection and then the state
of the art on unsupervised defect detection of known objects.

###  II-A Domain-generalized texture anomaly detection

Domain-generalized anomaly detection is an important topic for optimal
industrial process, since in specific industrial fields, the type of textures
often changes. The most obvious example is certainly fabric anomaly detection
where fabric can have different colors (red, blue, striped) and types (cotton,
polyester, silk, etc). The main objective is to detect defects on any type of
fabric without resorting to a time-consuming training. The feature extraction
from a pretrained classifier offers the most promising results with different
types of networks such as an episodic training [3], the use of extrinsic and
intrinsic aspects [4] and multiscale feature extractor with co-attention
modules [5].

###  II-B Unsupervised anomaly detection on known objects

More commonly, unsupervised anomaly detection deals with the problem of
detecting defects on an object or texture based on only good samples. In
industry or security scenarios, we often have a low rate of defects with a
vast number of different defect types which would lead to a time-consuming
annotation and a possibly non-pertinent classification if all the anomaly
types are not considered[6]. To tackle this question, several methods emerged
proposing different types of algorithms such as autoencoders [7] and
variational autoencoder variants [8] [9]. Another common way of detecting
anomalies is Generative Adversarial Networks (GAN) introduced by [10] adapted
to unsupervised anomaly detection such as Ano-GAN [11], G2D [12] and OCR-GAN
[13]. More recently, approaches using a pretrained classifier has been at the
heart of the research in industrial anomaly detection and offers outstanding
performance. There are three main feature extraction-based approaches:
normalizing flow, knowledge distillation and memory banks. The normalizing
flow approach consists of a flow training based on relevant features of good
samples from a pretrained network such as AlexNet [14], Resnet [15] or
efficient-net [16] trained on imageNet. Different strategies were used to
enhance performance, such as a 2D flow [17] or a cross-scale flow [18].
Another interesting approach is the use of a memory bank to extract relevant
information from different good samples and to use this memory bank to compare
and detect if there is an anomaly [19]. Finally, the concept of knowledge
distillation was adapted for unsupervised anomaly detection and localization
[2]. The idea is to train a student network based on the output features of a
teacher (already pretrained for a classification purpose) and on good samples.
The student will be able to reproduce teacher features on a good sample, but
will not be as precise for a defective sample. Several methods used this
principle with different strategies such as a multi-layer feature selection
[2], an asymmetric student teacher [20], a coupled-hypersphere-based feature
adaptation [21] and a mixed-teacher approach [22].

##  III Knowledge distillation generalization

The proposed model is based on the knowledge distillation framework, where a
pretrained network is used as a teacher and a student network is trained to
reproduce the teacher output on good samples. The student network is then
expected to not be able to reproduce teacher features on defective samples, a
property which is used to detect abnormal samples.
For domain generalization, we propose to train the student on different types
of textures and using many teachers to guarantee generalization. In order to
achieve this objective, we first constitute a new dataset based on fabric
datasets [23] which regroups different categories of textures with different
quality and homogeneity.

Figure 1:  Samples employed for the custom fabric dataset (extracted from the
fabrics dataset [23])

Then, to tackle the problem of texture domain generalization, we used a
specific student teacher architecture with different branches based on the
paradigm that each pre-trained classifier have a different bias towards
classification.
In terms of layer selection, the deeper a layer, the more the information
relates to the context and conversely, the shallower a layer, the more
information it contains on contours, edges, and colors. Based on different
layer configurations, we show that for the purpose of texture domain
generalization, mid-level features would be the best choice to combine texture
specific information such as contours and edges and a general vision of what a
texture is.
At least two classifiers are needed to attenuate each bias. We have used
Resnet18 and EfficientNet-b0 for computation time speed and meaningful
features.
To fully exploit each classifier information, we used a parallel architecture
which can be seen as a multiple teachers/multiple students architecture where
the training happen independently for each classifier, only the anomaly score
is calculated with the two networks outputs. Our framework is an adaptation of
MixedTeacher [22] with a different layer selection strategy. The first Resnet
layer is not used as its output features are too specific to training dataset
textures. We used the features of the three first residual blocks of Resnet18
and the last 2 convolutional blocks of efficientNet-b0. As in [22], we used a
reduced version of the Resnet18 model with a reduction of the block size and a
reduction of the dimension of each layer with an adaptive average pooling,
while we keep the same architecture for the EfficientNet part.

Given a training dataset of images without anomaly
${D=[I_{1},I_{2},...,I_{n}]}$, our goal is to extract the information of $L$
mid-level layers. For an image ${I_{k}}\in R^{w*h*c}$ where $w$ is the width,
$h$ the height, and $c$ the number of channel, the teacher outputs features
$F_{t}^{l}(I_{k})\in R^{w_{l}*h_{l}*c_{l}}$ and $F_{s}^{l}(I_{k})\in
R^{w_{l}/2*h_{l}/2*c_{l}/2}$ with $l>1$ and $F_{s}^{l}(I_{k})\in
R^{w_{l}*h_{l}*c_{l}}$ if $l=1$. The loss is obtained by applying the $l2$
distance of normalized feature vectors for each pixel of the feature map and
summing them. For the Resnet student part, we used an adaptive average pooling
layer on teacher features. The used layers are $l=\left\\{1,2,3\right\\}$ for
the Resnet part and $l=\left\\{5,6\right\\}$ for the EfficientNet part.

Pixel loss for the resnet part is defined in the following Eq.1:

$loss^{l}(I_{k})_{ij}=\frac{1}{2}\lVert
norm(AAP(F_{Resnet18}^{l}(I_{k}))_{ij})-norm(F_{s}^{l}(I_{k})_{ij})\rVert$ (1)

where AAP refers to Adaptive Average Pooling. For the EfficientNet part, pixel
loss is defined in the following Eq.2:

$loss^{l}(I_{k})_{ij}=\frac{1}{2}\lVert
norm(F_{EffNetb0}^{l}(I_{k})_{ij})-norm(F_{s}^{l}(I_{k})_{ij})\rVert$ (2)

For the layer $l$, the loss is defined as:

$loss^{l}(I_{k})=\frac{1}{w_{l}h_{l}}\sum_{i=1}^{w_{l}}\sum_{j=1}^{h_{l}}loss^{l}(I_{k})_{ij}$
(3)

and finally, for the total loss is written as:

$loss(I_{k})=\sum^{l}loss^{l}(I_{k})$ (4)

Figure 2:  Architecture of Resnet student teacher (left) and EfficientNet
student teacher (right)

##  IV Auto-learning process for industrial deployment

The previous part was presented in the context of industrial efficiency, where
it was not allowed to retrain for every new type/color of texture/fabric. The
objective of this section is to propose a general classifier for handling the
anomaly detection role while we gather enough images and train a specific
model for increased efficiency.
This section is divided in 3 parts: training and self-evaluation, recognition
of an already trained type of fabric, and a typical industrial use-case in
fabric industry.

###  IV-A Training and self-evaluation

Given the deployment constraints, we considered different criteria for the
choice of the student-teacher network architecture: (i) the inference and
training time, (ii) the performance and (iii) the robustness to defective
samples in the training set. We also considered the possibility of running the
process on several asynchronous defect detectors. The model Reduced Student
proposed in [22] is a good candidate. Thanks to its reduced architecture, we
can train a specific model in an acceptable time. To minimize the number of
potential defective samples in the training, we gathered the samples with
acceptable anomaly score from the domain-generalized model, i.e samples
classified as good samples. Based on a test-error approach, we determined the
optimal number of epochs (during specific training) where the specific model
becomes better than the domain-generalized one so that we can start using the
best model even if the complete training is not finished.
The self-evaluation part is based on two types of data: (i) the first type is
defective samples detected by the domain-generalized anomaly detector and (ii)
the second type is generated data with a procedure inspired by DRAEM [9]:
Perin noise and the texture database dtd [24]. We used the same approach to
generate non-absurd defects to self-evaluate our model.

Figure 3:  Custom defective samples generated with Perin noise

###  IV-B Already seen fabric recognition

To guarantee an automated anomaly detector without the help of an operator for
selecting an already-trained model, we propose an algorithm to precisely
recognize a fabric type already considered previously. For each trained model,
we save $x$ extracted features from the specific model on good samples reduced
using coreset subsampling introduced in PatchCore [19] to guarantee fast
computation. Each specific model is also saved in a model bank of $N$ models
and linked to its features in a feature bank. When we have to decide if the
fabric was already seen, we calculate the sample/model proximity by extracting
features from all trained specific models from the model bank, applying the
coreset subsampling and comparing these features to the $x$ features from the
feature bank of each specific model with cosine similarity distance as
described in equation 6. We then compute the intra-class proximity by
calculating the cosine similarity between the $x$ features of the same model
as reported in equation 7. The proximity score is defined as the absolute
value of the difference between the sample/model proximity and the intra-class
proximity. We finally make the decision by comparing the maximum proximity
score with a $similarityThreshold$. The threshold is chosen based on what is
known about the similarity between the fabric.

Figure 4:  Fabric Recognition process : The extraction of anomaly-free images
features in the testing phase stops when the $x$ number of elements is
reached, the model bank is the bank of all previously trained specific models

Even though it may seem laborious if the model bank is consistent, it is still
real-time deployable thanks to the inference speed of the reduced student
architecture proposed in the previous subsection and the coreset subsampling,
as we show in the experiment part. This is by far the most accurate method for
comparing a new piece of fabric with a previously seen one and, we believe, it
is still usable even in a specific case of thousands of specific models.
The cosine similarity formula is:

$sim(feat_{A},feat_{B})=\frac{feat_{A}.feat_{B}}{\left\|feat_{A}\right\|\left\|feat_{B}\right\|}$
(5)

with $feat_{A}$ and $feat_{B}$ the extracted features. The sample/model
proximity is defined as:

$prox_{sm}(S,Model)=\frac{1}{x}\sum_{i=1}^{x}sim(feat_{S},feat_{i})$ (6)

The intra-class proximity is defined as:

$prox_{ic}(Model)=\frac{1}{x(x-1)}\sum_{i=1}^{x}\sum_{j=1,i\neq
j}^{x}sim(feat_{i},feat_{j})$ (7)

The proximity score is :

$proxScore(S,Model)=abs(prox_{sm}(S,Model)-prox_{ic}(Model))$ (8)

And the already-seen decision is described as :

$\underset{i\in N}{\text{max}}(proxScore(S,Model_{i}))>similarityThreshold$
(9)

###  IV-C Typical use-case : fabric industry

To demonstrate the effectiveness of our method, we describe a typical real
defect analysis use case. In a vast majority of mid-range clothing industry,
the fabric is analyzed several times during the whole fabrication process by
operators that scroll the fabric and look for defects. This is a laborious job
and often distraction occurs resulting in a globally low detection percentage
of defects, not to mention the difficulty for the eyes to look at certain
fabric categories such as striped fabrics. Our automated process aims at
assisting the operator for the classification task and to speed up the
scrolling of the fabric. The operator is still needed since he has to install
the fabric roll on the machine and to verify the defect classification done by
the domain generalized model since the accuracy is still low for a full
automation process.
For every fabric roll, the process start with an identification of the fabric
to control. Two different cases may happen:
\- If this type of fabric has never been seen, a specific training is started
while still doing the anomaly detection with the domain-generalized model, we
may have to slow the scrolling of the fabric during the training depending on
the computational power. When the trained model becomes better than the
domain-generalized model, we used it instead, even if the training is not
completely finished. When the training is finished, we used the completely
trained model for anomaly detection while keeping some features of defective
samples for the recognition part.
\- If this type of fabric is already-seen, the specific trained model is used
for anomaly detection.
The process is fully automated and does not require any help from the operator
except for the activation or deactivation, which could be done also by
connecting the visiting machine with the central unit to send an activation
signal.

##  V Experiments

This section is divided into 3 parts: the analysis of the domain-
generalization model compared with state of the art for different training
configurations, the analysis of the training speed and inference speed of our
model and finally the estimation of the number of required epochs on a
specific training to outperform the domain-generalization algorithm.

###  V-A State-of-the-art comparison

For the evaluation of our model, we used two different databases for training.
For the “MVTEC” one, we trained the DG model on all good samples of MVTEC AD
textures except the one we are testing on to reproduce the evaluation protocol
of the other state-of-the-art papers. The “cotton” one is trained on the
custom fabric dataset presented in section III and was created for fabric
anomaly which explain the SOTA performances on carpet and leather. The results
are presented in table I.
For the training, we used stochastic gradient descent with a learning rate of
0.4 for 200 epochs with a batch size of 16. Both networks are pretrained on
ImageNet. We resized all the images to 256x256, keeping 80% for training and
20% for validation. We kept the checkpoint with the lowest validation loss.

TABLE I: AUC comparison between our method and existing ones on MVTEC AD

textures | Epi-FRC+[3] | EISNet+[4] | DGTSAD[5] | Ours(MVTEC) | Ours(Coton)
---|---|---|---|---|---
carpet | 0.916 | 0.982 | 0.943 | 0.985 | 0.996
leather | 1.000 | 1.000 | 1.000 | 0.991 | 0.996
wood | 0.941 | 0.979 | 0.962 | 0.999 | 0.948
tile | 0.951 | 0.851 | 0.994 | 0.965 | 0.964
grid | 0.725 | 0.728 | 0.730 | 0.937 | 0.944
mean | 0.907 | 0.909 | 0.918 | 0.975 | 0.969

Figure 5:  Outputs of detected anomaly with our domain-generalization model

As seen in table I, the 2 types of training offers approximately the same mean
AUC but the dataset “cotton” only contains one hundred images and is supposed
to be effective for fabric defects detection whereas it shows the best results
on both carpet, leather and grid which contains patterns [25] and are the most
fabric-like textures of the dataset MVTEC AD.
Comparing to the state-of-the-art methods, our approach obtains 0.057 AUC more
than the previous best model, which is an excellent improvement and is in the
way for closing the gap between domain generalization and classical anomaly
detection for textures.

###  V-B Inference speed

For the inference speed, all the following tests were done with a RTX 2080Ti.
The training parameters are the same as the previous part except for the
number of epochs where we limited it to 100. To perform these experiments, we
used an optimal algorithm for fast processing in batch of 8 which outperforms
the classic anomalib [26] in terms of inference speed on knowledge
distillation methods for anomaly detection (see table II). For the training
time, with 100 epochs and 100 images, we report 4 minutes and 20 seconds.

TABLE II: Inference speed

Category | Ours(Domain-Generalized) | Ours(Specific)
---|---|---
FPS | 333 | 1111
Latency (ms) | 3 | 0.9

###  V-C Necessary training epochs before model replacement

The main purpose of the whole method is to detect defects as precisely as
possible during the whole process. We need to estimate at which point the
specific model surpasses the DG model in terms of AUC. In order to find this
number, we tested our specific model at different epochs, and we compared
results to the DG model performance.

TABLE III: Epoch performance

textures | 10 epochs | 20 epochs | 30 epochs | 100 epochs | DG model
---|---|---|---|---|---
carpet | 0.915 | 0.920 | 0.963 | 1.0 | 0.996
leather | 0.949 | 0.974 | 0.990 | 0.997 | 0.996
wood | 0.988 | 0.989 | 0.995 | 0.996 | 0.948
tile | 0.986 | 0.987 | 0.991 | 0.987 | 0.964
mean | 0.959 | 0.967 | 0.984 | 0.995 | 0.969

Based on the mean results, we could consider using the specific model after 30
epochs of training i.e one minute and 30 seconds, but considering the scores
of carpet and leather which are the most fabric-like textures, it may be a
counter performance to switch to the specific model before the end of the
training if the DG model is close enough to real distribution.

##  VI Conclusion and future work

We proposed an industry-ready automation process for industrial defect
detection, especially deployable for fabric with fast inference results for
both domain-generalized and specific model. The outstanding capability of our
architecture to mutual aid between humans and artificial intelligence makes it
a great and reliable tool for visual inspection. Nevertheless, several
improvements could be considered. According to [6], semi-supervised can easily
surpass unsupervised anomaly detection even with a few annotated anomalies and
this could be applied to our method by using the DG detected anomalies to
train a semi supervised specific model to increase the detection performances.
Some automation improvement could be considered in terms of potential use of
the method, such as a direct software included in the fabric visiting machine
to monitor speed, stop and dysfunction mode more precisely.

## References

  * [1] Paul Bergmann, Michael Fauser, David Sattlegger and Carsten Steger “MVTec AD — A Comprehensive Real-World Dataset for Unsupervised Anomaly Detection” In _2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ Long Beach, CA, USA: IEEE, 2019, pp. 9584–9592 DOI: 10.1109/CVPR.2019.00982
  * [2] Guodong Wang, Shumin Han, Errui Ding and Di Huang “Student-Teacher Feature Pyramid Matching for Anomaly Detection” In _arXiv:2103.04257 [cs]_ , 2021 arXiv: http://arxiv.org/abs/2103.04257
  * [3] Da Li, Jianshu Zhang, Yongxin Yang, Cong Liu, Yi-Zhe Song and Timothy Hospedales “Episodic Training for Domain Generalization” In _2019 IEEE/CVF International Conference on Computer Vision (ICCV)_ Seoul, Korea (South): IEEE, 2019, pp. 1446–1455 DOI: 10.1109/ICCV.2019.00153
  * [4] Shujun Wang, Lequan Yu, Caizi Li, Chi-Wing Fu and Pheng-Ann Heng “Learning from Extrinsic and Intrinsic Supervisions for Domain Generalization” arXiv, 2020 arXiv: http://arxiv.org/abs/2007.09316
  * [5] Shang-Fu Chen, Yu-Min Liu, Chia-Ching Lin, Trista Pei-Chun Chen and Yu-Chiang Frank Wang “Domain-Generalized Textured Surface Anomaly Detection” arXiv, 2022 arXiv: http://arxiv.org/abs/2203.12304
  * [6] Songqiao Han, Xiyang Hu, Hailiang Huang, Mingqi Jiang and Yue Zhao “ADBench: Anomaly Detection Benchmark” arXiv, 2022 arXiv: http://arxiv.org/abs/2206.09426
  * [7] Shuang Mei, Yudan Wang and Guojun Wen “Automatic Fabric Defect Detection with a Multi-Scale Convolutional Denoising Autoencoder Network Model” In _Sensors_ 18.4, 2018, pp. 1064 DOI: 10.3390/s18041064
  * [8] Quoc Phong Nguyen, Kar Wai Lim, Dinil Mon Divakaran, Kian Hsiang Low and Mun Choon Chan “GEE: A Gradient-based Explainable Variational Autoencoder for Network Anomaly Detection” In _arXiv:1903.06661 [cs, stat]_ , 2019 arXiv: http://arxiv.org/abs/1903.06661
  * [9] Vitjan Zavrtanik, Matej Kristan and Danijel Skočaj “DRAEM – A discriminatively trained reconstruction embedding for surface anomaly detection” arXiv, 2021 arXiv: http://arxiv.org/abs/2108.07610
  * [10] Ian J. Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville and Yoshua Bengio “Generative Adversarial Networks” In _arXiv:1406.2661 [cs, stat]_ , 2014 arXiv: http://arxiv.org/abs/1406.2661
  * [11] Thomas Schlegl, Philipp Seeböck, Sebastian M. Waldstein, Georg Langs and Ursula Schmidt-Erfurth “f-AnoGAN: Fast unsupervised anomaly detection with generative adversarial networks” In _Medical Image Analysis_ 54, 2019, pp. 30–44 DOI: 10.1016/j.media.2019.01.010
  * [12] Masoud Pourreza, Bahram Mohammadi, Mostafa Khaki, Samir Bouindour, Hichem Snoussi and Mohammad Sabokrou “G2D: Generate to Detect Anomaly” event-place: Waikoloa, HI, USA In _2021 IEEE Winter Conference on Applications of Computer Vision (WACV)_ IEEE, 2021, pp. 2002–2011 DOI: 10.1109/WACV48630.2021.00205
  * [13] Yufei Liang, Jiangning Zhang, Shiwei Zhao, Runze Wu, Yong Liu and Shuwen Pan “Omni-frequency Channel-selection Representations for Unsupervised Anomaly Detection” In _arXiv:2203.00259 [cs]_ , 2022 arXiv: http://arxiv.org/abs/2203.00259
  * [14] Kaiming He, Xiangyu Zhang, Shaoqing Ren and Jian Sun “Deep Residual Learning for Image Recognition” arXiv, 2015 arXiv: http://arxiv.org/abs/1512.03385
  * [15] Alex Krizhevsky, Ilya Sutskever and Geoffrey E. Hinton “ImageNet classification with deep convolutional neural networks” In _Communications of the ACM_ 60.6, 2017, pp. 84–90 DOI: 10.1145/3065386
  * [16] Mingxing Tan and Quoc V. Le “EfficientNet: Rethinking Model Scaling for Convolutional Neural Networks” arXiv, 2020 arXiv: http://arxiv.org/abs/1905.11946
  * [17] Jiawei Yu, Ye Zheng, Xiang Wang, Wei Li, Yushuang Wu, Rui Zhao and Liwei Wu “FastFlow: Unsupervised Anomaly Detection and Localization via 2D Normalizing Flows” In _arXiv:2111.07677 [cs]_ , 2021 arXiv: http://arxiv.org/abs/2111.07677
  * [18] Marco Rudolph, Tom Wehrbein, Bodo Rosenhahn and Bastian Wandt “Fully Convolutional Cross-Scale-Flows for Image-based Defect Detection” In _arXiv:2110.02855 [cs]_ , 2021 arXiv: http://arxiv.org/abs/2110.02855
  * [19] Karsten Roth, Latha Pemula, Joaquin Zepeda, Bernhard Schölkopf, Thomas Brox and Peter Gehler “Towards Total Recall in Industrial Anomaly Detection” In _arXiv:2106.08265 [cs]_ , 2021 arXiv: http://arxiv.org/abs/2106.08265
  * [20] Marco Rudolph, Tom Wehrbein, Bodo Rosenhahn and Bastian Wandt “Asymmetric Student-Teacher Networks for Industrial Anomaly Detection” arXiv, 2022 arXiv: http://arxiv.org/abs/2210.07829
  * [21] Sungwook Lee, Seunghyun Lee and Byung Cheol Song “CFA: Coupled-hypersphere-based Feature Adaptation for Target-Oriented Anomaly Localization” arXiv, 2022 arXiv: http://arxiv.org/abs/2206.04325
  * [22] Simon Thomine, Hichem Snoussi and Mahmoud Soua “MixedTeacher : Knowledge Distillation for fast inference textural anomaly detection”, 2023
  * [23] Christos Kampouris, Stefanos Zafeiriou, Abhijeet Ghosh and Sotiris Malassiotis “Fine-Grained Material Classification Using Micro-geometry and Reflectance” Series Title: Lecture Notes in Computer Science In _Computer Vision – ECCV 2016_ 9909 Cham: Springer International Publishing, 2016, pp. 778–792 DOI: 10.1007/978-3-319-46454-1˙47
  * [24] Mircea Cimpoi, Subhransu Maji, Iasonas Kokkinos, Sammy Mohamed and Andrea Vedaldi “Describing Textures in the Wild” In _2014 IEEE Conference on Computer Vision and Pattern Recognition_ Columbus, OH, USA: IEEE, 2014, pp. 3606–3613 DOI: 10.1109/CVPR.2014.461
  * [25] Henry Y. T. Ngan, Grantham K. H. Pang and Nelson H. C. Yung “Automated fabric defect detection—A review” In _Image and Vision Computing_ 29.7, 2011, pp. 442–458 DOI: 10.1016/j.imavis.2011.02.002
  * [26] Samet Akcay, Dick Ameln, Ashwin Vaidya, Barath Lakshmanan, Nilesh Ahuja and Utku Genc “Anomalib: A Deep Learning Library for Anomaly Detection” arXiv, 2022 arXiv: http://arxiv.org/abs/2202.08341