Processing technique for optical retina imaging

Abstract

To provide a computer-implemented method for detecting a physiological response of the retina of an object eye to an optical stimulus.SOLUTION: There is provided a computer-implemented method for processing respective phase-components of a first OCT image and a second OCT image of a sequence of OCT images of a common portion of the retina of an eye acquired by a Fourier domain OCT imaging system after stimulation of the common portion by an optical stimulus, the common portion including a layer of the retina whose thickness has changed during acquisition of the sequence of OCT images, and determining an index of a position along an axial direction in the OCT images of a boundary of the layer, the method comprising: processing the phase-components of the first OCT image and the phase-components of the second OCT image (S10) to calculate a velocity-profile indicative of a distribution of velocities along the axis in the common portion of the retina; and determining the index on the basis of the calculated velocity-profile (S20).SELECTED DRAWING: Figure 3

Description

[Technical Field]

[0001] Exemplary aspects of the present specification relate generally to the field of optical coherence tomography (OCT), and more particularly to techniques for processing OCT data generated by a Fourier-domain OCT imaging system to generate optical retinal imaging (optoretinography) data indicative of the physiological response of the retina of a subject's eye to optical stimuli. [Background technology]

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

[0003] Depending on how the depth range is achieved, OCT imaging systems can be classified as time-domain OCT (TD-OCT) or Fourier-domain OCT (FD-OCT) (also called frequency-domain OCT). In TD-OCT, the optical path length of the reference arm of the imaging system's interferometer is varied in time during acquisition of a reflectance profile of the scattering medium (referred to herein as the "imaging target") being imaged by the OCT imaging system; the reflectance profile is commonly referred to as a "depth scan" or "axial scan" ("A-scan"). In FD-OCT, the spectral interferogram resulting from the interference of light in the reference arm with light in the sample arm of the interferometer at each A-scan position is Fourier transformed to simultaneously acquire all points along the depth of the A-scan without requiring changes to the optical path length of the reference arm. In FD-OCT, all backreflections from the sample are measured simultaneously, enabling much faster imaging than scanning the sample arm mirror in an interferometer. 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 the detector to measure all wavelengths simultaneously. In SS-OCT (also called time-encoded frequency-domain OCT), the light source is swept across a range of wavelengths, and the temporal output of the detector is converted into spectral interference.

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

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

[0006] OCT imaging systems can also be classified as phase-resolved, in which 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 often have a degree of phase stability that allows them to function as phase-resolved OCT imaging systems.

[0007] Optical retinal imaging (ORG) generally refers to the detection of the physiological response of the retina of the eye to optical stimuli (i.e., functional activity of the retina induced by light). ORG techniques involve noninvasive optical imaging of this physiological response of the retina. For example, OCT imaging systems can be used to image retinal neurons, which exhibit dimensional (size) changes in response to excitation by optical stimuli. These dimensional changes (typically changes in the length of the photoreceptor outer segment (OS) in the retina, i.e., the difference in depth between the inner-outer segment (IS / OS) junction and the cone outer segment tip (COST) of the cone photoreceptor) result in a change in the phase of the light waves returned from the eye that is large enough to be detectable by phase-resolved OCT imaging systems. This phase is highly sensitive to motion within tissue, and many phase-resolved OCT imaging systems can resolve displacements of less than 10 nm, which can be much smaller than the axial resolution of the system or the wavelength of the light used for imaging. Summary of the Invention

[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 second OCT image of a sequence of OCT images of a common portion of a retina of an eye, the common portion including a retinal layer (e.g., an outer segment (OS) of a retinal photoreceptor cell) whose thickness has changed in response to the optical stimulus during acquisition of the sequence of OCT images, acquired by a Fourier-domain optical coherence tomography (OCT) imaging system after stimulation of the common portion with an optical stimulus, to determine at least one of a first index of a first position along an axial (z-axis) direction in the OCT images of a first boundary of the layer and a second index of a second position along an axial direction in the OCT images of a second, different boundary of the layer. The method includes processing the phase components of the first OCT image and the second OCT image to calculate a velocity profile indicative of a distribution of velocities along the axis in the common portion of the retina, and determining at least one of the first index and the second index based on the calculated velocity profile.

[0009] At least one of the first index and the second index can be determined based on the calculated velocity profile by: determining, when the first index is determined, a first position in the velocity profile corresponding to a maximum value of velocity indicated by the velocity profile as the first index; and determining, when the second index is determined, a second position in the velocity profile corresponding to a minimum value of velocity indicated by the velocity profile as the second index. Thus, the velocity of 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 corresponding to a first boundary of the layer at a location on the retina to a minimum velocity at a second position in the velocity profile corresponding to a second boundary of the layer at a location on the retina, and the first index can indicate the first position in the velocity profile of the maximum velocity, and the second index can indicate the second position in the velocity profile of the minimum velocity.

[0010] At least one of the first index and the second index may be determined by processing the calculated velocity profile using a cluster analysis algorithm or a dimension reduction (e.g., component analysis) algorithm. The dimension reduction algorithm may include, for example, 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. When the dimension reduction algorithm includes a PCA algorithm, at least one of the first index and the second index may be determined using principal components determined by applying the PCA algorithm to the calculated velocity profile.

[0011] Both the first position in the velocity profile and the second position in the velocity profile can be determined by: calculating, for each combination of the candidate first position and the candidate second position in the velocity profile, an individual difference value between the individual velocity indicated by the velocity profile of the candidate first position in the combination and the individual velocity indicated by the velocity profile of the candidate second position in the combination, for all combinations of the candidate first position in the velocity profile and the candidate second position in the velocity profile; and identifying, as the first position in the velocity profile and the second position in the velocity profile, the candidate first position and the candidate second position of one combination among the plurality of combinations that has the largest calculated difference value.

[0012] Instead, 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 and the candidate second position in the speed profile, an individual difference value between the average of each speed indicated by the speed profile for a set of adjacent positions including the candidate first position in the combination and the average of each speed indicated by the speed profile for the set of adjacent positions including the candidate second position in the combination, for all combinations of the candidate first position in the speed profile and the candidate second position in 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 one combination among the plurality of combinations for which the calculated difference value is the largest.

[0013] The computer-implemented method of the first exemplary aspect or any of the exemplary implementations thereof described above may further include one of: determining a first indicator of a first position along the axial direction in the OCT image of the first boundary of the layer based on an amplitude component of at least one OCT image in the sequence of OCT images, wherein the second indicator is determined by identifying a second position in the velocity profile corresponding to a minimum value of the velocity indicated by the velocity profile; and determining a second indicator of a second position along the axial direction in the OCT image of the second boundary of the layer based on an amplitude component of at least one OCT image in the sequence of OCT images, wherein the first indicator is determined by identifying a first position in the velocity profile corresponding to a maximum velocity indicated by the velocity profile.

[0014] In some exemplary embodiments, the computer-implemented method may further include calculating a comparison value, which may be either (i) a value of an image quality metric calculated based on at least one of the amplitude components of the first OCT image and the amplitude components of the second OCT image, or (ii) a value of an image similarity metric calculated based on the first OCT image and the second OCT image, providing a measure of similarity between the images. The computer-implemented method may further include comparing the comparison value to 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 at least one of the determined first and second indices 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 at least one of the determined first and second indices is unreliable. The comparison value may be a maximum value of cross-correlation calculated between the first and second OCT images. Other similarity measures that may alternatively be used may be, for example, sum of squared differences, mutual information, normalized mutual information, and Kullback-Leibler distance.

[0015] In some other exemplary embodiments, the computer-implemented method may further include calculating, for each of the plurality of different sets of two or more OCT images in the sequence of OCT images, a respective comparison value, the respective comparison value being either (i) a respective value of an image quality metric calculated based on the amplitude component of at least one OCT image in the set, or (ii) a respective value of an image similarity metric calculated based on the amplitude component of at least two OCT images in the set, the respective value providing a measure of similarity between the images. The method may further include comparing each comparison value to a threshold to determine whether the comparison value is greater than or equal to the threshold. For each of the plurality of different sets of OCT images for which the calculated comparison value is determined to be greater than or equal to the threshold, the phase components of the OCT images in the set are processed to generate a respective velocity profile indicating a distribution. Furthermore, for each of the plurality of different sets of OCT images for which the calculated comparison value is determined not to be greater than the threshold, a respective velocity profile is generated that indicates zero velocity at all positions along the axial direction of the velocity profile. A concatenation of the calculated velocity profiles is then generated such that the concatenation of the velocity profiles indicates how the distribution changes over time, and each portion of the concatenation of the velocity profiles 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 index) is integrated to generate data indicative of the change in optical path length at the position along the axial direction over time. 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 set of the one or more sets of equal consecutive values ​​with one or more estimated values ​​calculated based on adjacent values ​​adjacent to the set of equal consecutive values ​​in the generated data.

[0016] In yet another exemplary embodiment, the computer-implemented method may further include calculating, for each of the plurality of different sets of two or more OCT images in the sequence of OCT images, an individual comparison value, which may be either (i) an individual value of an image quality metric calculated based on the amplitude component of at least one OCT image in the set, or (ii) an individual value of an image similarity metric providing a measure of similarity between the images, calculated based on the amplitude component of at least two OCT images in the set. In these embodiments, the method may further include comparing each comparison value to a threshold to determine whether the comparison value is greater than or equal to the threshold. Some 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 OCT image sets for which the calculated comparison value is determined to be greater than or equal to the threshold, the phase components of the OCT images in the set are processed to generate an individual velocity profile representing the distribution. For each of the plurality of different OCT image sets for which the calculated comparison value is determined not to be greater than or equal to the threshold, an individual velocity profile representing the distribution is generated (estimated) based on the individual velocity profiles calculated for one or more other sets of the plurality of different OCT image sets. A concatenation of the generated velocity profiles is then generated such that the concatenation of the velocity profiles indicates how the distribution changes over time. Portions of the concatenation of velocity profiles having the same axial position are integrated to generate data indicative of the change in optical path length over time at that axial position.

[0017] In further exemplary embodiments, the computer-implemented method may further include calculating, for each of the plurality of different sets of two or more OCT images in the sequence of OCT images, an individual comparison value, which may be either (i) an individual value of an image quality metric calculated based on the amplitude component of at least one OCT image in the set, or (ii) an individual value of an image similarity metric providing a measure of similarity between the images, calculated based on the amplitude component of at least two OCT images in the set. In these exemplary embodiments, each comparison value is compared to a threshold to determine whether the comparison value is greater than or equal to the threshold. Some 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 OCT image sets for which the calculated comparison value is determined to be greater than or equal to the threshold, the phase components of the OCT images in the set are processed to generate an individual velocity profile representing the distribution. For each of the plurality of different OCT image sets for which the calculated comparison value is determined not to be greater than or equal to the threshold, an individual velocity profile representing the distribution is generated based on the individual velocity profiles calculated for one or more other sets of the plurality of different OCT image sets. A concatenation of the generated velocity profiles is then generated such that the concatenation of the velocity profiles indicates how the distribution changes over time, and each portion of the concatenation of the velocity profiles having the same position along the axial direction is integrated to generate data indicative of the change in optical path length over time at a position along the axial direction.

[0018] The individual comparison value calculated for each set may be the individual maximum value of the calculated cross-correlations between at least two OCT images in the set. Other similarity measures that may alternatively be used may be, for example, sum of squared differences, mutual information, and normalized mutual information.

[0019] In an exemplary embodiment in which a linkage is generated, an individual index of the change in optical path length over time at each of the first axial position and at least one of the second axial position indicated by at least one of the first index and the second index may be determined by integrating respective portions of the linkage at the at least one of the first axial position and the second axial position.

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

[0021] Additionally, in any of the foregoing, the computer-implemented method may further include, prior to calculating the velocity profile, processing a phase component of the first OCT image and a phase component of the second OCT image to compensate for bulk motion of the common portion of the retina during acquisition of the sequence of OCT images by the Fourier domain OCT imaging system after stimulation of the common portion by the optical stimulus.

[0022] According to a second exemplary aspect of the present specification, 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 of the exemplary implementations and embodiments thereof described above. The computer program may be stored on a non-transitory computer-readable storage medium (e.g., a computer hard disk or a CD, etc.) or may be carried by a computer-readable signal.

[0023] According to a third exemplary aspect of the present disclosure, there is also provided a data processing device configured to process a phase component of a first OCT image and a phase component of a second OCT image of a sequence of OCT images of a retina of an eye, the OCT image sequence being acquired by a Fourier domain optical coherence tomography (OCT) imaging system after stimulation of the common portion with an optical stimulus, the OCT image sequence including a retinal layer whose thickness has changed in response to the optical stimulus during acquisition of the sequence of OCT images, to determine at least one of a first indicator of a first axial position in the OCT image of a first boundary of the layer and a second indicator of a second axial position in the OCT image of a second boundary of the layer. The data processing device is configured to perform the method of the first exemplary aspect or any of the above-described exemplary implementations and embodiments thereof. The data processing device 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 of the above-described exemplary implementations and embodiments thereof. Alternatively, the data processing device may be implemented in non-programmable hardware, such as an ASIC, an FPGA, or other integrated circuit configured to perform any of these methods.

[0024] According to a fourth exemplary aspect of the present specification, there is also provided a Fourier domain OCT imaging system including the data processing device of the third exemplary aspect described above.

[0025] Exemplary embodiments will now be described in detail, by way of non-limiting example only, with reference to the accompanying drawings in which: like reference numbers appearing in different drawings may indicate identical or functionally similar elements, unless otherwise indicated, and in which: [Brief explanation of the drawings]

[0026] [Figure 1] FIG. 1 is a schematic diagram of a system including a Fourier domain OCT imaging system 30 and a data processing apparatus 100, according to an exemplary embodiment of the present disclosure. [Figure 2]FIG. 2 is a schematic diagram of an exemplary implementation of the data processing device 100 of the exemplary embodiment in programmable signal processing hardware. [Figure 3] FIG. 3 is a flow diagram illustrating a method according to an exemplary embodiment for processing the phase component of an OCT image from a sequence of OCT images of a common portion of the retina acquired after optical stimulation of the common portion to identify layer boundaries. [Figure 4A] FIG. 4A shows a sinusoidal bulk vibration imposed on the oscillatory motion of two model retinal layers. [Figure 4B] Figure 4B shows the effect of removing bulk motion on the change over time in the phase angle of light reflected from two model retinal layers. [Figure 5] Figure 5 shows the ORG contributors for two model retinal layers set apart from each other and not overlapping. [Figure 6] Figure 6 shows the magnitude of the model retinal ORG response when no compensation for bulk motion is provided. [Figure 7A] FIG. 7A shows the phase of the model retinal ORG response when no compensation for bulk motion is provided. [Figure 7B] FIG. 7B shows the phase of the ORG response when the bulk motion is compensated. [Figure 8] FIG. 8 shows the PCA components obtained by applying the PCA algorithm in MATLAB® to the ORG response calculated for the bulk motion compensated case. [Figure 9] Figure 9 shows the absolute values ​​of the principal component coefficients returned by the PCA algorithm. [Figure 10] FIG. 10 shows the ORG contributors of two retinal layers in a variant where the ORG contributors of the two retinal layers partially overlap each other. [Figure 11] FIG. 11 shows the PCA components obtained by applying the PCA algorithm to the ORG response calculated using the ORG contributors of FIG. 10 for the bulk motion compensated case. [Figure 12]FIG. 12 shows the absolute values ​​of the principal component coefficients returned by the modified PCA algorithm. [Figure 13] FIG. 13 is a flow diagram illustrating how the data processing device 100 of an exemplary embodiment can determine whether the (possibly) determined first indicator ZB-1 and / or second indicator ZB-2 are reliable or unreliable and generate an indicator indicating the determined reliability / unreliability. [Figure 14A] FIG. 14A shows the first pair of highly correlated B-scans. [Figure 14B] FIG. 14B shows the calculated cross-correlation of the B-scans of FIG. 14A. [Figure 15A] FIG. 15A shows a second pair of highly uncorrelated B-scans. [Figure 15B] FIG. 15B shows the calculated cross-correlation of the B-scans of FIG. 15A. [Figure 16] FIG. 16 is a flow diagram illustrating a process by which the data processing device 100 of the second exemplary embodiment herein can process pairs of adjacent B-scans in a sequence of B-scans to generate ORG data indicative of the response of the retina to an applied light stimulus. [Figure 17] FIG. 17 shows an example of a concatenation of velocity profiles that may be generated by a data processing device according to process S180 of FIG. [Figure 18] FIG. 18 shows a plot of the change in optical path length ΔOPL over time calculated by the data processing device of the second exemplary embodiment for selected positions within the velocity profile of a concatenation of velocity profiles. [Figure 19] FIG. 19 shows a plot of an exemplary ORG signal generated by the data processing device 100 of the exemplary embodiment herein. [Figure 20] FIG. 20 shows the result of smoothing the ORG signal of FIG. [Figure 21]FIG. 21 is a flow diagram illustrating a process by which a data processing device of a second exemplary embodiment herein can process adjacent OCT images in a sequence of OCT images to generate ORG data indicative of the response of the retina to an applied light stimulus. [Figure 22] FIG. 22 is a flow diagram illustrating a process by which the data processing device of the third exemplary embodiment can process adjacent OCT images in a sequence of OCT images to generate ORG data indicative of the response of the retina to an applied light stimulus. DETAILED DESCRIPTION OF THE INVENTION

[0027] To analyze the OCT response, different layers of the retina are often selected, and the change in optical path length (OPL) between the layers of interest is calculated. When OCT images in the form of B-scans are processed, a flattening algorithm is typically used to distort the B-scan so that retinal layers are at a consistent depth in the distorted B-scan. Retinal layers are then typically found by manual selection or by automatically locating the depth location of the intensity peak in the B-scan that corresponds to the layer. Manual techniques are undesirable because they require operator input, potentially introducing biases that slow the analysis, and require operator knowledge. While automatic detection of intensity peaks typically has a high success rate, it can be problematic when the OCT signal is weak and different layers are not clearly defined, especially in scans with poor axial resolution, scans around the retina, or when the layer of interest has weak intensity (e.g., the inner ganglion layer).

[0028] The OCT data processing methods described herein can solve some or all of these problems and do not rely on identifying OCT signal intensity peaks, which can be difficult when the OCT data has insufficient axial resolution or is noisy. At least some of the OCT data processing methods described herein may enable effective layer segmentation and identification of OCT data from FD-OCT systems with insufficient axial resolution, which may not be sufficient to reveal COST, rod outer segment (ROST), and / or retinal pigment epithelium (RPE) intensity peaks. Furthermore, they may be advantageous in enabling extraction of ORG responses from different regions of the retina, such as ganglion cells, which typically provide weak OCT signal intensity. At least some of the OCT data processing methods described herein also enable finding conventional optical path length responses between photoreceptor layers using conventional point-scanning FD-OCT systems that use sequential registration and do not employ adaptive optics or a tracking system.

[0029] (First Exemplary Embodiment) 1 is a schematic diagram of a data processing device 100 according to a first exemplary embodiment. The data processing device 100 is arranged to process composite OCT data generated by a Fourier-domain OCT (FD-OCT) imaging system. More specifically, the data processing device 100 is arranged to process a phase component (phase information) of a first OCT image and a phase component of a second OCT image of a sequence of OCT images of a common portion of a retina of an eye 20 acquired by a phase-resolved FD-OCT imaging system 30 after stimulation of the common portion of the retina with an optical stimulus 40 generated by an optical stimulus source 45, such as a light-emitting diode (LED). The sequence of OCT images 10 is acquired by the FD-OCT imaging system 30 as the thickness of a layer L of the common portion of the retina changes in response to the optical stimulus 40, and may form part of a wider sequence of OCT images acquired by the FD-OCT imaging system 30 that also includes an OCT image acquired before application of the optical stimulus 40.

[0030] The FD-OCT imaging system 30 may be a swept-source OCT (SS-OCT) system, as in this exemplary embodiment. However, the FD-OCT imaging system 30 need not be provided in this form and may take alternative forms, such as spectral-domain OCT (SD-OCT). More generally, the exemplary embodiment may be provided as any form of phase-resolved FD-OCT imaging system capable of generating complex OCT data, i.e., a Fourier transform of respective spectral interferograms (interference spectra) representing complex A-scan information obtained for each scan position at which OCT measurements are made during a scan. Such complex OCT data encodes phase information from acquired OCT measurements that can be used by the data processing device 100 described herein to identify boundaries of one or both retinal layers of the eye 20.

[0031] The FD-OCT imaging system 30 may include well-known components, including a scanning system, a photodetector, OCT data processing hardware, and a light beam generator (not shown). The scanning system may be configured to perform one-dimensional and / or two-dimensional point scans of a light beam across the retina and collect light scattered by the retina during the point scans. Thus, the scanning system is configured to acquire A-scans at respective scan locations distributed across the surface of the retina by sequentially irradiating the scan locations with the light beam, one scan location at a time, and collecting at least a portion of the light scattered by the retina at each scan location. The scanning system may acquire OCT images in the form of repeated B-scans by performing point scans using a linear scan pattern that follows a set of overlapping scan lines during the scan. In this exemplary embodiment, the scanning system is configured to acquire repeated B-scans by performing point scans; however, in other embodiments, the scanning system may be configured to acquire repeated B-scans by performing line scans using hardware 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 a C-scan by performing a point scan or line scan using techniques well known to those skilled in the art, or by using a full-field setup.

[0032] The first OCT image described above may be in the form of a first OCT B-scan 10-1, as in this exemplary embodiment, and the second OCT image may be in the form of a second OCT B-scan 10-2, as shown in FIG. 1 . However, the forms of the first and second OCT images are not particularly limited; for example, the first and second OCT images may each be an OCT C-scan. The first B-scan 10-1 and second B-scan 10-2 of a common portion of the retina may be adjacent to each other in the sequence of OCT B-scans 10, as in this exemplary embodiment, or may be separated by n intervening B-scans (n≧1, but small enough so that the B-scans have a sufficiently high phase correlation with each other due to the lack of significant relative movement of the eye 20 with respect to the FD-OCT imaging system 30). Furthermore, as described in more detail below, the first B-scan 10-1 may be processed as one of a first set of adjacent B-scans in a sequence 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 a sequence of B-scans 10, and the first and second sets of adjacent OCT images may or may not have one or more OCT images in common.

[0033] Layer L, as in this exemplary embodiment, may be the photoreceptor outer segment (OS) of the retina, which contains rod and cone cells and functions to convert absorbed visible light signals into changes in membrane potential. The photoreceptor OS is particularly well suited for ORG measurements. The photoreceptor cells therein are elongated and behave like optical waveguides, providing relatively strong reflections from each end (i.e., the IS / OS junction and COST). The phase change between these reflections is used to quantify the optical path length change in the OS. However, layer L need not include (or be limited to) the photoreceptor OS, but instead may be another layer or layers of the retina that undergo a measurable change in thickness when stimulated, such as the ganglion cell layer / inner plexiform layer (GCL / IPL). The GCL and IPL still provide detectable reflections, albeit very weak, especially if steps are taken to suppress motion artifacts (see, e.g., C. Pfaffle et al., "Simultaneous functional imaging of neuronal and photoreceptor layers in living human retina," Optic Letters, Vol. 44, No. 23, pp. 5671-5674 (December 1, 2019)).

[0034] As will be 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 generates a first index Z of a first position along the axial direction (z-axis) z in the B-scan 10 of a first boundary (edge) B-1 of a layer L of the retina. B-1 , and a second index Z of a second position along the axial direction z in the B-scan 10 of the second boundary B-2 of the layer L. B-2 The axial direction (or z-axis direction) in a B-scan is whichever direction in the B-scan the elements of the component A-scans are aligned in. If the layer L is a photoreceptor cell OS, as in this exemplary embodiment, 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. The first index Z B-1and the second index Z B-2 Each of the indexes Z can provide a separate row index for the row of the B-scan 10 where the separate boundary is located, as in this exemplary embodiment. B-1 and Z B-2 may instead be scaled versions of these indices (e.g., if the depth of the IS / OS junction and COST from the retinal surface is required).

[0035] As will be explained in more detail below, the data processing device 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 (e.g., as shown in FIG. 1 ) to calculate a tissue velocity profile P, the variation of which with position in the velocity profile comprising values ​​of tissue velocity indicative of the distribution of velocities from point to point along the axial direction z in a common portion of the retina. Thus, the velocity profile P may be, as in this exemplary embodiment, 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., of a calculated change in phase or optical path length) are associated with corresponding values ​​of position. The data processing device 100 calculates a first index Z based on the calculated velocity profile P. B-1 and the second index Z B-2 The method is arranged to determine one or both of:

[0036] The data processing device 100 may be provided in any suitable form, for example as programmable signal processing hardware 200 of the type shown schematically in Figure 2. The programmable signal processing hardware 200 receives a B-scan 10 (or another form of OCT image, such as a C-scan) from the FD-OCT imaging system 30 and (optionally) displays a first index Z B-1 and the second index Z B-2and / or a graphical representation thereof. The communications interface 210 may also output a reliability indicator, which will be described in more detail below. The signal processing hardware 200 further comprises a processor (e.g., a central processing unit, CPU, and / or a graphics processing unit, GPU) 220, a working memory 230 (e.g., random access memory), and an instruction store 240 that stores a computer program 245 containing computer-readable instructions that, when executed by the processor 220, cause the processor 220 to perform various functions of the data processing apparatus 100 described herein. The working memory 230 stores information used by the processor 220 during execution of the computer program 245. The instruction store 240 may comprise a ROM (e.g., in the form of an electrically erasable programmable read-only memory (EEPROM) or flash memory) pre-loaded with computer-readable instructions. Alternatively, the instruction store 240 may comprise RAM or a similar type of memory, and the computer readable instructions of the computer program 245 may be input from a computer program product, such as a non-transitory computer readable storage medium 250 in the form of a CD-ROM, DVD-ROM, etc., or a computer readable signal 260 carrying the computer readable instructions. In either case, the computer program 245, when executed by the processor 220, causes the processor 220 to perform the functions of the data processing apparatus 100 described herein. Thus, the data processing apparatus 100 of the exemplary embodiment comprises a computer processor 220, which, when executed by the processor 220, causes the processor 220 to process 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, and to calculate a first index Z based on the calculated velocity profile P. B-1 and / or the second index Z B-2 and a memory 240 storing computer readable instructions for determining

[0037] However, it should be noted that the data processing apparatus 100 may alternatively be implemented with 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 a combination of such non-programmable hardware with programmable signal processing hardware 200, as described above with reference to FIG. 2.

[0038] The data processing device 100 may be provided as a stand-alone product or as part of a system 1000 comprising an optical stimulus source 45 arranged to provide an optical stimulus 40 and an FD-OCT imaging system 30 arranged to acquire a sequence of OCT images 10 of a common portion of a retina of the eye 20 after stimulation of the common portion by the optical stimulus 40. The data processing device 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 sequence of OCT images 10 acquired by the FD-OCT imaging system 30 using techniques described herein.

[0039] FIG. 3 is a flow diagram illustrating a method according to an exemplary embodiment, in which phase components of B-scans 10-1 and 10-2 from a sequence of B-scans 10 of a common portion of the retina acquired after optical stimulation of the common portion are processed by a data processing device 100 to identify one or both of a first boundary B-1 and a second boundary B-2 of layer L, which may be an IS / OS junction and a COST, respectively, as in this exemplary embodiment.

[0040] 3, the data processing device 100 processes the phase components of the first B-scan 10-1 and the second B-scan 10-2 to calculate a velocity profile P that indicates the distribution of velocities v along the axial direction z within the common portion of the retina. Thus, the data processing device 100 calculates the velocity profile P, which includes velocity v values ​​(or, for example, values ​​of phase change or values ​​of optical path length change ΔOPL) whose variation with position in the velocity profile indicates the distribution of velocities in the portion of the retina (within the common portion) at each point on an axially aligned line. The velocity profile P is thus a two-dimensional data structure in which values ​​of velocity v (or values ​​of a variable indicative thereof, such as ΔOPL or phase change Δφ) are associated with corresponding values ​​at that position, and this data structure can be visualized in the form of a plot of velocity v (or Δφ or ΔOPL) as a function of position z, as shown schematically in FIG.

[0041] Prior to calculating the velocity profile P, the data processing device 100, as in this exemplary embodiment, may process the phase components of the first OCT image 10-1 and the second OCT image 10-2 to compensate for bulk motion of the common portion of the retina during acquisition of the sequence of OCT images 10 by the Fourier domain OCT imaging system 30 after stimulation of the common portion by optical stimulation 40. Compensating for bulk motion has been found to increase the reliability of the velocity-based layer segmentation techniques described herein.

[0042] Figure 4A shows a sinusoidal bulk vibration applied to oscillatory motion of two model retinal layers. In this example, the bulk motion is modeled by adding a phase exp(iφ), where φ represents a fraction of the wavelength (half-wavelength motion is 2π for double-pass reflection). Figure 4A shows how the phase angle of light reflected from a first model retinal layer (Layer 1) and a second model retinal layer (Layer 2) changes over time while the layers are subjected to a common bulk motion as well as their individual vibrations (i.e., the plots labeled "Layer 1 + Bulk" and "Layer 2 + Bulk" in Figure 4A). Figure 4A also shows how the phase angle of light reflected from either model retinal layer changes over time while the layer is subjected to only bulk motion (i.e., the plot labeled "Bulk" in Figure 4A). Figure 4B shows the effect of removing bulk motion on the change in the phase angle of light reflected from the two model retinal layers over time.

[0043] The MATLAB® code used to generate the plots of FIGS. 4A and 4B is as follows:

[0044] t=(0:1000);

[0045] A = 0.95;

[0046] B=0.8;

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

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

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

[0050] Layer1totalMovement=Layer1Movement. * RetinaBulkMovement;

[0051] Layer2totalMovement=Layer2Movement. * RetinaBulkMovement;

[0052] %%%%%%%%%%%%%%%%%

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

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

[0055] subplot(211);

[0056] plot(angle(Layer1totalMovement));hold on;plot(angle(Layer2totalMovement));plot(angle(RetinaBulkMovement));legend(’Layer1+Bulk’,’Layer2+Bulk’,’Bulk’);

[0057] xlim([300 400]);

[0058] subplot(212);

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

[0060] xlim([300 400]);

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

[0062] Referring again to FIG. 3, in process S20, the data processing device 100 calculates a first index Z based on the calculated velocity profile P. B-1 and the second index Z B-2 Therefore, in contrast to the conventional OCT signal amplitude-based layer segmentation commonly used in ORG (whether performed manually by a user or automatically), velocity-based layer segmentation is used in this exemplary embodiment.

[0063] We now describe an example of a process by which the data processing device 100 can perform S10 of Figure 3, which is based on the velocity-based ORG technique described in Kari V. Vienola et al., "Velocity-based optoretinography for clinical applications," Optica 9, 1100-1108 (2022), the entire contents of which are incorporated herein by reference.

[0064] 3, the data processing device 100 can first flatten the first B-scan 10-1 and the second B-scan 10-2 so that the IS / OS and COST reflections are at substantially the same height for each A-scan in each B-scan. Next, the data processing device 100 can align the second B-scan 10-2 with the first B-scan 10-1. Next, 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) can be unwrapped in the temporal dimension to minimize the magnitude of the phase difference between the data sets of the two B-scans 10-1 and 10-2. After unwrapping and processing the phase components of the first OCT image 10-1 and the second OCT image 10-2 to compensate for bulk motion of the retina during imaging, the difference between corresponding phase values ​​of the B-scans 10-1 and 10-2 is calculated for each spatial location specified by the corresponding coordinate pair, and the difference between the acquisition times of the B-scans is then used to calculate the instantaneous velocity of the spatial location. These instantaneous velocities can, as in this exemplary embodiment, be averaged across the lateral dimension (i.e., along the x-axis) to provide a measure of the instantaneous depth dependence of velocity along the z-axis, i.e., a one-dimensional velocity profile P of the type shown schematically in FIG. 1 , which shows the distribution of velocity v along the axis (z-axis) in the common portion of the retina covered by the first and second B-scans 10-1 and 10-2. The B-scan amplitude can also be averaged across the lateral dimension (x-axis) (if desired) to provide a measure of the instantaneous depth dependence of backscatter. It should be noted that averaging of the lateral dimension is not necessary; instead, a two-dimensional velocity profile may be generated showing the distribution of velocity v along both the axial (z-axis) and lateral (x-axis) directions in the common portion of the retina covered by the first B-scan 10-1 and the second B-scan 10-2.

[0065] Once the boundaries of the OS are identified as described below, their respective velocities can be extracted and the difference between them provides the velocity of contraction / extension of the OS at the time B-scans 10-1 and 10-2 were acquired by the FD-OCT imaging system 30, which can be used to generate ORG data indicative of the response of the intersecting retina to an applied stimulus.

[0066] In process S20 of FIG. 3, the data processing device 100 calculates a first index Z B-1 The maximum velocity v in the velocity profile P is max The first position z in the velocity profile P corresponds to max By determining the first index Z of the first boundary B-1 of the layer L along the axial direction z in the B-scan 10, B-1 Here, the maximum value can be a global maximum of the velocity in the velocity profile P, or a maximum value in a predefined portion of the velocity profile P where the OS is expected to be located. The predefined portion can be defined relative to the upper surface of the retina, for example, from the known physiology of the eye. Similarly, the data processing device 100 determines the second index Z in process S20 of FIG. B-2 The minimum velocity value v in the velocity profile P is min A second position z in the velocity profile P corresponds to min By determining the second index Z of the second boundary B-2 of the layer L along the axial direction z in the B-scan 10, B-2 Here, the minimum value can be a global minimum of the velocity in the velocity profile P, or a minimum value within a predefined portion where the OS is expected to be located.

[0067] Therefore, the velocity v of the portion of the velocity profile P corresponding to the photoreceptor OS in the retina is the velocity v of the first position z in the velocity profile P corresponding to the IS-OS junction at a given location on the retina. max Maximum speed at v max from a second location z in the velocity profile P corresponding to COST at that location on the retina. minThe minimum velocity at v min The first index Z B-1 is the maximum speed v max First position z in the velocity profile P max and the second index Z B-2 is the minimum speed v min At the second position z in the velocity profile P min Shows.

[0068] The position(s) of one or both boundaries of the OS or other layers L of the retina that change in response to an applied optical stimulus (e.g., the position of the rod outer segment tip (ROST) or retinal pigment epithelium (RPE)) can be determined using the velocity-based approach for segmentation described herein. For example, the velocity profile P may be determined such that z>z, which corresponds to the change in rod length in response to the applied stimulus. min , there may be further maxima followed by further minima going along the z-axis of FIG. 1. In this exemplary embodiment, the data processing apparatus 100 first calculates all combinations of candidate first and second positions (z i ,z j ) at the first position z of the candidate in the velocity profile P i and the second position z of the candidate in the velocity profile P j For each combination of and, the combination (z i ,z j ) the first position z of the candidate in i Individual velocities v at i and the second position of the candidate z j Individual velocities v at j Difference D ij At the first position z in the velocity profile P, calculate the individual values ​​of max and the second position z min Next, the data processing device 100 determines both the calculated difference D ij The first position z of one of the candidate combinations that maximizes the value of i and the second position z of the candidate j the first position z max and the second position z min Identify as.

[0069] In a variation of this exemplary embodiment, the data processing device 100 first selects all combinations of candidate first and second positions (z i ,z j ) at the first position z of the candidate in the velocity profile P i and the second position z of the candidate in the velocity profile P j For each combination of , the first position z of the candidate in the combination i Five adjacent positions z centered on z (which may otherwise include the first position) i-2 , z i-1 , z i , z i+1 and z i+2 Each velocity v in the set i-2 , v i-1 , v i , v i+1 and v i+2 and the second position z of the candidate in the combination j Five adjacent positions z centered on z (which may otherwise include the second position) j-2 , z j-1 , z j , z j+1 and z j+2 Each velocity v in the set j-2、 v j-1、 v j、 v j+1 and v j+2 The first position z in the velocity profile P is calculated by calculating the individual values ​​of the difference between the average max and the second position z min In this exemplary embodiment, the number of adjacent positions is five, but it will be understood that this is given only as an example and that in other embodiments there may be a different number of adjacent positions (typically two or more). The data processing apparatus 100 then determines, among the calculated difference values, the calculated difference D ij The first position z of one of the candidate combinations that maximizes the value of i and the second position z of the candidate j the first position z max and the second position z min Identify as: Difference Dij Finding the maximum of gives the maximum change in optical path length and therefore indicates where the neuron experiences the greatest change due to stimulation.

[0070] In a further variant of the present exemplary embodiment, the data processing device 100 processes the calculated velocity profile P using a cluster analysis algorithm or a dimension reduction algorithm to obtain the first index Z B-1 and / or the second index Z B-2 These techniques for processing the velocity profile P can be used to identify changes in retinal layers that may be obscured by other layers and / or anatomical features. The dimensionality reduction algorithm may be one of several different types known to those skilled in the art and may be based on machine learning (ML) techniques. Dimensionality reduction algorithms may include, for example, principal component analysis (PCA), independent component analysis (ICA), linear discriminant analysis (LDA), or non-negative matrix factorization (NMF) algorithms.

[0071] When the dimension reduction algorithm is a PCA algorithm, as in this modification, the first index Z B-1 or the second index Z B-2 One or both of the above can be determined using the principal components determined by applying a PCA algorithm to the calculated velocity profile P. To better understand how this can be done, some illustrative examples of how PCA can be used in MATLAB to identify data elements of two model retinal layers that contribute to the ORG response (hereinafter, these data elements are referred to as "ORG contributors"), localize these ORG contributors, identify the layers to which they belong, and thereby define the retinal layers that exhibit the ORG response will be described with reference to Figures 5 to 12.

[0072] In these examples, for simplicity, a single A-scan is considered, and ORG contributors are provided at multiple locations within the A-scan for each of the two layers. In the first of these examples, the ORG contributors for the two retinal layers are spaced apart and do not overlap, as shown in Figure 5, where the x-axis represents the contributor's position within the A-scan and the y-axis represents the contributor's intensity. The MATLAB® code used to generate the plot of Figure 5 is as follows:

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

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

[0075] figure;

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

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

[0078] The phases of model retinal layer 1 and model retinal layer 2 are then assigned to the ORG contributors shown in Figure 5. The magnitude of the (complex) model retinal ORG response when no compensation for bulk motion is provided is shown in Figure 6, and the phase of the ORG response is shown in Figure 7A. The MATLAB® code used to generate Figures 6 and 7A is as follows:

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

[0080] figure;

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

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

[0083] Figure 7B shows the phase of the ORG response when the bulk motion is compensated. The MATLAB code used to generate Figure 7B is as follows:

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

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

[0086] PCA is then used to automatically detect the major ORG contributors and identify their locations within the A-scan. Applying MATLAB's PCA algorithm to the (complex) ORG response calculated for the bulk motion compensated case (the phase of the ORG response is shown in FIG. 7B) results in two PCA components, as shown in FIG. 8. The MATLAB code used to generate FIG. 8 is as follows:

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

[0088] stem(LATENT2);

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

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

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

[0092] In Figure 9, the ORG contributor associated with layer 1 (labeled "retinal element PCA1") appears to the right of the ORG contributor associated with layer 2 (labeled "retinal element PCA2"), which contrasts with Figure 5, where the ORG contributor associated with layer 1 (labeled "retinal element 1") appears to the left of the ORG contributor associated with layer 2 (labeled "retinal element 2"). The ORG contributor values ​​in Figure 9 are also slightly different from those in Figure 5. Nevertheless, the PCA algorithm has identified where the main contributors come from and grouped them.

[0093] In the above example, there is no spatial overlap of the ORG contributors associated with layer 1 and layer 2, as can be seen in Figure 5. However, in a more realistic situation, two ORG contributors may be collected at one location within an A-scan, and we will now describe, with reference to Figures 10-12, the application of PCA to a variation of the above example in which the ORG contributors of two retinal layers partially overlap each other.

[0094] Figure 10 shows the ORG contributors for the two retinal layers in the variation described above. Again, the x-axis represents the contributor's position within the A-scan, and the y-axis represents the contributor's intensity. The MATLAB® code used to generate Figure 10 is as follows:

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

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

[0097] figure();

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

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

[0100] PCA is then used to automatically detect the major ORG contributors and identify their locations within the A-scan. Again, two PCAs are obtained by applying the PCA algorithm in MATLAB to the ORG response calculated for the bulk motion compensated case, as shown in Figure 11. The MATLAB code used to generate Figure 11 is as follows:

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

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

[0103] RetinaORGCompb_t=transpose(RetinaORGCompb);

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

[0105] stem(LATENTb);

[0106] The absolute values ​​of the principal component coefficients returned by the PCA algorithm are plotted in Figure 12. The MATLAB® code used to generate Figure 12 is as follows:

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

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

[0109] As can be seen in Figure 12, even though the layers are shown in reverse order and with slightly different values ​​(as in Figure 9), the PCA algorithm again identified where the major contributors came from and grouped them.

[0110] In addition to, or alternatively to, determining the location of both boundaries of layer L using the velocity-based approach for segmentation described herein, data processing device 100 may provide functionality to determine the location of one of the boundaries of layer L using a conventional OCT signal amplitude-based approach and the remaining boundary using a velocity-based approach, as described below.

[0111] In a variation of an exemplary embodiment having such functionality, the data processing device 100 may include a first position z along the axis z in the B-scan 10. max The first index Z B-1 Determine this first position z maxis the location of the first boundary B-1 of layer L. This determination is based on the amplitude component of at least one B-scan in the sequence of B-scans 10, preferably one of B-scans 10-1 and 10-2 from which velocity profile P is generated as described above. If layer L is a photoreceptor cell OS, as in this case, 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, each of which is associated with a distinct reflectivity peak or maximum in the amplitude of the OCT signal in the B-scan of the sequence of B-scans 10 recorded by the FD-OCT imaging system 30. The data processing device 100 can determine this amplitude peak in a predefined region of the B-scan expected to include the IS / OS junction and COST using any suitable image processing technique known to those skilled in the art, for example, by finding a peak in a moving average of OCT signal amplitude values ​​of one or more A-scans of the B-scan taken along the axial direction of the B-scan. The data processing device 100 then determines the minimum velocity v of the portion in velocity profile P. min The second position z of the velocity profile P corresponds to min A second index Z (of the position of the IS / OS junction along the axial z axis in the B-scan of the OS) is obtained by identifying B-2 The data processing device 100 is further arranged to first determine a position z in the velocity profile P that corresponds to the position of the determined amplitude peak. peak and z peak and all combinations of candidate positions (z peak ,z i ) candidate position z in the velocity profile P i For each combination of and, the combination (z peak ,z i ) in z peak Individual velocities v at peak and candidate position z i Individual velocities v at i Difference D peak_i At the second position z in the velocity profile P, calculate the individual values ​​of min Next, the data processing device 100 identifies the calculated difference D peak_iThe position z of one of the candidate combinations that maximizes the value of i the second position z min Identify as.

[0112] In a further variation of the exemplary embodiment having the functionality described above, the data processing device 100 may be configured to measure a second position z along the axis z in the B-scan 10. min The second index Z B-2 Determine this second position z minは、 The location of the second boundary B-2 of layer L. This determination is based on the amplitude component of at least one B-scan in the sequence of B-scans 10, preferably one of B-scans 10-1 and 10-2 from which velocity profile P is generated as described above. In the present case where layer L is a photoreceptor cell OS, the second boundary B-2 of layer L is a COST, and the first boundary B-1 of layer L is an IS / OS junction, each of which is associated with a distinct reflectivity peak or maximum in the amplitude of the OCT signal in the B-scan of the sequence of B-scans 10 recorded by the FD-OCT imaging system 30. The data processing device 100 can determine this amplitude peak within a predefined region of the B-scan expected to include the IS / OS junction and COST using any suitable image processing technique known to those skilled in the art, for example, by finding a peak in a moving average of OCT signal amplitude values ​​of one or more A-scans of the B-scan taken along the axial direction of the B-scan. The data processing device 100 then determines the maximum velocity v of the portion within velocity profile P. max The first position z of the velocity profile P corresponds to max The first index Z (of the position along the axial z axis in the COST B-scan of the OS) is identified by B-1 The data processing device is further arranged to first determine a position z in the velocity profile P that corresponds to the position of the determined amplitude peak. peak and z peak and all combinations of candidate positions (z peak ,z i ) candidate position z in the velocity profile P i For each combination of peak ,z i) in z peak Individual velocities v at peak and candidate position z i Individual velocities v at i Difference D peak_i At the first position z in the velocity profile P, calculate the individual values ​​of max Next, the data processing device 100 identifies the calculated difference D peak_i The position z of one of the candidate combinations that maximizes the value of i the first position z max Identify as.

[0113] The data processing device 100 of the present exemplary embodiment or any of the variants thereof described above may (possibly) determine a first indicator Z B-1 and / or the second index Z B-2 The data processing device 100 may be further arranged to determine whether the first B-scan 10-1 and the second B-scan 10-2 are reliable or unreliable and to generate an indicator indicative of the determined reliability / unreliability. The reliability of the velocity-based OS segmentation performed by the data processing device 100 is 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 in many amplitude-based approaches to segmentation) and how closely the B-scans are correlated with each other. The maximum calculated cross-correlation between two B-scans is an indicator of how closely the B-scans are correlated, i.e., the first indicator Z determined by the data processing device 100. B-1 and / or the second index Z B-2 The data processing device 100 may communicate an indicator of the determined reliability to a user (e.g., by displaying a message or other type of graphic on the screen informing the user of the determined reliability / unreliability). Additionally or alternatively, the data processing device 100 may (possibly) calculate the first indicator Z determined by the data processing device 100. B-1 and / or the second index Z B-2The indicator associated with the data may be stored. The stored indicator may be used in subsequent data processing operations, as described below.

[0114] More specifically, the data processing apparatus 100 of the present exemplary embodiment or any of its variants described above may be further arranged to carry out the process shown schematically in FIG.

[0115] In process S30 of FIG. 13, 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 can be a value of an image quality metric or attribute calculated based on at least one of the amplitude component of the first B-scan 10-1 and the amplitude component of the second B-scan 10-2. The image quality metric (attribute) can take the form of dynamic range, contrast, signal-to-noise ratio, or a sharpness function value from the B-scans. A measure of success in intensity-based segmentation of the B-scans is another example of an image quality metric that can be used. Alternatively, the comparison value can be a 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 first B-scan 10-1 and the second B-scan 10-2. The image similarity metric can take the form of a maximum value of cross-correlation between the first B-scan 10-1 and the second B-scan 10-2, as in this exemplary embodiment. Other image similarity metrics that may alternatively be used include sum of squared differences, mutual information, normalized mutual information, and Kullback-Leibler distance, among others.

[0116] f * The cross-correlation between two complex functions f(t) and g(t) of a real variable t, denoted by g, is given by f * g=fbar(t)*g(t)(in the formula, * denotes convolution, and fbar(t) is the complex conjugate of f(t).

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

[0118] For comparison, Figure 15A shows two different B-scans. This dissimilarity can be caused by factors such as eye movement, leading to poor image registration. For the exemplary B-scan shown in Figure 15A, the maximum calculated cross-correlation between the B-scans is much lower than for the B-scan shown in Figure 14A, approximately 1 x 10, as shown in Figure 15B. 4 is.

[0119] 13, in process S40, the data processing device 100 compares the comparison value with a threshold value to determine whether the comparison value is greater than or equal to the threshold value. The threshold value is determined by the data processing device 100 (and possibly by the Z value determined by the data processing device 100) while taking into account the maximum cross-correlation value calculated for the source B-scans. B-1 and / or Z B-2 , the segmentation results from the boundary(s) of layer L, as indicated by the value of Z, can be established by comparing the layer segmentation resulting from inspection of the source B-scans by the user. In this way, pairs of B-scans for which the maximum cross-correlation calculated between them is below a certain threshold, can be used to determine the layer segmentation results from the boundary(s) of layer L, as indicated by the value of Z, as indicated by the value of Z, as indicated by the value of Z. B-1 and / or Z B-2 tends to yield unreliable values ​​for Z, but any pair of B-scans for which the maximum cross-correlation calculated between them is greater than or equal to that threshold will (possibly) B-1 and / or Z B-2 It can be determined that the eigenvalues ​​tend to yield reliable values ​​for .

[0120] For example, the B-scans shown in FIG. 14A compare favorably with the results of manual segmentation of the OS (in some cases) performed by inspection of these B-scans. B-1 and / or Z B-2 , and the B-scans shown in FIG. 15A (in some cases) yielded values ​​of Z that were significantly different from the results of manual segmentation of the OS performed by inspection of these B-scans. B-1 and / or Z B-2 In this case, the threshold is set to the maximum value of the cross-correlation shown in FIGS. 14B and 15B (i.e., approximately 3×10 4 and 1 x 10 4 ), e.g., about 1.4×10 as shown by the horizontal black line in FIG. 15B. 4 It may be appropriate to set this to a value of . B-scans resulting in a maximum cross-correlation value below this threshold may be deemed not to allow reliable segmentation to be performed.

[0121] If it is determined in process S40 that the comparison value is equal to or greater than the threshold value ("YES" in S50), the data processing apparatus 100 performs a process S60 in FIG. 13 to calculate the first index Z determined in process S20 in FIG. B-1 and / or the second index Z B-2 Generate a first indicator that the is trustworthy.

[0122] If it is determined in process S40 that the comparison value is not greater than or equal to the threshold value ("No" in S50), the data processing apparatus 100 performs a process S70 in FIG. 13 to calculate the first index Z determined in process S20 in FIG. B-1 and / or the second index Z B-2 Generate a second indicator that the is untrusted.

[0123] (Second Exemplary Embodiment) Above, a single B-scan pair including a first B-scan 10-1 and a second B-scan 10-2 is processed by the data processing device 100 to generate a first index Z B-1 and / or the second index Z B-2However, the described data processing operations may be used to process multiple pairs of B-scans, or C-scans, as in this exemplary embodiment, or more generally, multiple sets of two or more B-scans or C-scans.

[0124] FIG. 16 is a flow diagram illustrating a process by which the data processing device 100 of this exemplary embodiment processes pairs of adjacent B-scans 10-1 and 10-2 in a sequence of B-scans 10 in sequence (e.g., B-scans 10-1 and 10-2 followed by B-scans 10-3 and 10-4 in the sequence of B-scans 10 shown in FIG. 1 ) to generate ORG data indicative of the response of a common portion of the retina to an applied light stimulus.

[0125] Optical retinal imaging focuses on minute changes in retinal layers. The process of B-scan alignment is often used to account for eye movements that may otherwise obscure the ORG response. Alignment involves spatially shifting the two scans relative to each other so that they are at the point of best similarity, as defined by the maximum value of the cross-correlation function defined above. This maximum value of the cross-correlation function provides a good indication of how closely the two B-scans match. If they do not match closely, these scans may produce non-representative results, which may affect the final ORG results. B-scans with poor image quality can also negatively impact ORG results. Therefore, it can be beneficial to avoid including velocity profiles from such B-scans in the analysis, as described here.

[0126] In process S110 of FIG. 16, the data processing device 100 calculates a separate comparison value for each of a plurality of different pairs of consecutive B-scans 10-1, 10-2 in the sequence 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, dynamic range, contrast, signal-to-noise ratio, or a sharpness function value from the B-scans. A measure of success in intensity-based segmentation of the B-scans 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 similarity between the images, the image similarity metric value being calculated based on the amplitude components of the B-scans in the pair. The image similarity metric may take the form of a maximum value of the cross-correlation between the first B-scan 10-1 of the pair and the second B-scan 10-2 of the pair. Other image similarity metrics that may alternatively be used include sum of squared differences, mutual information, normalized mutual information, and Kullback-Leibler distance, among others.

[0127] 16, the data processing device 100 compares each comparison value with a threshold value to determine whether the comparison value is greater than or equal to the threshold value. In this exemplary embodiment, some pairs of B-scans result in comparison values ​​greater than or equal to the threshold value, and some other pairs of B-scans result in comparison values ​​less than the threshold value.

[0128] If the calculated comparison value is determined to be greater than or equal to the threshold value in process S120 of FIG. 16 ("Yes" in S130 of FIG. 16), the data processing device 100 processes the phase components of the pair of B-scans 10-1 and 10-2 in S140 of FIG. 16 to generate an individual velocity profile P showing the distribution, as described above in connection with process S10 of FIG. 3.

[0129] On the other hand, if it is determined in process S120 of FIG. 16 that the calculated comparison value is not greater than or equal to the threshold value (“No” in S130 of FIG. 16), the data processing device 100 generates, in process S150 of FIG. 16, a velocity profile P that shows zero velocity v at all positions along the axial direction z within the velocity profile P, i.e., a null velocity profile.

[0130] Following processes S140 and S150 of Figure 16, the data processing apparatus 100 determines in process S160 whether all pairs of adjacent B-scans in the sequence of B-scans 10 have been processed. If not ("No" in S160), the data processing apparatus 100 selects the next pair of adjacent B-scans in the sequence of B-scans 10 (i.e., the pair of B-scans adjacent to the last pair of B-scans processed up to that point in processing), and the selected pair of B-scans is processed as described above, starting from process S110 of Figure 16. However, if the data processing apparatus 100 determines in process S160 that all pairs of adjacent B-scans in the sequence of B-scans 10 have already been processed ("Yes" in S160 of Figure 16), processing proceeds to S180 of Figure 16, where the data processing apparatus 100 concatenates the generated velocity profiles P such that the concatenation of the velocity profiles shows how the distribution changes over time.

[0131] FIG. 17 shows an example of a concatenation 300 of velocity profiles P generated by data processing apparatus 100 in process S180 of FIG. 16. Each velocity profile P extends along the y-axis of FIG. 17 (labeled "OCT depth layer"), and velocity profiles calculated for adjacent pairs of B-scans in the sequence of B-scans 10 are arranged along the x-axis of FIG. 17 (labeled "time (ms)"). FIG. 17 shows the change in optical path length ΔOPL (providing an indication of velocity v) with position in each velocity profile P derived from B-scans acquired after application of an optical stimulus at time t=200 ms. The change in optical path length ΔOPL varies with position (i.e., along the y-axis of FIG. 17) in each velocity profile from a minimum ΔOPL of approximately −300 (arbitrary units) to a maximum ΔOPL of approximately +250 (arbitrary units). Although there is some variation between the velocity profiles P in the distribution of ΔOPL with position, especially between the velocity profiles for times between 200 ms and 300 ms (immediately after application of the light stimulus at time t = 200 ms), each velocity profile nevertheless has a global maximum near the OCT depth layer 11 associated with IS / OS and a global minimum near the OCT depth layer 20 associated with COST.

[0132] 16, in process S190, the data processing device 100 integrates each portion of the concatenation of velocity profiles P having the same position along the axial direction z in the concatenation of velocity profiles to generate ORG data indicative of the change in optical path length over time at that position along the axial direction z. In process S190 of FIG. 16, the data processing device 100 processes each velocity profile P in the set of concatenated velocity profiles by selecting a portion (segment or data element) of the velocity profile P located at a predetermined position along the axial direction z of the velocity profile P, and can then integrate the selected (co-located) portion by calculating a cumulative sum (or total) of the velocity values ​​in these portions, which indicates how the optical path length of the optical path of the OCT sample beam that ends at a position in the retina corresponding to the predetermined position of the velocity profile P changes over time.

[0133] First index Z B-1 A measure of the change in optical path length over time at a first position along the axial direction z, denoted by Z, can be determined by integrating (i.e., calculating a cumulative sum or total) each portion of the linkage at the first position along the axial direction z. Additionally or alternatively, a second measure Z B-2 An index of the change in optical path length over time at the second location along the axial direction z, denoted by: may be determined by integrating (i.e., calculating a cumulative sum or total) each portion of the coupling at the second location along the axial direction z.

[0134] Figure 18 shows a plot of the change in optical path length ΔOPL over time calculated by the data processing device 100 for selected positions within the velocity profile of a concatenation of velocity profiles. The plot in Figure 18 shows how ΔOPL at a common position varies from one B-scan to the next in a sequence of B-scans 10.

[0135] Figure 19 shows a plot of the cumulative change in optical path length over time, known as the ORG signal (or ORG data). The ORG signal is obtained by calculating the cumulative sum of the optical path length changes for selected locations within a concatenation of velocity profiles. As shown in Figure 19, the ORG signal becomes negative immediately 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 also shows the ORG signal having several plateaus, which, as noted above, are caused by the displacement of the velocity profile calculated using a B-scan that does not result in a sufficiently high maximum cross-correlation value with a velocity profile indicating zero velocity (i.e., a null velocity profile).

[0136] 16 may include one or more sets of equal consecutive values, and the data processing device 100 may be configured 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 ​​calculated based on neighboring values ​​of the set of equal consecutive values ​​in the generated ORG data. The data processing device 100 may perform this smoothing, for example, by using a moving median or moving average, in which the median or average value of a specified number of points on either side of the missing data point(s) is calculated and then assigned to the missing point(s). Alternatively, the data processing device 100 may, for example, detect each set of equal consecutive values ​​in the ORG data and replace the value of each detected set with a corresponding estimated value obtained by interpolating (for example, linearly) between a first value in the ORG data adjacent to the first value of the detected set of equal consecutive values ​​and a second value in the ORG data adjacent to the last value of the detected set of equal consecutive values.

[0137] The results of smoothing the ORG signal using a moving average are shown in Figure 20. Before smoothing is applied, the majority of the data represents zero OPL change, as shown in Figure 19, and therefore generally under-represents the actual values ​​that would be observed. Using the methods described herein, the equivalent values ​​are replaced with estimated values ​​(shown as circles in Figure 20) that more accurately represent the change in OPL.

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

[0139] Furthermore, although the data processing device 100 of this exemplary embodiment is described as processing pairs of B-scans, it may more generally be configured to process sets of three or more B-scans in a sequence of B-scans 10, where the B-scans may be consecutive B-scans in the sequence separated from one another by one or more intervening B-scans that do not form part of the set, or individual B-scans in the sequence. Each set of (e.g., five) B-scans used in a loop of the process of FIG. 16 may be selected by sliding a window selection function along the sequence of B-scans 10, where, in each loop of the process, the window selection function may, for example, select all B-scans in a windowed portion of the sequence of B-scans 10, or every other B-scan in the windowed portion, and the selected sets may have one or more B-scans in common. In such another exemplary embodiment, a modified version of the process described above with reference to FIG. 16, based on the technique described in Kari V. Vienola et al., "Velocity-based optoretinography for clinical applications," Optica 9, pp. 1100-1108 (2022), is as follows:

[0140] In a modification of process S110 of FIG. 16, the data processing device 100 calculates a comparison value for one set of a plurality of different sets of consecutive B-scans in the sequence of B-scans 10. The comparison value can take any of the different forms described above. Then, as in process S120 of FIG. 16, the data processing device 100 compares the comparison value with a threshold to determine whether the comparison value is greater than or equal to the threshold. Some sets of B-scans result in comparison values ​​greater than or equal to the threshold, while some other sets of B-scans result in comparison values ​​less than the threshold. If the calculated comparison value is determined to be greater than or equal to the threshold, the data processing device 100 processes the phase components of at least some of the B-scans in the set to generate respective velocity profiles P indicative of the distribution, according to a modification of process S140 of FIG. 16. Here, the data processing device 100 may first flatten the B-scans of the set so that the IS / OS reflections and the COST reflections are at substantially the same height for each A-scan in each B-scan. Next, the data processing device 100 may align the B-scans with each other. The phase data cube θ(x,z,t) is then calculated by subtracting |θ(x p ,z q ,t r )-θ(x p ,z q ,t r-1 )|, θ(x p ,z q ,t r ) where t r and t r-1 represents a continuous phase B-scan. This step involves scanning a pair of spatial coordinates (x p ,z q) for each scan. The rate of phase change is calculated for each coordinate pair by performing a least-squares linear fit to t, giving Δθ / Δt(x,z) rad / s. From this, the instantaneous velocity at each spatial location can be calculated as Δz / Δt(x,z)=Δθ / Δt(x,z)·λ / 4πn', where λ is the wavelength of the OCT light being 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 a measure of the instantaneous depth dependence of velocity along the z-axis, i.e., a one-dimensional velocity profile P of the type shown schematically in Figure 1, which shows the distribution of velocity v along an axis (z-axis) in the common portion of the retina covered by the repeated B-scans. The B-scan amplitude can also be averaged in the lateral dimension (x-axis) (if desired) to obtain a measure of the instantaneous depth dependence of backscatter. It should be noted that averaging in the lateral dimension is not required; instead, a two-dimensional velocity profile may be generated that shows the distribution of velocities v along both the axial (z-axis) and lateral (x-axis) directions in the common portion of the retina covered by the repeated B-scans. On the other hand, if it is determined that the calculated comparison value is not greater than or equal to the threshold, the data processing device 100 generates a velocity profile P that shows zero velocity v at all positions along the axial direction z within the velocity profile P, i.e., a null velocity profile, in process S150 of Figure 16. The data processing device 100 of the other exemplary embodiments described above will perform further processes as described above with reference to Figure 16, but the references to a "pair" of B-scans in processes S160 and S170 of Figure 16 are replaced with a "set" of B-scans as described above.

[0141] More generally, the data processing apparatus 100 may be arranged to carry out the process described below with reference to FIG.

[0142] 21, the data processing device 100 calculates an individual comparison value for each of a plurality of different sets of two or more OCT images in the sequence of OCT images 10. The individual comparison values ​​may be individual values ​​of an image quality metric calculated based on the amplitude component of at least one OCT image in the set, or individual values ​​of an image similarity metric calculated based on the amplitude components of at least two OCT images in the set, providing a measure of similarity between the images. The comparison values ​​may take any of the exemplary forms described above.

[0143] In process S220 of FIG. 21, the data processing device 100 compares each comparison value with a threshold value to determine whether the comparison value is greater than or equal to the threshold value.

[0144] For each of the multiple different OCT image sets for which the calculated comparison value is determined to be equal to or greater than the threshold, the data processing device 100 processes the phase components of the OCT images in the set to generate an individual velocity profile showing the distribution in process S230 of FIG. 21 .

[0145] Furthermore, for each set of multiple different OCT images for which it is determined that the calculated comparison value is not greater than or equal to the threshold, the data processing device 100 generates, in process S240 of FIG. 21, an individual velocity profile P that shows zero velocity at all positions along the axial direction z of the velocity profile P.

[0146] In process S250 of FIG. 21, the data processing apparatus 100 generates a concatenation of the generated velocity profiles such that the concatenation of the velocity profiles shows how the distribution changes over time.

[0147] Finally, in process S260 of FIG. 21, the data processing device 100 integrates each portion of the concatenation of velocity profiles having the same position along the axial direction to generate data showing the change in optical path length over time at the position along the axial direction.

[0148] (Third Exemplary Embodiment)

[0149] In the second exemplary embodiment described above, the data processing device 100 can reduce degradation in the quality of the generated ORG data by using the comparison values ​​calculated for a set of two or more B-scans in the sequence of B-scans 10 to set a velocity profile exhibiting zero velocity throughout (i.e., a null velocity profile) for use in generating the ORG data and to use in place of any velocity profile derived from a set in which one or more of the component B scans have unacceptable image quality or in which two or more component B scans have insufficient similarity, thereby smoothing the ORG data and optionally suppressing any resulting artifacts in the ORG data. However, the data processing device 100 may alternatively be arranged to reduce degradation in the quality of the generated ORG data by replacing any velocity profile derived from a set in which one or more of the component B scans have unacceptable image quality or in which two or more component B scans have insufficient similarity with an estimated (non-null) velocity profile, as in the present exemplary embodiment, as will now be described in more detail with reference to FIG. 22 .

[0150] Processes S310 to S330, S350, and S360 in FIG. 22 are the same as processes S210 to S230, S250, and S260 in FIG. 21, respectively, which were described in detail above.

[0151] In process S340 of FIG. 22 , the data processing device 100 generates, for each of a plurality of different sets of B-scans for which it has determined that the calculated comparison value is not greater than or equal to the threshold value (i.e., less than the threshold value), an individual estimated velocity profile P indicating the distribution of velocities v along the axial direction (z-axis) in the common portion of the retina covered by the B-scans. Each estimated velocity profile is based on a respective one or more velocity profiles calculated according to process S330 of FIG. 22 and adjacent estimated velocity profiles in the concatenation of the velocity profiles generated in process S350 of FIG. 22 . In other words, for each set of two or more unreliable B-scans for which it has determined that the calculated comparison value is not greater than or equal to the threshold value, the data processing device 100 of the present exemplary embodiment generates an individual velocity profile P based on two or more adjacent reliable sets of B-scans (for which it has determined that the calculated comparison value is greater than or equal to the threshold value) instead of the given set of two or more B-scans. In the present exemplary embodiment, the ORG data generated in process S360 of FIG. 22 does not have a plateau as shown in FIG. 19 , and therefore smoothing of the ORG data may not be necessary.

[0152] Some of the above exemplary embodiments are summarized in the following numbered clauses E1-E18.

[0153] E1. A phase component of a first OCT image 10-1 and a phase component of a second OCT image 10-2 of a sequence of OCT images 10 of a retina of an eye 20, acquired by a Fourier domain optical coherence tomography (OCT) imaging system 30 after stimulation of the intersection with an optical stimulus 40, the OCT image 10 including a retinal layer L whose thickness changed in response to the optical stimulus 40 during acquisition of the sequence of OCT images 10, are processed to obtain a first index Z of a first position along an axial direction z in the OCT image 10 of a first boundary B-1 of the layer L. B-1 and a second index Z of a second position along the axial direction z in the OCT image 10 of the second boundary B-2 of the layer L. B-2 a data processing apparatus 100 arranged to determine at least one of processing the phase components of the first OCT image 10-1 and the second OCT image 10-2 to calculate a velocity profile P indicative of the distribution of velocities v along the axial direction z at the common portion of the retina; Based on the calculated velocity profile, the first indicator Z B-1 and the second index Z B-2 The data processing apparatus 100 is further arranged to determine at least one of:

[0154] E2. The data processing device 100 calculates the first index Z B-1 and the second index Z B-2 at least one of the following based on the calculated velocity profile: First index Z B-1 is determined, the maximum value of the velocity v indicated by the velocity profile P max The first position z in the velocity profile P corresponds to max The first index Z B-1 and Second index Z B-2 is determined, the minimum value of the velocity v indicated by the velocity profile P min A second position z in the velocity profile P corresponds to min the second index Z B-2 The data processing device 100 according to E1 is arranged to determine as and to determine by.

[0155] E3. The data processing device 100 processes the calculated velocity profile P using a cluster analysis algorithm or a dimension reduction algorithm to obtain a first index Z B-1 and / or the second index Z B-2 a data processing device 100 according to E1 or E2, arranged to determine at least one of:

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

[0157] E5. The dimension reduction algorithm includes a PCA algorithm, and the first index Z B-1 and the second index Z B-2 The data processing device 100 according to E4, wherein at least one of the following is determined using principal components determined by applying a PCA algorithm to the calculated velocity profile P:

[0158] E6. The data processing device 100 calculates the velocity profile P at a first position z max and a second position z in the velocity profile P min Both of these determining, for each combination of a candidate first position in the velocity profile P and a candidate second position in the velocity profile P, by calculating an individual value of the difference between the individual velocity indicated by the velocity profile (P) of the candidate first position in the combination and the individual velocity indicated by the velocity profile (P) of the candidate second position in the combination, for all combinations of a candidate first position and a candidate second position in the velocity profile P; The first position and the second position of the candidate for one combination of the plurality of combinations that has the largest calculated difference value are determined as the first position z in the velocity profile. max and a second position z in the velocity profile P min and a data processing device 100 by E2, arranged to identify the data as the

[0159] E7. The data processing device 100 calculates the velocity profile at a first position z max and a second position z in the velocity profile P min Both of these calculating, for each combination of the candidate first position in the velocity profile P and the candidate second position in the velocity profile P, a separate value of the difference between the average of each velocity indicated by the velocity profile P for a set of adjacent positions including the candidate first position in the combination and the average of each velocity indicated by the velocity profile P for a set of adjacent positions including the candidate second position in the combination, for all combinations of the candidate first position and the candidate second position in the velocity profile P; The first position and the second position of the candidate for one combination of the plurality of combinations that has the largest calculated difference value are determined as the first position z in the velocity profile. max and a second position z in the velocity profile P min and a data processing device 100 by E2, arranged to identify the data as the

[0160] E8. The data processing device 100 Based on the amplitude component of at least one OCT image in the sequence of OCT images 10, a first position z along the axial direction z in the OCT image of the first boundary B-1 of the layer L is determined. max The first index Z B-1 and determining a second indicator Z B-2 is a second position z in the velocity profile P that corresponds to the minimum value of the velocity exhibited by the velocity profile P. min and Based on the amplitude component of at least one OCT image in the sequence of OCT images 10, a second position z along the axial direction z in the OCT image of the second boundary of the layer L is determined. min The second index Z B-2 determining a first index Z B-1 is the maximum velocity v indicated by the velocity profile P. max The first position z in the velocity profile P corresponds to max and the data processing device 100 according to any of E1 to E7 is further arranged to perform one of:

[0161] E9. The data processing device 100: a value of an image quality metric calculated based on at least one of the amplitude component of the first OCT image 10-1 and the amplitude component of the second OCT image 10-2; and a value of an image similarity metric calculated based on the amplitude component of the first OCT image 10-1 and the amplitude component of the second OCT image 10-2, the value providing a measure of similarity between the images; Calculate a comparison value that is one of comparing the comparison value to a threshold value to determine whether the comparison value is greater than or equal to the threshold value; generating a first indicator that indicates that at least one of the determined first indicator and second indicator is reliable if the comparison value is determined to be equal to or greater than the threshold value; A data processing device 100 according to any of E1 to E8, further arranged to generate a second indicator indicating that at least one of the determined first indicator and second indicator is unreliable if it is determined that the comparison value is not greater than or equal to the threshold value.

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

[0163] E11. The data processing device is For each set of a plurality of different sets of two or more OCT images 10-1, 10-2 in the sequence of OCT images 10: a discrete value of an image quality metric calculated based on the amplitude component of at least one OCT image in the set; and a distinct value of an image similarity metric that provides a measure of similarity between the images, calculated based on the amplitude components of at least two OCT images in the set; Calculate a separate comparison value, which is one of comparing each comparison value to a threshold to determine whether the comparison value is greater than or equal to the threshold; For each of the plurality of different sets of OCT images for which the calculated comparison value is determined to be equal to or greater than the threshold value, processing the phase components of the OCT images in the set to generate a distinct velocity profile indicating the distribution; generating, for each of the plurality of different OCT image sets for which it is determined that the calculated comparison value is not greater than or equal to the threshold, a separate velocity profile P that exhibits zero velocity at all locations along the axial direction z within the velocity profile P; generating a concatenation 300 of the generated velocity profiles P such that the concatenation 300 of the velocity profiles shows how the distribution changes over time; A data processing device 100 according to any of E1 to E8, further arranged to integrate respective portions of concatenations 300 of velocity profiles P having the same position along the axial direction z to generate data indicative of optical path length change over time at a position along the axial direction z.

[0164] E12. The generated data contains one or more sets of equal consecutive values; the data processing apparatus 100 is further arranged to smooth the generated data by replacing one or more values ​​in one set of the one or more sets of equal consecutive values ​​with one or more estimated values ​​calculated based on adjacent values ​​adjacent to the set of equal consecutive values ​​in the generated data, Data processing device 100 by E11.

[0165] E13. The data processing device is For each set of a plurality of different sets of two or more OCT images 10-1, 10-2 in the sequence of OCT images 10: a discrete value of an image quality metric calculated based on the amplitude component of at least one OCT image in the set; and a distinct value of an image similarity metric that provides a measure of similarity between the images, calculated based on the amplitude components of at least two OCT images in the set; Calculate a separate comparison value, which is one of comparing each comparison value to a threshold to determine whether the comparison value is greater than or equal to the threshold; For each of the plurality of different sets of OCT images for which the calculated comparison value is determined to be equal to or greater than the threshold value, processing the phase components of the OCT images in the set to generate a distinct velocity profile indicating the distribution; For each of the plurality of different OCT image sets for which it is determined that the calculated comparison value is not greater than or equal to the threshold value, generating an individual velocity profile indicating the distribution based on the individual velocity profile calculated for each of one or more other sets of the plurality of different OCT image sets; generating a concatenation 300 of the generated velocity profiles such that the concatenation 300 of the velocity profiles shows how the distribution changes over time; A data processing device 100 according to any of E1 to E8, further arranged to integrate respective portions of concatenated velocity profiles 300 having the same position along the axial direction to generate data indicative of optical path length change over time at position along the axial direction.

[0166] E14. A data processing device 100 according to any of E11 to E13, wherein the data processing device 100 is arranged to calculate, as an individual comparison value for each set, an individual maximum value of the cross-correlation between at least two OCT images in the set.

[0167] E15. The data processing device 100 calculates the first index Z B-1 and the second index Z B-2 a data processing device 100 according to any of E11 to E14, arranged to determine an individual indicator of the change in optical path length over time at each of the first position along the axial direction z and at least one of the second positions along the axial direction z, as indicated by at least one of:

[0168] E16. The data processing device 100 according to any of E1 to E15, wherein the retinal layer includes the outer segment of the photoreceptor cells.

[0169] E17. A data processing device 100 according to any of E1 to E16, wherein the data processing device 100 is further arranged to process the phase component of the first OCT image 10-1 and the phase component of the second OCT image 10-2 prior to calculation of the velocity profile P to compensate for bulk motion of the common portion of the retina during acquisition of the sequence of OCT images 10 by the Fourier domain OCT imaging system 30 after stimulation of the common portion by the optical stimulus 40.

[0170] E18. An optical stimulus source 45 arranged to provide the optical stimulus 40; a Fourier domain optical coherence tomography imaging system 30 arranged to acquire a sequence of OCT images 10 of the common portion of the retina of the eye 20 after stimulation of the common portion with an optical stimulus 40; a data processing device 100 according to any of E1 to E17 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 a sequence of OCT images 10 acquired by a Fourier domain optical coherence tomography imaging system 30; The system 1000 comprises:

[0171] In the foregoing description, exemplary aspects have been described with reference to several exemplary embodiments. Accordingly, the present specification should be considered illustrative rather than restrictive. Similarly, the diagrams shown in the drawings that highlight the functionality and advantages of exemplary embodiments are presented for illustrative purposes only. The architecture of the exemplary embodiments is sufficiently flexible and configurable so that it can be utilized in ways other than those shown in the accompanying figures.

[0172] Some aspects of the examples presented herein, such as the processing methods described with reference to Figures 3, 13, 16, 21, and 22, may be provided as a computer program or software, e.g., one or more programs having instructions or sequences of instructions contained in or stored on an article of manufacture, such as a machine-accessible or machine-readable medium, instruction store, or computer-readable storage device, which, in one example embodiment, may be non-transitory. The program or instructions on the non-transitory machine-accessible or machine-readable medium, instruction store, or computer-readable storage device may be used to program a computer system or other electronic device. Machine or computer-readable media, instruction stores, and storage devices may include, but are not limited to, floppy diskettes, optical disks, and optical-magnetic disks, or other types of media / machine-readable media / instruction stores / storage devices suitable for storing or transmitting electronic instructions. The techniques described herein are not limited to any particular software configuration. They may be applicable in any computing or processing environment. As used herein, the terms "computer-readable," "machine-accessible medium," "machine-readable medium," "instruction store," and "computer-readable storage device" are intended to include any medium capable of storing, encoding, or transmitting instructions or sequences of instructions for execution by a machine, computer, or computer processor that cause the machine / computer / computer processor to perform any one 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 producing 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 and produce a result.

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

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

[0175] Some implementations include software stored on one or more computer-readable media, instruction store(s), or storage device(s) that controls both the system hardware and the system and enables the system or microprocessor to interact with a human user or other mechanism utilizing the results of the exemplary embodiments described herein. Such software can include, but is not limited to, device drivers, operating systems, and user applications. Finally, such computer-readable media or storage device(s) further include software for performing exemplary aspects of the present invention, as described above.

[0176] The programming and / or software of the system includes software modules for performing the procedures described herein. In some exemplary embodiments herein, the modules include software, while in other exemplary embodiments herein, the modules include hardware or a combination of hardware and software.

[0177] While various exemplary embodiments of the present invention have been described above, it should be understood that they are presented by way of example, not limitation. Various changes in form and detail will be apparent to those skilled in the art. Therefore, the present invention should not be limited by any of the above-described exemplary embodiments, but should be defined only in accordance with the following claims and their equivalents.

[0178] While this specification contains details of many specific embodiments, these should not be construed as limitations on the scope of any invention or what may be claimed, but rather as descriptions of features unique to the particular embodiments described herein. Certain features described herein 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 separately in multiple embodiments or in any suitable subcombination. Furthermore, even if features are described above as acting in a particular combination and initially claimed as such, one or more features from a claimed combination may, in some cases, be deleted from the combination, and the claimed combination may also be directed to a subcombination or variations of the subcombination.

[0179] In certain circumstances, multitasking and parallel processing may be advantageous. Furthermore, the separation of various components in the above-described embodiments should not be understood as requiring such separation in all embodiments, and it should be understood that the described program components and systems may generally be integrated together in a single software product or packaged in multiple software products.

[0180] Having now described several exemplary embodiments and implementations, it should be apparent that the foregoing is presented by way of example, not limitation. In particular, while many of the examples presented herein involve particular combinations of device or software elements, those elements may be combined in other ways to achieve the same purpose. Operations, elements, and features discussed only in connection with one embodiment are not intended to be excluded from a similar role in other embodiments or implementations.

Claims

1. 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) of a sequence of OCT images (10) of a retina of an eye (20) acquired by a Fourier domain optical coherence tomography (OCT) imaging system (30) after stimulation of the common portion with an optical stimulus (40), the OCT image including a layer (L) of the retina that has changed in thickness in response to the optical stimulus (40) during acquisition of the sequence of OCT images (10) to determine at least one of a first indicator (ZB-1) of a first location along an axial direction (z) in the OCT image (10) of a first boundary (B-1) of the layer (L) and a second indicator (ZB-2) of a second location along the axial direction (z) in the OCT image (10) of a second boundary (B-2) of the layer (L), the method comprising: processing (S10) the phase components of the first OCT image (10-1) and the second OCT image (10-2) to calculate a velocity profile (P) indicative of a distribution of velocities (v) along the axial direction (z) in the common portion of the retina; and determining (S20) the at least one of the first index (ZB-1) and the second index (ZB-2) based on the calculated velocity profile (P).

