Axial positioning error suppression method in three-dimensional microscopic trajectory tracking

Abstract

The invention discloses an axial positioning error suppression method in three-dimensional microscopic trajectory tracking. The invention relates to the technical field of particle positioning error suppression, in particular to an axial positioning error suppression method in three-dimensional microscopic trajectory tracking. The method comprises the following steps: acquiring a three-dimensional motion track of a particle; gaussian fitting in three orthogonal motion directions is carried out; axial positioning error suppression of Fourier transform high-order filtering is carried out; axial positioning error suppression is carried out on the basis of a Winsorizing method; carrying out axial positioning error suppression on Gaussian distribution statistical parameters; and obtaining a three-dimensional motion track of the particle target with the same random jitter degree in the horizontal direction and the axial direction. According to the method, the problem that in the three-dimensional tracking process of the particles, due to the fact that the axial positioning precision is different from the horizontal direction due to the limitation of a characterization mode, the three-dimensional motion track of the particles is distorted is solved. The method is especially suitable for a high-signal-to-noise-ratio system. The method is suitable for various three-dimensional particle tracking technologies based on out-of-focus imaging.

Description

Technical Field

[0001] The present invention relates to the technical field of particle positioning error suppression, and in particular to a method for suppressing axial positioning errors in three-dimensional microscopic trajectory tracking. Background Art

[0002] With the continuous expansion of the application scope of microorganisms and drug particles in the fields of medical treatment and industrial production, and the rapid development of related scientific research technologies, the microscopic characterization of particles, microbubbles or microorganisms is no longer limited to the two-dimensional level. Three-dimensional particle tracking analysis based on microscopic imaging is gradually becoming an important real-time microscopic characterization method.

[0003] The three-dimensional particle tracking analysis technologies that have been commercialized on the market include: fluorescence microscopy based on defocus imaging, digital holographic microscopy, total internal reflection microscopy, etc. These technologies mainly rely on the two-dimensional optical imaging sensor plane (such as CMOS or CCD sensors of scientific digital cameras, photomultiplier tube arrays, etc.) to obtain the three-dimensional real-time motion trajectory of the target sample. It is usually necessary to use different imaging and positioning principles in the axial direction than in the horizontal direction, which further leads to inconsistent positioning accuracy in the axial and horizontal directions. This inconsistency in positioning accuracy will not only cause trajectory distortion, but also may introduce significant differences in the data analysis of the axial and horizontal directions of the trajectory. For example, the analysis of mean square displacement (MSD) may be affected. Especially in samples with low exposure time, high frame rate and high signal-to-noise ratio conditions where the real signal is relatively weak, the differences and effects will be more obvious. Commonly used trajectory smoothing methods may also cause similar significant differences. At present, there is still a lack of effective solutions on the market to solve the trajectory distortion problem caused by inconsistent positioning errors in the axial and horizontal directions. Summary of the invention

[0004] The present invention proposes a method for suppressing axial positioning errors in three-dimensional microscopic trajectory tracking. The method is a three-dimensional motion trajectory data optimization method based on statistical laws. It utilizes the isotropy of Brownian motion of particles and the Gaussian distribution characteristics of instantaneous velocity, and combines Gaussian fitting, Fourier transform filtering, Winsorizing method and other non-smooth trajectory processing methods to suppress the trajectory distortion caused by the different positioning accuracy of the microscopic characterization device in the horizontal direction and the axial direction. This error suppression method can make the suppressed trajectory closer to the actual situation, and effectively avoid introducing significant differences in the horizontal direction and the axial direction in the data processing of the suppressed trajectory, and obtain the three-dimensional motion trajectory of the particle target with the same random jitter degree in the horizontal direction and the axial direction. The present invention solves the problem that in the process of three-dimensional tracking of particles, the axial positioning accuracy is different from that in the horizontal direction due to the limitation of the characterization method, thereby causing the three-dimensional motion trajectory of the particles to be distorted. It is particularly suitable for high signal-to-noise ratio systems caused by low exposure time, high frame rate and small samples. It is suitable for a variety of three-dimensional particle tracking technologies based on defocus imaging such as interference, fluorescence, dark field and digital holography.

