1. Introduction

Scanning laser ophthalmoscopy (SLO) is a powerful tool for imaging the retina of the human eye. Images are acquired raster-wise at rates typically in the range of tens of frames per second. Since the eye is usually moving during image acquisition, these images are usually displaced randomly with respect to each other, requiring image registration. Registration offers the ability to analyse the retinal motion, to build a mosaic image from different retinal images, to increase signal-to-noise from several images and to achieve video stabilization. However, since the images are acquired as a raster scan they are often distorted by fixational eye movements that occur during the acquisition. These distortions, which can be severe with saccades, make the registration of sets of such images non-trivial. (One approach to study the eye movements relies on measuring these distortions using blood vessels [1]: vertically oriented blood vessels are required to measure horizontal saccades and vice-versa.)

A number of methods have been proposed to register retinal images and the approach taken depends on both the acquisition method and to a large extent on the angular size of the image. For large fields of view, it is important to consider the field curvature arising from the geometrical properties of the retina, The key point becomes then to properly model the retinal geometry (see [2]). In contrast, over small fields of view (1 or 2 degrees), the distortions resulting from the geometrical properties of the retina are negligible and the image registration problem is dominated by fixational eye movements.

Despite current knowledge about the eye movement during image acquisition [3], there is no consensus for an appropriate model for the deformation field to map small field-of-view SLO images and a variety of alternate approaches have been proposed. Patch-based cross-correlation methods are often used to remove translational motion of SLO images with a reduced field of view [4,5]. More generally, rotational motion has been modelled using the map-seeking circuit algorithm [6, 7]. Recently, the use of a second-order polynomial to model the field that maps two SLO images has been proposed [8], with the transformation estimated by tracking stable point features in the images. A similar procedure was previously used by Nourrit et al. but with the points tracked manually [9]. Nourrit et al. estimated transformation models that included a range of global and local mapping strategies.

Our approach here is to devise a deformation field using a few limited assumptions about the geometry and motion of the eye and the principles of the imaging acquisition system. In general the eye may rotate around three axes (horizontal, vertical, and torsional) during the image acquisition, which combined with the non-planar geometry of the retina results in a wide range of image deformations that the registration procedure must deal with. However, we can make a number of simplifying assumptions:

The deformations can then be described completely by rigid translation of each line in the image.

Note that methods using cross correlation techniques [4, 5] are implicitly based on these assumptions but applied to a strip comprising several lines of the image rather than to an individual line. However, these rigid translation assumptions are not always valid for a strip, especially in the case of a saccade, and this leads to failure of the registration process during such motion. Failure to cope with saccades is particularly problematic if the purpose of the study is to analyse eye motion. On the other hand, applying cross-correlation to single lines is complicated both by the poorer signal-to-noise that results from fewer samples in the correlation and by the difficulty estimating displacements in the direction perpendicular to the line.

Our proposed solution is to model the displacement vector field using hierarchical basis functions. This approach is of great interest in the image registration community [10–12] since it enables estimation of the transformation at successively finer scales making the approach more robust. Here we propose modelling the displacement vector field using a multiscale parametric representation (the B-spline functions as suggested in [11]), while constraining each line of the reference image to be mapped by a pure translation.

2. Method

2.1. Problem formalization

The goal of registration is to find a geometrical mapping h that maps a floating image I_float onto a reference image I_ref such that the two images are aligned. In the ideal case, the mapping is such that any point on the retina shares the same coordinates on both the reference image and the registered image.

Since the mapping is continuous, sample points in the floating image do not necessarily correspond to samples in the reference image and interpolation is required to evaluate the floating image on the grid defined by the samples on the reference image. The mapping could be implemented either from floating image coordinates into reference image coordinates or vice-versa. The second of these (reference to floating) is generally preferred since it permits standard regular sample interpolation. In this case interpolation of the floating image is made using the regular sample grid of the floating image itself, rather than using the corresponding irregularly-spaced samples that this grid maps onto in the reference coordinate system. The registered image I_reg at coordinate vector s is finally obtained as

2.2. Properties of the mapping between the reference image and the floating image

We write the mapping from the reference to the floating image coordinates as

