Automatic method of analysis and measurement of additional parameters of corneal deformation in the Corvis tonometer

Abstract

Introduction

The method for measuring intraocular pressure using the Corvis tonometer provides a sequence of images of corneal deformation. Deformations of the cornea are recorded using the ultra-high-speed Scheimpflug camera. This paper presents a new and reproducible method of analysis of corneal deformation images that allows for automatic measurements of new features, namely new three parameters unavailable in the original software.

Material and method

The images subjected to processing had a resolution of 200 × 576 × 140 pixels. They were acquired from the Corvis tonometer and simulation. In total 14000 2D images were analysed. The image analysis method proposed by the author automatically detects the edge of the cornea and sclera fragments. For this purpose, new methods of image analysis and processing proposed by the author as well as those well-known, such as Canny filter, binarization, median filtering etc., have been used. The presented algorithms were implemented in Matlab (version 7.11.0.584 - R2010b) with Image Processing toolbox (version 7.1 -R2010b) using both known algorithms for image analysis and processing and those proposed by the author.

Results

Owing to the proposed algorithm it is possible to determine three parameters: (1) the degree of the corneal reaction relative to the static position; (2) the corneal length changes; (3) the ratio of amplitude changes to the corneal deformation length. The corneal reaction is smaller by about 30.40% compared to its static position. The change in the corneal length during deformation is very small, approximately 1% of its original length. Parameter (3) enables to determine the applanation points with a correlation of 92% compared to the conventional method for calculating corneal flattening areas. The proposed algorithm provides reproducible results fully automatically within a few seconds/per patient using Core i7 processor.

Conclusions

Using the proposed algorithm, it is possible to measure new, additional parameters of corneal deformation, which are not available in the original software. The presented analysis method provides three new parameters of the corneal reaction. Detailed clinical studies based on this method will be presented in subsequent papers.

Introduction

Current advances in technology enable to measure intraocular pressure using a number of methods. These are the methods of non-contact and impression applanation tonometry [1, 2]. One such type is the Corvis tonometer which allows for the quantitative and qualitative (visual) assessment of biomechanical properties of the cornea. Based on a sequence of images of corneal deformation, the tonometer measures corneal thickness, deformation amplitude, applanation length, corneal deformation speed and intraocular pressure [3]. A schematic sequence of corneal deformation images and selected parameters are shown in Figure  1. In recent years, these parameters have been the subject of many papers and comparisons both among themselves as well as among other types of tonometers [4–47].

Figure 1

Measurement method with the use of the Corvis tonometer and the obtained results of a sequence of corneal deformation images. The cornea is subjected to an air puff in the Corvis tonometer. Consequently, every 230 μs an image of the corneal deformation in the line arranged on the main axis is obtained. Particular important deformation phases are shown on the left in the image a). Reactions of the cornea and the eyeball are shown on the right in the image b).

For example, there is a group of publications related to the analysis of corneal biomechanical parameters for the Ocular Response Analyzer (ORA) [4–7] and the impact of the patient’s age [8–12], glaucoma [13–18] and wound healing [19] on the results obtained. Also predictive numerical simulation of corneal biomechanical parameters [20] and intraocular pressure measured in vivo [21] have been analysed. A significant number of papers is devoted to corneal hysteresis and its association with glaucoma damage [22], pachymetry [23–25] or hysteresis measured in children and healthy patients using the Reichert ocular response analyzer [26, 27]. Corneal hysteresis was also analysed for porcine eyes or open-angle glaucoma [28–30]. A significant part of the papers deals with the relationship between keratoconus and its biomechanical properties [31–33]. The situation is similar to assisted laser in situ keratomileusis [34, 35], or the ultrastructure of the corneal stroma [36]. Therefore corneal biomechanical properties [37–41] are measured in various ways with different types of tonometers [42, 43] for various types of diseases [44, 47].

