∎

11institutetext: Nguyen Hieu Thao22institutetext: Delft Center for Systems and
Control, Delft University of Technology, 2628CD Delft, The Netherlands.
Department of Mathematics, School of Education, Can Tho University, Can Tho,
Vietnam. 22email<EMAIL_ADDRESS><EMAIL_ADDRESS>
Oleg Soloviev33institutetext: Delft Center for Systems and Control, Delft
University of Technology, 2628CD Delft, The Netherlands. Flexible Optical
B.V., Polakweg 10-11, 2288 GG Rijswijk, The Netherlands. 33email:
<EMAIL_ADDRESS>
Jacques Noom44institutetext: Delft Center for Systems and Control, Delft
University of Technology, 2628CD Delft, The Netherlands. 44email:
<EMAIL_ADDRESS>
Michel Verhaegen55institutetext: Delft Center for Systems and Control, Delft
University of Technology, 2628CD Delft, The Netherlands. 55email:
<EMAIL_ADDRESS>
# Subpixel image reconstruction using nonuniform defocused images ††thanks:
This project has received funding from the ECSEL Joint Undertaking (JU) under
grant agreement No. 826589. The JU receives support from the European Union’s
Horizon 2020 research and innovation programme and Netherlands, Belgium,
Germany, France, Italy, Austria, Hungary, Romania, Sweden and Israel.

Nguyen Hieu Thao  Oleg Soloviev  Jacques Noom and Michel Verhaegen

###### Abstract

This paper considers the problem of reconstructing an object with high-
resolution using several low-resolution images, which are degraded due to
nonuniform defocus effects caused by angular misalignment of the subpixel
motions. The new algorithm, indicated by the Superresolution And Nonuniform
Defocus Removal (SANDR) algorithm, simultaneously performs the nonuniform
defocus removal as well as the superresolution reconstruction. The SANDR
algorithm combines non-sequentially the nonuniform defocus removal method
recently developed by Thao _et al._ and the least squares approach for
subpixel image reconstruction. Hence, it inherits global convergence from its
two component techniques and avoids the typical error amplification of multi-
step optimization contributing to its robustness. Further, existing
acceleration techniques for optimization have been proposed that assure fast
convergence of the SANDR algorithm going from rate $\mathcal{O}(1/k)$ to
$\mathcal{O}(1/k^{2})$ compared to most existing superresolution (SR)
techniques using the gradient descent method. An extensive simulation study
evaluating the new SANDR algorithm has been conducted. As no algorithms are
available to address the combined problem, in this simulation study we
restrict the comparison of SANDR with other SR algorithms neglecting the
defocus aberrations. Even for this case the advantages of the SANDR algorithm
have been demonstrated.

###### Keywords:

Superresolution, Image reconstruction, Computational imaging, Deconvolution,
Inverse problems

##  1 Introduction

Image-based quality control is one of the important tools used during the
manufacturing and end quality checks in semiconductor HUANG20151 , automotive
Zhou2019 , and many other industries. For Industry 4.0, requiring a fully
automatized quality checks, the spatial resolution (size of the smallest
feature that can be inspected) is the key factor that affects the overall
efficiency and throughput of the control tool. For high-quality imaging
systems used in these tools, the spatial resolution is defined as the quotient
between the pixel size and the magnification, and thus for a higher
resolution, either a smaller pixel size or a larger magnification is required.
Larger magnification corresponds to smaller field of view (FOV in Fig. 1),
often it is desirable to have a smaller pixel size. However, there are
technological and design limits to the magnification and the smallest pixels
that can be manufactured and/or used in these tools and thus a computational
approach to increasing spatial resolution provides an interesting alternative.
The reconstruction of an object with high-resolution from several low-
resolution (LR) images capturing the object at subpixel-offset positions,
called the SuperResolution (SR) problem, has been studied for many decades
PelKerSch87 ; UrGro92 . A number of solution approaches have been proposed for
the SR problem, including direct methods KimBosVal90 ; KimSu93  and iterative
algorithms SauAll87 ; NguMilGol01 ; FarRobElaMil04 ; SroCriFlu07 ; TakKan08 ;
LiHuGaoTaoNin10 ; ZhaYuaSheLi11 ; ZhaLiShiLin11 ; LagGhaHakRag16 ;
WanLinDenAn17 ; HuaSunYanFanLin17 .

Superresolution reconstruction is possible if the LR images are registered for
different subpixel-offset positions of the object. In practice, shifting the
object at subpixel scale can be a major challenge to the SR problem and gives
rise to a number of important questions that need to be addressed. Camera
shake and motion blur induced by the shifts have been analyzed in BasBlaZis96
; KanMil13 . Inaccuracy of subpixel registration has been considered in
LeeKan03 ; TakMilProEla09 . Ideally, the shifting process should not cause any
variations in the object orientation with respect to the camera. However, this
is not always the case in practice and such deviations cause undesirable
deterioration of the data images and thus the reconstruction. Imprecise
displacements with respect to the optical axis would introduce defocus blurs
in the acquired images. More challenging, the shifting process can induce
rotational movements of the object causing nonuniform defocus effects in the
data. To the best of our knowledge, the latter challenge has not been
considered in the literature of SR by subpixel motions.

In this paper, we consider the problem of reconstructing a Superresolution
Image using Nonuniform Defocused images, called the SIND problem. Our
consideration was primarily motivated by the inspection of wafers in
semiconductor industry and the basic hypotheses are mainly inspired by its
practical context, but the resulting solution is also scalable for similar
applications of computer vision. As an alternative to shifting the object, LR
images can be registered using multiple cameras whose optical axes are
typically at different directions towards the object. This also results in
_nonuniform defocus blurs_ in the acquired data, and the SIND problem covers
this challenge as a special case with known and fixed blurs.

Solution approaches to the SIND problem should address three main tasks,
including estimation of nonuniform defocus models, removal of nonuniform
defocus effects, and reconstruction of an SR image. Assuming that the
nonuniform defocus models have been estimated, this paper is devoted to the
last two tasks. More specifically, we propose a new algorithm to
simultaneously perform both Superresolution reconstruction And Nonuniform
Defocus Removal (SANDR). The SANDR algorithm combines the nonuniform defocus
removal method recently developed in ThaOleJacMic21  and the least squares
approach VerVer07  for subpixel image reconstruction but not in a sequential
manner. Hence, it inherits global convergence from its two component
techniques and avoids the typical error amplification of multi-step
optimization contributing to its robustness. Further, existing acceleration
techniques for optimization BecTeb09  have been proposed that assure fast
convergence of the SANDR algorithm going from rate $\mathcal{O}(1/k)$ to
$\mathcal{O}(1/k^{2})$ compared to most existing SR techniques using the
gradient descent method, where $k$ is the number of iterations.

