1. Introduction

The advent of the green fluorescence protein (GFP) transgenic mouse model with the glial fibrillary acidic protein (GFAP) promoter has significantly advanced in vivo retinal imaging of GFP labeled astrocytes [1]. It provides valuable insight to the dynamic changes in astrocyte morphology during development and in response to physiological and pathological conditions.

The images of these fluorescently labeled cells in the optic disc of the mouse are captured via a scanning laser ophthalmoscope (SLO) where the mouse retina is scanned with a laser beam from a point source and the fluorescence emitted is then collected by a photodetector. However, the in vivo imaging of the mouse retina does not achieve diffraction limited resolution due to the wavefront aberration induced by the eye.

Despite the higher system complexity and cost incurred, the use of adaptive optics in mouse retinal imaging has reduced higher order aberration [2, 3]. However, the poor SNR level of these images, brought about by the thermal and shot noise from the detector and its amplifier, is another important factor which adversely affects image resolution. This is due to the application of image averaging on the sequence of successively acquired retinal images where the noise level may be reduced but at the expense of losing image resolution. The respiratory motion of the mouse results in significant changes to the morphology and location of the astrocytes between successive retinal images such that the averaged composite image appears diffused and occasionally contains undesirable image artifacts.

Various denoising methods are available [4, 5, 6, 7] but, for lack of a better method and also due to its effective noise reducing capability [8], image averaging has become the de facto technique for filtering noise from image sequences. Several modifications to the image averaging method have however been suggested such as the alignment of images prior to averaging [9, 10, 11, 12, 13, 14, 15, 16, 17] or alternatively the rejection of some images which are significantly different from that of a predefined reference image [18, 19, 20, 21, 22].

Image alignment requires the identification of landmark points such as the bifurcation of vessel branches and the establishment of point to point correspondence between landmark points in successive images via automated techniques [15, 17, 23]. However, unlike fluorescence angiography images, the high noise levels and the poor image resolution in our endogenous GFAP–GFP fluorescence images makes it a challenging task to even define these landmark points manually let alone in an automated manner. Even the alignment of these images is error prone since the astrocytes in our images exhibit a far greater degree of freedom compared to the vessels in the fluorescence angiography images which can be easily aligned using a rigid body transformation. More powerful alignment methods are therefore required but ironically this requires a larger number of landmark points which we cannot possibly identify from our images [12]. Conversely, alignment methods which do not rely on landmark points have not been convincing [19]. Here, we propose a novel frame averaging method based on a pixel rank matching criterion which we will henceforth refer to as “averaging via rank matching” (ARM). The two main advantages of ARM are that it achieves significantly higher SNR and cellular resolution compared to conventional frame averaging methods. This is attributed to the noise robustness of the rank matching criterion and the fact that it leverages on the a-priori estimates of the underlying cellular profile obtained by deconvolving the individual retinal images [24] from the sequence. Section 2.4.4 describes this in detail. This is a significant finding as it enables for the first time, the in vivo analysis on the role of astrocytes in the pathogenesis of various diseases such as diabetic retinopathy, Alzheimer’s disease and Parkinson’s disease. In addition, ARM can also be easily applied in the averaging of images acquired from other modalities such as MRI and confocal microscopy.

2. Materials and methods

2.1. Transgenic GFAP–GFP mice

The generation and genotyping of the transgenic GFAP–GFP mice are as previously described by Zhuo et al. [1]. We used adult mice (8–10 weeks old) in the FVB/N background. Animal husbandry was provided by the Biological Resource Centre in Biopolis, Singapore whereas the experimental protocol covering the current study was approved by the Institutional Animal Care and Use Committee (IACUC).

2.2. Preparation of transgenic mice for retinal imaging

Mice were anaesthetized by intra-peritoneal (i.p.) injections with 0.15ml/10g body weight of Avertin (1.5% 2,2,2-tribromoethanol; T48402) purchased from Sigma-Aldrich (St. Louis, MO, USA), and their pupils dilated with a drop of 0.5% Cyclogyl® sterile ophthalmic solution (cyclopentolate hydrochloride, Alcon®, Puurs, Belgium). Custom-made PMMA hard contact lenses (from Cantor & Nissel, Northamptonshire, UK) were used to avoid dehydration of the cornea. Careful eye examination prior to scanning laser ophthalmoscope imaging ruled out the presence of any corneal or lens opacities.

2.3. Scanning laser ophthalmoscope (SLO) imaging

