A. Murat Tekalp. “Image and Video Restoration.”
2000 CRC Press LLC. <>.
ImageandVideoRestoration
A.MuratTekalp
UniversityofRochester
53.1Introduction
53.2Modeling
Intra-FrameObservationModel
•
MultispectralObserva-
tionModel
•
MultiframeObservationModel
•
Regularization
Models
53.3ModelParameterEstimation
BlurIdentification
•
EstimationofRegularizationParameters
•
EstimationoftheNoiseVariance
53.4Intra-FrameRestoration
BasicRegularizedRestorationMethods
•
RestorationofIm-
agesRecordedbyNonlinearSensors
•
RestorationofImages
DegradedbyRandomBlurs
•

useofparticulardegradationmodels.Ontheotherhand,superresolutionreferstoestimatingan
imageataresolutionhigherthanthatoftheimagingsensor.Imagesequencefiltering(restoration
andsuperresolution)becomesespeciallyimportantwhenstillimagesfromvideoaredesired.This
isbecausetheblurandnoisecanbecomeratherobjectionablewhenobservinga“freeze-frame”,
althoughtheymaynotbevisibletothehumaneyeattheusualframerates.Sincemanyvideosignals
encounteredinpracticeareinterlaced,weaddressthecasesofbothprogressiveandinterlacedvideo.
c

1999byCRCPressLLC
Theproblemofimagerestorationhassparkedwidespreadinterestinthesignalprocessingcommu-
nityoverthe past20or 30years. Becauseimage restorationis essentiallyanill-posed inverseproblem
whichisalsofrequentlyencounteredinvariousotherdisciplinessuchasgeophysics,astronomy,med-
ical imaging, and computer vision, the literature that is related to image restoration is abundant. A
concisediscussion of early results can befound inthebooks by Andrewsand Hunt[1] andGonzalez
and Woods [2]. More recent developments are summarized in the book by Katsaggelos [3], and re-
view papers byMeinel[4], Demoment [5], Sezan and Tekalp[6], andKaufmanand Tekalp[7]. Most
recently, printinghigh-quality still images from video sources has become an important application
for multi-frame restoration and superresolution methods. An in-depth coverage of video filtering
methods can be found in the book D igital Video Processing by Tekalp [8]. This chapter summarizes
key results in digital image and video restoration.
53.2 Modeling
Every image restoration/superresolution algorithm is based on an observation model, which relates
the observed degraded image(s) to the desired “ideal” image, and possibly a regularization model,
whichconveystheavailablea priori informationabouttheideal image. Thesuccessofimage restora-
tion and/or superresolution depends on how good the assumed mathematical models fit the actual
application.
53.2.1 Intra-Frame Observation Model
Letthe observed and ideal imagesbe sampled on the same 2-D lattice . Then, the observedblurred
and noisy image can be modeled as
g = s(Df ) + v

,m
2
)
f
(
n
1
− m
1
,n
2
− m
2
)


+ v
(
n
1
,n
2
)
(53.2)
where d(m
1
,m
2
) and S
d

1
,n
2
)
d
(
n
1
,n
2
; m
1
,m
2
)
f
(
m
1
,m
2
)


+ v
(
n
1
,n
2




g
1
.
.
.
g
K



, f
.
=



f
1
.
.
.
f
K



, v

.
.
.
.
.
.
.
.
D
K1
··· D
KK



is an N
2
K × N
2
K matrix representing the multispectral blur operator. In most applications, D is
block diagonal, indicating no inter-band blurring.
53.2.3 Multiframe Observation Model
Supposeasequenceofblurredandnoisyimagesg
k
(n
1
,n
2
),k = 1, ,L, correspondingtomultiple
shots (from different angles) of a static scene sampled on a 2-D lattice or frames (fields) of video

)
∈
S
d
(
n
1
,n
2
;k
)
d
k
(
n
1
,n
2
; m
1
,m
2
)
f
(
m
1
,m
2
)

2
) withrespecttothe hig h resolution
grid, and whether there is additional optical (out-of-focus, motion, etc.) blur. Because the relative
positionsoflow-andhigh-resolution pixelsingeneralvary byspatialcoordinates,the discretesensor
PSFisspace-vary ing. Thesupportofthespace-varyingPSFisindicatedbytheshadedareainFig.53.1,
wheretherectangle depictedbysolidlinesshowsthesupport of a low-resolutionpixeloverthe high-
resolutionsensorarray. The shaded region corresponds to the area swept bythelow-resolutionpixel
due to motion duringthe aperture time [8].
c

