Image distortion real-time correction and real-time registration method based on high-speed laser scanning

Abstract

The invention discloses an image distortion real-time correction and real-time registration method based on high-speed laser scanning, and relates to the field of ophthalmology medical imaging and image processing. The method is suitable for various laser scanning systems based on the resonance galvanometer, and scanning lag and image distortion caused by sine scanning characteristics of the resonance galvanometer can be effectively corrected. According to the strip-based image real-time registration method, large-range full-frame motion and small-range intra-frame motion in the imaging process are accurately corrected, strip registration is carried out in the horizontal direction and the vertical direction, intra-frame gaps are avoided, and the continuity of the interior of the image is ensured. And meanwhile, a rolling average method is adopted to display a registration result in real time, so that the signal-to-noise ratio of a single-frame image is remarkably improved, and powerful support is provided for accurately observing a retina structure, dynamic and static characteristics of cells and hemodynamic changes.

Description

Technical Field

[0001] The present invention relates to the field of ophthalmic medical image processing and is applicable to various types of high-speed laser scanning imaging systems. Specifically, the present invention proposes an accurate correction method for the scanning lag and image distortion caused by the nonlinear scanning curve of the resonant galvanometer, and develops a strip-based real-time image registration algorithm, which effectively removes motion artifacts and significantly improves the signal-to-noise ratio of the real-time output image, laying the foundation for further quantitative analysis of weak signals at the cellular level. Background Art

[0002] As one of the main sensory organs of organisms, the eye is an important way to obtain information from the outside world. As the core component of the visual system, the retina is formed by a variety of optic nerve cells arranged according to specific rules. These optic nerve cells help the brain form visual perception by converting light signals into bioelectric signals and transmitting them to the visual center. However, with the emergence of various diseases, the normal physiological structure and function of the retina will be affected, and these diseases are often difficult to detect in the early stages. Therefore, it is very important to judge and prevent the type of disease by observing the morphology of the retina.

[0003] Scanning Laser Ophthalmoscopy (SLO) is a non-contact fundus imaging technology that usually uses a dual galvanometer galvanometer to perform two-dimensional point scanning imaging on the retina. Although the linear scanning characteristics of the galvanometer galvanometer enable it to obtain retinal images, the imaging speed is slow, which is not conducive to real-time observation of dynamic changes in the fundus. For example, with an image resolution of 250×250 pixels and a sampling rate of 250kHz, the frame rate may be limited to 2Hz due to the retracement time of the galvanometer. In order to achieve faster imaging speed and higher sampling density, researchers usually choose a resonant galvanometer instead of a galvanometer in the fast axis direction for scanning. For example, an 8kHz resonant galvanometer and a sampling rate of 25MHz can achieve a frame rate of 16Hz at an image resolution of 1562×1024 pixels. However, the scanning mode of the resonant galvanometer is sinusoidal motion. When the sampling rate is fixed, the image will produce a distortion of "wide in the middle and narrow on both sides" during the scanning process. In addition, due to the limitations of the drive circuit accuracy and mechanical processing technology, the resonant galvanometer will inevitably produce scanning lag in the periodic motion, that is, there is a certain time difference when the galvanometer passes the same position during the forward and reverse scanning process. In order to obtain a true and accurate image, the scanning lag and image distortion need to be effectively corrected.

[0004] In non-contact retinal imaging technology, eye movements that cannot be subjectively controlled (such as tremors, drifts, and tiny eye movements) have a significant impact on image quality. These movements can cause changes in the position, angle, and focal length of the image, thereby generating motion artifacts and reducing image quality. In human eye experiments, the fixation time of the eyeball usually does not exceed 3 seconds. In small animal in vivo experiments, although anesthesia can reduce uncontrollable body movements, eye displacement caused by breathing is still inevitable. In addition, the low signal-to-noise ratio of a single-frame image in the experiment will hinder researchers from observing the area of ​​interest, thereby reducing experimental efficiency. Therefore, it is urgent to develop a real-time registration method to correct eye movements, eliminate motion artifacts, and improve the quality of single-frame images to ensure the accuracy and reliability of experimental results. Summary of the invention

[0005] In view of the problems existing in the prior art, the present invention proposes a method for image distortion correction and real-time registration. This method is applicable to various laser scanning systems based on resonant galvanometers, and can effectively correct the scanning lag and image distortion caused by the sinusoidal scanning characteristics of the resonant galvanometer. The strip-based real-time image registration method accurately corrects the large-scale full-frame motion and small-scale intra-frame motion during the imaging process. By performing strip registration in the horizontal and vertical directions respectively, the appearance of intra-frame gaps is avoided, ensuring the continuity within the image. At the same time, the rolling average method is used to display the registration results in real time, which significantly improves the signal-to-noise ratio of a single-frame image, providing strong support for the accurate observation of retinal structure, dynamic and static characteristics of cells, and hemodynamic changes.