[0005] The present invention is implemented through the following technical solutions:

[0006] A method for suppressing axial positioning errors in three-dimensional microscopic trajectory tracking. The method comprises the following steps:

[0007] S1. Obtaining the three-dimensional motion trajectory of the particle using a microscopic characterization device, including a three-dimensional peak search sub-step and a three-dimensional scattered point connection sub-step;

[0008] S2, performing positioning optimization on the peak search parameters of the three orthogonal motion directions of the three-dimensional motion trajectory based on Gaussian fitting in each direction, to obtain the three-dimensional motion trajectory after positioning optimization;

[0009] S3, axial positioning error suppression by Fourier transform high-order filtering of abnormal high-frequency unidirectional repetitive jitter in the three-dimensional motion trajectory;

[0010] S4, differentiating the three-dimensional motion trajectory to obtain the axial instantaneous velocity, performing a normality test on the distribution of the instantaneous velocity, and performing axial positioning error suppression based on the Winsorizing method on the three-dimensional motion trajectory;

[0011] S5. Perform Gaussian distribution fitting on the instantaneous velocities in three directions of the three-dimensional motion trajectory, correct the abnormal distribution caused by the positioning error according to the ratio of the standard deviation of the distribution in each direction, and perform axial positioning error suppression on the axial instantaneous velocity distribution using the statistical parameters of the Gaussian distribution;

[0012] S6. Obtain a three-dimensional motion trajectory of the particle target with the same degree of random jitter distribution for horizontal and axial trajectory positioning.

[0013] Furthermore, in the step S1, a microscopic characterization device is used to obtain the three-dimensional motion trajectory of the particle, and the microscopic characterization device has different axial positioning accuracy and horizontal positioning accuracy for the particle.

[0014] Furthermore, in step S2, the light intensity as the peak-seeking parameter is used to confirm the peak-seeking location by a Gaussian distribution fitting method, and the all-directional positioning is optimized. The fitting calculation formula is as follows:

[0015]

[0016] Among them, I is the peak-finding parameter within a certain range around the original three-dimensional motion trajectory, which is the light intensity here; I0 is the background noise level of the peak-finding parameter; z0 is the original three-dimensional peak-finding position; the parameter to be fitted A is the peak height parameter; the parameter to be fitted σ is the standard deviation of the peak-finding parameter, which reflects the distribution width of the peak-finding parameter; the parameter to be fitted z G is the three-dimensional position after Gaussian distribution fitting.

[0017] Furthermore, in step S3, Fourier transform high-order filtering is performed on the three-dimensional motion trajectory to be suppressed, and the calculation formula is as follows:

[0018]

[0019] Where z[n] is the original time domain input; N is the length of the input vector; Z[k] is the transformed frequency domain output; k is the index of the current frequency component, ranging from 1 to N; the inverse Fourier transform formula after filtering is as follows:

[0020]

[0021] Where Z[k] is the frequency domain input; N is the length of the input vector; z[n] is the time domain output; and n is the index of the current time component, ranging from 1 to N.

[0022] Furthermore, each three-dimensional motion trajectory should include three-dimensional coordinate information and defocus distance of a certain particle (including solid particles, micro-nano bubbles, micro-nano vacuoles or microorganisms) at several time nodes within a certain time period.

[0023] Furthermore, the maximum defocus distance of the three-dimensional motion trajectory should be within the depth of field allowed by the microscopic characterization device.

[0024] Furthermore, the positioning accuracy in the horizontal direction and the axial direction of the microscopic characterization device for obtaining the three-dimensional motion trajectory of the particle can be considered not to change with time in a single or a single series of three-dimensional motion trajectory characterization experiments.

[0025] Furthermore, in step S5, the method for suppressing the axial positioning error of the Gaussian distribution statistical parameters is as follows: for the three-dimensional motion trajectory processed by the above steps, the instantaneous speed in the horizontal direction and the axial direction is obtained, and the instantaneous speed in the three directions is fitted with Gaussian distribution, and the fitting formula is as follows:

[0026]

