ORG retina multilayer optical physiological function signal extraction method based on image intensity analysis

Abstract

The invention provides an ORG retina multilayer optical physiological function signal extraction method based on image intensity analysis, and belongs to the field of image processing. According to the method, the functions of registration alignment, signal conversion, information extraction and the like are integrated, processing and quantitative analysis of obtained retina images are achieved, the situation that the amplitude of different interested retina layers changes along with the time sequence after photostimulation is visually explained, and relevant data support is provided for further analysis.

Description

Technical Field

[0001] The present invention is mainly aimed at the fields of basic and clinical research in ophthalmology. Specifically, it processes the time-series three-dimensional retinal images acquired by the imaging system, quantitatively extracts the dynamics and amplitude of the backscattered bands of interest in the retina, realizes precise and human-computer interactive data processing, and provides data support for further quantitative analysis. Background Art

[0002] In modern society, computers, mobile phones and other electronic devices have become indispensable tools in people's daily lives. Bad eye habits have led to frequent problems such as dry eyes and myopia. In addition, as my country gradually enters an aging society, a series of health problems have emerged. The incidence of blinding diseases such as cataracts, macular degeneration, and diabetic eye disease in the elderly population continues to increase. These diseases usually have no obvious symptoms in the early stages, and it is not until the middle and late stages that irreversible visual impairment or even blindness occurs. Therefore, early detection and intervention of retinal diseases are crucial to preventing visual impairment and protecting eye health.

[0003] At present, there are three commonly used retinal optical imaging technologies in clinical practice: Fundus Camera, Scanning Laser Ophthalmoscopy (SLO) and Optical Coherence Tomography (OCT). Among them, OCT technology is widely used in clinical examinations of various ophthalmic diseases with higher resolution. The human retina is relatively complex, including: Nerve Fiber Layer NFL (Nerve Fiber Layer), Inner Plexiform Layer IPL (Inner Plexiform Layer), Outer Plexiform Layer OPL (Outer Plexiform Layer), External Limiting Membrane ELM (External Limiting Membrane), Inner Segment / Outer Segment Junction ISOS (Inner Segment / Outer Segment Junction), Retinal Outer Segment Tips ROST (Retinal Outer Segment Tips), Retinal Pigment Epithelium-Bruch's Membrane RPE-BrM (Retinal Pigment Epithelium-Bruch's Membrane), Choriocapillaris CC (Choriocapillaris), etc.

[0004] Due to the complexity of the human retina and the high cost of research, mouse models with a high degree of similarity to the human retina have become an indispensable research object. Among them, the sample profile along the depth direction obtained by the OCT system is A-scan, and the two-dimensional cross-section composed of a series of closely adjacent A-scans is B-scan. Combining continuous B-scans can generate 3D OCT images of the retina. Time-series three-dimensional mouse retinal images can record the dynamic changes of the retina in real time, reflecting the changing states of different layers of the retina at different time points.

[0005] However, noise and artifact problems are difficult to avoid during image acquisition. Especially in biological experiments, due to the slight movement of the eyeball, rigid or nonlinear differences may occur between frames. In order to achieve more accurate qualitative comparison, registration becomes a challenging task in retinal OCT imaging. Summary of the invention

[0006] In view of the problems existing in the prior art, the present invention provides a method for extracting signals of multi-layer photophysiological function signals of ORG retina based on image intensity analysis. The method integrates functions such as registration alignment, signal conversion, and information extraction, realizes the processing and quantitative analysis of the acquired retinal images, and visualizes the changes in amplitude of different retinal layers of interest over time after light stimulation, providing relevant data support for further analysis.

[0007] The technical solution of the present invention:

[0008] A method for extracting multi-layer photophysiological function signals of ORG retina based on image intensity analysis, the steps are as follows:

[0009] Step 1: Preprocess the image

[0010] The OCT images were flipped vertically to ensure that the image orientation was consistent with the subsequent data processing standards for more accurate signal extraction.

[0011] Step 2: Initialize the parameter values

[0012] The parameter values ​​include the number of B-scans in each group M b , the number of frames of median filtering μ f , the average number of pixels in the horizontal direction μ l , the actual length C of one pixel in micrometers f , axial resolution R.

[0013] Step 3: Preliminary image alignment

[0014] (3-1) Averaging multiple adjacent A-scans: μ lThe preprocessed OCT image is divided into blocks in the transverse direction, and the average value of all A-scan positions in each block is calculated.