As, to our knowledge, no algorithms are available to address the SIND problem,
we demonstrate the advantages of the SANDR algorithm over other SR algorithms
neglecting the defocus aberrations, see Sect. 4.1. It is important to mention
that the Projected Gradient (PG) and the so-called Sequential Minimization
(SM) algorithms reported along with SANDR in the numerical section are also
considered for the SIND problem for the first time. Hence, comparing the SANDR
algorithm with them is not a goal of this paper.

##  2 Problem formulation

###  2.1 Superresolution by subpixel motions

a)
b)

Figure 1: a) Imaging scheme used in the problem formulation. An imaging lens
with magnification $Q$ creates image $o(x,y)$ of some planar object, which is
registered by a detector with pixel size $s$ to obtain $M$ sampled images
$o_{m}[i,j]$ (green, blue), with introduced subpixel offsets $v_{m}$ in each
of them. The superresolution problem is to restore $O[i,j]$ representing
$o(x,y)$ sampled with a finer grid (orange). b) Side view of the imaging
scheme with an example of misalignment of the object and detector planes
creating position-dependent defocus blur.

Let $o_{m}$ be the LR images, created by an imaging system, that sample image
$o$ of some planar object, see Fig. 1. Let each $o_{m}$ be registered with
some subpixel offset $v_{m}$, with coordinates expressed in pixels of $o_{m}$.
We have the following sampling:

$o_{m}\sim T_{v_{m}}(o),\;(m=1,2,\ldots,M),$ (1)

where $T_{v}$ denotes the translation by a vector $v$.

Let $O$ be the SR image to be reconstructed. The ratio between the sizes of
$O$ and $o_{m}$ is called the superresolution factor and denoted by $\tau$. In
this paper, the images are assumed to be square and the superresolution factor
is the same in both row and column directions for the sake of brevity. As the
shifts are measured in pixels of $o_{m}$, their coordinates with respect to
$O$ should be scaled up by $\tau$. Then we also have the following sampling:

$T_{\tau v_{m}}(O)\sim T_{v_{m}}(o),\;(m=1,2,\ldots,M).$ (2)

The combination of (1) and (2) leads to the following superresolution model:

$o_{m}\simeq D_{\tau}\circ T_{\tau v_{m}}(O),\;(m=1,2,\ldots,M),$ (3)

where $D_{\tau}$ is the downsampling operator with rate $\tau$, see Sect. 2.5
for the definition.

######  Remark 1 (external blurs)

The imaging model (3) can be extended as follows IraPel91 ; FarRobElaMil03 ;
FarRobElaMil03.2 ; ParParKan03 ; SroCriFlu07 ; TakKan08 ; LiHuGaoTaoNin10 ;
ZhaYuaSheLi11 ; ZhaLiShiLin11 ; LokSolSavVdo11 ; LagGhaHakRag16 ;
WanLinDenAn17 ; HuaSunYanFanLin17 :

$o_{m}\simeq D_{\tau}\circ H_{m}\circ T_{\tau v_{m}}(O),\;(m=1,2,\ldots,M),$

where $H_{m}$ denote the external blurs often modelled as isoplanatic
convolutions and assumed to be known. For our target application in wafer
inspection, external blurs are not so relevant and thus left for brevity
though they do not add major challenge to the problem under consideration.
Instead, we handle the more challenging anisoplanatic blurs induced by angular
misalignment of the shifts as detailed in the next section.

###  2.2 Nonuniform defocus effects

In practice, the subpixel shifts can be accomplished by moving either the
sample or the detector chip. In both cases, some angular misalignment can be
introduced, which can be difficult or costly (e.g., in terms of time/overall
throughput) to eliminate completely. Figure 1b shows an example of a
misaligned sample; a similar picture could be drawn for a misaligned detector,
where geometrical distortions might also appear.111In this paper, we consider
the geometrical distortions to be negligible compared to the position-
dependent blur. For simplicity, we do not discriminate between the object and
detector misalignment. Depending on the particular realisation and on the
optical magnification of the system, this might create presence of position-
dependent defocus blur in the image, which, as we show later, prevents the
direct application of existing superresolution algorithms.

We consider the challenge that the displacement process induces undesirable
rotational movements of the object and the acquired images are degraded by
nonuniform defocus blurs. In this case, the theoretical LR images $o_{m}$ and
the measured ones $i_{m}$ are related by

$i_{m}\simeq B_{m}(o_{m}),\;(m=1,2,\ldots,M),$ (4)

where the blur operators $B_{m}$ will be detailed shortly.

Let $o_{m}$ situate in $N$ defocus zones denoted by $D_{n}$
$(n=1,2,\ldots,N)$, for each of which the Point Spread Function (PSF) is
modelled using the Fourier transform Goo05 :

$p_{n}=\left|\mathcal{F}\left(A\cdot{\rm e}^{{\rm
j}d_{n}Z_{2}^{0}}\right)\right|^{2}\quad(n=1,2,\ldots,N),$ (5)

where the amplitude, product and square operations are elementwise,
$\mathcal{F}$ is the two-dimensional Fourier transform, $A$ is the binary mask
representing the camera aperture,222We assume that the diffraction-limited PSF
is not resolved by the camera pixels. ${\rm j}=\sqrt{-1}$ is the imaginary
unit, $d_{n}$ is the (directional) distance from $D_{n}$ to the focal plane,
and $Z_{2}^{0}$ is the Zernike polynomial of order two and azimuthal frequency
zero.

We make use of the following model of nonuniform defocus blurs, whose physical
relevance has been demonstrated, e.g., in ThaOleJacMic21 :

$B_{m}(o_{m})=\sum_{n=1}^{N}\left(\mu_{mn}\cdot
o_{m}\right)*p_{n},\;(m=1,2,\ldots,M),$ (6)