Under the aforementioned assumption that the deformation is described by translations of the lines, u(s) is constant for a horizontal line and we can write

 

Fig. 1 The mapping h and the associated displacement vector field u. A horizontal line in the reference is mapped by rigid translation to a shifted horizontal line in the floating image coordinates.

Download Full Size | PDF

2.3. Modelling of the displacement vector field

Each of the x and y displacements u_x(.), and u_y(.) can be represented with a one-dimensional B-spline of order 1. An advantage of using a low-order B-spline over other representations such as polynomial, harmonic functions, radial basis functions or wavelets is that they have a small overlap. This reduces the interdependency between the parameters which makes the minimization problem easier to solve [11]. Additionally, B-splines have been shown to have very good approximation properties [11].

A B-spline of order 1 defines an approximation of a function using a triangle function basis as illustrated in Fig. 2. At a given scale l, the function is represented as a weighted sum of 2^l+1 – 1 triangular basis functions ${\varphi }_i^l$. These basis functions are obtained by translation and scaling of the triangle function ϕ(t) = {1 – 2|t| for |t| < 1; 0 elsewhere}. Then, for the spline representation of the deformation field at scale l we have

 

Fig. 2 Modelling with one dimensional B-splines of order 1: (a) the triangle basis function ${\varphi }_i^l(.)$, (b) the three basis functions used at scale l = 1, and (c) an example of a function which can be modelled at scale l = 1.

Download Full Size | PDF

Note that we have defined the functions $u_x^l(.)$ and $u_y^l(.)$ on a support [−1, 1] that is twice as wide as the vertical support of the reference image, which is taken to be Ω_x = [−1/2, 1/2]. Using a larger support than the region of interest is necessary to allow the represented functions to have non-zero values at the boundary of Ω_x. However, this also means that at higher scales, some of the spline elements fall completely outside Ω_x and so have no influence on the registration mapping.

2.4. Optimization

The optimization cost function uses a similarity criterion to quantify how well the images match. Here we use a measure of the inter-image difference between the reference image and the mapped floating image for the similarity criterion,

Non-rigid registration is known to be an ill-posed problem and the following regularization term can be added to the cost function to favour smoothness [10],

Care needs to be taken to prevent the optimization procedure from becoming trapped in local minima of the cost function. We use a hierarchical optimization approach, in which optimization starts at the coarse scales, and then iterates through successively finer scales. We begin by estimating an initial translation between the images (which is added to the displacement vector for all subsequent optimization) before estimating the displacement vector using B-splines. We start at scale l = 0, and then move successively through finer scales until l = l_max. The solution obtained at each scale is used as the starting point for optimizing at the next scale (note that the scalable property of the splines means that the displacement vector at each scale l can be represented exactly at the next scale l + 1). Note also that the final resulting displacement vector, including the initial translation can be completely represented in the final B-spline field We also use a hierarchical representation of the images, starting with reduced resolution images of I_ref and of I_float at the lowest scales and increasing the resolution throughout the algorithm.

At each scale l, the optimization parameters are the set of spline weights for the displacement vector. Joint estimation of parameters is non-trivial, thus we use a sequential optimization procedure. We optimize the two parameters associated to one basis function at a time, keeping the rest constant. The result of a sequential optimization procedure in general can depend on the order of visitation of the parameters. However, the finite support of the spline basis functions allows us to split the basis functions into two sets (one comprising the odd indexed, the other the even indexed parameters). Within each set, the basis functions are completely disjoint, therefore the visitation order does not matter. Our procedure is to iterate alternately between the odd and the even sets of basis functions. The detailed procedure is described in Algorithm (1).

Algorithm 1. Blockwise optimization procedure

View Table | View all tables in this article

It is convenient to ensure monotonicity in the vertical transformation (which preserves the scan line ordering) by requiring that whenever the deformation field parameters $a_{x;i}^l$ are modified, they satisfy

This approach would be sufficient in most cases, but the optimization can get caught in local minima when the correct solution is far from the initialization (mainly at coarse scales in the presence of large eye movements or saccades). Therefore, to avoid such local minima, after optimization at a given scale a multi-start algorithm is employed on the region with the poorest fit (identified by the local cost function over the domain of each basis function). If this results in modification of the deformation parameters, further optimization in the adjacent, overlapping regions is carried out using Algorithm (1). The complete optimization procedure at scale l is finally described in Algorithm (2).