[0006] The technical solution of the present invention:

[0007] A real-time correction and registration method for image distortion based on high-speed laser scanning, the steps are as follows:

[0008] Step 1: Correct the lag during scanning

[0009] (1-1) The standard grid plate is imaged. The odd-numbered rows of pixels in the image are collected by forward scanning of the resonant galvanometer, and the even-numbered rows of pixels are collected by reverse scanning of the resonant galvanometer. The average intensity of the odd and even rows of the image is calculated respectively, and each average intensity is binarized to obtain the intensity-pixel position curves corresponding to the odd and even rows.

[0010] (1-2) The discrete points representing the boundaries of a single grid are extracted from the intensity-pixel position curve, and a cubic polynomial is fitted to the points to obtain the actual scanning displacement-pixel position curves of the galvanometer in the forward and reverse directions, respectively. The fitting formula is:

[0011] H θ (x) = θ 1 x+θ 2 x 2 +θ3 x 3

[0012] Among them, x is the pixel position corresponding to the discrete point, H θ (x) is the actual displacement of the galvanometer corresponding to the discrete point, θ 1 ,θ 2 ,θ 3 They represent the coefficients of the linear, quadratic, and cubic terms in the fitting curve respectively.

[0013] (1-3) The actual scanning displacement-pixel position curves in the forward and reverse scanning processes are fitted respectively to obtain a set of hysteresis curves. By calculating the energy consumption area of ​​the hysteresis curve and combining it with the actual displacement, the average number of pixel misalignments due to scanning lag can be estimated. The energy consumption area calculation formula is:

[0014]

[0015] Among them, A(x) is the energy consumption area, H θ+ (x), H θ- (x) are the actual scanning displacement-pixel position curves in the forward and reverse directions, respectively.

[0016] Calculate the average pixel misalignment number ΔP when the galvanometer moves to the same position based on the energy consumption area avg , the calculation formula is as follows:

[0017]

[0018] (1-4) The row pixels corresponding to the forward scanning of the galvanometer in the fixed image are calculated according to the average pixel misalignment number ΔP avg Adjust the corresponding row pixel position in the reverse scan to correct the image misalignment caused by scanning lag.

[0019] Step 2: Correct image distortion caused by scanning curve

[0020] (2-1) For the image after lag correction, recalculate the average intensity of the entire frame and perform binarization to obtain a single-line intensity-pixel position curve. Use the method in step (1-2) to fit the scanning curve again, and by solving the nonlinear equation, obtain the pixel position x of the linear model in the uniform scanning mode. l The corresponding pixel position x in the actual scanning curve r .

[0021]

[0022] H T (x l )=H θI (x r ), x r ∈[0,w-1]

[0023] Among them, H T (x l ) is the linear model in the uniform scanning mode, H θI (x r ) is the actual scanning displacement-pixel position curve, and w is the number of pixels contained in each row of the image.

[0024] (2-2) Based on the correspondence between the linear model calculated in step (2-1) and the pixel position in the actual scanning curve, the image after lag correction is interpolated using the inverse distance weighted interpolation algorithm. The specific formula is as follows:

[0025]

[0026]

[0027] y={0,1,2…h-1}

[0028] in, Indicates rounding down. Indicates rounding up, ω floor (x l ) and ω ceil (x l ) respectively represent the image after lag correction. and The weight at the position, I corr (x l , y) is the image after distortion correction, and h is the number of pixels in each column of the image. Step 3: Gaussian filling of the image boundary

[0029] (3-1) The image I after lag correction and distortion correction is corr (x l , y) expands to all sides, and initializes the expanded area with boundary pixels to obtain a new image I e (x l , y), the specific operations are as follows:

[0030]

[0031] Among them, δ is the number of expanded pixels.

[0032] (3-2) In order to make the pixel intensity of the extended area drop steadily and approach 0 at the boundary, the weights of each part can be generated with the help of Gaussian function, and a smooth transition can be achieved by adjusting the size of σ. The specific formula is as follows:

[0033]

[0034] (3-3) The one-dimensional Gaussian function can be used as the weight of the expanded rows and columns at the boundary of the original image. The formula is as follows:

[0035]

[0036] g T (k) = g(k), g L (k) = g T (k) T , at this time k={0,1,2…δ-1}

[0037] g B (k) = g(k+δ+1), g R (k) = g B (k) T , at this time k={0,1,2…δ-1}

[0038] Among them, σ is the standard deviation, which can control the distribution state of the Gaussian function. T (k), g B (k), g L (k), g R (k) represents the weights of the extended areas at the upper, lower, left, and right boundaries, respectively.

