OSA's Digital Library

Optics Express

Optics Express

  • Editor: C. Martijn de Sterke
  • Vol. 19, Iss. 5 — Feb. 28, 2011
  • pp: 4157–4169
« Show journal navigation

Scanning ophthalmoscope retinal image registration using one-dimensional deformation fields

S. Faisan, D. Lara, and C. Paterson  »View Author Affiliations


Optics Express, Vol. 19, Issue 5, pp. 4157-4169 (2011)
http://dx.doi.org/10.1364/OE.19.004157


View Full Text Article

Acrobat PDF (1297 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

We present a new, robust and automated method for registering sequences of images acquired from scanning ophthalmoscopes. The method uses a multi-scale B-spline representation of the deformation field to map images to each other and an hierarchical optimization method. We applied the method to video sequences acquired from different parts of the retina. In all cases, the registration was successful, even in the presence of large distortions from microsaccades, and the resulting deformation fields describe the fixational motion of the eye. The registration reveals examples of dynamic photoreceptor behaviour in the sequences.

© 2011 Optical Society of America

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

1. M. Stetter, R. A. Sendtner, and G. T. Timberlake, “A novel method for measuring saccade profiles using the scanning laser ophthalmoscope,” Vis. Res. 36, 1987–1994 (1996). [PubMed]

]: 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

2. A. Can, C. V. Stewart, B. Roysam, and H. L. Tannenbaum, “A Feature-Based, Robust, Hierarchical Algorithm for Registering Pairs of Images of the Curved Human Retina,” IEEE Trans. Pattern Anal. Mach. Intell. 24(3), 347–364 (2002)

]). 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

3. S. Martinez-Conde, S. L. Macknik, and D. H. Hube, “The role of fixational eye movements in visual perception,” Nat. Rev. Neurosci. 5, 229–240 (2004). [PubMed]

], 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

4. J. B. Mulligan, “Recovery of motion parameters from distortions in scanned images,” Proceedings of the NASA Image Registration Workshop (IRW97), NASA Goddard Space Flight Center, MD, 1997.

,5

5. S. B. Stevenson and A. Roorda, “Correcting for miniature eye movements in high resolution scanning laser ophthalmoscopy,” in Ophthalmic Technologies XV , edited by Fabrice Manns, Per Soderberg, and Arthur Ho, Proc. of SPIE Vol. 5688A, pp. 145–151 (2005).

]. More generally, rotational motion has been modelled using the map-seeking circuit algorithm [6

6. C. R. Vogel, D. W. Arathorn, A. Roorda, and A. Parker, “Retinal motion estimation in adaptive optics scanning laser ophthalmoscopy,” Opt. Express 14(2), 487–497 (2006). [PubMed]

, 7

7. D. W. Arathorn, Q. Yang, C. R. Vogel, Y. Zhang, P. Tiruveedhula, and A. Roorda, “Retinally stabilized cone-targeted stimulus delivery,” Opt. Express 15(21), 13731–13744 (2007) [PubMed]

]. Recently, the use of a second-order polynomial to model the field that maps two SLO images has been proposed [8

8. H. Li, J. Lu, G. Shi, and Y. Zhang, “Tracking features in retinal images of adaptive optics confocal scanning laser ophthalmoscope using KLT-SIFT algorithm,” Biomed. Opt. Express 1, 31–40 (2010).

], 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

9. V. Nourrit, J. M. Bueno, B. Vohnsen, and P. Artal, “Nonlinear registration for scanned retinal images: application to ocular polarimetry,” Appl. Opt. 47(29), 5341–5347 (2008). [PubMed]

]. 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:
  • For a small field of view, the non-planar geometry of the retina can be neglected. Then rotation of the eye about the horizontal and vertical axes can be treated as a rigid translation of the image.
  • We shall assume that the fixational eye motion during a single horizontal line scan can be neglected. A single horizontal scan is captured in tens of microseconds, and even at the peak velocity of a microsaccade, the line shift in this time is of the order of a fraction of a pixel. A similar assumption is made in [5

    5. S. B. Stevenson and A. Roorda, “Correcting for miniature eye movements in high resolution scanning laser ophthalmoscopy,” in Ophthalmic Technologies XV , edited by Fabrice Manns, Per Soderberg, and Arthur Ho, Proc. of SPIE Vol. 5688A, pp. 145–151 (2005).

    ].
  • We shall assume that rotation about the torsional axis can be neglected. The rotation around the torsional axis, which is known as ocular torsion, may give rise to complicated image distortions. However, these ocular torsions are neglected in most papers, leading nevertheless to good results.

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