[0015] (3-2) For the current small block, the average value of the A-scan position of the previous small block is used as the reference signal, and the similarity of each A-scan is calculated frame by frame in the frame number or time direction. The cross-correlation function for calculating similarity is:

[0016]

[0017] Where n is the depth position or pixel index of the A-scan to be processed in the current small block, k is the offset of the A-scan to be processed in the current small block to the reference signal, x(n) is the signal strength value of the A-scan of the current small block at the depth position n, y(nk) is the reference signal strength value of the corresponding depth position of the previous small block, R xy (k) is the cross-correlation function.

[0018] The relative displacement is determined by finding the maximum value of similarity, and each A-scan in the small block to be processed is aligned. The displacement function for alignment is:

[0019]

[0020] Where Δ is the difference between the reference A-scan and the current A-scan, n is the depth position or pixel index of the A-scan to be processed, and x ′ (n) is the signal strength value at depth position n in the A-scan after displacement.

[0021] By translating the A-scan up and down and padding the bottom or top with zeros, the A-scan in each small block is aligned.

[0022] (3-3) To suppress noise and enhance image quality, M b The B-scan is grouped in the frame direction with the step size of μ. f Perform median filtering for the step size:

[0023] If μ f is an odd number, take the signal strength value of the middle B-scan:

[0024]

[0025] If μ f If it is an even number, take the average signal strength of the two middle B-scans:

[0026]

[0027] Among them, μ f is the number of frames of median filtering, is the signal intensity value of the B-scan at the middle position of the step, and MedD(t) is the average signal intensity value of the B-scan after averaging.

[0028] After obtaining a smoother B-scan signal through median filtering, the B-scan data in each block are horizontally spliced, and finally the three-dimensional OCT data is converted into a two-dimensional image.

[0029] Step 4: Fine alignment of images (4-1) Rough alignment between blocks: Manually select a layer from all blocks as the reference layer (generally the ELM layer is selected as the reference layer), and perform cyclic shift on the B-scans of different blocks according to the selected reference layer. The formula for cyclic shift is:

[0030] A shifted (i) = A((i-ε-1)modN+1)

[0031] Where i is the position of the manually selected reference layer in the B-scan in the block, ε is the offset of the position of the manually selected reference layer in the B-scan in the current small block to the position of the manually selected reference layer in the previous small block, N is the number of small blocks, mod represents the modulus operation, A shifted (i) is the i-th element of the array after translation.

[0032] (4-2) Displacement calculation: By manually selecting the upper and lower limits of the reference layer as the B-scan alignment standard area, find the position of the maximum value in the alignment standard area in each column, compare each column of B-scan with the maximum value of the alignment standard area of ​​the first column of B-scan, and calculate the displacement:

[0033] S(j)=L(1)-L(j)

[0034] Where L(1) is the row where the maximum value of the first column B-scan is located, L(j) is the row where the maximum value of the j-th column is located, and S(j) is the number of pixels that the j-th column B-scan needs to be shifted compared to the first column B-scan.

[0035] (4-3) Fine alignment: Finely align the image by performing linear displacement up and down using the displacement calculated in step (4-2) so that the peaks are aligned. The displacement function for alignment is:

[0036] A ′ (i,j)=A(iS(j),j)

[0037] Among them, A ′(i, j) is the intensity value of the new B-scan signal after translation. If S(j) is a positive value, the B-scan signal is translated upward. If S(j) is a negative value, the B-scan signal is translated downward.

[0038] Step 5: Converting from raw OCT images to processed signals

[0039] (5-1) Sorting and noise reduction: The finely aligned data is sorted in the time dimension, and the average value of the middle c blocks is calculated to reduce the impact of noise and improve data quality.

[0040] (5-2) PCHIP difference: The sorted and denoised OCT signals are converted from original coordinates to physical coordinates, and the PCHIP interpolation method is applied to make the signal curve smoother, thereby improving the signal quality.

[0041] (5-3) Gaussian fitting: The Gaussian fitting formula is defined as:

[0042]

[0043] Where x is the pixel position of A-scan, x c is the center position of the Gaussian function, and σ is the standard deviation.

[0044] Furthermore, when When , the function value is half of the maximum value. Then:

[0045]

[0046] Taking the natural logarithm and solving the equation gives:

[0047]

[0048] Here, FWHM is the half-maximum width.

[0049] (5-4) In order to find the peak values ​​of each layer more accurately, Lucy-Richardson deconvolution is used to restore the original signal. In the deconvolution model:

[0050]

[0051] Among them, f r (x) is the estimated result of the rth iteration of the A-scan signal strength at position x, g(x) is the measured A-scan signal, h(x) is the point spread function PSF, which is the Gaussian function f(x), and * represents the convolution operation.

[0052] (5-5) Calculate the deconvolution residual: Convolve the restored original signal with the PSF, reintroduce the blurred or distorted part, and calculate the difference between the convolution result and the interpolated data in step (5-2).

[0053] Step 6: Layered display and analysis

[0054] (6-1) Moving average: For each data point y i Perform a moving average:

[0055]

[0056] Among them, x j is the signal strength value of the current A-scan at depth position j after deconvolution, y i is the signal strength value at depth position i after smoothing, and m represents half the size of the moving window.

[0057] (6-2) Layered processing of column data: Manually select the boundaries on both sides of the peaks in the NFL, IPL, OPL, ELM, ISOS, ROST, and RPE-BrM layers, and manually select the boundaries on both sides of the valleys in the CC layer. Initialize the movement amount and starting point, calculate the maximum value of the peak or valley, calculate the movement amount and update the boundary, and use the coefficient of the previous Gaussian fitting result as the new starting point for fitting. Among them, the NFL, IPL, OPL, ELM, ISOS, and CC layers use single-peak Gaussian fitting, and the ROST and RPE-BrM layers use double-peak Gaussian fitting;

[0058] Single-peak Gaussian fitting: Set the initial value and upper and lower bounds of the peak height, center position, and standard deviation in the Gaussian model respectively, and use the nonlinear least squares method to perform a single-peak Gaussian fitting on the data to ensure that the fitting result is closer to the actual signal characteristics. The selected single-peak Gaussian fitting model formula is:

[0059]

[0060] Where a is the height of the peak; b is the center position of the peak in this layer; c is the standard deviation, which determines the width of the peak.

[0061] Bimodal Gaussian fitting: Set the initial values ​​and upper and lower bounds of the peak height, center position, and standard deviation in the Gaussian model respectively, and use the nonlinear least squares method to perform bimodal Gaussian fitting on the data to obtain the characteristic parameters of the two main peaks of the signal to ensure that the fitting result is closer to the actual signal characteristics. The corresponding formula is:

[0062]

[0063] where a 1 is the height of the first peak, a 2is the height of the second peak; b 1 is the center position of the first peak, b 2 is the center position of the second peak; c 1 is the width of the first peak, c 2 is the width of the second peak.

[0064] (6-3) Obtaining fitting data: Based on the fitting results, the peak and valley intensity, position and standard deviation of different layers are obtained over time after light stimulation, and the multi-layer photophysiological function signal of the ORG retina is extracted, laying the foundation for further data analysis.

[0065] The beneficial effects of the present invention are as follows: the present invention pre-processes the retinal photophysiological function signals, takes Matlab as the core language, effectively suppresses the influence of noise through median filtering, Gaussian fitting, deconvolution and other algorithms, and performs quantitative analysis on different layers of the retina, intuitively sees the changes in the backscattered signals, provides data support for further data analysis, and provides an important reference basis for the early detection, early diagnosis and treatment of ophthalmic diseases. BRIEF DESCRIPTION OF THE DRAWINGS

[0066] Figure 1 It is a schematic diagram of the process of the present invention.

[0067] Figure 2 Figure 1 shows the functional response of the retina to visible light stimulation, where (a) is a schematic diagram of the spatiotemporal processing of OCT B-scan data. The symbol t represents time, k = 1, ..., n (usually n = 15-20) identifies the B-scan subregions used for analysis; (b) is a separate schematic diagram of OCTB-scan data, and (c) is an overlay of the average A-scan curves at different time points, with the color bar indicating time); the black dashed rectangle indicates the section used for analysis of rod length and scattering, and the cyan inset provides a cartoon model of a single rod drawn approximately to scale.

[0068] Figure 3(a) shows the image after preliminary alignment.

[0069] Figure 3(b) shows the image after fine alignment.

[0070] Figure 4(a) is a schematic diagram before the averaging operation.

[0071] Figure 4(b) is a schematic diagram after the averaging operation.

[0072] Figure 5(a) is a schematic diagram of the rough alignment operation between blocks.

[0073] Figure 5(b) is a schematic diagram of the precise alignment operation between blocks.

[0074] Figure 6(a) is a schematic diagram of ISOS hierarchical point selection.

[0075] Figure 6(b) is a schematic diagram of the ISOS stratification results.

[0076] Figure 6(c) is a schematic diagram of the layered point selection of ROST-RPE.

[0077] Figure 6(d) is a schematic diagram of the ROST stratification results.

[0078] Figure 6(e) is a schematic diagram of the RPE stratification results. DETAILED DESCRIPTION

[0079] The technical solution of the present invention will be further described below in conjunction with specific examples and drawings.

[0080] Example: Extraction of multi-layer photophysiological function signals from the retina of black mice

[0081] Mice are an important model for studying human retinal diseases and physiological mechanisms. The two have sufficient similarities in anatomical structure, function and gene expression. Through in-depth research on the mouse retina, scientists can obtain key biological information and potential therapeutic targets about retinal diseases, provide a theoretical basis for the study of human retinal diseases, and provide a reference for the development of clinical treatment methods.

[0082] C57BL / 6 mice were used as experimental models, and the method of the present invention was combined with Appdesigner of Matlab language to build a GUI operation interface. The OCT system used a scanning rate of 120 Hz for data acquisition. During the experiment, light stimulation was at the 30th second, and the system saved an image every 10 frames, and finally obtained 4800 images, with a total scanning time of 400 seconds.

[0083] Step 1: Preprocess the image

[0084] The OCT images were flipped vertically to ensure that the image orientation was consistent with the subsequent data processing standards for more accurate signal extraction.

[0085] Step 2: Initialize the parameter values

[0086] Set the number of B-scans per group M b =4, the number of frames of median filtering μ f 3, the horizontal average number of pixels μ l = 20, the actual length C per pixel expressed in microns f The axial resolution R of the system is 0.992 and 2.9.

[0087] Step 3: Preliminary image alignment

[0088] (3-1) Averaging multiple adjacent A-scans: The pre-processed OCT image of size 446*668*4800 is averaged in the lateral direction with a step size μ l =20 for block division, and a total of 668 / 20=33 small blocks. For each small block, the average signal strength of 20 adjacent A-scans in each block is calculated to obtain the A-scan average signal of the small block.

[0089] (3-2) For the current small block, the average value of the A-scan position of the previous small block is used as the reference signal, and the similarity of each A-scan is calculated frame by frame in the frame number or time direction. The cross-correlation function for calculating similarity is:

[0090]

[0091] Where n is the depth position or pixel index of the A-scan to be processed in the current small block, k is the offset of the A-scan to be processed in the current small block to the reference signal, x(n) is the signal strength value of the A-scan of the current small block at the depth position n, y(nk) is the reference signal strength value of the corresponding depth position of the previous small block, R xy (k) is the cross-correlation function.

[0092] The relative displacement is determined by finding the maximum value of similarity, and each A-scan in the small block to be processed is aligned. The displacement function for alignment is:

[0093]

[0094] Where Δ is the difference between the reference A-scan and the current A-scan, n is the depth position or pixel index of the A-scan to be processed, and x ′ (n) is the signal strength value at depth position n in the A-scan after displacement.

[0095] By translating the A-scan up and down and padding the bottom or top with zeros, the A-scan in each small block is aligned.

[0096] (3-3) To suppress noise and enhance image quality, M b =4 is the step size, 4800 B-scans are divided into 1200 blocks in the frame direction, each block contains four consecutive B-scans. f =3 is the step size for median filtering:

[0097] μ f =3 is an odd number, and the signal strength values ​​of the 4 B-scans of the middle block at each depth position (446 depth points) are taken:

[0098] MedD(t)=D 2 (t)

[0099] If μ f If the value is an even number, take the average signal intensity of the B-scans of the two middle groups:

[0100]

[0101] Among them, μ f is the number of frames of median filtering, is the signal intensity value of the B-scan at the middle position of the step, and MedD(t) is the average signal intensity value of the B-scan after averaging.

[0102] After obtaining a smoother B-scan signal through median filtering, the OCT image with a size of 446*668*4800 is converted into an OCT image with a size of 446*33*400. Each piece of B-scan data is horizontally spliced, and finally the three-dimensional OCT data is converted into a two-dimensional image with a size of 446*13200.

