A linear shift-invariant preprocessing technique is described which requires no specific knowledge of the image parameters and which is sufficiently general to allow the effective radius of the composite imaging system to be minimized while constraining other system parameters to remain within specified limits

Mcgillem, C. D.

Riemer, T. E.

English

NASA Technical Reports Server

f
' 022573
N 7 3 2 2 6 3 4
CONSTRAINED OPTIMIZATION OF IMAGE RESTORATION FILTERS
T,E, Riemer and C.P. McGillem
Purdue University, W, Lafayette, Indiana
Abstract
A linear shift-invariant preprocessing technique is described which re-
quires no specific knowledge of the image parameters and which is sufficiently
general to allow the effective radius of the composite imaging system to be
minimized while constraining other system parameters to remain within speci-
fied limits.
c
https://ntrs.nasa.gov/search.jsp?R=19730013907 2020-03-23T02:19:58+00:00Z
Introduction. The finite size of the point-spread function of a practical
imaging system causes the output to be a two-dimensional spatial smearing or
blurring of the original image. Methods of preprocessing the image to reduce
this blurring have been widely studied but no fully satisfactory processor has
as yet been developed. The preprocessing method proposed here, although new
in formulation, can be thought of as a generalization and extension of previous
work by Smith [Ij.] and Stuller [5]. This method accomodates the problem of non-
circular ly symmetric imaging system point-spread functions, provides for con-
trolled extent of the preprocessing filter to minimize distortion due to tran-
sients resulting from truncation errors and edge effects, can be used with
various kinds of system noise, and can be readily extended to provide constraint
of other system parameters. The analysis presented here relates to a line-
scanner system although it is applicable in principle to many other system
configurat ions.
Analysis of Preprocessing Technique. Referring to Fig. 1 and noting that a bar
over a variable denotes a two-dimensional spatial vector, a basis for choosing
the preprocessing filter point-spread function, h (v), in order to increase the
resolution of the composite imaging system is to minimize the radius of gyration,
R , of the composite imaging system point-spread function, g(v); given by
r" - 2|v| g(v) dv
R/
J° dv
where the composite system point-spread function is expressed
,.»
g(v) = 1 hr(z) hb(v-z) dz . (2)
-co
The magnitude of R may be minimized by setting the denominator equal to a con-
_ o
stant and minimizing the numerator with respect to h ('). The function |v| in
(l) is a penalty weighting function which causes h (v) to decrease as v increases.
Additional flexibility in controlling the shape of g(v) is obtained by using
a more general penalty weighting function w(v) which increases-at an arbitrarily
controlled rate outside a specified range of v. With this penalty function the
numerator and denominator of (l) can be written as the functionals
' w(v) g2(v) dv .(3)
oo
Ii- A2^ <*v - K! •- GO
-oo
Two additional constraints are considered essential for most imaging sys-
tems. The first limits the size, and to some extent the shape, of the prepro-
cessing filter point-spread function, h (v), so that truncation errors are mini-
mized. This constraint is
I2 - J s(v) hr2(v) dv = K2 (5)
-00
where s(v) is a penalty weighting function similar in form to w(v). The second
additional constraint limits the total noise power in the preprocessed image
as follows
J JRnn(i'^ \^ hr^> <^ dv = K (6)
-t - *o oo
where it is assumed that the noise source is a stationary ergodic random process
with autocorrelation function R (z-v).
nnv '
The required point-spread function is found by minimizing, with respect to
h (•)» the augmented functional formed by adding to (3) the constraints 00, (5)
and(6) each of which has been multiplied by an appropriate Lagrange multiplier [1,2]
The required minimization is obtained by converting to the spatial frequency
domain, setting the gradient of the augmented functional equal to zero and solving
the resulting equation subject to the constraints of (U)-(6). The gradient of
the augmented functional is as follows
-|H, (f ,f )j + 2\,L (f ,f ) H (f ,f ) = 0I1 DV x y'1 ^^ n^  x y ' J r ^ x y '
where H,(')» W(- ), S('), H (• ) and J(* ) are the two-dimensional Fourier trans-
forms of hb(-), w('), s('). hr(-) and Rnn(') respectively and X^ \2 and \ are
Lagrange multipliers associated with constraint equations (U)-(6) respectively.
The frequency domain representation of the constraint equations is
— |H (f f }| IH (f f )l df df = K
"1 JWJW bv x' y ' < rv x' y" x y~ 1
00 CO CO CO
_eo oo
Jo " I" 1" ? (f »f )IH (f »f )| df df * K. • (10)3 J J Innv x' y/! r x' y'1 x y 3 v '
—00—00
In order to proceed further an appropriate form for the penalty functions
w(v) and s(v) must be chosen. Several criteria must be satisfied by these
functions: they must increase rapidly away from the origin; they must be Fourier
transformable; and w(v) must be convex to assure the existence of a global
minimum of (3). A convenient representation centered at the origin is
2k 2k
w(v) = w (x)w (y) HC-5-) wx+ c KU-2-) wy+ c V .
 f...v v I v uv I I ^v ' u\r i 111 ix y L xw wxj (. yw wyj \,±±)
The constants 0 < c ,c < 1 are introduced to assure that the overall weighting
wx wy
function will be everywhere greater than zero thus avoiding singular behavior
of the solution. The constants k and k are positive integers. A similar
wx wy
expression may be written for s(v) = s (x)s (y). The shape of w(v) along one
x y
dimension is shown in Fig. 2. For the following analysis it is assumed that
k , k , k , and k are all unity. This assumption does not limit the general
wx wy sx sy
applicability of (y)-(lO) but is made to illustrate the method of solution of
these equations.
With the above assumption (7) can be put into the following form [3]
(f ,f )., 1). 2
X
 ^ 4—2- 5 H, (f ,f )H (f ,f ) -bit y c -*— H. (f ,f )H (f ,f )2^2 ) 2 2 b x x ' y / r v x y 7 3 w w . 2 b x x y ' r v x ' y '