Note that methods using cross correlation techniques [4

4. J. B. Mulligan, “Recovery of motion parameters from distortions in scanned images,” Proceedings of the NASA Image Registration Workshop (IRW97), NASA Goddard Space Flight Center, MD, 1997.

, 5

5. S. B. Stevenson and A. Roorda, “Correcting for miniature eye movements in high resolution scanning laser ophthalmoscopy,” in Ophthalmic Technologies XV , edited by Fabrice Manns, Per Soderberg, and Arthur Ho, Proc. of SPIE Vol. 5688A, pp. 145–151 (2005).

] 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

10. V. Noblet, C. Heinrich, F. Heitz, and J.-P. Armspach, “3-D deformable image registration: a topology preservation scheme based on hierarchical deformation models and interval analysis optimization,” IEEE Trans. Image Process. 14, 553–566 (2005). [PubMed]

12

12. H. Lester and S. R. Arridge, “A survey of hierarchical non-linear medical image registration,” Pattern Recogn. 32(1), 129–149 (1999).

] 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

11. J. Kybic and M. Unser, “Fast parametric elastic image registration,” IEEE Trans. Image Process. 12, 1427–1442 (2003).

]), 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 Ifloat onto a reference image Iref 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 Ireg at coordinate vector s is finally obtained as
Ireg(s)=Ifloath(s)Ifloat(h(s)).
(1)
The registration of the floating image Ifloat onto the reference image Iref can be stated as the following constrained optimization problem:
θopt=argminθΘE(Iref,Ifloathθ),
(2)
where θ represents the parameters that define the transformation hθ, Θ is the set of admissible parameters and E(.,.) is the similarity criterion, or cost function.

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
h(s)=s+u(s),
(3)
where s is a point on the reference image and u is the displacement vector field (representing the deformation). Note that if u is the null function, h corresponds to the identity transformation.

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
u(s)=u(sx,sy)=u(sx)=[ux(sx),uy(sx)]T,
(4)
where ux(.), and uy(.) are translations in the vertical (x), and horizontal (y) axes (see Fig. 1). Notice that these two functions depend only on sx: this means that the displacement field u is a function of one dimension.

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.

2.3. Modelling of the displacement vector field

Each of the x and y displacements ux(.), and uy(.) 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

11. J. Kybic and M. Unser, “Fast parametric elastic image registration,” IEEE Trans. Image Process. 12, 1427–1442 (2003).

]. Additionally, B-splines have been shown to have very good approximation properties [11

11. J. Kybic and M. Unser, “Fast parametric elastic image registration,” IEEE Trans. Image Process. 12, 1427–1442 (2003).

].

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 2l+1 – 1 triangular basis functions ϕil. 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
ul(sx)=(uxl(sx)uyl(sx))=(iax;ilϕil(sx)iay;ilϕil(sx)),withϕil(x)=2l/2ϕ(2lxi)
(5)
where al={ax,il,ay,il}i=(2l1)(2l1) are the spline weight parameters for the displacement vector.

Fig. 2 Modelling with one dimensional B-splines of order 1: (a) the triangle basis function ϕil(.), (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.

Note that we have defined the functions uxl(.) and uyl(.) 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,
E(sim)(h)=1N(h)1sΩd;h(s)Ωf(Iref(s)Ifloath(s)),
(6)
where Iref and Ifloat are the reference and floating image (with gray-levels linearly scaled to have the same mean and variance). The summation is over those image pixels in the reference image (the discrete space Ωd) that map onto the domain of the floating image (Ω). N(h) is the number of terms in the summation, which depends on the image mapping h. Note that whereas a normalized cross-correlation criterion would use a modulus-squared, here we use the absolute value for f which is less sensitive to outliers.

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

10. V. Noblet, C. Heinrich, F. Heitz, and J.-P. Armspach, “3-D deformable image registration: a topology preservation scheme based on hierarchical deformation models and interval analysis optimization,” IEEE Trans. Image Process. 14, 553–566 (2005). [PubMed]

],
E(reg)(h)=Ωx(ux(x)x)2+(uy(x)x)2dx.
(7)
The resulting minimization cost function is
E=(1λ)E(sim)+λCE(reg),
(8)
where λ is a weighting factor of the regularization term and C a scaling factor to make E(sim) and E(reg) comparable. In fact here we present results both using with and without a regularisation term in the cost function.

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

table-icon
View This Table
| View All Tables

