- Email: [email protected]
Least absolute deviation (LAD) image matching
PHOTOGRAMMETRY & REMOTESENSING ELSEVIER
ISPRS Journal of Photogrammetry & Remote Sensing 51 (1996) 223-229
Least absolute deviation (LAD) image matching M . E Calitz, H. Rtither* Department of Surveying and Geodetic Engineering, University of Cape Town, Cape Town, South Africa Received 13 June 1995; accepted 16 October 1995
Abstract The robust estimator properties of the L~-norm or least absolute deviation (LAD) is shown to provide better subpixel matching accuracy in the presence of outlier points than the least squares method widely employed for image matching applications. Two LAD algorithms are compared with each other and with the least squares (LS) method and the iteratively reweighted least squares (IRLS) method. Results indicate that the Barrodale-Roberts LAD algorithm can be used advantageously in conjunction with or in place of the IRLS and LS algorithms.
1. Introduction Most digital image matching strategies rely on techniques for initial localization of matching areas or points to pixel accuracy. A method to obtain matching to subpixel accuracy is then applied to the candidate matching points followed by a final acceptance or rejection of the match. The most widely accepted technique for matching to subpixel accuracy is the iteratively reweighted least squares (IRLS) grey-scale correlation method (Frrstner, 1984; Gruen, 1985). As shown in a previous paper (Calitz and RUther, 1994), the Lt-norm or least absolute deviation (LAD) can be used to improve the accuracy of the estimation problem in the presence of outlier points in the data. The LAD algorithm used in the latter investigation was based on the method of Barrodale and Young (1966), hereafter designated by (BY). This algorithm was found to be subject to numerical instabilities and tended to converge rather slowly. * Corresponding author. Elsevier Science B.V. SSDI 0 9 2 4 - 2 7 1 6 ( 9 5 ) 0 0 0 1 0 - 0
As an improvement on this algorithm the Barrodale and Roberts (1973, 1974), hereafter designated as BR, implementation was adapted to the correlation problem and found to converge stably. The following investigation compares the BR and BY algorithms with the least squares (LS) and IRLS methods and attempts to establish the conditions under which the BR algorithm is a more suitable alternative to the IRLS method. 2. Theory Digital images can be regarded as discrete twodimensional grey-scale distribution functions. Greyscale correlation matching selects a template window in one image and attempts to find the corresponding patch in another image which exhibits the greatest similarity. Usual measures of similarity are the grey-scale correlation and the grey-scale difference between the windows. In Lp-norm matching (0 < p < c~), with least squares being the L2-norm (Watson, 1980), the template f (x, y) is assumed to map into the patch g (x, y) by means of an affine
M.E Calitz, H. RCither/ISPRS Journal of Photogrammetry & Remote Sensing 51 (1996) 223-229
224
transformation (Gruen, 1985): X = all
is the parameter adjustment vector. The 'observations' are
--}-a12Xt -F-a21j
y = bll + bl2x' + b21y'
(1)
where x, y and x', y' are the respective coordinates in the template and patch; and aij are the affine transformation parameters. Lp-norm matching attempts to minimize the following expression involving the residuals between the template and the patch:
(EIrlp)I/P:=~ min
(2)
li ---- f
(xi,
Yi) --
g (xi, Yi)
and the rows of the design matrix are given by Ai -~
(gxl, gxiXi'
gxl Yi, gYi' gYl Xi' gyiYi)
for the ith point in the patch. The least squares solution with p = 2 is described in Gruen (1985). At each iteration the solution for x is given by x = (ATWA) -1 ATWl
(5)
where r = f (x, y) - g (x', y') is the residual corresponding to the template point (x, y). The point (x, y) is related to the patch coordinates by Eq. 1. Eq. 2 is essentially a non-linear minimization problem. As most matching strategies find approximate positions for the matching points, and assuming that the grey-scale distribution function is continuous in the local area of the matching points, Eq. 2 can be linearized by means of a first-order Taylor series expansion and reformulated in terms of the differentials of the affine transformation coefficients.
where A is the m x 6 matrix with each row given by Ai above; 1 is the m vector of li; and W is an m x m diagonal weight matrix. The weight matrix in this context is different to the one derived from observation quality. For the pure least squares (LS) method the weight matrix becomes the identity matrix. For the IRLS method the diagonal terms of the weight matrix are usually given as (Kubik et al., 1987)
( E IrlP) I/p = ( E
where ~b(ri) is an estimator function of the residuals. Commonly used M-estimators are the Huber, Hampel, Andrews and Biweight estimators (Li, 1985). s is an arbitrarily small constant to avoid division by zero. At each iteration the weights are recalculated according to the current values of the residuals. The Lp-norm estimators for p ~ (1, 2) can be simulated using IRLS with weights in the form
f (x, y ) - [g (x, y)
+ gx dall + gxX' dal2 + gxY' dazl + gy db11 + gyX' dbl2 + gyy' db21] p)l/p
(3)
where Og (x, y)
gx --
Ox
oqg (x, y)
'
gY =
Oy
The above linearized form of the minimization problem is an application of the Ganss-Newton method for solving non-linear problems, and is solved iteratively. The affine parameters are assigned initial values, and are updated at each iteration. The affine parameter corrections daij and dbij are calculated at each iteration until convergence is reached. Convergence is assumed when all the adjustments are less than a predefined value. Eq. 3 can be restated as
(EIFiIP)'/P=(EIAix-lilP)UP
(4)
where the summation is over all points in the patch. x T = (dall, dal2, da21, dbll, dbl2, db21)
dp(ri) Wi -- r~ + e
1 Wi -- Iril 2-p -~ e"
(6)
(7)
With p -- 1 Eq. 7 corresponds to the LAD estimator. The above approach of using the algorithmic solution of the least squares problem to simulate robust estimators is used in most applications due to its simplicity. No changes to the algorithmic structure are required other than altering the weights. The theory of least squares estimation is well established with regard to evaluating the results of the estimation procedure. The principal justification for attempting to use robust estimation techniques is to reduce the influence of outlier points on the results of the estimation.
M.E Calitz, H. Riither/lSPRS Journal of Photogrammetry & Remote Sensing 51 (1996) 223-229
(a)
(b)
Fig. 1. A template (left) and patch (right) containing the roof point of a building.
Least squares tends to become less reliable in the presence of outliers arising from noise between the template and patch. Outliers can also arise due to the inadequacy of the affine transformation assumption at discontinuities in the image due to sharp changes in height. Examples in aerial images would be the change in height from ground level to the roof of a building and the steep height difference at cliff faces. There are many examples in close range applications arising from the relatively large distance changes at the edges of objects in the scene. Different affine mapping parameters would apply to the various piecewise continuous regions that arise. These discontinuities occur at important matching primitive points such as corners and edges and are crucial to the eventual digital elevation model produced. An alternate mapping function between the template and patch can be derived from the collinearity condition equations and can be expressed as al + a2x' + a3y' 1 + a4x' + asy' bl + b2x' + b3y' (8) Y = l + b4x' + bs y' The coefficients in Eq. 8 are a function of the relative orientation between the two imaging devices and the distance to the objects in the scene. In this formulation a ten-parameter estimation problem needs to be solved. Using the IRLS with the Gauss-Newton formulation of the problem leads to unstable results due to the sensitivity of the mapping to the values of the denominators in Eq. 8. An alternative non-linear technique is required that retains higher-order terms. Analogous to the affine mapping, the above function will be piecewise continuous. Due to the numerical instability and the spurious accuracy that may X~
225
result from using an alternate mapping, the affine transformation seems to be justified. Fig. 1 shows a template and patch containing the roof point of a building. The affine mapping for the roof area will be different to that of the ground area. In solving for the affine parameters using the above-described formulation of the matching problem, the values determined will be those which shape the patch so that the residuals are minimized in a least squares sense. The effective mapping resulting from this model will then depend on the relative number of points in each piecewise continuous region. According to approximation theory (Watson, 1980), the ordinary least squares method will tend to solve for the mean mapping over the region. LAD, on the other hand, will solve for the median mapping over the region, which will correspond to the affine parameters of the continuous region with more than 50% of the points in the window. This would suggest that if the central points of the template and patch which are to be matched fall in or on the boundary of the largest continuous region then LAD matching should be more accurate. The roof point in Fig. 1 falls on the boundary of the ground area, which is the largest continuous area. A previous investigation (Calitz and Rtither, 1994), where artificially generated outliers were introduced into images, showed that LAD estimation was more precise than ordinary LS. IRLS seemed to give equally good results but the number of iterations required is greater than that for the BY algorithm used for LAD matching. Extensive comparisons of IRLS with linear programming algorithms such as BY and BR have been conducted by Armstrong and Frome (1976) (Bloomfield and Steiger, 1983). They conclude that the BR algorithm is generally faster and numerically more stable than the IRLS method. The numerical experiments described below tend to support this conclusion. For LAD matching the minimization problem described by Eq. (2) needs to be solved with p = 1. Let ri = ui - vi with ui, vi >_ O. ui and vi are non-negative slack variables introduced to recast the equations as a primal linear programming problem as follows: E
(ui + vi) ::~ rain
(9)
226
M.E Calitz, H. Riither / ISPRS Journal of Photogrammetry & Remote Sensing 51 (1996) 223-229 1 0,90.80.70.6-
~ 0.5~ 0.4~ 0.3" Z 0.2-
1
• BR t~ OLSM_
+ IRLS *BY
0.1.
ITERATION No.
Fig. 2, Convergence of the algorithms when matching a point with itself.
subject to the constraints li : A i X -I- Ui -- Ui.
(10)
The BY (Barrodale and Young, 1966) or BR (Barrodale and Roberts, 1974) algorithm is used to solve for x at each iteration. 3. N u m e r i c a l e x p e r i m e n t s
Numerical experiments were undertaken to compare the accuracy and convergence behaviour of the
four algorithms, viz. the - ordinary least-squares method (OLSM), - iteratively re-weighted least squares method (IRLS), - Barrodale-Roberts Ll-norm method (BR), - Barrodale-Young Ll-norm method (BY). The results of these experiments act as indicators of the likely general performance. All the algorithms took approximately the same time per iteration, with the BY algorithm being marginally slower. Fig. 2 shows the convergence behaviour of the algorithms when matching a point with itself. As expected all converge uniformly although the IRLS technique exhibits oscillations, and takes more iterations than the other methods. Twenty targets were matched from the left image to the fight image in Fig. 3, to compare the subpixel localization of target centres. The centre of gravity of the targets was chosen as the 'true' centre, although this is not necessarily the case due to radiometric gradients across the targets, random radiometric noise and the threshold value chosen to isolate the targets. The average results for the twenty targets are shown in Tables 1 and 2. For these matching experiments a resampled template is used with the centre of gravity of the target at its centre. The initial approximating patch in the right image is similarly based on the centre of gravity of
Fig. 3. Stereopair to test target centering using L-norm matching.
M.F. Calitz, H. Riither / ISPRS Journal of Photogrammetry & Remote Sensing 51 (1996) 223-229
Table 1 Average number of interations and pixel deviation e for the various L-norm algorithms at various tolerance values. Average taken over 20 targets matchedfrom Fig. 3
Table 3 Average number of iterations and pixel deviation e for the various L-norm algorithms at various tolerance values. Average taken over 10 selectedpoints from Fig. 4 Tolerances
Tolerances 0.1
OLSM IRLS BY BR
227
0.01
0.1
0.001
e
lter
e
Iter
E
Iter
0.08 0.21 0.15 0.12
4.0 5.9 4.3 3.5
0.09 0.13 0.10 0.10
5.7 14.6 8.4 6.4
0.09 0.10 0.11 0.11
14.0 28.9 15.8 8.9
Table 2 Average differences in pixels between the matching points determined by each method Tolerance
1-2
1-3
1--4
2-3
2-4
3-4
0.1 0.01 0.001
0.48 0.22 0.18
0.08 0.11 0.12
0.10 0.15 0.14
0.47 0.21 0.19
0.47 0.23 0.16
0.07 0.07 0.05
1 = OLSM, 2 = IRLS, 3 = BR, 4 = BY.
the corresponding target. A window size of 15 x 15 pixels was used. The initial matching therefore commences from the best possible approximating point. The iterations terminate when the affine parameter adjustments are less than the indicated tolerance value. The results indicate that the ordinary least squares method gives the best results in terms of precision and speed at the low tolerance of 0.1. At a tolerance of 0.001 the difference in the pixel deviations between the various methods are within 0.2 of a pixel. At the higher tolerance the BR algorithm converges significantly faster than the other methods. The IRLS algorithm appears to simulate the L1norm well at tolerance values of 0.001 or less, but it takes considerably more iterations, making it less efficient than the BR or BY method. Due to the absence of clearly defined outlier points, the ordinary least squares method would be expected to give the most precise results, as confirmed by the tests. A similar test was performed using the images shown in Fig. 4. Ten comer points were selected to pixel accuracy from each image. The results of this comparison are shown in Table 3. Clearly all the methods give poor results, with the OLSM and BR
OLSM IRLS BY BR
0.01
0.001
e
Iter
e
lter
e
Iter
4.26 11.0 6.94 3.52
24.7 17.1 30.4 5.6
4.75 9.36 10.56 4.52
38.1 39.5 45.9 24.9
4.78 5.55 8.25 5.02
45.0 50.0 46.4 40.8
methods being the most stable whereas the IRLS and BY methods are the most unstable. The accuracy in all cases is very poor. The relationship between the number of iterations and the tolerance level is similar to that of the previous experiment. The poor results can be ascribed to a number of factors, one of which is the low contrast at discontinuities. Contrast enhancement of the images only improved the results of a few points, but the overall accuracy remained low. The significant effects of the changed perspective views at certain comers also contribute to the poor results. This experiment indicates that subpixel areabased techniques are only effective in uniformly varying grey-scale regions and in the absence of discontinuities. The affine mapping reshaping model is generally not appropriate in the presence of discontinuities. A further experiment was conducted to illustrate the circumstances under which the Ll-norm algorithms give better results than the OLSM method. Outlier points were added to the 20 right image targets in Fig. 2 that were used in the target test described above. The target diameter was about eight pixels and the outliers were in the form of circles with radii of two pixels of a significantly different grey-scale value to the original targets. The results of this experiment are shown in Tables 4 and 5. The accuracy of all methods is affected by the presence of the outliers, but the pixel deviation of the OLSM is about one pixel greater than that of the other methods. The pixel deviations of the Ll-norm methods is less than or equal to 0.25 of a pixel at a tolerance of 0.001. These results indicate that the BR algorithm is
228
M.F. Calitz, H. Riither/ lSPRS Journal of Photogrammetry & Remote Sensing 51 (1996) 223-229
Fig. 4. Stereopair of house used to compare l-norm algorithms as a real image test case. Table 4 Average number of interations and pixel deviation e for the various L-norm algorithms at various tolerance values. Average taken over 20 targets matched from Fig. 3, with outliers added to the right image targets Tolerances 0.1 OLSM IRLS BY BR
0.01
0.001
Table 5 Average difference in pixels between the matching points determined by each method Tolerance
1-2
1-3
1-4
2-3
2-4
3-4
0.1 0.01 0.001
0.73 0.48 0.88
0.28 0.33 0.33
0.30 0.31 0.34
0.63 0.38 0.73
0.64 0.38 0.75
0.11 0.07 0.07
e
Iter
e
Iter
e
Iter
1 = OLSM, 2 = IRLS, 3 = BR, 4 = BY.
0.93 0.57 0.33 0.32
4.1 8.2 6.2 4.6
1.59 0.26 0.27 0.30
26.3 21.1 26.5 7.7
1.57 0.20 0.24 0.25
42.4 39.7 44.1 26.9
4. Conclusions
clearly the most efficient algorithm for subpixel matching in the presence of significant outliers. The number of iterations required is considerably less than that for the other methods while the accuracy is similar to that of the IRLS and BY algorithms. The unstable nature of the IRLS method is brought out by the rather large differences shown in Table 5. Although the matching points are close to the centre of gravity (Table 4), they can still be further from the matching positions found by the BR and BY methods. This can occur if many of the IRLS matching points are on the diametrically opposite side of the centre of gravity to the BR and BY matching points.
The main conclusions that can be drawn from the numerical experiments are: (1) The Barrodale-Roberts (BR) least absolute deviation (LAD) algorithm is as stable as the ordinary least-squares method (OLSM). (2) It is computationally faster and gives similar results to the OLSM in the absence of outliers. (3) All the Ll-norm methods investigated give significantly more accurate results than the O L S M in the presence of outliers. (4) When matching t h e comers of buildings between images with considerably different perspective views all the methods give poor results. Area-based methods are not appropriate in these circumstances and a feature-based approach must be adopted to
M.F. Calitz, H. Riither / ISPRS Journal of Photogrammetry & Remote Sensing 51 (1996) 223-229 locate the c o m e r points to subpixel accuracy. Further investigations are required in these situations.
Acknowledgements The research in digital i m a g e m a t c h i n g techniques is f u n d e d by the F o u n d a t i o n for Research and Dev e l o p m e n t (FRD) and the University of Cape Town. T h e i r contribution is gratefully acknowledged. The images o f buildings shown in Fig. 4 were provided by the E T H in Switzerland.
References Armstrong, R.D. and Frome, E.L., 1976. A comparison of two algorithms for absolute deviation curve fitting. J. Am. Stat. Assoc., 71: 328-330. Barrodale, I. and Roberts, F.D.K., 1973. An improved algorithm for discrete II linear approximation. SIAM J. Numer. Anal., 10(5): 839-848. Barrodale, I. and Roberts, F.D.K., 1974. Algorithm 478: Solution
229
of an overdetermined system of equations in the lj norm. Commun. ACM, 17(6): 319-320. Barrodale, I. and Young, A., 1966. Algorithms for best LI and Lo~ linear approximations on a discrete set. Numerische Mathematik, 8: 295-306. Bloomfield, P. and Steiger, W.L., 1983. Least Absolute Deviations. Birkh~iuser,Boston, Basel, Stuttgart. Calitz, M.F. and Rtither, H., 1994. L-Norm methods in areabased matching, IAPRS. Frrsmer, W., 1984. Quality assessment of object location and point transfer using digital image correlation techniques. IAPRS, Vol. XXV, A3a, Commission III, pp. 197-219. Gruen, A., 1985. Adaptive least-squares correlation - a powerful matching technique. S. Afr. J. Photogramm. Remote Sensing Cartography, 14(3): 175-187. Kubik, K., Merchant, D. and Schenk, T., 1987. Robust regression in photogrammetry. Photogramm. Eng. Remote Sensing, 53(2): 167-169. Li, Gouying, 1985. Robust regression. In: Hoaglin, D.C.. Mosteller and E, Tukey, J.W. (Editors), Exploring Data. Tables, Trends, and Shapes. John Wiley and Sons, New York, pp. 281-344. Watson, G.A., 1980. Approximation Theory and Numerical Methods. John Wiley and Sons, Chichester.
ISPRS Journal of Photogrammetry & Remote Sensing 51 (1996) 223-229
Least absolute deviation (LAD) image matching M . E Calitz, H. Rtither* Department of Surveying and Geodetic Engineering, University of Cape Town, Cape Town, South Africa Received 13 June 1995; accepted 16 October 1995
Abstract The robust estimator properties of the L~-norm or least absolute deviation (LAD) is shown to provide better subpixel matching accuracy in the presence of outlier points than the least squares method widely employed for image matching applications. Two LAD algorithms are compared with each other and with the least squares (LS) method and the iteratively reweighted least squares (IRLS) method. Results indicate that the Barrodale-Roberts LAD algorithm can be used advantageously in conjunction with or in place of the IRLS and LS algorithms.
1. Introduction Most digital image matching strategies rely on techniques for initial localization of matching areas or points to pixel accuracy. A method to obtain matching to subpixel accuracy is then applied to the candidate matching points followed by a final acceptance or rejection of the match. The most widely accepted technique for matching to subpixel accuracy is the iteratively reweighted least squares (IRLS) grey-scale correlation method (Frrstner, 1984; Gruen, 1985). As shown in a previous paper (Calitz and RUther, 1994), the Lt-norm or least absolute deviation (LAD) can be used to improve the accuracy of the estimation problem in the presence of outlier points in the data. The LAD algorithm used in the latter investigation was based on the method of Barrodale and Young (1966), hereafter designated by (BY). This algorithm was found to be subject to numerical instabilities and tended to converge rather slowly. * Corresponding author. Elsevier Science B.V. SSDI 0 9 2 4 - 2 7 1 6 ( 9 5 ) 0 0 0 1 0 - 0
As an improvement on this algorithm the Barrodale and Roberts (1973, 1974), hereafter designated as BR, implementation was adapted to the correlation problem and found to converge stably. The following investigation compares the BR and BY algorithms with the least squares (LS) and IRLS methods and attempts to establish the conditions under which the BR algorithm is a more suitable alternative to the IRLS method. 2. Theory Digital images can be regarded as discrete twodimensional grey-scale distribution functions. Greyscale correlation matching selects a template window in one image and attempts to find the corresponding patch in another image which exhibits the greatest similarity. Usual measures of similarity are the grey-scale correlation and the grey-scale difference between the windows. In Lp-norm matching (0 < p < c~), with least squares being the L2-norm (Watson, 1980), the template f (x, y) is assumed to map into the patch g (x, y) by means of an affine
M.E Calitz, H. RCither/ISPRS Journal of Photogrammetry & Remote Sensing 51 (1996) 223-229
224
transformation (Gruen, 1985): X = all
is the parameter adjustment vector. The 'observations' are
--}-a12Xt -F-a21j
y = bll + bl2x' + b21y'
(1)
where x, y and x', y' are the respective coordinates in the template and patch; and aij are the affine transformation parameters. Lp-norm matching attempts to minimize the following expression involving the residuals between the template and the patch:
(EIrlp)I/P:=~ min
(2)
li ---- f
(xi,
Yi) --
g (xi, Yi)
and the rows of the design matrix are given by Ai -~
(gxl, gxiXi'
gxl Yi, gYi' gYl Xi' gyiYi)
for the ith point in the patch. The least squares solution with p = 2 is described in Gruen (1985). At each iteration the solution for x is given by x = (ATWA) -1 ATWl
(5)
where r = f (x, y) - g (x', y') is the residual corresponding to the template point (x, y). The point (x, y) is related to the patch coordinates by Eq. 1. Eq. 2 is essentially a non-linear minimization problem. As most matching strategies find approximate positions for the matching points, and assuming that the grey-scale distribution function is continuous in the local area of the matching points, Eq. 2 can be linearized by means of a first-order Taylor series expansion and reformulated in terms of the differentials of the affine transformation coefficients.
where A is the m x 6 matrix with each row given by Ai above; 1 is the m vector of li; and W is an m x m diagonal weight matrix. The weight matrix in this context is different to the one derived from observation quality. For the pure least squares (LS) method the weight matrix becomes the identity matrix. For the IRLS method the diagonal terms of the weight matrix are usually given as (Kubik et al., 1987)
( E IrlP) I/p = ( E
where ~b(ri) is an estimator function of the residuals. Commonly used M-estimators are the Huber, Hampel, Andrews and Biweight estimators (Li, 1985). s is an arbitrarily small constant to avoid division by zero. At each iteration the weights are recalculated according to the current values of the residuals. The Lp-norm estimators for p ~ (1, 2) can be simulated using IRLS with weights in the form
f (x, y ) - [g (x, y)
+ gx dall + gxX' dal2 + gxY' dazl + gy db11 + gyX' dbl2 + gyy' db21] p)l/p
(3)
where Og (x, y)
gx --
Ox
oqg (x, y)
'
gY =
Oy
The above linearized form of the minimization problem is an application of the Ganss-Newton method for solving non-linear problems, and is solved iteratively. The affine parameters are assigned initial values, and are updated at each iteration. The affine parameter corrections daij and dbij are calculated at each iteration until convergence is reached. Convergence is assumed when all the adjustments are less than a predefined value. Eq. 3 can be restated as
(EIFiIP)'/P=(EIAix-lilP)UP
(4)
where the summation is over all points in the patch. x T = (dall, dal2, da21, dbll, dbl2, db21)
dp(ri) Wi -- r~ + e
1 Wi -- Iril 2-p -~ e"
(6)
(7)
With p -- 1 Eq. 7 corresponds to the LAD estimator. The above approach of using the algorithmic solution of the least squares problem to simulate robust estimators is used in most applications due to its simplicity. No changes to the algorithmic structure are required other than altering the weights. The theory of least squares estimation is well established with regard to evaluating the results of the estimation procedure. The principal justification for attempting to use robust estimation techniques is to reduce the influence of outlier points on the results of the estimation.
M.E Calitz, H. Riither/lSPRS Journal of Photogrammetry & Remote Sensing 51 (1996) 223-229
(a)
(b)
Fig. 1. A template (left) and patch (right) containing the roof point of a building.
Least squares tends to become less reliable in the presence of outliers arising from noise between the template and patch. Outliers can also arise due to the inadequacy of the affine transformation assumption at discontinuities in the image due to sharp changes in height. Examples in aerial images would be the change in height from ground level to the roof of a building and the steep height difference at cliff faces. There are many examples in close range applications arising from the relatively large distance changes at the edges of objects in the scene. Different affine mapping parameters would apply to the various piecewise continuous regions that arise. These discontinuities occur at important matching primitive points such as corners and edges and are crucial to the eventual digital elevation model produced. An alternate mapping function between the template and patch can be derived from the collinearity condition equations and can be expressed as al + a2x' + a3y' 1 + a4x' + asy' bl + b2x' + b3y' (8) Y = l + b4x' + bs y' The coefficients in Eq. 8 are a function of the relative orientation between the two imaging devices and the distance to the objects in the scene. In this formulation a ten-parameter estimation problem needs to be solved. Using the IRLS with the Gauss-Newton formulation of the problem leads to unstable results due to the sensitivity of the mapping to the values of the denominators in Eq. 8. An alternative non-linear technique is required that retains higher-order terms. Analogous to the affine mapping, the above function will be piecewise continuous. Due to the numerical instability and the spurious accuracy that may X~
225
result from using an alternate mapping, the affine transformation seems to be justified. Fig. 1 shows a template and patch containing the roof point of a building. The affine mapping for the roof area will be different to that of the ground area. In solving for the affine parameters using the above-described formulation of the matching problem, the values determined will be those which shape the patch so that the residuals are minimized in a least squares sense. The effective mapping resulting from this model will then depend on the relative number of points in each piecewise continuous region. According to approximation theory (Watson, 1980), the ordinary least squares method will tend to solve for the mean mapping over the region. LAD, on the other hand, will solve for the median mapping over the region, which will correspond to the affine parameters of the continuous region with more than 50% of the points in the window. This would suggest that if the central points of the template and patch which are to be matched fall in or on the boundary of the largest continuous region then LAD matching should be more accurate. The roof point in Fig. 1 falls on the boundary of the ground area, which is the largest continuous area. A previous investigation (Calitz and Rtither, 1994), where artificially generated outliers were introduced into images, showed that LAD estimation was more precise than ordinary LS. IRLS seemed to give equally good results but the number of iterations required is greater than that for the BY algorithm used for LAD matching. Extensive comparisons of IRLS with linear programming algorithms such as BY and BR have been conducted by Armstrong and Frome (1976) (Bloomfield and Steiger, 1983). They conclude that the BR algorithm is generally faster and numerically more stable than the IRLS method. The numerical experiments described below tend to support this conclusion. For LAD matching the minimization problem described by Eq. (2) needs to be solved with p = 1. Let ri = ui - vi with ui, vi >_ O. ui and vi are non-negative slack variables introduced to recast the equations as a primal linear programming problem as follows: E
(ui + vi) ::~ rain
(9)
226
M.E Calitz, H. Riither / ISPRS Journal of Photogrammetry & Remote Sensing 51 (1996) 223-229 1 0,90.80.70.6-
~ 0.5~ 0.4~ 0.3" Z 0.2-
1
• BR t~ OLSM_
+ IRLS *BY
0.1.
ITERATION No.
Fig. 2, Convergence of the algorithms when matching a point with itself.
subject to the constraints li : A i X -I- Ui -- Ui.
(10)
The BY (Barrodale and Young, 1966) or BR (Barrodale and Roberts, 1974) algorithm is used to solve for x at each iteration. 3. N u m e r i c a l e x p e r i m e n t s
Numerical experiments were undertaken to compare the accuracy and convergence behaviour of the
four algorithms, viz. the - ordinary least-squares method (OLSM), - iteratively re-weighted least squares method (IRLS), - Barrodale-Roberts Ll-norm method (BR), - Barrodale-Young Ll-norm method (BY). The results of these experiments act as indicators of the likely general performance. All the algorithms took approximately the same time per iteration, with the BY algorithm being marginally slower. Fig. 2 shows the convergence behaviour of the algorithms when matching a point with itself. As expected all converge uniformly although the IRLS technique exhibits oscillations, and takes more iterations than the other methods. Twenty targets were matched from the left image to the fight image in Fig. 3, to compare the subpixel localization of target centres. The centre of gravity of the targets was chosen as the 'true' centre, although this is not necessarily the case due to radiometric gradients across the targets, random radiometric noise and the threshold value chosen to isolate the targets. The average results for the twenty targets are shown in Tables 1 and 2. For these matching experiments a resampled template is used with the centre of gravity of the target at its centre. The initial approximating patch in the right image is similarly based on the centre of gravity of
Fig. 3. Stereopair to test target centering using L-norm matching.
M.F. Calitz, H. Riither / ISPRS Journal of Photogrammetry & Remote Sensing 51 (1996) 223-229
Table 1 Average number of interations and pixel deviation e for the various L-norm algorithms at various tolerance values. Average taken over 20 targets matchedfrom Fig. 3
Table 3 Average number of iterations and pixel deviation e for the various L-norm algorithms at various tolerance values. Average taken over 10 selectedpoints from Fig. 4 Tolerances
Tolerances 0.1
OLSM IRLS BY BR
227
0.01
0.1
0.001
e
lter
e
Iter
E
Iter
0.08 0.21 0.15 0.12
4.0 5.9 4.3 3.5
0.09 0.13 0.10 0.10
5.7 14.6 8.4 6.4
0.09 0.10 0.11 0.11
14.0 28.9 15.8 8.9
Table 2 Average differences in pixels between the matching points determined by each method Tolerance
1-2
1-3
1--4
2-3
2-4
3-4
0.1 0.01 0.001
0.48 0.22 0.18
0.08 0.11 0.12
0.10 0.15 0.14
0.47 0.21 0.19
0.47 0.23 0.16
0.07 0.07 0.05
1 = OLSM, 2 = IRLS, 3 = BR, 4 = BY.
the corresponding target. A window size of 15 x 15 pixels was used. The initial matching therefore commences from the best possible approximating point. The iterations terminate when the affine parameter adjustments are less than the indicated tolerance value. The results indicate that the ordinary least squares method gives the best results in terms of precision and speed at the low tolerance of 0.1. At a tolerance of 0.001 the difference in the pixel deviations between the various methods are within 0.2 of a pixel. At the higher tolerance the BR algorithm converges significantly faster than the other methods. The IRLS algorithm appears to simulate the L1norm well at tolerance values of 0.001 or less, but it takes considerably more iterations, making it less efficient than the BR or BY method. Due to the absence of clearly defined outlier points, the ordinary least squares method would be expected to give the most precise results, as confirmed by the tests. A similar test was performed using the images shown in Fig. 4. Ten comer points were selected to pixel accuracy from each image. The results of this comparison are shown in Table 3. Clearly all the methods give poor results, with the OLSM and BR
OLSM IRLS BY BR
0.01
0.001
e
Iter
e
lter
e
Iter
4.26 11.0 6.94 3.52
24.7 17.1 30.4 5.6
4.75 9.36 10.56 4.52
38.1 39.5 45.9 24.9
4.78 5.55 8.25 5.02
45.0 50.0 46.4 40.8
methods being the most stable whereas the IRLS and BY methods are the most unstable. The accuracy in all cases is very poor. The relationship between the number of iterations and the tolerance level is similar to that of the previous experiment. The poor results can be ascribed to a number of factors, one of which is the low contrast at discontinuities. Contrast enhancement of the images only improved the results of a few points, but the overall accuracy remained low. The significant effects of the changed perspective views at certain comers also contribute to the poor results. This experiment indicates that subpixel areabased techniques are only effective in uniformly varying grey-scale regions and in the absence of discontinuities. The affine mapping reshaping model is generally not appropriate in the presence of discontinuities. A further experiment was conducted to illustrate the circumstances under which the Ll-norm algorithms give better results than the OLSM method. Outlier points were added to the 20 right image targets in Fig. 2 that were used in the target test described above. The target diameter was about eight pixels and the outliers were in the form of circles with radii of two pixels of a significantly different grey-scale value to the original targets. The results of this experiment are shown in Tables 4 and 5. The accuracy of all methods is affected by the presence of the outliers, but the pixel deviation of the OLSM is about one pixel greater than that of the other methods. The pixel deviations of the Ll-norm methods is less than or equal to 0.25 of a pixel at a tolerance of 0.001. These results indicate that the BR algorithm is
228
M.F. Calitz, H. Riither/ lSPRS Journal of Photogrammetry & Remote Sensing 51 (1996) 223-229
Fig. 4. Stereopair of house used to compare l-norm algorithms as a real image test case. Table 4 Average number of interations and pixel deviation e for the various L-norm algorithms at various tolerance values. Average taken over 20 targets matched from Fig. 3, with outliers added to the right image targets Tolerances 0.1 OLSM IRLS BY BR
0.01
0.001
Table 5 Average difference in pixels between the matching points determined by each method Tolerance
1-2
1-3
1-4
2-3
2-4
3-4
0.1 0.01 0.001
0.73 0.48 0.88
0.28 0.33 0.33
0.30 0.31 0.34
0.63 0.38 0.73
0.64 0.38 0.75
0.11 0.07 0.07
e
Iter
e
Iter
e
Iter
1 = OLSM, 2 = IRLS, 3 = BR, 4 = BY.
0.93 0.57 0.33 0.32
4.1 8.2 6.2 4.6
1.59 0.26 0.27 0.30
26.3 21.1 26.5 7.7
1.57 0.20 0.24 0.25
42.4 39.7 44.1 26.9
4. Conclusions
clearly the most efficient algorithm for subpixel matching in the presence of significant outliers. The number of iterations required is considerably less than that for the other methods while the accuracy is similar to that of the IRLS and BY algorithms. The unstable nature of the IRLS method is brought out by the rather large differences shown in Table 5. Although the matching points are close to the centre of gravity (Table 4), they can still be further from the matching positions found by the BR and BY methods. This can occur if many of the IRLS matching points are on the diametrically opposite side of the centre of gravity to the BR and BY matching points.
The main conclusions that can be drawn from the numerical experiments are: (1) The Barrodale-Roberts (BR) least absolute deviation (LAD) algorithm is as stable as the ordinary least-squares method (OLSM). (2) It is computationally faster and gives similar results to the OLSM in the absence of outliers. (3) All the Ll-norm methods investigated give significantly more accurate results than the O L S M in the presence of outliers. (4) When matching t h e comers of buildings between images with considerably different perspective views all the methods give poor results. Area-based methods are not appropriate in these circumstances and a feature-based approach must be adopted to
M.F. Calitz, H. Riither / ISPRS Journal of Photogrammetry & Remote Sensing 51 (1996) 223-229 locate the c o m e r points to subpixel accuracy. Further investigations are required in these situations.
Acknowledgements The research in digital i m a g e m a t c h i n g techniques is f u n d e d by the F o u n d a t i o n for Research and Dev e l o p m e n t (FRD) and the University of Cape Town. T h e i r contribution is gratefully acknowledged. The images o f buildings shown in Fig. 4 were provided by the E T H in Switzerland.
References Armstrong, R.D. and Frome, E.L., 1976. A comparison of two algorithms for absolute deviation curve fitting. J. Am. Stat. Assoc., 71: 328-330. Barrodale, I. and Roberts, F.D.K., 1973. An improved algorithm for discrete II linear approximation. SIAM J. Numer. Anal., 10(5): 839-848. Barrodale, I. and Roberts, F.D.K., 1974. Algorithm 478: Solution
229
of an overdetermined system of equations in the lj norm. Commun. ACM, 17(6): 319-320. Barrodale, I. and Young, A., 1966. Algorithms for best LI and Lo~ linear approximations on a discrete set. Numerische Mathematik, 8: 295-306. Bloomfield, P. and Steiger, W.L., 1983. Least Absolute Deviations. Birkh~iuser,Boston, Basel, Stuttgart. Calitz, M.F. and Rtither, H., 1994. L-Norm methods in areabased matching, IAPRS. Frrsmer, W., 1984. Quality assessment of object location and point transfer using digital image correlation techniques. IAPRS, Vol. XXV, A3a, Commission III, pp. 197-219. Gruen, A., 1985. Adaptive least-squares correlation - a powerful matching technique. S. Afr. J. Photogramm. Remote Sensing Cartography, 14(3): 175-187. Kubik, K., Merchant, D. and Schenk, T., 1987. Robust regression in photogrammetry. Photogramm. Eng. Remote Sensing, 53(2): 167-169. Li, Gouying, 1985. Robust regression. In: Hoaglin, D.C.. Mosteller and E, Tukey, J.W. (Editors), Exploring Data. Tables, Trends, and Shapes. John Wiley and Sons, New York, pp. 281-344. Watson, G.A., 1980. Approximation Theory and Numerical Methods. John Wiley and Sons, Chichester.