2.5. Desinusoiding and Registration

During the acquisition of a line, the speed of the scanner may vary. This is the case for our images which were acquired using a resonant scanner with a sinusoidal horizontal scan speed. Desinusoiding is a key step before registration. Indeed, if the reference and the floating images are not properly desinusoided, the assumption according to which lines of the reference image need only a rigid translation to match those of the floating image is no longer valid. It is possible to estimate or measure the desinusoiding transformation h_des and its inverse ${\mathbf{h}}_{\text{des}}^{-1}$ for the instrument. Then the reference image can be mapped to the floating image with the transformation ${\mathbf{h}}_{\text{des}}\circ \mathbf{h}\circ {\mathbf{h}}_{\text{des}}^{-1}$, where h is again the registration transformation that satisfies the rigid translation properties we require. To avoid any extra interpolations (which can introduce artefacts) the desinusoiding transformations are incorporated directly so that I_float is warped according to ${\mathbf{h}}_{\text{des}}\circ \mathbf{h}\circ {\mathbf{h}}_{\text{des}}^{-1}$ in the optimization method. At the end, the reference image and the registered image can be warped into the same (desinusoided) coordinate space: the desinusoided reference image can be computed as I_ref ○ h_des and the desinusoided registered image can be computed as I_float ○ h_des ○ h.

Algorithm 2. Whole optimization procedure at scale l

View Table | View all tables in this article

3. Results

We applied the registration to 26 SLO video sequences each comprising 32 images of size 256 × 256 pixels, which corresponded to 1 and 2 degree fields on the retina using a modified Heidelberg HRT. The horizontal scanning motion is made at a round trip frequency of 4kHz, hence the time to acquire a single image is 32 ms. The frame rate is approximately 40Hz. The acquisition time of each line is approximately 60 μs.

The registration method was applied down to a scale corresponding to one line per basis function. At the finest scale the number of parameters to describe the deformation field was 256×2 (vertical and horizontal components). For each of the image sequences registration was performed both with and without a regularization term (equivalent to setting λ = 0 in Eq. (8)). Since there was no need to favour smoothness at coarser scales, λ was always set to 0 at the beginning of the registration for scales with functions with domain larger or equal than 16 pixels.

3.1. Video stabilization

We applied the registration without regularization to stabilization of video sequences, each containing 32 images. For each sequence, all images were registered onto a reference image (which was chosen manually).

Figure 3 shows frames from three representative video sequences both before and after registration. The video sequences can be seen online. Note that the cone photoreceptor mosaic can be clearly observed in all images, except for those from the fovea (c). In all unregistered sequences two clear motions can be observed: a relatively slow two-dimensional translational drift throughout the sequence, and vertical compression and expansion on different portions of the images. This is very clear on the sequence of the fovea. Both effects originate from ocular movements. The first sequence (a and b) also contains a large microsaccadic motion, which is also apparent in the distortion of the black regions around the registered image (b). In the registered sequences, these motions are no longer observable. The photoreceptor mosaic becomes stable, and the features in the retinal image do not move. The foveal image is also registered well, even if the photoreceptor mosaic could not be resolved by the optics of the instrument. Note that these sequences also show some distortion at the top edge of the images. This is in fact due to the non-uniform motion of the vertical scanner. Although not necessary for the registration, a desinusoiding procedure could be used in the vertical direction to remove such scanner distortion from the images. (This is in contrast to horizontal scan distortions which are relevant to our method of registration and for which reason we apply horizontal desinusoiding).

 

Fig. 3 Sample images from three SLO video sequences of the retina (all 1 degree field): ∼1.5 deg eccentricity (a) raw sequence ( Media 1) and (d) stabilized sequence ( Media 2); centre of the fovea (c) raw ( Media 3) and (d) stabilized ( Media 4); and ∼ 2 deg eccentricity (e) raw ( Media 5) and (f) stabilized ( Media 6). These correspond respectively to sequences 1, 14 and 25 in Fig. 5 (Note the black regions in the registered sequences correspond to regions of the reference image that were not contained in the floating image due to the eye motion).