where $*$ is the two-dimensional convolution, $\mu_{mn}$ are the mask
functions of $o_{m}$ defined by: for $n=1,2,\ldots,N$,

$\mu_{mn}[i,j]=\begin{cases}1&\text{if }\;o_{m}[i,j]\in D_{n},\\\
0&\text{otherwise}.\end{cases}$

This paper considers planar objects and $d_{n}$ take the following form
(ThaOleJacMic21, , Sect. IIC,):

$d_{n}=d(n_{0}-n),\quad(n=1,2,\ldots,N),$ (7)

where $d$ is the Depth of Focus (DoF) and $n_{0}$ is the _focal position_. The
number of defocus zones $N$ defines the degree of defocus in an image and can
be different for each $o_{m}$. For brevity, it is taken the same in this
paper.

The combination of (5), (6) and (7) yields the following blur model:
($m=1,2,\ldots,M$)

$B_{m}(o_{m})=\sum_{n=1}^{N}\left(\mu_{mn}\cdot
o_{m}\right)*\left|\mathcal{F}\left(A\cdot{\rm e}^{{\rm
j}d(n_{0}-n)Z_{2}^{0}}\right)\right|^{2}.$ (8)

###  2.3 The SIND problem

We consider the problem of reconstructing the Superresolution Image $O$ from
the Nonuniform Defocused images $i_{m}$ according to the relations (3) and
(4), where $B_{m}$ are given by (8). It is referred to as the SIND problem.

###  2.4 Optimization formulations

Combining (3) and (4) yields the imaging model:

$i_{m}=B_{m}\circ D_{\tau}\circ T_{\tau
v_{m}}(O)+w_{m},\quad(m=1,2,\ldots,M),$ (9)

where $B_{m}$ are given by (8), and $w_{m}$ represent the discrepancies
between the theoretical and the measured data, e.g., due to noise and model
deviations.

In this paper, $w_{m}$ are assumed to be independent zero-mean random
variables with jointly Gaussian-distributed entries.333This assumption is
ubiquitous and it does not rule out the case of Poisson noise as the latter
can be well approximated by a Gaussian distribution in view of the central
limit theorem provided that the image is registered with a sufficiently large
number of photon counts. For each $m=1,2,\ldots,M$, let $\mathcal{W}_{m}$ be
the covariance matrix of $\mathtt{Vec}\left(w_{m}\right)$, where
$\mathtt{Vec}$ denotes the vectorization operator. Then the maximum-likelihood
approach, e.g., (VerVer07, , Sect. 4.5.5,), applied to (9) leads to the
following minimization problem:

$\min_{O}\;\quad f(O)+\mathcal{G}(O),$ (10)

where $\mathcal{G}$ is the regularization capturing the physical attributes of
$O$ (see Sect. 3.1), and $f$ represents the data fidelity given by

$\displaystyle f(O)$
$\displaystyle=\sum_{m=1}^{M}R_{m}(O)^{T}\mathcal{W}_{m}^{-1}\,R_{m}(O),$ (11)

where $R_{m}$ ($m=1,2,\ldots,M$) are the fitting residual errors for the
(blurred) LR images:

$\displaystyle R_{m}(O)$ $\displaystyle=\mathtt{Vec}\left(B_{m}\circ
D_{\tau}\circ T_{\tau v_{m}}(O)-i_{m}\right).$

######  Remark 2 (sequential optimization)

The residual error in (3) is mainly due to the inaccuracy of the subpixel
shifts while the one in (4) is more related to measurement noise and model
deviations of (8). When the latter is less severe than the former,444This is
relevant to wafer inspection, where the camera is high-quality while
inexactness of the subpixel shifts poses the major challenge. one can also
address (4) and (3) sequentially via the following two-step optimization:

$\min_{O}\;\sum_{m=1}^{M}S_{m}(O)^{T}{\Sigma}_{m}^{-1}\,S_{m}(O)+\mathcal{G}(O),$
(12)

where for $m=1,2,\ldots,M$,

$\displaystyle S_{m}(O)=\mathtt{Vec}\left(D_{\tau}\circ T_{\tau
v_{m}}(O)-\widehat{o}_{m}\right),$
$\displaystyle\widehat{o}_{m}\in\arg\min_{o_{m}}Q_{m}(o_{m})^{T}\mathcal{E}_{m}^{-1}\,Q_{m}(o_{m})+\mathcal{H}\left(o_{m}\right),$
(13) $\displaystyle Q_{m}(o_{m})=\mathtt{Vec}\left(B_{m}(o_{m})-i_{m}\right).$

In the above, $\Sigma_{m}$ and $\mathcal{E}_{m}$ are respectively the
covariance matrices representing the noise in (4) and (3), and $\mathcal{H}$
is the regularization capturing the physical attributes of $o_{m}$.
Sequentially minimizing (13) and (12) gives rise to the so-called SM algorithm
(see Sect. 3.2), which suffers the typical error amplification of multi-step
optimization compared to the proposed solution method for solving (10), see
Sect. 4.3&4.4.

###  2.5 Downsampling operators

The downsampling operator with integer rate $\tau$ is given by

$D_{\tau}([u])=\frac{1}{\tau^{2}}\,\mathtt{conv}_{\tau}\left(u,\mathbf{1}_{\tau}\right),\quad(\forall\,u),$
(14)

where $\mathtt{conv}_{\tau}$ denotes the bivariate convolution operation with
striking sizes $\tau\times\tau$,555The terminology is standard in the field of
convolutional neural networks. and $\mathbf{1}_{\tau}$ is the all-ones matrix
of size $\tau\times\tau$. The striking sizes define the size reduction in row
and column directions. $D_{\tau}$ produces only the average intensity value of
every $\tau\times\tau$-block and hence it is not invertible without additional
information of $u$.

##  3 Solution approaches

Solution approaches to the SIND problem should address three main tasks,
including estimation of nonuniform defocus models, removal of nonuniform
defocus effects, and reconstruction of an SR image. Assuming that the
nonuniform defocus models have been estimated, this paper is devoted to the
last two tasks. We first discuss regularization schemes for the SIND problem.

###  3.1 Regularization functions

