 |
Description  |
|
|
BACKGROUND
The present invention relates generally to image processing, and more
particularly, to automatic subarea selection methods for use in
registering images.
In the prior an, selection of image subareas is based on scalar measures of
"edge density" or "entropy", such as is described in a paper entitled
"Spatial Registration of Multispectral Video Imagery" by Anuta, published
in SPIE J., Vol. 7, pp. 168-175, September 1969. In conventional practice,
there is no systematic way of considering subarea position and its impact
on registration accuracy. No general method is provided to limit the
uncertainty of registration parameters. No general method is provided for
handling complementary measurements because measurement quality indicators
are based on scalar quantities.
Accordingly, it is an objective of the present invention to provide methods
for use in registering images, and in particular, automatic subarea
selection methods that improve image registration between related images.
SUMMARY OF THE INVENTION
The present invention is a computer implemented method for automatically
selecting subareas from a reference image, such as a SAR image, for
example, such that registration accuracy is optimized between a current
image and the reference image. Optimal subarea selection reduces on-line
computation and reference data storage required for registration
effectiveness. The present method minimizes a cost function, corresponding
to the predicted total mean squared registration error (MSE), between the
images. The total MSE is predicted in terms of position and predicted
measurement covariance of each candidate subarea. The measurement
covariance for each subarea is predicted from local image statistics.
Combinatorial optimization procedures select a predetermined number of
subareas to minimize total MSE. The use of error prediction techniques and
the associated method for synthesizing optimal sets of subareas have not
been used in conventional image registration methods.
The selection criteria includes a vector quality (accuracy) metric for each
subarea. It is important to include a vector measure of subarea quality, a
covariance matrix rather than a scalar quality measure, to appropriately
handle measurements that are of good quality, say, in x but not in y. The
covariance matrix defines an error ellipse associated with a given
measurement. Previous subarea selection approaches have used scalar
indications of local measurement quality. The Anuta method referenced in
the Background section uses a measure of edge density to select subareas
for image registration.
The selection criteria includes subarea position (jointly with subarea
accuracy). For any transformation more general than pure translation,
scale factor (magnification) or rotation for example, registration
effectiveness depends on the relative positions of each measurement in
addition to the accuracy of the measurements themselves. For example,
moderately accurate measurements adequately dispersed in the image produce
better results than highly accurate measurements clustered in one comer of
the image. The present method jointly considers subarea layout and
covariance and therefore can distinguish between these cases. It can even
handle the situation where the addition of a near-degenerate but
well-positioned measurement is more effective than the addition of a
highly accurate but poorly-positioned measurement. The selection criteria
adapts to registration error model order. Also, the selection criteria
adapts to relative uncertainty in registration error parameters.
Experiments performed using the present invention indicate that the method
reduces on-line computation by a factor of 2 to 3 without loss of
accuracy. The present method considers the direction of measurement errors
in addition to their magnitude so that multiple locally-degenerate
measurements can be selected if they, together, yield a full-rank result
(complementary measurements). It explicitly considers the effect of
subarea (measurement) position and weighs this against the intrinsic
accuracy of the measurement. It systematically includes the prior
uncertainty into the selection process, preventing any parameter from
wrongly dominating the selection process. Also, it provides a general
method for balancing the competing requirements of different registration
parameter, for example, translation versus scale.
Thus, the present invention provides improvements over existing image
registration methods proposed to support change detection and other
multiple image analysis applications. The advantages provided by the
present invention are most significant for larger fields-of-view resulting
in higher order correction requirements. The present invention may be
adapted for use with military reconnaissance, attack, and counter-stealth
applications, and commercial remote sensing and mapping applications.
BRIEF DESCRIPTION OF THE DRAWINGS
The various features and advantages of the present invention may be more
readily understood with reference to the following detailed description
taken in conjunction with the accompanying drawings, wherein like
reference numerals designate like structural elements, and in which:
FIG. 1 shows an image processing system employing an image registration
method in accordance with the principles of the present invention;
FIG. 2 shows a flow diagram of the image registration method of the present
invention;
FIG. 3 shows a flow diagram illustrating the calculation of the total mean
square registration error performed in the image registration method of
FIG. 2;
FIG. 4 shows a flow diagram illustrating the calculation of the information
matrix performed in the image registration method of FIG. 2;
FIG. 5 shows an image with superimposed subareas, including rectangles that
show the extent of the subareas that is useful in illustrating the
advantages of the present invention; and
FIG. 6 shows an image that illustrates the importance of subarea layout and
of coupling subarea layout, vector quality measures, and an uncertainty
model, as is achieved by the present invention.
DETAILED DESCRIPTION
FIG. 1 shows an image processing system 10 employing an image registration
method 11 in accordance with the principles of the present invention. The
image registration method 11 operates off-line, and on reference (image)
data that is typically stored in memory, such as a hard disk, for example.
The reference data may be calibrated or uncalibrated image data. The image
registration method 11 is comprised of an uncertainty model 14 and a
subarea selection procedure 15. A more detailed description of the image
registration method 11 is provided below with reference to FIG. 2. The
uncertainty model 14 determines the relative uncertainty in locating the
actual position of an area of interest that is used in image registration
of live imagery. The subarea selection procedure 15 provides as its output
selected subareas 18 (position data) of the reference image that have been
determined to be optimal in registering the image. The uncertainty data is
combined with position data and is applied to an on-line processor 12 of
the system 10.
The on-line processor 12 operates on live imagery 21. The live imagery 21
is processed by a subarea extraction procedure 22 which uses a "window"
derived from the uncertainty data from the uncertainty model 14 and the
position data from the subarea selection procedure 15. The window (shown
as a dashed box around the selected subarea 18) comprises the image data
that has been extracted by the subarea extraction procedure 22. The
selected subareas 18 and the extracted subareas 26 are then correlated by
a multi-subarea correlation procedure 23. The output of the multi-subarea
correlation procedure 23 is applied to an image warping (resampling)
procedure 24 along with the live imagery 21 which registers a current
image to the reference image. Alternatively, and as is shown by the dashed
arrow identified as steering commands 26, the output of the multi-subarea
correlation procedure 23 may be directly used as inputs to a guidance
system, for example, that are used to guide a vehicle. The registered
current image is then output 25 to a display, printer or change detection
procedure, for example.
A detailed flow diagram of the image registration method 11 in accordance
with the principles of the present invention is shown in FIG. 2. The image
registration method 11 comprises an initialization procedure 30 and a
combinatorial optimization procedure 35.
By way of summary, in the initialization procedure 30, a covariance matrix
31 for each candidate subarea (i.e. candidate local measurement) is
predicted. Given this matrix 31 and the positions of candidate subareas 18
selected from the computed matrices by the subarea selection procedure 15,
the covariance of a distortion model parameter vector (A) is predicted to
produce the covariance matrix 31. The covariance of the distortion model
parameter vector A comprises the covariance of polynomial error
coefficients, wherein covariance propagation weights are derived from a
polynomial error model order and an image size parameter. For polynomial
distortion models, this covariance is independent of A and is therefore
computed from the reference image 13 alone without actually estimating A.
The final criterion for optimization is the predicted sum comprising the
calculation of a local mean squared registration error (MSE) 39 for each
pixel in the image. This is obtained by propagating the covariance of A
through to each pixel using the distortion model. Computation of the total
MSE objectively scales the different covariance terms according to their
impact on registration accuracy at the pixel level, avoiding any need for
ad hoc scaling of parameter covariance terms. A combinatorial optimization
procedure searches for the combination of subareas which results in the
lowest MSE calculation. The output of this optimization process is a list
of subarea positions.
To synthesize an optimal set of subareas, the impact on total MSE of each
candidate subarea 18 is calculated. The impact of eliminating a given
subarea generally depends on the other candidate subareas 18, and hence a
combinatorial optimization scheme is used to synthesize the best subset.
Two procedures are described below that comprise the image registration
method 11. The first starts with a maximal set of subareas and eliminates
those that minimize the increase in total MSE. The second starts with a
minimal set of subareas and then adds the subareas that maximize the
reduction of total MSE. The second method has been reduced to practice and
compared against selection based on minimizing the MSE of local
measurements alone and also a random selection of subareas.
The control structure for a subarea selection or extraction procedure 22
employed in the the image registration method 11 is the combinatorial
optimization procedure 35. Prior subarea selection methods have not used
combinatorial optimization. Data initialization 30 takes place prior to
the combinatorial optimization procedure 35.
Referring to FIG. 2, the first step in the subarea selection method 15 is
to calculate a local measurement quality (covariance) matrix 31. The input
of this process is the reference image 13, denoted by I.sub.1, the output
is a 2-by-2 matrix for each candidate subarea 18 indicating position
uncertainty (e.g. covariance matrix) or accuracy (e.g. the inverse
covariance, or information matrix). The inverse covariance C.sub.i.sup.-1,
referred to as the information matrix 31, is computed directly from local
image statistics, namely second moments of the spatial gradient of the
image. The computation is shown in detail in equations (7) through (10)
below.
The next step is to presort the subarea list 32 according to local quality.
This function takes all subareas for which the information matrix 31 has
been calculated and then sorts the list in decreasing order of information
(quality) or in increasing order of MSE (accuracy). This is a scalar
quantity for each subarea and does not take into account interaction with
other subareas. The output of this process is an ordered list of subarea
positions and information matrices 31. For example, the subareas may be
presorted in increasing order of local MSE given by the trace (sum of
diagonal elements) of the covariance matrix
tr(C.sub.i) (1)
The next step is to generate an initial subarea list 33. The input to this
process is the ordered list from the presorting step. This step selects
the best K subareas to provide the best initial list based on local scalar
criteria. Because the list 33 is sorted, the list 33 is comprised of the
first K entries. The output is a set of K pointers which specify the index
of each initial subarea in the full ordered list 33.
A control procedure 36 for the combinatorial optimization procedure 35
determines which subareas to add and/or delete from a current list in
order to form a new list, it compares the total MSE for the current
subarea list with that of the new list, and it determines whether the
stopping criterion has been met, i.e. whether to continue searching for a
better combination of subareas. Two realizations of this control procedure
37, namely steepest descent and steepest ascent are described below. Given
these two methods, the application of a more complicated combinatorial
optimization scheme, namely backtrack programming, is routine to those
skilled in the art.
The next step is to generate new list of candidate subareas 37. This step
creates a new list of subareas by adding or deleting K.sub.inc temporary
candidate subareas from the current list. It starts with the candidate
additions/deletions with the best local MSE and continues down the list
each time it is invoked.
The next step is to calculate total MSE 38 (optimization criterion). The
input to this step is the new list of candidate subareas 37. The output is
the total MSE, whose calculation is specified below.
The final step is to compare total MSE values 39. This step evaluates total
MSE for each new list of candidate subareas 19. This allows the control
process to decide which new list to accept.
The heart of the combinatorial optimization procedure 35 is the
optimization criterion, namely total MSE 38. Referring to FIG. 3, it shows
a flow diagram illustrating the calculation of the total mean square
registration error performed in the image registration procedure 35 of
FIG. 2. The inputs to the calculation of the total mean square
registration error procedure 38 are comprised of the list of subarea
positions 18, the local information matrix 31 (or covariance) associated
with each subarea 18, the prior covariance of the polynomial coefficients
17C, and image size 17B.
Polynomial error model order 17A. The preferred implementation is a set of
binary flags that indicate which polynomial coefficients are zero (or
known a priori) and which are nonzero (or unknown) would be estimated by
subsequent image registration processing. This model helps determine the
mapping from local MSE to total MSE 38.
Prior covariance of the polynomial coefficients 17C. The prior covariance
is a matrix B.sub.0 which indicates the uncertainty in polynomial error
model coefficients prior to application of the image registration
procedure. This is similar to the polynomial error model flags but it
includes the degree of uncertainty, instead of just a binary indicator.
For example, in a second order polynomial model, one might know from
operational limits on viewing geometry that the uncertainty in the
quadratic terms are much smaller than that of all lower order terms. One
can prevent the criterion function (equation (17) below) from being
inappropriately dominated by the effects of the quadratic terms by
specifying the relative degree of uncertainty using B.sub.0.
Calculation of covariance propagation weights 63. The inputs to this
process are the polynomial error model order 61 and the image size 62. The
calculations may be performed off-line given these inputs. The covariance
propagation weights 63, denoted by r.sub.jk, for the jkth polynomial term,
determine the mapping 65 from the covariance of polynomial coefficients 64
to total MSE 34, 38. The weights 63 are calculated by summing central
moments of image coordinates over the entire image. The combination of
moments used are determined from the polynomial model order 61. The
calculation method is detailed in equations (22) through (35) below.
Calculate covariance of polynomial coefficients 64. The inputs to this
procedure include the information matrix 31 for each candidate subarea 18
and all of the polynomial terms x.sub.i.sup.l y.sub.i.sup.m corresponding
to the position of the center (x.sub.i, y.sub.i) of each candidate subarea
18. The output B is the expected covariance matrix 64 for the polynomial
error coefficients. These polynomial coefficients are to be estimated by a
subsequent image registration procedure 12. The covariance matrix 64
indicates how accurate those coefficients are expected to be. Details of
the calculation are shown in equations (11) through (15) below. A modified
calculation which incorporates the a priori covariance B.sub.0 is
specified in equation (37) below.
Calculate total mean square registration error 65. The inputs to this
process are the covariance propagation weights r.sub.jk and the covariance
of polynomial coefficients B. This functions maps the covariance matrix 34
to a scalar quantity (total MSE) through a linear combination. The
rationale and method are shown in equations (16) through (21).
Optimization criterion calculation. Image registration requires estimation
of the coordinate transformation f that aligns two images. Assume that the
two images can be modeled as follows:
I.sub.1 (x)=I(x)+n.sub.1 (x) (2)
and
I.sub.2 (x)=I(f(x))+n.sub.2 (x), (3)
where the function I represents the repeatable, or noise-free, component of
the images and n.sub.1 and n.sub.2 are observation noise processes that
are independent of I. The total MSE criterion to be minimized is an
estimate of
##EQU1##
where .OMEGA. is the entire image, E denotes an ensemble average, and f
is the estimate of f obtained given the selected set of subareas 18.
For the formulation that follows, assume the coordinate transformation f
is estimated by first "measuring" its value x'.sub.i at several points
x.sub.i, in the image and then fitting a model to those measurements using
a weighted-least-squares method. This method minimizes the weighted sum of
squared residuals (or local measurement errors) given by
.epsilon..sub.i = f(x.sub.i)- x.sub.i. (5)
Finally, assume that x'.sub.i is estimated by maximizing the following
cross-correlation function
##EQU2##
where .OMEGA..sub.i is some window around x.sub.i.
Estimating the measurement information matrix. Then, the first step in
predicting MSE is to estimate the quality of each local measurement. Let
C.sub.i be the 2 by 2 covariance matrix for the ith measurement (subarea).
The inverse covariance C.sub.i.sup.-1, referred to as the information
matrix, is used and is computed directly from local image statistics. For
convenience, that development is summarized below.
An estimate of the measurement information matrix is given by
C.sub.i.sup.-1 =D.sub.i (2.sigma..sub.i.sup.2 D.sub.i +A.sub.w
b.sigma..sub.n.sup.4 I).sup.-1 D.sub.i, (7)
where .sigma..sub.n is an estimate of the standard deviation of the
additive noise (assumed equal for n.sub.1 and n.sub.2), D.sub.i is the
Hessian term with respect to x'.sub.i of the noise-free part of g.sub.i,
A.sub.w is the integral of the correlation window, I is the identity
matrix, and b<2. This estimation is easily accomplished during the process
of correlating I.sub.1 and I.sub.2 using a finite difference approximation
to the derivatives of the cross-correlation function. However, additional
assumptions are required in performing this estimation if only one image
I.sub.1 is observed, as is desired for off-line automatic subarea
selection.
Indirect calculation of the Hessian term. Given a single image, a
reasonable estimate of the Hessian term is given by
##EQU3##
where the second term compensates for an unknown bias due to the noise
covariance in I.sub.i.sup.x and I.sub.i.sup.y, w is the correlation
window, * denotes convolution in x=(x,y), and the superscript .sup.(x)
denotes the derivative with respect to x.
Assume that n.sub.1 is nearly white, and that the partial derivatives in
(8) are approximated using central differences as follows:
I.sub.i.sup.x (x,y)=1(I.sub.1 (x+1,y)-I.sub.1 (x-1, y)) (9)
and similarly for I.sub.i.sup.y. Then the noise terms from I.sub.i.sup.x
and I.sub.i.sup.y are approximately uncorrelated, yielding:
##EQU4##
The noise variance .sigma..sub.n.sup.2 is more difficult to estimate. For
single-look SAR imagery, for example, the a priori value .sigma..sub.n0
=5.57 dB may be used. For multi-look SAR imagery obtained by filtering and
downsampling higher resolution imagery, the image noise bandwidth should
be accounted for in calculating the variance of the multi-look pixels. For
other sensors, it is more desirable to estimate the noise bias from the
data using order statistics taken on the terms w * [I.sub.i.sup.x ].sup.2
and w * [I.sub.i.sup.y ].sup.2.
For spatially separable w, the convolutions may be performed separably,
requiring computation proportional to the linear width of w instead of the
area of w. For uniform windows, the convolutions may be performed
recursively with computation proportional to the image size. Depending on
the mechanization of the synthesis procedure and on system architecture,
the most efficient computational approach is to evaluate (10) only at the
centers of candidate subareas 18.
FIG. 4 shows a flow diagram illustrating the calculation of the information
matrix 31 performed in the image registration method 11 of FIG. 2. The
procedure for calculating the measurement information matrix 31 is as
follows. Calculate derivative images I.sub.i.sup.x and I.sub.i.sup.y, step
41. Calculate or look up b and A.sub.w. Estimate .sigma..sub.n, step 43.
For each Hessian term: (a) calculate product of derivatives, e.g.
[I.sub.i.sup.x ].sup.2, step 44, (b) filter the product with window w,
step 45, and (c) downsample by W/2, step 47, where W is the equivalent
width of w. For each point in downsampled Hessian images: (a) convert the
Hessian term to information matrix using equation (7), step 50.
Estimating parameter covariance. Let A denote the vector of registration
parameters that define a particular realization of the transformation f
and A denote its weighted-least-squares estimate given measurements
obtained from a set S of selected subareas. The covariance of A is known
given S and each C.sub.i.sup.-1 and may be predicted using the estimates
C.sub.i.sup.-1.
For compactness, let
B=covA, (11)
and define the following column vectors
X=col(x.sub.i, i.epsilon.S) (12)
where each x.sub.i denotes the centroid of the image subarea used in
obtaining the ith measurement, and
F=col(f(x.sub.i), i.epsilon.S), (13)
and let the block diagonal measurement covariance matrix be given by
C=diag (C.sub.i, i.epsilon.S). (14)
Then,
B=(.gradient. F.sup.T C.sub.i.sup.-1 F).sup.-1 =(.gradient. f(
x.sub.i).sup.T C.sub.i.sup.-1 .gradient. f(x.sub.i)).sup.-1(15)
where .gradient. is the gradient with respect to A and ()T denotes matrix
transposition.
Optimization criterion. The objective is to select the set S of subareas 18
such that the following optimization criterion is minimized:
##EQU5##
where .OMEGA. is a set of coordinates covering the entire image.
To relate J to covA,
##EQU6##
with equality holding for polynomial models.
Simplification for polynomial models. If f can be modeled as a polynomial,
then it is linear in A and its gradient (with respect to A) .gradient.F
does not depend on A. That is, f is a polynomial in (x.sub.i, y.sub.i)
of the form
##EQU7##
p is the length of A (i.e. the total number of parameters a in the model)
and p.sub.x is the number of parameters pertaining to the x-component
f.sub.x of the model. The corresponding terms in the gradient .gradient. F
are
##EQU8##
Therefore, equation (17) holds with equality and B can b estimated
independent of A, allowing subarea selection using off-line computation
on reference data alone.
After some manipulation, the optimization criterion becomes
##EQU9##
where b.sub.ij is the (i,j)th element of B, p is the length of A (i.e.
the number of parameters a in the model),
##EQU10##
and (f.sub.x, f.sub.y) are the x and y components of f.
The term r.sub.ij, in effect, weights the parameter covariance terms
according to their total contribution to pixel registration covariances.
Those weights only need to be computed once for a specified image size and
model order. The b.sub.ij 's are recomputed for each candidate S.
A better understanding of the present invention may be gained from example,
and several examples are presented below.
Example 1: Affine f, N-by-N Image
The components of f may be written as
f.sub.x =a.sub.1 +a.sub.2 x+a.sub.3 y (23)
f.sub.y =a.sub.4 +a.sub.5 x+a.sub.6 y. (24)
Then the partial derivatives are
##EQU11##
After collecting terms, the optimization criterion becomes:
J=(b.sub.11 +b.sub.44)N.sup.2 +(b.sub.22 +b.sub.55) x.sup.2 +(b.sub.33
+b.sub.66) y.sup.2 +2(b.sub.12 +b.sub.42) x+2(b.sub.13 +b.sub.46).sup. y(
30)
where
##EQU12##
and similarly for y.
Example 2: Pre-computed Weights for .OMEGA. Symmetric about the Origin,
Unit Sample Spacing, Affine f
Symmetry implies:
x= y=0 (33)
and unit sample spacing implies:
##EQU13##
so that
J=(b.sub.11 +b.sub.44)N.sup.2 +(b.sub.22 +b.sub.33 +b.sub.55
+b.sub.66)N.sup.4 /12. (35)
Equation (35) only involves the diagonal terms of B. Care should be taken
to use the same units for the x's used in computing B and those used for
pre-computing the inner sum.
Use of a priori parameter estimates. Suppose that we normally have initial
estimates of registration parameters and could quantify the accuracy of
those estimates or, at least bound their uncertainty. For example, in a
second order polynomial model, one might know from operational limits on
viewing geometry that the covariance of the quadratic terms are much
smaller than that of all lower order terms. Then the criterion function
(17) may inappropriately be dominated by the effects of the covariance for
the quadratic terms. This can be systematically prevented by including an
a priori covariance in the formulation, as described below.
Let the initial estimate be A.sub.0 with covariance B.sub.0. Then a
Bayesian estimate A of the parameters A is given by
A=(B.sub.0.sup.-1 +B.sup.-1).sup.-1 (B.sub.0.sup.-1 A.sub.0 +B.sup.-1
A).sup.-1 =(B.sub.0.sup.-1 +B.sup.-1).sup.-1 (B.sub.0.sup.-1 A.sub.0
+.gradient. F.sup.T C.sup.-1 X').sup.-1 (36)
where X' is the actual measurement vector.
The same subarea selection procedure applies but with a simple
modification: use of the covariance matrix:
B=cov A=(B.sub.0.sup.-1 +.gradient. F.sup.T C.sup.-1 .gradient. F).sup.-1(
37)
in place of B.
Example 3
Suppose that the covariance matrices are diagonal, B.sub.0 =diag(b.sub.0k,
k=1, . . . , p) and B=diag(b.sub.kk ; k=1, . . . , p). Then (32) becomes
J=(b.sub.11 +b.sub.44)N.sup.2 +(b.sub.22 +b.sub.33 +b.sub.55
+b.sub.66)N.sup.4 12 (38)
with
##EQU14##
This places a limit on the influence of the kth parameter according to the
following inequality:
b.sub.kk <b.sub.0k. (40)
Synthesis of optimal subarea sets. Given the means for computing J, the
problem is to synthesize sets of subareas that minimize J. The effect on J
of eliminating or adding a subarea can be dependent on all other subareas,
via .gradient. F. Therefore combinatorial optimization is used for
synthesizing the best S. If there are a total of K.sub.tot possible
subareas from which to select the best K, a brute-force search would
require
K.sub.tot !(K.sub.tot K)!K! (41)
evaluations of J. For typical values of K.sub.tot =1024 and K=16, this
would require over 10.sup.34 evaluations of J. Brute-force optimization is
clearly not feasible.
However, efficient search procedures can dramatically reduce the
computational requirements. Backtrack programming may provide a method to
solve for the optimum subset as is described in "Backtrack Programming",
by S. W. Golomb et al., in J. ACM, Vol. 12, pp. 516-524, October 1965. Two
suboptimal solution approaches are described below that can reduce the
cost to O(K.sup.2.sub.tot) or O(K.sub.tot) evaluations of J. They are
suboptimal because they do not include backtracking in the search
procedure.
Procedure 1
Start with maximal set of subareas and eliminate those that minimize the
increase in total MSE.
Initialization is accomplished by: (a) calculating downsampled information
matrix images, (b) calculating K.sub.tot =number of points in downsampled
image. (c) look up or calculating covariance propagation weights r, (d)
constructing S.sub.Ktot, t={1, . . . , K.sub.tot }, and (e) calculating
x.sub.k =position in original image of the kth point from the downsampled
image. For k=K.sub.tot through K in increments of K.sub.inc : (a) for
every active subarea i.epsilon.S.sub.k ; calculate J(S.sub.k -i), and (b)
deactivate the K.sub.inc active subareas with largest J(S.sub.k -i).
If K.sub.inc =1, J(S.sub.k -i) is recalculated from scratch each time and
brute force Gaussian elimination is used for calculating B, this procedure
requires O(p.sup.3 K.sub.tot .sup.3) calculations where p is the number of
registration parameters. Accelerated calculation where K.sub.inc =kB 2,
i.e. half of the active subareas are eliminated at each step, reduces this
to O(p.sup.3 K.sub.tot.sup.2 log K.sub.tot. Recursive calculation of
J(S.sub.k -i) reduces calculation by a factor of about K.sub.tot.
Efficient matrix inversion reduces this by a factor of p and, combined
with the above savings, results in O(p.sup.3 K.sub.tot.sup.3) log
K.sub.tot) computational cost.
Procedure 2
Start with a minimal set of subareas and then add the subareas that
maximize the reduction of total MSE. Initialization is accomplished by:
(a) calculate downsampled information matrix images, (b) calculate
K.sub.tot =number of points in downsampled image, (c) look up or calculate
covariance propagation weights r.sub.ij, (d) calculate x.sub.k =position
in original image of the kth point from the downsampled image, (e) For
every subarea i.epsilon.{1, . . . , K.sub.tot }: calculate local MSE
tr(C.sub.i), and (f) activate the K.sub.min .about.pB 2 subareas with
smallest tr(Ci). For k=K.sub.min through K in increments of K.sub.inc :
(a) for every inactive subarea i.epsilon.S.sub.k : calculate J(S.sub.k
+i), and (b) activate the K.sub.inc active elements with smallest
J(S.sub.k +i).
Example 4
N-by-N SAR image with N.sub.looks looks averaged. If w is a rectangular
window with area
A.sub.w =W.sup.2 (42)
with
##EQU15##
The information images, after decimation, contain less than
K.sub.tot =[2N/W].sup.2 (44)
samples, depending on implementation.
Example 5
Design for 512-by-512 SAR image with 4 looks averaged. If
W=32 (45)
A.sub.w =32.sup.2.sub.-- 1024 (46)
K.sub.tot [2.times.512/32].sup.2 =1024 (47)
and the desired number of selected subareas is
K=16. (48)
Procedure 1 may be computationally prohibitive with the above parameters
if, for example, K.sub.inc =1. The previously-mentioned acceleration
strategy selects a decreasing sequence where the number of active subareas
is halved at each iteration. Specifically,
K.sub.inc ={512, 256,128, 64, 32, 16}. (49)
Experimental Results
Procedure 2 was implemented and applied to SAR imagery. Four different s | | |