[0103] Step 4: Fine-tune image alignment

[0104] (4-1) Rough alignment between blocks: The ELM layer in all blocks is manually selected as the reference layer, and the B-scans of different blocks are cyclically shifted according to the selected ELM layer, as shown in Figure 5(a). The formula for cyclic shift is:

[0105] A shifted (i) = A((i-ε-1)modN+1)

[0106] Where i is the position of the manually selected reference layer in the B-scan of the block, ε is the offset of the position of the ELM layer in the B-scan of the current small block to the position of the ELM layer in the previous small block, N is the number of small blocks 33, mod represents the modulo operation, A shifted (i) is the i-th element of the array after translation.

[0107] (4-2) Calculation of displacement: By manually selecting the upper and lower limits of the ELM layer as the B-scan alignment standard area, find the position of the maximum value in the alignment standard area in each column, compare each column of B-scan with the maximum value of the alignment standard area of ​​the first column of B-scan, and calculate the displacement:

[0108] S(j)=L(1)-L(j)

[0109] Where L(1) is the row where the maximum value of the first column B-scan is located, L(j) is the row where the maximum value of the j-th column is located, and S(j) is the number of pixels that the j-th column B-scan needs to be shifted compared to the first column B-scan.

[0110] (4-3) Fine alignment: Finely align the images by performing linear displacement up and down using the displacement calculated in step (4-2) so that the peaks are aligned as shown in Figure 5(b). The displacement function for alignment is:

[0111] A ′ (i,j)=A(iS(j),j)

[0112] Among them, A ′ (i, j) is the intensity value of the new B-scan signal after translation. If S(j) is a positive value, the B-scan signal is translated upward. If S(j) is a negative value, the B-scan signal is translated downward.

[0113] Step 5: Converting from raw OCT images to processed signals

[0114] (5-1) Sorting and noise reduction: The finely aligned data is sorted in the time dimension, and the average of the middle five blocks is calculated to reduce the impact of noise and improve data quality.

[0115] (5-2) PCHIP difference: The sorted and denoised OCT signals are converted from original coordinates to physical coordinates, and the PCHIP interpolation method is applied to make the signal curve smoother, thereby improving the signal quality.

[0116] (5-3) Gaussian fitting: The Gaussian fitting formula is defined as:

[0117]

[0118] Where x is the pixel position of A-scan, x c is the center position of the Gaussian function, and σ is the standard deviation.

[0119] Furthermore, when When , the function value is half of the maximum value. Then:

[0120]

[0121] Taking the natural logarithm and solving the equation gives:

[0122]

[0123] The axial resolution of the system is R=2.9 as the half-maximum width FWHM.

[0124] (5-4) In order to find the peak values ​​of each layer more accurately, Lucy-Richardson deconvolution is used to restore the original signal. In the deconvolution model:

[0125]

[0126] Among them, f r (x) is the estimated result of the rth iteration of the A-scan signal strength at position x, g(x) is the measured A-scan signal, h(x) is the point spread function PSF, which is the Gaussian function f(x), and * represents the convolution operation, which is iterated 10 times in total.

[0127] (5-5) Calculate the deconvolution residual: Convolve the restored original signal with the PSF, reintroduce the blurred or distorted part, and calculate the difference between the convolution result and the interpolated data in step (5-2).

[0128] Step 6: Layered display and analysis

[0129] (6-1) Moving average: For each data point y i Perform a moving average:

[0130]

[0131] Among them, x j is the signal strength value of the current A-scan at depth position j after deconvolution, y i is the signal strength value at depth position i after smoothing, and m represents half the size of the moving window, which is 25.

[0132] (6-2) ISOS layer processing column data: Manually select the boundaries on both sides of the peak in the ISOS layer, as shown in Figure 6(a). Initialize the movement amount and starting point, calculate the peak maximum, calculate the movement amount and update the boundary, and use the coefficients of the previous Gaussian fitting result as the new starting point for fitting.

[0133] Single-peak Gaussian fitting: Set the initial value and upper and lower bounds of the peak height, center position, and standard deviation in the Gaussian model respectively, and use the nonlinear least squares method to perform a single-peak Gaussian fitting on the data to ensure that the fitting result is closer to the actual signal characteristics. The selected single-peak Gaussian fitting model formula is:

[0134]

[0135] Where a is the height of the peak; b is the center position of the peak in this layer; c is the standard deviation, which determines the width of the peak.

[0136] The fitting results are shown in Figure 6(b).

[0137] Processing column data of ROST and RPE layers: Manually select the left boundary of the ROST layer peak and the right boundary of the RPE layer peak, and perform a double-peak Gaussian fitting, as shown in Figure 6(c).

[0138] Bimodal Gaussian fitting: Set the initial values ​​and upper and lower bounds of the peak height, center position, and standard deviation in the Gaussian model respectively, and use the nonlinear least squares method to perform bimodal Gaussian fitting on the data to obtain the characteristic parameters of the two main peaks of the signal to ensure that the fitting result is closer to the actual signal characteristics. The corresponding formula is:

[0139]

[0140] where a 1 is the height of the first peak, a 2 is the height of the second peak; b 1 is the center position of the first peak, b 2 is the center position of the second peak; c 1 is the width of the first peak, c 2 is the width of the second peak.

[0141] The fitting results are shown in Figure 6(d) and Figure 6(e).

[0142] (6-3) Obtaining fitting data: Based on the fitting results, the peak intensity, position and width of the ISOS layer, ROST layer and RPE layer are respectively obtained as they change with time after light stimulation, and the multi-layer photophysiological function signal of the ORG retina is extracted, laying the foundation for further data analysis.

Claims