1999 by CRC Press LLC
FIGURE 53.1: Illustration of the discrete system PSF.
Note that the model (53.5) is invalid in case of occlusion. That is, each observed pixel (n
1
,n
2
; k)
canbe expressedas alinearcombinationofseveraldesiredhig h-resolutionpixels(m
1
,m
2
),provided
that (n
1
,n
2
; k) is connected to (m
1
,m
2

,m
2
)
∈
S
c
c
(
m
1
,m
2
)
f
(
n
1
− m
1
,n
2
− m
2
)
+ w
(
n
1
,n
2

Blur identification refers to estimation of both the support and parameters of the PSF {d(n
1
,n
2
) :
(n
1
,n
2
) ∈ S
d
}. It is a crucialelement of image restoration because the quality of restored images is
highly sensitiveto errors in thePSF [14]. An early approach to blur identification has been based on
the assumption that the original scene contains an ideal point source, and that its spread (hence the
PSF) can be determined from the observed image. Rosenfeld and Kak [15] show that the PSF can
alsobedeterminedfromanidealline source. These approachesare oflimiteduse in practicebecause
a scene, in general, doesnot contain an ideal point or line source and the observation noise may not
allow the measurement of a useful spread.
Models for certain types of PSF can be derived using principles of optics, if the source of the
blur is known [7]. For example, out-of-focus and motion blur PSF can be parameterized with afew
parameters. Further,theyarecompletelycharacterizedbytheirzerosinthefrequency-domain. Power
spectrumandcepstrum(Fouriertransformofthelogarithmofthepowerspectrum)analysismethods
have been successfully applied in many cases to identify the location of these zero-crossings [ 16, 17].
Alternatively, Chang et al. [18] proposed a bispectrum analysis method, which is motivated by the
fact that bispectrum is not affected, in principle, by the observation noise. However, the bispectral
method requires much more data than the method based on the power spectrum. Note that PSFs,
whichdonothavezerocrossingsinthe frequencydomain(e.g., GaussianPSFmodeling atmospheric
turbulence), cannot be identified by these techniques.
Yetanotherapproachforbluridentificationisthemaximumlikelihood(ML)estimationapproach.
TheMLapproachaimstofind those parameter values (including, in pr inciple, theobservationnoise

important role in defining constraints used in some of the restoration algorithms.
53.4 Intra-Frame Restoration
Westartbyfirstlookingatsomebasicregularizedrestorationstrategies,inthecaseofanLSIblurmodel
withnopointwisenonlinearity. Theeffectofthenonlinearmappings(.)isdiscussedinSection53.4.2.
Methods that allow PSFs with a random components are summarized in Section 53.4.3. Adaptive
restoration for ringing suppression and blind restoration are covered in Sections 53.4.4 and 53.4.5,
respectively. Restoration of multispectral images and space-varyingblurred images are addressed in
Sections 53.4.6 and 53.4.7, respectively.
53.4.1 Basic Regularized Restoration Methods
When the mapping s(.) is ignored, it is evident from Eq. (53.1) that image restoration reduces to
solving a set of simultaneous linear equations. If the matrix D is nonsingular (i.e., D
−1
exists) and
the vector g lies in the column space of D (i.e., there is no observation noise), then there exists a
uniquesolutionwhichcanbefoundbydirect inversion(alsoknown as inversefiltering). In practice,
however,wealmostalwayshaveanunderdetermined(duetoboundarytruncationproblem[14])and
inconsistent(due to observation noise) setof equations. In this case,we resort to a minimum-norm
least-squaressolution. A least squares (LS) solution (notunique when the columns of D arelinearly
dependent) minimizes the norm-square of the residual
J
LS
(f )
.
=||g − Df ||
2
(53.7)
LS solution(s) with the minimum norm (energy) is (are) generally known as pseudo-inverse solu-
tion(s) (PIS).
Restorationbypseudo-inversionis oftenill-posed owingtothepresenceofobservationnoise[14].
This follows because the pseudo-inverse operator usually has some very large eigenvalues. For ex-

i
and u
i
are the eigenvectors of D
T
D and DD
T
, respectively,
andR isthe rankofD. Clearly,reciprocation ofzerosingular-valuesisavoidedsince thesummation
runs to R, the rank of D. Under the assumption that D is block-circulant (corresponding to a
circular convolution), the PIS computed through Eq. (53.8) is equivalent to the frequency domain
c

1999 by CRC Press LLC
pseudo-inverse filtering
D
+
(u, v) =

1/D(u, v) if D(u, v) = 0
0 if D(u, v) = 0
(53.9)
where D(u, v) denotes the frequency response of the blur. This is because a block-circulant matrix
can be diagonalized by a 2-D discrete Fourier transformation (DFT) [2].
Regularization of the PIS can then be achieved by truncating the singular value expansion (53.8)
to eliminate all terms corresponding to small λ
i
(which are responsible for the noise amplification)
at the expense of reduced resolution. Truncation str ategies are generally ad-hoc in the absence of
additional information.

f
k
+ RD
T

g − Df
k


(53.11)
whereC is a nonexpansiveconstraint operator, i.e., ||C(f
1
) − C(f
2
)||≤||f
1
− f
2
||, to guarantee
theconvergenceoftheiterations. ApplicationofEq.(53.11)toimagerestorationhasbeenextensively
studied (see [31, 32] and the references therein).
Constrained Least Squares Method
Regularizedimagerestorationcanbeformulatedasaconstrainedoptimizationproblem,where
a functional ||Q(f )||
2
of the image is minimized subject to the constraint ||g − Df ||
2
= σ
2
.Here

D
H
g (53.13)
c

1999 by CRC Press LLC
where
H
stands for Hermitian (i.e., complex-conjugate and transpose). The parameter γ =
1
α
(the
regularization parameter) must be such that the constraint ||g − Df ||
2
= σ
2
is satisfied. It is often
computed iteratively [2]. A sufficient condition for the uniqueness of the CLS solution is that Q
−1
exists. For space-invariant blurs, the CLS solution canbe expressed in the frequency domain as [34]
ˆ
F (u, v) =
D
∗
(u, v)
|D(u, v)|
2
+ γ |L(u, v)|
2
G(u, v) (53.14)

stands for the
powerspectrum of the ideal image. The powerspectrumofthe ideal image is usually estimated from
a prototype. It can be easily seen that the CLS estimate (53.14) reduces to the Wiener estimate by
setting |L(u, v)|
2
= σ
2
v
/|P (u, v)|
2
and γ = 1.
A Kalman filter determines the causal (up to a fixed lag) LMMSE estimate recursively. It is based
on a state-space representation of the image and observ ation models. In the first step of Kalman
filtering, a prediction of the present state is formed using an autoregressive (AR) image model and
the previous state of the system. In the second step, the predictions are updated on the basis of the
observed image data to form the estimate of the present state. Woods and Ingle [39] applied 2-D
reduced-updateKalmanfilter (RUKF)toimagerestoration, wherethe updateislimited toonly those
state variables in a neighborhood of the present pixel. The main assumption here is that a pixel is
insignificantly correlated with pixels outside a certain neighborhood about itself. More recently, a
reduced-ordermodelKalmanfiltering(ROMKF),wherethestatevectoristruncatedtoasizethatison
the order of the image modelsupport has beenproposed [40]. Other Kalmanfiltering formulations,
including higher-dimensional state-space models to reduce the effective size of the state vector, have
been reviewed in [7]. The complexity of higher-dimensional state-space model based formulations,
however, limits their pr actical use.
Maximum A posteriori Probability Method
Themaximumaposterioriprobability(MAP)restorationmaximizestheaposterioriprobability
density function (pdf) p(f |g), i.e., the likelihood of a realization of f being the ideal image given
the observed data g. Through the application of the Bayes rule, we have
p(f |g) ∝ p(g|f )p(f )
(53.16)

(53.17)
whereR
v
denotes the covariancematrix of the noise process. Unlike the LMMSEmethod, theMAP
method uses complete pdf information. However, if both the image and noise are assumed to be
homogeneous Gaussian random fields, the MAP estimate reduces to the LMMSE estimate, under a
linear observation model.
Trusselland Hunt [10] used non-stationarya prioripdf models, andproposeda modifiedform of
thePicarditerationtosolvethe nonlinear maximizationproblem. They suggestedusingthe variance
of the residualsignal as a criterionforconvergence. Geman and Geman[11] proposedusing a Gibbs
randomfield modelforthea prioripdfoftheidealimage. Theyusedsimulatedannealingprocedures
to maximize Eq. (53.16). It should be noted that the MAP procedures usually require significantly
more computation compared to, for example, the CLS or Wiener solutions.
Maximum Entropy Method
A number of maximum entropy(ME) approacheshavebeen discussed in theliterature,which
vary in the way that the ME principle is implemented. A common feature of all these approaches,
however, is their computational complexity. Maximizing the entropy enforces smoothness of the
restored image. (In the absence of constraints, the entropy is highest for a constant-valued image).
One importantaspect of the ME approach is that the nonnegativity constraint isimplicitly imposed
on the solution because the entropy is defined in terms of the logarithm of the intensity.
Frieden was the first to apply the ME principle to image restoration [41]. In his formulation, the
sum of the entropy of the image and noise, given by
J
ME1
(f ) =−

i
f(i)ln f(i)−

i


i
f(i)= K
.
=

i
g(i) (53.23)
on the restored image. The optimization problem is solved using an ascent algorithm. Trussell [43]
showed that in the case of a prior distribution defined in terms of the image entropy, the MAP
solution is identical to the solution obtained by this ME formulation. Other ME formulations were
also proposed [44, 45]. Note that all ME methods are nonlinear in nature.
c

1999 by CRC Press LLC
Set-Theoretic Methods
Inset-theoreticmethods,firstanumberof“constraintsets”aredefinedsuchthattheirmembers
are consistent withthe observations and/or some a priori information about the ideal image. A set-
theoretic estimate of the ideal image is then defined as a feasible solution satisfying all constraints,
i.e., any member of the intersection of the constraint sets. Note that set-theoretic methods are, in
general, nonlinear.
Set-theoretic methods vary according to the mathematical properties of the constraint sets. In
the method of projections onto convex sets (POCS),the constraint sets C
i
are closed and convex in
an appropriate Hilbert space H. Given the sets C
i
, i = 1, ,M, and their respective projection
operators P
i

Image sensors and media may have nonlinear characteristics that can be modeled by a pointwise
(memoryless) nonlinearity s(.). Common examples are photographic film and paper, where the
nonlinear relationship between the exposure (intensity)and the silver densitydeposited on the film
orpaperis specifiedbya“d − loge”curve. The modelingofsensornonlinearitieswasfirstaddressed
by Andrews and Hunt [1]. However, it was not generally recognized that results obtained by taking
the sensor nonlinearity into account may be far more superior to those obtained by ignoring the
sensor nonlinearity, until the experimental work of Tekalp and Pavlovi
´
c[54, 55].
ExceptfortheMAPapproach,noneofthealgorithmsdiscussedaboveareequippedtohandlesensor
nonlinearity in a straightforward fashion. A simple approach would be to expand the observation
modelwiths(.) intoitsTaylorseriesaboutthemeanoftheobservedimageandobtainanapproximate
(linearized)model,whichcanbeusedwithanyoftheabovemethods[1]. However, theresultsdonot
show significant improvement over those obtained by ignoring the nonlinearity. The MAP method
is capable of taking the sensor nonlinearity into account directly. A modified Picard iteration was
proposedin[10], assuming both the image and noise are Gaussian distributed, which is given by
ˆ
f
k+1
=
¯
f
k
+ R
f
D
T
S
b
R

the ideal and blurred images in the exposure domain. However, the additive noise in the density
domain manifests itself as multiplicative noise in the exposure domain. To this effect, Tekalp and
Pavlovi
´
c[54] derive an LMMSE deconvolution filter in the presence of multiplicative noise under
certainassumptions. Theirresultsshowthataccountingforthesensornonlinearitymaydramatically
improve restoration results [54, 55].
53.4.3 Restoration of Images Degraded by Random Blurs
Basic regularized restoration methods (reviewed in Section 53.4.1) assume that the blur PSF is a
deterministic function. A more realistic model may be
D =
¯
D + D
(53.26)
where
¯
D is the deterministic part(known or estimated) of the blur operator andD stands for the
random component. Random component may represent inherent random fluctuations in the PSF,
for instance due to atmospheric turbulence or random relative motion, or it may model the PSF
estimation error.
A naive approach would be to employ the expected value of the blur operator in one of the
restoration algorithms discussed above. The resulting restoration, however, may be unsatisfactor y.
Slepian [56] derived the LMMSE estimate, which explicitly incorporated the randomcomponent of
the PSF. The resulting Wiener filter requires the a priori knowledge of the second order statistics of
theblurprocess. Wardetal. [57, 58] also proposedLMMSEestimators. Combettesand Trussell[59]
addressedrestorationofrandomblurswithintheframeworkofPOCS,wherefluctuations inthePSF
are reflected in the bounds defining the residual constraint sets. The method of total least squares
(TLS) has beenusedin the mathematicsliteraturetosolvea set oflinearequations withuncertainties
in the system matr ix. The TLSmethod amounts to findingthe minimum perturbations on D and g
to make the system of equations consistent. A variation of this principle has been applied to image

blurred images [63] and image recovery from Fourier phase information [64]. Lagendijk et al. [63]
appliedtheE-M algorithmtoblindimage restoration,whichalternatesbetweenMLparameteriden-
tification and minimum mean square error image restoration. Chen et al. [64] employed the POCS
method to estimatethe Fourier magnitude of the ideal image fromtheFourierphase of the observed
blurred image by assuming a zero-phase blur PSF so that the Fourier phase of the observed image is
undistorted. Both methods require the PSF to be real and symmetric.
53.4.6 Restoration of Multispectral Images
Atrivial solutionto multispectral image restoration, when there is no inter-bandblurring, may be to
ignore the spectral correlations among different bands and restore each band independently, using
oneofthealgor ithmsdiscussedabove. However,algorithmsthatareoptimalforsingle-bandimagery
may no longer be so when applied to individual spectral bands. For example, restoration of the red,
green, and blue bands of a color image independently usually results in objectionable color shift
artifacts.
To this effect, Hunt and Kubler [65] proposed employing the Karhunen-Loeve (KL) transfor m
to decorrelate the spectral bands so that an independent-band processing approach can be applied.
However, because the KL transform is image dependent, they then recommended using the NTSC
YIQ transformation as a suboptimum but easy-to-use alternative. Experimental evidence shows
that the visual quality of restorations obtained in the KL, YIQ, or another luminance-chrominance
domain are quite similar [65]. In fact, restoration of only the luminance channel suffices in most
cases. Thismethodappliesonlywhenthereisnointer-bandblurring. Further,oneshouldrealizethat
the observation noise becomes correlated with the image under a non-orthogonal transformation.
Thus,filteringbasedon the assumptionthattheimage and noiseareuncorrelatedisnot theoretically
founded in the YIQ domain.
Recent efforts in multispectral image restoration are concentrated on making total use of the
inherent correlations between the bands [66, 67]. Applying the CLS filter expression (53.13)tothe
observation model (53.4) with Q
H
Q = R
−1
f


R
f ;11
··· R
f ;1K
.
.
.
.
.
.
.
.
.
R
f ;K1
··· R
f ;KK





, and R
v
.
=




T
j
} and R
v;ij
.
= E{v
i
v
T
j
}, i, j = 1, 2, ,K denote the inter-band, cross-
correlation matrices. Note that if R
f ;ij
= 0 for i = j, i, j = 1, 2, ,K, then the multiframe
c

1999 by CRC Press LLC
estimate becomes equivalent to stacking the K single-frame estimates obtained independently.
Directcomputationof
ˆ
f through Eq. (53.27) requiresinversionof a N
2
L × N
2
L matrix. Because
the blur PSF is not necessarily the same in each band and the inter-band correlations are not shift-
invariant,thematricesD, R
f
,andR
v

v
is diagonal with all
diagonalentriesequaltoσ
2
v
. EstimationofthemultispectralcorrelationmatrixR
f
canbeperformed
by either the periodogram method or 3-D AR modeling [68]. Sezan and Trussell [69] show that the
multispectral Wiener filter is highly sensitive to the cross-power spectral estimates, which contain
phase information. Other multispectral restoration methods include Kalman filtering approach of
TekalpandPavlovi
´
c[67],leastsquaresapproachesofOhyamaetal.[70]andGalatsanosetal.[71],and
set-theoreticapproachof SezanandTrussell[23, 69] whoproposedmultispectral image constraints.
53.4.7 Restoration of Space-Varying Blur red Images
Inprinciple,allbasicregularizationmethodsapplytotherestorationofspace-varyingblurredimages.
However, because Fourier transforms cannot be utilized to simplify large matrix operations (such
as inversion or singular value decomposition) when the blur is space-varying, implementation of
some ofthese algorithms may be computationally formidable. There exist three distinct approaches
to attack the space-variant restoration problem: (1) sectioning, (2) coordinate transformation, and
(3) direct approaches.
The main assumption in sectioning is that the blur is approximately space-invariant over small
regions. Therefore, a space-varying blurred image can be restored by applying the well-known
space-invariant techniques to local image regions. Trussell and Hunt [73] propose using iterative
MAP restoration within rectangular, overlapping regions. Later, Trussell and Fo gel proposed using
a modified Landweber iteration [21]. A major drawback of sectioning methods is generation of
artifacts at the region boundaries. Overlapping the contiguous regions somewhat reduces these
artifacts, but does not completely suppress them.
Most space-varying PSF vary continuously from pixel to pixel (e.g., relative motion with acceler-

y :|r
(y)
(n
1
,n
2
)|≤δ
0

(53.28)
and
r
(y)
(
n
1
,n
2
)
.
= g
(
n
1
,n
2
)
−

(

2
)
(53.29)
is the residual at pixel (n
1
,n
2
) associated with y, which denotes an arbitrary member of the set.
The quantityδ
0
is ana priori bound reflecting the statistical confidence with which the actual image
is a member of the set C
n
1
,n
2
. Since r
(f )
(n
1
,n
2
) = v(n
1
,n
2
), the bound δ
0
is determined from
the statistics of the noise process so that the ideal image is a member of the set within a certain

=











x
(
i
1
,i
2
)
+
r
(x)
(
n
1
,n
2
)
−δ
0

(x)
(
n
1
,n
2
)
>δ
0
x
(
i
1
,i
2
)
if − δ
0
≤ r
(x)
(
n
1
,n
2
)
≤ δ
0
x
(

1
,o
2
)
h
(
n
1
,n
2
; i
1
,i
2
)
if r
(x)
(
n
1
,n
2
)
< −δ
0
(53.30)
The algorithm starts with an arbitrary x(i
1
,i
2

F
k
(u, v) =
S
f ;k
(u, v)
N

i=1
S
∗
f ;i
(u, v)D
∗
i
(u, v)G
i
(u, v)
N

i=1
|S
f ;i
(u, v)D
i
(u, v)|
2
+ σ
2
v

methods.)
Severalearlymotion-compensatedmethodsarein theformoftwo-stageinterpolation-restoration
algorithms[79,80]. Theyarebasedonthepremisethatpixelsfromallobservedframescanbemapped
backontoa desiredframe, based on estimatedmotiontrajectories, toobtainanupsampledreference
frame. However, unless we assume global translational motion, the upsampled reference frame is
nonuniformly sampled. Inordertoobtaina uniformly spaced upsampled image,interpolationonto
a uniform sampling grid needs to be performed. Image restoration is subsequently applied to the
upsampled image to remove the effect of the sensor blur. However, these methods do not use an
accurate image formation model, and cannot remove aliasing ar tifacts.
Motion-compensated (multiframe) superresolution methods that are based on the model (53.5)
canbeclassifiedasthosethataimtoeliminate(1)aliasingonly,(2)aliasingandLSIblurs,and(3)alias-
ingand space-varyingblurs. Inaddition, someofthese methodsaredesignedforglobaltr anslational
motiononly,while otherscanhandle space-varyingmotion fieldswithocclusion. Multiframesuper-
resolution was first introduced by Tsai and Huang [81] who exploited the relationship between the
continuous and discrete Fourier transforms of the undersampled frames to remove aliasing errors,
in the special case ofglobal motion. Their formulation has been extendedby Kim et. al. [82]to take
into account noise and blur in the low-resolution images, byposing theproblem in the least squares
sense. A further refinement by Kim and Su [83] allowed blurs that are different for each frame of
c

1999 by CRC Press LLC
low-resolution data, by using a Tikhonov regularization. However, the resulting algorithm did not
treat the formationof blur due to motion or sensor size, and suffers from convergence problems.
Inspection of the model (53.5) suggests that the superresolution problem can be stated in the
spatio-temporal domain as the solution of a set of simultaneous linear equations. Suppose that
the desired high-resolution frames are M × M, and we have L low-resolution observations, each
N × N . Then, from Eq. (53.5), we can set up at most L × N × N equations in M
2
unknowns to
reconstruct a particular hig h-resolution frame. These equations are linearly independent provided

intraframe restoration of space-varying blurred images. In this case, we define a different closed,
convex set for each observed low-resolutionpixel (n
1
,n
2
,k)(which can be connected to the desired
frame i by a motion trajectory)as
C
n
1
,n
2
;i,k
=

x
i
(
m
1
,m
2
)
:|r
(x
i
)
k
(
n

,n
2
)
−
M−1

m
1
=0
M−1

m
2
=0
x
i
(
m
1
,m
2
)
h
ik
(
m
1
,m
2
; n

(m
1
,m
2
) willbe. Ingeneral, itis desirable thatL>M
2
/N
2
.
Note,however,thatthePOCSmethodgeneratesareconstructedimagewithanynumberL ofavailable
frames. Thenumber L isjust anindicatorofhowlargethefeasibleset ofsolutionswillbe. Ofcourse,
the size of the feasible set can be further reduced by employing other closed, convex constraints in
the form of statistical or structural image models.
c

1999 by CRC Press LLC
53.6 Conclusion
Atpresent,factorsthatlimitthesuccessofdigitalimagerestorationtechnologyincludelackofreliable
(1) methods for blur identification, especially identification of space-variant blurs, (2) methods to
identify imaging system nonlinearities, and (3) methods to deal with the presence of artifacts in
restored images. Our experience with the restoration of real-life blurred images indicates that the
choice of a par ticular regularization strategy (filter) has a small effect on the qualit y of the restored
images as long as the parameters of the degradation model, i.e., the blur PSF and the SNR, and any
imagingsystemnonlinearityisproperlycompensated. Propercompensationofsystemnonlinearities
also plays a sig nificant role in blur identification.
References
[1] Andrews, H.C. and Hunt, B.R., Digital Image Restoration, Prentice-Hall, Englewood Cliffs,
NJ, 1977.
[2] Gonzales, R.C. and Woods, R.E.,
Digital Image Processing, Addison-Wesley, MA, 1992.

IEEE Trans. Sign. Proc., SP-39(3), 683–697, 1991.
[14] Tekalp, A.M. and Sezan, M.I., Quantitative analysis of artifacts in linearspace-invariant image
restoration,
Multidim. Syst. and Signal Proc., 1(1), 143–177, 1990.
[15] Rosenfeld, A. and Kak, A.C.,
Dig ital Picture Processing, Academic, New York, 1982.
[16] Gennery, D.B., Determination ofoptical transfer function by inspection of frequency-domain
plot,
J. Opt. Soc. Am., 63(12), 1571–1577, 1973.
[17] Cannon, M., Blind deconvolution of spatially invariant image blurs with phase,
IEEE Trans.
Acoust. Speech Sig. Proc.,
ASSP-24(1), 58–63, 1976.
[18] Chang, M.M., Tekalp, A.M. and Erdem, A.T., Blur identification using the bispectr um,
IEEE
Trans. on Sign. Proc.,
ASSP-39(10), 2323–2325, 1991.
[19] Lagendijk, R.L., Tekalp, A.M. and Biemond, J., Maximum likelihood image and blur identifi-
cation: A unifying approach,
Opt. Eng., 29(5), 422–435, 1990.
[20] Pavlovi
´
c, G. and Tekalp, A.M., Maximum likelihood parametric blur identification based on a
continuous spatial domain model,
IEEE Trans. Image Proc., 1(4), 496–504, 1992.
c

1999 by CRC Press LLC
[21] Trussell, H.J. and Fogel, S., Identification and restoration of spatially variant motion blurs in
sequential images,

[31] Biemond, J., Lagendijk, R.L. and Mersereau, R.M., Iterative methods for image deblurring,
Proc. IEEE, 78(5), 856–883, 1990.
[32] Katsaggelos, A.K., Iterative image restoration algorithms,
Opt. Eng., 28(7), 735–748, 1989.
[33] Tikhonov, A.N. and Arsenin, V.Y.,
Solutions of Ill-Posed Problems, V. H. Winston and Sons,
Washington, D.C., 1977.
[34] Hunt, B.R., The application of constrained least squares estimation to image restoration by
digital computer,
IEEE Trans. Comput., C-22(9), 805–812, 1973.
[35] Miller, K., Least squares method for ill-posed problems with a prescribed bound,
SIAM J.
Math. Anal.,
1, 52–74, 1970.
[36] Lagendijk, R.L., Biemond, J. and Boekee, D.E., Regularized iterative image restoration with
ringing reduction,
IEEE Trans. Acoust. Speech Sig. Proc., 36(12), 1874–1888, 1988.
[37] Zhou, Y.T., Chellappa,R.,Vaid,A.andJenkins,B.K.,Imagerestorationusinganeuralnetwork,
IEEE Trans. Acoust. Speech Sig. Proc., ASSP-36(7), 1141-1151, 1988.
[38] Yeh,S.J.,StarkH.andSezan,M.I., Hopfield-type neural networks: theirset-theoreticformula-
tions as associative memories, classifiers, and their application to image restoration, in
Digital
Image Restoration,
Katsaggelos, A. Ed., Springer Verlag, Berlin, 1991.
[39] Woods, J.W. and Ingle, V.K., Kalman filtering in two-dimensions-further results,
IEEE Trans.
Acoust. Speech Sig. Proc.,
ASSP-29, 188–197, 1981.
[40] Angwin, D.Land Kaufman,H., Imagerestorationusingreducedordermodels,
Sig. Processing,

IEEE Trans. Acoust.
Speech Sig. Proc.,
ASSP-32(4), 201-212, 1984.
[50] Youla,D.C.,Generalizedimagerestorationbythemethodofalternatingorthogonalprojections,
IEEE Trans. Circuits Syst., CAS-25(9), 694–702, 1978.
[51] Youla, D.C. and Velasco, V., Extensions of a result on the synthesis of signalsin the presence of
inconsistent constraints,
IEEE Trans. Circuits Syst., CAS-33(4), 465–467, 1986.
[52] Stark, H., Ed.,
Image Recovery: Theory and Application, Academic, Florida, 1987.
[53] Civanlar,M.R.andTrussell,H.J.,Digitalimagerestorationusingfuzzysets,
IEEETrans.Acoust.
Speech Sign. Proc.,
ASSP-34(8), 919-936, 1986.
[54] Tekalp, A.M. and Pavlovi
´
c, G., Image restoration with multiplicative noise: Incorporating the
sensor nonlinearity,
IEEE Trans. Sign. Proc., SP-39, 2132–2136, 1991.
[55] Tekalp,A.M. and Pavlovi
´
c, G., Digital restoration ofimages scanned fromphotographic paper,
J. Electronic Imaging, 2, 19–27, 1993.
[56] Slepian, D., Linear least squares filtering ofdistorted images,
J. Opt. Soc. Am., 57(7), 918–922,
1967.
[57] Ward, R.K. and Saleh, B.E.A., Deblurringrandom blur,
IEEE Trans. Acoust. Speech Sig. Proc.,
ASSP-35(10), 1494–1498, 1987.
[58] Quan, L. and Ward, R.K., Restoration of randomly blurred images by the Wiener filter,

[67] Tekalp, A.M. and Pavlovi
´
c, G., Multichannel image modeling and Kalman filtering for multi-
spectral image restoration,
Signal Process., 19, 221-232, 1990.
c

1999 by CRC Press LLC
[68] Ozkan, M.K., Erdem, A.T., Sezan, M.I. and Tekalp, A.M., Efficient multiframeWiener restora-
tion of blurred and noisy image sequences,
IEEE Trans. Image Proc., 1(4), 453–476, 1992.
[69] Sezan, M.I. and Trussell, H.J., Use of
a priori knowledge in multispectral image restoration,
Proc. IEEE ICASSP’89, Glasgow, Scotland, 1429–1432, 1989.
[70] Ohyama,N., Yachida,M., Badique, E., Tsujiuchi,J.and Honda,T.,Least-squaresfilter forcolor
image restoration,
J. Opt. Soc. Am., 5, 19–24, 1988.
[71] Galatsanos, N.P., Katsaggelos, A.K., Chin, R.T. and Hiller y, A.D., Least squares restoration of
multichannel images,
IEEE Trans. Sign. Proc., SP-39(10), 2222–2236, 1991.
[72] Ozkan, M.K., Tekalp, A.M. and Sezan, M.I., POCS-basedrestoration of space-varying blurred
images,
IEEE Trans. Image Proc., 3(3), 450–454, 1994.
[73] Trussell, H.J. and Hunt, B.R., Image restoration of space-variant blurs by sectioned methods,
IEEE Trans. Acoust. Speech Sig. Proc., ASSP-26(6) 608–609, 1978.
[74] Robbins, G.M. andHuang, T.S.,Inversefilteringforlinear shift-variantimaging systems,
Proc.
IEEE,
60(7), 1972.
[75] Sawchuck, A.A., Space-variant image restoration by coordinate transformations,

CVGIP: Graphical Models
and Image Proc.,
53, 231–239, 1991.
[85] Irani, M. and Peleg, S., Motion analysis for image enhancement: Resolution, occlusion and
transparency,
J. Vis. Comm. Image Rep., 4, 324–335, 1993.
[86] Stark, H.andOskoui,P.,High-resolutionimagerecoveryfromimageplanearraysusingconvex
projections,
J. Opt. Soc. Am., A 6, 1715–1726, 1989.
[87] Patti, A., Sezan, M.I. and Tekalp, A.M., Superresolution video reconstruction with arbitrar y
sampling lattices and nonzero aperture time,
IEEE Trans. Image Process., 6(8), 1064–1076,
1997.
[88] Schultz, R.R. and Stevenson, R.L., Extractionof high-resolution framesfrom video sequences,
IEEE Trans. Image Process., 5(6), 996–1011, 1996.
c

1999 by CRC Press LLC