- Email: [email protected]
Ultimate resolution: a mathematical framework
Ultramlcroscopy 47 (1992) 298-306 North-Holland
Ultimate resolution: a mathematical framework A v a n d e n Bos Department of Apphed Physws, Delft Unwerstty of Technology, P 0 Box 5046, 2600 GA Delft, The Netherlands Recewed at Editorial Office 8 January 1992
A new resolution hmlt is proposed Classical hmlts, such as Raylelgh's, are a measure of the width of the component intensity distributions to be resolved The proposed hmlt, on the other hand, is expressed in the errors m the observed intensity distributions It describes how the resolvabdlty of the components depends on these errors It is shown how this new limit can be computed for specified component functions at a given distance In addition, a test for the resolvabdlty of components from a given set of observations is presented
I. Introduction T h e p u r p o s e of this p a p e r is to propose a new definition of optical or e l e c t r o n - o p t l c a l resolution This definition is i n t e n d e d to be objective a n d o p e r a t i o n a l In both respects, it is different from the classical Raylelgh definition [1] widely in use in optical a n d electron-optical l i t e r a t u r e [2-5] Raylelgh's d e f i n i t i o n is based on p r e s u m e d limitations to the resolving capablhtles of the h u m a n visual system Since Raylelgh's days, visual inspection has b e e n s u p p l e m e n t e d with intensity m e a s u r i n g i n s t r u m e n t a t i o n a n d digital c o m p u t i n g facilities T h e resulting resolving capabilities exceed those of the h u m a n visual system This jUStlties the q u e s t i o n as to what the limitations to resolution are u n d e r these new conditions Raylelgh's hmlt to the resolution of two superimposed a n d overlapping c o m p o n e n t Intensity distributions is, to a certain extent, a m e a s u r e of the width of the c o m p o n e n t s No m e n t i o n is m a d e of the p r e s e n c e of errors in the observations of the c o m b i n e d intensity d i s t r i b u t i o n In the definition p r o p o s e d in this paper, these errors are central G i v e n the amplitudes, the locations, the m e a s u r e m e n t points a n d the f u n c t i o n a l description of the c o m p o n e n t distributions, the set of all errors, or equivalently that of the error-dist u r b e d observations, is c o m p u t e d for which the
c o m p o n e n t s are resolvable T h u s the f o r m u l a t e d hmlt to resolution is the b o u n d a r y of this set This definition is also o p e r a t i o n a l since it has b e e n f o u n d to be relatively easy to decide on which side of this b o u n d a r y a given set of observations is located So, o n e can d e t e r m i n e beforeh a n d w h e t h e r or not the c o m p o n e n t s can be resolved from the available observations T h e described error limit to resolution is based on a p a r a m e t e r e s t i m a t i o n approach [6,7] It is supposed that the t w o - c o m p o n e n t f u n c t i o n a l m o d e l is fitted with respect to both locations a n d both a m p l i t u d e s to e r r o r - c o r r u p t e d observations T h e n the resolution of the locations or, equivalently, of the c o m p o n e n t functions achieved is u l t i m a t e since it is b a s e d o n detailed a priori knowledge a b o u t the class of functions conc e r n e d T h u s formulated, the p r o b l e m at h a n d is an example of fitting wetghted-sum models A w e i g h t e d - s u m m o d e l of o r d e r K is described by
g( x , or, fl) = a l h ( X, ill) +
+.Kh( x, 13~), (1)
with flk 4: tim If k ~ m, where a a n d fl are the vectors of the linear p a r a m e t e r s a k a n d the n o n linear p a r a m e t e r s fig, while g ( x , a, fl) a n d h ( x , /3~) are functions of the i n d e p e n d e n t variable x T h e n u m b e r K is called the order of ( l )
0304-3991/92/$05 00 © 1992 - Elsevier Science Pubhshers B V All rights reserved
A van den B o s /
Ulttmate resolution a mathemattcalframework
Weighted-sum models are, for example, used to model multi-exponential decay The parameters a and /3 are then the amplitudes and the decay constants, respectively For weighted-sum models, the computation of the error limit to resolution of the/3 k has been extensively studied [8-11] In this p a p e r the model (1) will be extended to component functions h(x, y, 13k) of two independent variables x and y to cover optical and electron-optical applications An electron-optical application is described in ref [12] Furthermore, the non-linear parameters will, without exception, be locattons So, the component functions are described by
a~h(x -/3xk, y - / 3 y ~ ) ,
(2)
where a~ is the peak value a n d / 3 k = (/3xk, /3yk) IS the location of the k-th component Also, only second-order, that is, two-component models will be considered to keep the notation simple The analysis for models with one independent variable shows that extension to models of order higher than two IS straightforward The outline of this p a p e r is as follows In the next section, the Raylelgh resolution limit, and limits derived from it, are briefly reviewed In section 3, the incompatibility of the assumptions underlying the Raylelgh criterion with the possibilities of modern m e a s u r e m e n t facilities is discussed The alternative resolution limit is presented m section 4 Its quantification is the subject of section 5 A numerical example is given in section 6 Conclusions are drawn in section 7
2. The classical Raylelgh definition Rayleigh introduces his resolution limit in ref [1] H e considers two equal, superimposed, incoherent, sinc-square-shaped intensity distributions at a distance such that the center of the one d~stribut~on coincides with the first zero of the other Then the combined intensity distribution has two maxima, symmetrically located about a central relative minimum At this distance, the ratio of the value at the minimum to that at the m a m m a is ~ 0 81 Rayleigh states " W e may con-
299
sider this to be about the limit of closeness at which there could be any decided appearance of resolution" Thus defined, the limit clearly refers to visual inspection Later on, the Rayleigh limit has been generallzed to intensity distributions other than sincsquare [2] This generalization implies that the limit for these other distributions is also chosen as the distance from the center to the first zero of the component distribution concerned A difficulty with this choice is that it does not apply if the component distribution has no zero in the neighborhood of its center This is the case with spectroscopic lines as Gauss, Lorentz, and Volgt distributions [13] It appears, therefore, more appropriate to take the generalized limit as the distance for which the ratio of the value at the central minimum to that at the maxima on either side is equal to 0 81 Again other proposals are found m the literature For a summary, see ref [14] For the purposes of this paper, the most important conclusion with respect to these limits is that they are in the first place measures of the geometrical width of the intensity distributions
3. Fundamental difficulties with the Rayleigh definition Under modern conditions, there are difficulties connected with the Raylelgh criterion These difficulties, which are both conceptual and operational, are discussed in this section The Rayleigh resolution limit is, as we saw, based on the presumed capablhtles of the human visual system This leads to conceptual difficulties To see this, suppose that visual inspection is replaced by intensity m e a s u r e m e n t Then the resolution limit is set by the instrumental capablhties to measure differences in Intensity at various points of the combined intensity distribution With a hypothetical, perfect m e a s u r e m e n t system the obvious limit would then be the distance between the component distributions for which both maxlma and the minimum in between just coincide This limit is the so-called Sparrow limit [14] However, in Rayleigh's considerations, the functional form of the component distributions is
300
A van den Bos /Ulttmate resolution a mathemattcalframework
exactly known So, the limit could be further decreased by measuring the combined intensity distribution and numerically fitting a pair of comp o n e n t distributions to it Since Raylelgh does not mention any m e a s u r e m e n t errors, there would be no obvious limit to resolution at all O n e might object that the purpose of the Raylelgh resolution limit is not to establish an ultimate limit, but to c o m p a r e the resolving capablhtles of different instrumental set-ups used for the same purpose However, if the above procedure of intensity m e a s u r e m e n t and subsequent model fitting would be applied to each of the combined intensity distributions p r o d u c e d by the different instruments, the results would be the same unlimited resolution in all cases The reason why these theoretical considerations are not valid in practice is that the assumptions m a d e cannot be met These assumptions have been the following In the first place, it has been assumed that there are no errors m the observations of the c o m b i n e d intensity distribution So, there is no noise (non-systematic error) present in the intensity m e a s u r e m e n t s In the second place, there are no systematic errors in the sense of deviations from the supposed functional form of the c o m p o n e n t functions fitted to the intensity m e a s u r e m e n t s N o r are there systematic contributions to the observations not included in the model fitted, such as trends In prachce, such systematic and non-systematic errors are always present So, if there is a limit to resolution, the cause must be these errors since, as we have seen, there is no obvious hmlt m their absence T h e question then arises how the errors influence the resolvablllty and how resolution can be defined in terms of the errors This will be the subject of sections 4 and 5 of this paper A further difficulty with the Raylelgh criterion is that it is not, and has not been intended to be, operational in the following sense It does not tell us w h e t h e r the available e r r o r - c o r r u p t e d observations can reasonably be split into two c o m p o n e n t distributions or not This is particularly relevant to observations of the c o m b i n e d intensity of comp o n e n t distributions at a distance approximately equal to or smaller than the Sparrow limit U n d e r
these conditions, it would be extremely helpful if we could test the observations with respect to resolvablhty A proposal for such a test is discussed in section 5 of this paper Finally, the Rayleigh criterion is not applicable to c o h e r e n t imaging [14] Raylelgh's incoherent image model consists of the sum of squares of functtons For c o h e r e n t images described by the square of the sum of these funcuons [14], the author has derived results analogous to the ones presented in this paper These results will not be included but will be reported later Anyway, they show that the ideas presented here are not confined to strictly incoherent ~maglng
4. An alternative definItmn Suppose that h(x, y) describes the shape of the c o m p o n e n t distributions as a function of the spatial i n d e p e n d e n t variables x and y Also suppose that the peak value h(0, 0) of h(x, y) is normalized to one For example, h(x, y) may be a point image or a dlffracnon pattern T h e n the c o m b i n e d intensity distribution is described by
g ( x , y, ~, /3) = a l h ( x - f i x , ,
Y-fly,)
+a2h(x-flx2, y -/3y2),
(3)
where a = ( a l a 2 ) T are the amphtudes and /3 = (/3xl/3x2/3yl/3y2) f the peak locattons of the components, respectively In what follows, a reparameterlzed version of (3) will be used
g ( x , y, a, A, fl) = a[Ah(x-/3xl,
y -/3,,)
+(1-A)h(x-/3x2,
y-/3,2)],
(4)
with A = a l / ( a l + a 2) and a = a l + a 2 So, since a l and a 2 will be supposed to be strictly positive, A satisfies 0 < A < 1 Let us next assume that a n u m b e r of N error-disturbed observations w. have been m a d e on g(x, y, a, A, [!) described by
wn = g . ( c ~ , A, /3) + v,,,
n = 1,
,N,
(5)
where gn(a, A,/3) = g(x n, y.,, a, A,/3) and the v,, are the errors Systematic errors are defined as the expectanon E(vn) of the cn Non-systematic
A van den Bos /Ulttmate resolunon a mathemaucalframework errors are defined as v n - E ( v n ) The (x~, y~), n = 1, ,N are the measurement pomts Then the parameters (a, A, /3) of the component functions may be estimated from the w. by fitting the model g(x, y, a, :, b) with respect to (al : ] bY)T to the observatxons m the sense of a suitable criterion of goodness of fit, where a and l are variables corresponding to a and A, respectively, and b is defined as (b~l bx2 by I by2) x Next define the devtauon dn(a, :, b) of the model from the nth observation wn as
d . ( a , :, b) = w . - g . ( a ,
:, b)
(6)
The usual cnterm of goodness of fit are functions of all d~(a, l, b) and are to be minimized with respect to (a] : I bT) T Throughout, the criterion will be taken as the least-squares criterion
J2(a, : , O ) = Y',d2(a, :, b),
n=a,
,N,
n
(7) where the subscript 2 refers to the fact that a two-component, that is, a second-order model is fitted to the observations The choice of the least-squares criterion is hardly restrictive The analysis to follow can also be carried out for other criteria [15] However, the least-squares criterion has the advantage that it produces relatively s~mple expressions Next suppose that the locations (/3xl,/3y 1) and (/3~2, /3y2) are so close that h(x-/3~1, y - / 3 y l) and h ( x - f l x 2 , Y-/3y2) substantially overlap That is, their Euchdean distance A is smaller than their width Also suppose, for the moment, that all v~ are equal to zero Then the combined distribution g(x, y, o~, A, ~ ) becomes increasingly similar to one single-component distribution a h ( x - b~, y - by) with location (b~, by) somewhere in between (/3~1, ~yl ) and (~x2, ~y2 )' and a-~ a~ + o~2 as A decreases Now let this single component be fitted to the observations Then the deviations concerned are described by
d,,(a, b) = w. - a h ( x . - b., y. - by),
(8)
with b = (b x by) r So, the least-squares criterion is described by
Jl(a, b)
=
E d 2 ( a , b), n
n = 1,
,i,
(9)
301
where the subscript 1 refers to the fact that now a smgle-component, that is, a first-order model is fitted Let the least-squares solution for this problem be (4, bx, /~x) Then, under the assumptions made, ([~x, /~y) will be somewhere m between the closely located exact peak locations (/3xl, /3y 1) and (/3x2, /3y2), and d = a I + a 2, as explained above Notice that in the p a r a m e t e r space (a] : I bT) w of J2(a, :, b) the first-order criterion Jl(a, b) is the intersection of the planes bx~ = bx2 and by 1 = by2 with J2(a, :, b) If next the errors are non-zero and at least of the same order of magnitude as the first-order deviations (eq (8)) at (4, bx, /~y) for exact observations, the dlstmgmshlng of the component distributions m the observations may, mtultwely, become very difficult or perhaps impossible Notice that this will happen for smaller and smaller errors as A is smaller So, the question arises what the consequences are for the model-fittmg solutions (4, :, /~) if a two-component model g(x, y, a, : , b) is fitted to the observations (eq (5)) containing errors of this magmtude To answer this question the structure of the criterion J2(a, :, b) will now be analyzed It is observed that m the parameterlzaUon (al : I bX) T of J2(a, :, b), the first-order solution (4, /Sx, /~y) is represented by the points
where t' is an arbitrary scalar The following properties of the points (10) are relevant to the purposes of this p a p e r • Substitution of the points (10) for (a, :, b) m J2(a, :, b) shows that at the points (10) the value of this criterion is, not surprisingly, for all t equal to Jl(a, /~) This Is, by definition, the minim u m criterion value of Jl(a, b) If : is varied, the criterion value J2(a, :, b) at the points (10) remains equal to Jl(fi, b) • For any t', the points (10) are stattonary pomts of J2(a, :, b) with respect to (al bT) T Stationary points of a function are defined as points where the gradient is equal to zero, The statlonanty may easily be proved by first derwmg the expressions for the gradient of J2(a, :, b), defined by eq (7), with respect to (a I bY) v At the points (10), the
302
A van den Bos / Ultzmate resolutton a mathematwal Jramework
resulting elements of the gradient are then seen to be equivalent to the elements of the gradient of J l ( a , b) w~th respect to ( a l b T ) T at (4, /~x, /~) T h e y are, therefore, equal to zero • The points (10) are either saddle points or minima A stationary point ~s a saddle point ff there exist both one or m o r e d~rectlons in which it ~s a m i n i m u m and one or m o r e d~rectlons m which It is a maximum So, for a function of two variables It is a m o u n t a i n pass In the directions m the planes bxl = bx2 and by I = by 2, the point (10) is a m i n i m u m This is because the mtersect10n of these planes with Jz(a, l, b) is the crxten o n Jl(a, b) and (t~, bx, by) is the m l m m u m of the latter criterion The directions not in these planes are b ~ - b x e and by 1 - b y 2 If the point (10) is a maxxmum in any direction in the (b~l b~2, b y ! - b y 2 ) c o o r d i n a t e plane, it is a saddle point If not, it ~s a m l m m u m In section 5 of th~s paper, it will also be shown that for a particular set of observations the points (10) are e~ther all minima or all saddle points, ~rrespect~ve of the value of t' on (0,1) The consequences of the points (10) being either saddle points or minima will now be d~scussed F~rst, let the point (10) be a saddle point T h e n there are d~rect~ons m the (b~ - bx2, by~ by 2) coordinate plane m which the point (10) is a maximum Then, in any case, in the neighborhood and on either s~de of the point (10) there are points w~th b~i - b~2 ~ 0 a n d / o r by~ - by 2 4=- 0 at which the criterion value ~s smaller than at the point (10) More precisely, ~t can be shown that u n d e r these conditions, on e~ther side of the point (10) a m l m m u m is found T h e one miram u m ~s absolute and represents the solution The other ~s relatxve Consequently, there exists a solution for the parameters of the two-component model such that the solutions for the locations (/~1, bye) and (bx2, by2) are d~stlnct In any case, this happens ff the observations are exact T h e n (/~i, /~y~) = (/3x~, /3y~) =g (/3~2, /3y 2) = (/~2, /~y~) O n the other hand, ff the point (10) ~s a m i n i m u m and this m i n i m u m ~s absolute, / ~ = ~x2 = / ~ and ^/~y~ So, the locations (bx~, /)yl) and (bx2, by 2) c o i n c i d e Thus, in sp~te of the fact that the model g(x, y, a, A, /3) underlying the observations is a t w o - c o m p o n e n t model, the
y/~y2=/~y
model estimated from the observations is a onec o m p o n e n t model This may h a p p e n If the errors in the observations are of the order of m a g m t u d e of or larger than the deviations of the best fitting o n e - c o m p o n e n t model from exact observations This coincidence of the solution for/3x~ w~th that for/3x2 and of the solution for/3y~ with that for /3y 2 can be shown to be a change of structure of J2(a, t, b) around the points (10) The errors m the observations have c h a n g e d the two mamma, one saddle-point structure into a simple onem l m m u m structure A precise mathematical explanation of this structural change can be given on the basis of catastrophe theory [16,17] but this ~s outside the scope of this p a p e r O n the bas~s of these considerations we propose the following definition for the resolvablhty of c o m p o n e n t functions from error-d~sturbed observations "Two-component functtons are resoh'-
able from the at'adable observattons if the soluttons for thetr locatton parameters are dtstmct, they are not resolvable tf these soluttons exactly comctde"
5. ObJective quantification of resoivabihty In the prewous section we addressed the possibility of errors transforming the saddle point (10) into a m i n i m u m In th~s section, the errors will be c o m p u t e d which actually do so Generally, at a m i n i m u m all elgenvalues of the matrix of second-order derivatives of a function, the so-called Hessian matrix, are positive At a saddle point, one or more of these elgenvalues are negative Therefore, a saddle point ~s transformed into a minimum at an upward zero-crossmg of the smallest elgenvalue of the Hessian matrix So, all that is r e q m r e d to see when the point (10) becomes a minimum, xs to c o m p u t e all sets of errors v~, ,v u corresponding with such zero-crossings E a c h of these sets ~s represented by one point in the Euclidean N space of the errors All these points together form a hypersurface This hypersurface separates the error sets for whxch the point (10) is a, possibly local, minim u m from those for which It ~s a saddle point The hypersurface may also be shown to be ~dent~cal with or very close to another hypersurface
A van den Bos /Ulttmate resolutton a mathemattcal framework called btfurcatton set [9] The bifurcation set separates the error sets for which the point (10) IS the only, and therefore the absolute, minimum from those for which it is not Thus defined, the bifurcation set separates the error sets for which the least-squares solutions for the locations coincide from those for which the solutions are distinct So, in a simulation experiment, the coinciding solutions for the locations may be seen to spilt into two distinct ones once the errors generated cross the bifurcation set in the direction of the origin In what follows, the hypersurface where the point (10) becomes a minimum will be taken as an approximation of the bifurcation set This approxamatlon is justified by the closeness of both hypersurfaces To slmphfy the notation, the vector t is introduced, defined by t - - (a b~, b.2 by i by2) T
(11)
Next, the elements of the 5 > 5 Hessian matrix H 2 of J2(a, [, b) with respect to t are computed at the stationary point (10) These elements are described by O2J2(t ) (12)
c~t, Otj
at the point (10) with t, j = 1, ,5 Their computation is straightforward and will not be presented here Unfortunately, inspection of the resultmg expressions for the elements (12) does not yield much insight In the eigenvalue structure of H 2 This can be substantially improved by transforming t into
t R = Rt,
(13)
R=
oo 1-t
0
0
0
0
t'
1-
0 0
1 0
-1 0
0 1
-
(14)
So, the new coordinates tR are described by
tR=(a
[b~l+(1-1)b~2
bxl -- bx2
b,1 - bs2) T
In these coordinates, the Hessian m a t r i x /"/2 at the point (10) becomes
R-TH2R -1
(16)
Detailed computations, again left out, then reveal that at the point (10)
R-TH2R -1 = d l a g ( H 1 P ) ,
(17)
where P is a 2 x 2 matrix with elements
Pkm = - - 2 g ( 1 - - t ' ) a ~ d . ( ~ ) n
02h"(bx' /~) ~bkObm , (18)
where b I = b x and b z = by Furthermore, the detailed computations show that in these expressions H 1 is the Hessian matrix of Jl(a, b x, by) with respect to (a b x by) T and d.(?) represents the deviation w n - a hn(b ~, by), both evaluated at (& /~x, /~y) Notice that thus H 2, the Hessian m a m x , of J2(a, {, b) with respect to t at the stationary point (10) is completely determined by the first-order solution (a, bx, /~,) Since this solution is a minimum, the matrix H1 is positive definite at this point This agrees with the considerations of the previous section So, only the elgenvalues of the matrix P in eq (17) need be considered Hence, a necessary condition for the errors to belong to the bifurcation set IS that they make the smallest eIgenvalue of P just vanish That is, for these errors det(P), being the product of the latter elgenvalues, must be equal to zero under the condition that PI~ +P22, being the sum of the elgenvalues, is strictly positive So, this condition is described by d e t ( P ) =PIlP22 -P~2 = 0, with Pll +P22 > 0 (19)
where R is the 5 × 5 matrix lO 0 t
303
gby l + ( 1 - { ) b v 2 (15)
Notice that if this condition is satisfied for any t', it is satisfied for all g on (0,1) It is observed that the bifurcation set is defined by the three equations defining the firstorder solution (&/3~,/;v) and the equation det(P) =ptlp22-p22=O These constitute four equatlons in the N errors in the bifurcation set and tl, ~ and /~v So, after hypothetical ehmlnatlon of tl, b~ and by, one equation in the N errors results This is the equation of the bifurcation set in
304
A van den Bos / Ulnmate resolunon a mathemancal framework
terms of the errors Therefore, xt is a hypersurface, dlwdlng the Euchdean space of the errors into two regions In the one region, including the origin, the solutions for the locations are distract In the other, they coincide In this sense, the bifurcation set constitutes the limit to resoluUon From these considerations the following concluslons may be drawn In the first place, an objective and quanUfiable criterion for resolvablllty has been estabhshed It says that the component functions can only be resolved if the errors in the observations are on the same side of the bifurcation set as the origin The bifurcation set is defined by the class of functions to which the component funcuons underlying the observations belong, by their locations and amplitudes, and by the measurement points In the second place, the proposed resolvabfllty criterion is operational in the following sense Suppose that a set of two-component observations is available Then these observations can be tested with respect to resolvablhty by fitting a one-component model of the same family to them Next, the resulting solution (a, b~, by) is substituted in eq (18) to compute the elements P11, P~2 and P22 If these satisfy
resolve the components This illustrates that modelling errors may have an adverse effect on resolvabxhty Secondly, analogous to the resolvability theory developed for the fitting of functions of one independent variable [8-11], the results of this paper can be extended to include the fitting of models consisting of more than one component and, If needed, systematic contrlbutlons, such as trends
6. Numerical example The purpose of this numerical example is to show the coincidence of the solutions for the locations under the influence of the errors and to demonstrate the accuracy of the resolvablhty criterion (20) For that purpose, In a simulation experiment a number of 100 observations are generated described by
{ ,[ ( x , , - ~ x 2 ) 2 + ( y . - 1 3 ~ 2 ) -,])
+0~ 2 exp - g
(21) Pll > 0
and PIlP2z-p122 > 0,
(20)
the m a m x P ~s posmve defimte and no resolution is possible Otherwise, it is Notice that this test does not require the two-component model to be actually fitted to the observations The purpose is not to determine the distance of the components The purpose is to find out whether the components can be resolved For that the fitting of the one-component model suffices Also notice that the condmon (20) is satisfied for all C on (0,1) whenever it is satisfied for any [ on this interval With respect to th~s test for resolvablhty, two remarks may be made Firstly, assume that the fitted two-component model does not belong to the same family of functions as the component functions underlying the observations Then this may, possibly in combination with the addmve errors vn, also result In poSltlVe definiteness of the matrix P and, hence, in the lmposslbihty to
wlth a 1 = 0 6 , a 2 = 0 4 , ~ x l = 19, ~ x 2 = 2 0 , ~,1 = 19, 13y2 = 18 The measurement points are described by (x n, Y n ) = ( 0 4 p , 0 4 q ) , p, q = 0, ,9 The errors generated satisfy un = K[C 0 + C~(x n - 1 8)],
(22)
with C 0 = - 0 0 0 1 and C x = 0 0 0 0 1 In a n u m b e r of simulations, the error gain K is stepwlse increased from zero on In each step, the two-component model with { = 0 6
-/)
and the one-component model
A uan den Bos /Ulttmate resolutton a mathemattcalframework
20t
305
zero-crossing of d e t ( P ) It is emphasized once m o r e that the c o m p u t a t i o n of d e t ( P ) and P~ll only requtres the o n e - c o m p o n e n t solution (tl, bx, /~y) to be known
l°i[det'p'/ -~0
.o01U h----~
F~g 1 Results of the numerical example of section 6 Soluuons for locahons computed from error-corrupted observations (top) and resolvabd~tycriterion (m~ddle and bottom) as a function of the error gain are fitted to the observations T h e resulting solutions are (fi, bxl, bye, bx2, by2) and (fi, bx, ~ ) , respectwely T h e t w o - c o m p o n e n t solutions bx~, /~yl, b~2 and /~y2 as a function of K are shown m fig 1 (top) To test the c o n d m o n s (20), from the o n e - c o m p o n e n t solution (d, / ~ , / ~ ) the quantltmS Pl~ and det(/z') are c o m p u t e d for the same values of K These are shown m fig 1, b o t t o m and middle, respectwely T h e solutions /~l and /~2 and the solutions by~ and /~y2 for the locations of the t w o - c o m p o n e n t model are seen to s~multaneously coincide m pa~rs for K --- 2 and stay coincident ff K further increases A simple calculation shows that for ~ = 2 the absolutely largest error is equal to - 2 3 6 × 10 -~ This gwes an impression of the m a g m t u d e of the errors revolved T h e solutions for other values of {, not shown m the figure, have been found to coincide for the same value of ~ Th~s agrees with the theory of section 5 which predicts that coincidence does not dep e n d on { F u r t h e r m o r e , the graphs of det(R) and p ~ d e m o n s t r a t e that the criterion (20) accurately indicates when the coincidence of the twoc o m p o n e n t solutions occurs N o h c e that, as expected, th~s coincidence occurs at an upward
7. Conclusions A new deflmtlon of the concept resolution has b e e n p r o p o s e d In the classical Raylelgh defimt~on, the hm~t to resolution of observed component intensity d~strlbutlons is based on the shape of the c o m p o n e n t s and the p r e s u m e d capabdltms of the h u m a n visual system T h e p r o p o s e d deflmtlon is based on the errors in the observations of the mtens~ty d~stnbutlons T h e p r o p o s e d defimt~on ~s objectwe N o assumphons about the hum a n visual system are n e e d e d Which errors make resolution impossible and which do not, only depends on the type of c o m p o n e n t function underlying the observations, the locations and the amphtudes of the components, the m e a s u r e m e n t points, and the type of c o m p o n e n t function fitted to the observations Also, different from the classical definlhon, the p r o p o s e d deflmtlon is apphcable to c o m p o n e n t s having different a m p h t u d e s It is, m addition, operational in the sense that the experimenter can relatively easdy test his observations wxth respect to c o m p o n e n t resolvabdlty
References [1] J W Strutt, Baron Raylelgh, in Scmntlfic Papers by John Wdham Strutt, Baron Raylelgh, Vol 1, 1869-1881 (Cambridge Unwerslty Press, Cambridge, 1899) [2] M Born and E Wolf, Principles of Optics (Pergamon, Oxford, 1980) [3] S G L~pson and H Llpson, Optical Physms (Cambridge Unwerslty Press, Cambridge, 1981) [4] P R Buseck, J M Cowley and L Eyrmg, Hlgh-Resolut~on Transmms~on Electron Microscopy and Assocmted Techmques (Oxford University Press, Oxford, 1988) [5] L Relmer, Scanning Electron Mmroscopy (Springer, Berhn, 1985) [6] A van den Bos, Parameter Estimation, m Handbook of Measurement Scmnce, Ed P Sydenham (Wdey, Chlchester, 1982) p 331 [7] J Norton, An Introduction to Identification (Academic Press, London, 1986) [8] A van den Bos, Automatlca 16 (1980) 487
306
A van den Bos /Ulttmate resolution a mathemattcal framework
[9] A van den Bos, IEE Proc D 128 (1981) 109 [10] A van den Bos, J Opt Soc Am A4 (1987) 1402 [11] A van den Bos, IEEE Trans Autom Control AC-128 (1983) 1118 [12] D van Dyck, to be submitted [13] P A Jansson, Ed, Deconvolution with Apphcatlons in Spectroscopy (Academic Press, Orlando, 1984)
[14] A Papouhs, Systems and Transforms with Applications in Optics (McGraw-Hill, New York, 1968) [15] A van den Bos, Automatica 24 (1988) 803 [16] T Poston and I N Stewart, Catastrophe Theory and its Applications (Pitman, London, 1980) [17] P T Saunders, An Introduction to Catastrophe Theory (Cambridge University Press, Cambridge, 1980)
Ultimate resolution: a mathematical framework A v a n d e n Bos Department of Apphed Physws, Delft Unwerstty of Technology, P 0 Box 5046, 2600 GA Delft, The Netherlands Recewed at Editorial Office 8 January 1992
A new resolution hmlt is proposed Classical hmlts, such as Raylelgh's, are a measure of the width of the component intensity distributions to be resolved The proposed hmlt, on the other hand, is expressed in the errors m the observed intensity distributions It describes how the resolvabdlty of the components depends on these errors It is shown how this new limit can be computed for specified component functions at a given distance In addition, a test for the resolvabdlty of components from a given set of observations is presented
I. Introduction T h e p u r p o s e of this p a p e r is to propose a new definition of optical or e l e c t r o n - o p t l c a l resolution This definition is i n t e n d e d to be objective a n d o p e r a t i o n a l In both respects, it is different from the classical Raylelgh definition [1] widely in use in optical a n d electron-optical l i t e r a t u r e [2-5] Raylelgh's d e f i n i t i o n is based on p r e s u m e d limitations to the resolving capablhtles of the h u m a n visual system Since Raylelgh's days, visual inspection has b e e n s u p p l e m e n t e d with intensity m e a s u r i n g i n s t r u m e n t a t i o n a n d digital c o m p u t i n g facilities T h e resulting resolving capabilities exceed those of the h u m a n visual system This jUStlties the q u e s t i o n as to what the limitations to resolution are u n d e r these new conditions Raylelgh's hmlt to the resolution of two superimposed a n d overlapping c o m p o n e n t Intensity distributions is, to a certain extent, a m e a s u r e of the width of the c o m p o n e n t s No m e n t i o n is m a d e of the p r e s e n c e of errors in the observations of the c o m b i n e d intensity d i s t r i b u t i o n In the definition p r o p o s e d in this paper, these errors are central G i v e n the amplitudes, the locations, the m e a s u r e m e n t points a n d the f u n c t i o n a l description of the c o m p o n e n t distributions, the set of all errors, or equivalently that of the error-dist u r b e d observations, is c o m p u t e d for which the
c o m p o n e n t s are resolvable T h u s the f o r m u l a t e d hmlt to resolution is the b o u n d a r y of this set This definition is also o p e r a t i o n a l since it has b e e n f o u n d to be relatively easy to decide on which side of this b o u n d a r y a given set of observations is located So, o n e can d e t e r m i n e beforeh a n d w h e t h e r or not the c o m p o n e n t s can be resolved from the available observations T h e described error limit to resolution is based on a p a r a m e t e r e s t i m a t i o n approach [6,7] It is supposed that the t w o - c o m p o n e n t f u n c t i o n a l m o d e l is fitted with respect to both locations a n d both a m p l i t u d e s to e r r o r - c o r r u p t e d observations T h e n the resolution of the locations or, equivalently, of the c o m p o n e n t functions achieved is u l t i m a t e since it is b a s e d o n detailed a priori knowledge a b o u t the class of functions conc e r n e d T h u s formulated, the p r o b l e m at h a n d is an example of fitting wetghted-sum models A w e i g h t e d - s u m m o d e l of o r d e r K is described by
g( x , or, fl) = a l h ( X, ill) +
+.Kh( x, 13~), (1)
with flk 4: tim If k ~ m, where a a n d fl are the vectors of the linear p a r a m e t e r s a k a n d the n o n linear p a r a m e t e r s fig, while g ( x , a, fl) a n d h ( x , /3~) are functions of the i n d e p e n d e n t variable x T h e n u m b e r K is called the order of ( l )
0304-3991/92/$05 00 © 1992 - Elsevier Science Pubhshers B V All rights reserved
A van den B o s /
Ulttmate resolution a mathemattcalframework
Weighted-sum models are, for example, used to model multi-exponential decay The parameters a and /3 are then the amplitudes and the decay constants, respectively For weighted-sum models, the computation of the error limit to resolution of the/3 k has been extensively studied [8-11] In this p a p e r the model (1) will be extended to component functions h(x, y, 13k) of two independent variables x and y to cover optical and electron-optical applications An electron-optical application is described in ref [12] Furthermore, the non-linear parameters will, without exception, be locattons So, the component functions are described by
a~h(x -/3xk, y - / 3 y ~ ) ,
(2)
where a~ is the peak value a n d / 3 k = (/3xk, /3yk) IS the location of the k-th component Also, only second-order, that is, two-component models will be considered to keep the notation simple The analysis for models with one independent variable shows that extension to models of order higher than two IS straightforward The outline of this p a p e r is as follows In the next section, the Raylelgh resolution limit, and limits derived from it, are briefly reviewed In section 3, the incompatibility of the assumptions underlying the Raylelgh criterion with the possibilities of modern m e a s u r e m e n t facilities is discussed The alternative resolution limit is presented m section 4 Its quantification is the subject of section 5 A numerical example is given in section 6 Conclusions are drawn in section 7
2. The classical Raylelgh definition Rayleigh introduces his resolution limit in ref [1] H e considers two equal, superimposed, incoherent, sinc-square-shaped intensity distributions at a distance such that the center of the one d~stribut~on coincides with the first zero of the other Then the combined intensity distribution has two maxima, symmetrically located about a central relative minimum At this distance, the ratio of the value at the minimum to that at the m a m m a is ~ 0 81 Rayleigh states " W e may con-
299
sider this to be about the limit of closeness at which there could be any decided appearance of resolution" Thus defined, the limit clearly refers to visual inspection Later on, the Rayleigh limit has been generallzed to intensity distributions other than sincsquare [2] This generalization implies that the limit for these other distributions is also chosen as the distance from the center to the first zero of the component distribution concerned A difficulty with this choice is that it does not apply if the component distribution has no zero in the neighborhood of its center This is the case with spectroscopic lines as Gauss, Lorentz, and Volgt distributions [13] It appears, therefore, more appropriate to take the generalized limit as the distance for which the ratio of the value at the central minimum to that at the maxima on either side is equal to 0 81 Again other proposals are found m the literature For a summary, see ref [14] For the purposes of this paper, the most important conclusion with respect to these limits is that they are in the first place measures of the geometrical width of the intensity distributions
3. Fundamental difficulties with the Rayleigh definition Under modern conditions, there are difficulties connected with the Raylelgh criterion These difficulties, which are both conceptual and operational, are discussed in this section The Rayleigh resolution limit is, as we saw, based on the presumed capablhtles of the human visual system This leads to conceptual difficulties To see this, suppose that visual inspection is replaced by intensity m e a s u r e m e n t Then the resolution limit is set by the instrumental capablhties to measure differences in Intensity at various points of the combined intensity distribution With a hypothetical, perfect m e a s u r e m e n t system the obvious limit would then be the distance between the component distributions for which both maxlma and the minimum in between just coincide This limit is the so-called Sparrow limit [14] However, in Rayleigh's considerations, the functional form of the component distributions is
300
A van den Bos /Ulttmate resolution a mathemattcalframework
exactly known So, the limit could be further decreased by measuring the combined intensity distribution and numerically fitting a pair of comp o n e n t distributions to it Since Raylelgh does not mention any m e a s u r e m e n t errors, there would be no obvious limit to resolution at all O n e might object that the purpose of the Raylelgh resolution limit is not to establish an ultimate limit, but to c o m p a r e the resolving capablhtles of different instrumental set-ups used for the same purpose However, if the above procedure of intensity m e a s u r e m e n t and subsequent model fitting would be applied to each of the combined intensity distributions p r o d u c e d by the different instruments, the results would be the same unlimited resolution in all cases The reason why these theoretical considerations are not valid in practice is that the assumptions m a d e cannot be met These assumptions have been the following In the first place, it has been assumed that there are no errors m the observations of the c o m b i n e d intensity distribution So, there is no noise (non-systematic error) present in the intensity m e a s u r e m e n t s In the second place, there are no systematic errors in the sense of deviations from the supposed functional form of the c o m p o n e n t functions fitted to the intensity m e a s u r e m e n t s N o r are there systematic contributions to the observations not included in the model fitted, such as trends In prachce, such systematic and non-systematic errors are always present So, if there is a limit to resolution, the cause must be these errors since, as we have seen, there is no obvious hmlt m their absence T h e question then arises how the errors influence the resolvablllty and how resolution can be defined in terms of the errors This will be the subject of sections 4 and 5 of this paper A further difficulty with the Raylelgh criterion is that it is not, and has not been intended to be, operational in the following sense It does not tell us w h e t h e r the available e r r o r - c o r r u p t e d observations can reasonably be split into two c o m p o n e n t distributions or not This is particularly relevant to observations of the c o m b i n e d intensity of comp o n e n t distributions at a distance approximately equal to or smaller than the Sparrow limit U n d e r
these conditions, it would be extremely helpful if we could test the observations with respect to resolvablhty A proposal for such a test is discussed in section 5 of this paper Finally, the Rayleigh criterion is not applicable to c o h e r e n t imaging [14] Raylelgh's incoherent image model consists of the sum of squares of functtons For c o h e r e n t images described by the square of the sum of these funcuons [14], the author has derived results analogous to the ones presented in this paper These results will not be included but will be reported later Anyway, they show that the ideas presented here are not confined to strictly incoherent ~maglng
4. An alternative definItmn Suppose that h(x, y) describes the shape of the c o m p o n e n t distributions as a function of the spatial i n d e p e n d e n t variables x and y Also suppose that the peak value h(0, 0) of h(x, y) is normalized to one For example, h(x, y) may be a point image or a dlffracnon pattern T h e n the c o m b i n e d intensity distribution is described by
g ( x , y, ~, /3) = a l h ( x - f i x , ,
Y-fly,)
+a2h(x-flx2, y -/3y2),
(3)
where a = ( a l a 2 ) T are the amphtudes and /3 = (/3xl/3x2/3yl/3y2) f the peak locattons of the components, respectively In what follows, a reparameterlzed version of (3) will be used
g ( x , y, a, A, fl) = a[Ah(x-/3xl,
y -/3,,)
+(1-A)h(x-/3x2,
y-/3,2)],
(4)
with A = a l / ( a l + a 2) and a = a l + a 2 So, since a l and a 2 will be supposed to be strictly positive, A satisfies 0 < A < 1 Let us next assume that a n u m b e r of N error-disturbed observations w. have been m a d e on g(x, y, a, A, [!) described by
wn = g . ( c ~ , A, /3) + v,,,
n = 1,
,N,
(5)
where gn(a, A,/3) = g(x n, y.,, a, A,/3) and the v,, are the errors Systematic errors are defined as the expectanon E(vn) of the cn Non-systematic
A van den Bos /Ulttmate resolunon a mathemaucalframework errors are defined as v n - E ( v n ) The (x~, y~), n = 1, ,N are the measurement pomts Then the parameters (a, A, /3) of the component functions may be estimated from the w. by fitting the model g(x, y, a, :, b) with respect to (al : ] bY)T to the observatxons m the sense of a suitable criterion of goodness of fit, where a and l are variables corresponding to a and A, respectively, and b is defined as (b~l bx2 by I by2) x Next define the devtauon dn(a, :, b) of the model from the nth observation wn as
d . ( a , :, b) = w . - g . ( a ,
:, b)
(6)
The usual cnterm of goodness of fit are functions of all d~(a, l, b) and are to be minimized with respect to (a] : I bT) T Throughout, the criterion will be taken as the least-squares criterion
J2(a, : , O ) = Y',d2(a, :, b),
n=a,
,N,
n
(7) where the subscript 2 refers to the fact that a two-component, that is, a second-order model is fitted to the observations The choice of the least-squares criterion is hardly restrictive The analysis to follow can also be carried out for other criteria [15] However, the least-squares criterion has the advantage that it produces relatively s~mple expressions Next suppose that the locations (/3xl,/3y 1) and (/3~2, /3y2) are so close that h(x-/3~1, y - / 3 y l) and h ( x - f l x 2 , Y-/3y2) substantially overlap That is, their Euchdean distance A is smaller than their width Also suppose, for the moment, that all v~ are equal to zero Then the combined distribution g(x, y, o~, A, ~ ) becomes increasingly similar to one single-component distribution a h ( x - b~, y - by) with location (b~, by) somewhere in between (/3~1, ~yl ) and (~x2, ~y2 )' and a-~ a~ + o~2 as A decreases Now let this single component be fitted to the observations Then the deviations concerned are described by
d,,(a, b) = w. - a h ( x . - b., y. - by),
(8)
with b = (b x by) r So, the least-squares criterion is described by
Jl(a, b)
=
E d 2 ( a , b), n
n = 1,
,i,
(9)
301
where the subscript 1 refers to the fact that now a smgle-component, that is, a first-order model is fitted Let the least-squares solution for this problem be (4, bx, /~x) Then, under the assumptions made, ([~x, /~y) will be somewhere m between the closely located exact peak locations (/3xl, /3y 1) and (/3x2, /3y2), and d = a I + a 2, as explained above Notice that in the p a r a m e t e r space (a] : I bT) w of J2(a, :, b) the first-order criterion Jl(a, b) is the intersection of the planes bx~ = bx2 and by 1 = by2 with J2(a, :, b) If next the errors are non-zero and at least of the same order of magnitude as the first-order deviations (eq (8)) at (4, bx, /~y) for exact observations, the dlstmgmshlng of the component distributions m the observations may, mtultwely, become very difficult or perhaps impossible Notice that this will happen for smaller and smaller errors as A is smaller So, the question arises what the consequences are for the model-fittmg solutions (4, :, /~) if a two-component model g(x, y, a, : , b) is fitted to the observations (eq (5)) containing errors of this magmtude To answer this question the structure of the criterion J2(a, :, b) will now be analyzed It is observed that m the parameterlzaUon (al : I bX) T of J2(a, :, b), the first-order solution (4, /Sx, /~y) is represented by the points
where t' is an arbitrary scalar The following properties of the points (10) are relevant to the purposes of this p a p e r • Substitution of the points (10) for (a, :, b) m J2(a, :, b) shows that at the points (10) the value of this criterion is, not surprisingly, for all t equal to Jl(a, /~) This Is, by definition, the minim u m criterion value of Jl(a, b) If : is varied, the criterion value J2(a, :, b) at the points (10) remains equal to Jl(fi, b) • For any t', the points (10) are stattonary pomts of J2(a, :, b) with respect to (al bT) T Stationary points of a function are defined as points where the gradient is equal to zero, The statlonanty may easily be proved by first derwmg the expressions for the gradient of J2(a, :, b), defined by eq (7), with respect to (a I bY) v At the points (10), the
302
A van den Bos / Ultzmate resolutton a mathematwal Jramework
resulting elements of the gradient are then seen to be equivalent to the elements of the gradient of J l ( a , b) w~th respect to ( a l b T ) T at (4, /~x, /~) T h e y are, therefore, equal to zero • The points (10) are either saddle points or minima A stationary point ~s a saddle point ff there exist both one or m o r e d~rectlons in which it ~s a m i n i m u m and one or m o r e d~rectlons m which It is a maximum So, for a function of two variables It is a m o u n t a i n pass In the directions m the planes bxl = bx2 and by I = by 2, the point (10) is a m i n i m u m This is because the mtersect10n of these planes with Jz(a, l, b) is the crxten o n Jl(a, b) and (t~, bx, by) is the m l m m u m of the latter criterion The directions not in these planes are b ~ - b x e and by 1 - b y 2 If the point (10) is a maxxmum in any direction in the (b~l b~2, b y ! - b y 2 ) c o o r d i n a t e plane, it is a saddle point If not, it ~s a m l m m u m In section 5 of th~s paper, it will also be shown that for a particular set of observations the points (10) are e~ther all minima or all saddle points, ~rrespect~ve of the value of t' on (0,1) The consequences of the points (10) being either saddle points or minima will now be d~scussed F~rst, let the point (10) be a saddle point T h e n there are d~rect~ons m the (b~ - bx2, by~ by 2) coordinate plane m which the point (10) is a maximum Then, in any case, in the neighborhood and on either s~de of the point (10) there are points w~th b~i - b~2 ~ 0 a n d / o r by~ - by 2 4=- 0 at which the criterion value ~s smaller than at the point (10) More precisely, ~t can be shown that u n d e r these conditions, on e~ther side of the point (10) a m l m m u m is found T h e one miram u m ~s absolute and represents the solution The other ~s relatxve Consequently, there exists a solution for the parameters of the two-component model such that the solutions for the locations (/~1, bye) and (bx2, by2) are d~stlnct In any case, this happens ff the observations are exact T h e n (/~i, /~y~) = (/3x~, /3y~) =g (/3~2, /3y 2) = (/~2, /~y~) O n the other hand, ff the point (10) ~s a m i n i m u m and this m i n i m u m ~s absolute, / ~ = ~x2 = / ~ and ^/~y~ So, the locations (bx~, /)yl) and (bx2, by 2) c o i n c i d e Thus, in sp~te of the fact that the model g(x, y, a, A, /3) underlying the observations is a t w o - c o m p o n e n t model, the
y/~y2=/~y
model estimated from the observations is a onec o m p o n e n t model This may h a p p e n If the errors in the observations are of the order of m a g m t u d e of or larger than the deviations of the best fitting o n e - c o m p o n e n t model from exact observations This coincidence of the solution for/3x~ w~th that for/3x2 and of the solution for/3y~ with that for /3y 2 can be shown to be a change of structure of J2(a, t, b) around the points (10) The errors m the observations have c h a n g e d the two mamma, one saddle-point structure into a simple onem l m m u m structure A precise mathematical explanation of this structural change can be given on the basis of catastrophe theory [16,17] but this ~s outside the scope of this p a p e r O n the bas~s of these considerations we propose the following definition for the resolvablhty of c o m p o n e n t functions from error-d~sturbed observations "Two-component functtons are resoh'-
able from the at'adable observattons if the soluttons for thetr locatton parameters are dtstmct, they are not resolvable tf these soluttons exactly comctde"
5. ObJective quantification of resoivabihty In the prewous section we addressed the possibility of errors transforming the saddle point (10) into a m i n i m u m In th~s section, the errors will be c o m p u t e d which actually do so Generally, at a m i n i m u m all elgenvalues of the matrix of second-order derivatives of a function, the so-called Hessian matrix, are positive At a saddle point, one or more of these elgenvalues are negative Therefore, a saddle point ~s transformed into a minimum at an upward zero-crossmg of the smallest elgenvalue of the Hessian matrix So, all that is r e q m r e d to see when the point (10) becomes a minimum, xs to c o m p u t e all sets of errors v~, ,v u corresponding with such zero-crossings E a c h of these sets ~s represented by one point in the Euclidean N space of the errors All these points together form a hypersurface This hypersurface separates the error sets for whxch the point (10) is a, possibly local, minim u m from those for which It ~s a saddle point The hypersurface may also be shown to be ~dent~cal with or very close to another hypersurface
A van den Bos /Ulttmate resolutton a mathemattcal framework called btfurcatton set [9] The bifurcation set separates the error sets for which the point (10) IS the only, and therefore the absolute, minimum from those for which it is not Thus defined, the bifurcation set separates the error sets for which the least-squares solutions for the locations coincide from those for which the solutions are distinct So, in a simulation experiment, the coinciding solutions for the locations may be seen to spilt into two distinct ones once the errors generated cross the bifurcation set in the direction of the origin In what follows, the hypersurface where the point (10) becomes a minimum will be taken as an approximation of the bifurcation set This approxamatlon is justified by the closeness of both hypersurfaces To slmphfy the notation, the vector t is introduced, defined by t - - (a b~, b.2 by i by2) T
(11)
Next, the elements of the 5 > 5 Hessian matrix H 2 of J2(a, [, b) with respect to t are computed at the stationary point (10) These elements are described by O2J2(t ) (12)
c~t, Otj
at the point (10) with t, j = 1, ,5 Their computation is straightforward and will not be presented here Unfortunately, inspection of the resultmg expressions for the elements (12) does not yield much insight In the eigenvalue structure of H 2 This can be substantially improved by transforming t into
t R = Rt,
(13)
R=
oo 1-t
0
0
0
0
t'
1-
0 0
1 0
-1 0
0 1
-
(14)
So, the new coordinates tR are described by
tR=(a
[b~l+(1-1)b~2
bxl -- bx2
b,1 - bs2) T
In these coordinates, the Hessian m a t r i x /"/2 at the point (10) becomes
R-TH2R -1
(16)
Detailed computations, again left out, then reveal that at the point (10)
R-TH2R -1 = d l a g ( H 1 P ) ,
(17)
where P is a 2 x 2 matrix with elements
Pkm = - - 2 g ( 1 - - t ' ) a ~ d . ( ~ ) n
02h"(bx' /~) ~bkObm , (18)
where b I = b x and b z = by Furthermore, the detailed computations show that in these expressions H 1 is the Hessian matrix of Jl(a, b x, by) with respect to (a b x by) T and d.(?) represents the deviation w n - a hn(b ~, by), both evaluated at (& /~x, /~y) Notice that thus H 2, the Hessian m a m x , of J2(a, {, b) with respect to t at the stationary point (10) is completely determined by the first-order solution (a, bx, /~,) Since this solution is a minimum, the matrix H1 is positive definite at this point This agrees with the considerations of the previous section So, only the elgenvalues of the matrix P in eq (17) need be considered Hence, a necessary condition for the errors to belong to the bifurcation set IS that they make the smallest eIgenvalue of P just vanish That is, for these errors det(P), being the product of the latter elgenvalues, must be equal to zero under the condition that PI~ +P22, being the sum of the elgenvalues, is strictly positive So, this condition is described by d e t ( P ) =PIlP22 -P~2 = 0, with Pll +P22 > 0 (19)
where R is the 5 × 5 matrix lO 0 t
303
gby l + ( 1 - { ) b v 2 (15)
Notice that if this condition is satisfied for any t', it is satisfied for all g on (0,1) It is observed that the bifurcation set is defined by the three equations defining the firstorder solution (&/3~,/;v) and the equation det(P) =ptlp22-p22=O These constitute four equatlons in the N errors in the bifurcation set and tl, ~ and /~v So, after hypothetical ehmlnatlon of tl, b~ and by, one equation in the N errors results This is the equation of the bifurcation set in
304
A van den Bos / Ulnmate resolunon a mathemancal framework
terms of the errors Therefore, xt is a hypersurface, dlwdlng the Euchdean space of the errors into two regions In the one region, including the origin, the solutions for the locations are distract In the other, they coincide In this sense, the bifurcation set constitutes the limit to resoluUon From these considerations the following concluslons may be drawn In the first place, an objective and quanUfiable criterion for resolvablllty has been estabhshed It says that the component functions can only be resolved if the errors in the observations are on the same side of the bifurcation set as the origin The bifurcation set is defined by the class of functions to which the component funcuons underlying the observations belong, by their locations and amplitudes, and by the measurement points In the second place, the proposed resolvabfllty criterion is operational in the following sense Suppose that a set of two-component observations is available Then these observations can be tested with respect to resolvablhty by fitting a one-component model of the same family to them Next, the resulting solution (a, b~, by) is substituted in eq (18) to compute the elements P11, P~2 and P22 If these satisfy
resolve the components This illustrates that modelling errors may have an adverse effect on resolvabxhty Secondly, analogous to the resolvability theory developed for the fitting of functions of one independent variable [8-11], the results of this paper can be extended to include the fitting of models consisting of more than one component and, If needed, systematic contrlbutlons, such as trends
6. Numerical example The purpose of this numerical example is to show the coincidence of the solutions for the locations under the influence of the errors and to demonstrate the accuracy of the resolvablhty criterion (20) For that purpose, In a simulation experiment a number of 100 observations are generated described by
{ ,[ ( x , , - ~ x 2 ) 2 + ( y . - 1 3 ~ 2 ) -,])
+0~ 2 exp - g
(21) Pll > 0
and PIlP2z-p122 > 0,
(20)
the m a m x P ~s posmve defimte and no resolution is possible Otherwise, it is Notice that this test does not require the two-component model to be actually fitted to the observations The purpose is not to determine the distance of the components The purpose is to find out whether the components can be resolved For that the fitting of the one-component model suffices Also notice that the condmon (20) is satisfied for all C on (0,1) whenever it is satisfied for any [ on this interval With respect to th~s test for resolvablhty, two remarks may be made Firstly, assume that the fitted two-component model does not belong to the same family of functions as the component functions underlying the observations Then this may, possibly in combination with the addmve errors vn, also result In poSltlVe definiteness of the matrix P and, hence, in the lmposslbihty to
wlth a 1 = 0 6 , a 2 = 0 4 , ~ x l = 19, ~ x 2 = 2 0 , ~,1 = 19, 13y2 = 18 The measurement points are described by (x n, Y n ) = ( 0 4 p , 0 4 q ) , p, q = 0, ,9 The errors generated satisfy un = K[C 0 + C~(x n - 1 8)],
(22)
with C 0 = - 0 0 0 1 and C x = 0 0 0 0 1 In a n u m b e r of simulations, the error gain K is stepwlse increased from zero on In each step, the two-component model with { = 0 6
-/)
and the one-component model
A uan den Bos /Ulttmate resolutton a mathemattcalframework
20t
305
zero-crossing of d e t ( P ) It is emphasized once m o r e that the c o m p u t a t i o n of d e t ( P ) and P~ll only requtres the o n e - c o m p o n e n t solution (tl, bx, /~y) to be known
l°i[det'p'/ -~0
.o01U h----~
F~g 1 Results of the numerical example of section 6 Soluuons for locahons computed from error-corrupted observations (top) and resolvabd~tycriterion (m~ddle and bottom) as a function of the error gain are fitted to the observations T h e resulting solutions are (fi, bxl, bye, bx2, by2) and (fi, bx, ~ ) , respectwely T h e t w o - c o m p o n e n t solutions bx~, /~yl, b~2 and /~y2 as a function of K are shown m fig 1 (top) To test the c o n d m o n s (20), from the o n e - c o m p o n e n t solution (d, / ~ , / ~ ) the quantltmS Pl~ and det(/z') are c o m p u t e d for the same values of K These are shown m fig 1, b o t t o m and middle, respectwely T h e solutions /~l and /~2 and the solutions by~ and /~y2 for the locations of the t w o - c o m p o n e n t model are seen to s~multaneously coincide m pa~rs for K --- 2 and stay coincident ff K further increases A simple calculation shows that for ~ = 2 the absolutely largest error is equal to - 2 3 6 × 10 -~ This gwes an impression of the m a g m t u d e of the errors revolved T h e solutions for other values of {, not shown m the figure, have been found to coincide for the same value of ~ Th~s agrees with the theory of section 5 which predicts that coincidence does not dep e n d on { F u r t h e r m o r e , the graphs of det(R) and p ~ d e m o n s t r a t e that the criterion (20) accurately indicates when the coincidence of the twoc o m p o n e n t solutions occurs N o h c e that, as expected, th~s coincidence occurs at an upward
7. Conclusions A new deflmtlon of the concept resolution has b e e n p r o p o s e d In the classical Raylelgh defimt~on, the hm~t to resolution of observed component intensity d~strlbutlons is based on the shape of the c o m p o n e n t s and the p r e s u m e d capabdltms of the h u m a n visual system T h e p r o p o s e d deflmtlon is based on the errors in the observations of the mtens~ty d~stnbutlons T h e p r o p o s e d defimt~on ~s objectwe N o assumphons about the hum a n visual system are n e e d e d Which errors make resolution impossible and which do not, only depends on the type of c o m p o n e n t function underlying the observations, the locations and the amphtudes of the components, the m e a s u r e m e n t points, and the type of c o m p o n e n t function fitted to the observations Also, different from the classical definlhon, the p r o p o s e d deflmtlon is apphcable to c o m p o n e n t s having different a m p h t u d e s It is, m addition, operational in the sense that the experimenter can relatively easdy test his observations wxth respect to c o m p o n e n t resolvabdlty
References [1] J W Strutt, Baron Raylelgh, in Scmntlfic Papers by John Wdham Strutt, Baron Raylelgh, Vol 1, 1869-1881 (Cambridge Unwerslty Press, Cambridge, 1899) [2] M Born and E Wolf, Principles of Optics (Pergamon, Oxford, 1980) [3] S G L~pson and H Llpson, Optical Physms (Cambridge Unwerslty Press, Cambridge, 1981) [4] P R Buseck, J M Cowley and L Eyrmg, Hlgh-Resolut~on Transmms~on Electron Microscopy and Assocmted Techmques (Oxford University Press, Oxford, 1988) [5] L Relmer, Scanning Electron Mmroscopy (Springer, Berhn, 1985) [6] A van den Bos, Parameter Estimation, m Handbook of Measurement Scmnce, Ed P Sydenham (Wdey, Chlchester, 1982) p 331 [7] J Norton, An Introduction to Identification (Academic Press, London, 1986) [8] A van den Bos, Automatlca 16 (1980) 487
306
A van den Bos /Ulttmate resolution a mathemattcal framework
[9] A van den Bos, IEE Proc D 128 (1981) 109 [10] A van den Bos, J Opt Soc Am A4 (1987) 1402 [11] A van den Bos, IEEE Trans Autom Control AC-128 (1983) 1118 [12] D van Dyck, to be submitted [13] P A Jansson, Ed, Deconvolution with Apphcatlons in Spectroscopy (Academic Press, Orlando, 1984)
[14] A Papouhs, Systems and Transforms with Applications in Optics (McGraw-Hill, New York, 1968) [15] A van den Bos, Automatica 24 (1988) 803 [16] T Poston and I N Stewart, Catastrophe Theory and its Applications (Pitman, London, 1980) [17] P T Saunders, An Introduction to Catastrophe Theory (Cambridge University Press, Cambridge, 1980)