[0039] The two-dimensional Gaussian function can be used as the weight of the expanded area at the four corners of the original image. The formula is as follows:

[0040] G(k x , k y )=exp[-(1 / 2){((k x -δ) 2 +(k y -δ) 2 ) / σ 2}]k x , k y ={0, 1, 2…δ-1}

[0041] G TL (k x , k y )=G(k x , k y )

[0042] G TR (k x , k y )=G(k x +δ+1,k y )

[0043] G BL (k x , k y )=G(k x , ky +δ+1)

[0044] G BR (k x , k y )=G(k x +δ+1,k y +δ+1)

[0045] Among them, G TL (k x , k y ), G TR (k x , k y ), G BL (k x , k y ), G BR (k x , k y ) represent the weights of the expanded areas at the upper left, upper right, lower left, and lower right corners of the expanded image.

[0046] (3-4) By multiplying the weights with the pixel values ​​of the corresponding areas of the image, we can obtain the Gaussian-filled image and reduce the discontinuity of the boundaries.

[0047] The images mentioned in the subsequent steps must first be processed in steps one to three.

[0048] Step 4: Perform phase correlation calculation on the image after Gaussian filling

[0049] (4-1) Before performing phase correlation calculation, a reference frame needs to be selected as the benchmark for the entire registration sequence. The first frame image acquired after the registration starts is selected as the reference image, and subsequent images are registered as target images. The specific formula is as follows:

[0050] r(x, y) = t(x-Δx, y-Δy)

[0051] R(u,v)=T(u,v)e -2π(uΔx+vΔy)i

[0052]

[0053]

[0054] Among them, r(x, y) and R(u, v) represent the reference frames in the time domain and frequency domain respectively, t(x, y) and T(u, v) represent the target frames in the time domain and frequency domain respectively, and the cross power spectrum H(u, v) is calculated and inverse FFT is performed on it. The impulse function δ is obtained, and the target frame displacement (Δx, Δy) can be obtained by determining the coordinate position of the impulse function.

[0055] (4-2) By rearranging the pixel positions in the image, moving the four corner elements to the center, and combining Gaussian filtering to effectively remove false peaks, retain a single true peak, and further optimize the registration results. The template formula for rearranging pixel positions and calculating normalized Gaussian filtering is as follows:

[0056] I rearrange (x, y) = fftshift( I e(x,y))

[0057]

[0058] Among them, I e (x, y) is the image after Gaussian filling, I rearrange (x, y) is the rearranged image, W Nor (x, y) is the normalized Gaussian template, sum(w(x, y)) is the sum of the values ​​at each position of the Gaussian template, n is the length and width of the template, u x 、u y is the template center value, and σ is the standard deviation.

[0059] Step 5: Rough whole frame registration

[0060] After Gaussian filtering, usually only one clear peak value (true peak) is retained. According to the position of the true peak, the actual displacement of the image can be calculated and corrected accordingly. The specific actual displacement calculation formula is as follows:

[0061] x offset =Δx-W mid

[0062] y offset =Δy-H mid

[0063] Among them, x offset ,y offset The distribution is the actual displacement distance in the horizontal and vertical directions, Δx and Δy are the maximum peak positions detected after Gaussian filtering, and W mid , H mid are the midpoints of the image width and height, respectively.

[0064] Step 6: Fine strip registration

[0065] (6-1) The whole frame of the registered image is divided into strips of equal size in the horizontal and vertical directions.

[0066] (6-2) Repeat steps 4 to 6, use horizontal and vertical strips instead of the entire frame image for re-registration, and calculate the peak position in each strip to obtain the displacement between strips. The specific formula is as follows:

[0067] x s_offset =Δx hs -W mid

[0068] y s_offset =Δy vs -H mid

[0069] Among them, x s_offset ,y s_offset are the horizontal and vertical displacements of the strip, Δx hs is the horizontal coordinate of the maximum peak position in the horizontal strip, Δy vs The ordinate of the maximum peak position in the vertical strip. The calculated strip displacement is applied to the horizontal or vertical strips of the target frame to correct the slight displacement within the frame.

[0070] The beneficial effects of the present invention are as follows: the scanning lag and image distortion caused by the sinusoidal scanning characteristics of the resonant galvanometer are effectively corrected by using cubic polynomial fitting and inverse distance weighted algorithm interpolation. The registration of the whole frame and the strip is achieved through multiple operations such as Gaussian filling, phase correlation calculation, pixel rearrangement, Gaussian filtering, strip segmentation and rolling average, which can effectively correct the large-scale movement of the image while carefully adjusting the internal micro-displacement. The registration of the horizontal and vertical strips can ensure the continuity and accuracy of the frame content and avoid the appearance of gaps within the frame. The real-time display of the rolling average result significantly improves the signal-to-noise ratio of the single-frame image, enhances the visualization and stability of the image, and helps researchers accurately locate the area of ​​interest.

[0071] In addition, the present invention can be applied to post-processing algorithms to generate high-quality image results by quickly registering tens of thousands of frames of images, laying the foundation for subsequent analysis of tiny signals. BRIEF DESCRIPTION OF THE DRAWINGS

[0072] Figure 1 It is a flow chart of the present invention.

[0073] Figure 2 Schematic diagram of hysteresis correction. (a) is the binary curve of image intensity-pixel position of odd rows (forward scanning) and even rows (reverse scanning). Red represents odd rows (forward scanning) and blue represents even rows (reverse scanning). (b) is the actual displacement-pixel position curve during the scanning process of the resonant galvanometer. The blue line represents the forward scan; the green line represents the reverse scan; the red dotted line represents the fitting curve after misalignment correction; the yellow dotted line represents the linear model under the uniform scanning mode; the dots indicate the discrete points corresponding to each type of curve. The enlarged area can more clearly observe the hysteresis phenomenon.

[0074] Figure 3The distortion correction comparison diagram is shown in Figure 2. (a) The left side is the image before distortion correction, and the right side is the binary curve of image intensity-pixel position. (b) The left side is the image after distortion correction, and the right side is the binary curve of image intensity-pixel position.

[0075] Figure 4 Schematic diagram of Gaussian filling. (a) is the theoretical diagram of Gaussian filling. (b) is the effect diagram of Gaussian filling.

[0076] Figure 5 Schematic diagram of the change of image peak value during the registration process. (a) The effect after phase correlation calculation. (b) The effect after Gaussian filling. (c) The effect after pixel reorganization and Gaussian filtering.

[0077] Figure 6 Schematic diagram of the image registration process. (a) Rough whole-frame registration diagram. (b) Schematic diagram of horizontal segmentation strips. (c) Schematic diagram of vertical segmentation strips.

[0078] Figure 7 Comparison of human (reflected light) and mouse (fluorescent) retina before and after registration. (a) Original image of human eye. (b) Image after rolling average of human eye (16 frames). (c) Average result after human eye registration (100 frames). (d) Original image of mouse retina. (e) Image after rolling average of mouse retina (16 frames). (f) Average result after mouse retina registration (1000 frames). DETAILED DESCRIPTION

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

[0080] Example: Image distortion correction and real-time registration display during human eye reflected light imaging and mouse fluorescence imaging.

[0081] The present invention proposes a technical solution for real-time correction and registration to address the problems of image distortion and motion artifacts that occur during retinal imaging of humans and mice. Compared with other fundus imaging methods, SLO has multiple imaging modes that can achieve reflected light and fluorescence imaging. However, the experimental requirements for fluorescence imaging of the human eye are high and accompanied by certain risks. It must be completed under the accompaniment of medical staff. As a key biological model for basic scientific research and clinical ophthalmology research, mice have a pupil numerical aperture that is more than twice that of the human eye, have higher optical performance, and are more likely to observe dynamic changes in cells after fluorescent labeling. Therefore, mice have become the best choice for fluorescence imaging experiments.

[0082] In the specific experiment, the present invention uses human eyes as the experimental object for reflected light imaging, and mice as the experimental object for fluorescent imaging. The core algorithm is implemented by C++, and the user interface (GUI) is built in combination with Qt Creator. The SLO system acquires images at an acquisition rate of 16Hz. In order to ensure the consistency of pixel resolution in the horizontal and vertical directions of the image, the original image of 1562×1024pixels is interpolated to 1024×1024pixels. When storing data, the system sets corresponding storage intervals according to different experimental objects. The gaze duration of the human eye is usually less than 3s, so 100 frames are selected as a storage unit; for mouse experiments conducted under gas anesthesia, 1000 frames are selected as a storage unit.

[0083] Step 1: Correct the lag during scanning

[0084] (1-1) Image a standard black and white grid plate, where the spacing of the white grid is 0.25 mm and the spacing of the black grid is 0.05 mm. In order to calculate the degree of scanning lag, it is necessary to extract discrete points representing the actual displacement and pixel position in the odd rows (forward scanning) and even rows (reverse scanning) of the image respectively to fit the actual scanning curve. However, the interference of various types of noise during the data acquisition process will make it impossible to accurately distinguish the black and white boundaries of the grid in the image. Therefore, it is necessary to average the intensity of each pixel and perform binarization processing for subsequent calculations, such as Figure 2 As shown in (a) in .

[0085] (1-2) 40 discrete points were extracted in the forward and reverse scanning cycles for fitting the curves, e.g. Figure 2 As shown in (b), the forward and reverse scanning curves are:

[0086] H θ+ (x) = (8.153 × 10 -5 )x+(7.198×10 -6 )x 2 -(3.08×10 -9 )x 3

[0087] H θ- (x) = (2.805 × 10 -4 )x+(7.05×10 -6 )x 2 -(3.065×10 -9 )x 3

[0088] Among them, H θ+ (x) is the actual displacement-pixel position curve of forward scanning, H θ-(x) is the actual displacement-pixel position curve of reverse scanning, and x is the pixel position.

[0089] (1-3) The imaging area contains 19 black grids and 20 white grids. The maximum displacement of the scanning mirror is 5.95 mm, and the actual number of pixels in the image is 1562 pixels. The energy consumption area A(x) of the hysteresis curve is calculated as follows:

[0090] A(x)=∫ 0 1562 |H θ+ (x)-H θ- (x)|dx

[0091] Where A(x) = 77.0399 mm 2 , the average number of pixel misalignments is estimated based on the area, and the formula is as follows:

[0092]

[0093] The final calculated average number of pixel misalignments ΔP avg ≈13 pixels.

[0094] (1-4) Fix the pixel positions of odd rows in the image and calculate the ΔP obtained in step (1-3) avg Adjusting the pixel positions of even-numbered rows effectively corrects image misalignment caused by scanning lag.

[0095] Step 2: Correct image distortion caused by scanning curve

[0096] (2-1) The scanning hysteresis has been corrected in step 1, so there is no need to calculate the odd and even rows of the image separately. The intensity of the whole frame image that has completed the hysteresis correction is averaged in the vertical direction, and the single-line intensity-pixel position curve is obtained by binarization. 40 discrete points are extracted and steps (1-2) are repeated to fit the actual displacement curve:

[0097] H θ (x) = (8.186 × 10 -5 )x+(7.196×10 -6 )x 2 -(3.078×10 -9 )x 3

[0098] By solving the equation, the pixel position x of the linear model in the uniform scanning mode can be obtained l The actual scanning curve pixel position x r The corresponding relationship between the linear model H in the uniform scanning mode T (x) is as follows:

[0099]

[0100] H T (x i )=H θI (x r ), x r ∈[0, 1561]

[0101] Among them, H T (x i ) is the linear model in the uniform scanning mode, H θI (x r ) is the actual scanning displacement-pixel position curve.

[0102] (2-2) Based on the pixel position correspondence calculated in step (2-1), the inverse distance weighted interpolation algorithm is used to interpolate and correct the image after lag correction. The formula is as follows:

[0103]

[0104]

[0105]

[0106] y={0,1,2…1023}

[0107] in, Indicates rounding down. Indicates rounding up, ω floor (x i ) and ω ceii (x i ) respectively represent the image after lag correction. and The weight at the position, I corr (x i , y) is the image after distortion correction.

[0108] To ensure the consistency of the pixel resolution in the horizontal and vertical directions of the image, the above interpolation algorithm is used again to interpolate the original image of 1562×1024 pixels to 1024×1024 pixels. Figure 3 As shown in FIG. 4 , by comparing the images before and after distortion correction, it can be found that the method is effective.

[0109] Step 3: Gaussian filling of image boundaries

[0110] (3-1) Figure 4 As shown in (a), the image I after lag and distortion correction is completed corr (x i, y) expands to all sides, and initializes the expanded area with boundary pixels to obtain a new image I e (x i , y), the specific operations are as follows:

[0111]

[0112] Among them, δ=5 is the number of pixels expanded, w=1024 pixels is the image width, and h=960 pixels is the image height. In the image registration process, there is no need to register the retrace area, so the image height involved in the subsequent processing will be subtracted from the height of the retrace area.

[0113] (3-2) The weights of each part are generated by using the Gaussian function, and σ = 2.15 is set to make the intensity of the edge pixels of the image decrease smoothly:

[0114]

[0115] (3-3) Calculate the one-dimensional and two-dimensional Gaussian functions as the weights of the image expansion area. The one-dimensional Gaussian function is used as the weight of the expanded rows and columns at the boundary of the original image. The formula is as follows:

[0116]

[0117] gT(k)=g(k),g L (k) = g T (k) T , at this time k={0,1,2…δ-1}

[0118] g B (k) = g(k+δ+1), g R (k) = g B (k) T , at this time k={0,1,2…δ-1}

[0119] Among them, g T (k), g B (k), g L (k), g R (k) represents the weights of the extended areas at the upper, lower, left, and right boundaries, respectively.

[0120] The two-dimensional Gaussian function is used as the weight of the expanded area at the four corners of the original image. The formula is as follows:

[0121] G(k x , k y )=exp[-(1 / 2){((k x -δ) 2 +(k y -δ) 2 ) / σ2}]k x , k y ={0, 1, 2…δ-1}

