Mathematical image assembly in a scanning-type microscope

Abstract

A method and apparatus for imaging a specimen using a scanning-type microscope, by irradiating a specimen with a beam of radiation using a scanning motion, and detecting a flux of radiation emanating from the specimen in response to the irradiation, in the first sampling session {S.sub.1} of a set {S.sub.n}, gathering data from a first collection of sparsely distributed sampling points {P.sub.1} of set {P.sub.n}. A mathematical registration correction is made to compensate for drift mismatches between different members of the set {P.sub.n}, and an image of the specimen is assembled using the set {P.sub.n} as input to an integrative mathematical reconstruction procedure.

Claims

1. A method of accumulating an image of a specimen using a scanning-type microscope, comprising the following steps: providing a beam of radiation that is directed from a source through an illuminator so as to irradiate the specimen; providing a detector for detecting a flux of radiation emanating from the specimen in response to said irradiation; causing said beam to undergo scanning motion relative to a surface of the specimen, and recording an output of the detector as a function of scan position, in a first sampling session S.sub.1, gathering detector data from a first collection P.sub.1 of sampling points distributed sparsely across the specimen, the collection P.sub.1 comprising fewer than all the sampling points in a sampling grid; repeating this the procedure of gathering detector data from subsequent collections of sampling points so as to accumulate a set {P.sub.n} of such collections, gathered during an associated set {S.sub.n} of sampling sessions, each set with a cardinality N>1; assembling an image of the specimen by using the set {P.sub.n} as input to an integrative mathematical reconstruction procedure, wherein, as part of said assembly process, a mathematical registration correction is made to compensate for drift mismatches between different members of the set {P.sub.n}.

2. A method according to claim 1, wherein: each member P.sub.n of the set {P.sub.n} is used to mathematically reconstruct a corresponding sub-image I.sub.n; said mathematical registration correction is used to align the members of the sub-image set {I.sub.n}; a combined image is mathematically composed from said aligned sub-image set.

3. A method according to claim 1, wherein: prior to reconstruction, said mathematical registration correction is used to align the members of the collection set {P.sub.n}; a composite image is mathematically reconstructed from said aligned collection set.

4. A method according to claim 1, wherein different members of the set {P.sub.n} have different associated sparse distributions of sampling points across the specimen.

5. A method according claim 1, wherein at least one member P.sub.n of the set {P.sub.n} comprises a sparse distribution of sampling points that is not arranged on a regular grid.

6. A method according to claim 1, wherein correction is made for lower-order drift mismatches selected from the group consisting of displacement, rotation, and combinations hereof.

7. A method according to claim 1, wherein correction is made for higher-order drift mismatches selected from the group consisting of skew, shear, scaling, and combinations hereof.

8. A method according to claim 1, wherein the positions of sampling points in at least one given collection P.sub.n are at least partially elected on the basis of an analysis of at least some scan information obtained from at least one previous sampling sessions S.sub.n.

9. A method according to claim 8, wherein; in a given sampling session S.sub.n, sampling points in the associated collection P.sub.n. are visited sequentially while scanning out a line-by-line pattern on the specimen; along a given line L.sub.j in said line-by-line pattern, the positions of sampling points are elected using detection results obtained in scanning a previous line L.sub.i in said line-by-line pattern.

10. A method according to claim 1, wherein, in at least one sampling session S.sub.n, at least some of the sampling points in the associated collection P.sub.n are located below said surface of the specimen.

11. A method according to claim 1, wherein the set {P.sub.n} is accumulated using a plurality of beams of radiation.

12. A method according to claim 1, wherein said radiation comprises charged particles and said microscope comprises a charged-particle microscope.

13. A method according to claim 12, wherein said charged-particle microscope is selected from the group consisting of a Scanning Electron Microscope and a Scanning Transmission Electron Microscope.

14. A method according to claim 1, wherein said radiation comprises photons and said microscope comprises a confocal microscope.

15. A scanning-type microscope, comprising: a specimen holder, for holding a specimen; a source, for producing a beam of radiation; an illuminator, for directing said beam so as to irradiate said specimen; a detector, for detecting a flux of radiation emanating from the specimen in response to said irradiation; beam deflectors, for causing said beam to undergo scanning motion relative to a surface of the specimen; a controller, for recording an output of said detector as a function of scan position, wherein the controller stores instructions which can be invoked to execute the following steps: in a first sampling session S.sub.1, gathering detector data from a first collection P.sub.1 of sampling points distributed sparsely across the specimen, the collection P.sub.1 comprising fewer than all the sampling points in a sampling grid; repeating the procedure of gathering detector data from subsequent collections of sampling points so as to accumulate a set {P.sub.n} of such collections, gathered during an associated set {S.sub.n} of sampling sessions, each set with a cardinality N>1; assembling an image of the specimen by using the set {P.sub.n} as input to an integrative mathematical reconstruction procedure; and as part of said assembly process, making a mathematical registration correction to compensate for drift mismatches between different members of the set {P.sub.n}.

