METHOD FOR AUTOMATIC SHAPE QUANTIFICATION OF AN OPTIC NERVE HEAD

Abstract

The invention relates to a method and a computer program for automatic shape quantification of an optic nerve head from three-dimensional image data (1) acquired with optical coherence tomography, comprising the steps of: a) Providing (100) three-dimensional image data (1) of the retina, the image data comprising at least a portion of the optic nerve head, wherein the image data comprises pixels with associated pixel values; b) In the three-dimensional image data (1) identifying (200, 300) anatomic portions of the optic nerve head, the anatomic portions comprising a retinal pigment epithelium (RPE) portion (3) and an inner limiting membrane (ILM) portion (2); c) Determining an RPE polygon mesh (30) for a lower boundary of the retinal pigment epithelium portion (3), wherein the RPE polygon mesh (30) extends along the lower boundary of the retinal pigment epithelium portion (3); d) Determining an ILM polygon mesh (20) for the inner limiting membrane portion (2), wherein the ILM polygon mesh (20) extends along the inner limiting membrane portion (2); e) Determining a morphologic parameter (10) of the optic nerve head from the RPE polygon mesh (30) and the ILM polygon mesh (20); f) Displaying the morphologic parameter (10) of the optic nerve head and/or a representation of at least a portion of the RPE polygon mesh (30) and/or a representation of at least a portion of the ILM polygon mesh (20).

Claims

1. A method for automatic shape quantification of an optic nerve head from three-dimensional image data (1) acquired with optical coherence tomography, comprising the steps of: a) Providing (100) three-dimensional image data (1) of the retina, the image data comprising at least a portion of the optic nerve head, wherein the image data comprises pixels with associated pixel values; b) In the three-dimensional image data (1) identifying (200, 300) anatomic portions of the optic nerve head, the anatomic portions comprising a retinal pigment epithelium (RPE) portion (3) and an inner limiting membrane (ILM) portion (2); c) Determining an RPE polygon mesh (30) for a lower boundary of the retinal pigment epithelium portion (3), wherein the RPE polygon mesh (30) extends along the lower boundary of the retinal pigment epithelium portion (3); d) Determining an ILM polygon mesh (20) for the inner limiting membrane portion (2), wherein the ILM polygon mesh (20) extends along the inner limiting membrane portion (2); e) Determining a morphologic parameter (10) of the optic nerve head from the RPE polygon mesh (30) and the ILM polygon mesh (20); f) Displaying the morphologic parameter (10) of the optic nerve head and/or a representation of at least a portion of the RPE polygon mesh (30) and/or a representation of at least a portion of the ILM polygon mesh (20).

2. Method according to claim 1, wherein the RPE polygon mesh (30) and the ILM polygon mesh (20) each comprise vertices and faces (22, 32) and edges (23, 33), wherein for each face (32) of the ILM polygon mesh (20) a corresponding vertex of the RPE polygon mesh (30) is determined (400), such that a correspondence between the ILM polygon mesh (20) and RPE polygon mesh (30) is determined.

3. Method according to claim 1, wherein in the image data (1) openings in the Bruch's membrane (BMO) (43) are automatically identified, wherein from the identified Bruch's membrane openings (43) a set of points (40) is generated representing the Bruch's membrane openings, wherein the set of points (40) encloses a Bruch's membrane opening region (42).

4. Method according to claim 3, wherein an oval (41) is fitted to the points (40) representing the Bruch's membrane openings, wherein the oval (41) encloses the Bruch's membrane opening region (42) and wherein the oval (41) is comprised by a plane extending orthogonally to an A-scan direction of the image data (1).

5. Method according to claim 4, wherein the oval (41) is transformed to a circle by an affine transformation.

6. Method according to claim 3, wherein for the RPE polygon mesh (30) a central RPE region (34) and for the ILM polygon mesh (20) a central ILM region (24) is determined (500), wherein the respective central region (24, 34) is enclosed by an area corresponding to a projection of the oval (41) onto the respective polygon mesh (20, 30), wherein the projection is along the A-scan direction of the image data (1).

7. Method according to claim 1, wherein an optic nerve head cup volume (7) is determined, wherein the optic nerve head cup volume (7) is calculated from the ILM and RPE polygon mesh (20, 30), wherein the optic nerve head cup volume (7) is the volume enclosed by the portion of the ILM polygon mesh (25) that extends below the RPE polygon mesh (30).

8. Method according to claim 1, wherein an optic nerve head cup surface (8) is calculated from the RPE polygon mesh (30) and the ILM polygon mesh (20), wherein the optic nerve head cup surface (8) corresponds to the surface portion of the ILM polygon mesh (25) that extends below the RPE polygon mesh (30).

9. Method according to claim 1, wherein a bending energy (6) is determined for the inner limiting membrane portion (2) from the ILM polygon mesh (20), wherein the bending energy (6) corresponds to a magnitude of curvature of the ILM polygon mesh (20), wherein the bending energy (6) is determined only for the central ILM region (25).

10. Method according to claim 1, wherein the RPE polygon mesh (30) and/or the ILM polygon mesh (20) is a triangulated surface, wherein for each face (22) of the ILM triangulated surface (20) a corresponding vertex of the RPE triangulated surface (30) is determined, such that a correspondence between the ILM triangulated surface (20) and the RPE triangulated surface (30) is established.

11. Method according to claim 1, wherein prior to the determination of the RPE polygon mesh (30) for the retinal pigment epithelium portion (3) said portion is smoothed by means of a thin plate spline method, wherein the smoothing is applied to the lower boundary of the retinal pigment epithelium or wherein the smoothing is applied only to the lower boundary of the retinal pigment epithelium, wherein the RPE polygon mesh (30) is determined for the lower boundary of the smoothed retinal pigment epithelium portion.

12. Method according to claim 1, wherein prior to the determination of the RPE polygon mesh (30) for the lower boundary of the retinal pigment epithelium portion (2) said portion (2) is interpolated by means of a thin plate spline method for the Bruch's membrane opening region (42), wherein the RPE polygon mesh (30) is determined for the interpolated lower boundary of the retinal pigment epithelium portion.

13. Method according to claim 1, wherein a transformation is determined under which retinal pigment epithelium portion (2) becomes planar, wherein the image data (1) and/or the ILM polygon mesh (20), the RPE polygon mesh (30) and the Burch membrane openings are transformed by the same transformation.

14. A computer program comprising instructions to cause a computer to execute the steps of the method of claim 1.

Description

[0105] In the following, a description of exemplary embodiments is given by means of examples and figure description. It is shown in

[0106] FIG. 1 a flow diagram depicting selected processing steps of the method according to the invention;

[0107] FIG. 2 a B-scan with highlighted processing steps according to the invention;

[0108] FIG. 3 a RPE surface with motion artefacts and the same RPE surface after smoothing.

[0109] FIG. 4 representation of the RPE and ILM polygon meshes as well as the points representing the Bruch's membrane openings;

[0110] FIG. 5 Set of points representing the Bruch's membrane openings and the fitted ellipse enclosing the Bruch's membrane opening region;

[0111] FIG. 6 RPE polygon mesh with projected ILM polygon mesh for determining the corresponding faces and vertices;

[0112] FIG. 7 a visual representation of different morphologic parameters;

[0113] FIG. 8 a B-scan of the ONH, where the Bruch's membrane is visible close to the ONH; and

[0114] FIG. 9 a representation of various morphological parameters.

[0115] The invention relates to a fully automatic 3D, i.e. three-dimensional shape analysis of the optic nerve head (ONH) region. According to some embodiments it is possible to calculate novel three-dimensional shape parameters characterizing the ONH and to provide robust and reliable three-dimensional morphologic parameters that describe different aspects of the various shapes of the ONH.

[0116] In the following, the procedure to compute several three-dimensional morphologic, particularly shape parameters of the ONH including the preprocessing of the ILM surface and the RPE surface and a correspondence between them, is demonstrated.

[0117] FIG. 1 shows a flow diagram, where an OCT volume scan is the input of the method, i.e. the image data provided 100 to the method comprise an OCT volume scan.

[0118] Next, the RPE surface, representing the Bruch's membrane is segmented 200 together with the Bruch's membrane opening (BMO) points. The term “RPE surface” is used interchangeably with the term “RPE polygon mesh” representing the lower boundary of the RPE portion.

[0119] It is noted that in the context of the specification the lower boundary of the RPE is identical to the Bruch's membrane layer.

[0120] In a subsequent or parallel step, the ILM surface, also referred to as the ILM polygon mesh in the context of the current specification is determined 300 from the image data.

[0121] These three anatomic portions, represented by the RPE polygon mesh, the ILM polygon mesh and the Bruch's membrane openings can be used to represent the ONH three-dimensional shape and serve as inputs for further shape analysis of the ONH.

[0122] For each face of the ILM polygon mesh, a corresponding vertex of the RPE polygon mesh is determined 400. Furthermore, a central region is determined 500 for both polygon meshes by means of a region enclosed by the BMO points, such that at least one morphologic parameter can be estimated and displayed to a user of the method.

[0123] The ILM and the RPE surfaces will be represented as M.sub.ilm and M.sub.rpe, respectively and the BMO points are denoted by P. In this example, the ILM and the RPE surfaces are triangulated manifold surfaces, i.e. specific polygon meshes, and can be written in terms of a set of vertices and faces (triangles):

M.sub.ilm={V.sub.ilm,F.sub.ilm} and M.sub.rpe={V.sub.rpe,F.sub.rpe}

[0124] The ILM and RPE surface can have different numbers of vertices and number of triangles, wherein n.sub.ilm and m.sub.ilm represent the numbers of vertices and faces in the ILM surfaces. Analogously, n.sub.rpe and m.sub.rpe denote the size of V.sub.rpe and F.sub.rpe respectively. The BMO points are represented as:

P={p.sub.i∈R.sup.3;i=1 . . . n.sub.p}, where n.sub.p is the number of the BMO points.

[0125] The volume scans of the ONH scans can be obtained with a spectral-domain OCT (Heidelberg Spectralis SDOCT, Heidelberg Engineering, Germany) using a custom protocol with 145 B-scans, focusing the optic nerve head with a scanning angle of 15°×15° and a resolution of 384 A-scans per B-scan. The spatial resolution in x direction is approximately 12.6 μm, in axial direction approximately 3.9 μm and the distance between two B-scans approximately 33.5 μm.

[0126] FIG. 2 shows the Bruch's membrane surface and BMO points detection exemplified using a single B-scan. FIG. 2(a) shows one B-scan of the original image data (gray arrow shows blood vessels and the shadow artifacts these produce, the delimited region is part of the ONH disc. FIG. 2 (b) shows the same B-scan after applying smoothing and intensity normalization. FIG. 2 (c) shows approximated ISOS Junction points (dotted line) after removing outliers. As a convention in the computation, the outliers detected at this step are set to have in the axial (z) direction the coordinate equal to 1. FIG. 2 (d) shows the same B-scan with an approximated upper boundary of the RPE portion, i.e. upper boundary RPE points (dotted line). FIG. 2 (e) shows the same B-scan with the smoothed and interpolated lower boundary of the RPE portion 30 (dotted line). FIG. 2(f) shows the same B-scan with the identified BMO points 40 (white dots); even in the presence of blood vessels the BMO points 40 are still detected. 384 pixels represent approximately 4402.80 μm and 226 pixels approximately 881.40 μm.

[0127] In the following a detailed example is given for generating the RPE polygon mesh 30 for the Bruch's membrane 44.

Determining the Lower Boundary of the RPE

[0128] The Bruch's membrane 44 represents the termination of the retina and is therefore an important parameter in morphometric computations. According to this example, the RPE polygon mesh 30 is determined such that it represents the Bruch's membrane. One way to segment the Bruch's membrane is taught in [1]. Several preprocessing steps are performed that are commonly employed in image data from OCT. Consider I(q.sub.xy) the intensity of a pixel q.sub.xy. In a first step, a Gaussian smoothing filter (σ=5 pixels isotropic with kernel size=(10 μm×14 μm) is applied on each B-scan separately. The smoothing by the Gauss filter not only reduces speckle noise present in most OCT data, but also facilitates the approximation of the two most hyperintense layers, namely the retinal nerve fiber layer (RNFL) and the RPE.

[0129] In a second step, in order to address varying intensities in the image data, a contrast rescaling on each slice (B-scan) is performed. Contrast inhomogeneities can occur in form of a B-scan having regions with different illumination or as several B-scans of the same volume with very different intensity ranges. Specifically, a histogram-based amplitude normalization method [2] is used to map the signals in the original image linearly between the pixel values [0; 1] using as low cut off the first 66.sup.th percentiles and as high cut off the 99.sup.th percentile on the histogram of the B-scan where the sampled column (A-scan) is located.

[0130] FIG. 2(b) shows a B-scan of the volume data that has been smoothed and normalized according to the above procedure. FIG. 2(a) shows the same B-scan with its original grey values.

Upper Boundary RPE Approximation

[0131] In order to approximate the upper boundary of the RPE portion 3, first the ILM portion 2 is approximated as an upper boundary for the segmenting method. At each A-scan, the first pixel from top (i.e. on the side where the vitreous 51 is located (see e.g. FIG. 2) in the smoothed and normalized volume data, I.sub.SN is selected, that is higher than ⅓ of the maximum value in the B-scan containing the A-scan. This provides a set of initial estimate points for the ILM, denoted by ILM.sub.init. Next, the upper boundary of RPE is approximated. First, the image derivative, ∇I.sub.SN, of each B-scan (vertical gradient) using a Sobel kernel is calculated. Along each A-scan, starting from the ILM.sub.init set, several intermediate steps for the RPE upper boundary approximation are performed. Inner and outer segment Junction regions (ISOS) are determined by finding the first set of points p, shown in FIG. 2(c):

p.sub.xy={∇I.sub.SN(qi.sub.xy)<∇I.sub.SN(s.sub.xy)/3},s.sub.xy∈ILM.sub.init.

[0132] Starting from the set p.sub.xy, the RPE upper boundary is approximated:

RPE.sub.upper=max(∇I.sub.SN({tilde over (p)}.sub.xy)),

wherein only the points below IS/OS are considered:

{tilde over (p)}.sub.xy∈{∥s.sub.i.sub.xy−p.sub.i.sub.xy∥<60μm}.

[0133] At this stage the input is a list of points that belong to the upper boundary of RPE in each B scan. This list comprises among the points correctly positioned at the upper RPE boundary, also several outliers, especially in the presence of shadows cast by blood vessels, as well as at the region of the optic nerve head. In order to remove the outliers that lie in the upper or lower third of each B-scan, the gradient of a line consisting of upper RPE points reflecting the RPE upper boundary can be determined, and the mean value of said upper RPE points can be computed from coordinates that most likely belong to the correct upper RPE points. These coordinates represent RPE boundary points from the largest part of the gradient line between outliers (outliers in the gradient are considered to be >40 μm). The first seed point is then detected as the one closest to the mean value. Starting from this seed a q.sub.seed, outliers are iteratively removed from RPE.sub.upper (points where ∥q.sub.seed−q.sub.new∥>70 μm). Similarly, outliers from the last third of the B-scan are removed. The resulting point set of one B-scan is shown by the white points in FIG. 2(d). The points removed from the RPE.sub.upper roughly estimate the ONH region, as well as the BMO area. Note that the ILM can have a very complex topology, while other retinal layers are missing in this area. A mask of the ONH from the removed RPE.sub.upper part by fitting an ellipse to its contour, A.sub.ONH is created.

RPE Lower Boundary Detection

[0134] The lower boundary of the RPE is denoted with RPE.sub.lower and is identical to the Bruch's membrane. Points with the largest negative gradient below RPE.sub.upper, closest to the RPE.sub.upper are selected as the lower boundary of the RPE (i.e. if several minimum points have similar values, the point with the smallest distance to the corresponding RPE.sub.upper is selected). Using only the maximum gradient values leads to spurious points along each surface.

[0135] Correction of these errors is done by applying a cubic smoothing spline with a high smoothing parameter. Note that in case of presence of blood vessels, large regions of missing coordinates for the RPE might occur and the cubic spline can present deviations from the desired smooth contour.

[0136] Finally, to account for motion artifacts in consecutive B-scans, but also for the natural curvature of the retina, an efficient two-stage thin-plate spline fitting (TPS) process is performed, which improves the approach proposed by [3] without making use of the orthogonal scans presented in the work of [4]. First, TPS least-square approximation is performed. The number of control points used is determined by the size of the surface along each axial dimension. At this stage, the number is set to 25% in the slow scan direction, 15%, respectively, and the control points are evenly distributed along each direction. This enables the TPS, in combination with a smoothing parameter α (c.f. [5]) set to 0.85, to create a more smoothed surface that keeps the curvature of the retina without being influenced from motion artifacts especially along the slow axis. Particularly values of α=[0.70; 0.85] provide consistent results. Extreme grid points in the original surface defined as mean+standard deviation in local non-overlapping neighborhoods of 10×10 grid points of the TPS surface are removed. Then, the actual TPS fitting similar is applied. The choice of parameters at this second step is strongly influenced by the fact that outliers at a previous stage have been removed. Specifically, for grid points 20% in the slow scan direction are used, 10% in the fast scan direction, with smoothing parameter 0.45. Consistent results were obtained for α=[0.30; 0.50]. The strategy for is to obtain a TPS closer to the data in the grid points, while smoothing the artifacts present in the position of the detected RPE.sub.lower points especially at the presence of blood vessels, or in the close vicinity of the approximated ONH region. Both stages are done on the RPE.sub.lower without including A.sub.ONH. FIG. 3(a) shows the original RPE surface with typical artifacts in the in-between B-scans direction. These are corrected after applying our TPS approach, while keeping the shape of the original surface. The result is presented in FIG. 3(b).

[0137] The result of the RPE lower boundary 3, i.e. the Bruch's membrane, after performing the TPS is shown in FIG. 2(e).

[0138] Once also the ILM surface is determined, a correspondence between the RPE polygon mesh and the ILM polygon mesh is established. For further analysis it is necessary to find the corresponding points between these two surfaces. This process is demonstrated in the following. Vertices in the RPE surface corresponding to each face (f.sub.i∈F.sub.ilm|i=0, . . . , m.sub.ilm−1) of the ILM surface are computed. In general, the RPE surface, represented here as a function graph: M.sub.rpe: R.sup.2.fwdarw.R, has a less complex structure compared to the ILM. In the OCT scanner, the number of the A-scans (x-direction) and the number of the B-scans (y-direction) are fixed, which creates a regular XY-grid as a domain for the RPE graph function. Therefore, the index of each vertex of the RPE surface can be computed using the numbers of x-lines (vertical lines) and y-lines (horizontal lines) and the sampling size in both directions, denoted by ϵ.sub.x and ϵ.sub.y respectively. The numbers of x-lines and y-lines are computed using the following equation:

[00001] $xline = \frac{x_{\max} - x_{\min}}{ϵ_{x}} + 1, yline = \frac{y_{\max} - y_{\min}}{ϵ_{y}} + 1,$

where x.sub.max, x.sub.min, y.sub.max, and y.sub.min are the bounding values of the RPE surface in x and y directions. For each face f.sub.i∈F.sub.ilm, the vertex of the RPE surface onto the XY-plane is computed, which approximates the position of the corresponding A-scan and B-scan in the volume scan. Let us consider that c.sub.i represents the centroid of the face f.sub.i. To compute the corresponding vertex in the RPE surface, the face f.sub.i is projected onto the corresponding XY-plane, {tilde over (c)}.sub.i represents the projected centroid. The terms {tilde over (c)}.sub.xi and {tilde over (c)}.sub.yi are the corresponding x and y coordinates. The x-index (i.sub.x), the y-index (i.sub.y), and the vertex (i) index are computed using xline and yline in the RPE surface for the face f.sub.i using the following equation:

[00002] $i_{x} = .Math. \frac{{\tilde{c}}_{x_{i}} - x_{\min} + ϵ_{x} / 2}{{.Math.}_{x}} .Math., i_{y} = .Math. \frac{{\tilde{c}}_{y} - y_{\min} + ϵ_{y} / 2}{ϵ_{y}} .Math., i = i_{x} + i_{y} .Math. xline,$

where ┌ . . . ┐ represents the ceil function and i denotes the corresponding vertex in the RPE surface. For an accurate computation, the neighborhood of vertex i of the RPE surface is checked and the corresponding vertex is computed as:

{tilde over (v)}.sub.i={tilde over (v)}.sub.j∈Ω.sub.i|min|{tilde over (c)}.sub.i−{tilde over (v)}.sub.j|,

where Ω.sub.i represents the 3×3 neighborhood (at XY-plane) of vertex i. The term {tilde over (v)}.sub.i corresponds to the projection of vertex v.sub.i∈V.sub.rpe onto XY-plane. Finally, we get the set

C={v.sub.i∈ custom-character .sup.3|i=0, . . . ,m.sub.ilm−1}

which represents the set of RPE surface vertices corresponding to each face in F.sub.ilm. A visual representation of the correspondence computation is shown in FIG. 6, where the RPE's regular XY-grid is shown in blue and the projected ILM's vertices and edges are painted in red.

BMO Points Computation

[0139] BMO is the termination of the Bruch's membrane (BM) layer, i.e. the lower boundary of the RPE, and serves as a stable zero reference plane for ONH quantification. Thus, the BMO is an important parameter in the detection of ONH morphologic parameters. A challenge in BMO detection is the correct identification of these points 40, especially in the presence of shadows caused by blood vessels, or the border tissue of Elsching 50—a structure similar to the BM. The BMO points 40 are segmented in the three-dimensional image data directly without the use of a two-dimensional projection image in the XY-plane, as know from the state of the art.

[0140] For BMO determination the image data are flattened. This step refers to the translation of all A-scans such that a chosen boundary in the volume is flat. In this example, the retina is aligned to the smoothed RPE.sub.lower. The alignment facilitates the volume reduction process, as well as the differentiation of BMO from other tissue.

[0141] The end-points of the rough ONH area, A.sub.ONH provide the starting points for BMO points detection. On each B-scan, the starting points are updated with new BMO points candidates if they meet the following conditions: 1) have minimum value in the 2D Morlet filtered image, 2) d(p.sub.new,p.sub.BMO-seed)<30 μm, and 3) the curvature in a neighboring region of 5 voxels is almost zero to avoid including the tissue of Elsching. In case the BMO points detected in the left and right part of one B-scan overlap, the BMO starting or end region previously defined by A.sub.ONH are updated accordingly. An example of a pair (left and right) detected BMO points is shown in FIG. 2(f).

[0142] Due to blood vessels around the ONH, noise components and three-dimensional OCT scan patterns, the BMO points are non-uniform and noisy as shown in FIG. 5(a). To remove these artifacts, an ellipse is fitted to the BMO points. Another key parameter in the ONH shape analysis is the center of the BMO points. This is computed as the barycenter of the all BMO points. FIG. 5 shows that the ellipse fitting is not only removing the noise but also increases the data points uniformly.

[0143] For the ONH shape analysis, the region inside the BMO points is of special interest since BMO points represent the optic disc margin. To segment this region, the elliptic representation of the BMO points in R.sup.2, i.e. the fitted ellipse, along with the barycenter of the BMO points are used. First, the centers of the ILM and the RPE surfaces corresponding to the center {tilde over (p)}.sub.c of the BMO points are computed.