SR methods often minimize a cost function consisting of data fidelity and
regularization IraPel91 ; NguMilGol01 ; FarRobElaMil03 ; FarRobElaMil03.2 ;
ParParKan03 ; FarRobElaMil04 ; SroCriFlu07 ; TakKan08 ; LiHuGaoTaoNin10 ;
ZhaYuaSheLi11 ; ZhaLiShiLin11 ; LokSolSavVdo11 ; LagGhaHakRag16 ;
WanLinDenAn17 ; HuaSunYanFanLin17 . Data fidelity is typically a norm of the
residual between the theoretical and the measured data while regularization is
driven by the _a priori_ known physical attributes of the solution. The latter
pertains to each particular application and is the main difference between
existing SR techniques. Total variation and Tikhonov regularization were
considered in, e.g., NguMilGol01 ; FarRobElaMil03 . The Bilateral Total
Variation (BTV) was introduced in FarRobElaMil03  and later adapted in
FarRobElaMil04 ; LiHuGaoTaoNin10 ; ZhaLiShiLin11 ; LagGhaHakRag16 ;
WanLinDenAn17 . In LagGhaHakRag16  BTV was used in combination with the
Laplace operator while in LiHuGaoTaoNin10  it was used in combination with
another regularization to enhance the consistence of the gradient variation.

In this paper, the images are assumed to have intensities in $[0,1]$, and the
set of matrices satisfying this constraint is denoted by $\Omega$. This
constraint is easy to handle, but essential for the success of our proposed
algorithms, where acceleration optimization mechanisms are exploited. Its
effectiveness has been widely known in the literature of deconvolution, see,
e.g., WilSolPozVdoVer17 ; ThaOleJacMic21 . There are several approaches to
this constraint, e.g., the penalty approaches using the associated distance
function or its square. In this paper, we make use of the _indicator function_
VA :

$\iota_{\Omega}(x)=\begin{cases}0&\text{ if }x\in\Omega,\\\ \infty&\text{
otherwise}.\end{cases}$ (15)

In our simulation results, this constraint is a precise regularization and
hence its advantages over the other schemes are clearly observed, see Sect.
4.1.

###  3.2 The proposed algorithms

In view of Remark 2, the SIND problem can be addressed by solving (13) and
(12) sequentially. For each $m=1,2,\ldots,M$, (13) is the single-frame
nonuniform defocus removal problem recently studied in ThaOleJacMic21 . Hence,
it can be solved by the algorithm proposed in that paper, where its challenges
including the typical ill-posedness were also discussed and global convergence
of the proposed algorithm was also established. The main challenge of (12) is
that the downsampling operator $D_{\tau}$ is not invertible, in particular,
closed-form solutions for it are not available. We propose to apply the
regularization (15) and make use of the fast proximal gradient method
introduced in BecTeb09 , often known as FISTA, for solving (12).

The algorithm resulted from this sequential approach will be referred to as
the Sequential Minimization (SM) algorithm for the SIND problem. However, we
chose to skip its details for the sake of brevity. The main advantages of the
SM algorithm include its simplicity and the parallelism of (13) while its
major disadvantage is the typical error amplification of multi-step
optimization.

To overcome the drawback of SM, we next propose a new algorithm to
simultaneously handle both Superresolution reconstruction And Nonuniform
Defocus Removal (SANDR). The SANDR algorithm combines the nonuniform defocus
removal method developed in ThaOleJacMic21  and the least squares approach
VerVer07  for subpixel image reconstruction but not in a sequential manner.
Hence, it inherits global convergence from its two component techniques and
avoids the typical error amplification of multi-step optimization contributing
to its robustness, see Sects. 4.3&4.4. Making use of the acceleration
techniques for optimization of FISTA assures fast convergence of the SANDR
algorithm going from rate $\mathcal{O}(1/k)$ to $\mathcal{O}(1/k^{2})$
compared to most existing SR techniques using the gradient descent method.

For simplicity, the noise covariance matrices $\mathcal{W}_{m}$ in (11) are
taken to be the identity matrix in the sequel. The repetitive term
$(m=1,2,\ldots,M)$ following the subscript $m$ will be omitted for brevity.

In the sequel, $U_{\tau}$ will denote a right inverse of the downsampling
operator $D_{\tau}$ defined in (14), i.e., $D_{\tau}\circ
U_{\tau}=\operatorname{Id}$, the identity mapping.666$U_{\tau}$ is not unique
and in general $D_{\tau}\circ U_{\tau}\neq U_{\tau}\circ D_{\tau}$. $U_{\tau}$
can be understood as a numerical upsampling operator, and in our numerical
results, it is taken to be the interpolation with block constant values.
Recall that the translation by a vector with integer coordinates
$(v_{x},v_{y})$, is given by

$T_{(v_{x},v_{y})}(u)(r,c)=u(r-v_{x},c-v_{y}),\quad(\forall\,u),$ (16)

where $(r,c)$ are the row-column coordinates of the pixels.

######  Algorithm 1 (the SANDR algorithm)

_Input:_ $i_{m}$ – LR images, $B_{m}$ – blur operators, $\lambda$ – stepsize,
$t^{(0)}$ – initial acceleration stepsize, $K$ – number of iterations, and
$\varepsilon>0$.
_Initialization:_ $X^{(0)}=O^{(0)}=\frac{1}{M}\sum_{m=1}^{M}T_{(-\tau
v_{m})}(U_{\tau}(i_{m}))$.
_Iteration process_ : given $X^{(k)}$, $O^{(k)}$, $t^{(k)}$

$\displaystyle G_{m}^{(k)}=T_{(-\tau v_{m})}\circ U_{\tau}\left(\nabla
f_{m}\left(D_{\tau}\circ T_{\tau v_{m}}\left(O^{(k)}\right)\right)\right),$
$\displaystyle X_{m}^{(k+1)}=P_{\Omega}\left(O^{(k)}-\lambda
G_{m}^{(k)}\right),$ $\displaystyle
X^{(k+1)}=\frac{1}{M}\sum_{m=1}^{M}X_{m}^{(k+1)},$ $\displaystyle
t^{(k+1)}=\frac{1+\sqrt{1+4{t^{(k)}}^{2}}}{2},$ $\displaystyle
O^{(k+1)}=X^{(k+1)}+\frac{t^{(k)}-1}{t^{(k+1)}}\left(X^{(k+1)}-X^{(k)}\right).$

_Stopping criteria_ : $k>K$ or

$\sum_{m=1}^{M}\left\|G_{m}^{(k)}\right\|>\sum_{m=1}^{M}\left\|G_{m}^{(k-1)}\right\|+\varepsilon.$
(17)

_Output:_ $\widehat{O}=P_{\Omega}\left(O^{(\mathtt{end})}\right)$.

In Algorithm 1, $P_{\Omega}$ is the projection operator associated with
$\Omega$ and the functions $f_{m}$ are given by

$f_{m}(x)=\frac{1}{2}\left\|B_{m}(x)-i_{m}\right\|^{2},\;(m=1,2,\ldots,M).$

##  4 Numerical simulations

As explained in Sect. 2.2, a higher degree of defocus in an (LR) image
corresponds to a larger number of defocus zones and smaller supports (nonzero
entries) of the mask functions and vice versa. To simplify simulation of
random defocus levels in LR images, we chose to fix these parameters, but
consider the DoF $d$ in (7) as the single parameter quantifying the defocus in
each image, called the blur coefficient of the image in the sequel. It is
important to mention that our choice for convenience does not contradict the
fact that DoF is a fixed physical parameter of the camera because
underestimation of DoF does not introduce model deviations.777It only costs
computational time as the number of defocus zones increases accordingly. The
larger the blur coefficient is, the more the defocus blur in the image.

Simulation data is generated according to the forward imaging model (9).
Except for the analysis regarding the number of input images in Sect. 4.5,
each data set consists of four images corresponding to the shift vectors
$v_{1}=(0,0)$, $v_{2}=(1/2,0)$, $v_{3}=(0,1/2)$ and $v_{4}=(1/2,1/2)$. Half of
the images contain defocus blur varying in the vertical direction and half in
the horizontal direction. Unless otherwise specified, the common parameters
are as in Table 1.

Table 1: Parameters used in Sect. 4. $M$ is the number of LR images, $\mu$ – size of the supports of mask functions in pixel rows/columns, $\rho$ – PSF size (square), $\tau$ – SR factor, $\lambda$ – stepsize, and $t^{(0)}$ – initial acceleration stepsize. Parameter | $M$ | $\mu$ | $\rho$ | $\tau$ | $\lambda$ | $t^{(0)}$
---|---|---|---|---|---|---
Value | 4 | 3 | 11 | 2 | 1 | 1

Except for the noise analysis in Sect. 4.4, the data is corrupted with Poisson
noise using the MATLAB function $\mathtt{imnoise}$. The quality of SR
reconstruction is measured by the Root Mean Square (RMS) error of the restored
SR image relative to the ideal one:
${\|\widehat{O}-O\|}\big{/}{\left\|O\right\|}$. The stopping criterion (17) is
not implemented as it is not so relevant for simulations.

As no algorithms are available to address the SIND problem, we can only
demonstrate its advantages over other SR methods neglecting the defocus
effects. It is important to mention that the Projected Gradient (PG) and the
Sequential Minimization (SM) algorithms are also first considered for the SIND
problem, and hence comparing the SANDR algorithm with them is not a goal of
this section. Instead, their own advantages and disadvantages in various
problem settings will be of our primary interest.

###  4.1 Comparison to known SR methods

Figure 2: LR image (left) and the ideal SR (right). The ROIs are shown in Fig.
3 for visual comparison of different SR methods.

Figure 3: SR images obtained by L1-BTV, L1-BTV-L, L2-BTV and SANDR are shown
together with an LR image and the ideal SR. Only the ROIs are shown for
clarity. The relative RMS error of each ROI is also reported. The SANDR
algorithm outperforms the other SR methods.

Most existing SR methods minimize a cost function consisting of data fidelity
and regularization using the gradient descent method FarRobElaMil04 ;
SroCriFlu07 ; TakKan08 ; LiHuGaoTaoNin10 ; ZhaYuaSheLi11 ; ZhaLiShiLin11 ;
LagGhaHakRag16 ; WanLinDenAn17 ; HuaSunYanFanLin17 . Data fidelity is
typically the (weighted) $L_{p}$-norm ($1\leq p\leq 2$) of the residual
between the theoretical and the measured data while regularization is driven
by the _a priori_ known physical attributes of the solution, see Sect. 3.1. In
this section, we compare the SANDR algorithm with three existing SR methods
minimizing (1) the $L_{1}$-norm with bilateral total variation (L1-BTV)
FarRobElaMil04 , (2) the $L_{1}$-norm with BTV and Laplace operator (L1-BTV-L)
LagGhaHakRag16 , and (3) the $L_{2}$-norm with BTV (L2-BTV) WanLinDenAn17 .
Each iteration of the algorithms is additionally followed by a projection on
the constraint $\Omega$ to improve their performance, especially in terms of
stability. Note that without defocus effects, the SM and the SANDR algorithms
coincide.

Figure 4: Relative RMS errors of the SR images obtained by L1-BTV, L1-BTV-L,
L2-BTV and SANDR are shown in iterations. The SANDR algorithm is superior to
the other methods in both convergence speed and accuracy.

Figure 2 shows an LR image (left) and the ideal SR (right). The SR images
obtained by L1-BTV, L1-BTV-L, L2-BTV and SANDR are shown in Fig. 3 together
with an LR image and the ideal SR. Only the ROIs are shown for clarity. The
relative RMS error of each ROI is also reported. The SANDR algorithm clearly
outperforms the other methods both visually and in terms of RMS errors. This
is further explained in Fig. 4, where the RMS errors are shown in iterations.
The SANDR algorithm is far superior to the others in both convergence speed
and accuracy. Faster convergence is due to the acceleration feature of SANDR
while higher accuracy can be explained by the fact that $\Omega$ is a precise
regularization in this simulation problem. Ripple behaviours of L1-BTV and
L1-BTV-L in Fig. 4 can be explained by the step-size being larger than the
distance from the iteration to a local minimum. This phenomenon is more likely
to happen to $L_{1}$-norm cost functions since their gradient includes the
sign function, which does not depend on the residual gap of the current
iteration.888An advantage of $L_{1}$-norm cost functions is their potential to
suppress outliers. Gradually decreasing the stepsize is a possible remedy for
this issue, however, we chose not to distract the reader further in that
direction because there does not exist a unified recipe for such tasks while
the methods are not applicable to the SIND problem.

###  4.2 Convergence properties