16. The scanning-type microscope of claim 15 wherein the stored instructions include instructions for correction of lower-order drift mismatches selected from the group consisting of displacement, rotation, and combinations hereof.

17. The scanning-type microscope of claim 15 wherein the stored instructions include instructions for correction of higher-order drift mismatches selected from the group consisting of skew, shear, scaling, and combinations hereof.

18. The scanning-type microscope of claim 15, wherein said radiation comprises charged particles and said microscope comprises a charged-particle microscope.

19. The scanning-type microscope of claim 15, wherein said radiation comprises photons and said microscope comprises a confocal microscope.

20. The scanning-type microscope of claim 15 wherein the stored instructions include instructions for, in at least one sampling session S.sub.n, at least some of the sampling points in the associated collection P.sub.n are located below said surface of the specimen.

Description

(1) The invention will now be elucidated in more detail on the basis of exemplary embodiments and the accompanying schematic drawings, in which:

(2) FIG. 1 renders a longitudinal cross-sectional elevation of a scanning-type microscope in which an embodiment of the current invention can be carried out.

(3) FIGS. 2A and 2B schematically depict certain aspects of a conventional method of image accumulation in a scanning-type microscope.

(4) FIGS. 3A and 3B schematically depict certain aspects of an embodiment of a method of image accumulation in a scanning-type microscope according to the current invention.

(5) In the Figures, where pertinent, corresponding parts are indicated using corresponding reference symbols. It should be noted that, in general, the Figures are not to scale.

(6) Embodiment 1

(7) FIG. 1 is a highly schematic depiction of an embodiment of a scanning-type microscope 1 that lends itself to use in conjunction with the current invention; the depicted microscope is a STEM (i.e. a TEM, with scanning functionality) but, in the context of the current invention, it could just as validly be a SEM, confocal microscope, scanning ion microscope, etc. In the Figure, within a vacuum enclosure 2, an electron source 4 (such as a Schottky gun, for example) produces a beam of electrons that traverse an electron-optical illuminator 6, serving to direct/focus them onto a chosen region of a (substantially planar) specimen S. This illuminator 6 has an electron-optical axis 8, and will generally comprise a variety of electrostatic/magnetic lenses, (scan) deflectors, correctors (such as stigmators), etc.; typically, it can also comprise a condenser system.

(8) The specimen S is held on a specimen holder 10 than can be positioned in multiple degrees of freedom by a positioning device (stage) 12; for example, the specimen holder 10 may comprise a finger that can be moved (inter alia) in the XY plane (see the depicted Cartesian coordinate system). Such movement allows different regions of the specimen S to be irradiated/imaged/inspected by the electron beam traveling along axis 8 (in theZ direction) (and/or allows scanning motion to be performed, as an alternative to beam scanning). An optional cooling device 14 is in intimate thermal contact with the specimen holder 10, and is capable of maintaining the latter at cryogenic temperatures, e.g. using a circulating cryogenic coolant to achieve and maintain a desired low temperature.

(9) The focused electron beam traveling along axis 8 will interact with the specimen S in such a manner as to cause various types of stimulated radiation to emanate from the specimen S, including (for example) secondary electrons, backscattered electrons, X-rays and optical radiation (cathodoluminescence). If desired, one or more of these radiation types can be detected with the aid of detector 22, which might be a combined scintillator/photomultiplier or EDX (Energy-Dispersive X-Ray Spectroscopy) detector, for instance; in such a case, an image could be constructed using basically the same principle as in a SEM. However, alternatively or supplementally, one can study electrons that traverse (pass through) the specimen S, emerge from it and continue to propagate (substantially, though generally with some deflection/scattering) along axis 8. Such transmitted electrons enter an imaging system (combined objective/projection lens) 24, which will generally comprise a variety of electrostatic/magnetic lenses, deflectors, correctors (such as stigmators), etc. In normal (non-scanning) TEM mode, this imaging system 24 can focus the transmitted electrons onto a fluorescent screen 26, which, if desired, can be retracted/withdrawn (as schematically indicated by arrows 28) so as to get it out of the way of axis 8. An image of (part of) the specimen S will be formed by imaging system 24 on screen 26, and this may be viewed through viewing port 30 located in a suitable portion of the wall 2. The retraction mechanism for screen 26 may, for example, be mechanical and/or electrical in nature, and is not depicted here.

(10) As an alternative to viewing an image on screen 26, one can instead make use of electron detector D, particularly in STEM mode. To this end, adjuster lens 24 can be enacted so as to shift the focus of the electrons emerging from imaging system 24 and re-direct/focus them onto detector D (rather than the plane of retracted screen 26: see above). At detector D, the electrons can form an image (or diffractogram) that can be processed by controller 50 and displayed on a display device (not depicted), such as a flat panel display, for example. In STEM mode, an output from detector D can be recorded as a function of (X,Y) scanning beam position on the specimen S, and an image can be constructed that is a map of detector output as a function of X,Y. The skilled artisan will be very familiar with these various possibilities, which require no further elucidation here.