^2
y
W
2c -^r H, (f ,f )H (f ,f ) + 87/Sc 2y 2c c H, (f ,f )H (f ,f )!•
w wx
 f 2 bv x' y' rv x' y' w 7w wx wy bv x y' rv x y J
y
XP r >U 9 P *2
+ —. o o 1—^ ^ H (f , f ) - UTTy c -2-^  H - (f , f )
8T/Sc y ^af af y syaf r x y
s s x y x
2
- UTT2x 2c -S-s H (f ,f ) + 8/x 2y 2c c H (f ,f )}
s sx..f 2 rv x yy s 7s sx sy rx x y'J
y
4[2X 1 lH b(f x , f y) | 2 +2X 3J n n(f x , f y)]H r(f x , f y) = 0. (12)
Equation (l2) is a two-dimensional fourth-order partial differential equation
which, when solved simultaneously with (8), (9), and (lO), specifies the spatial
frequency spectrum of the required preprocessing filter.
Image Collecting Systems with Separable Point-Spread Functions. For many
physically realizable image collecting systems, certain assumptions may be made
which greatly simplify the form of (12). A large class of physically realizable
image collecting system point-spread functions may be modelled as separable
point-spread functions, from which it follows that
H. ( £ , £ ) * H. (f ) H. (£) (13)b ^ x ' y ' b x v x / b y y / X J /
H (f ,£ ) - H (f ) H (f ) (110
r x y rx x ryx y7 '
(f ,f ) = J (f ) T (f ) . (15)x
 x y' -3-nnxx x -Lnnyx y' * •"nn  nx x -nn  y
With these assumptions, it can be shown that the solution for the two-
dimensional point-spread function of h (• ) defined by (12) and constraint equa-
tions (8)-(lO) reduces to two sets of equations [3]: one set specifying the
x-axis component of h (•)> the other the y-axis component. Since these equa-
tion sets are identical except for appropriate changes in x- and y-axis component
variables, only the equation set for the x-axis components will be given.
This equation set is composed of the second-order ordinary differential equation
2H* (f ) H' (f )x 2
HrVfx> + 22 S 2 *i*<fx>x X
 IIL (f )] x 2 2 rx x1
 bxv x'1 os 2x w