Figure 5: LR image with defocus varying in the vertical direction (left) and
the unblurred one (right). Figure 6: RMS errors of the SR images obtained by
PG, SM and SANDR are shown in iterations. The algorithms exhibit convergence
properties and without acceleration, PG (red) converges slower than SM (black)
and SANDR (blue). The RMS error with 150 PG iterations is 0.79%, corresponding
to about 30 iterations of SANDR.

In this section, we demonstrate convergence properties of the SANDR algorithm
along with the PG and SM methods. We consider LR images of size $330\times
330$ pixels with $110$ defocus zones and blur coefficients randomly taken in
the interval $[0.001,0.06]$. The other parameters are as in Table 1. One of
the LR images with defocus effects varying in the vertical direction and its
unblurred version are shown in Fig. 5.

In Fig. 6 the RMS errors of the SR images obtained by PG, SM and SANDR are
shown in iterations. The algorithms exhibit convergence properties and without
acceleration, PG (red) converges slower than SM (black) and SANDR (blue). The
RMS error with 150 iterations of PG is 0.79%, corresponding to about 30
iterations of SANDR. In Fig. 7 the ROIs of the SR images obtained by 5, 15, 50
and 150 iterations of PG, SM and SANDR are shown in comparison with the ones
of an LR image and the ideal SR. The relative RMS error of each ROI is also
shown. Note that the RMS error of each unblurred LR image is around 13.83%.

a)
b)
c)
d)
e)

Figure 7: ROIs of the SR images obtained by a) 5, b) 15, c) 50, d) 150
iterations of PG, SM and SANDR are shown in comparison with e) the ones of an
LR image and the ideal SR. The RMS error of each ROI is also shown. The RMS
error of each unblurred LR image is around 13.83%.

###  4.3 Solvability analysis

As explained at the beginning of Sect. 4, blur coefficients quantify the
defocus in the LR images. The larger they are, the more challenging the SIND
problem is. In this section, we analyse the solvability of PG, SM and SANDR
with respect to this parameter. For each experiment, four blur coefficients of
the LR images are randomly taken in the interval $[0.001,d_{\mathtt{max}}]$
with $d_{\mathtt{max}}$ ranging from 0.06 to 0.3. The other parameters are as
in Table 2.999It is a trade-off between the computational complexity (number
of iterations) and the restoration accuracy, in view of Figs. 4 and 6 we chose
to run $50$ iterations for each experiment.

Table 2: Parameters used in Sect. 4.3. $d_{\mathtt{max}}$ is the upper bound of (random) blur coefficients, $N$ – # defocus zones, $n_{0}$ – focal position, $(l,w)$ – size of LR images, and $K$ – # iterations. Parameter | $d_{\mathtt{max}}$ | $N$ | $n_{0}$ | $l=w$ | $K$
---|---|---|---|---|---
Value | 0.06 – 0.3 | 55 | 28 | 165 | 50
Figure 8: Solvability analysis of PG, SM and SANDR with respect to the
distortion level of data images quantified by the blur coefficients. One
hundred experiments are reported for each value of $d_{\mathtt{max}}$ ranging
from 0.06 to 0.3. The restoration errors increase for larger
$d_{\mathtt{max}}$. SM works best for $d_{\mathtt{max}}$ up to 0.06, but it
quickly becomes problematic for $d_{\mathtt{max}}$ from 0.12 due to its
sequential optimization. SANDR is effective for blur coefficients up to 0.3
and outperforms PG and SM (for $d_{\mathtt{max}}\geq 0.12$) in both accuracy
and consistency. The superiority becomes more significant for larger values of
$d_{\mathtt{max}}$.

For each value of $d_{\mathtt{max}}$, one hundred experiments with PG, SM and
SANDR are reported in Fig. 8, where the RMS errors of the obtained SR images
with respect to the ideal one are presented. The restoration errors increase
for larger values of $d_{\mathtt{max}}$. SM works best for $d_{\mathtt{max}}$
up to 0.06, but it quickly becomes problematic for $d_{\mathtt{max}}$ from
0.12 due to its sequential optimization. SANDR is effective for blur
coefficients up to 0.3. It outperforms PG and SM (for $d_{\mathtt{max}}\geq
0.12$) in both accuracy and consistency, and the superiority becomes more
significant for larger $d_{\mathtt{max}}$. Higher accuracy is reflected by its
smaller average restoration errors while more consistency is indicated by its
smaller variances of the errors. To visualize the _blur coefficient_
parameter, the PSFs for the 30th and 50th defocus zones (counted from the
focal position) with blur coefficient 0.06 are shown in Fig. 9. Recall that
the distortion level of an image is proportional to the product of the blur
coefficient and the zone position in view of (7).

Figure 9: PSFs for the 30th (left) and 50th (right) defocus zones with blur
coefficient 0.06.

###  4.4 Noise analysis

We analyse the influence of noise on the performance of PG, SM and SANDR. Five
levels of Gaussian noise ranging from 45 to 65 dB (decibels) are considered.
Recall that the signal-to-noise ratio (SNR) in decibels is defined by:
$\text{SNR}=10\ln\left({P}/{P_{0}}\right)$, where $P$ and $P_{0}$ are the
powers of the signal and noise, respectively. To visualize the noise, an LR
image with SNR 45 dB is shown in Fig. 12 together with its residual relative
to the noiseless one.

Figure 10: Noise analysis of PG, SM and SANDR. One hundred experiments are
reported for each SNR from 45 to 65 dB. The reconstruction is more accurate
for higher SNR. SM works best for SNR from 55 dB, but it quickly becomes
problematic for SNR decreasing from 50 dB. Since SM and SANDR are acceleration
variants of PG, they are less robust than PG. Figure 11: Experiments similar
to those reported in Fig. 10 but with $d_{\mathtt{max}}=0.12$ instead of
$0.09$ show that SM deteriorates much faster than PG and SANDR for larger blur
coefficients.