(11) Note that the controller (computer processor) 50 is connected to various illustrated components via control lines (buses) 50. This controller 50 can provide a variety of functions, such as synchronizing actions, providing setpoints, processing signals, performing calculations, and displaying messages/information on a display device (not depicted). Needless to say, the (schematically depicted) controller 50 may be (partially) inside or outside the enclosure 2, and may have a unitary or composite structure, as desired. The skilled artisan will understand that the interior of the enclosure 2 does not have to be kept at a strict vacuum; for example, in a so-called Environmental STEM, a background atmosphere of a given gas is deliberately introduced/maintained within the enclosure 2.

(12) When an image of a specimen S is accumulated using a scanning-type microscope such as the subject 1 of FIG. 1, such accumulation occurs on a pixel-by-pixel basis, achieved by scanning the employed imaging beam relative to the specimen S. With reference to the concept of a scan grid G introduced and explained above (see FIG. 2A), such scanning conventionally causes the imaging beam to sequentially observe (and gather imaging data from) every cell C in the scan grid Gi.e. the beam performs 100% observation of the scan grid G. However, in the current invention, a radically different approach is employed, as will now be elucidated in more detail with reference to FIGS. 2 and 3.

(13) Embodiment 2

(14) FIGS. 2A and 2B schematically depict certain aspects of a conventional method of image accumulation in a scanning-type microscope (e.g. of a type as depicted in FIG. 1, or of an alternative type). In this context, FIG. 2A depicts a scan grid G of a type as alluded to above, which is an imaginary mathematical grid/matrix superimposed upon the (XY plane of the) specimen S and containing an array of juxtaposed sampling cells (pixels, sampling points) C; as here depicted, the grid G is orthogonal in nature, though this is not limiting, and other grid geometries (such as polar) could also be conceived. In conventional scanning microscopy, this entire scan grid G is filled because, in tracing out a scan path on the specimen S, the scanning beam sequentially observes (i.e. gathers data from) every cell C in the grid G. If grey shading is used to depict a cell C that is observed (measured) in this manner, then this prior-art situation is represented in FIG. 2B by the fact that the whole grid G is shaded grey. Since there is a certain dwelling time associated with the collection of data from each cell C (e.g. determined by the operating mechanism of detector(s) D and/or 22 in FIG. 1), such a scenario will obviously entail quite a large cumulative (i.e. summed) dwelling time to observe the whole grid G. This cumulative dwelling time will here be denoted by T.sub.G.

(15) Turning now to FIGS. 3A and 3B, these correspond (in broad terms) to FIG. 2B, except in that they depict certain aspects of an embodiment of an alternative, inventive method of image accumulation in a scanning-type microscope. In accordance herewith: (i) Data accumulation (from the specimen S) is now broken down into a number of distinct sampling sessions S.sub.1, S.sub.2, . . . , S.sub.i, . . . , performed in succession. Together, these form a set {S.sub.n} of sampling sessions. (ii) In each such sampling session S.sub.i, instead of visiting all of the cells C in grid G, only a relatively sparse collection P.sub.i, of cells C is visited, representing a relatively thinly-populated subset of all the cells C in G. This situation is schematically illustrated in FIG. 3A for a sampling session S.sub.1, whereby the depicted sprinkling/scattering of gray-shaded cells C represents a first (non-regular) collection P.sub.1, of data cells (pixels, sampling points) that are observed during measurement session S.sub.1. Because P.sub.1, is relatively sparse, the cumulative (summed) dwelling time T.sub.S1, associated with sampling session S.sub.1 will be relatively short. For example, P.sub.1, might have a sparsity (filling factor compared to the entire grid G) of the order of about 2%, whence T.sub.S10.02T.sub.G; this value is not limiting, and greater or smaller sparsities are fully consistent with the current invention (see, in this context, the end of item (v) below).

(16) (iii) FIG. 3B depicts the filling geometry of the grid G for a second (non-regular) collection P.sub.2 of data cells (pixels, sampling points) that are observed during a second measurement session S.sub.2; just as there is a set {S.sub.n} of sampling sessions, there is an associated set {P.sub.n} of sampling point collections, whereby a collection P.sub.i, is gathered during a corresponding session S.sub.i. As in FIG. 3A/item (ii) above, sparse collection P.sub.2 represents a relatively thinly-populated (grey-shaded) subset of all the cells C in G. Once again, because P.sub.2 is relatively sparse, the cumulative (summed) dwelling time T.sub.S2 associated with sampling session S.sub.2 will be relatively short; for example, just as in item (ii) above, P.sub.2 might have a sparsity (filling factor compared to the entire grid G) of the order of about 2%, whence T.sub.S20.02T.sub.G (once again, this value is given purely as an example). In general, in should be noted that, for any given pair of members P.sub.i, P.sub.j, in {P.sub.n}: It can typically be assumed that P.sub.iP.sub.j: in other words, the members of {P.sub.n} are generally mutually disparate (though, in principle, they may have some points/cells/pixels in common). Such mutual disparity can be usefully exploited in the integrative mathematical reconstruction procedure that is used by the current invention to assemble an image from {P.sub.n}, in that, as set forth above, it leads to greater cumulative coverage of the specimen by sampling points (with, as a consequence, a greater confidence level in the final reconstructed image). P.sub.i, P.sub.j do not have to have the same number of data cells (sampling points), i.e. they may have different sparsities.

(17) (iv) According to the invention, the cardinality N (size) of the set {S.sub.n} is a matter of choice, and can be selected in accordance with various factors, such as desired cumulative measurement time and/or imaging sharpness, specimen fragility, etc. In various experiments, the inventors used a whole scala of different values of Nvarying from as little as 2 to as many as 256 (which values are quoted here for purposes of example only, and are not intended to be limiting vis--vis the scope of the appended claims). Depending (inter alia) on the chosen value of N, the cumulative dwelling time T.sub.c=T.sub.Sn (for all N sampling sessions combined) may or may not exceed T.sub.G. For instance: In certain cases, T.sub.c will be less than T.sub.G, thus entailing an immediately evident throughput advantage relative to the prior-art. In other cases, although T.sub.c might not be less than T.sub.G, there will nevertheless be a throughput advantage. This is because, in contemporary scanning microscopy, use if often made of averaging techniques whereby several (e.g. a plurality n.sub.G) full-grid images are obtained and then averaged, so as to mitigate noise effects. In such a scenario, the time required to assemble an averaged image will be n.sub.GT.sub.G; consequently, if T.sub.c<n.sub.GT.sub.G, the invention will offer a throughput advantage relative to such averaging approaches. Even if there is no such throughput advantage, subdividing the imaging procedure into N component sampling sessions in accordance with the invention generally allows sharper imaging results to be obtained than in the prior art (inter alia because of the inter-frame registration correction referred to above and in (vi) below).

(18) (v) Using {P.sub.n} as a basis, an image can be assembled according to the invention using the aforementioned integrative reconstruction procedure. As part of this procedure, the various members of {P.sub.n} will (ultimately) be combined/integrated/hybridized into a composite data set P.sub.c. Depending (inter alia) on choices previously made in steps (i)-(iv), this composite data set P.sub.c may, in principle, have any of a range of possible sparsity values (filling factors compared to 100% coverage of the cells C in grid G). In many instances, P.sub.c will be relatively sparse (e.g. of the order of about 20%), but, despite such sparsity, the invention nevertheless allows a satisfactory image to be mathematically reconstructed. With due regard to points (i)-(iv) above, one can, for example, choose a desired target value for the sparsity of P.sub.c (e.g. 25%), and then correspondingly pick the cardinality N and sparsity of each component collection P.sub.n so as to arrive at this target value (making allowance for possible overlap/redundancy of sampling points within {P.sub.n}).

(19) (vi) As set forth above, the operation in step (v) will have an associated registration correction, which may be performed before, during or after said integrative reconstruction procedure. In this regard, one may, for example, adopt a Type I or Type II approach as discussed above.

(20) More details of the mathematical reconstruction procedure employed by the current invention will be given in the Embodiments that now follow.

(21) Embodiment 3

(22) As already set forth above, the current invention performs a mathematical registration correction to compensate for drift mismatches between different members of the set {P.sub.n}. The general principles of such a registration correction can be elucidated in more detail as follows, whereby the term set will be used to refer to a collection D of data points/pixels acquired for imaging purposes. In particular: When used in the context of a Type I approach (see item (I) above, and item (A) below), the term refers to a reconstructed sub-image I.sub.n. When used in the context of a Type II approach (see item (II) above, and item (B) below), the term refers to a raw collection P.sub.n of sampling points.

(23) One can now distinguish between the following two situations.

(24) (A)

(25) When registering a first set D.sub.1 with a second set D.sub.2, a typical alignment algorithm performs the following tasks:

(26) Regard D.sub.2 as being the result of applying a transformation T to D.sub.1, and start with an initial estimate of this transformation T. Formulate a cost function J (T(D.sub.1), D.sub.2), e.g. as a sum of squared differences, correlation, etc. Solve a local minimization problem

(27) $\underset{T}{\min }J\left(T\left(D_1\right),D_2\right).$

(28) These steps are repeated until convergence occurs, which can, for example, be detected when J no longer decreases substantially. At each step, a pixel-to-pixel comparison is used in the evaluation of the cost function, and J can be typically expressed as:
J(T(D.sub.1), D.sub.2)=(T(D.sub.1)(x,y), D.sub.2 (x,y))dx dy (1a)
where (.,.) is a local set similarity measure (e.g. an I.sup.P norm (.Math..sub.p), a correlation, an inter-pixel mutual information measure, etc.). Because one typically assumes a continuous function for the transformation T (e.g. rotation, scaling, shear, etc.), when T(D.sub.1)(x,y) is evaluated, interpolation can be used to compute an estimate from an original discrete image grid (full regular scan grid G).

(29) (B)

(30) Using the elucidation set forth in (A) above, one can extend the described registration approach to sparse image datasets by comparing a transformed image data point to the nearest one (x.sub.*, y.sub.*) in the target image. This results in the following reformulation of expression (1a):