In the case of the Corvis tonometer, the comparative analysis carried out in the literature (e.g. Smedowski el at [3]) applies only to the parameters available in the device (intraocular pressure, pachymetry, applanation 1 time, applanation 1, length, applanation 1 velocity, applanation 2 time, applanation 2 length, applanation 2 velocity, highest concavity time, peak distance, radius, deformation amplitude maximum). The corneal deformation image analysis enables to determine a significantly larger number of interesting parameters than those available in the original software [48–50]. Few papers have been published so far in this area. These include the papers of Tejwani et al. [51] and Koprowski et al. [52–54]. The first one ( [51]) refers to spectral analysis and comparison of the ORA with the Corvis tonometer. The second and third papers ( [52, 54]) present the analysis of new features (e.g. reaction of the eyeball) obtained from the image analysis method proposed by the authors. In this paper, a new analysis method of the cornea images acquired from the Corvis tonometer is presented. It allows for automatic measurements of new, not yet published, features of the cornea - unavailable in the original software. These three new automatically designated parameters are: the degree of the corneal reaction relative to the static position, the corneal length changes, the ratio of amplitude changes to the corneal deformation length. They are described in detail in the following sections.

Material

As part of the paper, the correctness of the algorithm operation was tested on sequences of corneal deformation images obtained every 230 μs with a fixed M × N × I resolution of 200 × 576 × 140 pixels. The data of about 100 eyes were obtained from real data (open-access medical image repositories) and simulation data from the Corvis tonometer. Simulation data are derived from a corneal deformation generator, specially designed for this purpose, which provides any number of corneal deformations without patient’s participation. There were 140 2D images in each measurement, which in total gave 14000 2D images for analysis. The analysis presented in the following part of the paper enabled to enter image sequences in the source recording format *.cst or as a file *.avi.

Method

The new method of data analysis consists of two stages. (1) Image pre-processing- the data from the Corvis tonometer was subjected to analysis which enabled reconstruction of the cornea shape changes and separation of the eyeball reaction and corneal deformation; (2) The corneal reaction was analysed, which resulted in three new parameters.

The presented algorithms were implemented in Matlab (Version 7.11.0.584 - R2010b) with Image Processing toolbox (Version 7.1 -R2010b) using both known algorithms for image analysis and processing and those proposed by the author.

Pre-processing

The image pre-processing stage is partly known from previous publications of the author [52–54]. As mentioned earlier and in paper [52], the input images L _GRAY (m,n,i) (where m-row m∈(1,M), n-column n∈(1,N), and i – another 2D image i∈(1,I)) acquired from the Corvis tonometer had an M × N × I resolution of 200 × 576 × 140 pixels - Figure  2. Pixels for this type of image matrix symbols are numbered from one, first rows (m) and then columns (n) and the image number in the sequence (i). As a result, corneal deformation is shown in the three-dimensional graph (Cartesian coordinate system) in the reverse form relative to the image L _GRAY (m,n,i) visible in the Corvis tonometer. Input images were acquired directly in the *.cst format. First, a sequence of images L _GRAY (m,n,i) underwent median filtering with a mask h ₁ sized M _h1 ×N _h1 ×I _h1  = 3×3×3 pixels. The mask size was chosen arbitrarily, taking into consideration the size of artefacts and distortions that enter the optical path. Next, the filtered image L _M (m,n,i) was subjected to further preliminary transformations. These transformations were developed and their repeatability and accuracy were corrected in previous papers [52, 54]. These include:

Figure 2

Example of an input image L _GRAY ( m,n,i ). (where m-row m∈(1,M), n-column n∈(1,N), for i = 45 pixel) acquired from the Corvis tonometer (resolution M × N = 200 × 576 pixels) and the outer corneal contour L _T (n,i = 45) (red).

• detection of the outer edge of the cornea (using the Canny method [47]) - result L _C (m,n,i),

• morphological close operation,

• corneal contour L _T (n,i) (Figure  2),

• division of the contour L _T (n,i) into the sum of components: L _TD (n,i) – constant component of the corneal shape for t = 0 (for i = 1); L _TR (n,i) – corneal deformation and L _TO (n,i) – reaction of the eyeball (which is shown schematically in Figure  1).