We employed a commercially available scanning laser ophthalmoscope, the Heidelberg Retina Angiograph 2, HRA 2 (Heidelberg Engineering, Dossenheim, Germany), for retinal imaging of the mice. The imaging system was adapted for the optics of the mouse eye by replacing the 30° focal lens with a 55° wide angle objective. This reduces the laser beam diameter to 1.7mm and in turn allows more light to be coupled into the small mouse pupil. Maximum retinal irradiance is approximately 2.0mWcm^-2 and, therefore, lies below the limits established by the American National Standards Institute. Since our primary interest is in the endogenous GFP signal in the mouse retina, we operate the HRA 2 in the fluorescence mode where the 488nm Argon laser provides the excitation light. The emission light is collected, with the use of the 500nm barrier filter, by a photodetector. The collected signal is then fed to a frame grabber interfaced to a computer at a frame rate of 5Hz where each image sequence comprises 45 single image frames. Each image is defined by 1536×1536 pixels with an intensity resolution of 8 bits/pixel. The estimated size of a focused spot on the retina i.e. resolution limit is approximately 10µm. By Nyquist sampling criterion, the sampling density must be at least 5.0µm/pixel i.e. twice the bandwidth frequency. Since we assume that 1° of visual angle subtended at the mouse eye corresponds to about 31µm of its retina [25] and our field of view is about 55.9°, the actual sampling density works out to be approximately $\frac{31\times 55.9}{1536}\mu m/\mathrm{pixel}=1.13\mu m/\mathrm{pixel}$ .

2.4. Proposed frame averaging algorithm

Here, we aim to obtain a composite fluorescence image of the mouse retina, with improved SNR and cellular resolution, from a sequence of relatively noisy and poorly resolved retinal fluorescence images. This is achieved by first aligning the images such that the between-frame motion artifacts are reduced and then averaging the aligned images. To this end, the proposed ARM algorithm comprises four major sections: (i) Detecting landmark points, (ii) Establishing point-to-point correspondence, (iii) Aligning images via affine transformation and (iv) Aligning images via the novel rank matching technique.

2.4.1. Detecting landmark points

The high noise levels in the image sequence prevent us from accurately detecting the landmark points. The image sequence can be formally defined as I_i, i=1,2, …,N where the subscript i identifies a particular image in the sequence of N images and N=45. In ARM, we propose an automated method of selecting landmark points from the a-priori estimates J_i, i=1,2, …,N of the underlying fluorescence signal in each retinal image from the original sequence. The a-priori estimates are obtained by deconvolving each image I_i (Huygens Essential Software, Scientific Volume Imaging BV) [24] in the sequence. The deconvolution software computes a theoretical point spread function (PSF), resembling a 2-D Gaussian profile in the x-y plane, based on the ophthalmoscope device settings. Also, it automatically estimates the image SNR and uses this as a regularization parameter to control the sharpness of the deconvolved image. In doing so, photon noise is effectively removed. Next, we apply a gray-scale morphological operation [26] to locate the intensity maximas (peak intensities) in each J_i estimate. Naturally, the locations of the intensity maximas in the I_i images will be identical to those in the corresponding J_i images. These intensity maximas serve as landmark points where their displacement over successive images in the sequence can be used to align these images. Their localized “spot-like” neighborhoods in J_i makes them effective landmark points since displacements along any direction, in the x-y plane, can be accurately detected.

2.4.2. Establishing point-to-point correspondence between adjacent images

Classical frame averaging methods establish correspondence between the landmark points in the images I_i, from a sequence, to those in the reference image I ₁ from the same sequence. However, establishing this correspondence can be impossible in some cases since the fluorescence signal in image I_i of the sequence may have changed significantly in terms of intensity and location from the corresponding signal in the reference image I ₁. ARM establishes point-to-point correspondence between successive images I _i-1 and I_i in the sequence such that each pair of successive images in the sequence has its own set of corresponding landmark points. To this end, the normalized cross-correlation measure (NCC) [23]

is applied as shown in Eq. 1 where J_i,u,v is the (2W_c+1)×(2W_c+1) neighborhood of the i^th deconvolved image centered around one of its landmark points (u,v) and W_c is the width of the neighborhood, J̄_i,u,v is the mean value of this neighborhood, J _{i-1,u′,v′} is the (2W_c+1)×(2W_c+1) neighborhood of the (i-1)^th deconvolved image centered around one of its landmark points (u′,v′) and J̄ _{i-1,u′,v′} is the mean value of this neighborhood. The notation 〈·〉 denotes the dot operator whereas |·| denotes the magnitude operator. We again use the a-priori estimates J_i since they are cleaner and better resolved compared to the I_i images.