[0122] G TL (k x , k y )=G(k x , k y )

[0123] G TR (k x , k y )=G(k x +δ+1,k y )

[0124] G BL (k x , k y )=G(k x , k y +δ+1)

[0125] G BR (k x , k y )=G(k x +δ+1,k y +δ+1)

[0126] Among them, G TL (k x , k y ), G TR (k x , k y ), G BL (k x , k y ), G BR (k x , k y ) represent the weights of the expanded areas at the upper left, upper right, lower left, and lower right corners of the expanded image respectively.

[0127] (3-4) Multiply the calculated weight by the corresponding extended area in the image, such as Figure 4 As shown in (b), a Gaussian filled image with smoothly decreasing edge pixel intensity is obtained.

[0128] Step 4: Perform phase correlation calculation on the image after Gaussian filling

[0129] (4-1) Since the present invention involves real-time registration, the first frame image acquired after the registration is started is selected as the reference image, and the subsequent frames are selected as the target frames to perform phase correlation calculation:

[0130] r(x, y) = t(x-Δx, y-Δy)

[0131] R(u,v)=T(u,v)e -2π(uΔx+v)i

[0132]

[0133] Among them, r(x, y) and R(u, v) are the reference frames in the time domain and frequency domain, t(x, y) and T(u, v) are the target frames in the time domain and frequency domain respectively. The cross power spectrum H(u, v) is calculated in the frequency domain and the inverse FFT transform is performed on it. The impulse function δ is obtained, and the coordinate position of the impulse function is the displacement (Δx, Δy) of the target frame.

[0134] (4-2) Various noises, spectrum leakage, edge discontinuities and other realistic factors in the image may cause multiple peaks in the impulse function, resulting in incorrect registration results. Although Gaussian filling of the image boundary can effectively improve this problem, false peaks may still exist at the four corners of the image. In order to improve the robustness and accuracy of the registration method, the image I after Gaussian filling is e The resolution of (x, y) is increased to 1034×970 pixels, and the coordinates of the center of the image are (518, 486). The image is divided into four parts along the center line, and the pixel positions are rearranged to obtain a new image I rearrange (x, y), the formula is as follows:

[0135] I rearrange (x, y) = fftshift(I e (x, y)

[0136] When performing fftshift processing, image I e The upper left and lower right areas, and the lower left and upper right areas of (x, y) are swapped to move the four corners of the image to the center. Subsequently, the rearranged image is low-pass filtered using a Gaussian template with a standard deviation of σ = 3 and a template size of n = 13 × 13 to minimize the false peak value. The normalized Gaussian template calculation formula is as follows:

[0137]

[0138] Among them, w(x, y) is the normalized Gaussian filter template, sum(w(x, y)) is the sum of the values ​​at each position of the Gaussian template, and u x ,u y The center position of the template.

[0139] like Figure 5 As shown in (a) in the figure, an image that has only been calculated through phase correlation will be affected by various noises, spectrum leakage, and edge discontinuities, which will cause the false peaks in the image to cover up the real peaks. Figure 5 As shown in (b) in Figure 3, after Gaussian filling, although the false peaks in the center of the image are effectively suppressed, there are still false peaks in the four corners. Figure 5 As shown in (c) in the figure, the pixel positions are rearranged and Gaussian filtering is applied to effectively remove the false peaks, retain the single true peak, and significantly improve the stability of the results.

[0140] Step 5: Rough whole frame registration

[0141] The actual displacement of the entire frame is calculated based on the peak position to correct a wide range of eye movements, such as Figure 6 The specific displacement calculation formula is as follows:

[0142] x offset =Δx-W mid

[0143] y offset =Δy-H mid

[0144] Among them, x offset ,y offset The distribution is the actual displacement distance in the horizontal and vertical directions, Δx and Δy are the maximum peak positions detected after Gaussian filtering, and W mid =518, H mid =486 are the midpoints of the image width and height respectively.

[0145] Step 6: Fine strip registration

[0146] (6-1) The whole frame of the registered image is divided into four strips of 1024×240 pixels in the horizontal direction, as follows: Figure 6 As shown in (b) of . It is divided into 4 strips of 256×960 pixels in the vertical direction, as shown in Figure 6 As shown in (c) in .

[0147] (6-2) Repeat steps 4 to 6, use horizontal and vertical strips instead of the entire frame image to re-register, and calculate the displacement between strips. The specific formula is as follows:

[0148] x s-offset =Δx hs -W mid

[0149] y s_offset =Δy vs -H mid

[0150] Among them, x s_offset ,y s_offset are the horizontal and vertical displacements of the strip, Δxhs is the horizontal coordinate of the maximum peak position of the horizontal strip, Δy vs is the ordinate of the maximum peak position of the vertical strip. The horizontal strip of the target frame is adjusted according to the calculated displacement to correct the micro-motion within the frame.

[0151] Step 7, result display

[0152] (7-1) Calculate the average result of the intensity of the recorded 100 frames of human eye reflected light images and 1000 frames of mouse retinal fluorescence images. The specific process is as follows: establish two arrays, one of which is used as a counting board to count the number of times the pixel intensity at each position of the image is not 0 after the time series registration is completed, and the other is used to accumulate the pixel intensity at the same position after the registration. The real average image is obtained by calculation, such as Figure 7 As shown in (c) and (f) in .

[0153] M(x, y) = S(x, y) / A(x, y)

[0154] Among them, M(x, y) is the registered average image, S(x, y) is the sum array, A(x, y) is the count array, if A(x, y) = 0, then M(x, y) = 0.

[0155] (7-2) Rolling average display: Calculate the rolling average result of the registered image and display it in real time at a frame rate of 16 Hz.

[0156]

[0157] Among them, M i (x,y),A i (x, y) is the i-th group rolling average result and count array, I sj (x, y) is the jth frame image after strip registration, and k=16 is the number of frames contained in each group. Figure 7 As shown in (a) and (b), (d) and (e), the real-time rolling average results significantly improve the signal-to-noise ratio of a single frame image, and intuitively display the structure of retinal blood vessels, the static distribution of cells, and their dynamic changes.

Claims

1. A real-time correction and registration method for image distortion based on high-speed laser scanning, characterized in that Its real-time nature allows it to be used for high-speed distortion correction and image registration during data acquisition. The steps are as follows: Step 1: Correct the lag during scanning (1-1) Imaging a standard grid plate, the odd-numbered rows of pixels of the image are collected by forward scanning of a resonant galvanometer, and the even-numbered rows of pixels are collected by reverse scanning of a resonant galvanometer, and the average intensity of the odd-numbered rows and the even-numbered rows of the image are calculated respectively, and each average intensity is binarized to obtain the intensity-pixel position curves corresponding to the odd-numbered rows and the even-numbered rows; (1-2) The discrete points representing the boundaries of a single grid are extracted from the intensity-pixel position curve, and a cubic polynomial is fitted to the points to obtain the actual scanning displacement-pixel position curves of the galvanometer in the forward and reverse directions. The fitting formula is: H θ (x)=θ1x+θ2x 2 +θ3x 3 Among them, x is the pixel position corresponding to the discrete point, H θ (x) is the actual displacement of the galvanometer corresponding to the discrete point, θ1, θ2, and θ3 represent the coefficients of the first-order, second-order, and third-order terms in the fitting curve, respectively; (1-3) The actual scanning displacement-pixel position curves in the forward and reverse scanning processes are fitted respectively to obtain a set of hysteresis curves. By calculating the energy consumption area of ​​the hysteresis curve and combining it with the actual displacement, the average number of pixel misalignments due to scanning lag can be estimated. The energy consumption area calculation formula is: Among them, A(x) is the energy consumption area, H θ+ (x), H θ- (x) are the actual scanning displacement-pixel position curves in the forward and reverse directions, respectively; Calculate the average pixel misalignment number ΔP when the galvanometer moves to the same position based on the energy consumption area avg , the calculation formula is as follows: (1-4) The row pixels corresponding to the forward scanning of the galvanometer in the fixed image are calculated according to the average pixel misalignment number ΔP avg Adjust the corresponding row pixel position in the reverse scan to correct the image misalignment caused by scanning lag; Step 2: Correct image distortion caused by scanning curve (2-1) For the image after lag correction, recalculate the average intensity of the entire frame image, and perform binarization processing to obtain a single-line intensity-pixel position curve; use the method of step (1-2) to fit the scanning curve again, and obtain the pixel position x of the linear model in the uniform scanning mode by solving the nonlinear equation. l The corresponding pixel position x in the actual scanning curve r ; H T (x l )=H θI (x r ),x r ∈[0,w-1] Among them, H T (x l ) is the linear model in the uniform scanning mode, H θI (x r ) is the actual scanning displacement-pixel position curve, w is the number of pixels contained in each row of the image; (2-2) Based on the correspondence between the linear model calculated in step (2-1) and the pixel position in the actual scanning curve, the image after lag correction is interpolated using the inverse distance weighted interpolation algorithm. The specific formula is as follows: y={0,1,2…h-1} in, Indicates rounding down. Indicates rounding up, ω floor (x l ) and ω ceil (x l ) respectively represent the image after lag correction. and The weight at the position, I corr (x l , y) is the image after distortion correction, h is the number of pixels contained in each column of the image; Step 3: Gaussian filling of image boundaries (3-1) The image I after lag correction and distortion correction is corr (x l , y) expands to all sides, and initializes the expanded area with boundary pixels to obtain a new image I e (x l , y), the specific operations are as follows: Among them, δ is the number of expanded pixels; (3-2) In order to make the pixel intensity of the extended area drop steadily and approach 0 at the boundary, the weights of each part can be generated with the help of Gaussian function, and a smooth transition can be achieved by adjusting the size of σ. The specific formula is as follows: (3-3) The one-dimensional Gaussian function can be used as the weight of the expanded rows and columns at the boundary of the original image. The formula is as follows: g T (k)=g(k),g L (k)=g T (k) T , evttak={0,1,2…δ-1} g B (k)=g(k+δ+1),g R (k)=g B (k) T , evttak={0,1,2…δ-1} Among them, σ is the standard deviation, which can control the distribution state of the Gaussian function; g T (k), g B (k), g L (k), g R (k) represents the weights of the extended areas at the upper, lower, left, and right boundaries, respectively; The two-dimensional Gaussian function can be used as the weight of the expanded area at the four corners of the original image. The formula is as follows: G(k x ,k y )=exp[-(1 / 2){((k x -δ) 2 +(k y -δ) 2 ) / σ 2 }] k x ,k y ={0,1,2…δ-1} G TL (k x ,k y )=G(k x ,k y ) G TR (k x ,k y )=G(k x +δ+1,k y ) G BL (k x ,k y )=G(k x ,k y +δ+1) G BR (k x ,k y )=G(k x +δ+1,k y +d+1) Among them, G TL (k x , k y ), G TR (k x , k y ), G BL (k x , k y ), G BR (k x , k y ) represent the weights of the expanded areas at the upper left, upper right, lower left, and lower right corners of the expanded image respectively; (3-4) Multiplying the weights by the pixel values ​​of the corresponding areas of the image will yield a Gaussian-filled image, thus reducing the discontinuity of the boundaries. The images mentioned in the subsequent steps must first be processed through steps one to three; Step 4: Perform phase correlation calculation on the image after Gaussian filling (4-1) Before performing phase correlation calculation, a reference frame needs to be selected as the benchmark for the entire registration sequence. The first frame image acquired after the registration starts is selected as the reference image, and subsequent images are registered as target images. The specific formula is as follows: r(x, y) = t(x-Δx, y-Δy) R(u,v)=T(u,v)e -2π(uΔx+vΔy)i Among them, r(x, y) and R(u, v) represent the reference frames in the time domain and frequency domain respectively, t(x, y) and T(u, v) represent the target frames in the time domain and frequency domain respectively, and the cross power spectrum H(u, v) is calculated and inverse FFT is performed on it. Obtain the impulse function δ, and by determining the coordinate position of the impulse function, the target frame displacement (Δx, Δy) can be obtained; (4-2) By rearranging the pixel positions in the image, moving the four corner elements to the center, and combining Gaussian filtering to effectively remove false peaks, retain a single true peak, and further optimize the registration results; the template formula for rearranging pixel positions and calculating normalized Gaussian filtering is as follows: I rearrange (x,y)=fftshift(I e (x,y)) Among them, I e (x, y) is the image after Gaussian filling, I rearrange (x, y) is the rearranged image, w Nor (x, y) is the normalized Gaussian template, sum(w(x, y)) is the sum of the values ​​at each position of the Gaussian template, n is the length and width of the template, u x 、u y is the template center value, σ is the standard deviation; Step 5: Rough whole frame registration After Gaussian filtering, usually only one clear true peak value is retained; according to the position of the true peak, the actual displacement of the image can be calculated and corrected accordingly; the specific actual displacement calculation formula is as follows: x offset =Δx-W mid y offset =Δy-H mid Among them, x offset ,y offset The distribution is the actual displacement distance in the horizontal and vertical directions, Δx and Δy are the maximum peak positions detected after Gaussian filtering, and W mid , H mid are the midpoints of the image width and height, respectively; Step 6: Fine strip registration (6-1) Divide the whole frame of the registered image into strips of equal size in the horizontal and vertical directions; (6-2) Repeat steps 4 to 6, use horizontal and vertical strips instead of the entire frame image for re-registration, and calculate the peak position in each strip to obtain the displacement between strips. The specific formula is as follows: x s_offset =Δx hs -W mid y s_offset =Δy vs -H mid Among them, x s_offset ,y s_offset are the horizontal and vertical displacements of the strip, Δx hs is the horizontal coordinate of the maximum peak position in the horizontal strip, Δy vs is the ordinate of the maximum peak position in the vertical strip; the calculated strip displacement is applied to the horizontal or vertical strips of the target frame to correct the slight displacement within the frame.

Images