Processing techniques for optical retinal imaging

Abstract

The invention relates to processing techniques for optical retinal imaging. A computer-implemented method for processing respective phase components of a first OCT image and a second OCT image of a series of OCT images of a common portion of the retina acquired by a Fourier domain OCT imaging system after stimulation of the common portion by optical stimulation, the common portion includes a layer of the retina whose thickness varies during acquisition of a series of OCT images to determine an indication of a position of a boundary of the layer in an axial direction in the OCT images, the method comprises: processing a phase component of the first OCT image and a phase component of the second OCT image to calculate a velocity profile indicative of a distribution of velocities in an axial direction within the common portion of the retina; and determining the indication based on the calculated velocity profile.

Description

field

[0001] Example aspects of the present invention generally relate to the field of optical coherence tomography (OCT), and specifically to techniques for processing OCT data generated by a Fourier domain OCT imaging system to generate optical retinal imaging (ORG) data that is indicative of the physiological response of the retina of a subject's eye to optical stimulation. background

[0002] Optical coherence tomography (OCT) is an imaging technique based on low-coherence interferometry that is widely used to acquire high-resolution two-dimensional and three-dimensional images of optically scattering media such as biological tissue.

[0003] Depending on how depth ranging is achieved, OCT imaging systems can be classified as time-domain OCT (TD-OCT) or Fourier domain OCT (FD-OCT) (also known as frequency domain OCT). In TD-OCT, the optical path length of the reference arm of the interferometer of the imaging system varies with time during the acquisition of the reflectivity profile of the scattering medium (referred to herein as the "imaging target") imaged by the OCT imaging system, and the reflectivity profile is often referred to as a "depth scan" or "axial scan" ("A scan"). In FD-OCT, the spectral interferogram produced by the interference between the light in the reference arm of the interferometer and the light in the sample arm at each A scan position is Fourier transformed to simultaneously acquire all points along the depth of the A scan without any change in the optical path length of the reference arm. FD-OCT can allow much faster imaging than scanning of the sample arm mirror in the interferometer because all back reflections from the sample are measured simultaneously. Two common types of FD-OCT are spectral domain OCT (SD-OCT) and swept source OCT (SS-OCT). In SD-OCT, a broadband light source delivers many wavelengths to the imaging target, and a spectrometer is used as a detector to measure all wavelengths simultaneously. In SS-OCT (also known as time-coded frequency-domain OCT), the light source is swept across a range of wavelengths, and the temporal output of the detector is converted into spectral interferometry.

[0004] OCT imaging systems can also be classified as point scanning (also called "point detection" or "scanning point"), line scanning, or full-field, depending on how the imaging system is configured to acquire OCT data at various locations on an imaging target. A point scanning OCT imaging system acquires OCT data by scanning a focused sample beam across the surface of an imaging target, typically along a single line (which may be straight, or alternatively curved to define a circle or spiral) or a set of (usually substantially parallel) lines on the imaging target surface, and acquiring an axial depth profile (A-scan) for each of a plurality of points along the line, one point at a time, to create OCT data comprising a one-dimensional or two-dimensional array of A-scans representing a two-dimensional (i.e., B-scan) or three-dimensional (i.e., C-scan or volume scan) reflectance profile of the sample.

[0005] Line scan OCT imaging systems acquire OCT data by scanning with focused light across the surface of an imaging target. The measured reflectance from the imaging target is used to generate OCT data comprising a two-dimensional reflectance distribution (i.e., a B-scan) of the sample. By scanning with focused light across multiple locations on the imaging target, OCT data comprising a three-dimensional reflectance distribution (i.e., a C-scan or volume scan) of the sample can be obtained. Typically, the focused light is straight and is scanned in a direction perpendicular to it, although in some cases it can be curved and the scanning direction adjusted accordingly. Full-field OCT imaging systems acquire OCT data by projecting a light beam onto the imaging target to acquire OCT data comprising a three-dimensional reflectance distribution (i.e., a C-scan or volume scan) of the sample.

[0006] OCT imaging systems can also be classified as phase-resolved, where both the intensity and phase of light reflected from the imaging target are measured as a function of axial depth. Modern FD-OCT imaging systems typically possess a degree of phase stability, which allows them to be used as OCT imaging systems with phase-resolved capabilities.

[0007] Optical retinal imaging (ORG) generally refers to detecting the physiological response of the retina of the eye to optical stimulation (i.e., the functional activity of the retina induced by light). ORG technology includes non-invasive optical imaging of this physiological response of the retina. For example, an OCT imaging system can be used to image retinal neurons that exhibit size (size) changes in response to excitation of an optical stimulus. These size changes (usually changes in the length of the outer segment (OS) of the photoreceptors in the retina, which is the depth difference between the junction of the inner and outer segments (IS / OS) of the cone photoreceptors and the outer segment tip (COST) of the cone) cause phase changes in the light waves returned from the eye that are large enough to be detected by an OCT imaging system with phase resolution capability. Phase is very sensitive to movement in tissue, and many OCT imaging systems with phase resolution capability are able to resolve displacements of less than 10 nm, which may be much smaller than the axial resolution of the system or the wavelength of light used in imaging. Overview

[0008] According to a first exemplary aspect of the present disclosure, a computer-implemented method is provided for processing phase components (information) of a first OCT image and a phase component of a second OCT image of a series of OCT images of a common portion of a retina of an eye acquired by a Fourier-domain optical coherence tomography (OCT) imaging system after the common portion of the retina is stimulated by an optical stimulus, wherein the common portion includes a layer of the retina whose thickness changes in response to the optical stimulus during acquisition of the series of OCT images (such as the outer segments (OS) of photoreceptor cells in the retina) to determine at least one of a first indication of a first location of a first boundary of the layer along an axial (z-axis) direction in the OCT image and a second indication of a second location of a different second boundary of the layer along the axial direction in the OCT image. The method includes: processing the phase components of the first OCT image and the phase components of the second OCT image to calculate a velocity profile indicating a distribution of velocity within the common portion of the retina along the axial direction; and determining at least one of the first indication and the second indication based on the calculated velocity profile.

[0009] At least one of the first indication and the second indication can be determined based on the calculated velocity profile by: determining, in the case of determining the first indication, a first position in the velocity profile corresponding to a maximum value of the velocity indicated by the velocity profile as the first indication; and determining, in the case of determining the second indication, a second position in the velocity profile corresponding to a minimum value of the velocity indicated by the velocity profile as the second indication. Thus, the velocity in a portion of the velocity profile corresponding to a layer in the retina can vary from a maximum velocity at a first position in the velocity profile to a minimum velocity at a second position in the velocity profile, the first position in the velocity profile corresponding to a first boundary of the layer at a location on the retina, the second position in the velocity profile corresponding to a second boundary of the layer at the location on the retina, and the first indication can indicate the maximum velocity at the first position in the velocity profile, and the second indication can indicate the minimum velocity at the second position in the velocity profile.

[0010] At least one of the first indication and the second indication may be determined by processing the calculated velocity profile using a cluster analysis algorithm or a dimensionality reduction (e.g., component analysis) algorithm. For example, the dimensionality reduction algorithm may include one of a principal component analysis (PCA) algorithm, an independent component analysis (ICA) algorithm, a linear discriminant analysis (LDA) algorithm, or a non-negative matrix factorization (NMF) algorithm. In the case where the dimensionality reduction algorithm includes a PCA algorithm, at least one of the first indication and the second indication may be determined using principal components determined by applying the PCA algorithm to the calculated velocity profile.

[0011] Both the first position in the speed distribution map and the second position in the speed distribution map can be determined in the following manner: for each combination of the candidate first position in the speed distribution map and the candidate second position in the speed distribution map among all combinations of the candidate first position and the candidate second position in the speed distribution map, calculating a corresponding value of the difference between the corresponding speed for the candidate first position indicated by the speed distribution map and the corresponding speed for the candidate second position indicated by the speed distribution map in the combination; and identifying the candidate first position and the candidate second position of the following combination among multiple combinations as the first position in the speed distribution map and the second position in the speed distribution map: for the combination among the multiple combinations, the calculated value of the difference is the largest among multiple calculated values ​​of the difference.

[0012] Alternatively, both the first position in the speed profile and the second position in the speed profile are determined by: calculating, for each combination of the candidate first position in the speed profile and the candidate second position in the speed profile among all combinations of the candidate first position and the candidate second position in the speed profile, a corresponding value of the difference between an average value of corresponding speeds for a set of neighboring positions including the candidate first position in the combination and an average value of corresponding speeds for the set of neighboring positions including the candidate second position in the combination indicated by the speed profile; and identifying, as the first position in the speed profile and the second position in the speed profile, the candidate first position and the candidate second position of the following combination among a plurality of combinations: for the combination among the plurality of combinations, the calculated value of the difference is the largest among the plurality of calculated values ​​of the difference.

[0013] The computer-implemented method of the above-mentioned first example aspect or any example implementation thereof may also include one of the following: determining a first indication of a first position of a first boundary of the layer in an axial direction in an OCT image based on an amplitude component of at least one OCT image in an OCT image of a series of OCT images, wherein the second indication is determined by identifying a second position in the velocity distribution map corresponding to a minimum value of the velocity indicated by the velocity distribution map; and determining a second indication of a second position of a second boundary of the layer in the axial direction in the OCT image based on an amplitude component of at least one OCT image in an OCT image of a series of OCT images, wherein the first indication is determined by identifying the first position in the velocity distribution map corresponding to a maximum value of the velocity indicated by the velocity distribution map.

[0014] In some example embodiments, the computer-implemented method may further include calculating a comparison value, which is (i) a value of an image quality metric calculated based on at least one of an amplitude component of the first OCT image and an amplitude component of the second OCT image, or (ii) a value of an image similarity metric that provides a measure of the degree of similarity between the images, the image similarity metric calculated based on the first OCT image and the second OCT image. The computer-implemented method may further include comparing the comparison value with a threshold to determine whether the comparison value is equal to or greater than the threshold. If the comparison value is determined to be equal to or greater than the threshold, a first indicator may be generated indicating that the determined at least one of the first indication and the second indication is reliable. If the comparison value is determined not to be equal to or greater than the threshold, a second indicator may be generated indicating that the determined at least one of the first indication and the second indication is unreliable. The comparison value may be a maximum value of the calculated cross-correlation between the first OCT image and the second OCT image. For example, other similarity metrics that may be used alternatively include sum of squared differences, mutual information, normalized mutual information, and Kullback Leibler distance.

[0015] In some other example embodiments, the computer-implemented method may further include: for each of a plurality of different sets of two or more OCT images in a series of OCT images, calculating a corresponding comparison value, the comparison value being (i) a corresponding value of an image quality metric calculated based on an amplitude component of at least one OCT image in the set, or (ii) a corresponding value of an image similarity metric that provides a measure of similarity between the images, the value of the image similarity metric being calculated based on the amplitude components of at least two OCT images in the set. The method may further include comparing each comparison value with a threshold to determine whether the comparison value is equal to or greater than the threshold. For each of the plurality of different sets of OCT images, for sets where the calculated comparison value is determined to be equal to or greater than the threshold, processing the phase components of the OCT images in the set to generate a corresponding velocity profile indicating a distribution. Furthermore, for each of the plurality of different sets of OCT images, for sets where the calculated comparison value is determined to be not greater than the threshold, generating a corresponding velocity profile indicating zero velocity at all locations along the axial direction in the velocity profile. A concatenation of the calculated velocity profiles is then generated such that the concatenation indicates how the profile changes over time, and corresponding portions of the concatenation of velocity profiles are integrated to generate data indicating how the optical path length changes over time at positions along the axial direction, the portions having the same position along the axial direction (e.g., one of the first position and the second position indicated by the generated first or second indication). The generated data may include one or more sets of equal consecutive values, and the method may further include smoothing the generated data by replacing one or more values ​​in at least one of the one or more sets of equal consecutive values ​​with one or more estimated values, the one or more estimated values ​​being calculated based on neighboring values ​​adjacent to one of the sets of equal consecutive values ​​in the generated data.

[0016] In yet other example embodiments, the computer-implemented method may further include, for each of a plurality of different sets of two or more OCT images in the series of OCT images, calculating a corresponding comparison value, the comparison value being (i) a corresponding value of an image quality metric calculated based on an amplitude component of at least one OCT image in the set, or (ii) a corresponding value of an image similarity metric that provides a measure of similarity between the images, the value of the image similarity metric being calculated based on the amplitude components of at least two OCT images in the set. In these example embodiments, the method further includes comparing each comparison value to a threshold to determine whether the comparison value is equal to or greater than the threshold. Some of the comparison values ​​may be determined to be greater than or equal to the threshold, while other comparison values ​​may be determined to be less than the threshold. For each of the plurality of different sets of OCT images, for each set for which the calculated comparison value is determined to be equal to or greater than the threshold, the phase components of the OCT images in the set are processed to generate a corresponding velocity profile indicating a distribution. For each of the plurality of different sets of OCT images, for a set in which the calculated comparison value is determined not to be greater than or equal to the threshold value, a corresponding velocity profile indicating a distribution is generated (estimated) based on a corresponding velocity profile calculated for each of one or more other sets of the plurality of different sets of OCT images. A cascade of the generated velocity profiles is then generated such that the cascade of velocity profiles indicates how the distribution changes over time. Corresponding portions of the cascade of velocity profiles (the portions having the same position along the axial direction) are integrated to generate data indicating a change in the optical path length at the position along the axial direction over time.

[0017] In further example embodiments, the computer-implemented method may further include calculating, for each of a plurality of different sets of two or more OCT images in the series of OCT images, a corresponding comparison value, the comparison value being (i) a corresponding value of an image quality metric calculated based on an amplitude component of at least one OCT image in the set, or (ii) a corresponding value of an image similarity metric that provides a measure of similarity between the images, the value of the image similarity metric being calculated based on the amplitude components of at least two OCT images in the set. In these example embodiments, each comparison value is compared to a threshold to determine whether the comparison value is equal to or greater than the threshold. Some of the comparison values ​​may be determined to be greater than or equal to the threshold, while other comparison values ​​may be determined to be less than the threshold. For each of the plurality of different sets of OCT images, for each set for which the calculated comparison value is determined to be equal to or greater than the threshold, the phase components of the OCT images in the set are processed to generate a corresponding velocity profile indicating a distribution. For each of the plurality of different sets of OCT images, for a set in which the calculated comparison value is determined not to be greater than or equal to the threshold value, a corresponding velocity profile indicating the distribution is generated based on the corresponding velocity profiles calculated for each of one or more other sets of the plurality of different sets of OCT images. A cascade of the generated velocity profiles is then generated such that the cascade of velocity profiles indicates how the distribution changes over time, and corresponding portions of the cascade of velocity profiles (the portions having the same position along the axial direction) are integrated to generate data indicating how the optical path length at the position along the axial direction changes over time.

[0018] The respective comparison value calculated for each set may be the respective maximum value of the calculated cross-correlations between at least two OCT images in the set. For example, other similarity metrics that may be used instead are sum of squared differences, mutual information, and normalized mutual information.

[0019] In an example embodiment of generating a cascade, a corresponding indication of the change in optical path length over time at each of at least one of a first position along the axial direction and a second position along the axial direction indicated by at least one of the first indication and the second indication can be determined in the following manner, namely, by integrating the corresponding part of the cascade at at least one of the first position along the axial direction and the second position along the axial direction.

[0020] In any of the foregoing, the layer of the retina may include the outer segments of photoreceptor cells.

[0021] In addition, in any of the foregoing cases, the computer-implemented method may also include processing the phase component of the first OCT image and the phase component of the second OCT image before calculating the velocity distribution map to compensate for the bulk motion of the common portion during the acquisition of a series of OCT images by the Fourier domain OCT imaging system after the common portion of the retina is stimulated by optical stimulation.

[0022] According to a second exemplary aspect of the present invention, there is also provided a computer program comprising computer-readable instructions that, when executed by a processor, cause the processor to perform the method of the first exemplary aspect or any exemplary implementation and embodiment thereof set forth above. The computer program may be stored on a non-transitory computer-readable storage medium (e.g., a computer hard disk or CD) or may be carried by a computer-readable signal.

[0023] According to a third exemplary aspect of the present disclosure, a data processing apparatus is further provided. The apparatus is configured to process a phase component of a first OCT image and a phase component of a second OCT image of a series of OCT images of a common portion of the retina of an eye acquired by a Fourier-domain optical coherence tomography (OCT) imaging system after the common portion of the retina is stimulated by optical stimulation, wherein the common portion includes a layer of the retina whose thickness changes in response to the optical stimulation during acquisition of the series of OCT images, to determine at least one of a first indication of a first position of a first boundary of the layer along an axial direction in the OCT image and a second indication of a second position of a second boundary of the layer along an axial direction in the OCT image. The data processing apparatus is arranged to perform the method of the first exemplary aspect or any exemplary implementation and embodiment thereof set forth above. The data processing apparatus may include a processor and a storage device storing computer-readable instructions that, when executed by the processor, cause the processor to perform the method of the first exemplary aspect or any exemplary implementation and embodiment thereof set forth above. Alternatively, the data processing apparatus may be implemented in non-programmable hardware, such as an ASIC, FPGA, or other integrated circuit configured to perform any of these methods.

[0024] According to a fourth exemplary aspect of the present disclosure, a Fourier domain OCT imaging system is further provided, which includes the data processing apparatus of the third exemplary aspect set forth above. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] Example embodiments will now be explained in detail, by way of non-limiting example only, with reference to the accompanying drawings described below.Unless otherwise indicated, like reference numerals appearing in different figures of the drawings may represent the same element or a functionally similar element.

[0026] Figure 1is a schematic diagram of a system including a Fourier domain OCT imaging system 30 and a data processing apparatus 100 according to an example embodiment herein.

[0027] Figure 2 is a schematic diagram of an example implementation of a data processing apparatus 100 of an example embodiment in programmable signal processing hardware.

[0028] Figure 3 is a flow chart illustrating a method by which phase components of OCT images from a series of OCT images of a common portion of the retina acquired after optical stimulation of the common portion are processed to identify boundaries of layers, according to an example embodiment.

[0029] Figure 4A Shown is the sinusoidal bulk oscillation added to the oscillatory movement of the two model retinal layers.

[0030] Figure 4B Shown is the effect of removing global motion on the temporal variation of the phase angle of light reflected from two model retinal layers.

[0031] Figure 5 The ORG contributors for two model retinal layers are shown separate and non-overlapping from each other.

[0032] Figure 6 Shown are the magnitudes of the model retinal ORG responses without providing compensation for global motion.

[0033] Figure 7A Shown are the phases of the model retinal ORG responses without providing compensation for global motion.

[0034] Figure 7B The phase of the ORG response is shown for the case where the overall motion has been compensated.

[0035] Figure 8 shows that by inserting MATLAB TM The PCA algorithm in is applied to the PCA components obtained from the ORG responses calculated for the case where the overall motion has been compensated.

[0036] Figure 9 The absolute values ​​of the principal component coefficients returned by the PCA algorithm are shown.

[0037] Figure 10 Shown are the ORG contributors of two retinal layers in a variation in which the ORG contributors of the two retinal layers partially overlap with each other.

[0038] Figure 11shows the case where the overall motion has been compensated by applying the PCA algorithm to the Figure 10 The tissue contributing factors are calculated as PCA components obtained from the ORG response.

[0039] Figure 12 The absolute values ​​of the principal component coefficients returned by the PCA algorithm in the variant are shown.

[0040] Figure 13 is a flow chart showing a method by which the data processing apparatus 100 of the exemplary embodiment may determine the determined first indication Z B-1 and / or second indication Z B-2 is reliable or unreliable (as the case may be), and generates an indicator indicating the determined reliability / unreliability.

[0041] Figure 14A The first pair of B-scans that are highly correlated are shown.

[0042] Figure 14B Shows the Figure 14A Calculate the cross-correlation of the B-scans.

[0043] Figure 15A A second pair of B-scans that are not highly correlated is shown.

[0044] Figure 15B Shows the Figure 15A Calculate the cross-correlation of the B-scans.

[0045] Figure 16 is a flow chart illustrating a process by which the data processing apparatus 100 of the second example embodiment herein may process pairs of adjacent B-scans in a series of B-scans to generate ORG data indicative of the response of the retina to applied light stimulation.

[0046] Figure 17 It shows that the data processing device can be used according to Figure 16 An example of a concatenation of velocity profiles generated by process S180.

[0047] Figure 18 A graph showing the change in optical path length over time ΔOPL is shown, which has been calculated by the data processing apparatus of the second example embodiment for selected positions in the velocity profiles in the cascade of velocity profiles.

[0048] Figure 19 A graph of an example ORG signal generated by the data processing apparatus 100 of the example embodiment herein is shown.

[0049] Figure 20 Shows smooth Figure 19The result of the ORG signal in .

[0050] Figure 21 is a flowchart illustrating a process by which a data processing apparatus according to a second example embodiment of the present invention may process adjacent OCT images in a series of OCT images to generate ORG data indicating a response of the retina to applied light stimulation.

[0051] Figure 22 is a flowchart illustrating a process by which the data processing apparatus of the third example embodiment may process adjacent OCT images in a series of OCT images to generate ORG data indicating a response of the retina to applied light stimulation. Detailed Description of Example Embodiments

[0052] To analyze the ORG response, different layers of the retina are typically selected and the changes in optical path length (OPL) between the layers of interest are calculated. In the case of processing OCT images in the form of B-scans, a flattening algorithm is typically used to distort the B-scan so that the retinal layers are located at a constant depth in the distorted B-scan. The layers are then typically found by manual selection or by automatically finding the depth positions of the intensity peaks in the B-scan corresponding to the retinal layers. The manual approach is undesirable because it requires input from the operator and can introduce bias and slow down the analysis, and also requires operator knowledge. Although automatic discovery of intensity peaks generally has a high success rate, this can be problematic when the OCT signal is weak and the different layers are not well defined, particularly in scans with low axial resolution, in scans of the retinal periphery, or when the layer of interest has weak intensity (e.g., the inner ganglion layer).

[0053] The OCT data processing methods described herein may address some or all of these issues and do not rely on identifying OCT signal intensity peaks, which may be difficult when the OCT data has poor axial resolution or is noisy. At least some of the OCT data processing methods described herein may allow layer segmentation to be efficiently performed on OCT data from an FD-OCT system with poor axial resolution, which may be insufficient to show, for example, COST, rod outer segment (ROST) and / or retinal pigment epithelium (RPE) intensity peaks, and allow such boundaries to be identified. In addition, they may be advantageous in allowing ORG responses to be extracted from different regions of the retina (such as ganglion cells) that typically provide weak OCT signal intensities. At least some of the OCT data processing methods described herein also allow the use of a conventional point scanning FD-OCT system with sequential registration and without adaptive optics or tracking systems to discover conventional optical path length responses between photoreceptor layers.

[0054] First exemplary embodiment

[0055] Figure 1 is a schematic diagram of a data processing apparatus 100 according to a first exemplary embodiment. The data processing apparatus 100 is arranged to process complex OCT data generated by a Fourier domain OCT (FD-OCT) imaging system. More specifically, the data processing apparatus 100 is arranged to process a phase component (phase information) of a first OCT image and a phase component of a second OCT image in a series of OCT images of a common portion of the retina acquired by an FD-OCT imaging system 30 having phase resolution capability after the common portion of the retina of the eye 20 is stimulated by an optical stimulus 40 generated by an optical stimulus source 45, such as a light emitting diode (LED). The series of OCT images 10 are acquired by the FD-OCT imaging system 30 when the thickness of a layer L in the common portion of the retina changes in response to the optical stimulus 40, and may form part of a wider series of OCT images acquired by the FD-OCT imaging system 30, the wider series of OCT images also including OCT images acquired before the application of the optical stimulus 40.

[0056] As in the present example embodiment, the FD-OCT imaging system 30 may be a swept source OCT (SS-OCT) system. However, the FD-OCT imaging system 30 need not be provided in this form and may, for example, take the alternative form of spectral domain OCT (SD-OCT). More specifically, the example embodiments may be provided as any form of FD-OCT imaging system with phase resolution capability that is capable of generating complex OCT data, i.e., a Fourier transform of individual spectral interference patterns (interference spectra) representing complex A-scan information obtained for each scan position at which an OCT measurement is taken during a scan. Such complex OCT data encodes phase information from the acquired OCT measurements that may be used by the data processing apparatus 100 as described herein to identify one or two boundaries of a layer of the retina of the eye 20.

[0057] The FD-OCT imaging system 30 may include well-known components, including a scanning system, a photodetector, OCT data processing hardware, and a beam generator (not shown). The scanning system may be arranged to perform one-dimensional and / or two-dimensional point scanning of a light beam across the retina and collect light that has been scattered by the retina during the point scanning. The scanning system is therefore arranged to acquire A-scans at various scanning locations distributed on the surface of the retina by sequentially illuminating the scanning locations with a light beam (one scanning location at a time) and collecting at least a portion of the light scattered by the retina at each scanning location. The scanning system may acquire OCT images in the form of repeated B-scans by performing point scanning using a linear scanning pattern, wherein a set of overlapping scan lines is followed during scanning. Although in this exemplary embodiment, the scanning system is arranged to acquire repeated B-scans by performing point scanning, in other exemplary embodiments, the scanning system may alternatively be arranged to acquire repeated B-scans by performing line scanning using hardware well known to those skilled in the art. The FD-OCT imaging system 30 may alternatively be arranged to acquire OCT images in the form of C-scans by performing point scanning or line scanning using techniques well known to those skilled in the art, or by employing a full-field setup.

[0058] like Figure 1 As shown, in this exemplary embodiment, the first OCT image described above may be in the form of a first OCT B-scan 10-1, and the second OCT image may be in the form of a second OCT B-scan 10-2. However, the form of the first and second OCT images is not limited thereto, and, for example, the first and second OCT images may be corresponding OCT C-scans. As in this exemplary embodiment, the first B-scan 10-1 and the second B-scan 10-2 of a common portion of the retina may be adjacent to each other in a series of OCT B-scans 10, or they may be separated by n intermediate B-scans, where n ≥ 1 but is sufficiently small so that the B-scans have a sufficiently high degree of phase correlation with each other due to the absence of significant relative movement of the eye 20 relative to the FD-OCT imaging system 30. Furthermore, as explained in more detail below, the first B-scan 10-1 may be processed as one of a first set of adjacent B-scans in the series of B-scans 10, and the second B-scan 10-2 may be processed as one of a second set of adjacent B-scans in the series of B-scans 10, wherein the first and second sets of adjacent OCT images may or may not have one or more OCT images in common.

[0059] As in this example embodiment, layer L may be the photoreceptor outer segments (OS) of the retina, which contain rods and cones and serve to convert absorbed visible light signals into changes in membrane potential. The photoreceptor OS is particularly well suited for ORG measurements; the photoreceptors therein are long and narrow, behave like optical waveguides, and provide relatively strong reflections from each end (i.e., the IS / OS junction and the COST). Phase changes between these reflections are used to quantify changes in the optical path length of the OS. However, layer L need not include (or be limited to) the photoreceptor OS, and may instead be another layer or layers of the retina that produce measurable thickness changes when stimulated, such as the ganglion cell layer / inner plexiform layer (GCL / IPL). The GCL and IPL will provide much weaker reflections that are still detectable, especially if steps are taken to suppress motion artifacts (see, e.g., C. et al., “Simultaneous functional imaging of neuronal and photoreceptor layers in living human retina,” Optics Letters, vol. 44, no. 23, pp. 5671–5674 (1 December 2019).

[0060] As described in more detail below, processing of the phase components of the first B-scan 10-1 and the second B-scan 10-2 by the apparatus 100 results in determination by the apparatus 100 of a first indication Z of a first position of a first boundary (edge) B-1 of a layer L of the retina along an axial direction (z-axis direction) z in the B-scan 10. B-1 and a second indication Z of a second position of a second boundary B-2 of the layer L along the axial direction z in the B-scan 10 B-2 One or both of the following. The axial direction (or z-axis direction) in a B-scan is the direction along which the elements of the A-scan, the components of the B-scan, are arranged in any B-scan. In the case where layer L is the photoreceptor OS, as in the present example embodiment, the first boundary B-1 of layer L is the IS / OS junction, and the second boundary B-2 of layer L is the COST. As in the present example embodiment, the first indicator Z B-1 and the second indication Z B-2 Each of the Z may provide a corresponding row index for the row of the B-scan 10 at which the corresponding boundary is located, although the Z B-1 and Z B-2 Scaled versions of these indices may be used instead (eg, where the depth of the IS / OS junction and the COST from the retinal surface are required).

[0061] As described in more detail below, the data processing apparatus 100 is arranged to process at least a portion of the phase component of the first B-scan 10-1 and at least a portion of the phase component of the second B-scan 10-2 to calculate a tissue velocity profile P comprising tissue velocity values ​​which, as a function of position in the velocity profile, are indicative of a distribution of velocities among points along the axial direction z in a common portion of the retina (e.g., as Figure 1 ). Thus, as in this example embodiment, the velocity profile P may be a two-dimensional data structure in which values ​​of tissue velocity v are associated with corresponding values ​​of position. The velocity profile P may more generally be a two-dimensional data structure in which values ​​indicative of velocity v (e.g., calculated changes in phase or optical path length) are associated with corresponding values ​​of position. The data processing apparatus 100 is arranged to determine the first indication Z based on the calculated velocity profile P. B-1 and the second indication Z B-2 One or two of .

[0062] The data processing device 100 may be provided in any suitable form, for example as Figure 2 The programmable signal processing hardware 200 is provided as schematically shown in FIG. The programmable signal processing hardware 200 includes a communication interface (I / F) 210 for receiving the B-scan 10 (or another form of OCT image, such as a C-scan) from the FD-OCT imaging system 30 and outputting a first indication Z B-1 and the second indication Z B-2The communication interface 210 may also output a reliability indicator, which will be described in more detail below. The signal processing hardware 200 also includes a processor 220 (e.g., a central processing unit CPU and / or a graphics processing unit GPU), a working memory 230 (e.g., a random access memory), and an instruction storage device 240 storing a computer program 245, which includes computer readable instructions that, when executed by the processor 220, cause the processor 220 to perform various functions of the data processing device 100 described herein. The working memory 230 stores information used by the processor 220 during execution of the computer program 245. The instruction storage device 240 may include a ROM (e.g., in the form of an electrically erasable programmable read-only memory (EEPROM) or flash memory) preloaded with computer readable instructions. Alternatively, the instruction storage device 240 may include RAM or a similar type of memory, and the computer-readable instructions of the computer program 245 may be input to the instruction storage device 240 from a computer program product (e.g., a non-transitory computer-readable storage medium 250 in the form of a CD-ROM, DVDROM, etc.) or a computer-readable signal 260 carrying the computer-readable instructions. In any case, when the computer program 245 is executed by the processor 220, the processor 220 performs the functions of the data processing device 100 described herein. Thus, the data processing device 100 of the example embodiment may include a computer processor 220 and a memory 240 storing computer-readable instructions that, when executed by the processor 220, cause the processor 220 to process the phase components of the first B-scan 10-1 and the phase components of the second B-scan 10-2 to calculate a velocity profile P, and determine a first indication Z based on the calculated velocity profile P. B-1 and / or second indication Z B-2 .

[0063] However, it should be noted that the data processing apparatus 100 may alternatively be implemented in non-programmable hardware (such as an ASIC, FPGA or other integrated circuit dedicated to performing the functions of the data processing apparatus 100 described herein), or in a combination of such non-programmable hardware and the above-referenced Figure 2 The system is implemented in a combination of programmable signal processing hardware 200 as described.

[0064] The data processing apparatus 100 may be provided as a standalone product or as part of a system 1000 comprising an optical stimulus source 45 arranged to provide optical stimulus 40 and an FD-OCT imaging system 30 arranged to acquire a series of OCT images 10 of a common portion of the retina of the eye 20 after stimulating the common portion by the optical stimulus 40. The data processing apparatus 100 is arranged to process a phase component of a first OCT image 10-1 and a phase component of a second OCT image 10-2 of the series of OCT images 10 acquired by the FD-OCT imaging system 30 using the techniques described herein.

[0065] Figure 3 is a flowchart illustrating a method according to an example embodiment, by which a data processing apparatus 100 processes phase components of B-scans 10-1 and 10-2 from a series of B-scans 10 of a common portion of the retina acquired after optical stimulation of the common portion of the retina to identify one or both of a first boundary B-1 and a second boundary B-2 of a layer L, such as in the present example embodiment, the first boundary B-1 and the second boundary B-2 may be an IS / OS junction and a COST, respectively.

[0066] exist Figure 3 In the process S10, the data processing device 100 processes the phase component of the first B-scan 10-1 and the phase component of the second B-scan 10-2 to calculate a velocity profile P indicating the distribution of the velocity v in the common portion of the retina along the axial direction z. Thus, the data processing device 100 calculates the velocity profile P including the values ​​of the velocity v (or the values ​​of the change in phase or optical path length, such as ΔOPL), the changes in the values ​​of the velocity v with position in the velocity profile indicating the distribution of the velocity v of each portion of the retina (in the common portion) located at the corresponding point on the line aligned with the axial direction. The velocity profile P is thus a two-dimensional data structure in which the values ​​of the velocity v (or the values ​​of the variables indicating the velocity v, such as ΔOPL or the change in phase) are calculated. ) is associated with a corresponding value of the position, and the data structure can be used to store the velocity v (or or ΔOPL) can be visualized in the form of a graph, such as Figure 1 Schematically shown in FIG.

[0067] Before calculating the velocity profile P, as in the present example embodiment, the data processing apparatus 100 may process the phase component of the first OCT image 10-1 and the phase component of the second OCT image 10-2 to compensate for global motion of the common portion of the retina during acquisition of a series of OCT images 10 by the Fourier domain OCT imaging system 30 after the common portion of the retina is stimulated by the optical stimulus 40. Compensating for global motion has been found to increase the reliability of the velocity-based layer segmentation technique described herein.

[0068] Figure 4A The overall oscillation of a sinusoidal waveform is shown added to the oscillatory movement of two model retinal layers. In this example, by adding the phase to model the overall motion, where Represents a fraction of a wavelength (in a two-way reflection, a half-wavelength shift is 2π). Figure 4A shows how the phase angle of light reflected from the first model retinal layer (layer 1) and the second model retinal layer (layer 2) varies with time while these layers undergo a common overall motion as well as their individual oscillations (i.e., Figure 4A Graphs labeled “Layer 1 + Whole” and “Layer 2 + Whole”). Figure 4A It also shows how the phase angle of light reflected from any model retinal layer changes with time while the layer is subjected only to bulk motion (i.e., Figure 4A ). Figure 4B Shown is the effect of removing global motion on the temporal variation of the phase angle of light reflected from two model retinal layers.

[0069] Used to generate Figure 4A and Figure 4B MATLAB plots in TM The code is as follows:

[0070] t=(0:1000);

[0071] A = 0.95;

[0072] B = 0.8;

[0073] RetinaBulkMovement=exp(j*cos(t*pi / sqrt(700)));

[0074] Layer1Movement=exp(j*cos(t*pi / sqrt(30)));

[0075] Layer2Movement=exp(j*cos(t*pi / sqrt(50)));

[0076] Layer1totalMovement = Layer1Movement * RetinaBulkMovement;

[0077] Layer2totalMovement = Layer2Movement * RetinaBulkMovement;

[0078] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[0079] Layer1MovementDeduced = angle(Layer1totalMovement * conj(RetinaBulkMovement));

[0080] Layer2MovementDeduced = angle(Layer2totalMovement * conj(RetinaBulkMovement));

[0081] subplot(211);

[0082] plot(angle(Layer1totalMovement)); hold

[0083] on; plot(angle(Layer2totalMovement)); plot(angle(RetinaBulkMovement)); legend('Layer1

[0084] +Bulk','Layer2+Bulk','Bulk');

[0085] xlim([300 400]);

[0086] subplot(212);

[0087] plot(Layer1MovementDeduced); hold on; plot(Layer2MovementDeduced); shg; hold off;

[0088] xlim([300 400]);

[0089] legend('Layer1','Layer2');

[0090] Refer again Figure 3In process S20, the data processing device 100 determines a first indication Z based on the calculated velocity distribution map. B-1 and the second indication Z B-2 Thus, in contrast to conventional OCT signal amplitude-based slice segmentation typically employed in ORG (whether performed manually by a user or automatically), velocity-based slice segmentation is used in this exemplary embodiment.

[0091] The data processing apparatus 100 may now be described as follows: Figure 3 An example of the process of S10 of , which is based on the velocity-based ORG technology described in Kari V. Vienola et al., “Velocity-based optoretinography for clinical applications,” Optica 9, pp. 1100-1108 (2022), the contents of which are incorporated herein by reference in their entirety.

[0092] exist Figure 3 In process S10, as in this exemplary embodiment, the data processing device 100 may first flatten the first B-scan 10-1 and the second B-scan 10-2 so that the IS / OS and COST reflections are located at substantially the same height for each A-scan in each B-scan. The data processing device 100 may then register the second B-scan 10-2 relative to the first B-scan 10-1. The phase data of the two B-scans for each spatial coordinate pair (i.e., phase data at the same (x, z) coordinate in the B-scan) may then be unfolded in the time dimension to minimize the magnitude of the phase difference between the datasets of the two B-scans 10-1 and 10-2. After unfolding and processing the phase components of the first and second OCT images 10-1 and 10-2 to compensate for global motion of the retina during imaging, the difference between the corresponding phase values ​​in the B-scans 10-1 and 10-2 is calculated for each spatial location specified by the corresponding coordinate pair, and then the instantaneous velocity for the spatial location is calculated using the difference between the acquisition times of the B-scans. As in this example embodiment, these instantaneous velocities can be averaged over the lateral dimension (ie, along the x-axis) to give an instantaneous, depth-dependent measure of velocity along the z-axis, ie Figure 1A one-dimensional velocity profile P of the type shown schematically in FIG1 indicates the distribution of velocity v in the axial direction (z-axis) within the common portion of the retina covered by the first B-scan 10-1 and the second B-scan 10-2. The B-scan amplitudes can also be averaged in the transverse dimension (x-axis) if desired to provide an instantaneous, depth-dependent measure of backscatter. It should be noted that averaging in the transverse dimension is not necessary, and a two-dimensional velocity profile can instead be generated that indicates the distribution of velocity v in both the axial direction (z-axis) and the transverse direction (x-axis) within the common portion of the retina covered by the first B-scan 10-1 and the second B-scan 10-2.

[0093] Once the boundaries of the OS are identified as described below, their corresponding velocities can be extracted, and the difference between them provides the contraction / elongation rate of the OS when the B-scans 10-1 and 10-2 are acquired by the FD-OCT imaging system 30, which can be used to generate ORG data indicating the response of the retina in the common portion to the applied stimulus.

[0094] exist Figure 3 In the process S20, the data processing device 100 can generate a maximum speed v in the speed distribution map P by max The first position z max Determined as the first indicator Z B-1 , determining a first indication Z of a first position of a first boundary B-1 of the layer L in the axial direction z in the B-scan 10 B-1 The maximum value here may be a global maximum value of the velocity in the velocity profile P or a maximum value within a predefined portion in which the OS of the velocity profile P is expected to be located. For example, the predetermined portion may be defined relative to the top surface of the retina based on the known physiology of the eye. Similarly, the data processing apparatus 100 may define the maximum value v in the velocity profile P by placing the maximum value v in the velocity profile P corresponding to the velocity in the velocity profile P. min The second position z min Determined as the second indication Z B-2 ,Sure Figure 3 The second indication Z of the second position of the second boundary B-2 of the layer L in the process S20 along the axial direction z in the B scan 10 B-2 As mentioned above, the minimum here may be a global minimum of the speed in the speed profile P or a minimum within a predefined portion in which the OS is expected to lie.

[0095] Thus, the velocity v in a portion of the velocity profile P corresponding to a photoreceptor OS in the retina is (typically monotonically) modulated from a first position z in the velocity profile P corresponding to the IS-OS junction at a given location on the retina. max The maximum speed v maxChange to the second position z of the COST corresponding to the location on the retina in the velocity distribution map P min The minimum velocity v at min , and the first indication Z B-1 Indicates maximum speed v max The first position z in the velocity profile P max , the second indicator Z B-2 Indicates minimum speed v min The second position z in the velocity profile P min .

[0096] The velocity-based segmentation methods described herein can be used to determine the location of one or both of the boundaries of the OS or other layers L of the retina that change in response to an applied optical stimulus (e.g., the location of the rod outer segment tips (ROSTs) or the retinal pigment epithelium (RPE)). min , the velocity profile P can have a Figure 1 In this example embodiment, the data processing apparatus 100 first determines the candidate first position z in the velocity profile P by first determining the maximum value of the z axis in the velocity profile P and then determining the minimum value of the z axis in the velocity profile P. i and the candidate second position z in the velocity profile P j All combinations of (z i ,z j ) for each combination of the candidate first position and the candidate second position in the calculation of the combination (z i ,z j ) in the candidate first position z i The corresponding velocity v i and candidate second position z j The corresponding velocity v j The difference between ij The corresponding value of the velocity distribution diagram P determines the first position z max and the second position z min Then, the data processing device 100 identifies the candidate first position z of the following combinations among the plurality of combinations: i and candidate second position z j As the first position z max and the second position z min : For this combination in multiple combinations, the calculated difference D ij The value of is the largest of the calculated differences.

[0097] In a variation of this exemplary embodiment, the data processing apparatus 100 first processes the velocity profile P for a candidate first position z i and the candidate second position z in the velocity profile P jAll combinations of (z i ,z j ) in each combination of the candidate first position and the candidate second position, calculate the combination of the candidate first position z i is the center (but may also include candidate first positions z i )'s five adjacent positions z i-2 、z i-1 、z i 、z i+1 and z i+2 The corresponding velocity v at the collection i-2 、v i-1 、v i 、v i+1 and v i+2 The average value of the combination and the candidate second position z j is the center (but may also include candidate second positions z j )'s five adjacent positions z j-2 、z j-1 、z j 、z j+1 and z j+2 The corresponding velocity v at the collection j-2 、v j-1 、v j 、v j+1 and v j+2 The corresponding difference between the average values ​​of max and the second position z min Although the number of adjacent positions is five in this example embodiment, it should be understood that this is given by way of example only, and in other embodiments there may be a different number of adjacent positions (typically, two or more). The data processing apparatus 100 then identifies candidate first positions z for the following combinations among the plurality of combinations: i and candidate second position z j As the first position z max and the second position z min : For this combination in multiple combinations, the calculated difference D ij The value of is the largest among the calculated differences. Find the difference D ij The maximum value of gives the maximum change in optical path length and thus indicates the position where the neuron undergoes the largest change due to the stimulation.

[0098] In another variation of this exemplary embodiment, the data processing device 100 is arranged to determine the first indication Z by processing the calculated velocity profile P using a cluster analysis algorithm or a dimensionality reduction algorithm. B-1 and / or second indication Z B-2. These methods of processing the velocity profile P can be used to identify changes associated with retinal layers that may be hidden by other layers and / or anatomical features. The dimensionality reduction algorithm can be one of several different types known to those skilled in the art and can be based on machine learning (ML) methods. For example, the dimensionality reduction algorithm can include a principal component analysis (PCA) algorithm, an independent component analysis (ICA) algorithm, a linear discriminant analysis (LDA) algorithm, or a non-negative matrix factorization (NMF) algorithm.

[0099] In the case where the dimensionality reduction algorithm is a PCA algorithm, as in this variant, the first indicator Z B-1 Or the second indication Z B-2 One or both of can be determined using the principal components that have been determined by applying the PCA algorithm to the calculated velocity profile P. To better understand how this is done, reference will now be made to Figures 5 to 12 Describe how to use MATLAB TM PCA was used in the present study to identify data elements of two model retinal layers that contribute to the ORG response (these data elements are referred to as "ORG contributors" below), and to locate these ORG contributors, as well as to identify to which layer they belong, thereby allowing for the classification of some illustrative examples of retinal layers that exhibit ORG responses.

[0100] In these examples, a single A-scan is considered for simplicity, and ORG contributors are provided at multiple locations for each of the two layers in the A-scan. In the first of these examples, the ORG contributors for the two retinal layers are set separately from each other and do not overlap, as shown in FIG. Figure 5 where the x-axis represents the location of the contributing factor within the A-scan and the y-axis represents the intensity of the contributing factor. Figure 5 MATLAB TM The code is as follows:

[0101] l1=transpose([0 0.15 0.4 0.25 0 0 0 0 0 0]);

[0102] l2=transpose([0 0 0 0 0 0 0.15 0.5 0.3 0]);

[0103] figure;

[0104] plot(l1);hold on;plot(l2);hold off;

[0105] legend('Retina element 1','Retina element 2')

[0106] The phases of model retinal layer 1 and model retinal layer 2 are then assigned to Figure 5 The ORG contributing factors are shown. For the case where no compensation for global motion is provided, the magnitude of the (complex) model retinal ORG response is Figure 6 The phase of the ORG response is shown in Figure 7A As shown in. For generating Figure 6 and Figure 7A MATLAB TM The code is as follows:

[0107] RetinaORG=l1*Layer1totalMovement+l2*Layer2totalMovement;

[0108] figure;

[0109] imagesc(abs(RetinaORG)); colormap gray; title('OCT amplitude')

[0110] imagesc(angle(RetinaORG));colormap parula;colorbar;

[0111] Figure 7B The phase of the ORG response is shown for the case where the overall motion has been compensated. Figure 7B MATLAB TM The code is as follows:

[0112] RetinaORGComp=RetinaORG.*conj(RetinaBulkMovement);

[0113] imagesc(angle(RetinaORGComp));colormap parula;colorbar;

[0114] PCA was then used to automatically detect the main ORG contributors and identify their localization in the A-scans. TM The PCA algorithm is applied to the case where the overall motion has been compensated (the phase of the ORG response is in Figure 7B The (complex) ORG response calculated as shown in Figure 2 yields two PCA components, such as Figure 8 As shown. Used to generate Figure 8 MATLAB TM The code is as follows:

[0115] [COEFF2,SCORE2,LATENT2]=pca((RetinaORGComp_t));

[0116] stem(LATENT2);

[0117] The absolute values ​​of the principal component coefficients (i.e., loadings) returned by the PCA algorithm are plotted on Figure 9 Used to generate Figure 9 MATLAB TM The code is as follows:

[0118] plot(abs(COEFF2(:,1:2)));

[0119] legend('Retina element PCA1','Retina element PCA2');

[0120] exist Figure 9 In the , the ORG contributors associated with layer 1 (labeled “retinal element PCA1”) appear to the right of the ORG contributors associated with layer 2 (labeled “retinal element PCA2”), which is consistent with Figure 5 In contrast, Figure 5 The ORG contributor associated with layer 1 (labeled "Retina Element 1") appears to the left of the ORG contributor associated with layer 2 (labeled "Retina Element 2"). Figure 9 ORG enabler values ​​in Figure 5 The values ​​in are also slightly different. However, the PCA algorithm has identified where the main contributing factors come from and grouped them.

[0121] In the example above, there is no spatial overlap between the ORG enablers associated with Tier 1 and Tier 2, as shown in Figure 5 However, in a more realistic scenario, positioning in the A-scan can collect both ORG contributors and will now refer to Figures 10 to 12 Describe the application of PCA to a variation of the above example, in which the ORG contributors of two retinal layers partially overlap with each other.

[0122] Figure 10 The ORG contributors for two retinal layers in the above variation are shown, where the x-axis again represents the location of the contributor within the A-scan and the y-axis represents the intensity of the contributor. Figure 10 MATLAB TM The code is as follows:

[0123] l1b=transpose([0 0.15 0.4 0.25 0 0 0 0 0 0]);

[0124] l2b=transpose([0 0 0 0.15 0.5 0.3 0 0 0 0]);

[0125] figure();

[0126] plot(l1b); hold on; plot(l2b); hold off;

[0127] legend('Retina element 1','Retina element 2')

[0128] PCA was then used to automatically detect the main ORG contributors and identify their localization in the A-scans. TM Applying the PCA algorithm to the ORG response calculated for the case where the overall motion has been compensated again produces two PCA components, such as Figure 11 As shown. Used to generate Figure 11 MATLAB TM The code is as follows:

[0129] RetinaORGb=l1b*Layer1totalMovement+l2b*Layer2totalMovement;

[0130] RetinaORGCompb=RetinaORGb.*conj(RetinaBulkMovement);

[0131] RetinaORGCompb_t=transpose(RetinaORGCompb);

[0132] [COEFFb,SCOREb,LATENTb]=pca((RetinaORGCompb_t));

[0133] stem(LATENTb);

[0134] The absolute values ​​of the principal component coefficients returned by the PCA algorithm are plotted on Figure 12 Used to generate Figure 12 MATLAB TM The code is as follows:

[0135] plot(abs(COEFFb(:,1:2)));

[0136] legend('Retina element PCA 1','Retina element PCA2');

[0137] As in Figure 12 Although the layers are shown in reverse order and have slightly different values ​​(as in Figure 9 ), but the PCA algorithm has again identified where the main contributing factors come from and has grouped them.

[0138] In addition to or as an alternative to using the velocity-based segmentation method described herein to determine the positions of two boundaries of layer L, the data processing apparatus 100 may provide functionality for determining the position of one of the boundaries of layer L using a conventional OCT signal amplitude-based method and determining the position of the remaining boundary using a velocity-based method, as will now be described.

[0139] In a variation of the exemplary embodiment having this functionality, the data processing device 100 determines a first position z along the axial direction z in the B-scan 10. max The first indication Z B-1 , the first position z max is the position of the first boundary B-1 of the layer L. This determination is based on an amplitude component of at least one B-scan in the series of B-scans 10, preferably based on one of the B-scans 10-1 and 10-2 from which the velocity profile P was generated as described above. In the case where the layer L is the photoreceptor OS, as in the present example, the first boundary B-1 of the layer L is the IS / OS junction and the second boundary B-2 of the layer L is the COST, each of which is associated with a maximum value or corresponding reflectivity peak in the amplitude of the OCT signal in one of the B-scans of the series of B-scans 10 recorded by the FD-OCT imaging system 30. The data processing device 100 can use any image-appropriate image processing technique known to those skilled in the art to determine this amplitude peak in a predefined region of the B-scan expected to contain the IS / OS junction and the COST, for example by finding a peak in a moving average of the amplitude values ​​of the OCT signal of one or more A-scans of the B-scan taken along the axial direction in the B-scan. The data processing device 100 is further arranged to determine the peak value in the velocity profile P by identifying a minimum velocity v in a portion of the velocity profile P that is in agreement with the minimum velocity v in the velocity profile P. min The corresponding second position z min To determine the second indication Z (the position of the IS / OS junction of the OS in the axial direction z in the B-scan) B-2 The data processing device 100 can first calculate the position z corresponding to the determined amplitude peak in the velocity distribution map P. peak and the candidate position z in the velocity profile P i All combinations of (z peak ,z i ) in the position z peak and candidate position zi For each combination of peak ,z i ) in z peak The corresponding velocity v peak and at candidate position z i The corresponding velocity v i The difference between peak_i The corresponding value of , identifies the second position z in the velocity profile P min Then, the data processing device 100 selects the candidate position z of the following combinations from among the multiple combinations: i Identified as the second position z min :For the combination in multiple combinations, the difference D peak_i The calculated value of is the largest among the calculated values ​​of difference.

[0140] In another variation of the exemplary embodiment having the above functionality, the data processing device 100 determines a second position z along the axial direction z in the B-scan 10. min The second instruction Z B-2 , the second position z min is the location of the second boundary B-2 of layer L. This determination is based on an amplitude component of at least one B-scan in the series of B-scans 10, preferably one of the B-scans 10-1 and 10-2 from which the velocity profile P was generated as described above. In the case where layer L is the photoreceptor OS, as in the present example, the second boundary B-2 of layer L is the COST, and the first boundary B-1 of layer L is the IS / OS junction, each of which is associated with a maximum value, or corresponding reflectivity peak, in the amplitude of the OCT signal in one of the B-scans of the series of B-scans 10 recorded by the FD-OCT imaging system 30. The data processing device 100 can use any image-appropriate image processing technique known to those skilled in the art to determine this amplitude peak in a predefined region of the B-scan expected to contain the IS / OS junction and the COST, for example by finding a peak in a moving average of the amplitude values ​​of the OCT signal of one or more A-scans of the B-scan taken along the axial direction in the B-scan. The data processing device 100 is then further arranged to determine the maximum velocity v in the velocity profile P by identifying a portion of the velocity profile P that corresponds to the maximum velocity v in the velocity profile P. max The corresponding first position z max To determine the first indication Zv (the position of the COST of the OS in the axial direction z in the B-scan) B-1 The data processing device 100 can first calculate the position z corresponding to the determined amplitude peak in the velocity distribution map P. peak and the candidate position z in the velocity profile P i All combinations of (z peak ,z i ) in the position zpeak and candidate position z i For each combination of peak ,z i ) in z peak The corresponding velocity v peak and at candidate position z i The corresponding velocity v i The difference between peak_i The corresponding value of , identifies the first position z in the velocity profile P max Then, the data processing device 100 selects the candidate position z of the following combinations from among the multiple combinations: i Identified as the first position z max :For the combination in multiple combinations, the difference D peak_i The calculated value of is the largest among the calculated values ​​of difference.

[0141] The data processing apparatus 100 of this exemplary embodiment or any of its variants described above may also be arranged to determine the first indication Z determined by the data processing apparatus 100. B-1 and / or second indication Z B-2 The method determines whether the velocity-based OS segmentation performed by the data processing apparatus 100 is reliable or unreliable (as the case may be), and generates an indicator indicating the determined reliability / unreliability. The reliability of the velocity-based OS segmentation performed by the data processing apparatus 100 will be affected by the image quality of the first B-scan 10-1 and the second B-scan 10-2 (although to a lesser extent than many amplitude-based segmentation methods) and how closely the B-scans are correlated with each other. The maximum value of the calculated cross-correlation between the two B-scans provides a good measure of how closely the B-scans are correlated, and therefore provides a first indicator Z determined by the data processing apparatus 100. B-1 and / or second indication Z B-2 The data processing device 100 may communicate the determined reliability indicator to the user (e.g., by displaying a message or other kind of graphic on the screen that informs the user of the determined reliability / unreliability). Additionally or alternatively, the data processing device 100 may store the first indicator Z that the data processing device 100 has determined. B-1 and / or second indication Z B-2 (As appropriate) The associated indicator. The stored indicator can be used in subsequent data processing operations, as described below.

[0142] More specifically, the data processing apparatus 100 of the present exemplary embodiment or any of its variations described above may also be arranged to perform Figure 13 The process is schematically shown in FIG.

[0143] exist Figure 13In process S30, the data processing device 100 calculates a comparison value based on the first B-scan 10-1 and the second B-scan 10-2. The comparison value may be a value of an image quality metric or attribute calculated based on at least one of the amplitude components of the first B-scan 10-1 and the amplitude components of the second B-scan 10-2. The image quality metric (attribute) may take the form of a dynamic range, contrast, signal-to-noise ratio, or sharpness function value derived from the B-scan. A measure of the success of intensity-based segmentation of the B-scan is another example of an image quality metric that may be used. Alternatively, the comparison value may be the value of an image similarity metric that provides a measure of the degree of similarity between the images, the image similarity metric being calculated based on the first B-scan 10-1 and the second B-scan 10-2. As in this example embodiment, the image similarity metric may take the form of a maximum value of the cross-correlation between the first B-scan 10-1 and the second B-scan 10-2. Other alternative image similarity metrics include sum of squared differences, mutual information, normalized mutual information, and Kullback-Leibler distance.

[0144] The cross-correlation between two complex functions f(t) and g(t) of a real variable t is denoted as f*g and is defined as where * denotes convolution, and is the complex conjugate of f(t).

[0145] The maximum value of the cross-correlation between two B-scans is a good measure of how closely the B-scans are related to each other. For example, Figure 14A The two B-scans shown in FIG are highly correlated with each other, and the cross-correlation calculated between these B-scans has a value of approximately 3·10 4 The peak value, such as Figure 14B shown.

[0146] For comparison, Figure 15A Two dissimilar B-scans are shown. This dissimilarity could be caused by factors such as eye movement, for example, and would result in poor image registration. Figure 14A For the B-scan shown in Figure 15A For the example B-scans shown in FIG, the maximum value of the calculated cross-correlation between the B-scans is much lower and is approximately 1·10 4 ,like Figure 15B shown.

[0147] Reference again Figure 13 In process S40, the data processing device 100 compares the comparison value with a threshold value to determine whether the comparison value is equal to or greater than the threshold value. The threshold value can be set by comparing the segmentation result from the data processing device 100 (i.e., the segmentation result of Z B-1 and / or Z B-2The boundary or boundaries of the layer L indicated by the value of (as the case may be) and the layer segmentation resulting from the user's inspection of the source B-scans, while taking into account the maximum cross-correlation value calculated for the source B-scans. In this way, it can be determined that pairs of B-scans between which the maximum value of the cross-correlation calculated is below a certain threshold tend to produce Z when the B-scans are processed by the data processing apparatus 100 described herein. B-1 and / or Z B-2 (as the case may be), and pairs of B-scans for which the maximum value of the cross-correlation calculated between them is equal to or higher than the threshold value tend to produce Z B-1 and / or Z B-2 (as the case may be) a reliable value.

[0148] For example, it can be found Figure 14A The B-scan shown in the Z B-1 and / or Z B-2 (as appropriate), which compares favourably with the results of manual segmentation of OS performed by examining these B-scans, while Figure 15A The B-scan shown in the Z B-1 and / or Z B-2 (as the case may be), which differs significantly from the results of manual segmentation of OS performed by examining these B-scans. In this case, it may be appropriate to set the threshold to Figure 14B and Figure 15B The maximum values ​​of the cross-correlations shown (i.e., approximately 3·10 4 and 1.10 4 ), for example, about 1.4·10 4 The value of Figure 15B This is shown by the horizontal black line in . B-scans that produce maximum cross-correlation values ​​below this threshold can be considered not allowing reliable segmentation to be performed.

[0149] In the case where it is determined in the process S40 that the comparison value is equal to or greater than the threshold value ("Yes" at S50), the data processing apparatus 100 Figure 13 In the process S60, a first indicator is generated, which indicates that Figure 3 The first indication Z determined in process S20 B-1 and / or second indication Z B-2 It is reliable.

[0150] In the case where it is determined in the process S40 that the comparison value is not equal to or greater than the threshold value (“No” at S50 ), the data processing apparatus 100 Figure 13 In the process S70, a second indicator is generated, which indicates that Figure 3 The first indication Z determined in process S20 B-1and / or second indication Z B-2 is unreliable.

[0151] Second exemplary embodiment

[0152] In the above, a single B-scan pair including the first B-scan 10 - 1 and the second B-scan 10 - 2 is processed by the data processing apparatus 100 to determine the first indication Z B-1 and / or second indication Z B-2 , and in some cases, also determining a first indicator indicating that the determination is reliable or a second indicator indicating that the determination is unreliable. However, the described data processing operations can be used to process multiple pairs of B-scans or C-scans, as in the present example embodiment, or more generally, multiple sets of two or more B-scans or C-scans.

[0153] Figure 16 is a flowchart showing a process by which the data processing apparatus 100 of this exemplary embodiment sequentially processes adjacent pairs of adjacent B scans 10-1 and 10-2 (eg, Figure 1 The series of B-scans 10 is shown as B-scans 10-1 and 10-2, followed by B-scans 10-3 and 10-4, etc., to generate ORG data indicative of the response of a common portion of the retina to an applied light stimulus.

[0154] Optical retinal imaging focuses on small changes in the retinal layers. The registration process of B-scans is often used to account for eye movements that would otherwise engulf the ORG response. Registration involves spatially shifting the two scans relative to each other so that they are located at the optimal similarity position defined by the maximum value of the cross-correlation function defined above. The maximum value of this cross-correlation function provides a good indication of how closely the two B-scans match. If they are not closely matched, the scans may result in unrepresentative results, affecting the final ORG results. B-scans with poor image quality can also adversely affect the ORG results. Therefore, it may be beneficial to avoid including velocity profiles from such B-scans in the analysis, as will now be explained.

[0155] exist Figure 16In process S110, the data processing device 100 calculates a corresponding comparison value for each of a plurality of different pairs of consecutive B-scans 10-1 and 10-2 in the series of B-scans 10. Each comparison value may be the value of an image quality metric calculated based on the amplitude component of at least one B-scan in the pair. The image quality metric (attribute) may take the form of, for example, a dynamic range, contrast, signal-to-noise ratio, or sharpness function value derived from the B-scan. A measure of the success of intensity-based segmentation of the B-scan is another example of an image quality metric that may be used. Alternatively, each comparison value may be the value of an image similarity metric that provides a measure of the degree of similarity between the images, the image quality metric being calculated based on the amplitude components of the B-scans in the pair. For example, the image similarity metric may take the form of the maximum value of the cross-correlation between the first B-scan 10-1 and the second B-scan 10-2 of the pair. Other image similarity metrics that may be used alternatively include sum of squared differences, mutual information, normalized mutual information, and Kullback-Leibler distance.

[0156] exist Figure 16 In process S120, the data processing apparatus 100 compares each comparison value with a threshold value to determine whether the comparison value is equal to or greater than the threshold value. In this exemplary embodiment, some pairs of B-scans produce comparison values ​​greater than or equal to the threshold value, and some other pairs of B-scans produce comparison values ​​less than the threshold value.

[0157] When Figure 16 The comparison value calculated in the process S120 is determined to be equal to or greater than the threshold value (in Figure 16 In the case of "Yes" at S130 of the data processing device 100, Figure 16 The phase components of a pair of B scans 10-1 and 10-2 are processed at S140, as described above with respect to Figure 3 The process S10 is described to generate a corresponding velocity distribution map P indicating the distribution.

[0158] On the other hand, when Figure 16 In the process S120 in the embodiment, it is determined that the calculated comparison value is not greater than or equal to the threshold value ( Figure 16 In the case of "No" in S130, Figure 16 In the process S150 , the data processing device 100 generates a velocity profile P, which indicates zero velocity v at all positions along the axial direction z in the velocity profile P, ie, a zero velocity profile.

[0159] exist Figure 16After processes S140 and S150 in the process, the data processing apparatus 100 determines in process S160 whether all adjacent pairs of B-scans in the series of B-scans 10 have been processed. If not ("No" at S160), the data processing apparatus 100 selects the next pair of adjacent B-scans in the series of B-scans 10 (i.e., a pair of B-scans adjacent to the last pair of B-scans that have been processed up to that point in the process), and then processes the data from the series of B-scans 10 as described above. Figure 16 The process S110 starts processing the selected pair of B-scans. However, if the data processing apparatus 100 determines in process S160 that all adjacent pairs of B-scans in a series of B-scans 10 have been processed (in Figure 16 If "Yes" is given at S160, the process proceeds to Figure 16 S180 , wherein the data processing device 100 concatenates the generated velocity profiles P such that the concatenation of the velocity profiles indicates how the distribution changes over time.

[0160] Figure 17 Shown in Figure 16 An example of a cascade 300 of velocity profiles P generated by the data processing device 100 in the process S180. Each velocity profile P is along Figure 17 The y-axis in FIG. 1 (labeled “OCT depth layer”) extends along the y-axis, and the velocity profiles calculated for adjacent pairs of B-scans in the series of B-scans 10 are plotted along the y-axis. Figure 17 The x-axis (labeled as "Time (ms)") in the diagram is arranged. Figure 17 The optical path length change ΔOPL (which provides an indication of velocity v) is shown as a function of position in each velocity profile P derived from a B-scan acquired after application of the optical stimulus at time t=200 ms. The optical path length change ΔOPL is shown as a function of position (i.e., along the Figure 17 The y-axis in FIG1 varies from a minimum ΔOPL of approximately −300 (arbitrary units) to a maximum ΔOPL of approximately +250 (arbitrary units). Although there is some variation in the distribution of ΔOPL with position between the velocity profiles P, particularly between the velocity profiles at times between 200 ms and 300 ms (shortly after application of the optical stimulus at time t=200 ms), each velocity profile still has a global maximum near OCT depth layer 11 associated with IS / OS and a global minimum near OCT depth layer 20 associated with COST.

[0161] Reference again Figure 16 In process S190, the data processing device 100 integrates corresponding parts of the velocity profile P in the cascade of velocity profiles, which parts have the same position along the axial direction z, to generate ORG data indicating a change in the optical path length at the position along the axial direction z over time. Figure 16 In process S190, the data processing device 100 can process each velocity distribution map P in the cascade set of velocity distribution maps by selecting a portion (segment or data element) of the velocity distribution map P set at a predetermined position along the axial direction z of the velocity distribution map P, and then integrating the selected (commonly located) portions by calculating the cumulative sum (or running total) of the velocity values ​​in these portions, which cumulative sum indicates how the optical path length of the optical path of the OCT sample beam ending at the position in the retina corresponding to the predetermined position in the velocity distribution map P changes with time.

[0162] By the first indication Z B-1 The indication of the change in optical path length over time at the first position along the axial direction z can be determined by integrating (i.e. calculating the cumulative sum or accumulated total) the corresponding parts of the cascade at the first position along the axial direction z. Additionally or alternatively, the indication of the change in optical path length over time at the first position along the axial direction z can be determined by integrating (i.e. calculating the cumulative sum or accumulated total) the corresponding parts of the cascade at the first position along the axial direction z. B-2 The indication of the change in optical path length over time at the second position along the axial direction z may be determined by integrating (ie calculating a cumulative sum or cumulative total) the corresponding portion of the cascade at the second position along the axial direction z.

[0163] Figure 18 A graph of the optical path length versus time ΔOPL is shown, which has been calculated by the data processing device 100 for selected positions in the velocity profiles in the cascade of velocity profiles. Figure 18 The graph in shows how the ΔOPL at a common location varies from one B-scan to the next in a series of B-scans 10 .

[0164] Figure 19 A graph showing the cumulative change in optical path length over time is referred to as an ORG signal (or ORG data). The ORG signal is obtained by calculating the cumulative sum of the optical path length changes for selected positions in the velocity profiles in the cascade of velocity profiles. Figure 19 As shown, the ORG signal becomes negative shortly after the stimulus is applied, then increases to become positive, and finally levels off at a value of approximately 250 nm in this example. Figure 19 It is also shown that the ORG signal has several plateaus, which are caused by replacing the velocity profile calculated using the B-scan with a velocity profile indicating zero velocity (i.e., a zero velocity profile) as described above, which did not produce a sufficiently high maximum cross-correlation value.

[0165] Therefore, in Figure 16The ORG data generated in process S190 may include one or more sets of equal consecutive values, and the data processing device 100 may be arranged to smooth the generated ORG data by replacing one or more values ​​in at least one of the one or more sets of equal consecutive values ​​with one or more estimated values, the one or more estimated values ​​being calculated based on adjacent values ​​adjacent to the set of equal consecutive values ​​in the ORG data. The data processing device 100 may perform this smoothing by using a moving median or moving average, for example, where the median or average of a specified number of points on either side of a missing data point is calculated and then assigned to the missing point. Alternatively, the data processing device 100 may, for example, detect each set of equal consecutive values ​​in the ORG data and replace the values ​​in each detected set with a corresponding estimated value obtained by (e.g., linearly) interpolating between a first value in the ORG data that is adjacent to the first of the equal consecutive values ​​in the detected set and a second value in the ORG data that is adjacent to the last of the equal consecutive values ​​in the detected set.

[0166] The result of smoothing the ORG signal using the moving average is as follows Figure 20 Before smoothing is applied, most of the data show that the OPL variation is zero, as shown in Figure 19 , and therefore generally underestimate the actual values ​​that will be observed. Using the method described herein, equivalent values ​​are replaced with estimated values ​​(e.g. Figure 20 ), which more accurately represents the changes in OPL.

[0167] Although the data processing device 100 processes adjacent pairs of B-scans in a series of B-scans 10 (i.e., one pair of B-scans adjacent to another pair of B-scans in the series) sequentially to generate the velocity profile P, the data processing device 100 may alternatively process adjacent pairs of B-scans in parallel to generate the velocity profile P, thereby significantly speeding up the processing. Furthermore, although pairs of B-scans adjacent to each other in a series of B-scans 10 are processed (i.e., consecutive B-scans in the series), the data processing device 100 may alternatively process pairs of B-scans in the series where, for example, the B-scans of each pair are separated from each other by one or more intermediate B-scans.

[0168] Furthermore, although the data processing apparatus 100 of this example embodiment has been described as processing pairs of B-scans, it may more generally be configured to process a set of more than two B-scans in a series of B-scans 10, which may be consecutive B-scans in the series or individual B-scans in the series separated from each other by one or more intermediate B-scans that do not form part of the set. Figure 165) B-scans used in a loop of the process, wherein in each loop of the process the window selection function selects, for example, all B-scans in a window portion of a series of B-scans 10 or every other B-scan in a window portion, and the selected sets may have one or more B-scans in common. In such other example embodiments, the above reference to Figure 16 A modified version of the described procedure will be as follows, which is based on the technique described in Kari V. Vienola et al., “Velocity-based optoretinography for clinical applications,” Optica 9, pp. 1100-1108, 2022.

[0169] exist Figure 16 In a modified form of process S110 in FIG. 1 , the data processing apparatus 100 calculates a comparison value for one of a plurality of different sets of consecutive B-scans in the series of B-scans 10. The comparison value may take any of the different forms described above. Then, as in Figure 16 In the process S120, the data processing apparatus 100 compares the comparison value with the threshold value to determine whether the comparison value is equal to or greater than the threshold value. Some sets of B scans produce comparison values ​​greater than or equal to the threshold value, and some other sets of B scans produce comparison values ​​less than the threshold value. In the case where the calculated comparison value is determined to be equal to or greater than the threshold value, the data processing apparatus 100 performs the comparison according to Figure 16 A modified form of process S140 processes the phase components of at least some of the B-scans in the set to generate a velocity profile P indicative of a distribution. Here, the data processing apparatus 100 may first flatten the B-scans of the set so that the IS / OS and COST reflections are located at substantially the same height for each A-scan in each B-scan. The data processing apparatus 100 may then register the B-scans relative to each other. The phase data cube θ(x, z, t) is then flattened by the phase data cube θ(x, z, t) at θ(x p ,z q ,t r ) by adding or subtracting 2π to minimize |θ(x p ,z q ,t r )-θ(x p ,z q ,t r-1 )| to expand in the time dimension, where t r and t r-1 represents a continuous phase B scan. For each spatial coordinate pair (x p ,z q ) performs this step. The phase change rate for each coordinate pair is calculated by performing a least squares linear fit with respect to t, given The unit is rad / s. Therefore, the instantaneous velocity of each spatial location can be calculated as where λ is the wavelength of the OCT light used, and n is the nominal refractive index of the eye 20. These instantaneous velocities can be averaged in the lateral dimension (i.e., along the x-axis) to give an instantaneous, depth-dependent measure of velocity along the z-axis, i.e. Figure 1 A one-dimensional velocity profile P of the type shown schematically in FIG is indicative of the distribution of velocities v in the axial direction (z-axis) within a common portion of the retina covered by repeated B-scans. The B-scan amplitudes may also be averaged over the lateral dimension (x-axis) if desired to give an instantaneous, depth-dependent measure of backscatter. It should be noted that averaging over the lateral dimension is not necessary, and a two-dimensional velocity profile may instead be generated that indicates the distribution of velocities v in both the axial direction (z-axis) and the lateral direction (x-axis) within a common portion of the retina covered by repeated B-scans. On the other hand, in the event that the calculated comparison value is determined to be not equal to or greater than the threshold value, as in Figure 16 In the process S150, the data processing device 100 generates a velocity profile P, which indicates zero velocity v at all positions along the axial direction z in the velocity profile P, i.e., a zero velocity profile. The data processing device 100 of the other exemplary embodiments mentioned above will perform the above-mentioned reference Figure 16 Another process described, but refer to Figure 16 The "pairs" of B-scans in processes S160 and 170 are replaced by the "sets" of B-scans discussed above.

[0170] The data processing apparatus 100 may more generally be arranged to perform what will now be referred to as Figure 21 Describe the process.

[0171] exist Figure 21 In process S210, the data processing apparatus 100 calculates a corresponding comparison value for each of a plurality of different sets of two or more OCT images in the series of OCT images 10. The corresponding comparison value may be a corresponding value of an image quality metric calculated based on an amplitude component of at least one OCT image in the set, or a corresponding value of an image similarity metric that provides a measure of the degree of similarity between images, the value of the image similarity metric being calculated based on the amplitude components of at least two OCT images in the set. The comparison value may take any of the above-described exemplary forms.

[0172] exist Figure 21 In process S220 , the data processing apparatus 100 compares each comparison value with a threshold value to determine whether the comparison value is equal to or greater than the threshold value.

[0173] For each set of the plurality of different sets of OCT images where the calculated comparison value is determined to be equal to or greater than the threshold value, the data processing apparatus 100 Figure 21 In process S230 , the phase components of the OCT images in the set are processed to generate corresponding velocity profiles indicating the distribution.

[0174] For each set of the plurality of different sets of OCT images where the calculated comparison value is determined not to be greater than or equal to the threshold value, the data processing apparatus 100 determines that: Figure 21 In the process S240 , a corresponding velocity profile P is generated, which indicates zero velocity v at all positions in the velocity profile P along the axial direction z.

[0175] exist Figure 21 In the process S250 , the data processing apparatus 100 generates a concatenation of the generated velocity profiles such that the concatenation of the velocity profiles indicates how the profile changes over time.

[0176] Finally, in Figure 21 In process S260, the data processing apparatus 100 integrates corresponding portions of the concatenation of velocity profiles, the portions having the same position along the axial direction, to generate data indicating a temporal change in the optical path length at the position along the axial direction.

[0177] [Third exemplary embodiment]

[0178] In the second example embodiment described above, the data processing apparatus 100 may reduce quality degradation of the generated ORG data by using a comparison value calculated for a set of two or more B-scans in a series of B-scans 10 to set the velocity profiles used to generate the ORG data and replacing any velocity profiles derived from a set in which one or more of the component B-scans had unacceptable image quality or in which two or more of the component B-scans had an insufficient degree of similarity with a velocity profile indicating zero velocity throughout (i.e., a zero velocity profile), and optionally suppressing any resulting artifacts in the ORG data by smoothing the ORG data. However, as in the present example embodiment, the data processing apparatus 100 may alternatively be arranged to reduce quality degradation of the generated ORG data by replacing any velocity profiles derived from a set in which one or more of the component B-scans had unacceptable image quality or in which two or more of the component B-scans had an insufficient degree of similarity with an estimated (rather than zero) velocity profile, as will now be referred to. Figure 22 Described in more detail.

[0179] Figure 22 The processes S310 to S330, S350 and S360 are respectively the same as those described in detail above. Figure 21 The processes S210 to S230, S250 and S260 are the same.

[0180] exist Figure 22 In the process S340, the data processing device 100 generates a corresponding estimated velocity profile P for each set of the plurality of different sets of B-scans for which the calculated comparison value is determined to be not greater than or equal to the threshold value (i.e., less than the threshold value), the estimated velocity profile P indicating the distribution of velocity v along the axial direction (z-axis) within the common portion of the retina covered by the B-scans. Each estimated velocity profile is based on the respective Figure 22 The process S330 calculates one or more velocity profiles and compares them with the velocity profiles in Figure 22 In other words, the data processing apparatus 100 of the present exemplary embodiment generates, for each unreliable set of two or more B-scans whose calculated comparison values ​​are determined not to be greater than or equal to the threshold value, a corresponding velocity profile P based on an adjacent reliable set of two or more B-scans whose calculated comparison values ​​are determined to be equal to or greater than the threshold value, rather than based on a given set of two or more B-scans. In the present exemplary embodiment, in Figure 22 The ORG data generated in process S360 does not have Figure 19 The plateau shown, and therefore smoothing of the ORG data may not be required.

[0181] Some of the example embodiments described above are summarized in the following numbered clauses E1 to E18.

[0182] E1. A data processing apparatus 100 arranged to process a phase component of a first OCT image 10-1 and a phase component of a second OCT image 10-2 in a series of OCT images 10 of a common portion of the retina of an eye 20 acquired by a Fourier domain optical coherence tomography (OCT) imaging system 30 after the common portion of the retina is stimulated by an optical stimulus 40, wherein the common portion includes a layer L of the retina whose thickness changes in response to the optical stimulus 40 during acquisition of the series of OCT images 10, to determine a first indication Z of a first position of a first boundary B-1 of the layer L along an axial direction z in the OCT images 10. B-1 and a second indication Z of a second position of the second boundary B-2 of the layer L along the axial direction z in the OCT image 10 B-2 At least one of the data processing apparatus 100 is arranged to:

[0183] processing the phase component of the first OCT image 10 - 1 and the phase component of the second OCT image 10 - 2 to calculate a velocity profile P indicating a distribution of velocity v along an axial direction z within a common portion of the retina; and

[0184] Determine the first indication Z based on the calculated velocity profile B-1 and the second indication Z B-2 At least one of .

[0185] E2. The data processing device 100 according to E1, wherein the data processing device 100 is arranged to determine the first indication Z based on the calculated speed profile by: B-1 and the second indication Z B-2 At least one of:

[0186] In determining the first indication Z B-1 In the case of the velocity profile P, the maximum value v corresponding to the velocity indicated by the velocity profile P is determined. max The first position z max As the first indication Z B-1 ;and

[0187] In determining the second indication Z B-2 In the case of , determining the minimum value v in the velocity profile P corresponding to the velocity indicated by the velocity profile P min The second position z min As the second indicator Z B-2 .

[0188] E3. The data processing device 100 according to E1 or E2, wherein the data processing device 100 is arranged to determine the first indication Z by processing the calculated velocity profile P using a cluster analysis algorithm or a dimensionality reduction algorithm B-1 and the second indication Z B-2 At least one of .

[0189] E4. The data processing apparatus 100 according to E3, wherein the dimensionality reduction algorithm comprises one of a principal component analysis (PCA) algorithm, an independent component analysis (ICA) algorithm, a linear discriminant analysis (LDA) algorithm, or a non-negative matrix factorization (NMF) algorithm.

[0190] E5. The data processing apparatus 100 according to E4, wherein the dimensionality reduction algorithm comprises a PCA algorithm, and the first indication Z B-1 and the second indication Z B-2 At least one of is determined using principal components that have been determined by applying a PCA algorithm to the calculated velocity profile P.

[0191] E6. The data processing device 100 according to E2, wherein the data processing device 100 is arranged to determine the first position z in the velocity profile P by: max and the second position z in the velocity distribution diagram P min :

[0192] calculating, for each combination of the candidate first position in the velocity profile P and the candidate second position in the velocity profile P among all combinations of the candidate first position and the candidate second position in the velocity profile P, a respective value of a difference between a respective speed indicated by the velocity profile (P) for the candidate first position and a respective speed indicated by the velocity profile (P) for the candidate second position in the combination; and

[0193] The candidate first position and candidate second position of the following combination among the plurality of combinations are identified as the first position z in the velocity profile max and the second position z in the velocity profile P min : For the combination among multiple combinations, the calculated value of the difference is the largest among the multiple calculated values ​​of the difference.

[0194] E7. The data processing device 100 according to E2, wherein the data processing device 100 is arranged to determine the first position z in the velocity profile by: max and the second position z in the velocity profile P min :

[0195] calculating, for each combination of the candidate first position in the speed profile P and the candidate second position in the speed profile P among all combinations of the candidate first position and the candidate second position in the speed profile P, a respective value of a difference between an average value of respective speeds indicated by the speed profile P for a set of neighboring positions including the candidate first position in the combination and an average value of respective speeds indicated by the speed profile P for a set of neighboring positions including the candidate second position in the combination; and

[0196] The candidate first position and candidate second position of the following combination among the plurality of combinations are identified as the first position z in the velocity profile max and the second position z in the velocity profile P min : For the combination among multiple combinations, the calculated value of the difference is the largest among the multiple calculated values ​​of the difference.

[0197] E8. The data processing apparatus 100 according to any one of E1 to E7, wherein the data processing apparatus 100 is further arranged to perform one of the following:

[0198] A first position z of a first boundary B-1 of the layer L along the axial direction z in the OCT image is determined based on the amplitude component of at least one OCT image in the series of OCT images 10. max The first indication Z B-1 , where the second indicator Z B-2By identifying a second position z in the velocity profile P corresponding to a minimum value of the velocity indicated by the velocity profile P min to determine; and

[0199] A second position z of the second boundary B-2 of the layer L along the axial direction z in the OCT image is determined based on the amplitude component of at least one OCT image in the series of OCT images 10. min The second instruction Z B-2 , where the first indicator Z B-1 By identifying the maximum value v of the velocity indicated by the velocity profile P max The corresponding first position z in the velocity distribution diagram P max to confirm.

[0200] E9. The data processing apparatus 100 according to any one of E1 to E8, wherein the data processing apparatus 100 is further arranged to:

[0201] Computes a comparison value, which is one of the following:

[0202] a value of an image quality metric calculated based on at least one of the magnitude component of the first OCT image 10 - 1 and the magnitude component of the second OCT image 10 - 2 ; and

[0203] Providing a value of an image similarity metric that is a measure of a degree of similarity between the images, the value of the image similarity metric being calculated based on an amplitude component of the first OCT image 10 - 1 and an amplitude component of the second OCT image 10 - 2 ;

[0204] comparing the comparison value with a threshold value to determine whether the comparison value is equal to or greater than the threshold value;

[0205] generating a first indicator indicating that at least one of the first indication and the second indication is reliable in a case where the comparison value is determined to be equal to or greater than the threshold value; and

[0206] In a case where it is determined that the comparison value is not equal to or greater than the threshold value, a second indicator is generated indicating that at least one of the first indication and the second indication is unreliable.

[0207] E10. The data processing apparatus 100 according to E9, wherein the comparison value is a maximum value of the calculated cross-correlations between the first OCT image 10 - 1 and the second OCT image 10 - 2 .

[0208] E11. The data processing apparatus 100 according to any one of E1 to E8, wherein the data processing apparatus is further arranged to:

[0209] For each of a plurality of different sets of two or more OCT images 10 - 1 , 10 - 2 in the series of OCT images 10 , a corresponding comparison value is calculated as one of:

[0210] a respective value of an image quality metric calculated based on an amplitude component of at least one of the OCT images in the set; and

[0211] providing a corresponding value of an image similarity metric that is a measure of a degree of similarity between the images, the value of the image similarity metric being calculated based on amplitude components of at least two of the OCT images in the set;

[0212] comparing each comparison value to a threshold value to determine whether the comparison value is equal to or greater than the threshold value;

[0213] for each of a plurality of different sets of OCT images, for a set in which the calculated comparison value is determined to be equal to or greater than a threshold, processing phase components of the OCT images in the set to generate a corresponding velocity profile indicative of a distribution;

[0214] for each of the plurality of different sets of OCT images, for the sets in which the calculated comparison value is determined not to be greater than or equal to the threshold, generating a corresponding velocity profile P indicating zero velocity at all positions along the axial direction z in the velocity profile P;

[0215] generating a cascade 300 of generated velocity profiles P such that the cascade 300 of velocity profiles indicates how the profile changes over time; and

[0216] Corresponding portions of the concatenation 300 of velocity profiles P having the same position along the axial direction z are integrated to generate data indicative of the optical path length change over time at a position along the axial direction z.

[0217] E12. The data processing apparatus 100 according to E11, wherein

[0218] The generated data includes one or more sets of equal consecutive values, and

[0219] The data processing apparatus 100 is further arranged to smooth the generated data by replacing one or more values ​​of one of the one or more sets of equal consecutive values ​​with one or more estimated values, the one or more estimated values ​​being calculated based on adjacent values ​​adjacent to the set of equal consecutive values ​​in the generated data.

[0220] E13. The data processing apparatus 100 according to any one of E1 to E8, wherein the data processing apparatus is further arranged to:

[0221] For each of a plurality of different sets of two or more OCT images 10 - 1 , 10 - 2 in the series of OCT images 10 , a corresponding comparison value is calculated as one of:

[0222] a respective value of an image quality metric calculated based on an amplitude component of at least one of the OCT images in the set; and

[0223] providing a corresponding value of an image similarity metric that is a measure of a degree of similarity between the images, the value of the image similarity metric being calculated based on amplitude components of at least two of the OCT images in the set;

[0224] comparing each comparison value to a threshold value to determine whether the comparison value is equal to or greater than the threshold value;

[0225] for each of a plurality of different sets of OCT images, for a set in which the calculated comparison value is determined to be equal to or greater than a threshold, processing phase components of the OCT images in the set to generate a corresponding velocity profile indicative of a distribution;

[0226] for each of the plurality of different sets of OCT images, for a set for which the calculated comparison value is determined not to be greater than or equal to the threshold, generating a respective velocity profile indicative of a distribution based on respective velocity profiles calculated for each of one or more other sets of the plurality of different sets of OCT images;

[0227] generating a cascade 300 of generated velocity profiles such that the cascade 300 of velocity profiles indicates how the profile changes over time; and

[0228] Corresponding portions of the concatenation 300 of velocity profiles are integrated to generate data indicative of the change in optical path length over time at positions along the axial direction, the portions having the same position along the axial direction.

[0229] E14. The data processing apparatus 100 according to any one of E11 to E13, wherein the data processing apparatus 100 is arranged to calculate a respective maximum value of the cross-correlation between at least two of the OCT images in each set as a respective comparison value for the set.

[0230] E15. The data processing device 100 according to any one of E11 to E14, wherein the data processing device 100 is arranged to B-1 and the second indication Z B-2A corresponding partial integration of the cascade 300 at at least one of the indicated first positions along the axial direction z and the second positions along the axial direction z determines a corresponding indication of the change in optical path length over time at each of the at least one of the first positions along the axial direction z and the second positions along the axial direction z.

[0231] E16. The data processing apparatus 100 according to any one of E1 to E15, wherein the layers of the retina comprise outer segments of photoreceptor cells.

[0232] E17. The data processing apparatus 100 according to any one of E1 to E16, wherein the data processing apparatus 100 is further arranged to:

[0233] Before calculating the velocity distribution map P, the phase component of the first OCT image 10-1 and the phase component of the second OCT image 10-2 are processed to compensate for the overall movement of the common portion during the acquisition of a series of OCT images 10 by the Fourier domain OCT imaging system 30 after the common portion of the retina is stimulated by the optical stimulation 40.

[0234] E18. A system 1000 comprising:

[0235] an optical stimulus source 45 arranged to provide the optical stimulus 40;

[0236] a Fourier domain optical coherence tomography imaging system 30 arranged to acquire a series of OCT images 10 of a common portion of the retina of the eye 20 after stimulation of the common portion by an optical stimulus 40; and

[0237] The data processing apparatus 100 according to any one of E1 to E17 is arranged to process phase components of a first OCT image 10 - 1 and a second OCT image 10 - 2 of a series of OCT images 10 acquired by a Fourier domain optical coherence tomography imaging system 30 .

[0238] In the foregoing description, example aspects have been described with reference to several example embodiments. Therefore, the description should be regarded as illustrative rather than restrictive. Similarly, the figures shown in the accompanying drawings, which highlight the features and advantages of the example embodiments, are presented for illustrative purposes only. The architecture of the example embodiments is sufficiently flexible and configurable that it can be utilized in ways other than those shown in the accompanying drawings.

[0239] In one example embodiment, some aspects of the examples given herein, such as reference Figure 3 、 Figure 13 、 Figure 16 、 Figure 21 and Figure 22The processing methods described herein may be provided as a computer program or software, such as one or more programs having instructions or instruction sequences, contained in or stored on an article of manufacture, such as a machine-accessible or machine-readable medium, instruction storage, or computer-readable storage device, each of which may be non-transitory. The program or instructions on the non-transitory machine-accessible medium, machine-readable medium, instruction storage, or computer-readable storage device may be used to program a computer system or other electronic device. Machine-readable or computer-readable media, instruction storage, and storage devices may include, but are not limited to, floppy disks, optical disks, magneto-optical disks, or other types of media / machine-readable media / instruction storage / storage devices suitable for storing or transmitting electronic instructions. The techniques described herein are not limited to any particular software configuration. They may be employed in any computing or processing environment. As used herein, the terms "computer-readable," "machine-accessible medium," "machine-readable medium," "instruction storage," and "computer-readable storage device" shall include any medium capable of storing, encoding, or transmitting instructions or instruction sequences for execution by a machine, computer, or computer processor, and causing the machine / computer / computer processor to perform any of the methods described herein. Furthermore, it is common in the art to refer to software in one form or another (e.g., program, procedure, process, application, module, unit, logic, etc.) as taking an action or causing a result. Such expressions are merely a shorthand way of stating that execution of the software by a processing system causes the processor to perform an action to produce a result.

[0240] Some or all of the functionality of the OCT data processing hardware 130 may also be implemented by preparing application specific integrated circuits, field programmable gate arrays, or by interconnecting an appropriate network of conventional component circuits.

[0241] The computer program product may be provided in the form of one or more storage media, instruction storage means, or storage devices having stored thereon or therein instructions that can be used to control or cause a computer or computer processor to perform any of the processes of the example embodiments described herein. The storage media / instruction storage means / storage devices may include, by way of example and without limitation, optical disks, ROM, RAM, EPROM, EEPROM, DRAM, VRAM, flash memory, flash memory cards, magnetic cards, optical cards, nanosystems, molecular memory integrated circuits, RAID, remote data storage / archiving / warehousing devices, and / or any other type of device suitable for storing instructions and / or data.

[0242] For storage on any one of one or more computer-readable media, instruction storage devices, or storage devices, some embodiments include hardware for controlling the system and software for enabling the system or microprocessor to utilize the results of the example embodiments described herein to interact with a human user or other mechanism. Such software may include, without limitation, device drivers, operating systems, and user applications. Ultimately, as described above, such computer-readable media or storage devices also include software for performing example aspects of the present invention.

[0243] Software modules for implementing the processes described herein are included in the system's programming and / or software. In some example embodiments herein, the modules include software, but in other example embodiments herein, the modules include hardware or a combination of hardware and software.

[0244] Although various exemplary embodiments of the present invention have been described above, it should be understood that they have been presented by way of example and not limitation. It will be apparent to those skilled in the relevant art that various changes in form and detail may be made. Therefore, the present invention should not be limited by any of the above exemplary embodiments, but should be defined only in accordance with the appended claims and their equivalents.

[0245] Although this specification contains many specific embodiment details, these should not be understood as limitations on the scope of any invention or what may be claimed, but rather as descriptions of features specific to the particular embodiments described herein. Certain features described in this specification in the context of separate embodiments may also be implemented in combination in a single embodiment. Conversely, various features described in the context of a single embodiment may also be implemented in multiple embodiments individually or in any suitable subcombination. Furthermore, although features may be described above as acting in a particular combination and even initially claimed as such, one or more features from the claimed combination may in some cases be deleted from the combination, and the claimed combination may be directed to subcombinations or variations of subcombinations.

[0246] In some cases, multitasking and parallel processing may be advantageous. In addition, the separation of the various components in the above embodiments should not be understood as requiring such separation in all embodiments, and it should be understood that the described program components and systems can generally be integrated together in a single software product or packaged into multiple software products.

[0247] Now that some illustrative embodiments and implementations have been described, it will be apparent that the foregoing embodiments are illustrative rather than restrictive and have been presented by way of example. In particular, although many of the examples presented herein involve specific combinations of devices or software elements, these elements can be combined in other ways to achieve the same objectives. Actions, elements, and features discussed in connection with only one embodiment are not intended to be excluded from similar roles in that embodiment or other embodiments.

Claims

1. A computer-implemented method for processing a phase component of a first OCT image (10-1) and a phase component of a second OCT image (10-2) in a series of OCT images (10) of a common portion of a retina of an eye (20) acquired by a Fourier domain optical coherence tomography (OCT) imaging system (30) after the common portion of the retina of the eye (20) is stimulated by an optical stimulus (40), wherein: The common portion includes a layer (L) of the retina, the thickness of the layer (L) changing in response to the optical stimulus (40) during acquisition of a series of the OCT images (10) to determine a first indication (Z) of a first position of a first boundary (B-1) of the layer (L) along an axial direction (z) in the OCT images (10). B-1 ) and a second indication (Z) of a second position of a second boundary (B-2) of the layer (L) along the axial direction (z) in the OCT image (10) B-2 ), the method comprising: processing (S10) a phase component of the first OCT image (10-1) and a phase component of the second OCT image (10-2) to calculate a velocity profile (P) indicating a distribution of velocity (v) within the common portion of the retina along the axial direction (z); and The first indication (Z) is determined (S20) based on the calculated velocity profile (P) B-1 ) and the second indication (Z B-2 ) at least one of .

2. The computer-implemented method of claim 1 , wherein: The first indication (Z B-1 ) and the second indication (Z B-2 ) is determined based on the calculated velocity profile (P) in the following manner: In determining the first indication (Z B-1 ), determining a first position in the velocity profile (P) corresponding to a maximum value of the velocity (v) indicated by the velocity profile (P) as a first indication (Z B-1 ); and In determining the second indication (Z B-2 ), determining a second position in the velocity profile (P) corresponding to a minimum value of the velocity (v) indicated by the velocity profile (P) as a second indication (Z B-2 ).

3. The computer-implemented method of claim 1 or claim 2, wherein: The first indication (Z B-1 ) and the second indication (Z B-2 ) is determined by processing the calculated velocity distribution map (P) using a cluster analysis algorithm or a dimensionality reduction algorithm.

4. The computer-implemented method of claim 3, wherein: The dimensionality reduction algorithm includes one of a principal component analysis algorithm, an independent component analysis algorithm, a linear discriminant analysis algorithm or a non-negative matrix decomposition algorithm.

5. The computer-implemented method of claim 4, wherein: The dimensionality reduction algorithm includes a principal component analysis algorithm, and the first indication (Z B-1 ) and the second indication (Z B-2 ) are determined using principal components that have been determined by applying the principal component analysis algorithm to the calculated velocity profile (P).

6. The computer-implemented method of claim 2, wherein: The first position in the velocity profile (P) and the second position in the velocity profile (P) are both determined by: calculating, for each combination of a candidate first position in the velocity profile (P) and a candidate second position in the velocity profile (P) among all combinations of candidate first positions and candidate second positions in the velocity profile (P), a respective value of a difference between a respective speed (v) indicated by the velocity profile (P) for the candidate first position and a respective speed (v) indicated by the velocity profile (P) for the candidate second position in the combination; and A candidate first position and a candidate second position of the following combination among multiple combinations are identified as the first position in the speed distribution map (P) and the second position in the speed distribution map (P): for this combination among the multiple combinations, the calculated value of the difference is the largest among the multiple calculated values ​​of the difference.

7. The computer-implemented method of claim 2, wherein: The first position in the velocity profile (P) and the second position in the velocity profile (P) are both determined by: calculating, for each combination of a candidate first position in the speed profile (P) and a candidate second position in the speed profile (P) among all combinations of candidate first positions and candidate second positions in the speed profile (P), a respective value of a difference between an average value of respective speeds indicated by the speed profile (P) for a set of neighboring positions including the candidate first position in the combination and an average value of respective speeds indicated by the speed profile (P) for a set of neighboring positions including the candidate second position in the combination; and The candidate first position and the candidate second position of the following combination among multiple combinations are identified as the first position in the speed distribution map and the second position in the speed distribution map (P): for this combination among the multiple combinations, the calculated value of the difference is the largest among the multiple calculated values ​​of the difference.

8. The computer-implemented method of any one of claims 1 to 7, further comprising one of the following: determining a first indication (Z) of a first position of the first boundary (B-1) of the layer (L) along the axial direction (z) in the OCT image (10) based on an amplitude component of at least one OCT image in a series of the OCT images (10); B-1 ), wherein the second indication (Z B-2 ) is determined by identifying said second position in said velocity profile (P) corresponding to a minimum value of the velocity indicated by said velocity profile (P); and determining a second indication (Z) of a second position of the second boundary (B-2) of the layer (L) along the axial direction (z) in the OCT image (10) based on an amplitude component of at least one OCT image in a series of the OCT images (10); B-2 ), wherein the first indication (Z B-1 ) is determined by identifying said first position in said velocity profile (P) corresponding to a maximum value of the velocity indicated by said velocity profile (P).

9. The computer-implemented method of any one of claims 1 to 8, further comprising: Compute a comparison value, which is one of the following: a value of an image quality metric calculated based on at least one of an amplitude component of the first OCT image (10-1) and an amplitude component of the second OCT image (10-2); and Providing a value of an image similarity metric that is a measure of the degree of similarity between the images, wherein the value of the image similarity metric is calculated based on the amplitude component of the first OCT image (10-1) and the amplitude component of the second OCT image (10-2); comparing the comparison value with a threshold value to determine whether the comparison value is equal to or greater than the threshold value; In the case where it is determined that the comparison value is equal to or greater than the threshold, a first indicator is generated, wherein the first indicator indicates that the first indication (Z B-1 ) and the second indication (Z B-2 ) is reliable; and In the case where it is determined that the comparison value is not equal to or greater than the threshold, a second indicator is generated, wherein the second indicator indicates that the first indicator (Z B-1 ) and the second indication (Z B-2 ) is unreliable.

10. The computer-implemented method of claim 9, wherein: The comparison value is a maximum value of the calculated cross-correlations between the first OCT image (10-1) and the second OCT image (10-2).

11. The computer-implemented method of any one of claims 1 to 8, further comprising: For each of a plurality of different sets of two or more OCT images (10-1, 10-2) in the series of OCT images (10), a corresponding comparison value is calculated as one of: a respective value of an image quality metric calculated based on an amplitude component of at least one of the OCT images in the set; and providing a corresponding value of an image similarity metric that is a measure of a degree of similarity between the images, the value of the image similarity metric being calculated based on amplitude components of at least two of the OCT images in the set; comparing each comparison value with a threshold value to determine whether the comparison value is equal to or greater than the threshold value; for each of a plurality of different sets of OCT images, for the set for which the calculated comparison value is determined to be equal to or greater than the threshold, processing phase components of the OCT images in the set to generate a corresponding velocity profile indicative of the distribution; generating, for each of a plurality of different sets of OCT images, for the set in which the calculated comparison value is determined not to be greater than or equal to the threshold, a corresponding velocity profile (P) indicating zero velocity (v) at all locations along the axial direction (z) in the velocity profile (P); generating a cascade (300) of the generated velocity profiles such that the cascade (300) of velocity profiles indicates how the profile changes over time; and Respective portions of the concatenation (300) of velocity profiles are integrated to generate data indicative of optical path length variations over time at positions along the axial direction, the portions having the same position along the axial direction.

12. The computer-implemented method of claim 11, wherein: The generated data includes one or more sets of equal consecutive values, and The method also includes smoothing the generated data by replacing one or more values ​​in one of the one or more sets of equal consecutive values ​​with one or more estimated values, wherein the one or more estimated values ​​are calculated based on adjacent values ​​adjacent to the set of equal consecutive values ​​in the generated data.

13. The computer-implemented method of any one of claims 1 to 8, further comprising: For each of a plurality of different sets of two or more OCT images (10-1, 10-2) in the series of OCT images (10), a corresponding comparison value is calculated as one of: a respective value of an image quality metric calculated based on an amplitude component of at least one of the OCT images in the set; and providing a corresponding value of an image similarity metric that is a measure of a degree of similarity between the images, the value of the image similarity metric being calculated based on amplitude components of at least two of the OCT images in the set; comparing each comparison value with a threshold value to determine whether the comparison value is equal to or greater than the threshold value; for each of a plurality of different sets of OCT images, for the set for which the calculated comparison value is determined to be equal to or greater than the threshold, processing phase components of the OCT images in the set to generate a corresponding velocity profile indicative of the distribution; for each of the plurality of different sets of OCT images, for the set for which the calculated comparison value is determined not to be greater than or equal to the threshold, generating a respective velocity profile indicative of the distribution based on respective velocity profiles calculated for each of one or more other sets of the plurality of different sets of OCT images; generating a cascade (300) of the generated velocity profiles such that the cascade (300) of velocity profiles indicates how the profile changes over time; and Respective portions of the concatenation (300) of velocity profiles are integrated to generate data indicative of optical path length variations over time at positions along the axial direction, the portions having the same position along the axial direction.

14. The computer-implemented method of any one of claims 11 to 13, wherein: The respective comparison value calculated for each set is a respective maximum value of the calculated cross-correlations between at least two of the OCT images in the set.

15. The computer-implemented method of any one of claims 11 to 14, wherein: By the first indication (Z B-1 ) and the second indication (Z B-2 ) is determined by integrating corresponding portions of the cascade (300) at at least one of the first position along the axial direction (z) and the second position along the axial direction (z).

16. The computer-implemented method of any one of claims 1 to 15, wherein: The layers of the retina include the outer segments of photoreceptor cells.

17. The computer-implemented method according to any one of claims 1 to 16, further comprising processing the phase component of the first OCT image (10-1) and the phase component of the second OCT image (10-2) before calculating the velocity distribution map (P) to compensate for the overall motion of the common portion during acquisition of a series of the OCT images (10) by the Fourier domain OCT imaging system (30) after the common portion of the retina is stimulated by the optical stimulus (40).

18. A computer program (245) comprising computer readable instructions which, when executed by a processor (220), cause the processor (200) to perform the method according to any one of claims 1 to 17.

19. A data processing device (100) arranged to process a phase component of a first OCT image (10-1) and a phase component of a second OCT image (10-2) in a series of OCT images (10) of a common portion of a retina of an eye (20) acquired by a Fourier domain optical coherence tomography (OCT) imaging system (30) after the common portion of the retina of the eye (20) is stimulated by an optical stimulus (40), wherein: The common portion includes a layer (L) of the retina, the thickness of the layer (L) changing in response to the optical stimulus (40) during acquisition of a series of the OCT images (10) to determine a first indication (Z) of a first position of a first boundary (B-1) of the layer (L) along an axial direction (z) in the OCT images (10). B-1 ) and a second indication (Z) of a second position of a second boundary (B-2) of the layer (L) along the axial direction (z) in the OCT image (10) B-2 ), the data processing device (100) being arranged to: processing (S10) a phase component of the first OCT image (10-1) and a phase component of the second OCT image (10-2) to calculate a velocity profile (P) indicating a distribution of velocity (v) within the common portion of the retina along the axial direction (z); and The first indication (Z) is determined (S20) based on the calculated velocity profile (P) B-1 ) and the second indication (Z B-2 ) at least one of .

Images