(31) $\begin{array}{cc}J\left(T\left(D_1\right),D_2\right)=\left(T\left(D_1\right)\left(x,y\right),D_2\left(x_{*},y_{*}\right)\right)xy\mathrm{where}\left(x_{*},y_{*}\right)=\arg {\min }_{u,v}.Math.\left(\begin{array}{c}u\\ v\end{array}\right)-\left(\begin{array}{c}x\\ y\end{array}\right).Math.,& \left(1b\right)\end{array}$

(32) (u,v) set of coordinates of D.sub.2 data points.

(33) If desired, one can limit candidate nearest points to those lying within a certain radius, using an appropriate distance threshold.

(34) It should be noted that such an approach may encompass a point sets registration technique such as the Iterative Closest Point (ICP) algorithm; see, in this context, the following Wikipedia link, for example: en.wikipedia.org/wiki/Iterative_closest_point
Embodiment 4

(35) Some general information on the mathematics of Compressive Sensing (Scanning/Sampling) can, for example, be gleaned from the following references: dsp.rice.edu/cs Emmanuel Cands, Compressive Sampling, Int. Congress of Mathematics 3, pp. 1433-1452, Madrid, Spain, 2006: www-stat.stanford.edu/candes/papers/CompressiveSampling.pdf Richard Baraniuk, Compressive Sensing, IEEE Signal Processing Magazine 24(4), pp. 118-121, Jul. 2007: dsp.rice.edu/sites/dsp.rice.edu/files/cs/baraniukCSlecture07.pdf Justin Romberg, Imaging via Compressive Sampling, IEEE Signal Processing Magazine 25(2), pp. 14-20, Mar. 2008: dsp.rice.edu/sites/dsp.rice.edu/files/cs/Imaging-via-CS.pdf

(36) Essentially, the goal of Compressive Scanning algorithms is the reconstruction of an original signal from compressed measurements thereof. The following elucidation will outline a general approach to such a reconstruction, from which (with the aid of the various references above) the skilled artisan will be able to implement the current invention.

(37) If x .sup.n is -sparse, which is defined as x.sub.0K<<n, one can characterize a sparse acquisition/measurement process by a measurement matrix .sup.mn(m<n).

(38) One can then express the attendant measurements as:
y=x (2)

(39) Literature references show that one can recover the sparse signal x by solving an l.sup.0-minimization problem:

(40) $\begin{array}{cc}\underset{x}{\min }{.Math.x.Math.}_0& \left(3\right)\end{array}$
such that y=x

(41) It has been shown that, if any set of 2 columns from are linearly independent, then the l.sup.0-minimization approach can perfectly recover the original vector x. Despite the fact that an l.sup.0-minimization technique can provide an accurate recovery of x, it is known that, due to the non-convexity of the l.sup.0norm, such reconstruction requires an exhaustive search over all possible combinations, so as to find the sparsest solution. To find a less computationally expensive approach to l.sup.0-minimization, there have been many efforts to develop alternative algorithms. One alternative is to replace an l.sup.0-minimization problem by an l.sup.1-minimization problem:

(42) $\begin{array}{cc}\underset{x}{\min }{.Math.x.Math.}_1& \left(4\right)\end{array}$
such that y=x

(43) If the l.sup.1-norm is assumed to be convex, then solving (4) is computationally feasible. Also, it is known from convex optimization that solving (4) is equivalent to solving the Linear Programming (LP) problem:

(44) $\begin{array}{cc}\underset{t}{\min }{l}^{T}t,& \left(5\right)\end{array}$
subject to txt and y=x where the vector inequality xt means that x.sub.it.sub.i for all i. An advantage of l.sup.1-minimization is the existence of proven numerical solvers. Additionally, this form of minimization has been shown to provide relatively simple conditions guaranteeing the accurate recovery of -sparse signals. These conditions can be formalized as the so-called Restricted Isometry Property (RIP) and the additional Incoherence Property (see mentioned references).

(45) It is worth mentioning that several possible variations on the previously mentioned algorithms take into account various noise models (deterministic noise, stochastic noise, etc.). Furthermore regularization techniques and Bayesian formulations can be used to stabilize convergence and embed prior knowledge.

(46) Despite its advantages, the complexity associated with the LP approach is cubic in the size of the original vector to be recovered (O(n.sup.3)), so that this approach tends to be impractical for large systems. An alternative, more computationally-tractable approach to finding the sparest solution of (2) is based on so-called greedy algorithms. Such algorithms iteratively find an approximation of the original signal and an associated support (defined as the index set of nonzero elements), either by sequentially identifying the support of the signal, or by refining the estimate of the signal gradually.

(47) Representative algorithms of this category include Orthogonal Matching Pursuit (OMP), Iterative Hard Thresholding (I HT), Subspace Pursuit (SP), and Compressive Sampling Matching Pursuit (CoSaMP) algorithms (which are set forth in more detail in the provided references).

(48) In particular, one well-known representative of the greedy approach familyOMP is attractive for its good performance and low computational complexity. The OMP algorithm iteratively estimates the signal x and its support. If the -sparse vector x is supported on T and if we define variables T.sub.k , X.sub.k and r.sub.k as, respectively, the estimated support, the estimated sparse signal, and the residual (r.sub.k=yx.sub.k) in the k-th iteration, then the OMP algorithm repeats the following steps until r.sub.k reaches zero or until a user-defined number of iterations has been reached (assuming initial values k=0, r.sub.0=y, T.sub.0=): Find the largest element in magnitude and the corresponding index t.sub.k among correlations between .sub.i(i-th column of ) and the residual r.sub.k1 generated in the (k1)-th iteration:
t.sub.k=argmax.sub.i|r.sub.k1,.sub.i|(6) Add index t.sub.k into the estimated support set:
T.sub.k=T.sub.k1{t.sub.k}(7) Solve the least squares (LS) problem:
x.sub.k=argmin.sub.supp(u)=T.sub.kyu.sub.2 (8) Update the residual of the k-th iteration as
r.sub.k=yx.sub.k (9)
Embodiment 5

(49) Some further mathematical considerations pertaining to the sparse image registration correction of the current invention will now be elucidated.

(50) Defining a Differentiable and Asymptotically Convex Point-Sets Registration Criterion

(51) As an alternative to the ICP algorithm described earlier, one can use a technique called Gaussian Fields Registration (GFR) to align the sparse image data points (see, for example, references [1], [2] below). This approach defines the registered position as one resulting in the maximum point-to-point overlap (or maximum proximity, in a relaxed form) between reference and transformed datasets.

(52) To derive the GFR criterion, one starts with a basic combinatorial Boolean criterion satisfying the maximum (point-to-point) overlap of two sparse image point-sets:

(53) $M={\left\{P_i\right\}}_{i=1\mathrm{.Math.}{N}_M}\mathrm{and}D={\left\{Q_j\right\}}_{j=1\mathrm{.Math.}{N}_D},$
which are registered using a transformation Tr*. Let us first assume a noiseless case (noise will be addressed later), and also assume that M and D have a maximum point-to-point overlap at the registered position. The ICP algorithm (previously alluded to) was based on this same assumption. Given these definitions, the following criterion (10) will have a global maximum at Tr*:

(54) $\begin{array}{cc}E\left(\mathrm{Tr}\right)=\underset{\underset{j=1\mathrm{.Math.}{N}_D}{i=1\mathrm{.Math.}{N}_M}}{.Math.}\left(d\left(\mathrm{Tr}\left(P_i\right),Q_j\right)\right)\mathrm{with}\left(t\right)=\{\begin{array}{cc}1& \mathrm{for}t=0\\ 0& \mathrm{otherwise}\end{array},& \left(10\right)\end{array}$
where d(P,Q) is a distance measure (e.g. Euclidean) between points. In addition to the sparse point locations, adding a quantity such as the associated image intensity to this criterion is straightforward, and requires just using a higher-dimensional representation of the datasets, where points are defined by both position and a vector of intensity/color attributes:

(55) $M={\left\{\left(P_i,S\left(P_i\right)\right)\right\}}_{i=1\mathrm{.Math.}{N}_M}\mathrm{and}D={\left\{\left(Q_j,S\left(Q_j\right)\right)\right\}}_{j=1\mathrm{.Math.}{N}_D}.$

(56) Given that the combinatorial criterion in (10) is not continuous with respect to the alignment transformations, it will be difficult to find the global maximum. To overcome this problem, one can use a smooth approximation of E(Tr) obtained using an analytical method known as Mollification (see, for example, reference [3] below, in which a similar approach is employed to regularize ill-posed problems with non-differentiable cost functions).

(57) An arbitrary non-differentiable function f(t) defined on .sup.d can be mollified by convolution with the Gaussian kernel

(58) ${}_{}\left(t\right)=\exp \left(\frac{-t^2}{{}^2}\right)$
as follows:

(59) 0$\begin{array}{cc}{f}_{}\left(t\right)=\left({}_{}*f\right)\left(t\right)={}_{}\exp \left(\frac{-{\left(t-s\right)}^2}{{}^2}\right)f\left(s\right)s.& \left(11\right)\end{array}$

(60) The resulting function (t) will satisfy

(61) $\underset{.fwdarw.0}{\lim }{f}_{}\left(t\right)=f\left(t\right)\mathrm{and}{f}_{}{C}^{}\left(\right).$
The transformation described in (11) is also known as the Gauss Transform. If one applies mollification to the criterion E(Tr) [see (10)], one obtains:

(62) $E_{}\left(\mathrm{Tr}\right)=\exp \left(-\frac{{\left(d\left(\mathrm{Tr}\left(P_i\right),Q_j\right)-s\right)}^2}{{}^2}\right)\left\{\underset{\underset{j=1\mathrm{.Math.}{N}_D}{i=1\mathrm{.Math.}{N}_M}}{.Math.}\left(d\left(\mathrm{Tr}\left(P_i\right),Q_j\right)\right)\right\}s$

(63) One can now define:

(64) $d_{\mathrm{ij}}=d\left(\mathrm{Tr}\left(P_i\right),Q_j\right)$$\mathrm{and}$$\begin{array}{c}{E}_{}\left(\mathrm{Tr}\right)=\exp \left(-\frac{{\left(d_{\mathrm{ij}}-s\right)}^2}{{}^2}\right)\left\{\underset{\underset{j=1\mathrm{.Math.}{N}_D}{i=1\mathrm{.Math.}{N}_M}}{.Math.}\left(d_{\mathrm{ij}}\right)\right\}s\\ =\underset{\underset{j=1\mathrm{.Math.}{N}_D}{i=1\mathrm{.Math.}{N}_M}}{.Math.}\exp \left(-\frac{{\left(d_{\mathrm{ij}}-s\right)}^2}{{}^2}\right)\left(d_{\mathrm{ij}}\right)s\\ =\underset{\underset{j=1\mathrm{.Math.}{N}_D}{i=1\mathrm{.Math.}{N}_M}}{.Math.}\exp \left(-\frac{s^2}{{}^2}\right)\left(d_{\mathrm{ij}}-s\right)s.\end{array}$

(65) Knowing that (d.sub.ijs) is non-zero only for s=d.sub.ij, the last integral will be simplified to:

(66) $\exp \left(-\frac{s^2}{{}^2}\right)\left(d_{\mathrm{ij}}-s\right)s=\exp \left(-\frac{d_{\mathrm{ij}}^2}{{}^2}\right)$
which leads to:

(67) $\begin{array}{cc}{E}_{}\left(\mathrm{Tr}\right)=\underset{\underset{j=1\mathrm{.Math.}{N}_D}{i=1\mathrm{.Math.}{N}_M}}{.Math.}\exp \left(-\frac{d_{\mathrm{ij}}^2}{{}^2}\right)=\underset{\underset{j=1\mathrm{.Math.}{N}_D}{i=1\mathrm{.Math.}{N}_M}}{.Math.}\exp \left(-\frac{d^2\left(\mathrm{Tr}\left(P_i\right),Q_j\right)}{{}^2}\right)& \left(12\right)\end{array}$

(68) The mollified criterion E(Tr) is a sum of Gaussians of distances between all pairs of reference and transformed data points. Deriving an analogy from physics, expression (12) can be viewed as the integration of a potential field generated by sources located at points in one of the datasets acting on targets in the other one. The effects of noise, affecting the spatial localization of the point sets, are addressed by relaxing the parameter to values near that of noise variance.

(69) The Gaussian registration criterion can now be extended to include measurement information (e.g. Backscatter intensity, emitted photons intensity . . . ) which is used in addition to the spatial location of the sparse points. This is done by extending the distance measure between points in the criterion as follows:

(70) $\begin{array}{cc}{E}_{,}\left(\mathrm{Tr}\right)=\underset{\underset{j=1\mathrm{.Math.}{N}_D}{i=1\mathrm{.Math.}{N}_M}}{.Math.}\exp \left(-\frac{{.Math.\mathrm{Tr}\left(P_i\right)-Q_j.Math.}^2}{{}^2}-{\left(S\left(\mathrm{Tr}\left(P_i\right)\right)-S\left(Q_j\right)\right)}^T{}^{-1}\left(S\left(\mathrm{Tr}\left(P_i\right)\right)-S\left(Q_j\right)\right)\right))& \left(13\right)\end{array}$
where . . . is the Euclidean distance, and the matrix , which is associated with the measurements vector S(.) is a diagonal matrix with positive elements, which extends the mollification to higher dimensions. Defining:
.sub..sup.ij(Tr)=exp((S(Tr(P.sub.i))S(Q.sub.j)).sup.T.sup.1(S(Tr(P.sub.i))S(Q.sub.j))))
the registration criterion becomes:

(71) $\begin{array}{cc}{E}_{,}\left(\mathrm{Tr}\right)=\underset{\underset{j=1\mathrm{.Math.}{N}_D}{i=1\mathrm{.Math.}{N}_M}}{.Math.}{}_{}^{\mathrm{ij}}\left(\mathrm{Tr}\right)\exp \left(-\frac{{.Math.\mathrm{Tr}\left(P_i\right)-Q_j.Math.}^2}{{}^2}\right)& \left(14\right)\end{array}$

(72) Given that the measurement vector is independent of the aligning transformations, the coefficients .sub..sup.ij will not depend on Tr .

(73) Optimizing the Criterion

(74) For various registration transformations, including rigid and affine models, the criterion E.sub.,(Tr) can be shown to be continuously differentiable. Furthermore, E.sub.,(Tr) will generally have a bell-shape in parameter space in the case of a mixture of closely packed Gaussians. Given this and the nature of the current datasets, one can assume a smooth convex behavior around the registered position. This allows for the use of a variety of powerful convex optimization techniques, such as the quasi-Newton algorithm: see, for example: en.wikipedia.org/wiki/Quasi-Newton_method

(75) The gradient of E.sub.,(Tr) with respect to a given registration parameter is expressed as:

(76) $\begin{array}{cc}\frac{E_{,}\left(\mathrm{Tr}\right)}{}=\underset{\underset{j=1\mathrm{.Math.}{N}_D}{i=1\mathrm{.Math.}{N}_M}}{.Math.}\frac{-2{}_{}^{\mathrm{ij}}}{{}^2}\frac{\mathrm{Tr}\left(P_i\right)}{}.Math.\left(\mathrm{Tr}\left(P_i\right)-Q_j\right)\exp \left(-\frac{{.Math.\mathrm{Tr}\left(P_i\right)-Q_j.Math.}^2}{{}^2}\right)& \left(15\right)\end{array}$

(77) The gradient expression and an approximation of the Hessian are used in the quasi-Newton scheme to update descent directions minimizing E.sub.,(Tr). In each descent direction, a line search routine is used to find the optimum. The procedure is iterated until convergence.

(78) Fast computation of the registration criterion

(79) Evaluating E.sub.,(Tr) at each iteration of the registration algorithm will have a relatively high computational cost of O(N.sub.MN.sub.D). A technique called the Fast Gauss Transform (FGT) (see, for example, references [4], [5] below) can be employed to speed up the process, leading to a computational complexity of only O(N.sub.M+N.sub.D). The FGT method uses the fact that calculations are only needed up to a given accuracy. For computing sums of the form:

(80) $S\left(t_i\right)=\underset{j=1}{\overset{N}{.Math.}}{f}_j\exp \left(-{\left(\frac{s_j-t_i}{}\right)}^2\right),i=1,\mathrm{.Math.},M,$
where T,? are the centers of the Gaussians known as sources and T,? are defined as targets, the following reformulation and expansion in Hermite series is employed:

(81) 0$\begin{array}{cc}\exp \left(\frac{-{\left(t-s\right)}^2}{{}^2}\right)=\exp \left(\frac{-{\left(t-s_0-\left(s-s_0\right)\right)}^2}{{}^2}\right)=\exp \left(\frac{-{\left(t-s_0\right)}^2}{{}^2}\right)\underset{n=0}{\overset{}{.Math.}}\frac{1}{n!}{\left(\frac{s-s_0}{}\right)}^n{H}_n\left(\frac{t-s_0}{}\right)& \left(16\right)\end{array}$

(82) where H.sub.n are the Hermite polynomials. These series converge rapidly and only few terms are needed for a given precision; therefore the new expression can be used to cluster several sources into one virtual source s.sub.0 with a linear cost for a given precision.

(83) The clustered sources can then be evaluated at the targets. In a case where the number of targets is also relatively large, Taylor series (17) can now be used to cluster targets together into a virtual center t.sub.0, further reducing the number of computations

(84) $\begin{array}{cc}\exp \left(-{\left(\frac{t-s}{}\right)}^2\right)=\exp \left(\frac{-{\left(t-t_0-\left(s-t_0\right)\right)}^2}{{}^2}\right)\underset{n=0}{\overset{p}{.Math.}}\frac{1}{n!}{h}_n\left(\frac{s-t_0}{}\right){\left(\frac{t-t_0}{}\right)}^n& \left(17\right)\end{array}$

(85) In (17), the Hermite functions h.sub.n(t) are defined by h.sub.n(t)=e.sup.t.sup.2H.sub.n(t). This can be shown to converge asymptotically to a linear computational complexity as the number of sources and targets increases. For further gains in speed, a variant of the FGT method, called the Improved Fast Gauss Transform (IFGT), can be employed (see, for example, reference [6] below), where a data-clustering scheme along with a multivariate Taylor expansion allows for further computational gains in datasets with high dimensions. In the current convex optimization scheme, computing the gradient of the Gaussian criterion is reduced to the computation of a weighted sum version similar to the criterion itself.

(86) Therefore the gradient can also be evaluated efficiently using FGT techniques.

(87) Some background information relating to certain of the mathematical concepts referred to above can, for example, be gleaned from the following literature sources: [1] F. Boughorbel, A. Koschan, B. Abidi, and M. Abidi, Gaussian Fields: a New Criterion for 3D Rigid Registration, Pattern Recognition 37 (7), pp. 1567-1571 (Jul. 2004). [2] F. Boughorbel, M. Mercimek, A. Koschan, and M. Abidi, A new method for the registration of three-dimensional point-sets: The Gaussian Fields Framework, Image and Vision Computing 28, pp. 124-137 (2010). [3] D. A. Murio, The Mollification Method and the Numerical Solution of III-Posed Problems, Wiley, N.Y. (1993). [4] A. Elgammal, R. Duraiswami, L. Davis, Efficient kernel density estimation using the Fast Gauss Transform with applications to color modeling and tracking, IEEE Trans. Pattern Analysis and Machine Intelligence 25 (11), pp. 1499-1504 (2003). [5] Greengard, J. Strain, The fast Gauss Transform, SIAM J. Scientific Computing 12 (1), pp. 79-94 (1991). [6] C. Yang, R. Duraiswami, N. A. Gumerov, L. Davis, Improved fast gauss transform and efficient kernel density estimation, Proc. Ninth Int. Conf. Computer Vision, IEEE, Nice, France, pp. 464-471 (2003).