1. A method for extracting multi-layer photophysiological function signals of ORG retina based on image intensity analysis, characterized in that: Here are the steps: Step 1: Preprocess the image Flip the OCT image vertically to ensure that the image orientation is consistent with the subsequent data processing standards for more accurate signal extraction; Step 2: Initialize the parameter values The parameter values ​​include the number of B-scans in each group M b , median filter frame number μ f , the average number of pixels in the horizontal direction μ l , the actual length C of one pixel in micrometers f , axial resolution R; Step 3: Preliminary image alignment (3-1) Averaging multiple adjacent A-scans: μ l The preprocessed OCT image is divided into blocks in the transverse direction, and the average value of all A-scan positions in each block is calculated; (3-2) For the current small block, the average value of the A-scan position of the previous small block is used as the reference signal, and the similarity of each A-scan is calculated frame by frame in the frame number or time direction; The cross-correlation function for calculating similarity is: Where n is the depth position or pixel index of the A-scan to be processed in the current small block, k is the offset of the A-scan to be processed in the current small block to the reference signal, and x(n) is the signal strength value of the A-scan of the current small block at the depth position n. y(nk) is the reference signal strength value of the previous small block corresponding to the depth position, R xy (k) is the cross-correlation function; The relative displacement is determined by finding the maximum value of similarity, and each A-scan in the small block to be processed is aligned; the displacement function for alignment is: Where Δ is the difference between the reference A-scan and the current A-scan, n is the depth position or pixel index of the A-scan to be processed, and x ′ (n) is the signal strength value at depth position n in the A-scan after displacement; By translating the A-scan up and down, and padding the bottom or top with zeros, the A-scan in each small block is aligned. (3-3) To suppress noise and enhance image quality, M b The B-scan is grouped in the frame direction with the step size; the grouped B-scan is divided into μ f Perform median filtering for the step size: If μ f is an odd number, take the signal strength value of the middle B-scan: If μ f If it is an even number, take the average signal strength of the two middle B-scans: Among them, μ f is the number of frames of median filtering, is the signal intensity value of the B-scan at the middle position of the step, and MedD(t) is the average signal intensity value of the B-scan after averaging; After obtaining a smoother B-scan signal through median filtering, the B-scan data in each block are spliced ​​horizontally, and finally the three-dimensional OCT data is converted into a two-dimensional image; Step 4: Fine-tune image alignment (4-1) Rough alignment between blocks: A layer among all blocks is manually selected as the reference layer, and the B-scans of different blocks are cyclically shifted according to the selected reference layer. The formula for cyclic shift is: A shifted (i)=A((i-ε-1)modN+1) Where i is the position of the manually selected reference layer in the B-scan in the block, ε is the offset of the position of the manually selected reference layer in the B-scan in the current small block to the position of the manually selected reference layer in the previous small block, N is the number of small blocks, mod represents the modulus operation, A shifted (i) is the i-th element of the array after translation; (4-2) Displacement calculation: By manually selecting the upper and lower limits of the reference layer as the B-scan alignment standard area, find the position of the maximum value in the alignment standard area in each column, compare each column of B-scan with the maximum value of the alignment standard area of ​​the first column of B-scan, and calculate the displacement: S(j)=L(1)-L(j) Where L(1) is the row with the maximum value of the first column B-scan, L(j) is the row with the maximum value of the jth column, and S(j) is the number of pixels that the jth column B-scan needs to be shifted compared to the first column B-scan; (4-3) Fine alignment: Finely align the image by performing up and down linear displacements based on the displacement calculated in step (4-2) so that the peaks are aligned; the displacement function for alignment is: A ′ (i,j)=A(i-S(j),j) Among them, A ′ (i, j) is the intensity value of the new B-scan signal after translation. If S(j) is a positive value, the B-scan signal is translated upward, and if S(j) is a negative value, the B-scan signal is translated downward; Step 5: Converting from raw OCT images to processed signals (5-1) Sorting and noise reduction: The finely aligned data is sorted in the time dimension, and the average value of the middle c blocks is calculated to reduce the impact of noise and improve data quality; (5-2) PCHIP difference: The sorted and denoised OCT signals are converted from original coordinates to physical coordinates, and the PCHIP interpolation method is applied to make the signal curve smoother, thereby improving the signal quality; (5-3) Gaussian fitting: The Gaussian fitting formula is defined as: Where x is the pixel position of A-scan, x c is the center position of the Gaussian function, σ is the standard deviation; Furthermore, when When , the function value is half of the maximum value; then: Taking the natural logarithm and solving the equation gives: Among them, FWHM is the half-maximum width; (5-4) In order to find the peak values ​​of each layer more accurately, Lucy-Richardson deconvolution is used to restore the original signal. In the deconvolution model: Among them, f r (x) is the estimated result of the rth iteration of the A-scan signal strength at position x, g(x) is the measured A-scan signal, h(x) is the point spread function PSF, which is the Gaussian function f(x), and * represents the convolution operation; (5-5) Calculate the deconvolution residual: Convolve the restored original signal with the PSF, reintroduce the blurred or distorted part, and calculate the difference between the convolution result and the interpolated data in step (5-2); Step 6: Layered display and analysis (6-1) Moving average: For each data point y i Perform a moving average: Among them, x j is the signal strength value of the current A-scan at depth position j after deconvolution, y i is the signal strength value at depth position i after smoothing, and m represents half the size of the moving window; (6-2) Layered processing of column data: Manually select the boundaries on both sides of the peaks in the NFL, IPL, OPL, ELM, ISOS, ROST, and RPE-BrM layers, and manually select the boundaries on both sides of the valleys in the CC layer; initialize the movement amount and starting point, calculate the maximum value of the peak or valley, calculate the movement amount and update the boundary, and use the coefficient of the previous Gaussian fitting result as the new starting point for fitting; Among them, the NFL, IPL, OPL, ELM, ISOS, and CC layers use single-peak Gaussian fitting, and the ROST and RPE-BrM layers use double-peak Gaussian fitting; Single-peak Gaussian fitting: Set the initial value and upper and lower bounds of the peak height, center position, and standard deviation in the Gaussian model respectively, and use the nonlinear least squares method to perform a single-peak Gaussian fitting on the data to ensure that the fitting result is closer to the actual signal characteristics. The selected single-peak Gaussian fitting model formula is: Where a is the height of the peak; b is the center position of the peak in this layer; c is the standard deviation, which determines the width of the peak; Bimodal Gaussian fitting: Set the initial values ​​and upper and lower bounds of the peak height, center position, and standard deviation in the Gaussian model respectively, and use the nonlinear least squares method to perform bimodal Gaussian fitting on the data to obtain the characteristic parameters of the two main peaks of the signal to ensure that the fitting result is closer to the actual signal characteristics. The corresponding formula is: Where a1 is the height of the first peak, a2 is the height of the second peak; b1 is the center position of the first peak, b2 is the center position of the second peak; c1 is the width of the first peak, c2 is the width of the second peak; (6-3) Obtaining fitting data: Based on the fitting results, the peak and valley intensity, position and standard deviation of different layers are obtained over time after light stimulation, and the multi-layer photophysiological function signal of the ORG retina is extracted, laying the foundation for further data analysis.

Images