It is convenient to ensure monotonicity in the vertical transformation (which preserves the scan line ordering) by requiring that whenever the deformation field parameters ax;il are modified, they satisfy
ax;ilax;i1l>123l/2andax;ilax;i+1l<123l/2.
(9)
In practice, the parameters are optimized initially without constraints using a simplex algorithm [13

13. V. Chvátal, Linear programming, (Freeman, 1983).

]. If the resulting parameters do not satisfy the conditions in Eq. (9), which is very unusual, a constrained optimization procedure is then used. Note that modifying the parameters ax,il and ay,il affects only the values of the functions uxl or uyl inside the domain of the corresponding spline element ϕil(.), which avoids costly re-evaluation of the whole cost function when changing these two parameters.

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 hdes and its inverse hdes1 for the instrument. Then the reference image can be mapped to the floating image with the transformation hdeshhdes1, 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 Ifloat is warped according to hdeshhdes1 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 Irefhdes and the desinusoided registered image can be computed as Ifloathdesh.

Algorithm 2. Whole optimization procedure at scale l

table-icon
View This Table
| View All Tables

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).

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

14. R. S. Jonnal, J. Rha, Y. Zhang, B. Cense, W. Gao, and D. T. Miller, “In vivo functional imaging of humancone photoreceptors,” Opt. Express 15(24), 16141–16160 (2007). [PubMed]

].

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.

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.

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 (2562). 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.)

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 (ux and uy). In the microsaccade near the bottom of Ifloat (extending approximately from line 160 to line 230) the horizontal component (uy) 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−1 (assuming the reference image is fixed). Martinez-Conde et al. refer to studies that indicate microsaccades up to 120° s−1 [3

3. S. Martinez-Conde, S. L. Macknik, and D. H. Hube, “The role of fixational eye movements in visual perception,” Nat. Rev. Neurosci. 5, 229–240 (2004). [PubMed]

]. 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 Ifloath 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.

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 ux 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.

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

6. C. R. Vogel, D. W. Arathorn, A. Roorda, and A. Parker, “Retinal motion estimation in adaptive optics scanning laser ophthalmoscopy,” Opt. Express 14(2), 487–497 (2006). [PubMed]

] 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.

Acknowledgments

The authors would like to acknowledge Ibolya Kepiro for providing some of the raw image sequences used in this study. This work was carried out with financial support of the Wellcome Trust.

References and links

1.

M. Stetter, R. A. Sendtner, and G. T. Timberlake, “A novel method for measuring saccade profiles using the scanning laser ophthalmoscope,” Vis. Res. 36, 1987–1994 (1996). [PubMed]

2.

A. Can, C. V. Stewart, B. Roysam, and H. L. Tannenbaum, “A Feature-Based, Robust, Hierarchical Algorithm for Registering Pairs of Images of the Curved Human Retina,” IEEE Trans. Pattern Anal. Mach. Intell. 24(3), 347–364 (2002)

3.

S. Martinez-Conde, S. L. Macknik, and D. H. Hube, “The role of fixational eye movements in visual perception,” Nat. Rev. Neurosci. 5, 229–240 (2004). [PubMed]

4.

J. B. Mulligan, “Recovery of motion parameters from distortions in scanned images,” Proceedings of the NASA Image Registration Workshop (IRW97), NASA Goddard Space Flight Center, MD, 1997.

5.

S. B. Stevenson and A. Roorda, “Correcting for miniature eye movements in high resolution scanning laser ophthalmoscopy,” in Ophthalmic Technologies XV , edited by Fabrice Manns, Per Soderberg, and Arthur Ho, Proc. of SPIE Vol. 5688A, pp. 145–151 (2005).

6.

C. R. Vogel, D. W. Arathorn, A. Roorda, and A. Parker, “Retinal motion estimation in adaptive optics scanning laser ophthalmoscopy,” Opt. Express 14(2), 487–497 (2006). [PubMed]

7.

D. W. Arathorn, Q. Yang, C. R. Vogel, Y. Zhang, P. Tiruveedhula, and A. Roorda, “Retinally stabilized cone-targeted stimulus delivery,” Opt. Express 15(21), 13731–13744 (2007) [PubMed]

8.

H. Li, J. Lu, G. Shi, and Y. Zhang, “Tracking features in retinal images of adaptive optics confocal scanning laser ophthalmoscope using KLT-SIFT algorithm,” Biomed. Opt. Express 1, 31–40 (2010).

9.

V. Nourrit, J. M. Bueno, B. Vohnsen, and P. Artal, “Nonlinear registration for scanned retinal images: application to ocular polarimetry,” Appl. Opt. 47(29), 5341–5347 (2008). [PubMed]

10.

V. Noblet, C. Heinrich, F. Heitz, and J.-P. Armspach, “3-D deformable image registration: a topology preservation scheme based on hierarchical deformation models and interval analysis optimization,” IEEE Trans. Image Process. 14, 553–566 (2005). [PubMed]

11.

J. Kybic and M. Unser, “Fast parametric elastic image registration,” IEEE Trans. Image Process. 12, 1427–1442 (2003).

12.