[0144] Let us consider that {tilde over (V)}.sub.ilm and {tilde over (V)}.sub.rpe represent projected sets of vertices onto XY-plane. Then, the indices of the centers of both surfaces can be computed as:

ilmc=j∈0, . . . ,n.sub.ilm−1|min|{tilde over (v)}.sub.j.sup.ilm−{tilde over (p)}.sub.c|

rpec=j∈0, . . . ,n.sub.rpe−1|min|{tilde over (v)}.sub.j.sup.rpe−{tilde over (p)}.sub.c|

where {tilde over (v)}.sub.j.sup.ilm∈{tilde over (V)}.sub.ilm and {tilde over (v)}.sub.j.sup.rpe∈{tilde over (V)}.sub.rpe. The vertices v.sub.ilmc∈V.sub.ilm and v.sub.rpec∈V.sub.rpe represent the center of the ILM and the RPE surfaces, respectively.

[0145] To compute the BMO regions in the ILM and the RPE surfaces, the ellipse is transformed into a circle using the affine transform. Let us consider {tilde over (p)}.sub.i.sup.e represents the fitted ellipse point to {tilde over (p)}.sub.i and {tilde over (p)}.sub.i.sup.c denotes the corresponding affine transformed point on a circle of radius r. This transformation reduces the complexity of the BMO region computation and increases the speed of the method, particularly when implemented as a computer program. Now, by using the circular representation of the BMO points, the BMO regions in both RPE and ILM surfaces is computed as:

Ω.sub.ilm={f.sub.j∈F.sub.ilm∥{tilde over (c)}.sub.j.sup.ilm−{tilde over (v)}.sub.ilmc|≤r},

Ω.sub.rpe={f.sub.j∈F.sub.rpe∥{tilde over (c)}.sub.j.sup.rpe−{tilde over (v)}.sub.rpec|≤r},

where r is the radius of the transformed circle and {tilde over (c)}.sub.j.sup.ilm, {tilde over (c)}.sub.j.sup.rpe represent the centroid of a face in the ILM and the RPE surfaces, respectively. The computation of the respective regions is done using the disk growing method.

[0146] With these processing steps of the image data, the central region of the ILM polygon mesh, the central region of the RPE polygon mesh as well as the BMO opening region are determined.

Determination of Morphologic Parameters 10

ONH Cup Volume 7

[0147] The ONH cup is defined as a segment of the ILM surface 20, i.e. the ILM polygon mesh 20, particularly the central ILM region 25, which extends below the RPE surface 30 as shown in FIGS. 7(a) and (b). Note that the cup is not present in every ONH volume scan. For example, in the case of a swollen ONH, the ILM surface is always above the RPE surface. To detect the presence of the ONH cup, for each face f.sub.i∈F.sub.ilm and compute its centroid c.sub.i. As mentioned above, each face f.sub.i of the ILM has the corresponding vertex v.sub.i ∈C in RPE. If c.sub.i.sup.z−v.sub.i.sup.z≥0 for all faces in ILM, then, there is no cup available in the ONH region. Otherwise, there is a cup. The terms c.sub.i.sup.z and v.sub.i.sup.z show the corresponding z-coordinates (height). Similarly, it is possible to compute the cup region:

Ω.sub.cup={f.sub.i∈f.sub.ilm|(c.sub.i.sup.z−v.sub.i.sup.z)≤0},

where Ω.sub.cup consists of faces (triangles) of ILM, which are below the RPE surface. As it can be seen from FIG. 7(b), the cup region is also a manifold surface with proper face normal orientation. To compute the volume of the cup accurately, the face normal information at each triangle of the region can be used.

[0148] The cup volume is computed using the following formula:

[00003] $V_{c u p} = \underset{(f_{i} \in Ω_{c u p})}{.Math.} A_{i} h_{i},$

where A.sub.i represents the area of a triangle, which is a projection of the face f.sub.i of the ILM surface onto the XY-plane and h.sub.i is the height with respect to the RPE surface. These variables are defined as:

A.sub.i=½({right arrow over (e)}.sub.0×{right arrow over (e)}.sub.1),h.sub.i=(c.sub.i.sup.z−v.sub.i.sup.z),

where {right arrow over (e)}.sub.0 and {right arrow over (e)}.sub.1 are the connected edges of the projected triangle. The cross product between the two edges will take care of the orientation of the corresponding face and enables a precise volume computation even in complex topological regions.

Central ONH Thickness (CONHT):

[0149] The CONHT is defined as the height difference between the center of the ILM and the RPE surfaces as shown in FIG. 7(c).

BMO Region Volume:

[0150] The BMO region volume is computed using the segmented ILM and RPE surfaces. For this reason the cup volume is separated from the BMO region volume such that it does not include the cup volume, if it exists. Then, the BMO region volume can be computed as:

Ω.sub.bmo=Ω.sub.ilm/Ω.sub.cup.

[0151] Similar to the ONH cup volume, the BMO region volume is determined using the similar formula:

[00004] $V_{bmo} = \underset{(f_{i} \in Ω_{bmo})}{.Math.} A_{i} h_{i},$

where A.sub.i is area of the face f, which belongs to the set Ω.sub.bmo.

[0152] Similarly, h.sub.i is also computed using the corresponding vertices in the RPE surface.

ONH Total Volume:

[0153] Similar to the ONH cup and the BMO region volumes, the ONH total volume is also computed using ILM and RPE surfaces. The total volume is computed from the circular region, with radius 1.5 mm, centered at v.sub.ilmc and v.sub.rpec for ILM and RPE surfaces respectively, as shown in FIG. 9(a). The circular regions for ILM and RPE surfaces are computed using the following equation:

Ω.sub.1.5mm.sup.ilm={f.sub.j∈F.sub.ilm∥{tilde over (c)}.sub.j.sup.ilm−{tilde over (v)}.sub.ilmc|≤1.5 mm},

Ω.sub.1.5mm.sup.rpe={f.sub.j∈F.sub.rpe∥{tilde over (c)}.sub.j.sup.rpe−{tilde over (v)}.sub.rpec|≤1.5 mm},

where Ω.sub.1.5mm.sup.ilm and Ω.sub.1.5mm.sup.rpe are the sets consisting of all faces within the 1.5 mm region of the ILM and RPE surfaces from their centers. Then, the total volume region Ω.sub.tv is calculated using:

Ω.sub.tv=Ω.sub.1.5 mm\Ω.sub.cup

where represents Ω.sub.tv the total volume region on the ILM surface as shown in FIG. 9(b), which is further employed for the total volume computation:

[00005] $V_{tv} = \underset{(f_{i} \in Ω_{t ν})}{.Math.} A_{i} h_{i},$

where A.sub.i is the area of the face f.sub.i ∈Ω.sub.tv and h.sub.i is the height w.r.t. corresponding vertex in the RPE surface.

ONH Annular Region Volume:

[0154] The ONH annular region represents the ONH outer region, see FIG. 9(c). On ILM surface this region is computed using the following equation:

Ω.sub.ab=Ω.sub.1.5 mm\Ω.sub.bmo

where Ω.sub.av consists of all the ILM surface faces which belong to the annular region of the ONH. The volume of the annular region is computed using the ILM and RPE surfaces correspondence.

[00006] $V_{av} = \underset{(f_{i} \in Ω_{av})}{.Math.} A_{i} h_{i} .$

[0155] The annular region volume helps to see the change in the outer region of the ONH volume in different cohorts.

Bending Energy 6:

[0156] The roughness on the ILM surface within the BMO region is an important parameter and commonly known as the bending energy on a manifold surface. The bending energy measures the fairness of a surface in terms of the curvature. In general, the outer region of the ILM surface is quite smooth and flat unlike the one inside the BMO, which has very complex topological structure. In this paper, we define the bending energy within the BMO region using the element-based normal voting tensor (ENVT). The ENVT exploits the orientation information (face normals) to compute a shape analysis operator at each face f.sub.i∈Ω.sub.bmo and is defined as:

[00007] $M_{i} = \frac{1}{{.Math.}_{f_{j} \in Ω_{bmo}} a_{j}} \underset{f_{j} \in Ω_{b m o}}{.Math.} a_{j} n_{j} .Math. n_{j}^{T},$

where n.sub.j represents the normal of face f.sub.i and n.sub.j.sup.T the transpose of n.sub.j. The term a.sub.j is the area of the face f.sub.i. To assure robustness against irregular sampling of the ILM surface, the above equation is weighted by the corresponding face area a.sub.j. The ENVT, M.sub.i is a symmetric and positive semidefinite matrix, so, it can be decomposed into its spectral components:

[00008] $M_{i} = {.Math.}_{k = 0}^{2} λ_{k}^{i} e_{k} e_{k}^{T},$

where λ.sub.k.sup.i={λ.sub.0.sup.i,λ.sub.1.sup.i,λ.sub.2.sup.i} are the eigenvalues of the eigenvectors and these eigenvalues are sorted in decreasing order (λ.sub.0.sup.i>λ.sub.1.sup.i>λ.sub.2.sup.i≥0). The corresponding eigenvector is denoted by e.sub.k. In general, the dominant eigenvalue λ.sub.0.sup.i has the corresponding eigenvector in the direction of the face normal and the remaining two eigenvectors will be aligned to the principle curvature direction on the ILM surface. On the planar region, only λ.sub.0.sup.i will be significant, on the edge region, λ.sub.0.sup.i and λ.sub.1.sup.i will be significant and at the corners, all of these eigenvalues are significant. Using the anisotropic properties of these eigenvalues, we define the bending energy inside the BMO region using the following Equation:

[00009] $E_{b} = \underset{f_{i} \in Ω_{bmo}}{.Math.} λ_{2}^{i} + λ_{3}^{i},$

where λ.sub.1.sup.i and λ.sub.2.sup.i are the two least dominant eigenvalues of the ENVT of the face f.sub.i. FIG. 7(d) shows how each face of the BMO region is colored based on the bending energy. The color is scaled from white (flat regions) to dark (sharp features).

BMO-MRW:

[0157] BMO-MRW, has been proposed by [6] as a valid alternative structural measure. It computes the minimum distance between the BMO points and the ILM surface. The average BMO-MRW, denoted by avg.sub.ni, is calculated as:

[00010] $a υ g_{mrw} = \frac{1}{n_{p}} {.Math.}_{i = 1}^{n_{p}} .Math. p_{i}^{mrw} - p_{i}^{e} .Math. .$

[0158] BMO-MRW surface area:

[0159] BMO-MRW surface area, BMO-MRA, is computed by taking the whole region defined by all BMO-MRW. The ellipse fitted BMO points P.sub.2D and the z-coordinates from P are combined and are represented as: P.sub.e={p.sub.i.sup.e∈R.sup.3|i=0, . . . , n.sub.p−1}.

[0160] For each point p.sub.i.sup.e∈P, a point p.sub.i.sup.mrw on the ILM surface is computed:

p.sub.i.sup.mrw=v.sub.j∈V.sub.ilm|min|v.sub.j−p.sub.i.sup.e|.

[0161] The MRW points P.sub.mrw={p.sub.i∈R.sup.3|i=0, . . . , n.sub.p−1} are lying on the ILM surface. A A quad surface is created using point sets P.sub.e and P.sub.mrw by introducing edges between the corresponding vertices in both point sets and connecting the neighbor points. The number of quad elements in the MRD surface is equal to the number of points in each point sets and is represented as: Q={q.sub.i|i=0, . . . , n.sub.p−1}. MRW-MRA is computed as:

[00011] $A_{mrw} = \underset{q_{i} \in Q}{.Math.} a_{q_{i}},$

where a.sub.qi represents the area of the quad q.sub.i.

BMO Area:

[0162] As shown in FIG. 5, BMO area (BMOA) represents the area under the fitted ellipse to the BMO points and is computed using the conic representation:

A.sub.bmo=πr.sub.1r.sub.2

where r.sub.1 and r.sub.2 are the major and minor axes of the fitted ellipse.

Experiment and Results

[0163] In order to evaluate the method according to the invention repeated-measurement reliability tests were performed, the ONH shape in healthy subjects was investigated, and tested in order determine if the method is able to detect differences in patients with diseases known to affect the ONH in form of swelling and atrophy.

[0164] In order to estimate the repeated-measurement reliability three repeated scans of each eye from ten healthy subjects were taken. These subjects were measured each in a time frame of a week, and then again in the following week. Table 1 shows the repeatability results. The method according to the invention scores highly at every parameter presented, with lowest intraclass correlation coefficient (ICC) of 0.905 for CONHT, and highest 0.998 for V.sub.cup. The ICC and confidence intervals were estimated using the variance components from a one-way ANOVA.

[0165] The method according to the invention, was also tested with several other scan protocols of the same device (ONH cube with 73 B-scans, scanning angle of 15°×15° and resolution 384 A-scans per B-scan, spatial resolution in x direction is ≈12.6 μm, in axial direction ≈3.9 μm and the distance between two B-scans≈61 μm, ONH star scan with 24 B-scans, scanning angle of 15°×15° and a resolution of 768 A-scans per B-scan, spatial resolution in x direction is ≈5.36 μm, in axial direction ≈3.9 μm) and the volumetric ONH-centered protocol acquired using Cirrus HD OCT (Carl Zeiss Meditec, Dublin, Calif.), which covers 6×6×2 mm.sup.3 region with 200×200×1024 voxels and obtained positive results.

[0166] BMO detection was validated and the RPE segmentation was checked. Five scans from the ones used in the repeatability testing were randomly selected. An experienced grader manually selected the BMO points. This resulted in a total amount of 488 B-scans with manually selected BMO points which corresponded to the number detected automatically.

TABLE-US-00001 TABLE 1 Repeatability test for the 3D parameters. Parameters ICC LCI UCI ONH Cup Vol.(V.sub.cup) (mm.sup.3) 0.998 0.996 0.999 Central ONH Thickness (CONHT) (mm) 0.905 0.813 0.958 BMO Region Vol (V.sub.bmo) (mm.sup.3) 0.993 0.986 0.997 ONH Total Vol. (V.sub.tv) (mm.sup.3) 0.995 0.989 0.998 ONH Annular Vol. (V.sub.av) (mm.sup.3) 0.983 0.965 0.993 Bending Energy (E.sub.b) 0.911 0.824 0.961 BMO-MRW Surface Area (A.sub.mrw) (mm.sup.2) 0.910 0.823 0.960 BMO-MRW (d.sub.mrw) (mm) 0.993 0.986 0.997 BMO Area (A.sub.bmo) (mm.sup.2) 0.989 0.976 0.995 Abbreviations: ICC - intra-class correlation coefficient, LCI - lower boundary of 95% confidence interval and UCI - upper boundary of 95% confidence interval.

[0167] Furthermore, the mean signed and unsigned error in the x axis was compared, as well as in the axial one (in z axis). If the automated method identified the BMO closer to the optic disc center, the sign of distance in the x-direction was positive. Similarly, if the automated BMO located below the manual BMO, the sign of distance in the z-direction was positive. Results are shown in Table 2.

TABLE-US-00002 TABLE 2 Mean (±SD) Mean (±SD) Mean (±SD) Mean (±SD) unsigned error (pixels) unsigned error (μm) signed error (pixels) signed error (μm) x axis 4.9098 (±4.9182) 61.8635 (±61.9693) −0.6107 (±6.9261) −7.6943 (±87.2693) z axis 2.8828 (±3.1801) 12.4024 (±11.2429) −0.3618 (±4.2790) −1.4110 (±16.6881) Mean unsigned and signed error in pixel and μm, for the x axis, and z axis between automatic (proposed) and manual segmentation.

Clinical Evaluation

[0168] In this section, the results of the method according to the invention for 248 OCT scans are presented, from three groups, 71 healthy control eyes (HC), 31 eyes of patients suffering from idiopathic intracranial hypertension (IIH), which causes swelling of the ONH (papilledema). We also included 146 eyes of patients with autoimmune central nervous system disorders (multiple sclerosis (MS) and nueromyelitis spectrum disorder (NMOSD)) and a history of optic neuritis (ON), an inflammatory optic neuropathy that damages the optic nerve leading to neuroaxonal degeneration.

[0169] In IIH patients, the ONH volume is increased and was shown to correlate with cerebrospinal fluid (CSF) pressure. The longitudinal analyses revealed that ONH volume measured by OCT decreased after the initial lumbar puncture and initiation of therapy with acetazolamide. Additionally, increased ONH volume was associated with lower visual acuity in IIH patients, which points out to the potential clinical relevance of the parameter.

[0170] ON is one of the most common initial clinical presentations of MS without any prior history of a demyelinating event. During the course of the disease, acute ON affects 50%, 70% of MS patients. After initial swelling due to edema in the acute phase of ON, retinal nerve fiber layer (RNFL) thickness decreased over the following 6 months. Optic neuritis (ON) is the first NMOSD-related clinical event in 55% of the patients, which causes severe structural damage to the optic nerve and retina with resulting functional impairment. Recurrent ONs in NMOSD give rise to severely thinned pRNFL and combined ganglion cell layer and inner plexiform layer (GCIP).

[0171] The results presented in Table 3 demonstrate that the method according to the invention successfully captures the differences. The only parameter showing no difference between groups is bending energy, E.sub.b. Although it was expected that in the IIH affected eyes the ONH shape inside the BMO to have a smoother convex shape, the data is still extremely heterogeneous. Thus, the bending energy reflects this extreme variability in the data.

[0172] Especially inside the BMO the topologies from one ONH to the other can be extremely different.