2. The at least one of the first index (ZB-1) and the second index (ZB-2) is determined based on the calculated velocity profile (P), When the first index (ZB-1) is determined, 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 the first index (ZB-1); When the second index (ZB-2) is determined, a second position in the velocity profile (P) corresponding to a minimum value of the velocity (v) indicated by the velocity profile (P) is determined as the second index (ZB-2). The computer-implemented method of claim 1 .

3. 3. The computer-implemented method of claim 1, wherein the at least one of the first index (ZB-1) and the second index (ZB-2) is determined by processing the calculated velocity profile (P) using a cluster analysis algorithm or a dimension reduction algorithm.

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

5. 5. The computer-implemented method of claim 4, wherein the dimensionality reduction algorithm includes a principal component analysis algorithm, and wherein the at least one of the first index (ZB-1) and the second index (ZB-2) is determined using principal components determined by applying the principal component analysis algorithm to the calculated velocity profile (P).

6. Both the first position in the velocity profile (P) and the second position in the velocity profile (P) calculating, for each combination of a candidate first position in the velocity profile (P) and a candidate second position in the velocity profile (P), an individual value of the difference between the individual velocity (v) indicated by the velocity profile (P) of the candidate first position in the combination and the individual velocity (v) indicated by the velocity profile (P) of the candidate second position in the combination, for all combinations of a candidate first position and a candidate second position in the velocity profile (P); and identifying a candidate first position and a candidate second position of one combination among a plurality of combinations that produces a maximum calculated difference value among the calculated difference values ​​as the first position in the velocity profile (P) and the second position in the velocity profile (P).

7. Both the first position in the velocity profile (P) and the second position in the velocity profile (P) calculating, for each combination of a candidate first position in the speed profile (P) and a candidate second position in the speed profile (P), an individual value of the difference between the average of each velocity indicated by the speed profile (P) for a set of adjacent positions including the candidate first position in the combination and the average of each velocity indicated by the speed profile (P) for a set of adjacent positions including the candidate second position in the combination, for all combinations of a candidate first position and a candidate second position in the speed profile (P); and identifying the candidate first position and the candidate second position of one combination among a plurality of combinations that produces the largest calculated difference value as the first position in the velocity profile (P) and the second position in the velocity profile (P).

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

9. a value of an image quality metric calculated based on at least one of the amplitude component of the first OCT image (10-1) and the amplitude component of the second OCT image (10-2); and an image similarity metric value, which is a measure of the similarity between the images, calculated based on the amplitude component of the first OCT image (10-1) and the amplitude component of the second OCT image (10-2); calculating a comparison value, the comparison value being one of: comparing the comparison value to a threshold value to determine whether the comparison value is greater than or equal to the threshold value; generating a first indicator that the at least one of the determined first indicator (ZB-1) and second indicator (ZB-2) is reliable if the comparison value is determined to be greater than or equal to the threshold value; and if it is determined that the comparison value is not greater than or equal to the threshold, generating a second indicator that the at least one of the determined first indicator (ZB-1) and second indicator (ZB-2) is unreliable. The computer-implemented method of any one of claims 1 to 8.

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

11. For each set of two or more different sets of OCT images (10-1, 10-2) in the sequence of OCT images (10), a discrete value of an image quality metric calculated based on an amplitude component of at least one of the OCT images in the set; and a distinct value of an image similarity metric that provides a measure of similarity between images, the distinct value being calculated based on amplitude components of at least two of the OCT images in the set; calculating a distinct comparison value, the distinct comparison value being one of: comparing each comparison value to a threshold value to determine whether the comparison value is greater than or equal to the threshold value; For each of the plurality of different OCT image sets for which the calculated comparison value is determined to be greater than or equal to the threshold, processing the phase components of the OCT images in the set to generate a respective velocity profile indicative of the distribution; generating, for each of the plurality of different OCT image sets for which it is determined that the calculated comparison value is not greater than or equal to the threshold, a separate velocity profile P that exhibits zero velocity (v) at all positions along the axial direction (z) within the velocity profile (P); generating a concatenation (300) of the generated velocity profiles such that the concatenation (300) shows how the distribution changes over time; and integrating each portion of the concatenation (300) of the velocity profiles having the same position along the axial direction (z) to generate data indicative of optical path length change over time at the position along the axial direction (z). The computer-implemented method of any one of claims 1 to 8.

12. the generated data includes one or more sets of equal consecutive values; the computer-implemented method further comprising 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 ​​calculated based on neighboring values ​​adjacent to the set of equal consecutive values ​​in the generated data.

12. The computer-implemented method of claim 11.

13. For each set of two or more different sets of OCT images (10-1, 10-2) in the sequence of OCT images (10), a discrete value of an image quality metric calculated based on an amplitude component of at least one of the OCT images in the set; and a distinct value of an image similarity metric that provides a measure of similarity between images, the distinct value being calculated based on amplitude components of at least two of the OCT images in the set; calculating a distinct comparison value, the distinct comparison value being one of: comparing each comparison value to a threshold value to determine whether the comparison value is greater than or equal to the threshold value; For each of the plurality of different OCT image sets for which the calculated comparison value is determined to be greater than or equal to the threshold, processing phase components of the OCT images in the set to generate a respective velocity profile indicative of the distribution; generating, for each of the plurality of different OCT image sets for which it is determined that the calculated comparison value is not greater than or equal to the threshold, an individual velocity profile indicative of the distribution based on the individual velocity profile calculated for each of one or more other sets of the plurality of different OCT image sets; generating a concatenation (300) of the generated velocity profiles P such that the concatenation (300) of velocity profiles shows how the distribution changes over time; and integrating each portion of the concatenation (300) of the velocity profiles having the same position along the axial direction (z) to generate data indicative of optical path length change over time at the position along the axial direction (z). The computer-implemented method of any one of claims 1 to 8.

14. 14. The computer-implemented method of claim 11, wherein the individual comparison value calculated for each set is an individual maximum value of calculated cross-correlations between at least two of the OCT images in the set.

15. 15. The computer-implemented method of claim 11, wherein an individual index of change in optical path length over time at each of the first position along the axial direction (z) and at least one of the second positions along the axial direction (z), indicated by the at least one of the first index (ZB-1) and the second index (ZB-2), is determined by integrating respective portions of the linkage (300) at the 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 layer (L) of the retina comprises outer segments of photoreceptor cells.

17. 17. The computer-implemented method of claim 1, further comprising, prior to calculation of the velocity profile (P), processing the phase component of the first OCT image (10-1) and the phase component of the second OCT image (10-2) to compensate for bulk motion of the common portion of the retina during acquisition of the sequence of OCT images (10) by a Fourier-domain OCT imaging system (30) after stimulation of the common portion by the optical stimulus (40).

18. When executed by a processor (220), the processor (220) A computer program (245) comprising computer readable instructions for carrying out the method of 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) of a sequence of OCT images (10) of a retina of an eye (20) acquired by a Fourier domain optical coherence tomography (OCT) imaging system (30) after stimulation of the common portion with an optical stimulus (40), the OCT image comprising a layer (L) of the retina whose thickness has changed in response to the optical stimulus (40) during acquisition of the sequence of OCT images (10) to determine at least one of a first indicator (ZB-1) of a first position along an axial direction (z) in the OCT image (10) of a first boundary (B-1) of the layer (L) and a second indicator (ZB-2) of a second position along the axial direction (z) in the OCT image (10) of a second boundary (B-2) of the layer (L), processing (S10) the phase components of the first OCT image (10-1) and the second OCT image (10-2) to calculate a velocity profile (P) indicative of a distribution of velocities (v) along the axial direction (z) within the common portion of the retina; A data processing device (100) further arranged to determine (S20) said at least one of said first index (ZB-1) and said second index (ZB-2) based on said calculated velocity profile (P).

Images