H. Lester and S. R. Arridge, “A survey of hierarchical non-linear medical image registration,” Pattern Recogn. 32(1), 129–149 (1999).

13.

V. Chvátal, Linear programming, (Freeman, 1983).

14.

R. S. Jonnal, J. Rha, Y. Zhang, B. Cense, W. Gao, and D. T. Miller, “In vivo functional imaging of humancone photoreceptors,” Opt. Express 15(24), 16141–16160 (2007). [PubMed]

OCIS Codes
(100.0100) Image processing : Image processing
(170.4470) Medical optics and biotechnology : Ophthalmology
(170.5810) Medical optics and biotechnology : Scanning microscopy
(170.5755) Medical optics and biotechnology : Retina scanning

ToC Category:
Medical Optics and Biotechnology

History
Original Manuscript: October 8, 2010
Revised Manuscript: December 17, 2010
Manuscript Accepted: December 19, 2010
Published: February 17, 2011

Virtual Issues
Vol. 6, Iss. 3 Virtual Journal for Biomedical Optics

Citation
S. Faisan, D. Lara, and C. Paterson, "Scanning ophthalmoscope retinal image registration using one-dimensional deformation fields," Opt. Express 19, 4157-4169 (2011)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-19-5-4157


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. M. Stetter, R. A. Sendtner, and G. T. Timberlake, "A novel method for measuring saccade profiles using the scanning laser ophthalmoscope," Vision Res. 36, 1987-1994 (1996). [PubMed]
  2. A. Can, C. V. Stewart, B. Roysam, and H. L. Tannenbaum, "A Feature-Based, Robust, Hierarchical Algorithm for Registering Pairs of Images of the Curved Human Retina," IEEE Trans. Pattern Anal. Mach. Intell. 24(3), 347-364 (2002).
  3. S. Martinez-Conde, S. L. Macknik, and D. H. Hube, "The role of fixational eye movements in visual perception," Nat. Rev. Neurosci. 5, 229-240 (2004). [PubMed]
  4. J. B. Mulligan, "Recovery of motion parameters from distortions in scanned images," Proceedings of the NASA Image Registration Workshop (IRW97), NASA Goddard Space Flight Center, MD, 1997.
  5. S. B. Stevenson, and A. Roorda, "Correcting for miniature eye movements in high resolution scanning laser ophthalmoscopy," in Ophthalmic Technologies XV, edited by Fabrice Manns, Per Soderberg, Arthur Ho, Proc. of SPIE Vol. 5688A, pp. 145-151 (2005).
  6. C. R. Vogel, D. W. Arathorn, A. Roorda, and A. Parker, "Retinal motion estimation in adaptive optics scanning laser ophthalmoscopy," Opt. Express 14(2), 487-497 (2006). [PubMed]
  7. D. W. Arathorn, Q. Yang, C. R. Vogel, Y. Zhang, P. Tiruveedhula, and A. Roorda, "Retinally stabilized conetargeted stimulus delivery," Opt. Express 15(21), 13731-13744 (2007). [PubMed]
  8. H. Li, J. Lu, G. Shi, and Y. Zhang, "Tracking features in retinal images of adaptive optics confocal scanning laser ophthalmoscope using KLT-SIFT algorithm," Biomed. Opt. Express 1, 31-40 (2010).
  9. V. Nourrit, J. M. Bueno, B. Vohnsen, and P. Artal, "Nonlinear registration for scanned retinal images: application to ocular polarimetry," Appl. Opt. 47(29), 5341-5347 (2008). [PubMed]
  10. V. Noblet, C. Heinrich, F. Heitz, and J.-P. Armspach, "3-D deformable image registration: a topology preservation scheme based on hierarchical deformation models and interval analysis optimization," IEEE Trans. Image Process. 14, 553-566 (2005). [PubMed]
  11. J. Kybic, and M. Unser, "Fast parametric elastic image registration," IEEE Trans. Image Process. 12, 1427-1442 (2003).
  12. H. Lester, and S. R. Arridge, "A survey of hierarchical non-linear medical image registration," Pattern Recognit. 32(1), 129-149 (1999).
  13. V. Chvátal, Linear programming, (Freeman, 1983).
  14. R. S. Jonnal, J. Rha, Y. Zhang, B. Cense, W. Gao, and D. T. Miller, "In vivo functional imaging of human cone photoreceptors," Opt. Express 15(24), 16141-16160 (2007). [PubMed]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Supplementary Material


» Media 1: AVI (654 KB)     
» Media 2: AVI (619 KB)     
» Media 3: AVI (654 KB)     
» Media 4: AVI (631 KB)     
» Media 5: AVI (654 KB)     
» Media 6: AVI (625 KB)     
» Media 7: AVI (47 KB)     

« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited