An OCTA image optimization method and device based on GPU real-time processing

Abstract

The present invention discloses an OCTA image optimization method and device based on GPU real-time processing. An optical coherence detection method is used to collect OCT scattering signals of scattering signal samples in three-dimensional space; a real-time blood flow optimization imaging method is based on GPU to construct an optimization classifier of a two-dimensional feature space between the inverse of the signal-to-noise ratio and the decorrelation coefficient of the OCT scattering signal, and the optimization classifier is used to perform real-time optimization classification of dynamic blood flow signals and static tissue signals, thereby optimizing and obtaining OCTA images; by constructing a data processing framework based on graphics processing unit (GPU) accelerated computing, high-speed processing and real-time imaging of OCTA data are realized; and the OCTA image quality is calculated in real time; data processing and imaging parameters are adjusted through current image quality feedback, and the OCTA image with the best quality is obtained in real time. The present invention realizes real-time processing and imaging of OCTA data and real-time self-optimization of OCTA blood flow projection image quality through GPU accelerated computing.

Description

Technical Field

[0001] The present invention generally relates to the field of biomedical imaging, and more particularly to real-time image processing and image quality optimization methods associated with Optical Coherence Tomography (OCT) and OCT Angiography (OCTA). Background Art

[0002] Hemodynamic response is an important parameter for measuring physiological function and pathological state, and optical coherence tomography (OCTA) technology can replace traditional exogenous fluorescent markers with endogenous blood flow movement to perform non-invasive, high-resolution, label-free three-dimensional imaging of microvascular networks in biological tissues, which can be used as an effective means to observe hemodynamic response. In clinical experiments, real-time OCTA imaging and high-quality blood flow angiography are reliable prerequisites for analyzing hemodynamic response.

[0003] In order to obtain OCTA blood flow images, it is usually necessary to perform repeated sampling (repeated A-line scans or B-frame scans) at a certain time interval at each spatial position of the biological tissue. The motion intensity at each signal is quantified by analyzing the temporal dynamics of the OCT scattering signal, and the blood flow signal and the static tissue signal are classified according to the quantified motion intensity. The OCTA blood flow classification reported so far is mainly based on the difference, variance or decorrelation calculation between adjacent A-line scans (or adjacent B-scan frames). Among them, the OCTA blood flow classification based on decorrelation calculation has a higher reliability of classification results in theory because it fully utilizes the statistical characteristics of multiple signals in the window. At the same time, since decorrelation measures the similarity between adjacent B-scan frames, it is less affected by changes in the overall light source intensity.

[0004] However, the quantification effect of decorrelation on motion contrast is significantly dependent on the noise level of the original OCT scattering signal. As the signal intensity decays (such as in deep tissue areas), random noise will gradually occupy the main component, and will also produce a large decorrelation value, resulting in decorrelation artifacts. The motion contrast generated based on the decorrelation operation cannot distinguish between the randomness of the noise and the decorrelation caused by the movement of red blood cells. Therefore, areas with weak signal-to-noise ratios are easily misjudged as blood flow signal areas, which seriously affects the contrast of blood flow images. The existing method based on the inverse SNR-decorrelation coefficient (ID) feature space uses the 3σ boundary of the static signal in the inverse SNR-decorrelation coefficient feature space as a fixed classification boundary. Although the ID-OCTA algorithm can remove most of the static areas, it also removes the dynamic signals that overlap with the static signals in the inverse SNR-decorrelation coefficient feature space, and cannot achieve accurate classification of dynamic signals. The quality of the generated OCTA images has room for further improvement. At the same time, the ID-OCTA method under the CPU framework has a slow processing speed because it cannot process all OCT signals in parallel. Summary of the invention

[0005] In view of the shortcomings of the prior art, the present invention proposes an OCTA image optimization method and device based on GPU real-time processing. In the method, a processing framework is constructed using GPU, and the quality parameters of the blood flow projection image are calculated and generated in real time by constructing an image quality evaluation module. The linear classifier threshold in the inverse signal-to-noise ratio-decorrelation coefficient space is adjusted according to the blood flow projection image quality feedback, thereby realizing real-time processing of OCTA signals and image quality optimization. The present invention can significantly accelerate the OCTA signal processing speed, optimize the quality of blood flow projection images, improve the contrast between blood flow and background, and improve vascular continuity.

[0006] The objective of the present invention is achieved through the following technical solutions:

[0007] 1. An OCTA image optimization method based on GPU real-time processing:

[0008] An optical coherence detection method for collecting OCT scattering signals of scattering signal samples in three-dimensional space;

[0009] A real-time blood flow optimization imaging method is proposed. Based on GPU, an optimization classifier is constructed in the two-dimensional feature space between the inverse of the signal-to-noise ratio and the decorrelation coefficient of the OCT scattering signal. The optimization classifier is used to perform real-time optimization classification of dynamic blood flow signals and static tissue signals, thereby optimizing the OCTA image.

[0010] The optical coherence detection method is one of the following:

[0011] Time domain OCT imaging method by scanning and changing the optical path length of the reference arm;

[0012] or a spectral domain OCT imaging method using a spectrometer to record spectral interference signals;

[0013] A swept frequency OCT imaging method that uses a detector to record spectral interference signals.

[0014] In the real-time blood flow optimization imaging method, the step of using GPU real-time optimization specifically includes:

[0015] S1. Using first-order autocovariance and zero-order autocovariance analysis to calculate and analyze the OCT scattering signal, the two characteristics of the inverse signal-to-noise ratio and decorrelation coefficient of each OCT scattering signal are obtained (21).

[0016] S2. Construct the inverse signal-to-noise ratio-decorrelation coefficient feature space of the OCT scattering signal (22). Determine the upper bound of the static tissue signal distribution and the lower bound of the dynamic blood flow signal distribution in the inverse signal-to-noise ratio-decorrelation coefficient feature space based on the multivariate time series model and the L'Aida criterion. Divide the inverse signal-to-noise ratio-decorrelation coefficient feature space into three regions according to the upper and lower bounds, namely, a dynamic region dominated by dynamic signals, a static region dominated by static signals, and a transition region where dynamic and static signals are mixed (23).

[0017] S3. Find an optimized classification line in the transition region so that the blood flow angiography image generated by the dynamic signal extracted by the classification line has the highest quality. Further separate the signal from the transition region according to the optimized classification line, and then combine it with the signal in the dynamic region to obtain a dynamic OCTA signal, and use the dynamic OCTA signal to generate an optimized blood flow angiography image (24).

[0018] The step S1 specifically comprises: calculating the inverse of the signal-to-noise ratio of each OCT scattering signal, calculating the decorrelation coefficient of the amplitude part of the OCT scattering signal obtained by scanning T different time points at the same spatial position or the OCT scattering signal including amplitude and phase, and obtaining two characteristics of the inverse of the signal-to-noise ratio and the decorrelation coefficient of each OCT scattering signal.

[0019] The step S2 specifically includes: obtaining the asymptotic distribution of the OCT scattering signal in the inverse signal-to-noise ratio-decorrelation coefficient feature space according to the multivariate time series model, obtaining the variance of the asymptotic distribution by simulation analysis, and determining the upper bound of the static signal distribution and the lower bound of the dynamic signal distribution.

[0020] The specific definition of the multivariate time series model is:

[0021] The collected OCT scattering signal X(s, t) is composed of the backscattered light signal A(s, t) of the biological tissue sample and the complex Gaussian white noise n(s, t), where the backscattered light signal is A(s, t) = αA d (s,t)+(1-α)A s (s, t), A d (s, t) and A s (s, t) are the dynamic and static components in the backscattered light signal, respectively. α adjusts the ratio of the two. s and t represent the three-dimensional space and time of the voxel, respectively. Generally, the backscattered light signal A and the noise n satisfy the following conditions: 1. n(s, t) is white noise that is independent of the time and space dimensions and has nothing to do with the backscattered light A(s, t); 2. In the time dimension, the dynamic part of the backscattered light A d (s, t) follows independent and identical distribution, the static part A s (s, t) is a constant vector.

[0022] Through multivariate time series model analysis, for the OCT scattering signal with a dynamic degree of α caused by motion, its asymptotic distribution in the inverse signal-to-noise ratio-decorrelation coefficient feature space satisfies:

[0023]

[0024] Where, → represents convergence, iSNR represents the inverse of the signal-to-noise ratio, D represents the decorrelation coefficient, and α represents the dynamic degree of the OCT scattering signal.

[0025] The variance of the asymptotic distribution obtained by simulation analysis is specifically: when the average kernel size for calculating the decorrelation coefficient is finite, the distribution variance σ of the static signal in the inverse signal-to-noise ratio-decorrelation coefficient feature space under α=0 is further obtained by simulation. s , satisfying the following formula:

[0026]

[0027] Among them, G is a constant obtained from simulation analysis, N is the average kernel size for calculating the decorrelation coefficient, N = S × T, S and T represent the three-dimensional spatial size and time size of the voxel respectively.

[0028] The upper bound of the static signal distribution and the lower bound of the dynamic signal distribution are determined as follows: the upper bound D of the static signal distribution is obtained according to the following formula: s1 and the lower bound D of the dynamic signal distribution d2 :

[0029]

[0030]

[0031] Where α0 is the lower limit of significant dynamics.

[0032] In S2, in the inverse signal-to-noise ratio-decorrelation coefficient feature space, the decorrelation value is greater than the upper bound D of the static signal distribution. s1 The area is the dynamic area R dominated by dynamic signals. d , the lower bound D between the dynamic signal distribution d2 and the upper bound D of the static signal distribution s1 The area between them is the transition area R where dynamic and static signals are mixed. i , the decorrelation value is less than the lower bound D of the dynamic signal distribution d2 The area is regarded as the static area R dominated by static signals. s .

[0033] The step S3 is specifically as follows:

[0034] S31, determining two dividing lines according to two classifier thresholds, respectively, extracting signals whose signal-to-noise ratio reciprocals are less than the two dividing lines corresponding to the two classifier thresholds to obtain two dynamic OCTA signals, and using the two dynamic OCTA signals to generate a blood flow angiography image respectively, thereby obtaining two blood flow angiography images (24);

[0035] S32, calculating the quality assessment parameter AQI of the two blood flow angiography images (24);

[0036] S33, compare the quality evaluation parameter AQI of the two blood flow angiography images (24) and iteratively update the classifier threshold, use the updated classifier threshold to return to step S31 to re-obtain the two blood flow angiography images, and then proceed to step S32 to re-obtain the quality evaluation parameter AQI of the two blood flow angiography images (24), thereby continuously repeating the steps iteratively until the quality evaluation parameter AQI is maximized, and then take the two classifier thresholds corresponding to the maximum quality evaluation parameter AQI;

[0037] S34. Calculate the average value of the two classifier thresholds corresponding to the maximum quality evaluation parameter AQI, segment the transition area with the dividing line determined by the average value, and use the area in the transition area where the decorrelation value is greater than the dividing line as the separated signal. Superimpose the separated signal and the signal in the dynamic area to obtain the dynamic OCTA signal.

[0038] GPU accelerated computing is used in the above process to achieve real-time optimization of the quality of blood flow angiography images.

[0039] In step S31, the upper boundary of the static tissue signal distribution and the lower boundary of the dynamic blood flow signal distribution obtained in step S2 are used as the two initial dividing lines, and the angle between the two initial dividing lines and the decorrelation coefficient coordinate axis in the signal-to-noise ratio inverse-decorrelation coefficient feature space is used as the two initial classifier thresholds.

[0040] The signals of the initial two dividing lines whose inverses of the signal-to-noise ratio are respectively less than the initial two classifier thresholds are extracted to obtain two dynamic OCTA signals respectively, and an initial blood flow angiography image is generated respectively using the two dynamic OCTA signals to obtain two initial blood flow angiography images.

[0041] In the step S32, the parameters of the blood vessel density and image contrast of the two blood flow angiography images (24) are calculated in real time using the GPU, and the product of the two is used as the quality evaluation parameter AQI of the blood flow projection image.

[0042] The blood vessel density and image contrast parameters of the blood flow angiography are calculated as follows:

[0043] Analyze the two-dimensional blood flow projection map A. Let the side length of the blood flow projection map be L. and Select four straight lines that are perpendicular to each other, calculate the gray intensity difference of each adjacent pixel on each straight line, and form a pixel intensity gradient curve. The number of intersections between the intensity gradient curve of the i-th straight line and the horizontal axis (pixel position) is V. i , and obtain the vascular density VD of the blood flow projection map:

[0044]

[0045] The 2D blood flow projection images present sharp grayscale intensity variations at the locations of blood vessels, and the frequency of such variations is captured by the proposed vessel density VD.

[0046] The image contrast of blood flow projection image is expressed by the root mean square contrast C RMS To express:

[0047]

[0048] Among them, N A Represents the total number of pixels in the blood flow projection map A, I (x，y) Represents the grayscale intensity of the pixel at the coordinate (x, y) in the blood flow projection map A, Represents the average gray intensity of blood flow projection image A.

[0049] Finally, the vascular density VD and the root mean square contrast C RMS The product of is the quality evaluation parameter AQI of the blood flow projection image A:

[0050] AQI=VD×C RMS .

[0051] In step S33, the classifier threshold is iteratively updated, specifically:

[0052] The two classifier thresholds are the first classifier threshold γ1 and the first classifier threshold γ2, respectively. The quality evaluation parameters AQI of the blood flow projection image obtained corresponding to the first classifier threshold γ1 and the first classifier threshold γ2 are the first quality evaluation parameter AQI1 and the second quality evaluation parameter AQI2, respectively, and the following judgment is performed:

[0053] If the first quality evaluation parameter AQI1<the second quality evaluation parameter AQI2, the first classifier threshold γ1 and the first classifier threshold γ2 of the next iteration are respectively set to and AQI2; otherwise, the first classifier threshold γ1 and the first classifier threshold γ2 of the next iteration are set to AQI1 and

[0054] In the step S33, after a certain iteration, if the difference between the first classifier threshold γ1 and the first classifier threshold γ2 is less than 1, it is considered that the quality evaluation parameter AQI is the largest at this time.

[0055] 2. Optical coherence angiography device based on GPU real-time processing:

[0056] An OCT optical coherence tomography detection device, used for collecting OCT scattering signals of scattering signal samples in three-dimensional space;

[0057] One or more GPU-accelerated signal processors are used to acquire and analyze the inverse signal-to-noise ratio and decorrelation coefficient of the OCT scattering signal, thereby constructing an optimized classifier in the two-dimensional feature space between the inverse signal-to-noise ratio and the decorrelation coefficient. The optimized classifier is used to perform real-time optimized classification of dynamic blood flow signals and static tissue signals, thereby optimizing the OCTA image.

[0058] In the specific implementation, a flow projection image quality evaluation module is also constructed internally, which is used for feedback adjustment and update of the threshold of the dynamic signal linear classifier in the inverse signal-to-noise ratio-decorrelation coefficient feature space.

[0059] The OCT optical coherence tomography device is one of the following:

[0060] Includes low coherence light source, interferometer and detector;

[0061] or including a low coherence light source, an interferometer, and a spectrometer;

[0062] Or include a swept wide spectrum light source, an interferometer and a detector.

[0063] The OCT optical coherence tomography detection device is selectively equipped with a visible light indicating device for indicating the position of the OCT detection beam and guiding the placement position of the detection target.

[0064] The visible light indicating device is mainly composed of a visible light indicating light source and a collimating lens.

[0065] The present invention realizes high-speed OCTA data processing and real-time imaging by constructing a data processing framework based on graphics processing unit (GPU) accelerated computing; and calculates OCTA image quality in real time; and adjusts data processing and imaging parameters through current image quality feedback to obtain OCTA images with optimal quality in real time.

[0066] The present invention realizes real-time processing and imaging of OCTA data and real-time self-optimization of OCTA blood flow projection image quality through GPU accelerated computing.

[0067] The beneficial effects and innovative features of the present invention are as follows:

[0068] The present invention calculates the inverse signal-to-noise ratio and decorrelation coefficient of the OCT scattering signal to form a two-dimensional inverse signal-to-noise ratio-decorrelation coefficient feature space, and combines the time series model to divide the inverse signal-to-noise ratio-decorrelation coefficient feature space into a dynamic region, a transition region, and a static region. A blood flow projection image quality evaluation module is designed to optimize and adjust the linear classifier threshold parameters in the inverse signal-to-noise ratio-decorrelation coefficient feature space based on vascular density and image contrast feedback, thereby extracting dynamic blood flow signals. A GPU computing framework is constructed to accelerate all signal processing processes in the process, thereby realizing real-time processing of OCTA signals and optimizing the quality of blood flow projection images. Compared with existing methods, the present invention can significantly accelerate the OCTA signal processing speed, improve the accuracy of dynamic signal classification, and optimize the OCTA image quality.

[0069] The image quality optimization method of the present invention is data-driven, starting from the blood vessel density and contrast of the data itself, and can adaptively optimize each set of data to its best image quality.

[0070] Compared with the prior art, the present invention has the following significant advantages:

[0071] 1. OCTA based on decorrelation calculation, due to the dependence between the decorrelation coefficient of the OCT scattering signal and the signal intensity, the decorrelation artifacts introduced by random noise in the low signal-to-noise ratio area cannot be distinguished from the decorrelation introduced by blood flow movement. The common solution is to set an empirical classifier threshold for intensity masking, which is equivalent to distinguishing the signal with the decorrelation classifier threshold in the intensity-decorrelation coefficient feature space, and then distinguishing it with a signal-to-noise ratio inverse classifier threshold. The more complex dependence between the signal decorrelation coefficient and the signal-to-noise ratio makes the boundary between the actual blood flow signal and other signals very different from the intensity classifier threshold straight line, resulting in misclassification. However, the classifier proposed in the present invention is based on the time series model analysis of the static signal in the signal-to-noise ratio inverse-decorrelation coefficient feature space, and has the advantage of signal-to-noise ratio adaptation. In addition, through the blood flow projection image quality evaluation module, the linear classifier threshold is adjusted by feedback, thereby self-optimizing the classification effect.

[0072] 2. The image quality optimization method based on blood flow projection image quality evaluation proposed in the present invention can adaptively adjust the classifier threshold parameters of the linear classifier in the inverse signal-to-noise ratio-decorrelation coefficient feature space, thereby effectively subclassifying the dynamic and static signals in the transition area, while suppressing static signals and effectively extracting more sufficient dynamic signals.

[0073] 3. Compared with the existing methods, the GPU-accelerated computing framework constructed by the present invention can realize real-time OCTA signal processing and image quality optimization, and its processing speed performance is significantly better than other methods; the image quality-based adaptive optimization method of the inverse signal-to-noise ratio-decorrelation coefficient feature space linear classifier is performed. After verification by a large number of samples, it is significantly better than the traditional method in terms of image signal-to-noise ratio, contrast-to-noise ratio, and vascular continuity parameters. BRIEF DESCRIPTION OF THE DRAWINGS

[0074] Figure 1 is a schematic diagram of the method of the present invention;

[0075] Figure 2 is a schematic diagram of the device of the present invention;

[0076] Figure 3 is a schematic diagram of a device according to an embodiment of the present invention;

[0077] Figure 4 This is a flowchart of a signal processing framework based on GPU acceleration according to an embodiment of the present invention;

[0078] Figure 5 A result diagram of a threshold value of a classification classifier for adjusting the image quality evaluation parameter feedback of an exemplary embodiment of the present invention;

[0079] Figure 6 This is a graph showing the experimental results of retinal blood flow imaging in living mice according to an exemplary embodiment of the present invention;

[0080] Figure 7 This is a comparison chart of the analysis of the mouse retinal angiography results of the present invention and the results of the traditional method.

[0081] In the figure: 1- an optical coherence detection method; 2- a real-time blood flow optimization imaging method; 21- calculate the two features of the inverse signal-to-noise ratio and the decorrelation coefficient of the OCT signal; 22- construct the inverse signal-to-noise ratio-decorrelation coefficient feature space of the OCT scattering signal; 23- calculate the upper and lower bounds of the static signal distribution, and divide the inverse signal-to-noise ratio-decorrelation coefficient feature space into three parts: dynamic area, transition area, and static area; 24- find the optimal classification line in the transition area to extract the dynamic signal and generate a blood flow angiography image. DETAILED DESCRIPTION

[0082] The specific implementation of the present invention will be described in detail below in conjunction with the accompanying drawings, which form a part of the present invention. It should be noted that these descriptions and examples are only exemplary and cannot be understood as limiting the scope of the present invention. The scope of protection of the present invention is defined by the attached claims, and any changes based on the claims of the present invention are within the scope of protection of the present invention.

[0083] Embodiments of the present invention are as follows:

[0084] To facilitate understanding of the embodiments of the present invention, each operation is described as a plurality of discrete operations, but the order of description does not represent the order of implementing the operations.

[0085] In this description, the sample measurement space is represented by a three-dimensional coordinate system of xyz based on the spatial direction. This description is only used to facilitate discussion and is not intended to limit the application of the embodiments of the present invention. Wherein: the depth z direction is the direction along the incident light axis; the xy plane is a plane perpendicular to the optical axis, where x is orthogonal to y, and x represents the OCT transverse fast scanning direction, and y represents the slow scanning direction.

[0086] The above-mentioned α, A, N, etc. represent variables, which are only used to facilitate discussion and are not intended to limit the application of the embodiments of the present invention, and can be any value such as 1, 2, 3, etc.

[0087] The method of the present invention is as follows Figure 1 As shown, in the signal acquisition part, OCT three-dimensional scanning imaging is performed on the tissue sample, and the same or adjacent spatial position is repeatedly sampled at T different time points 1.

[0088] In the signal classification part, a two-dimensional feature space is constructed by combining the scattering signal intensity and time decorrelation to achieve the classification of dynamic blood flow signals and static tissue background2.

[0089] The specific steps are: firstly, the mouse retinal sample is sampled in time and space by the collector, and the collected original spectral signal is reconstructed by OCT, including interpolation, dispersion correction, fast Fourier transform (FFT), and inter-frame correlation registration, as shown in Figure 4 The OCT part of the flowchart is shown in the figure. Then, the reconstructed OCT complex signal matrix is ​​subjected to OCTA processing, including stratification, constructing the inverse signal-to-noise ratio-decorrelation coefficient feature space classifier to extract blood flow signals, as shown in Figure 4 The image quality evaluation module is redesigned to evaluate the quality of the OCTA projection images generated in the previous step, and the classifier threshold of the linear classifier of the inverse signal-to-noise ratio-decorrelation coefficient feature space in the OCTA part is adjusted according to the evaluation parameter feedback, as shown in Figure 4 The blood flow projection image quality evaluation and linear ID classifier threshold optimization part of the flowchart are shown in the figure. Finally, the three-dimensional OCT / OCTA processing results are saved to generate a self-optimized OCTA blood flow projection map. The signal processing processes in the above steps are all carried out in the data processing framework based on GPU accelerated computing to achieve real-time OCTA signal processing and blood flow projection image quality optimization.

[0090] The calculation of the decorrelation coefficient used to construct the inverse signal-to-noise ratio-decorrelation coefficient feature space is based on analyzing the relative motion of the blood flow signal and the surrounding tissue to obtain the correlation degree information of each OCT scattering signal. The decorrelation operation includes calculating the decorrelation coefficient for the amplitude part of the OCT scattering signal obtained by scanning at T different time points, or for the complex OCT scattering signal containing amplitude and phase. The correlation calculation of the complex signal can theoretically obtain a higher motion contrast. The correlation coefficient between adjacent B scan frames (xz plane) can be calculated by the following formula:

[0091]

[0092] Among them, CorrMap i (z, x) represents the correlation coefficient of the voxel calculated at (z, x), i represents the i-th B scan frame, z represents the depth direction, x represents the line scan direction, A i (,) represents a complex signal, * represents the conjugate of a complex number, and is the complex intensity of two adjacent B scan frames. P and Q represent the size of the window in the z direction and x direction when calculating the correlation coefficient. Here, the window size is selected as 5×1 (z×x) pixels, and p and q represent the pixel positions in the corresponding window. The range of the calculated correlation coefficient is 0 to 1, and it can be converted into a decorrelation coefficient by subtracting 1 from the correlation coefficient. During the calculation process, each window corresponds to a decorrelation coefficient value. By moving the window to traverse the entire zx plane, the decorrelation values ​​of all positions in the entire plane can be obtained. The decorrelation values ​​of adjacent frames at the same scanning position are averaged as the result.

[0093] In OCT systems, the main source of noise is shot noise, which is considered to be approximately constant throughout the scanning volume and can be determined by scanning air data or by calibrating the system in advance.

[0094] The inverse signal-to-noise ratio iSNR of each voxel is then calculated using the following formula, which is defined as follows:

[0095]

[0096] Among them, s 2 is the noise level of the OCT system, and I represents the OCT signal intensity value of each voxel.

[0097] The noise level is obtained by calculating the intensity of the data blank area in the OCT scattering signal or by pre-calibrating the OCT system. So far, the calculation of the OCT signal decorrelation coefficient and the inverse of the signal-to-noise ratio has been completed 21.

[0098] The inverse of the OCT signal-to-noise ratio obtained by the above signal processing and the OCTA blood flow information obtained based on decorrelation calculation are combined to establish the ID two-dimensional feature space, and the OCT scattering signal is projected into the feature space 22. According to the results of the multivariate time series model analysis, for the OCT scattering signal caused by motion, its asymptotic distribution in the inverse of the signal-to-noise ratio-decorrelation coefficient feature space satisfies:

[0099]

[0100] Where, → represents convergence, iSNR represents the inverse of the signal-to-noise ratio, α represents the dynamic degree of the OCT scattering signal, and D represents the decorrelation coefficient.

[0101] When the average kernel size is limited, the distribution variance σ of the static signal (α = 0) in the inverse signal-to-noise ratio-decorrelation coefficient feature space is further obtained through simulation. s , satisfying the following formula:

[0102]

[0103] Among them, G≈1.5 is a constant obtained from simulation analysis, N=S×T is the Gaussian average kernel size, and S and T represent the three-dimensional spatial size and time size of the voxel, respectively.

[0104] Then, based on the multivariate time series model and simulation analysis, combined with empirical values, the upper bound D of the static signal distribution is defined in the inverse signal-to-noise ratio-decorrelation coefficient feature space. s1 and the lower bound D of the dynamic signal distribution d2 , defined as follows:

[0105]

[0106]

[0107] Where α0 is the lower limit of significant dynamics, and its value in this example is 0.4.

[0108] According to the definitions of equations (5) and (6), the SNR inverse-decorrelation coefficient feature space can be divided into three parts23: the static region R dominated by static signals s , dynamic region R dominated by dynamic signals d , the transition region R of the dynamic and static signal mixture t ,like Figure 5 As shown in (a).

[0109] When the dynamic signal in the inverse signal-to-noise ratio-decorrelation coefficient feature space is initially classified, the classifier threshold used is D s1 With D d2 , to ensure that in the subsequent iterative optimization process of the classifier threshold, the linear classifier is always in the transition region, so as to extract appropriate dynamic signals in the transition region. s1 With D d2 The classifier classifies the mixed signals in the transition area, extracts the dynamic signals therein, and generates a corresponding blood flow angiography image 24.

[0110] OCTA image quality evaluation parameters are constructed based on the vascular density and image contrast of blood flow angiography images, such as Figure 5 As shown in each two-dimensional blood flow diagram in (b), the two-dimensional blood flow projection diagram A is analyzed, and the side length of the blood flow projection diagram is L. and Select four mutually perpendicular straight lines at each position, calculate the gray intensity difference of each adjacent pixel on each straight line, and form a pixel intensity gradient curve, such as Figure 5 As shown in (c), the number of intersections between the intensity gradient curve of the ith straight line among the four straight lines and the horizontal axis (pixel position) is V i , the vascular density VD of the blood flow projection map can be obtained:

[0111]

[0112] The two-dimensional blood flow projection image shows sharp grayscale intensity changes at the blood vessel position. The frequency of such changes is captured by the given blood vessel density VD, which can reflect the number of blood vessels in the blood flow projection image.

[0113] The image contrast of blood flow projection image is expressed by the root mean square contrast C RMS To express:

[0114]

[0115] Among them, N A Represents the total number of pixels in the two-dimensional blood flow projection image A, I (x，y) Represents the grayscale intensity of the pixel at the coordinate (x, y), Represents the average gray intensity of the blood flow projection image.

[0116] The vascular density VD and the root mean square contrast C RMS The product of is the quality evaluation parameter AQI of the blood flow projection image:

[0117] AQI=VD×C RMS (9)

[0118] Figure 5 (b) shows the blood flow images generated by the same OCT data after classification by different linear ID classifiers γ, and their image quality AQI values. In the transition area of ​​the inverse signal-to-noise ratio-decorrelation coefficient feature space, a linear ID classifier corresponding to a certain classifier threshold can extract dynamic signals more completely and generate OCTA projections with the best image quality. Blood flow projections generated by different ID classifiers for the same data, with vascular density VD, root mean square contrast C RMS The relationship between the blood flow projection image quality evaluation parameter AQI and the classifier threshold of its corresponding classifier is as follows: Figure 5 As shown in (d), the vascular density VD increases with the increase of the classifier threshold γ. It first increases because more dynamic signals are extracted. After the dynamic signals are extracted, it fluctuates around its maximum value due to the influence of static tissue and noise. RMS As the classifier threshold γ increases, it first increases and then decreases. The blood flow projection image quality evaluation parameter AQI obtained by multiplying the two also tends to first increase and then decrease as the classifier threshold γ increases.

[0119] Under the same data conditions, the ID classifier threshold corresponding to the maximum blood flow projection image quality evaluation parameter AQI can be determined by an iterative search method based on . Specifically, let the two starting classifier thresholds in a certain iteration be γ1 and γ2, and the corresponding blood flow projection image quality evaluation parameter AQI are AQI1 and AQI2 respectively. If AQI1<AQI2, then the starting classifier threshold of the next iteration is with AQI2, otherwise with AQI1 Until the difference of the initial classifier threshold after a certain iteration is less than 1, the ID classifier threshold corresponding to the maximum blood flow projection image quality evaluation parameter AQI is found, thereby achieving self-optimization of OCTA image quality.

[0120] The inverse signal-to-noise ratio-decorrelation coefficient feature space is classified and the dynamic signal is extracted according to the classifier threshold optimized based on the blood flow projection image quality evaluation feedback. Gaussian filtering or Hessian filtering can be performed from the graphics perspective to further obtain the blood flow angiography results.

[0121] The above signal processing processes are all in Figure 4 The processing flow framework based on GPU accelerated computing is shown in the example. The GPU model used in the example is RTX 2080Ti (32GB video memory), which reaches a line processing speed of 133kHz. For a collector (1) with a line acquisition speed of 120kHz, real-time OCTA signal processing and image quality optimization can be achieved.

[0122] Figure 2The structure diagram of the acquisition device of the OCT blood flow imaging technology based on GPU real-time processing in the present invention is shown. The main structure of the low-coherence interference measurement part of the device is an interferometer, which is composed of 11 to 23, wherein the light emitted by the light source 11 is divided into two light beams by the beam splitter 12: one of the light beams enters the reference arm of the interferometer through a polarization controller 13, and is irradiated on the plane reflector 15 through the reference arm collimator 14; the other light beam enters the sample arm through another polarization controller 13, and is focused on the sample to be tested 21 through the collimator lens 16 and the scanning device optical path. In the scanning device optical path, the light beam is reflected by the two-dimensional scanning galvanometer group 17, 18, the "4f" lens group 54, 55 and the dichroic mirror 19, and then is focused on the sample to be tested 21 through the focusing objective lens 20. The design of the lens group 54, 55 is to ensure that the center of the light beam of the mirror surface of the two-dimensional scanning galvanometer and the center of the light beam of the reflection surface of the dichroic mirror are fixed and unchanged during scanning, so that the light beam in the OCT sample arm does not affect the imaging properties of the objective lens during scanning. Then, the light reflected from the reference arm and the sample arm interferes and is received by the interference signal detection device 22, which is then connected to the signal processor module and the calculation unit 23. For the optical fiber optical path, the polarization controller 13 is used to adjust the polarization state of the light beam to maximize the signal interference effect.

[0123] The specific implementation is also provided with a visible light indicating device, which includes a low-power visible light source 25, a collimating lens 24 and a filter 52. The visible light used for indicating passes through the collimating lens 24, the dichroic mirror 19 and the focusing objective lens 20 to the sample 21 to be tested.

[0124] According to different ways of detecting signals by low coherence interference, Figure 2 The acquisition device of an OCT blood flow imaging technology based on GPU real-time processing specifically includes:

[0125] 1) Time domain measurement device. The light source 11 uses broadband low-coherence light, the plane reflector 15 can move along the optical axis, and the interference signal detection device 22 is a point detector. By moving the plane reflector 15 to change the optical path of the reference arm, the interference signals of the two arms are detected by the point detector 22, and the low-coherence interference detection of the scattered signal in the z direction of a certain spatial depth is performed, thereby obtaining a sampling volume of the depth space dimension.

[0126] 2) Spectral domain measurement device. The light source 11 uses broadband low-coherence light, the plane reflector 15 is fixed, and the interference signal detection device 22 uses a spectrometer. The interference signal passes through the linear array camera in the spectrometer and records the interference spectrum at the same time. The Fourier analysis method is used to analyze the interference spectrum signal, and the scattering information in the depth z direction is obtained in parallel, so as to obtain a sampling volume in the depth space dimension.

[0127] 3) Sweep frequency measurement device. The light source 11 uses a swept frequency light source, the plane reflector 15 is fixed, and the interference signal detection device 22 uses a point detector. The point detector records the low-coherence interference spectrum of the swept frequency light source in time-sharing. The Fourier analysis interference spectrum signal is sampled, and the scattering information in the depth z direction is obtained in parallel, so as to obtain a sampling volume of the depth space dimension.

[0128] Figure 3 The figure shows an exemplary embodiment of the present invention. The feature space-based three-dimensional blood flow imaging system includes a broadband low-coherence light source 26, an optical circulator 27, a fiber coupler 28 with a splitting ratio of 50:50, a first polarization controller 29, a first fiber collimator 30, a focusing lens 36, a plane reflector 37, a second polarization controller 38, a second fiber collimator 39, a two-dimensional scanning galvanometer assembly 40 and 41, a dichroic mirror 42, a focusing objective lens 43, a third fiber collimator 45, a grating 46, a focusing lens 47, a high-speed linear array camera 48, a signal processor module and a computing unit 49, a visible light indicator light source 50, a collimator lens 51, and a "4f" lens group 56 and 57, wherein the broadband low-coherence light source 26 adopts a superluminescent diode light source with a central wavelength of 1325nm and a bandwidth of 100nm, the focusing objective lens 43 adopts an achromatic double glued lens with a focal length of 30mm, and the high-speed linear array camera 48 adopts a 2048-pixel single The invention discloses a linear array scanning camera composed of two elements; wherein the light emitted by the low coherence broadband light source 26 used in the device of the present invention passes through an optical circulator 27 and enters a fiber coupler 28 with a splitting ratio of 50:50, and the light emitted from the fiber coupler 28 is divided into two sub-beams: one beam of light is connected to the first fiber collimator 30 in the reference arm through an optical fiber through a first polarization controller 29, and irradiates the plane reflector 37 after passing through a collimating and focusing lens 36; the other beam of light is connected to the second fiber collimator 39 in the sample arm through an optical fiber through a second polarization controller 38, and after being collimated, it is reflected by two scanning galvanometers 40, 41, a "4f" lens group 56, 57 and a dichroic mirror 42, and then is focused by a focusing objective lens 43 onto a sample 44 to be measured, wherein the design of the lens groups 56, 57 is to ensure that the beam center of the mirror surface of the two-dimensional scanning galvanometer and the beam center of the reflection surface of the dichroic mirror are fixed and unchanged during scanning. The light reflected by the plane reflector 37 in the reference arm interferes with the light backscattered by the sample under test in the sample arm at the fiber coupler 28. The interference light is detected and recorded by the spectrometers 45 to 48, and then collected and processed by the signal processor module and the calculation unit 49. The spectrometer includes devices 45 to 48 connected in sequence, wherein the device 45 is a fiber coupler, 46 is a grating, and 47 is a converging lens, which focuses the light dispersed by the grating onto the linear array detector shown in 48.

[0129] The specific implementation is also provided with a visible light indicating device, which includes a visible light indicating light source 50 and a collimating lens 51 . The visible light indicating light source 50 emits visible light for indicating, which passes through the collimating lens 51 , the dichroic mirror 42 and the focusing objective lens 43 and then reaches the sample 44 to be tested.

[0130] Figure 6 The figure is the result of the retinal blood flow imaging experiment in living mice. The OCTA image optimization method and device based on GPU real-time processing proposed in the present invention can accurately stratify the structural layers of the retina and generate high-quality retinal blood flow projection images of each layer, including the full retinal layer (RE), superficial retinal layer (SVP), intermediate retinal layer (ICP), deep retinal layer (DCP) and choroid layer (CH).

[0131] Figure 7 The figure is a comparison chart of the analysis of the mouse retinal angiography results of the present invention and the results of the traditional method. The traditional method uses the ID-OCTA method that uses the upper limit of the static signal in the inverse signal-to-noise ratio-decorrelation coefficient feature space distribution as the classification classifier threshold. The classifier threshold adaptive optimization method based on image quality evaluation proposed in the present invention ( Figure 7 (a)) and ID-OCTA( Figure 7 (b)) can obtain lower background noise, and the blood flow signal extracted by the present invention is richer, has a higher contrast relative to the background, and has better vascular continuity, indicating that the present invention can effectively extract dynamic signals while suppressing noise signals in the transition area. In ID-OCTA, since the upper bound of the static signal in the inverse of the signal-to-noise ratio-decorrelation coefficient feature space distribution is used as the classification classifier threshold, while retaining the complete dynamic signal, the possible dynamic signal in the transition area is completely eliminated, resulting in the removal of part of the blood flow signal, affecting the continuity of the blood vessel. Figure 7 (c) It reflects that the present invention extracts more OCTA signals than the ID-OCTA method in the inverse signal-to-noise ratio-decorrelation coefficient feature space.

[0132] Furthermore, in order to quantify the image quality improvement brought by the classifier threshold adaptive optimization method based on image quality evaluation relative to the ID-OCTA method, the present invention characterizes the performance of the two methods based on the image signal-to-noise ratio SNR, the contrast-to-noise ratio CNR and the vascular continuity CON. The noise ratio CNR is calculated by the following formula:

[0133]

[0134] in, are the mean values ​​of the decorrelation coefficients of the vascular region and the background noise region, σ n is the standard deviation of the background noise area.

[0135] Vascular continuity CON is calculated by the following formula:

[0136]

[0137] Among them, N c is the number of connection points after the projection graph is skeletonized, N s is the number of all pixels after the skeletonization of the projection image.

[0138] Table 1 shows the quantitative evaluation results of the mouse retinal full-layer blood flow projection map for 10 sets of data. The parameters of SNR, CNR, and CON are described in the form of mean ± standard deviation, and Δ represents the improvement of the parameter value.

[0139] Table 1 Quantitative comparison of blood flow imaging performance between the present invention and the ID-OCTA blood flow classification method

[0140]

[0141] The quantitative analysis results in Table 1 show that the blood flow angiography results obtained by the classifier threshold adaptive optimization method based on image quality evaluation proposed in the present invention have better signal-to-noise ratio, contrast-to-noise ratio, and vascular continuity of the generated blood flow projection map than the ID-OCTA method. This comparison result fully demonstrates that the OCTA image optimization method based on GPU real-time processing involved in the present invention can improve the accuracy of blood flow signal classification, achieve effective enhancement of blood flow contrast and increase signal processing speed, and has outstanding and significant technical effects.

Claims

1. An OCTA image optimization method based on GPU real-time processing, characterized in that include: An optical coherence detection method (1) for collecting OCT scattering signals of scattering signal samples in three-dimensional space; A real-time blood flow optimization imaging method (2), which constructs an optimization classifier of a two-dimensional feature space between the inverse of the signal-to-noise ratio and the decorrelation coefficient of the OCT scattering signal based on a GPU, and uses the optimization classifier to perform real-time optimization classification of dynamic blood flow signals and static tissue signals, thereby optimizing and obtaining OCTA images; The real-time blood flow optimization imaging method (2) specifically comprises: S1. Using first-order autocovariance and zero-order autocovariance analysis to calculate and analyze the OCT scattering signal, the two characteristics of the inverse signal-to-noise ratio and decorrelation coefficient of each OCT scattering signal are obtained (21). S2, constructing the inverse signal-to-noise ratio-decorrelation coefficient feature space of the OCT scattering signal (22), determining the upper bound of the static tissue signal distribution and the lower bound of the dynamic blood flow signal distribution in the inverse signal-to-noise ratio-decorrelation coefficient feature space, and dividing the inverse signal-to-noise ratio-decorrelation coefficient feature space into three regions according to the upper bound and the lower bound, namely, a dynamic region dominated by dynamic signals, a static region dominated by static signals, and a transition region where dynamic and static signals are mixed (23); S3, finding an optimized classification line in the transition region, further separating the signal from the transition region according to the optimized classification line, combining it with the signal in the dynamic region to obtain a dynamic OCTA signal, and using the dynamic OCTA signal to generate an optimized blood flow angiography image (24); The step S3 is specifically as follows: S31, determining two dividing lines respectively according to two classifier thresholds, extracting signals whose signal-to-noise ratio reciprocals are respectively less than the two dividing lines to obtain two dynamic OCTA signals, and using the two dynamic OCTA signals to generate a blood flow angiography image respectively, to obtain two blood flow angiography images (24); S32, calculating the quality assessment parameter AQI of the two blood flow angiography images (24); S33, compare the quality evaluation parameter AQI of the two blood flow angiography images (24) and iteratively update the classifier threshold, use the updated classifier threshold to return to step S31 to re-obtain the two blood flow angiography images, and then proceed to step S32 to re-obtain the quality evaluation parameter AQI of the two blood flow angiography images (24), thereby continuously repeating the steps iteratively until the quality evaluation parameter AQI is maximized, and then take the two classifier thresholds corresponding to the maximum quality evaluation parameter AQI; S34. Calculate the average value of the two classifier thresholds corresponding to the maximum quality evaluation parameter AQI, segment the transition area with the dividing line determined by the average value, and use the area in the transition area where the decorrelation value is greater than the dividing line as the separated signal. Superimpose the separated signal and the signal in the dynamic area to obtain the dynamic OCTA signal.

2. The OCTA image optimization method based on GPU real-time processing according to claim 1, characterized in that: The optical coherence detection method (1) is one of the following: Time domain OCT imaging method by scanning and changing the optical path length of the reference arm; or a spectral domain OCT imaging method using a spectrometer to record spectral interference signals; A swept frequency OCT imaging method that uses a detector to record spectral interference signals.

3. The OCTA image optimization method based on GPU real-time processing according to claim 1, characterized in that: In S2, in the inverse signal-to-noise ratio-decorrelation coefficient feature space, the decorrelation value is greater than the upper bound D of the static signal distribution. s" The area is the dynamic area R dominated by dynamic signals. $ , the lower bound D between the dynamic signal distribution $％ and the upper bound D of the static signal distribution s" The area between them is the transition area R where dynamic and static signals are mixed. i , the decorrelation value is less than the lower bound D of the dynamic signal distribution $％ The area is regarded as the static area R dominated by static signals. s .

4. The OCTA image optimization method based on GPU real-time processing according to claim 1, characterized in that: In step S31, the upper boundary of the static tissue signal distribution and the lower boundary of the dynamic blood flow signal distribution obtained in step S2 are used as the two initial dividing lines, and the angle between the two initial dividing lines and the decorrelation coefficient coordinate axis in the signal-to-noise ratio inverse-decorrelation coefficient feature space is used as the two initial classifier thresholds.

5. The OCTA image optimization method based on GPU real-time processing according to claim 1, characterized in that: In the step S32, the parameters of the blood vessel density and image contrast of the two blood flow angiography images (24) are calculated in real time using the GPU, and the product of the two is used as the quality evaluation parameter AQI of the blood flow projection image.

6. The OCTA image optimization method based on GPU real-time processing according to claim 1, characterized in that: In step S33, the classifier threshold is iteratively updated, specifically: The two classifier thresholds are the first classifier threshold γ " With the first classifier threshold γ ％ , the first classifier threshold γ " With the first classifier threshold γ ％ The quality evaluation parameters AQI of the corresponding blood flow projection images are respectively the first quality evaluation parameter AQI " and the second quality evaluation parameter AQI ％ , and make the following judgments: If the first quality evaluation parameter AQI " <Second quality evaluation parameter AQI ％ , then the first classifier threshold γ of the next iteration " With the first classifier threshold γ ％ Set as with AQI ％ ; Otherwise, the first classifier threshold γ of the next iteration " With the first classifier threshold γ ％ Set as AQI " and 7. An optical coherence angiography device based on GPU real-time processing for implementing any of the methods of claims 1 to 6, comprising: An OCT optical coherence tomography detection device, used for collecting OCT scattering signals of scattering signal samples in three-dimensional space; One or more GPU-accelerated signal processors are used to acquire and analyze the inverse signal-to-noise ratio and decorrelation coefficient of the OCT scattering signal, thereby constructing an optimized classifier in the two-dimensional feature space between the inverse signal-to-noise ratio and the decorrelation coefficient. The optimized classifier is used to perform real-time optimized classification of dynamic blood flow signals and static tissue signals, thereby optimizing the OCTA image.

8. The optical coherence angiography device based on GPU real-time processing according to claim 7, characterized in that: The OCT optical coherence tomography device is one of the following: Includes low coherence light source, interferometer and detector; or including a low coherence light source, an interferometer, and a spectrometer; or includes a swept wide spectrum light source, an interferometer, and a detector; A visible light indicating device is selectively configured in the OCT optical coherence tomography detection device to indicate the position of the OCT detection beam and guide the placement position of the detection target; The visible light indicating device is mainly composed of a visible light indicating light source and a collimating lens.

Images