*-H (f ) H/'(f ) - Ulc (c +\, )|H. (f ) ! + ,2 | T?xN x' DV x' w L wx lx'1 t?xv x^' 3xrnn x sxJV „ /c \x < -- *" - - — 5 - 5 — - > 'VvA* / = 08 l 22 2
 J rx
 xIH . (£ ) | x - f X 0 x1
 DX X/( S 2X W
and the constraint equations
KxMXxM2 dfx
—3—5 iTn" (f ) H (f ) + H". (f ) H . (f ) I df
..2. 2 J L rrxv
 x rrx x' rixx xy rixx x' J xH« X —oo
+
 «„ f KxM2 «, ' S. (3S)
f?nn(O |H (f )|2 df = K (19)
V M.LL1L A. L A. X OA
•CO
where
H (f ) = H (f ) + J H . (f ) . (20)
rxv x rrx x' J rix x7 v '
Thus the spatial frequency transform of the x-axis component of h (•) is
obtained by the simultaneous solution of (l6) and the constraints of (17)-(19).
This system of equations can be solved by standard numerical procedures and
the functions h (•) and h (•) determined. The preprocessing operation thus
becomes two one-dimensional convolutions carried out sequentially rather than
a single two-dimensional convolution as would generally be the case.
Conclusions. The proposed method for specifying an image preprocessing point-
spread function is sufficiently general to allow a variety of penalty weighting
functions to be used for control of the effective radius of the point-spread
function while constraining other system parameters to remain within specified
limits.
Although (12) in conjunction with (8)-(lO) specify the general form of the
spatial frequency transform of the optimum preprocessing filter, it was shown
that in the case of a separable point-spread function the solution can be con-
siderably simplified. This same simplification also results for the case of
a radially symmetric point-spread function, however, the resulting processor
must perform a two-dimensional convolution in this case rather than two sequential
one-dimensional convolutions.
The simplified solution based upon separable point-spread functions can
be extended to a large class of practical systems by approximation techniques.
For example, by properly selecting the parameters of a two-dimensional Gaussian
function, an elliptical or circular point-spread function can be approximated
and preprocessing carried out as sequential one-dimensional convolutions; thus
substantially reducing the processing time.
\
Acknowledgement; This work was carried out at the Laboratory for Applications
of Remote Sensing, Purdue University under NASA grant NGL-15-005-112.
References:
1, Epstein, B., Linear Functional Analysis. W. B. Saunders, Co., 1970.
2, Franks, L. E., Signal Theory, Prentice-Hall, 1969.
3, McGillem, C.D., and T.E. Riemer, "A Linear Shift-Invariant Image Preprocessing
Technique for Multispectral Scanner Systems," LARS Information Note: 022673,
Laboratory for the Application of Remote Sensing, Purdue University.
U. Smith, H. A., "Improvement of the Resolution of a Linear Scanning Device,"
SIAM Journal on Applied Mathematics, January 1966, Vol. lU, No. 1.
5. Stuller, J. A., "An Algebraic Approach to Image Restoration Filter Design"
Computer Graphics and Image Processing, Dec. 1972, pp. 107-122.
6. Zemanian, A. H., Distribution Theory and Transform Analysis. McGraw-Hill,
1965-
i(V)
J^
Original
Image
hb(V)
Imaging
System
n(
\
•v/i_
V)
Noise
\? ,
-Ai/ x
+
hr(V)
Preprocessing
Filter
?(V)+nT(V£
^
Processed
Image
Figure 1. Block Diagram of Imaging and Preprocessing Systems.
l+c
k=2
vl V2 ^"V
Figure 2. One-Dimensional Penalty Weighting Function.
022573 Constrained Optimization of Image Restoration Filters.
TTE. Riemer and C.U^McGillenuA linear shift-
invariant preprocessing technique is described which
requires no specific knowledge of the image parameters
and which is sufficiently general to allow the effec-
tive radius of the composite imaging system to be
minimized while constraining other system parameters
to remain within specified limits.


Constrained optimization of image restoration filters

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server