[0027] Among them, v and v0 are the instantaneous speed in one direction before and after fitting, P is the probability value of the distribution, P0 is the probability reference value, w is the overall standard deviation of the Gaussian distribution, and B is the proportional coefficient. v0, P0, w, and B are all parameters to be fitted.

[0028] Compared with the prior art, the advantages of the present invention are:

[0029] The present invention proposes for the first time a method for suppressing the distortion of three-dimensional motion trajectory caused by system errors of microscopic characterization devices. Compared with the trajectory smoothing method of averaging the three-dimensional motion trajectory within a range, this method focuses on screening out data that are greatly affected by system errors and correcting these data, which will better retain the true speed and acceleration characteristics of the original three-dimensional motion trajectory. It is suitable for systems with high signal-to-noise ratios due to low exposure time, high frame rate and small samples, and the measurement results are more susceptible to system errors. It is suitable for a variety of particle three-dimensional motion trajectory tracking technologies based on interference, fluorescence, dark-field digital holography, etc. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Figure 1 It is a flow chart of a method for suppressing axial positioning errors in three-dimensional microscopic trajectory tracking according to the present invention;

[0031] Figure 2 is a comparison diagram of the typical three-dimensional motion trajectory of the polystyrene standard microspheres dispersed in water in Example 1 before and after inhibition;

[0032] Figure 3 This is a comparison chart of the MSD results of the typical trajectory of the polystyrene standard microspheres dispersed in water in Example 1 before and after suppression.

[0033] Figure 4 This is a comparison diagram of the three-dimensional motion trajectory of the typical trajectory of the polystyrene standard microspheres dispersed in the liquid crystal in Example 2 before and after inhibition. DETAILED DESCRIPTION

[0034] The following is an example of a method for suppressing the distortion of a three-dimensional motion trajectory caused by a system error of a characterization device proposed by the present invention, with reference to the accompanying drawings. However, the embodiments of the present invention are not limited thereto.

[0035] Example 1

[0036] A particle axial positioning optimization method based on error suppression, such as Figure 1 As shown, the following steps are included.

[0037] S1. Obtain the three-dimensional motion trajectory of the sample. Prepare microparticles, microbubbles or microorganism samples and observe them under a microscopic characterization device that can characterize their three-dimensional motion trajectory. Use an optical sensor to record a series of images of the microparticles in real time, containing three-dimensional motion trajectory information (the information here refers to kinematic information such as velocity, acceleration, mean square end displacement, etc.). The sample in this embodiment uses polystyrene standard microspheres with an average particle size of 0.96 μm dispersed in ultrapure water, and the concentration of the microspheres is 10 7 / mL. The three-dimensional microscopic characterization device uses a digital holographic microscope; the light source is an LED light source with a central wavelength of 365nm, a 100x objective lens, and the camera parameters are set as follows: the frame rate is set to 40fps, the exposure time is 0.025s, and the image resolution is 1024×1024; 2400 optical interference pattern grayscale images with a bit depth of 16bit are continuously shot, and the side length of the square field of view of each optical interference pattern is 66.56μm, and the side length of a single pixel after magnification is 65nm. The three-dimensional peak-finding sub-step and the three-dimensional scattered point connection sub-step are as follows: through a program based on the principle of light propagation, the three-dimensional light field of polystyrene standard microspheres with an average particle size of 0.96 μm in each frame is calculated and restored, and the peak is found in the three-dimensional light field by the maximum light intensity criterion to obtain the preliminary three-dimensional position information of the polystyrene standard microspheres; according to the principle of space-time consistency, the microsphere coordinate position information obtained in each frame and expressed in data form as three-dimensional scattered points is connected in series into a three-dimensional motion trajectory; in this embodiment, the following axial positioning error suppression will revolve around these trajectories.

[0038] S2. Perform directional positioning optimization based on Gaussian fitting on the peak-seeking parameter distribution of the three orthogonal motion directions of the three-dimensional motion trajectory. According to the guidance of the three-dimensional motion trajectory, the light intensity information of the three-dimensional light field within the radius of 32 pixels and the surrounding three-dimensional motion trajectory is obtained by using the original camera image and the three-dimensional light field reconstruction program based on the light propagation principle. The light intensity information is the peak-seeking parameter; the distribution curve of the light intensity is obtained in the X and Y horizontal directions and the axial direction with the three-dimensional coordinates given by the trajectory as the center; the peak-seeking positioning position is confirmed by the Gaussian distribution fitting method for the light intensity as the peak-seeking parameter, and the directional positioning is optimized. The distribution curve of the light intensity is fitted based on the least squares method using the following Gaussian distribution formula:

[0039]

[0040] Among them, I is the peak-finding parameter within a certain range around the original trajectory, which is the light intensity of the three-dimensional light field here; I0 is the background noise level of the three-dimensional light intensity, which is obtained by averaging the darker 80% light intensity values ​​in the entire three-dimensional light field; z0 is the three-dimensional coordinate position given by the trajectory, that is, the original three-dimensional peak-finding position, and the axial position is taken as an example here; the parameter to be fitted A is the peak height parameter; the parameter to be fitted σ is the standard deviation of the peak-finding parameter, which reflects the distribution width of the peak-finding parameter; the parameter to be fitted z G is the axial position of the microsphere in this frame after Gaussian distribution fitting. The steps for positioning optimization of Gaussian fitting in the two horizontal directions of X and Y are the same as above. After fitting the distribution curve of light intensity in the horizontal direction and the axial direction as described above, the center peak coordinates are obtained, which are the three-dimensional coordinates of the microparticle after positioning optimization. Repeat the above operation for each frame of the three-dimensional motion trajectory to obtain the three-dimensional motion trajectory of the microparticle after positioning optimization.

[0041] S3. Suppress the axial positioning error by Fourier transform high-order filtering of abnormal high-frequency unidirectional repetitive jitter in the three-dimensional motion trajectory. Since most particles, microbubbles or microorganisms will not perform long-term high-frequency jitter movement in a certain direction without external conditions, the same applies to polystyrene standard microspheres without external conditions. The occurrence of this kind of high-frequency jitter can be determined as regular abnormal axial positioning. Therefore, for the axial positioning result with the largest degree of jitter in the three-dimensional motion trajectory, the axial instantaneous velocity information is obtained by differentiating the axial position and performing a one-dimensional Fourier transform. The three-dimensional motion trajectory is processed by Fourier transform high-order filtering, and the formula used is as follows:

[0042]

[0043] Where z[n] is the original time domain input, which is the axial position at time n after Gaussian fitting positioning optimization; N is the length of the input vector. Since the difference method is used for calculation speed, the vector length here is 2399; Z[k] is the frequency domain output after Fourier transformation; k is the index of the current frequency component, ranging from 1 to 2399, and i is the imaginary unit;

[0044] For high-frequency data that reaches 50% of the highest frequency after transformation, the abnormal high-frequency peaks are eliminated. For most particles, regular high-frequency jitter is abnormal, and high-frequency data is generally close to 0. Here, the high-frequency data that is greater than 2 times the mean of all high-frequency data is eliminated; then a one-dimensional inverse Fourier transform is performed to obtain the axial instantaneous velocity of the Fourier transform high-order filter. The formula for the inverse Fourier transform after filtering is as follows:

[0045]

[0046] Among them, Z[k] is the frequency domain input, k is the frequency domain value; N is the length of the input vector, here the vector length is 2399; z'[n] is the time domain output; n is the index of the current time component, ranging from 1 to 2399, and i is the imaginary unit. After the transformation, the axial motion trajectory and the three-dimensional motion trajectory can be further obtained according to the following steps.

[0047] S4. Perform axial positioning error suppression based on the Winsorizing method. The Winsorizing method is a method that makes a set of data that does not follow the Gaussian distribution more consistent with the Gaussian distribution by correcting extreme values. Differentiate the three-dimensional motion trajectory to obtain the axial instantaneous velocity, extract the axial instantaneous velocity information of each frame of the axial data of the three-dimensional motion trajectory processed in the previous step, and perform a normality test on the distribution of these instantaneous velocities; if the result of the normality test is less than 0.05 (indicating that the current data is not significantly from the normal distribution population), further determine the 5% and 95% quantiles of the current instantaneous velocity distribution, and perform Winsorizing data processing. Winsorizing data processing is a signal processing method that clips outliers. It reduces the impact of possible false outliers by limiting the extreme values ​​in the statistical data, thereby transforming the statistical data; obtain the instantaneous velocity after axial positioning error suppression based on the Winsorizing method, integrate the velocity to obtain the axial motion trajectory, and then combine it with the horizontal motion trajectory to obtain the three-dimensional motion trajectory.

[0048] S5. Perform axial positioning error suppression using Gaussian distribution statistical parameters on the distribution of instantaneous velocities in the two horizontal directions and the axial direction. For the three-dimensional motion trajectory processed by the above steps, obtain the instantaneous velocities in the horizontal direction and the axial direction, and perform Gaussian distribution fitting on the instantaneous velocities in the three directions. The fitting formula is as follows:

[0049]

[0050] Among them, v and v0 are the instantaneous speed in one direction before and after fitting, P is the probability value of the distribution, P0 is the probability reference value, w is the overall standard deviation of the Gaussian distribution, and B is the proportional coefficient. v0, P0, w, and B are all parameters to be fitted. Here we mainly focus on the fitting parameter w.

[0051] The fitting is performed to obtain three groups of Gaussian distribution parameters in three directions, each group of parameters includes the mean and standard deviation of the instantaneous velocity distribution in that direction; if the standard deviation of the axial direction is still significantly different from the two standard deviations in the horizontal direction, that is, the means are significantly different at the 0.05 level in the variance analysis, then the standard deviation of the axial velocity is set equal to the average of the two standard deviations in the horizontal direction; in this way, the instantaneous velocity of each axial direction will change, and the changed axial instantaneous velocity will be inverted through the definite integral to obtain the axial trajectory after the axial positioning error suppression of the statistical parameters of the Gaussian distribution, and combined with the horizontal trajectory as a reference for positioning error suppression, the three-dimensional motion trajectory is obtained.

[0052] S6. Obtain the three-dimensional motion trajectory of the particle target corresponding to the original trajectory, with the same random jitter degree in the horizontal and axial directions and positioning error suppression. Select one of the typical trajectories, and the comparison diagram between this trajectory and the original trajectory is as follows: Figure 2 By comparison, it can be seen that the axial abnormal positioning jump in the original trajectory is effectively suppressed, and the processed trajectory is more consistent with the classic Brownian motion trajectory. Figure 3 The MSD results of the trajectory before processing show an obvious upward warping of the axial direction compared to the horizontal direction in the short time interval part. In the MSD results of the trajectory after processing, this upward warping is suppressed, which is more consistent with the property of Brownian motion being the same in all directions.

[0053] Example 2

[0054] The sample in Example 1 was changed to a liquid liquid crystal dispersed with polystyrene standard microspheres with an average particle size of 0.2 μm and a viscosity of 40 mPa·s. The concentration of the microspheres was 10 8 / mL.

[0055] The light source is an LED light source with a central wavelength of 505nm, a 100x objective lens, and the camera parameters are set as follows: the frame rate is set to 40fps, the exposure time is 0.025s, the image resolution is 1024×1024, and 12,000 optical interference pattern grayscale images with a bit depth of 16bit are continuously captured. The side length of the square field of view of each optical interference original image is 66.56μm, and the side length of the unit pixel after magnification is 65nm. The position information of the polystyrene standard microspheres in space is calculated and restored through a program based on the principle of light propagation, and the three-dimensional motion trajectory of the Brownian motion of the microspheres in the liquid crystal is obtained.

[0056] In the axial positioning error suppression of the three-dimensional motion trajectory, compared with Example 1, in step S5, the standard deviation of the axial velocity is not directly equal to the average of the two standard deviations in the horizontal direction, but the axial velocity is multiplied by a parameter, which comes from the ratio of the standard deviation of the axial velocity of the polystyrene standard microspheres to the two standard deviations in the horizontal direction under the same experimental conditions and defocus distance. The remaining experimental and data processing steps and the formulas and parameters used remain unchanged.