Download Full Size | PDF

Moreover, the registration reveals features in the sequences that were very difficult to observe in the raw sequences. It is possible to observe small dynamic motion in small localized parts of the stabilized video. Observe for example the regions marked on Fig. 3(f) in the corresponding video sequence. In regions A, and B, the marked photoreceptors appear to vascillate slightly from left to right (A) and up and down (B). Note that these motions are highly localized. The registration process performed rigid translations of complete horizontal lines, and given that the rest of the line does not move, it is unlikely that these local movements are artefacts of the registration process. Moreover, the motions typically last across several frames of the sequence, making it unlikely that they are artefacts of the acquisition or a result of eye tremors. This suggests that the motion in the images could in fact due to dynamic changes in the retina itself, perhaps due to cardio-vasculature motion, however this is not clear. Dynamic blood flow is clearly evident, particularly in the movie of Fig. 3(f) in the regions highlighted with a white curved line. In region C, we can observe that the framed photoreceptor fades out smoothly throughout the sequence. It seems unlikely that this is caused by focus variations as the rest of the image around it stays sharp and in focus. It may be that the length or the refractive index of that individual photoreceptor is changing rapidly, which would have an interference effect on the light that comes back from the eye, as suggested in [14].

Figure 4 shows temporal averages of the registered sequences of Fig. 3. The individual photoreceptors can be distinguished more clearly in the temporal average Fig. 4(b) than in the reference Fig. 4(a). The improvement in signal-to-noise in the temporal averages is further evidence of the registration accuracy. In sequences where noise was not so prevalent (Fig. 4(c) and 4(e)), the corresponding temporal averages (Fig. 4(d) and 4(f)) show no loss of detail such as would occur if the images were not accurately registered.

 

Fig. 4 Enlarged sections showing temporal averages of the registered sequences from Fig. 3. Left: (a), (c) and (e) show the reference frames from the sequences in Fig. 3. Right: (b), (d) and (f) show the temporal averages of the registered sequences.

Download Full Size | PDF

Figure 5 plots the standard deviation of the irradiance in the temporal averages as a measure of contrast for each of the 26 video sequences. In all cases the contrast for the registered sequences is comparable to that in the reference image. Although this measure provides a general indication of registration, caution should be excercised: noise averaging will reduce contrast in the temporal average even with good registration and the approach will not show if there are small areas or single frames that are poorly registered. The registered sequences were also assessed visually. Good registration was achieved over the whole field for every frame in all 26 sequences.

 

Fig. 5 Standard deviation of irradiance in the temporal average normalized to the standard deviation in the reference image for the 26 registered video sequences. Results are shown both for sequences registered with (blue) and without (red) regularization and for the unregistered sequences (black) for comparison.

Download Full Size | PDF

We observed that results obtained without regularization were very similar to those with regularization, the only difference being the presence of slightly larger high-frequency components in the non-regularized results (see Fig. 6 for example plots of the deformation fields both with and without regularization). This is perhaps unsurprising since the number of free parameters (2 × 256) is small compared with the number of image pixels (256²). Regularization becomes necessary when there is insufficient data in the image to guide the registration and, for example, non-rigid registration can easily become an ill-posed problem where homogeneous regions in an image have to be mapped. However, this is unlikely to happen when registering complete lines of retinal images which have fine structure everywhere in the image.

 

Fig. 6 Displacement fields during the microsaccade motion of sequence Fig. 3 (b) calculated both with and without regularization and the difference between the two. (Note, all units are in pixels and the fields are plotted only for the valid region where the reference and registered images overlap.)

Download Full Size | PDF

3.2. Microsaccade motion