Designation of the waveform L _TO (n,i) is possible owing to the analysis of the visible contour of the sclera at the border of the left and right image, which is shown schematically in Figure  3, i.e.:

$$L_{TO}\left(N,i\right)=\frac{1}{N_k}{\displaystyle \sum _{n=0}^{N_k}}\left(L_T\left(N-n,i\right)-L_{TD}\left(N-n,1\right)\right)$$

(1)

$$L_{TO}\left(1,i\right)=\frac{1}{N_k}{\displaystyle \sum _{n=1}^{N_k}}\left(L_T\left(n,i\right)-L_{TD}\left(n,1\right)\right)$$

(2)

Figure 3

Schematic diagram of the method for designating the boundary points of the reaction of the eyeball. The waveforms L _TO (N,i) and L _TO (1,i) determined on this basis, shown demonstratively on the left and right side, are the basis for determining the missing values for n∈(2,3,…,N-2,N-1).

The waveforms L _TO (N,i) and L _TO (1,i) designated based on the formulas (1) and (2) are the basis for determining the missing values for n∈(2,3,…,N-2,N-1). These values were determined based on linear interpolation. The results obtained, namely L _TD (n,i), L _TR (n,i), L _TO (n,i) and the summary waveform L _T (n,i), are shown in Figure  4. They are the basis for further processing steps presented hereafter.

Figure 4

Sample results obtained from the Corvis device. The waveform L _T (n,i) visible in figure a) is the response of the cornea and eyeball with a constant component to an air puff. On this basis and as a result of image pre-processing, L _TD (n,i), L _TR (n,i) and L _TO (n,i) were automatically separated. Figure b) shows the results obtained for n = N/2 in the form of a two-dimensional graph.

Processing

The results obtained, namely L _TD (n,i), L _TR (n,i), L _TO (n,i) and the summary waveform L _T (n,i), are subjected to further transformations, which allows for automatic determination of three new parameters:

• the degree of the corneal reaction relative to the static position,

• the corneal length changes,

• the ratio of amplitude changes to the corneal deformation length.

Determination of the degree of the corneal reaction relative to the static position requires measurement of the length L _S (n,i) and, on its basis, L _g (n,i). The value L _S (n,i) is the length of the cornea subjected to deformation. In contrast, L _g (n,i) is the surface area of the difference between corneal deformation and its inverted static shape. The measurement method is shown schematically in Figure  5. The value L _W (n,i) necessary for calculating L _S (n,i) is determined on the basis of L _TR (n,i), i.e.:

$$L_W\left(n,i\right)=\left\{\begin{array}{ccc}\hfill 1\hfill & \hfill \mathit{if}\hfill & \hfill L_{TR}\left(n,i\right)>p_r\hfill \\ \hfill 0\hfill & \hfill \mathit{other}\hfill \end{array}\right.$$

(3)

where p _r – binarization threshold set at the maximum possible amplitude of the noise, i.e. approximately two pixels [55].

Figure 5

Schematic diagram of the measurement method of the values L _S ( n , i ) and L _g ( n , i ). Part a) shows the measurement of the distance L _d as pick distance and the distance Ls as the corneal length where deformation relative to the static position is measured. Part b) shows the measured deformation value L _b as a measure of the difference in relation to the original position.

Next, after median filtering of the waveform (L _W (n,i)) with a filter mask h sized M _h  × N _h  = 3 × 3 pixels, the value L _S (n,i) is obtained. In the next stage, a straight line L _OS (the axis of symmetry) is drawn in the range for which L _S (n,i) = 1 for the extreme points c ₁ and c ₂ (Figure  5.) with the coordinates (m ₁ ,n ₁ ) and (m ₂ ,n ₂ ) respectively, i.e.:

$$L_{\mathit{OS}}\left(n,i\right)=\frac{\left(m_2-m_1\right)\cdot \left(n-n_1\right)}{\left(n_2-n_1\right)}+m_1$$

