Distortion correction method for scanning aberration of OCT (optical coherence tomography) galvanometer

Abstract

According to the distortion correction method for the OCT galvanometer scanning aberration, a pixel displacement vector caused by aberration in a three-dimensional space is decomposed into the sum of a displacement vector in the Z direction and a displacement vector in the XY plane, the aberration is compensated successfully, and space correction of a whole three-dimensional image is completed. A second step of rough correction and Z-direction correction: compensating the original data to complete rough correction; 3, fine correction and correction in an XY plane: performing fine correction on the OCT signal after coarse correction, namely performing phase compensation on aberration in the XY plane; and 4, obtaining corrected data: for the compensation of different depths, placing the plane mirror at the corresponding depth z0, and repeating the steps 2 and 3 to finish the aberration compensation of the OCT signal of the whole xyz space. According to the method, aberration caused by galvanometer scanning is removed, the imaging quality is improved, and the resolution of the three-dimensional coordinates in all directions xyz is improved.

Description

Technical Field

[0001] The present invention relates to the technical field of optical coherence tomography distortion correction, and in particular to a method for correcting the distortion of an OCT galvanometer scanning aberration. Background Art

[0002] OCT is a non-contact, high-precision 3D scanning imaging technology with a wide range of applications in medical and industrial inspection. This technology images a sample by measuring the back-reflected light at different depths. Specifically, each scanning point can image A points at its location, which are located at different depths of the sample. This data line of length A is called an Aline, as shown in the following figure. Figure 1 As shown, attached Figure 1 A, B, and C shown in the figure represent the three mutually perpendicular axes of the scan. In the industry, the depth direction is generally referred to as the A-direction scan, and the other two directions are the B- and C-direction scans. Then, by rotating the two galvanometer mirrors to change the scanning position, the entire B×C area is scanned. Since A points can be obtained each time the scan is performed, a cube data of size A×B×C can be obtained in the end. Through the image algorithm, the corresponding three-dimensional scan image of the sample can be obtained.

[0003] Due to the inherent characteristics of galvanometer scanning, the light rays reflected by the galvanometer are not parallel to each other. Therefore, the scanned asymmetric surface is a curved surface related to the scanning angle. Due to the superposition of the two galvanometers, this surface may be a saddle surface or a sphere, so the image is distorted accordingly. In order to obtain accurate sample information, we must correct this distortion. Summary of the Invention

[0004] The present invention aims to solve at least one of the technical problems existing in the prior art.

[0005] To this end, the present invention proposes a distortion correction method for OCT galvanometer scanning aberration, which removes the aberration caused by galvanometer scanning, improves imaging quality, and enhances the resolution of each direction of the three-dimensional coordinates xyz.

[0006] The method for correcting distortion caused by OCT galvanometer scanning aberrations according to an embodiment of the present invention decomposes the pixel displacement vector caused by the aberration in three-dimensional space into the sum of the displacement vector in the Z direction and the displacement vector in the XY plane, successfully compensates for the aberration, and completes the spatial correction of the entire three-dimensional image. The method includes the following steps:

[0007] Step 1: Obtain raw data: Place a plane mirror sample near the focal plane and use an OCT device to image the plane mirror sample to obtain a curved image in three-dimensional space. The distribution of the plane mirror surface in three-dimensional space is calculated based on the light intensity distribution. The formula is as follows:

[0008] Z(x,y)=g(maxIdx(V(x,y,z),z)) (1)

[0009] The meanings of the symbols in formula (1) are:

[0010] Z(x,y) represents the offset of each point in the (x,y) plane in the Z direction;

[0011] g represents median filtering;

[0012] maxIdx(V(x,y,z),z) indicates the position of the maximum grayscale value in the Z direction;

[0013] By fitting Z(x,y) with a polynomial, we can obtain the curved surface z′(x,y) of the plane mirror in the Z direction.

[0014] Step 2, coarse correction, Z direction correction: Compensate the original data to complete the coarse correction. The formula for coarse correction is as follows:

[0015] D′(x,y,z)=D(x,y,z)*exp(-1j*Z′(x,y)*2*pi) (2)

[0016] The meanings of the symbols in formula (2) are:

[0017] D′(x,y,z) represents the OCT signal after coarse correction;

[0018] D(x,y,z) represents the raw OCT signal before coarse correction at depth z;

[0019] j represents the imaginary unit;

[0020] Z′(x,y) represents the curved surface of the plane mirror in the Z direction;

[0021] pi stands for pi;

[0022] x, y, z represent the three coordinate axes in a three-dimensional rectangular coordinate system;

[0023] Step 3: Fine correction, correction in the XY plane: Perform fine correction on the coarsely corrected OCT signal, that is, perform phase compensation for the aberration in the XY plane. The formula for fine correction is as follows:

[0024] D″(x,y,z0)=D′(x,y,z0)*exp(-1j*Ф(x,y,z0)) (3)

[0025] The meanings of the symbols in formula (3) are:

[0026] D″(x,y,z0) represents the OCT signal after fine correction;

[0027] D′(x,y,z0) represents the OCT signal after coarse correction at depth z0;

[0028] Ф(x,y,z0) represents the compensation phase at depth z0;

[0029] z0 represents a certain depth on the Z axis;

[0030] Step 4: Obtain corrected data: For different depths of z0, place the plane mirror at different depths of z0 and repeat steps 2 and 3 to complete the aberration compensation of the OCT signal in the entire xyz space.

[0031] The beneficial effects of the present invention are that the method has a simple process, does not require complex optical and mechanical modeling to derive the image plane distribution formula, and does not require calibration and measurement of relevant parameters; the method can remove aberrations at all positions within the imaging depth, rather than just the aberrations at the depth where the calibrated plane mirror is located; and digital correction is used, without the need to add hardware, saving costs.

[0032] According to one embodiment of the present invention, the calculation formula of the compensation phase Φ(x, y) is as follows:

[0033]

[0034] The meanings of the symbols in formula (6) are:

[0035] p k Indicates the linear superposition coefficient corresponding to the compensation phase in the XY plane;

[0036] S k represents the k-th Zernike polynomial.

[0037] According to one embodiment of the present invention, the p k The updating method is simulated annealing.

[0038] According to one embodiment of the present invention, the p k The update method is gradient descent.

[0039] According to one embodiment of the present invention, the specific steps of solving the Zernike polynomial coefficients are as follows:

[0040] Step 1, coefficient initialization;

[0041] Step 2: Update the coefficients;

[0042] Step 3: Calculate the compensation phase from the coefficient;

[0043] Step 4: Phase compensation;

[0044] Step 5: Calculate the mean square error after compensation;

[0045] Step 6: Determine whether the mean square error is less than the threshold:

[0046] When the mean square error is judged to be less than the threshold, the correction ends;

[0047] When it is determined that the mean square error is not less than the threshold, return to step 2.

[0048] According to one embodiment of the present invention, in the fourth step, the calculation formula of the phase compensation is as follows:

[0049] D_new=D_old*exp(j*phi(x,y) (7)

[0050] The meanings of the symbols in formula (7) are:

[0051] D_new represents the data after compensation;

[0052] D_old represents the data before compensation;

[0053] j represents the imaginary unit, that is,

[0054] According to one embodiment of the present invention, the plane mirror is a plane reflecting mirror.

[0055] According to an embodiment of the present invention, in the fourth step, the difference between two adjacent depths is 1 pixel.

[0056] According to one embodiment of the present invention, in step 5, the calculation formula of the mean square error (MSE) is as follows:

[0057]

[0058] The meanings of the symbols in formula (8) are:

[0059] N represents the total number of pixels on the xy plane;

[0060] y i (x, y) represents the actual observation value;

[0061] Represents the predicted value obtained by the regression model.

[0062] Other features and advantages of the present invention will be described in the following description, and in part will become apparent from the description, or understood by practicing the present invention. The purpose and other advantages of the present invention are realized and obtained by the structures particularly pointed out in the description and the drawings.

[0063] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, preferred embodiments are given below and described in detail with reference to the accompanying drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0064] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are only some embodiments recorded in this application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0065] Figure 1 This is a schematic diagram of the galvanometer scanning of the OCT system;

[0066] Figure 2 is the image at depth near the focal plane of the system;

[0067] Figure 3 It is the left side view of the plane mirror sample;

[0068] Figure 4 It is the front side view of the plane mirror sample;

[0069] Figure 5 is the fitted image in the Z direction, i.e. z'(x,y);

[0070] Figure 6 This is the image after rough correction of the OCT signal of the plane mirror (top view);

[0071] Figure 7 is the distribution of the first 36 terms of the Zernike polynomials on the XY plane;

[0072] Figure 8 is an overall flow chart of the distortion correction method of the present invention;

[0073] Figure 9 This is a schematic diagram of solving the Zernike polynomial coefficients. DETAILED DESCRIPTION

[0074] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts shall fall within the scope of protection of the present invention.

[0075] The following describes in detail the distortion correction method for OCT galvanometer scanning aberrations according to an embodiment of the present invention, with reference to the accompanying drawings. This method is divided into two parts: coarse correction and fine correction. Coarse correction is used to correct large-scale pixel offsets in the Z direction, i.e., Z-direction aberrations; fine correction is used to correct pixel offsets within the XY plane, i.e., XY-plane aberrations. Generally speaking, pixel offsets in the XY plane are much smaller than pixel offsets in the Z direction.

[0076] The distortion correction method for OCT galvanometer scanning aberrations of the present invention regards the spatial aberration correction of a three-dimensional image as the displacement of pixels in three-dimensional space, and decomposes the displacement field into the superposition of a one-dimensional displacement field (in the Z direction) and a two-dimensional displacement field (in the XY plane). That is, by decomposing the pixel displacement vector caused by aberration in three-dimensional space into the sum of the displacement vector in the Z direction and the displacement vector in the XY plane, the aberration is successfully compensated, completing the spatial correction of the entire three-dimensional image.

[0077] It should be noted that the distortion correction method for OCT galvanometer scanning aberrations of the present invention mainly includes the following two points: 1. Coarse correction with a plane mirror (z direction), which can be measured at different depths, and can be measured at the depth of interest; 2. Fine correction (xy direction), using the Zernike polynomial surface for phase compensation. The entire search for the compensation phase is an iterative process.

[0078] The method for correcting the distortion of the OCT galvanometer scanning aberration specifically includes the following steps:

[0079] Step 1: Obtain raw data: Place a plane mirror sample near the focal plane and use the OCT device to image the plane mirror sample to obtain a curved image in three-dimensional space. The three views (cross-sections) of the image are as follows: Figure 2 、 Figure 3 、 Figure 4 As shown, Figure 2 It is an image at a depth near the focal plane of the system, and is the reflecting surface of a plane mirror. Due to the influence of aberration, the brightness distribution is uneven, and the overall distribution is saddle-shaped. The reason for the saddle shape is the deformation caused by the scanning aberration of the galvanometer. If there is no scanning aberration, it will present a plane with uniform brightness.

[0080] Figure 3 and Figure 4 This is a side view of the sample. It can be seen that the image formed by the plane mirror after passing through the OCT imaging system is bent due to the influence of the galvanometer scanning aberration.

[0081] Specifically, Figure 3 This is the left side view of the plane mirror sample. The white in the picture is the plane mirror, and the black is the area outside the plane mirror where no light is reflected, that is, the air.

[0082] Figure 4This is the front and side view of the plane mirror sample. The white in the figure is the plane mirror, and the black is the area outside the plane mirror where no light is reflected, that is, the air.

[0083] It should be noted that the reason why a standard plane mirror is placed in the system for measurement is because the plane mirror is absolutely horizontal. The plane mirror can be used as a reference, and then the changes reflected in the image after the reference passes through the system are calculated.

[0084] First, the distribution of the plane mirror surface in three-dimensional space is calculated through the light intensity distribution, as shown in formula (1). Formula (1) is as follows:

[0085] Z(x,y)=g(maxIdx(V(x,y,z),z)) (1)

[0086] The meanings of the symbols in formula (1) are:

[0087] Z(x,y) represents the offset of each point in the (x,y) plane in the Z direction;

[0088] g represents median filtering;

[0089] maxIdx(V(x,y,z),z) indicates the position of the maximum grayscale value in the Z direction, that is, the position of the maximum grayscale value on each Aline;

[0090] By fitting Z(x,y) with a polynomial, we can obtain the curved surface Z′(x,y) of the plane mirror in the Z direction, as shown in Figure 5 shown; among them, Figure 5 Different colors in represent different grayscale values; Figure 5 The horizontal axis represents pixels, and the vertical axis also represents pixels.

[0091] Taking a fifth-order polynomial as an example, the calculation formula for Z′(x,y) is as follows:

[0092] Z′(x,y)=p00+p10*x+p01*y+p20*x 2 +p11*x*y+p02*y^2+p30*x^3+p21*x^2*y+p12*x*y^2+p03*y^3+p40*x^4+p31*x^3*y+p22*x^2*y ^2+p13*x*y^3+p04*y^4+p50*x^5+p41*x^4*y+p32*x^3*y^2+p23*x^2*y^3+p14*x*y^4+p05*y^5

[0093] The meaning of each symbol in the formula is:

[0094] p00, p10, p01, p20, p11, p02, p30, p21, p12, p03, p40, p31, p22, p13, p04, p50, p41, p32, p23, p14, and p05 are all coefficients of the polynomial.

[0095] Step 2, coarse correction, Z direction correction: Compensate the original data using formula (2) to complete the coarse correction. The formula for coarse correction is as follows:

[0096] D′(x,y,z)=D(x,y,z)*exp(-1j*Z′(x,y)*2*pi) (2)

[0097] The meanings of the symbols in formula (2) are:

[0098] D′(x,y,z) represents the OCT signal after coarse correction;

[0099] D(x,y,z) represents the raw OCT signal before coarse correction at depth z;

[0100] j represents the imaginary unit, that is,

[0101] Z′(x,y) represents the curved surface of the plane mirror in the Z direction;

[0102] Pi stands for pi, which is a constant;

[0103] x, y, z represent the three coordinate axes in a three-dimensional rectangular coordinate system;

[0104] The roughly rectified image is Figure 6 As shown, it can be seen that Figure 2 The saddle-shaped shadow in the middle has been well removed, restoring the original appearance of the plane mirror sample.

[0105] Step 3: Fine correction, correction in the XY plane: Fine correction is performed on the OCT signal after the coarse correction, that is, phase compensation is performed on the aberration in the XY plane. In other words, the aberration in the XY plane can be achieved by phase compensation. The formula for fine correction is shown in formula (3), which is as follows:

[0106] D″(x,y,z0)=D′(x,y,z0)*exp(-1j*Ф(x,y,z0)) (3)

[0107] The meanings of the symbols in formula (3) are:

[0108] D″(x,y,z0) represents the OCT signal after fine correction;

[0109] D′(x,y,z0) represents the OCT signal after coarse correction at depth z0;

[0110] Ф(x,y,z0) represents the compensation phase at depth z0;

[0111] z0 represents a certain depth on the Z axis;

[0112] Formula (3) can complete the aberration compensation of the OCT signal at the depth z0.

[0113] Step 4: Obtain corrected data: To compensate for different depths, place the plane mirror at the corresponding depth z0 and repeat steps 2 and 3 to complete the aberration compensation of the OCT signal in the entire xyz space.

[0114] It should be noted that the unit of depth is pixels, so the difference between two adjacent depths is 1 pixel. By moving the relative position of the plane mirror and the reference arm, the plane mirror can be placed at different depths. Generally speaking, there are many z0s that need to be measured, depending on the number of sampling points in the depth direction, such as 512, 1024, etc. However, in this type of OCT system, the compensation phases measured at different depths are extremely close, so the compensation phase at a certain depth can be used to replace the compensation phases at all depths, so only one measurement is required.

[0115] To obtain the compensation phase Φ(x, y, z0) in formula (3), the Zernike polynomial fitting method can be used. The expression of the Zernike polynomial is as follows:

[0116]

[0117] The meanings of the symbols in formula (4) and formula (5) are as follows:

[0118] Z represents the (m,n)th term of the Zernike polynomial;

[0119] represents an intermediate variable;

[0120] r is the coordinate axis of polar coordinates;

[0121] θ is the coordinate axis of polar coordinates;

[0122] m is the coefficient of the Zernike polynomial, m>0;

[0123] n is the coefficient of the Zernike polynomial, n>0, and n≥m;

[0124] k represents an integer between 0 and (nm) / 2; n=0, 1, 2, ... (nn) even; wherein even represents an even number.

[0125] The first 36 terms of the Zernike polynomials are distributed on the XY plane as follows: Figure 7 shown. Figure 7 It is the image of different terms of Zernike polynomials. The compensation process is the process of weighted superposition of these 36 terms.

[0126] The compensation phase Φ(x,y) can be viewed as a linear combination of these 36 polynomials. The calculation formula for the compensation phase Φ(x,y) is as follows:

[0127]

[0128] The meanings of the symbols in formula (6) are:

[0129] p k Indicates the linear superposition coefficient corresponding to the compensation phase in the XY plane;

[0130] S k represents the k-th Zernike polynomial.

[0131] It should be noted that, for a certain depth, the compensation phase of the corresponding plane can be represented by Ф(x, y); for the entire three-dimensional cube, a certain depth z0 can be represented by Ф(x, y, z0).

[0132] The entire correction process is shown in the attached Figure 9 , where the update coefficient p k The method can adopt optimization methods such as simulated annealing or gradient descent, that is, p k The updating method of is simulated annealing or gradient descent, etc., which will not be described here.

[0133] The specific steps for solving the Zernike polynomial coefficients are as follows:

[0134] Step 1, coefficient initialization: for example, p1 = 1, p2 = 1, p3 = 1...p36 = 1;

[0135] Step 2: Update the coefficients: For example, p1 = 1.1, p2 = 1.3, p3 = 1.2 ... p36 = 0.9. Step 3: Calculate the compensation phase from the coefficients: Substitute the compensation phase into formula (6);

[0136] Step 4, Phase Compensation: The calculation formula for phase compensation is as follows:

[0137] D_new=D_old*exp(j*phi(x,y) (7)

[0138] The meanings of the symbols in formula (7) are:

[0139] D_new represents the data after compensation;

[0140] D_old represents the data before compensation;

[0141] j represents the imaginary unit, that is,

[0142] Step 5: Calculate the mean square error of the compensated form: MSE(fft2(D_new)), where MSE is the mean square error and fft2 is the two-dimensional Fourier transform;

[0143] The calculation formula of MSE is as follows:

[0144]

[0145] The meanings of the symbols in formula (8) are:

[0146] N represents the total number of pixels on the xy plane;

[0147] y i (x, y) represents the actual observation value;

[0148] represents the predicted value obtained by the regression model;

[0149] For a plane mirror, the MSE of an ideal image without aberration should be 0.

[0150] Step 6: Determine whether the mean square error is less than the threshold. The threshold T can be set according to the actual situation:

[0151] When the mean square error is judged to be less than the threshold, the correction ends;

[0152] When it is determined that the mean square error is not less than the threshold, return to step 2.

[0153] After completing the above coarse correction and fine correction, the compensated surface shape Z′(x, y) and the compensated phase Φ(x, y, z0) at different depths can be obtained respectively. When measuring other samples, the aberration correction of the system can be completed using formula (2) and formula (3) respectively to obtain the original appearance of the sample.

[0154] The above description is only a preferred specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any technician familiar with the technical field, within the technical scope disclosed by the present invention, who makes equivalent replacements or changes based on the technical solution and inventive concept of the present invention, should be covered by the scope of protection of the present invention.

Claims

1. A method for correcting the aberration of an OCT galvanometer scanning image, characterized in that: The pixel displacement vector caused by the aberration in three-dimensional space is decomposed into the sum of the displacement vector in the Z direction and the displacement vector in the XY plane. The aberration is successfully compensated and the spatial correction of the entire three-dimensional image is completed. The following steps are included: Step 1: Obtain raw data: Place a plane mirror sample near the focal plane and use an OCT device to image the plane mirror sample to obtain a curved image in three-dimensional space. The distribution of the plane mirror surface in three-dimensional space is calculated based on the light intensity distribution. The formula is as follows: Z(x,y)=g(maxIdx(V(x,y,z),z)) (1) The meanings of the symbols in formula (1) are: Z(x,y) represents the offset of each point in the (x,y) plane in the Z direction; g represents median filtering; maxIdx(V(x,y,z),z) indicates the position of the maximum grayscale value in the Z direction; By fitting Z(x,y) with a polynomial, we can obtain the curved surface Z′(x,y) of the plane mirror in the Z direction. Step 2, coarse correction, Z direction correction: Compensate the original data to complete the coarse correction. The formula for coarse correction is as follows: D′(x,y,z)=D(x,y,z)*exp(-1j*z′(x,y)*2*pi) (2) The meanings of the symbols in formula (2) are: D′(x,y,z) represents the OCT signal after coarse correction; D(x,y,z) represents the raw OCT signal before coarse correction at depth z; j represents the imaginary unit; Z′(x,y) represents the curved surface of the plane mirror in the Z direction; pi stands for pi; x, y, z represent the three coordinate axes in a three-dimensional rectangular coordinate system; Step 3: Fine correction, correction in the XY plane: Perform fine correction on the coarsely corrected OCT signal, that is, perform phase compensation for the aberration in the XY plane. The formula for fine correction is as follows: D″(x,y,z0)=D′(x,y,z0)*exp(-1j*Ф(x,y,z0)) (3) The meanings of the symbols in formula (3) are: D″(x,y,z0) represents the OCT signal after fine correction; D′(x,y,z0) represents the OCT signal after coarse correction at depth z0; Ф(x,y,z0) represents the compensation phase at depth z0; z0 represents a certain depth on the Z axis; Step 4: Obtain corrected data: To compensate for different depths, place the plane mirror at the corresponding depth z0 and repeat steps 2 and 3 to complete the aberration compensation of the OCT signal in the entire xyz space.

2. The method for correcting the aberration of an OCT galvanometer scanning image according to claim 1, wherein: The calculation formula of the compensation phase Ф(x,y) is as follows: The meanings of the symbols in formula (6) are: p k Indicates the linear superposition coefficient corresponding to the compensation phase in the XY plane; S k represents the k-th Zernike polynomial.

3. The method for correcting the aberration of an OCT galvanometer scanning image according to claim 2, wherein: The p k The updating method is simulated annealing.

4. The method for correcting the aberration of an OCT galvanometer scanning image according to claim 2, wherein: The p k The update method is gradient descent.

5. The method for correcting the aberration of an OCT galvanometer scanning image according to claim 2, wherein: The specific steps of solving the Zernike polynomial coefficients are as follows: Step 1, coefficient initialization; Step 2: Update the coefficients; Step 3: Calculate the compensation phase from the coefficient; Step 4: Phase compensation; Step 5: Calculate the mean square error after compensation; Step 6: Determine whether the mean square error is less than the threshold: When the mean square error is judged to be less than the threshold, the correction ends; When it is determined that the mean square error is not less than the threshold, return to step 2.

6. The method for correcting the aberration of an OCT galvanometer scanning image according to claim 5, wherein: In step 4, the calculation formula for phase compensation is as follows: D_new=D_old*exp(j*phi(x,y)) (7) The meanings of the symbols in formula (7) are: D_new represents the data after compensation; D_old represents the data before compensation; j represents the imaginary unit.

7. The method for correcting the aberration of an OCT galvanometer scanning image according to claim 1, wherein: The plane mirror is a plane reflecting mirror.

8. The method for correcting the aberration of an OCT galvanometer scanning image according to claim 1, wherein: In the fourth step, the difference between two adjacent depths is 1 pixel.

9. The method for correcting the aberration of an OCT galvanometer scanning image according to claim 5, wherein: In the fifth step, the calculation formula of the mean square error MSE is as follows: The meanings of the symbols in formula (8) are: N represents the total number of pixels on the xy plane; y i (x, y) represents the actual observation value; Represents the predicted value obtained by the regression model.

Images