TABLE-US-00003 TABLE 3 HC IIH ON GEE Parameters Mean (SD) Min-Max Mean (SD) Min-Max Mean (SD) Min-Max p(HC-IIH) p(HC-ON) p(IIH-ON) ONH Cup 0.044 (0.075) 0.000-0.370 0.004 (0.012) 0.000-0.049 0.035 (0.052) 0.000-0.256 0.0012 0.353 1.95E−07 Vol.(V.sub.cup) (mm.sup.3)) Central ONH 0.002 (0.25) −0.545-0.550 0.403 (0.340) −0.242-0.966 0.012 (0.231) −0.399-0.668 9.91E−06 0.747 4.22E−06 Thickness (CONHT) (mm) BMO 0.572 (0.21) 0.107-1.035 1.521 (1.022) 0.436-4.3098 0.444 (0.189) 0.105-1.050 5.02E−05 0.0012 2.66E−06 Region Vol (V.sub.bmo) (mm.sup.3) ONH 2.523 (0.250) 2.114-2.988 3.309 (0.799) 2.351-5.365 2.163 (0.295) 1.522-3.087 2.57E−05 5.28E−13 6.54E−10 Total Vol. (V.sub.tv) (mm.sup.3) ONH 1.906 (0.181) 1.504-2.268 1.78 (0.357) 0.572-2.200 1.683 (0.210) 1.162-2.118 0.146 2.19E−10 0.257 Annular Vol. (V.sub.av) (mm.sup.3) Bending 0.026 (0.009) 0.011-0.065 0.029 (0.013) 0.008-0.066 0.029 (0.010) 0.011-0.066 0.477 0.038 0.793 Energy (E.sub.b) BMO-MRW 1.791 (0.459) 0.763-2.789 2.975 (1.016) 1.439-5.351 1.414 (0.364) 0.615-2.765 1.35E−06 5.44E−06 3.94E−11 Surface Area (A.sub.mrw) (mm.sup.2) BMO-MRW 0.135 (0.054) 0.261-0.035 0.263 (0.105) 0.127-0.447 0.093 (0.040) 0.030-0.260 3.52E−07 2.69E−05 1.81E−12 (d.sub.mrw) BMO Area 1.906 (0.352) 1.158-2.646 2.682 (1.109) 1.057-5.709 1.861 (0.377) 1.030-3.326 0.0028 0.407 0.001 (A.sub.bmo) (mm.sup.2) Table 3: Analysis of all the 3D parameters defined for the HC and patient group. The last column shows the GEE analysis between the two groups. Abbreviations: HC—healthy controls. SD—standard deviation, Min—minimum value, Max—maximum value, GEE—generalized estimating equation models analysis accounting for the inter-eye/intra-subject dependencies, p—p value.

[0173] In FIG. 8 it is exemplary shown in which region the Bruch's membrane can be seen in a OCT scan. Close to the optic nerve head, the RPE portion 3 ends and the Bruch's membrane 44 (underlined by a white bar) is visible as a continuation of the RPE 3. At the end points of the Bruch's membrane 44 the Bruch's membrane openings 43 are visible that can be identified as well by the method according to the invention.

[0174] Moreover, the border tissue of Elsching 50 is also visible in this OCT scan. The vitreous 51 is also indicated in this OCT-scan.

[0175] With the proposed method that can be implemented in a computer for execution, it is possible to efficiently and accurately determine several morphologic parameters, particularly a bending energy that provides a novel morphologic parameter characterizing the optic nerve head.

REFERENCES

[0176] [1] E. M. Kadas, F. Kaufhold, C. Schulz, et al., 3D Optic Nerve Head Segmentation in Idiopathic Intracranial Hypertension, 262-267. Springer Berlin Heidelberg, Berlin, Heidelberg (2012). [0177] [2] C.-L. Chen, H. Ishikawa, G. Wollstein, et al., “Individual a-scan signal normalization between two spectral domain optical coherence tomography devices,” Investigative Ophthalmology & Visual Science 54(5), 3463-3471 (2013). [0178] [3] M. K. Garvin, M. D. Abramoff, R. Kardon, et al., “Intraretinal layer segmentation of macular optical coherence tomography images using optimal 3d graph search,” IEEE Transactions on Medical Imaging 27, 1495-1505 (2008). [0179] [4] B. J. Antony, A Combined Machine-learning and Graph-based Framework for the 3-d Automated Segmentation of Retinal Structures in Sd-oct Images. PhD thesis, Iowa City, Iowa, USA (2013). AAI3608177. [0180] [5] K. Rohr, H. S. Stiehl, R. Sprengel, et al., Point-based elastic registration of medical image data using approximating thin-plate splines, 297-306. Springer Berlin Heidelberg, Berlin, Heidelberg (1996). [0181] [6] A. S. C. Reis, N. O'Leary, H. Yang, et al., “Influence of clinically invisible, but optical coherence tomography detected, optic disc margin anatomy on neuroretinal rim evaluation,” Investigative Ophthalmology & Visual Science 53(4), 1852-1860 (2012).

METHOD FOR AUTOMATIC SHAPE QUANTIFICATION OF AN OPTIC NERVE HEAD

Assignee

Inventors

Cpc classification

Classification Explorer

A61B3/0058

HUMAN NECESSITIES

Classification Explorer

G06T7/62

PHYSICS

Classification Explorer

G06T7/0012

PHYSICS

Classification Explorer

A61B3/0025

HUMAN NECESSITIES

Classification Explorer

G06T2207/10101

PHYSICS

Classification Explorer

G06T17/20

PHYSICS

Classification Explorer

G06T7/155

PHYSICS

Classification Explorer

G06T7/12

PHYSICS

Classification Explorer

G06T2207/30041

PHYSICS

Classification Explorer

A61B3/1225

HUMAN NECESSITIES

Classification Explorer

A61B3/102

HUMAN NECESSITIES

Classification Explorer

G06T2210/41

PHYSICS

International classification

Classification Explorer

A61B3/00

HUMAN NECESSITIES

Classification Explorer

A61B3/10

HUMAN NECESSITIES

Classification Explorer

A61B3/12

HUMAN NECESSITIES

Classification Explorer

G06T17/20

PHYSICS

Classification Explorer

G06T7/155

PHYSICS

Abstract

Claims

Description