[0057] After the above image capture and data processing process, the three-dimensional motion trajectory of the particle corresponding to the original trajectory is obtained, with the positioning error suppressed and the same random jitter degree in the horizontal and axial directions. Select one of the typical trajectories, and the comparison chart between this trajectory and the original trajectory is as follows Figure 4 By comparison, it can be seen that the axial abnormal positioning jump in the original trajectory is effectively suppressed, and the processed trajectory can still reflect the anisotropic properties of the liquid crystal medium.

[0058] The above embodiments are preferred implementation modes of the present invention, but the implementation modes of the present invention are not limited to the above embodiments. Any other changes, modifications, substitutions, combinations, and simplifications that do not deviate from the spirit and principles of the present invention should be equivalent replacement methods and are included in the protection scope of the present invention.

Claims

1. A method for suppressing axial positioning errors in three-dimensional microscopic trajectory tracking, characterized in that: The following steps are involved: S1, obtaining the three-dimensional motion trajectory of the particle; S2, perform Gaussian fitting in three orthogonal motion directions; S3, performing Fourier transform high-order filtering to suppress axial positioning errors; S4, perform normality test and axial positioning error suppression based on Winsorizing method; S5, performing axial positioning error suppression using Gaussian distribution statistical parameters; S6. Obtain the three-dimensional motion trajectory of the particle target with the same random jitter degree in the horizontal direction and the axial direction.

2. The method for suppressing axial positioning errors in three-dimensional microscopic trajectory tracking according to claim 1, characterized in that: The microparticles include solid particles, micro-nano bubbles, micro-nano vacuoles and microorganisms.

3. The method for suppressing axial positioning errors in three-dimensional microscopic trajectory tracking according to claim 1, characterized in that: In step S1, obtaining the three-dimensional motion trajectory of the particle consists of a three-dimensional peak search sub-step and a three-dimensional scattered point connection sub-step.

4. The method for suppressing axial positioning errors in three-dimensional microscopic trajectory tracking according to claim 1, characterized in that: In step S2, the peak search parameters are fitted based on Gaussian distribution to confirm the peak search location and optimize the positioning in all directions. The fitting calculation formula is as follows: Where I is the peak-finding parameter within a certain range around the original trajectory; I0 is the background noise level of the peak-finding parameter; z0 is the original three-dimensional peak-finding position; the parameter to be fitted A is the peak height parameter; the parameter to be fitted σ is the standard deviation of the peak-finding parameter, which reflects the distribution width of the peak-finding parameter; the parameter to be fitted z G is the three-dimensional position after Gaussian distribution fitting.

5. The method for suppressing axial positioning errors in three-dimensional microscopic trajectory tracking according to claim 1, characterized in that: In step S3, the three-dimensional motion trajectory to be suppressed is subjected to Fourier transform high-order filtering, and the calculation formula is as follows: Where z[n] is the original time domain input; N is the length of the input vector; Z[k] is the transformed frequency domain output; k is the index of the current frequency component, ranging from 1 to N; the inverse Fourier transform formula after filtering is as follows: Where Z[k] is the frequency domain input; N is the length of the input vector; z[n] is the time domain output; and n is the index of the current time component, ranging from 1 to N.

6. The method for suppressing axial positioning errors in three-dimensional microscopic trajectory tracking according to claim 1, characterized in that: The Winsorizing method described in step S4 is a signal processing method for clipping outliers, which reduces the impact of possible false outliers by limiting the extreme values ​​in the statistical data, thereby transforming the statistical data.

7. The method for suppressing axial positioning errors in three-dimensional microscopic trajectory tracking according to claim 1, characterized in that: In step S5, the distribution of the instantaneous velocities in the two horizontal directions and the axial direction is subjected to axial positioning error suppression using Gaussian distribution statistical parameters. For the three-dimensional motion trajectory processed by the above steps, the instantaneous velocities in the horizontal direction and the axial direction are obtained, and Gaussian distribution fitting is performed on the instantaneous velocities in the three directions. The fitting formula is as follows: Among them, v and v0 are the instantaneous speeds in one direction before and after fitting, P is the probability value of the distribution, P0 is the probability reference value, w is the overall standard deviation of the Gaussian distribution, B is the proportional coefficient, and v0, P0, w, and B are all parameters to be fitted.

Images