(4)

Then each point of the corneal contour L _T (n,i) is reflected symmetrically with respect to the axis of symmetry L _OS (n,i) for which L _S (n,i) = 1, which gives L _K (n,i), i.e.:

$$L_K\left(n,i\right)=\left\{\begin{array}{ccc}\hfill 2\cdot L_{\mathit{OS}}\left(n,i\right)-L_T\left(n,i\right)\hfill & \hfill \mathit{if}\hfill & \hfill L_S\left(n,i\right)=1\hfill \\ \hfill L_T\left(n,i\right)-L_{TD}\left(n,i\right)\hfill & \hfill \mathit{other}\hfill \end{array}\right.$$

(5)

The result L _K (n,i) is shown schematically with the dotted line in Figure  5.b). The results L _K (n,i) and L _TR (n,i), after the removal of the constant component (the static corneal contour L _TD (n,i)), are shown in Figure  6.

Figure 6

Graphs L _K ( n , i ) and L _TR ( n , i ). Part a) shows three-dimensional graphs L _K (n,i) and L _TR (n,i) and their difference, whereas part b) shows two-dimensional graphs for the sample value i = 68 pixels.

Determination of the corneal length changes is the second calculated parameter and is directly linked to the results L _K (n,i) and L _TR (n,i). Changes in the corneal length are calculated within the range for which L _S (n,i) = 1. Three values are calculated, the length of the corneal deformation L _DDTR (i), the corneal length in a stable state covering the deformation area L _DDK (i) and the distance between the extreme points of deformation L _DDS (i) (Figure  5):

$$L_{DDK}\left(i\right)={\displaystyle \sum _{n=1}^{N-1}}\sqrt{1+{\left(L_K\left(n+1,i\right)-L_K\left(n,i\right)\right)}^2}$$

(6)

$$L_{\mathit{DDTR}}\left(i\right)={\displaystyle \sum _{n=1}^{N-1}}\sqrt{1+{\left(L_{TR}\left(n+1,i\right)-L_{TR}\left(n,i\right)\right)}^2}$$

(7)

$$L_{DDS}\left(i\right)={\displaystyle \sum _{n=1}^N}{L}_s\left(n,i\right)$$

(8)

The obtained measurement results are shown in Figure 7a and b. Additionally, Figure  7b shows the range of measurement accuracy calculated as the sum of ± LSB (Least Significant Bit) for each distance between adjacent points. It results from high sensitivity of formulas (6) and (7) to noise, the sum of relationships L _K (n + 1,i)-L _K (n,i) and L _TR (n,i)-L _TR (n,i). It is therefore impossible at this stage to clearly assess whether the cornea changed its length during deformation. Accordingly, the approach was modified and waveforms L _K (n,i) and L _TR (n,i) were approximated with a polynomial of degree 8, thus obtaining L _K2 (n,i) and L _TR2 (n,i). The degree of the polynomial was chosen on the basis of the performed measurements of the best match. The corneal length measurement, analogously to equations (6) and (7), was marked as L _AAK (i) and L _AATR (i):

$$L_{AAK}\left(i\right)=\left\{\begin{array}{ccc}\hfill {\displaystyle \sum _{n=1}^{N-1}}\sqrt{1+{\left(\frac{\partial L_{K2}\left(n,i\right)}{\partial n}\right)}^2}\hfill & \hfill \mathit{if}\hfill & \hfill L_{DDS}\left(i\right)>0\hfill \\ \hfill 0\hfill & \hfill \mathit{other}\hfill \end{array}\right.$$

(9)

$$L_{\mathit{AATR}}\left(i\right)=\left\{\begin{array}{ccc}\hfill {\displaystyle \sum _{n=1}^{N-1}}\sqrt{1+{\left(\frac{\partial L_{TR2}\left(n,i\right)}{\partial n}\right)}^2}\hfill & \hfill \mathit{if}\hfill & \hfill L_{DDS}\left(i\right)>0\hfill \\ \hfill 0\hfill & \hfill \mathit{other}\hfill \end{array}\right.$$

(10)

Figure 7

Graph of distance between the extreme points of deformation L _DDS ( n , i ) against the corneal reaction a), and a graph of corneal deformation length L _DDTR ( n , i ) and the corneal length in a stable state covering the deformation area - L _DDK ( n , i ) b).