Point-to-point correspondence between landmark points (u,v) and (u′,v′) is established if (u,v) lies within the (2W_c+1)×(2W_c+1) neighborhood of (u′,v′) and their NCC measure, ${\mathrm{NCC}}_{J_{i,u,v},J_{i-1,u\prime ,v\prime }}$ , is the largest exceeding an empirically predefined threshold, NCC_th. NCC_th may vary from one sequence to another due to various factors such as SNR levels but typically its value ranges from 0.7–0.9. NCC_th should be large enough so that spurious landmark points are avoided but at the same time it should be low enough such that valid landmark points are not excluded. W_c is empirically set to 11 since it was observed that the maximum displacement in the field of view remains within this vicinity for all images in the sequence. The point-to-point correspondence in the I_i images follows exactly that determined from the J_i images.

2.4.3. Linear alignment of images via affine transformation

ARM implements the linear alignment of images in a sequential manner and takes into consideration the translation, rotation, scaling and shearing effects between successive images in the sequence. We begin with the alignment of I ₂ to I ₁ resulting in the aligned image I^a ₂ followed by the alignment of I ₃ to I ₁ and the sequential alignment process continues in this way until the last image I_n in the sequence is aligned. The alignment process at every step can be formally expressed as follows where

the superscript a denotes that I_i has been aligned to I ₁ and f ₂ represents the affine transformation function determined from the point-to-point correspondence of the landmark points between I ₁ and I ₂. Similarly, f_i is determined from the point-to-point correspondence of the landmark points between adjacent images I _i-1 and I_i. The motivation here is that transformation functions f_i, i=1,2, …,N between adjacent images can be more accurately determined compared to those between any given image and a fixed reference image. The reason is that the motion artifacts between adjacent images tend to be lesser as compared to images that are further apart.

But the retina is a biological tissue and as such affine transformation alone may not accurately compensate the complex respiratory motion between successive images in the sequence. In the next section, we propose a non-linear alignment approach which overcomes this limitation.

2.4.4. Non-linear alignment of images via rank matching

In the heart of ARM lies a pixel rank matching criterion which defines our non-linear alignment process. The rank of a pixel denotes the order of that pixel in an image neighborhood assuming all pixels in that neighborhood have been sorted in ascending order based on their intensity values. If image A is to be aligned to image B, the criterion requires the pixel values in image A to be updated such that the rank of the center pixel in every local neighborhood of image A is identical to that of the center pixel from a corresponding neighborhood of image B. This is a reasonable requirement since the rank of both pixels should not differ despite there being variations in the average intensity level in their neighborhood. However, this assumption is only valid in the absence of noise which could significantly alter the pixel ranks.

In order to ensure the noise robustness of our rank matching criterion, we utilize the rank information obtained from the aligned a-priori estimates J^a_i, i=1,2, …,45 computed by subjecting J_i to the same transformation functions as in Eq. 2. As in Section 2.4.3, the non-linear alignment is implemented sequentially as shown in Eq. 3 where

where the RANK operation is performed on a local (2W_r+1)×(2W_r+1) neighborhood centered at (x,y) and W_r is the width. The pixel (x_o,y_o) is located within this neighborhood and the superscript r in I^a,r_i denotes that I^a_i has been non-linearly aligned such that the ranks of its pixels match those of J^ref _i-1 and i≥2. J^ref_i is defined as

where J^ref ₁=J ₁ and J^a,r_i is computed as follows:

Finally, the result of averaging i number of aligned images in the sequence is given by

where I^ref ₁=I ₁. Eqs. 3–6 are implemented in a cyclic fashion where for every cycle, i is incremented before proceeding to the next cycle where Eqs. 3–6 are recomputed and I^ref_N denotes the final averaged image where N is the total number of images in the sequence.

3. Results

3.1. Estimating the underlying fluorescence signal

Figure 1(a) is a typical image of the optic disc in the retina of a transgenic mouse where the fluorescence signal is from the GFP protein which, in this case, labels the glial cells. As observed, the image has a very low SNR and this makes it extremely difficult to detect any suitable landmark points and it is clearly inappropriate to use such images as reference images in rank matching. An a-priori estimate is therefore required to accurately estimate the underlying fluorescence signal. We have explored several methods of determining this estimate and found that the Huygens deconvolution method (Huygens Essential Software, Scientific Volume Imaging BV) [27] gives the best results. The deconvolved result of Fig. 1(a) is shown in Fig. 1(b) and as observed, the speckle noise is removed and the fluorescence intensity profile is characterized by various intensity maximas. Conversely, a 9×9 median filtering method of estimating the a-priori estimate fails to remove the speckle noise and instead distorts the image as observed in Fig. 1(c). Figure 1(d) shows the results of applying anisotropic diffusion [5] where the noise appears to be removed but at the expense of diffusing the underlying fluorescence signal resulting in the image appearing blurred overall.

 

Fig. 1. (a) Input image, I_i, of the GFP signal from the optic disc of the transgenic mouse (b)–(d) J_i obtained after subjecting I_i to (b) Huygens deconvolution (c) 9×9 median filtering and (d) speckle reducing anisotropic diffusion.

Download Full Size | PDF

3.1.1. Tuning of the deconvolution parameters

The Huygens deconvolution method requires certain device information such as the pinhole diameter, numerical aperture, excitation wavelength etc. but assumes that the device in question is a wide-field or confocal microscope and not a scanning laser ophthalmoscope. We circumvent this problem by defining parameters for a hypothetical confocal microscope system such that the ratio of its Nyquist critical sampling density over the actual sampling density is comparable to that of the HRA 2 system specified in Section 2.3 (5.00µm/pixel and 1.13µm/pixel respectively). To this end, we define the confocal microscope parameters as shown in Table 1 so that the Nyquist critical and actual sampling density are 70nm and 15.8nm respectively. This ensures that both systems have the same degree of oversampling i.e. $\frac{5.00}{1.13}=\frac{70.0}{15.8}$ . The theoretical point spread function for deconvolution was subsequently computed based on the parameters specified in Table 1.

Table 1. Parameters of imaginary confocal microscope for the Huygens deconvolution method

View Table

3.2. Fast and accurate tracking of landmark points

We apply the NCC measure in Section 2.4.1 to automatically select intensity maximas as landmark points for every a-priori deconvolved image J_i in a sequence. Although a large number of landmark points may have been detected in an image J_i, only a subset of these points correspond to those in the adjacent image J _i-1 since we require that ${\mathrm{NCC}}_{J_{i,u,v},J_{i-1,u\prime ,v\prime }}\ge 0.9$ and that the corresponding points lie within a 23×23 square neighborhood of each other (W_c=11 as in Section 2.4.2).

 

Fig. 2. Point-to-point correspondence between: First row-(a) J ₁ and (b) J ₂, Second row-(c) J ₂ and (d) J ₃, Corresponding points in the adjacent images carry the same number.

Download Full Size | PDF

Figure 2 illustrates this where Figs. 2(a) and (b) are a-priori estimates of the first and second images in a sequence whereas Figs. 2(c) and (d) are a-priori estimates of the second and third images in the same sequence. As observed, each pair has their own unique set of corresponding landmark points indicated by the labeled green and red dots.

3.3. Quantitative analysis of ARM

Here, we quantify the influence of the window size W_r and image SNR on the alignment accuracy of the rank-matching criterion. This is achieved by omitting the affine transformation step of Section 2.4.3 and directly applying the rank-matching criterion of Section 2.4.4 on a sequence of 45 noisy synthetic images where each 256×256 image comprises a 2-D “ball” shaped intensity profile with a diameter of approximately 36 pixels. The “ball” is vertically displaced by two pixels between adjacent images. The Huygens deconvolution method was applied to obtain the a-priori deconvolved images J_i from the sequence.

3.3.1. Influence of W_r on alignment accuracy

In each case, the final averaged image I^ref ₄₅ is compared against the ideal (noiseless) synthetic image I′ and the mean absolute error (MAE) is computed as $\mathrm{MAE}=\frac{1}{M}\sum _{i=1}^M\mid I_{45,i}^{\mathrm{ref}}-I_i^{\prime }\mid$ and the results shown in Fig. 3(g) where i and M denote the pixel entry and total number of image pixels respectively. It is evident here that the distortion is smaller for larger window sizes and that it is minimal for W_r=19. Fig. 3(a) shows a sample image from a noisy sequence with an SNR of 0dB whereas the final averaged results in Figs. 3(b)–(e) illustrates the influence of W_r on the alignment accuracy. We observe that the averaged results in Figs. 3(b)–(d) appear truncated. This is due to the insufficient window size of (2W_r+1)×(2W_r+1) which limits the validity of the rank matching operation to the peak intensity region of the “ball” profile as observed in Figs. 3(b)–(d) although the correctly aligned region increases for larger window sizes. Accurate alignment of the entire “ball” profile can only be obtained if the window is large enough to cover the entire “ball” and its displaced counterpart in the adjacent image i.e. 2W_r+1 ≃ Ball diameter+Between-frame pixel displacement. The window size requirement is met for Fig. 3(e) and as observed, the distortion is minimal. In contrast, the 2-D “ball’ profile appears totally distorted in the direct frame averaged result of Fig. 3(f) due to the vertical displacements between adjacent images in the sequence.

3.3.2. Influence of SNR on alignment accuracy

In Figs. 4(d)–(f), we observe the efficacy of the ARM method in reducing the noise levels while minimizing distortion to image sequences with different SNR levels as shown in Figs. 4(a)–(c). A window size of W_r=19 is used in all three cases. Significantly small error values of about 1.9% and 0.5% were achieved for SNR levels of -10dB and 0dB respectively. Even for a poor SNR level of -20dB, a reasonably low MAE of about 6.5% was achieved. The efficacy of ARM is attributed to the noise robustness of the deconvolution method and the proposed rank matching criterion.

3.4. Qualitative analysis of ARM on noisy optic disc images of transgenic mice

We have also evaluated the performance of ARM on noisy optic disc images of transgenic GFAP–GFP mice and found that the results corroborate the findings made in Section 3.3 on the synthetic images.

 

Fig. 3. Influence of W_r on the alignment accuracy: (a) Reference image, (b)–(f) Final averaged images corresponding to (b) W_r=4, (c) W_r=9, (d) W_r=14, (e) W_r=19, (f) Direct averaging and (g) Corresponding drop in mean absolute error (MAE) error values from (b) to (e).

Download Full Size | PDF

3.4.1. Influence of Wr on alignment accuracy

As W_r increases, the resultant improvement in the alignment accuracy can be observed from the final averaged images I^ref_N of Fig. 5 where N=45. Image artifacts are visible in I^ref_N of Fig. 5(a) due to the small window size of W_r=2. Although the underlying fluorescence profile is taking shape in Fig. 5(b), it appears that only the maxima and minima intensity regions of the underlying fluorescence have been aligned accurately resulting in the appearance of tiny fluorescence dots and narrow branches of blood vessels in the optic disc. This is consistent with the observation made in Section 3.3.1. Accurate results can only be achieved if W_r≃17 since the sum of the “spot” diameter in the J_i images and the displacement between adjacent retinal images is approximately 35 pixels (2W_r+1=35). This window size corresponds to the result shown in Fig. 5(d). A slight deterioration in signal resolution is observed when the larger window size of W_r=22 is applied in Fig. 5(e) whereas minor image artifacts are observable for the smaller window size of W_r=12.

 

Fig. 4. Influence of W_r on the alignment accuracy: (a) Reference image, (b)–(f) Final averaged images corresponding to (b) W_r=4, (c) W_r=9, (d) W_r=14, (e) W_r=19, (f) Direct averaging and (g) Corresponding drop in mean absolute error (MAE) error values from (b) to (e).

Download Full Size | PDF

 

Fig. 5. Final averaged image I^ref_N of a transgenic mouse retina from the ARM method (a) W_r=2, (b) W_r=7, (c) W_r=12 (d) W_r=17 and (e) W_r=22.

Download Full Size | PDF

3.4.2. Comparison against conventional averaging methods

We also find that the final averaged image I^ref_N of Fig. 6(a), processed using ARM, has a significantly better resolution compared to those obtained using direct frame averaging in Fig. 6(b) as well as the averaging module in the commercial Heidelberg Retinal Angiograph 2 (HRA 2) software in Fig. 6(c). In addition, the noise level in Fig. 6(a) appears to be notably lower than those in Figs. 6(b) and (c). The poor resolution of the results in Fig. 6(b) and (c) makes it extremely difficult if not impossible to perform morphometric and intensity analysis of the fluorescence spots in the optic disc region.

4. Conclusions

We have presented a fully automated frame averaging method of in vivo GFAP–GFP fluorescent images of the transgenic mouse retina. It is based on a non linear alignment method we call averaging via rank matching (ARM) which ensures a significantly higher resolution and SNR of the final averaged image compared to conventional frame averaging and motion compensated frame averaging methods. Rank matching provides a significant improvement to the alignment accuracy of standard motion compensation operations such as affine transformation as it is robust against noise and variations in the fluorescence intensity of astrocytes between successive images. The method also enables an automated, accurate and robust detection of landmark points from the peak intensity points of the fluorescence profiles. These landmark points are subsequently used to establish point-to-point correspondence between successive images rather than between a fixed reference frame and all other images in a sequence since successive images are more closely aligned. Other potential applications of ARM include the averaging of images acquired from other modalities such as confocal microscopy, MRI and PET.

 

Fig. 6. Qualitative comparison of the composite GFAP–GFP transgenic mouse retina images, from the same sequence of noisy images, processed using (a) the proposed ARM method, (b) Direct frame averaging and (c) HRA 2’s averaging module.

Download Full Size | PDF