For each SNR, one hundred experiments with $d_{\mathtt{max}}=0.09$ and the
other parameters as in Table 2 are reported in Fig. 10, where the RMS errors
of the SR images obtained by PG, SM and SANDR are presented. The
reconstruction is more accurate for higher SNR. SM works best for SNR from 55
dB, but it quickly becomes problematic for SNR decreasing from 50 dB. Its less
robustness against noise compared to SANDR is due to its two-step
optimization, see also Sect. 4.3. It is not a surprise that SM and SANDR are
less robust than PG because the former are acceleration variants of the latter
and there is a typical trade-off between robustness and convergence speed. In
view of Fig. 10, it is worth thinking about PG for the SIND problem with SNR
below 45 dB, but for higher SNR it is outperformed by the others. It is
important to recall that the conclusions drawn for the SM algorithm from Fig.
10 are valid only for $d_{\mathtt{max}}$ up to 0.09, which seems to be a limit
for it, see also Fig. 8. To demonstrate this point, we do similar experiments
but with slightly larger blur coefficients, $d_{\mathtt{max}}=0.12$ in place
of $0.09$. The results are summarized in Fig. 11, where SM deteriorates much
more than PG and SANDR in comparison with Fig. 10.

Figure 12: A noisy LR image (left) with SNR 45 dB and its residual (right)
with respect to the noiseless one.

###  4.5 Number of input images

The major practical challenge of SR by subpixel motions is to perform the
shifts accurately. Let us suppose that we are able to perform shifts at scale
$1/\tau$ pixel, where $1<\tau\in\mathbb{Z}$.101010In this study, inaccuracy of
subpixel registration is subsumed in noise. Then there are at most $\tau^{2}$
LR images and one cannot expect to gain a SR factor greater than $\tau$. In
this section, we briefly study the influence of the number of input images on
the quality of SR. We consider $\tau=4$ and construct the SR image using 2, 4,
8 and 16 LR images, respectively. In this experiment, $d_{\mathtt{max}}=0.06$
and the other parameters are as in Table 2.

The numerical results are summarized in Fig. 13, where only the ROIs and their
RMS errors are shown for brevity. It is clear that more input images result in
higher quality of the SR and the observation is consistent for PG, SM and
SANDR. SR images obtained with two LR images (the first row) already shows
improvement even in comparison with the unblurred LR images (the second in the
last row).

Figure 13: Influence of the number of input images on the SR images obtained
by PG, SM and SANDR. ROIs of the SR images obtained with 2, 4, 8 and 16 LR
images are shown together with their RMS errors. More input images result in
higher quality of SR.

###  4.6 Image cropping

Cropping the data images would introduce deviations to the imaging model (9).
This issue does not arise in the previous sections since the simulation object
there has almost constant intensity near the boundary. In this section, we
study the influence of image cropping on the performance of the PG, SM and
SANDR algorithms.

Table 3: Parameters used in Sect. 4.6. $d_{\mathtt{max}}$ is the upper bound of blur coefficients, $N$ – # defocus zones, $n_{0}$ – focal position, $(l,w)$ – size of LR images, and $K$ – # iterations. Parameter | $d_{\mathtt{max}}$ | $N$ | $n_{0}$ | $l=w$ | $K$
---|---|---|---|---|---
Value | 0.09 | 85 | 43 | 255 | 50
Figure 14: SR images obtained by PG, SM and SANDR are shown together with an
LR image (bottom left) and the ideal SR (bottom right). The restoration error
is smaller in the central region and becomes larger towards the boundary. The
SANDR algorithm is the most effective for this problem while the SM algorithm
is problematic due to high level of defocus effect. The SANDR and SM
algorithms suffer more boundary effects than PG. Figure 15: The ROIs marked in
Fig. 14 are shown for a better inspection of finer details. Their RMS errors
with respect to the ideal SR are also reported.

Four images are generated according to (9) with the parameters as in Table 3.
They are then windowed using Butterworth function to yield the LR images, one
of which is shown at the bottom left of Fig. 14. The SR images obtained by PG,
SM and SANDR are shown in Fig. 14 in comparison with an LR image and the ideal
SR. The SANDR algorithm is the most effective for this problem while the SM
algorithm is problematic due to high level of defocus effect as discussed in
Sect. 4.3. The SM and SANDR algorithms suffer more boundary effects than PG
since the former are more sensitive to noise than the latter as analyzed in
Sect. 4.4. To reduce the restoration errors near the boundary, the Butterworth
filter also need to be applied to every iteration update of the algorithms.
The RMS errors are computed for the central regions with 90% in radius of the
images. The reconstruction error is smaller in the central region and
increases towards the boundary. The ROIs are zoomed out in Fig. 15 for a
better inspection of finer details. The RMS error of each ROI is also
reported.

##  5 Concluding Remarks

We have investigated the problem of constructing an object with high-
resolution using several nonuniform defocused images, called the SIND problem.
Nonuniform defocus effects can arise in both standard techniques of data
registration, including the use of multiple cameras and moving the object.
However, the SIND problem has not been studied before. We have proposed the
efficient algorithm for SIND, called the SANDR algorithm, that can process
both subpixel image reconstruction and nonuniform defocus removal
simultaneously. Important theoretical and practical aspects of the SANDR
algorithm have been analyzed, including its global convergence, solvability,
noise robustness, dependence on the number of LR images, and sensitivity to
model deviations due to image croping. We have demonstrated advantages of the
SANDR algorithm over a number of existing superresolution methods without
considering defocus effects because the latter cannot handle this additional
challenge. Our consideration was primarily motivated by the inspection of
wafers in semiconductor industry, but the SANDR algorithm can be scalable for
similar applications of computer vision in Industry 4.0.

Funding. This project has received funding from the ECSEL Joint Undertaking
(JU) under grant agreement No. 826589. The JU receives support from the
European Union’s Horizon 2020 research and innovation programme and
Netherlands, Belgium, Germany, France, Italy, Austria, Hungary, Romania,
Sweden and Israel.