The sample results of changes in the values L _AAK (i) and L _AATR (i) are shown in Figure  8. The differences between the waveforms of corneal length changes L _AAK (i) and L _AATR (i) visible in Figure  8b) are in this case 4 ± 2 pixels.

Figure 8

Graph of corneal length during deformation and in a stable state L _AAK ( i ) and L _AATR ( i ) a) and the length difference L _AAK ( i )- L _AATR ( i ) b).

The ratio of the corneal deformation length to deformation amplitude

This parameter was previously calculated using equation (3) based on which the value L _S (n,i) and then L _DDS (i) are calculated (formula (8)). Thus, the ratio of the corneal deformation length to corneal deformation amplitude was calculated on the basis of L _DDS (i) and L _TR (n,i) as L _DDS/TR (i), i.e.:

$$L_{DDS/TR}\left(i\right)=\frac{L_{DDS}\left(i\right)}{\underset{n}{max}{L}_{TR}\left(n,i\right)}$$

(11)

assuming $\underset{n}{max}{L}_{TR}\left(n,i\right)\ne 0$.

A sample graph of L _DDS/TR (i) is shown in Figure  9. The visible two peaks designate two points for which the corneal deformation length is greatest in relation to the deformation amplitude. The impact of changes in binarization threshold p _r (formula (3)) on the obtained results is shown in Figure  9b) for p _r ∈(0.1,6).The results obtained using the presented new methods of automated analysis are discussed in the next section. The individual processing steps, discussed above, are shown as a block diagram in Figure 10.

Figure 9

Sample graph L _DDS/TR ( i ). The positions of the two peaks on the right and left side define the points of corneal applanation. Part a) shows the graph L _DDS/TR (i) with the graph L _TR (i) in the background, whereas part b) shows the impact of changes in the binarization threshold p _r on the measurement results obtained - the shape of the graph L _DDS/TR (i).

Figure 10

Block diagram of the various processing stages. In image pre-processing, data from the Corvis tonometer were analysed, which enabled to reconstruct changes in the shape of the cornea and separate the reaction of the eyeball. In the next stage of data processing, a detailed analysis of the response of the eyeball was performed.

Results