Several of the image sequences exhibited significant microsaccade motion. The registration process was nevertheless successful in all the cases tested. Figure 7 shows an example of registration (with regularization) of images with such motion. A two-frame movie ( Media 7) showing the reference and registered images of this example is also provided so that the reader can check that they are well-aligned. Note the severe distortion on the floating image and the associated distortion in the displacement field (u_x and u_y). In the microsaccade near the bottom of I_float (extending approximately from line 160 to line 230) the horizontal component (u_y) is particularly severe, with a horizontal shift of approximately 35 pixels over 70 lines - i.e., half a pixel shift per line scan. The black regions near the edges of the registered image, which correspond to regions that were outside the floating image, give a further indication of the severity of the motion that is nevertheless correctly registered here. This example corresponds to a microsaccade of approximate magnitude 16° s⁻¹ (assuming the reference image is fixed). Martinez-Conde et al. refer to studies that indicate microsaccades up to 120° s⁻¹ [3]. Although we did not observe any in our sequences, we note that our assumption of neglibible fixational eye motion during a line-scan would not be valid with such movements.

 

Fig. 7 Image registration with microsaccade motion. Floating, reference, and the registered floating image I_float ○ h are shown. The plots indicate the x and y components of the registration displacement field in pixel units for the overlap region. Media 7 shows alternately the reference and registered images to illustrate the alignment.

Download Full Size | PDF

There is also a strong vertical component to the displacement field near the top of the image, extending from line 10 to line 70 (see u_x in Fig. 7). In this case, eye motion is not completely responsible for the deformation: the vertical scan speed of this instrument is not constant near the top of the frame and small vertical motion of the eye can give larger displacement. The registration method has corrected for this motion as well as the fixation motion, illustrating the robustness of the approach. Although not necessary for the registration, a desinusoiding procedure could be used in the vertical direction to remove such scanner distortion from the images.

4. Conclusions and discussion

We have introduced a new and robust automated registration method for sequences of SLO images. The method represents the mapping between two images using a deformation field modelled with a multi-scale parametric B-spline function. Limited assumptions are made about the motion of the eye and the SLO acquisition: it is assumed that the scanning is fast enough to freeze the motion of the eye during the acquisition of a single line. This reduces the required mapping from a 2-D mapping between pixels to a 1-D mapping between scan-lines in the image. However, apart from limits on the maximum gradient of the deformation field (for convenience rather than strict necessity), no other restrictions are placed on its possible form. An hierarchical optimization method is then used to find the deformation field parameters. The process is entirely automated, apart from the initial selection of an image from the sequence to act as the reference. In our implementation, each registration took between 20 and 40 seconds using a single core of a 2.93GHz Intel Core 2 Duo CPU E7500 (depending on the complexity of the required deformation field). Our requirement was for a high fidelity post-process registration and we have not attempted to apply parallel programming techniques. However, we note that the use of B-splines should enable significant vectorization of the algorithm if required.

The registration algorithm was applied to 26 different video sequences of retinal images. In every case it was found to register all images in the sequence successfully. The stabilization of the video sequences revealed some examples of highly localized motion of photoreceptors and other features which were not readily apparent in the raw sequences. Some of these may be related to dynamic blood flow in small capillaries, but this requires further investigation.

Of particular note is that the registration was successful even in those sequences that exhibited significant distortion from microsaccades. Some of the examples shown included motion comparable in magnitude to the vertical scan speed. This is a capability that, to our knowledge, has not previously been reported and it suggests that the method could also prove particularly useful in studies of eye motion itself. Saccadic motion often leads to failures in other registration methods (such as patch-based correlation) and the corresponding data is lost. This is especially problematic if the purpose of the study is to analyse eye motion. In contrast, the relatively unrestricted nature of the deformation field used here means it is capable of representing such motion. Moreover, the high fidelity of the resultant stabilization indicates that the deformation field does indeed accurately capture the distortion introduced by these fixational eye movements. (Of course, estimating the motion of the eye from the displacement vector fields does not follow directly since the reference image itself is acquired while the eye is moving. However, a procedure as described in [6] can be used here).

We found that for many of the images, reasonable results could be obtained if the optimization were stopped before reaching the finest scale. The B-spline model is quite general and has not been specifically tailored for the expected image motion. By considering the constraints of the eye motion and its mechanics, it should be possible to derive a more compact, optimized model requiring fewer parameters. However, an advantage of the general B-spline approach is that the deformation fields are not influenced or limited by such prior knowledge (correct or otherwise). This is likely to be important if the purpose is to study the ocular motion itself.