## References

  * [1] B. Bascle, A. Blake, and A. Zisserman.  Motion deblurring and super-resolution from an image sequence.  In Bernard Buxton and Roberto Cipolla, editors, Computer Vision — ECCV ’96, pages 571–582, Berlin, Heidelberg, 1996. Springer Berlin Heidelberg.
  * [2] A. Beck and M. Teboulle.  A fast iterative shrinkage-thresholding algorithm for linear inverse problems.  SIAM J. Imaging Sci., 2(1):183–202, 2009.
  * [3] S. Farsiu, D. Robinson, M. Elad, and P. Milanfar.  Fast and robust super-resolution.  In Proceedings 2003 International Conference on Image Processing (Cat. No.03CH37429), volume 2, pages II–291, 2003.
  * [4] S. Farsiu, M. D. Robinson, M. Elad, and P. Milanfar.  Fast and robust multiframe super resolution.  IEEE Transactions on Image Processing, 13(10):1327–1344, 2004.
  * [5] Sina Farsiu, Dirk Robinson, Michael Elad, and Peyman Milanfar.  Robust shift and add approach to superresolution.  In Andrew G. Tescher, editor, Applications of Digital Image Processing XXVI, volume 5203, pages 121 – 130. International Society for Optics and Photonics, SPIE, 2003.
  * [6] J. W. Goodman.  Introduction to Fourier Optics.  Roberts & Company Publishers, 2005.
  * [7] Hidenori Takeshima and Toshimitsu Kaneko.  Image registration using subpixel-shifted images for super-resolution.  In 2008 15th IEEE International Conference on Image Processing, pages 2404–2407, 2008.
  * [8] S. Huang, J. Sun, Y. Yang, Y. Fang, and P. Lin.  Multi-frame super-resolution reconstruction based on gradient vector flow hybrid field.  IEEE Access, 5:21669–21683, 2017.
  * [9] Szu-Hao Huang and Ying-Cheng Pan.  Automated visual inspection in the semiconductor industry: A survey.  Computers in Industry, 66:1–10, 2015.
  * [10] Michal Irani and Shmuel Peleg.  Improving resolution by image registration.  CVGIP: Graphical Models and Image Processing, 53(3):231–239, 1991\.
  * [11] A. V. Kanaev and C. W. Miller.  Multi-frame super-resolution algorithm for complex motion patterns.  Opt. Express, 21(17):19850–19866, 2013.
  * [12] S.P. Kim, N.K. Bose, and H.M. Valenzuela.  Recursive reconstruction of high resolution image from noisy undersampled multiframes.  IEEE Transactions on Acoustics, Speech, and Signal Processing, 38(6):1013–1027, 1990.
  * [13] S.P. Kim and W.-Y. Su.  Recursive high-resolution reconstruction of blurred multiframe images.  IEEE Trans Image Process., 2(4):534–542, 1993.
  * [14] Amine Laghrib, Abdelghani Ghazdali, Abdelilah Hakim, and Said Raghay.  A multi-frame super-resolution using diffusion registration and a nonlocal variational image restoration.  Computers and Mathematics with Applications, 72(9):2535–2548, 2016\.
  * [15] Eun Sil Lee and Moon Gi Kang.  Regularized adaptive high-resolution image reconstruction considering inaccurate subpixel registration.  IEEE Transactions on Image Processing, 12(7):826–837, 2003.
  * [16] Xuelong Li, Yanting Hu, Xinbo Gao, Dacheng Tao, and Beijia Ning.  A multi-frame image super-resolution method.  Signal Processing, 90(2):405–414, 2010.
  * [17] Mikhail Loktev, Oleg Soloviev, Svyatoslav Savenko, and Gleb Vdovin.  Speckle imaging through turbulent atmosphere based on adaptable pupil segmentation.  Opt. Lett., 36(14):2656–2658, 2011.
  * [18] Nhat Nguyen, P. Milanfar, and G. Golub.  A computationally efficient superresolution image reconstruction algorithm.  IEEE Transactions on Image Processing, 10(4):573–583, 2001.
  * [19] Shmuel Peleg, Danny Keren, and Limor Schweitzer.  Improving image resolution using subpixel motion.  Pattern Recognition Letters, 5(3):223–226, 1987.
  * [20] R. T. Rockafellar and R. J. Wets.  Variational Analysis.  Grundlehren Math. Wiss. Springer-Verlag, Berlin, 1998.
  * [21] K. Sauer and J. Allebach.  Iterative reconstruction of bandlimited images from nonuniformly spaced samples.  IEEE Transactions on Circuits and Systems, 34(12):1497–1506, 1987\.
  * [22] F. Sroubek, G. Cristobal, and J. Flusser.  A unified approach to superresolution and multichannel blind deconvolution.  IEEE Transactions on Image Processing, 16(9):2322–2332, 2007.
  * [23] Sung Cheol Park, Min Kyu Park, and Moon Gi Kang.  Super-resolution image reconstruction: a technical overview.  IEEE Signal Processing Magazine, 20(3):21–36, 2003.
  * [24] Hiroyuki Takeda, Peyman Milanfar, Matan Protter, and Michael Elad.  Super-resolution without explicit subpixel motion estimation.  IEEE Transactions on Image Processing, 18(9):1958–1975, 2009.
  * [25] Nguyen Hieu Thao, Oleg Soloviev, Jacques Noom, and Michel Verhaegen.  Nonuniform defocus removal for image classification.  Manuscript submitted to IEEE Transactions on Image Processing in February 2021.
  * [26] Hanoch Ur and Daniel Gross.  Improved resolution from subpixel shifted pictures.  CVGIP: Graphical Models and Image Processing, 54(2):181–186, 1992\.
  * [27] Michel Verhaegen and Vincent Verdult.  Filtering and System Identification: A Least Squares Approach.  Cambridge University Press, 2007.
  * [28] Longguang Wang, Zaiping Lin, Xinpu Deng, and W. An.  Multi-frame image super-resolution with fast upscaling technique.  arXiv: Computer Vision and Pattern Recognition, 2017.
  * [29] Dean Wilding, Oleg Soloviev, Paolo Pozzi, Gleb Vdovin, and Michel Verhaegen.  Blind multi-frame deconvolution by tangential iterative projections (TIP).  Opt. Express, 25(26):32305–32322, 2017.
  * [30] Liangpei Zhang, Qiangqiang Yuan, Huanfeng Shen, and Pingxiang Li.  Multiframe image super-resolution adapted with local spatial information.  J. Opt. Soc. Am. A, 28(3):381–390, 2011.
  * [31] N. Zhao, C. Li, H. Shi, and C. Lin.  Multi-frame image super-resolution based on regularization scheme.  In 2011 International Conference on Control, Automation and Systems Engineering (CASE), pages 1–4, 2011.
  * [32] Qinbang Zhou, Renwen Chen, Bin Huang, Chuan Liu, Jie Yu, and Xiaoqing Yu.  An automatic surface defect inspection system for automobiles using machine vision methods.  Sensors, 19(3), 2019.