The proposed new algorithm allows for the automatic measurement of three new parameters of the corneal response during deformation using the Corvis tonometer. The calculation of the degree of the corneal reaction relative to the static position enabled automatic measurement of the difference between these waveforms, L _K (n,i = 68) and L _TR (n,i = 68), which was 30 pixels for n = 403 pixels (Figure  6). Based on the analysed 14000 2D images, the corneal reaction was always lower compared to the corresponding static position. For all the analysed images L _K (n,i = 68)-L _TR (n,i = 68) > 0 in each case. Due to the high sensitivity of the method (formulas (6)(7)), the analysis of changes in the corneal length L _DDK (i) and L _DDTR (i) must be preceded by filtering or interpolation. The differences between the waveforms of changes in the corneal length L _AAK (i) and L _AATR (i) visible in Figure  8b) are in this case 4 ± 2 pixels (about 1% in relation to its original length). These differences are different for various groups of images and on average they are a few pixels. In the analysis of this group of results, the measurement error and accuracy of obtaining the results L _K2 (i), L _TR2 (i) and, in particular, L _T (n,i), should be taken into account. Errors may arise in the initial phase of the algorithm operation in the form of not fully visible corneal contours. They replicate here introducing a substantial measurement error [56]. They also affect the results obtained from the calculation of the third parameter – the ratio of deformation to the corneal deformation length. This parameter designates two peaks (Figure  9). From a practical point of view, the position of the two peaks and the shape of the graph L _DDS/TR (i) between them are important. This fragment on the graph (Figure  9) runs stably for the approximate value of i∈(50,85). It means that the ratio of the deformation to deformation amplitude is constant. The cornea during deformation causes proportional changes in the deformation length. The position of the peaks (a ₁ , a ₂ ) is slightly dependent on the accuracy of the earlier analysis and especially on the adopted binarization threshold p _r . Tab. 1 shows the obtained results of changes in the position (relative to i) of the two peaks (a₁, a₂), for the adopted different binarization thresholds, i.e. for p _r ∈(0.1,6). According to the presented results, changes in peak positions are at a level of the resolution error (±1 pixel) for small values of the threshold p _r . In a further step, the correlation between the position of the peaks and applanation points, calculated in the standard software of the Corvis tonometer, was assessed. This correlation was 92%, which confirms the usefulness of the new discussed method for the automatic designation of applanation points. Typical analysis of applanation points may also be carried out in a conventional manner as a search for time instants (values i) for which corneal flattening occurs. This is a difficult task in terms of algorithmics due to the need to correct the angle of the cornea position relative to the tonometer or difficulties in the unambiguous determination of corneal flattening areas. Particular attention should be paid here to the calculation of the above errors that apply only to the discussed method. In practice, especially when calculating their values for a different type of a tonometer, for a different type of Scheimpflug camera, close attention should be paid to the accuracy and reproducibility of obtaining images L _GRAY (m,n,i). This also applies to data derived from simulation as well as, for example, data from the measured eye phantom. Diversity of technological parameters of the Corvis tonometer and the impact of other factors on the result and reproducibility of measurements should be taken into account: e.g. the impact of the patient’s head position or the effect of temperature and humidity. In extreme cases, these errors can skew each other or greatly increase the total measurement error.

Table 1 Results of position changes (relative to i ) of two peaks (applanation points) for the adopted different binarization thresholds, p _r ∈(0.1,6)

The presented new features enable to extend the range of possibilities for quantitative assessment of the corneal response and to find links with known features measured in the Corvis tonometer. The proposed algorithm provides reproducible results fully automatically using Core i7 10GB RAM within a few seconds per patient.

Critical summary

The paper presents a method for calculating three additional parameters of the quantitative assessment of corneal deformation. The presented algorithm for image analysis provides reproducible results in a fully automated manner. The presented algorithm is only one of several possible solutions to this task. In particular, the described profiled algorithm can be composed of other image analysis and processing methods [56]. For example, the presented problem of image analysis and processing can be solved using other methods [57, 58] derived from morphological transformations (open, close etc.) [59–61]. However, in any case, the algorithm must be profiled to a given application. Proper selection of the algorithm structure is an important part taking into account, on the one hand, full automation, and, on the other hand, large variation in images (both their quality as well as the parameters of the cornea itself). Due to full automation of the measurement, there is a need to profile the algorithm for a particular research problem and a group of images. On the other hand, the algorithm requires generalization in order to ensure its correct operation in different research institutions. Its final form is thus a compromise between these two elements.

For example, in paper [52], the authors present the division of corneal response into four different groups depending on the characteristics of their response to an air puff. New parameters that are not calculated in the Corvis tonometer are introduced there. For example, in papers [47, 62–65], conclusions and the whole analysis are only related to the parameters available from the Corvis tonometer. Paper [54] shows the correlation of the new parameters obtained from the Corvis tonometer with known parameters available in the proprietary software of the Corvis tonometer.

In subsequent papers, the author intends to:

• verify the repeatability of the three parameters obtained for the same patients measured with different types of the Corvis tonometer at different medical centres [56, 66],

• perform quantitative analysis of the impact of patient positioning on the results obtained [66].

Biomechanical parameters of the eye, measured in dynamic states, still represent a new area of interdisciplinary research combining engineering, medicine and computer science [67, 68].