CN113671682B - Frequency domain light source position accurate correction method based on Fourier laminated microscopic imaging - Google Patents

Abstract

In the frequency-domain light source position accurate correction method based on Fourier stack microscopic imaging of the present invention, (1) aiming at the problem of inaccurate correction position, an independent frequency-domain position error model is introduced for each LED in the LED array for precise correction ; (2) Aiming at the problem that the correction process needs to be involved in the image reconstruction process and the time consumption is too long, a high-resolution intensity image collected by a high-magnification objective lens is introduced as a reference image to simulate the Fourier stack microscopic imaging technology. During the image acquisition process, a virtual low-resolution intensity image corresponding to each LED is generated. Combining the virtual and actual collected images, a loss function is established, and the optimal frequency-domain position parameters corresponding to each LED are searched by a two-dimensional particle swarm algorithm. Therefore, the present invention uses the high-resolution intensity image under the high-magnification objective lens and the frequency-domain position error model to realize the accurate correction of the frequency-domain position of the LED light source based on the principle of Fourier stack imaging, and has the advantages of simple implementation, high robustness, and accurate position Etc.

Description

A method for accurate correction of light source position in frequency domain based on Fourier stack microscopy

Technical Field

The present invention relates to the technical field of quantitative phase computational microscopy, and in particular to a method for accurately correcting the position of a frequency domain light source based on Fourier stack microscopy.

Background Art

Traditional optical microscopy technology forms image contrast by recording light intensity information, opening the door to the microscopic world for humans. However, for some transparent biological samples, traditional optical microscopy technology that only relies on light intensity measurement values for imaging seems powerless. In 1930, Frits Zernike invented a high-contrast imaging technology for phase objects. Since then, quantitative phase imaging has gradually become a major method for label-free imaging of transparent biological samples. Unlike traditional optical microscopy, quantitative phase imaging forms high-contrast quantitative images based on the optical path delay introduced by the sample, and provides objective phase measurements of target features. Since the phase is proportional to the product of the refractive index and thickness of the sample, these measurements can further reveal the deep structural characteristics of the sample and provide corresponding conditions for subsequent judgment and decision-making. At present, quantitative phase imaging technology has been widely studied and applied in biological observation, cell detection, drug screening, trace element research, precision materials and other fields.

In order to break some inherent limitations in traditional optical microscopy technology, such as the mutual constraints between system resolution and objective field of view, system diffraction limit conditions, etc., quantitative phase computational microscopy technology came into being. As an emerging quantitative phase computational microscopy technology, Fourier ptychographic microscopy (FPM) can simultaneously achieve wide field of view, high resolution, multi-modality, digital refocusing, aberration correction, quantitative phase imaging and other functions by effectively combining synthetic aperture, phase retrieval, optimization theory and other computational imaging related technologies. Compared with traditional digital optical microscopy systems, this technology only needs to use a programmable LED array to replace the traditional reflective or active light source to complete the modification of the system hardware part, which is cheap and easy to implement.

The implementation of Fourier stacking microscopy technology is divided into two processes: image acquisition process and image reconstruction process. For the image acquisition process, the LEDs in the array are lit one by one under a low-power objective lens as the illumination source. When a single LED at different positions in the array is lit to obliquely illuminate a two-dimensional thin sample, the spectrum equivalent to the sample undergoes a relative translation effect in the frequency domain. This position translation transformation introduced by different oblique illuminations enables the objective lens to collect sample spectrum information outside the original cutoff frequency and capture a series of low-resolution intensity images on the camera surface. Later, during the image reconstruction process, these intensity images will iteratively update the sample's spectrum information in the frequency domain and eventually synthesize into a new spectrum. After inverse Fourier transform, the reconstructed spectrum information can be converted into high-resolution sample complex amplitude information. Similar to the concept of synthetic aperture, Fourier stacking microscopy technology improves the equivalent cutoff frequency of the system, thereby simultaneously achieving wide field of view, high resolution and quantitative phase imaging under low-power microscopes.

As the light source in the Fourier stacking microscope, the position parameters of each LED in the LED array determine the update position of each sub-spectrum in the frequency domain, and its accuracy directly affects the quality of the final image reconstruction. Therefore, in the actual system, accurate correction of each LED position is an important part of realizing Fourier stacking microscopy imaging technology. In 2016, Sun et al. first established a spatial global error model to describe the position error of the LED array, and combined the simulated annealing algorithm with the image reconstruction process to complete the position correction of the LED. The relevant experimental results show that this method effectively and automatically corrects the LED position error and reduces the demand for light source position accuracy in the construction of the experimental system. In 2017, Pan et al. combined the embedded pupil function recovery algorithm (EPRY) on this basis to further improve the robustness of the correction algorithm. However, these algorithms are all based on a unified spatial correction model, and the error between each LED and the actual position is characterized by a consistency parameter. In addition, the simulated annealing algorithm used is very easy to fall into a local optimal solution rather than a global optimal solution. Based on the above two analysis points, the existing correction algorithms still have the risk of inaccurate correction. On the other hand, the correction process based on the simulated annealing algorithm must participate in the image reconstruction process, which makes the entire correction process appear very lengthy.

Summary of the invention

The present invention aims to provide a method for accurately correcting the position of a frequency domain light source based on Fourier stack microscopy, so as to solve the problems of inaccurate correction position and excessively long correction process time in the existing correction algorithm based on Fourier stack microscopy technology.

The present invention provides a method for accurately correcting the position of a frequency domain light source based on Fourier stack microscopy imaging, comprising the following steps:

Step S100, FPM acquisition and reconstruction with position error under low-power lens: LEDs in the programmable LED array are illuminated one by one under the low-power objective lens, and a series of low-resolution intensity images of the sample at the corresponding positions of different LEDs are captured on the camera surface; then the captured low-resolution intensity images are reconstructed in the frequency domain using a reconstruction algorithm to obtain reconstructed high-resolution intensity and phase complex amplitude with position error;

Step S200, high-resolution intensity image acquisition and regional registration under a high-power objective lens: after replacing the low-power objective lens with a high-power objective lens and completing the focusing operation, a high-resolution intensity image of the sample under the high-power objective lens is acquired; then, the low-resolution intensity image and the high-resolution intensity image are matched using an image registration algorithm and the regions representing the same information in the two images are framed, and the regions representing the same information in the high-resolution intensity image are interpolated to obtain a registered high-resolution intensity image;

Step S300, generating a virtual low-resolution intensity image of the corresponding position of the LED: using the high-resolution intensity with position error and the phase complex amplitude and the registered high-resolution intensity image to form a new spectrum diagram of the sample; according to the ideal frequency domain position of each LED in the programmable LED array, obtaining the sub-spectrum information corresponding to the ideal frequency domain position of each LED from the new spectrum diagram, and obtaining a series of virtual low-resolution intensity images by inverse Fourier transforming the sub-spectrum information;

Step S400, searching for position parameters of independent LEDs in the array: based on the fact that the translation effect of the corresponding position of the LED in the spatial domain is equivalent to the translation effect of the aperture position in the frequency domain, a frequency domain position error model is introduced for each LED; thus, the ideal frequency domain position of the current LED is taken as the search center, the search area is specified, and the optimal position parameter search is performed using a two-dimensional particle swarm algorithm; when the difference between the generated virtual low-resolution intensity image and the currently actually collected low-resolution intensity image is the smallest, this position is considered to be the actual frequency domain position corresponding to the LED, and the position correction of the LED is completed.

Furthermore, the method for collecting FPM with position error under a low-power microscope in step S100 includes the following sub-steps:

Step S101, using a programmable LED array with a specification of m×n as a light source; assuming that the interval between adjacent LEDs in the programmable LED array is d, the ideal position P _i,j of the LED in the i-th row and j-th column (1≤i≤m, 1≤j≤n) is expressed as:

P _i,j = (x _i,j ,y _i,j ) = (i·d,j·d)

Where x _i,j and y _i,j represent the two-dimensional coordinate component values of the horizontal plane where the LED is located in the airspace;

Step S102, assuming that the light source is far enough from the sample, the illumination light wave emitted from the LED can be regarded as a monochromatic plane wave. Therefore, when a certain LED is used to illuminate the sample, the corresponding wave vector k _i,j is expressed as:

与

分别表示波矢量k_i,j在x与y方向上的波矢分量值；λ表示LED辐射的中心波长值；h表示样本距离当前照明样本的LED的垂直距离；In the formula,

and

Respectively represent the wave vector component values of wave vector k _i,j in the x and y directions; λ represents the central wavelength value of LED radiation; h represents the vertical distance between the sample and the LED currently illuminating the sample;

Step S103, normalize the intensity value of the incident light field of the sample to 1, then the spectrum emitted from the sample is expressed as:

F{e(r)}＝F{o(r)exp(ik _i,j r)}＝O(kk _i,j )

Wherein, F{} represents the two-dimensional Fourier transform, e(r) represents the complex amplitude emitted from the sample, o(r) represents the complex amplitude modulation function of the sample space, that is, the complex amplitude function of the sample; o(k) represents the spectrum of the sample in the frequency domain; O(kk _i,j ) represents that the spectrum emitted by the sample is equivalent to the spectrum center of the sample being translated to the wave vector k _i,j ;

Step S104, the complex amplitude emitted from the sample is low-pass filtered in the frequency domain by the pupil function P(k) of the low-power objective lens. The low-pass filtering process is expressed as:

Gi _,j (k)＝O(kki _,j )P(k)

Where G _i,j (k) represents the sample spectrum after passing through the low-power objective lens;

表示为：Step S105: the sample spectrum is propagated to the image plane for imaging and is captured by the camera and converted into a low-resolution intensity image.

It is expressed as:

Where, g _i,j (r) represents the complex amplitude of the sample on the image plane, and F ^-1 {} represents the two-dimensional inverse Fourier transform;

Step S106, lighting up the LEDs in the programmable LED array one by one, and executing steps S102 to S105 for each lit LED, thereby capturing a series of low-resolution intensity images of the sample at positions corresponding to different LEDs.

Furthermore, the reconstruction method in step S100 includes the following sub-steps:

Step S111, spectrum initialization, that is, the result of interpolation of the low-resolution intensity image collected at the center position of the programmable LED array is used as the initial value of the intensity part of the reconstructed sample spectrum, and the zero phase is selected as the initial value of the relevant phase, and then combined with the pupil function of the low-power objective lens, the initial value of the reconstructed sample spectrum is solved as:

表示重建样本频谱的初始值；u与v表示二维频域中的坐标分量；R()表示相关插值过程，I_center(r)表示LED阵列中心位置的LED照明样本时采集的低分辨强度图像，k_center表示LED阵列中心位置的LED照明样本时对应的波矢量；

表示初始位相；In the formula,

represents the initial value of the reconstructed sample spectrum; u and v represent the coordinate components in the two-dimensional frequency domain; R() represents the related interpolation process, I _center (r) represents the low-resolution intensity image collected when the LED is illuminated at the center of the LED array, and k _center represents the wave vector corresponding to the LED illumination sample at the center of the LED array;

represents the initial phase;

Step S112, obtaining the sub-aperture low-resolution complex amplitude, that is, using the principle of translation effect to intercept the corresponding sub-aperture sample spectrum from the current sample spectrum, and then using two-dimensional inverse Fourier transform to obtain the sub-aperture low-resolution complex amplitude estimation value under the corresponding LED:

表示更新位置序列(i,j)对应子孔径的样本频谱；(u_i,j,v_i,j)表示位置序列(i,j)在频域中对样本的平移效应，P(u+u_i,j,v+v_i,j)表示平移效应结果；

表示第n次迭代对应位置序列(i,j)的子频谱复振幅估计值；Where n represents the current number of iterations, 1≤n≤n _max , and n _max represents the set maximum number of iterations.

represents the sample spectrum of the sub-aperture corresponding to the updated position sequence (i, j); (u _{i, j} , vi _{, j} ) represents the translation effect of the position sequence (i, j) on the sample in the frequency domain, and P(u+u _{i, j} ,v+vi _{, j} ) represents the translation effect result;

Represents the estimated value of the complex amplitude of the sub-spectrum corresponding to the position sequence (i, j) at the nth iteration;

Step S113, updating the sub-spectrum complex amplitude estimation value:

The amplitude part of the sub-spectrum complex amplitude estimate is replaced by the low-resolution intensity image actually collected at the position sequence (i, j), and the phase part is kept unchanged. The update formula of the sub-spectrum complex amplitude estimate is:

表示更新后的第n次迭代对应位置序列(i,j)的子频谱复振幅估计值；In the formula,

Represents the updated sub-spectrum complex amplitude estimate of the position sequence (i, j) corresponding to the nth iteration;

Step S114, using two-dimensional inverse Fourier transform to obtain the corresponding sub-aperture spectrum information for the sub-spectrum complex amplitude estimation value of the nth iteration corresponding position sequence (i, j):

表示对应

的子孔径频谱信息，F{}表示二维傅里叶变换；In the formula,

Indicates the corresponding

Sub-aperture spectrum information, F{} represents two-dimensional Fourier transform;

Afterwards, the overall sample spectrum is updated according to the sub-aperture positions corresponding to the position sequence (i, j):

表示更新后的样本频谱；In the formula,

represents the updated sample spectrum;

Step S115, updating the sub-aperture spectrum information corresponding to all position sequences {(i,j)1≤i≤m,1≤j≤n} according to steps S112 to S114, and obtaining the n-th iterative updated sample spectrum as the reconstruction result O ⁿ (u,v);

其中，每一次迭代循环的初始值为上一次循环更新后的样本频谱结果；Step S116, iterate the steps S112 to S115 until the reconstruction result converges or the loop ends, and obtain the final reconstruction result, that is, the reconstructed sample spectrum with position error, which is expressed as

The initial value of each iteration cycle is the sample spectrum result updated in the previous cycle;

并将样本高分辨复振幅

分为强度部分与位相部分：Step S117, the reconstructed sample spectrum with position error is subjected to a two-dimensional inverse Fourier transform to obtain a reconstructed sample high-resolution complex amplitude with position error.

The sample high-resolution complex amplitude

It is divided into intensity part and phase part:

表示重建的带位置误差的样本高分辨复振幅，

表示重建的带位置误差的样本高分辨复振幅中的强度部分，

表示重建的带位置误差的样本高分辨复振幅中的位相部分。In the formula,

represents the reconstructed sample high-resolution complex amplitude with position error,

represents the intensity part of the reconstructed sample high-resolution complex amplitude with position error,

Represents the phase part of the reconstructed sample high-resolution complex amplitude with position error.

Furthermore, the method for determining whether the reconstruction result converges in step S116 is:

The loss judgment function corresponding to the nth iteration cycle is:

In the formula, when ^En reaches the minimum value or ( ^En - ^En-1 )/En ^-1 is less than the set threshold, the reconstruction result is considered to have converged.

Furthermore, step S200 includes the following sub-steps:

Step S201, performing a coarse matching of the display area of the low-resolution intensity image and the high-resolution intensity image through a correlation matching operation to obtain a coarse matching area of the two images;

Step S202, performing a binarization operation on the rough matching area, extracting the coordinates of the top, bottom, left and right four points used to describe the boundary of the area, and the area selected by the four points is identified as the area representing the same information;

Step S203, based on the FPM reconstructed image resolution, the regions representing the same information in the high-resolution intensity image are interpolated to obtain a registered high-resolution intensity image I ^hr .

Furthermore, step S300 includes the following sub-steps:

位相部分保持不变，经二维傅里叶变换后获得虚拟频谱O^vir：Step S301, using the registered high-resolution acquisition image I ^hr to replace the intensity part of the high-resolution intensity and phase complex amplitude with position error

The phase part remains unchanged, and the virtual spectrum ^Ovir is obtained after two-dimensional Fourier transform:

Step S302, following the FPM acquisition process in step S100, sub-aperture sampling is performed according to the ideal frequency domain position (u+u _i,j ,v+vi _,j ) corresponding to the position sequence (i,j):

表示位置序列(i,j)虚拟的采样子孔径频谱；In the formula,

represents the virtual sampling sub-aperture spectrum of the position sequence (i, j);

Step S303: Perform a two-dimensional Fourier transform on the virtual sampling sub-aperture spectrum to obtain a corresponding virtual low-resolution intensity map.

Step S304, repeat steps S302 to S303 until all position sequence LEDs are sampled, generating a series of virtual low-resolution intensity images

Furthermore, step S400 includes the following sub-steps:

Step S401, introduce an independent frequency domain position error model (Δu _i,j , Δv _i,j ) for each position sequence (i,j). Then, the corresponding virtual low-resolution intensity image after the frequency domain position error model is introduced is expressed as:

Where Δu _i,j and Δv _i,j represent the position error components of the position sequence (i, j) in the two-dimensional frequency domain;

Step S402, using the two-dimensional correlation function Corr2 to evaluate the difference between the virtual low-resolution intensity image corresponding to the position sequence (i, j) and the low-resolution intensity image actually acquired in step S100:

表示位置序列(i,j)对应的虚拟低分辨强度图像的均值，

表示位置序列(i,j)对应的步骤S100实际采集的低分辨强度图像的均值；In the formula,

represents the mean of the virtual low-resolution intensity image corresponding to the position sequence (i, j),

represents the mean value of the low-resolution intensity image actually acquired in step S100 corresponding to the position sequence (i, j);

Step S403, using a two-dimensional particle swarm algorithm (PSO) to search for the best frequency domain position parameters of the position sequence (i, j) in the frequency domain:

(Δu _i,j ,Δv _i,j )=argmin(fitness)

fitness＝1-Corr2for(i,j)

Where fitness represents the fitness function describing the particle swarm, i.e., the loss function. When the fitness function value of all particles converges to the minimum value, it means that the difference between the generated virtual low-resolution intensity image and the currently actually collected low-resolution intensity image is the smallest. At this time, (Δu _i,j ,Δv _i,j ) is the frequency domain correction parameter of the position sequence (i,j). Then, the frequency domain position parameters of the position sequence (i,j) are corrected using the frequency domain correction parameters to obtain the corrected optimal frequency domain position parameters.

Furthermore, step S404 further includes:

Step S404, using the corrected optimal frequency domain position parameters, perform FPM acquisition and reconstruction again according to the method of step S100, that is, obtain the corrected and reconstructed high-resolution intensity and phase complex amplitude with position error.

In summary, due to the adoption of the above technical solution, the beneficial effects of the present invention are:

In the method for accurately correcting the position of a frequency domain light source based on Fourier stack microscopy imaging of the present invention, (1) in order to address the problem of inaccurate correction position, an independent frequency domain position error model is introduced for each LED in the LED array for accurate correction; (2) in order to address the problem of excessive time consumption caused by the need to participate in the image reconstruction process during the correction process, a high-resolution intensity image collected by a high-power objective lens is introduced as a reference image to simulate the image acquisition process of the Fourier stack microscopy imaging technology, and a virtual low-resolution intensity image corresponding to each LED is generated. Combining the virtual and actual collected images, a loss function is established, and a two-dimensional particle swarm algorithm is used to search for the optimal frequency domain position parameters corresponding to each LED. Therefore, the present invention utilizes the high-resolution intensity image under the high-power objective lens and the frequency domain position error model to realize accurate correction of the frequency domain position of an LED light source based on the principle of Fourier stack imaging, which has the advantages of simple implementation, high robustness, and accurate position.

BRIEF DESCRIPTION OF THE DRAWINGS

In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings in the embodiments will be briefly introduced below. It should be understood that the following drawings only show certain embodiments of the present invention and therefore should not be regarded as limiting the scope. For ordinary technicians in this field, other related drawings can be obtained based on these drawings without paying creative work.

FIG. 1 is a flow chart of a method for accurately correcting the position of a frequency domain light source based on Fourier stack microscopy according to an embodiment of the present invention.

DETAILED DESCRIPTION

In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Generally, the components of the embodiments of the present invention described and shown in the drawings here can be arranged and designed in various different configurations.

Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the invention claimed for protection, but merely represents selected embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention.

Example

As shown in FIG1 , this embodiment proposes a method for accurately correcting the position of a frequency domain light source based on Fourier stack microscopy imaging, comprising the following steps:

Step S100, FPM acquisition and reconstruction with position error under low-power objective: light up the LEDs in the programmable LED array one by one under the low-power objective, and capture a series of low-resolution intensity images of the sample at the corresponding positions of different LEDs on the camera surface; then use the reconstruction algorithm to complete the reconstruction of the captured low-resolution intensity image in the frequency domain to obtain a high-resolution intensity and phase image of the sample with position error;

Specifically:

(1) The method of FPM acquisition with position error under low-power microscope includes the following sub-steps:

Step S101, using a programmable LED array with a specification of m×n as a light source; assuming that the interval between adjacent LEDs in the programmable LED array is d, the ideal position P _i,j of the LED in the i-th row and j-th column (1≤i≤m, 1≤j≤n) is expressed as:

P _i,j = (x _i,j ,y _i,j ) = (i·d,j·d)

Where x _i,j and y _i,j represent the two-dimensional coordinate component values of the horizontal plane where the LED is located in the airspace;

Step S102, assuming that the light source is far enough from the sample, the illumination light wave emitted from the LED can be regarded as a monochromatic plane wave. Therefore, when a certain LED is used to illuminate the sample, the corresponding wave vector _{k i,j} is expressed as:

与

分别表示波矢量k_i,j在x与y方向上的波矢分量值；λ表示LED辐射的中心波长值；h表示样本距离当前照明样本的LED的垂直距离；In the formula,

and

Respectively represent the wave vector component values of wave vector k _i,j in the x and y directions; λ represents the central wavelength value of LED radiation; h represents the vertical distance between the sample and the LED currently illuminating the sample;

Step S103, normalize the intensity value of the incident light field of the sample to 1, then the spectrum emitted from the sample is expressed as:

F{e(r)}＝F{o(r)exp(ik _i,j r)}＝O(kk _i,j )

Wherein, F{} represents the two-dimensional Fourier transform, e(r) represents the complex amplitude emitted from the sample, o(r) represents the complex amplitude modulation function of the sample space, that is, the complex amplitude function of the sample; o(k) represents the spectrum of the sample in the frequency domain; O(kk _i,j ) represents that the spectrum emitted by the sample is equivalent to the spectrum center of the sample being translated to the wave vector k _i,j ;

Step S104, the complex amplitude emitted from the sample is low-pass filtered in the frequency domain by the pupil function P(k) of the low-power objective lens. The low-pass filtering process is expressed as:

Gi _,j (k)＝O(kki _,j )P(k)

Where G _i,j (k) represents the sample spectrum after passing through the low-power objective lens;

表示为：Step S105: the sample spectrum is propagated to the image plane for imaging and is captured by the camera and converted into a low-resolution intensity image.

It is expressed as:

Where, g _i,j (r) represents the complex amplitude of the sample on the image plane, and F ^-1 {} represents the two-dimensional inverse Fourier transform;

Step S106, lighting up the LEDs in the programmable LED array one by one, and executing steps S102 to S106 for each lit LED, thereby capturing a series of low-resolution intensity images of the sample at positions corresponding to different LEDs.

(2) The reconstruction method includes the following sub-steps:

Step S111, spectrum initialization:

In order to make the iterative solution algorithm have a faster convergence speed, it is necessary to first give the initial value of the solution to the reconstructed spectrum. Generally, normal incidence is selected, and the result of the interpolation of the low-resolution intensity image collected at the center of the programmable LED array is used as the initial value of the intensity part of the reconstructed sample spectrum, and the zero phase is selected as the initial value of the relevant phase. Combined with the pupil function of the low-magnification objective lens, the initial value of the solution to the reconstructed sample spectrum is:

表示重建样本频谱的初始值；u与v表示二维频域中的坐标分量；R()表示相关插值过程，I_center(r)表示LED阵列中心位置的LED照明样本时采集的低分辨强度图像，k_center表示LED阵列中心位置的LED照明样本时对应的波矢量；

表示初始位相，一般取0；In the formula,

represents the initial value of the reconstructed sample spectrum; u and v represent the coordinate components in the two-dimensional frequency domain; R() represents the related interpolation process, I _center (r) represents the low-resolution intensity image collected when the LED is illuminated at the center of the LED array, and k _center represents the wave vector corresponding to the LED illumination sample at the center of the LED array;

Indicates the initial phase, usually 0;

Step S112, obtaining a sub-aperture low-resolution complex amplitude estimation value:

For the incident angle determined by the wave vector k _i,j , the corresponding sub-aperture spectrum information is intercepted from the current sample spectrum using the principle of translation effect, and then the two-dimensional inverse Fourier transform is used to obtain the sub-aperture low-resolution complex amplitude estimation value under the corresponding LED:

表示更新位置序列(i,j)对应子孔径的样本频谱；(u_i,j,v_i,j)表示位置序列(i,j)在频域中对样本的平移效应，P(u+u_i,j,v+v_i,j)表示平移效应结果；

表示第n次迭代对应位置序列(i,j)的子频谱复振幅估计值；Where n represents the current number of iterations (1≤n≤n _max ), n _max represents the set maximum number of iterations,

represents the sample spectrum of the sub-aperture corresponding to the updated position sequence (i, j); (u _{i, j} , vi _{, j} ) represents the translation effect of the position sequence (i, j) on the sample in the frequency domain, and P(u+u _{i, j} ,v+vi _{, j} ) represents the translation effect result;

Represents the estimated value of the complex amplitude of the sub-spectrum corresponding to the position sequence (i, j) at the nth iteration;

Step S113, updating the sub-spectrum complex amplitude estimation value:

The amplitude part of the sub-spectrum complex amplitude estimate is replaced by the low-resolution intensity image actually collected at the position sequence (i, j), and the phase part is kept unchanged. The update formula of the sub-spectrum complex amplitude estimate is:

Step S114, updating the sub-aperture region in the spectrum:

The sub-spectrum complex amplitude estimate of the updated n-th iteration corresponding to the position sequence (i, j) is transformed using the two-dimensional inverse Fourier transform to obtain the corresponding sub-aperture spectrum information:

Afterwards, the overall sample spectrum is updated according to the sub-aperture positions corresponding to the position sequence (i, j):

表示更新后的样本频谱；该式表示，在进行一次子频谱复振幅后，整体样本频谱只对位置序列(i,j)子孔径部分更新频谱信息，其余部分则保持不变。In the formula,

Represents the updated sample spectrum; this formula indicates that after a sub-spectrum complex amplitude is performed, the overall sample spectrum only updates the spectrum information of the sub-aperture part of the position sequence (i, j), and the rest remains unchanged.

Step S115, updating the sub-aperture spectrum information corresponding to all position sequences {(i,j)|1≤i≤m,1≤j≤n} according to steps S112 to S114, and obtaining the n-th iterative updated sample spectrum as the reconstruction result O ⁿ (u,v);

其中，每一次迭代循环的初始值为上一次循环更新后的样本频谱结果。Step S116, iterate the steps S112 to S115 until the reconstruction result converges or the loop ends, and obtain the final reconstruction result, that is, the reconstructed sample spectrum with position error, which is expressed as

The initial value of each iteration cycle is the sample spectrum result updated in the previous cycle.

The method for judging the convergence of reconstruction results is:

The loss judgment function corresponding to the nth iteration cycle is:

In the formula, when ^En reaches the minimum value or ( ^En - ^En-1 )/En ^-1 is less than the set threshold, the reconstruction result is considered to have converged.

并将样本高分辨复振幅

分为强度部分与位相部分：Step S117, the reconstructed sample spectrum with position error is subjected to a two-dimensional inverse Fourier transform to obtain a reconstructed sample high-resolution complex amplitude with position error.

The sample high-resolution complex amplitude

It is divided into intensity part and phase part:

表示重建的带位置误差的样本高分辨复振幅，

表示重建的带位置误差的样本高分辨复振幅中的强度部分，

表示重建的带位置误差的样本高分辨复振幅中的位相部分。In the formula, ~ represents the relevant results with position error, that is,

represents the reconstructed sample high-resolution complex amplitude with position error,

represents the intensity part of the reconstructed sample high-resolution complex amplitude with position error,

Represents the phase part of the reconstructed sample high-resolution complex amplitude with position error.

Step S200, high-resolution intensity image acquisition and regional registration under a high-power objective lens: after replacing the low-power objective lens with a high-power objective lens and completing the focusing operation, a high-resolution intensity image of the sample under the high-power objective lens is acquired; then, the low-resolution intensity image and the high-resolution intensity image are matched using an image registration algorithm and the regions representing the same information in the two images are framed, and the regions representing the same information in the high-resolution intensity image are interpolated to obtain a registered high-resolution intensity image;

Due to the difference in field of view, the display area of the high-resolution intensity image acquired under the high-power objective is a sub-area of the display area of the low-resolution intensity image captured under the low-power objective. In order to match the field of view of the two images and characterize the same information area, the high-resolution intensity image needs to be registered, which specifically includes the following sub-steps:

Step S201, performing a coarse matching of the display areas of the low-resolution intensity image and the high-resolution intensity image through a correlation matching operation, and obtaining a coarse matching area of the two images, that is, an approximate area range of the two images;

Step S202, performing a binarization operation on the rough matching area, extracting the coordinates of the top, bottom, left and right four points used to describe the boundary of the area, and the area selected by the four points is identified as the area representing the same information;

Step S203, based on the FPM reconstructed image resolution, the regions representing the same information in the high-resolution intensity image are interpolated to obtain the registered high-resolution intensity image I ^hr . It can be seen that since the present invention interpolates the regions representing the same information, in the subsequent position correction process, only the information in the regions representing the same information in the high-resolution intensity image is correlated.

Step S300, generating a virtual low-resolution intensity image of the corresponding position of the LED: using the high-resolution intensity with position error and the phase complex amplitude and the registered high-resolution intensity image to form a new spectrum diagram of the sample; according to the ideal frequency domain position of each LED in the programmable LED array, obtaining the sub-spectrum information corresponding to the ideal frequency domain position of each LED from the new spectrum diagram, and obtaining a series of virtual low-resolution intensity images by Fourier inverse transforming the sub-spectrum information; specifically comprising the following sub-steps:

位相部分保持不变，经二维傅里叶变换后获得虚拟频谱O^vir：Step S301, using the registered high-resolution acquisition image I ^hr to replace the intensity part of the high-resolution intensity and phase complex amplitude with position error

The phase part remains unchanged, and the virtual spectrum ^Ovir is obtained after two-dimensional Fourier transform:

Step S302, following the FPM acquisition process in step S100, sub-aperture sampling is performed according to the ideal frequency domain position (u+u _i,j ,v+vi _,j ) corresponding to the position sequence (i,j):

表示位置序列(i,j)虚拟的采样子孔径频谱；In the formula,

represents the virtual sampling sub-aperture spectrum of the position sequence (i, j);

Step S303: Perform a two-dimensional Fourier transform on the virtual sampling sub-aperture spectrum to obtain a corresponding virtual low-resolution intensity map.

Step S304, repeat steps S302 to S303 until all position sequence LEDs are sampled, generating a series of virtual low-resolution intensity images

Step S400, search for position parameters of independent LEDs in the array: based on the fact that the translation effect of the corresponding position of the LED in the spatial domain is equivalent to the translation effect of the aperture position in the frequency domain, a frequency domain position error model is introduced for each LED; thus, the current ideal frequency domain position of the LED is taken as the search center, the search area is specified, and the optimal position parameter search is performed using a two-dimensional particle swarm algorithm; when the difference between the generated virtual low-resolution intensity image and the currently actually collected low-resolution intensity image is the smallest, it is considered that this position is the actual frequency domain position corresponding to the LED, and the position correction of the LED is completed; specifically, the following sub-steps are included:

Step S401, in order to accurately correct the position error of each LED, the present invention introduces an independent frequency domain position error model (Δu _i,j , Δv _i,j ) for each position sequence (i,j). Then, the corresponding virtual low-resolution intensity image after the introduction of the frequency domain position error model is expressed as:

Where Δu _i,j and Δv _i,j represent the position error components of the position sequence (i, j) in the two-dimensional frequency domain;

Step S402, using the two-dimensional correlation function Corr2 to evaluate the difference between the virtual low-resolution intensity image corresponding to the position sequence (i, j) and the low-resolution intensity image actually acquired in step S100:

表示位置序列(i,j)对应的虚拟低分辨强度图像的均值，

表示位置序列(i,j)对应的步骤S100实际采集的低分辨强度图像的均值。In the formula,

represents the mean of the virtual low-resolution intensity image corresponding to the position sequence (i, j),

Represents the mean value of the low-resolution intensity image actually acquired in step S100 corresponding to the position sequence (i, j).

Step S403, using a two-dimensional particle swarm algorithm (PSO) to search for the best frequency domain position parameters of the position sequence (i, j) in the frequency domain:

(Δu _i,j ,Δv _i,j )=argmin(fitness)

fitness＝1-Corr2 for(i,j)

Where fitness represents the fitness function describing the particle swarm, i.e., the loss function. When the fitness function value of all particles converges to the minimum value, it means that the difference between the generated virtual low-resolution intensity image and the currently actually collected low-resolution intensity image is the smallest. At this time, (Δu _i,j ,Δv _i,j ) is the frequency domain correction parameter of the position sequence (i,j). Then, the frequency domain correction parameter is used to correct the frequency domain position parameter (coordinate component in the two-dimensional frequency domain) of the position sequence (i,j) to obtain the corrected optimal frequency domain position parameter.

Step S404, using the corrected optimal frequency domain position parameters, perform FPM acquisition and reconstruction again according to the method of step S100, that is, obtain the corrected and reconstructed high-resolution intensity and phase complex amplitude with position error.

So far, in the method for accurate correction of the frequency domain light source position based on Fourier stack microscopy imaging of the present invention, (1) in order to solve the problem of inaccurate correction position, an independent frequency domain position error model is introduced for each LED in the LED array for accurate correction; (2) in order to solve the problem that the correction process needs to participate in the image reconstruction process and consumes too much time, a high-resolution intensity image collected by a high-power objective lens is introduced as a reference image to simulate the image acquisition process of Fourier stack microscopy imaging technology, and generate a virtual low-resolution intensity image corresponding to each LED. Combining the virtual and actual collected images, a loss function is established, and a two-dimensional particle swarm algorithm is used to search for the optimal frequency domain position parameters corresponding to each LED. Therefore, the present invention utilizes the high-resolution intensity image under the high-power objective lens and the frequency domain position error model to realize accurate correction of the frequency domain position of the LED light source based on the principle of Fourier stack imaging, which has the advantages of simple implementation, high robustness, and accurate position.

The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. For those skilled in the art, the present invention may have various modifications and variations. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims (2)

1. A method for accurate correction of the position of a frequency domain light source based on Fourier stack microscopy imaging, characterized in that it comprises the following steps: Step S100, FPM acquisition and reconstruction with position error under low-power lens: LEDs in the programmable LED array are illuminated one by one under the low-power objective lens, and a series of low-resolution intensity images of the sample at the corresponding positions of different LEDs are captured on the camera surface; then the captured low-resolution intensity images are reconstructed in the frequency domain using a reconstruction algorithm to obtain reconstructed high-resolution intensity and phase complex amplitude with position error; Step S200, high-resolution intensity image acquisition and regional registration under a high-power objective lens: after replacing the low-power objective lens with a high-power objective lens and completing the focusing operation, a high-resolution intensity image of the sample under the high-power objective lens is acquired; then, the low-resolution intensity image and the high-resolution intensity image are matched using an image registration algorithm and the regions representing the same information in the two images are framed, and the regions representing the same information in the high-resolution intensity image are interpolated to obtain a registered high-resolution intensity image; Step S300, generating a virtual low-resolution intensity image of the corresponding position of the LED: using the high-resolution intensity with position error and the phase complex amplitude and the registered high-resolution intensity image to form a new spectrum diagram of the sample; according to the ideal frequency domain position of each LED in the programmable LED array, obtaining the sub-spectrum information corresponding to the ideal frequency domain position of each LED from the new spectrum diagram, and obtaining a series of virtual low-resolution intensity images by inverse Fourier transforming the sub-spectrum information; Step S400, search for position parameters of independent LEDs in the array: based on the fact that the translation effect of the corresponding position of the LED in the spatial domain is equivalent to the translation effect of the aperture position in the frequency domain, a frequency domain position error model is introduced for each LED; thus, the current ideal frequency domain position of the LED is taken as the search center, the search area is specified, and the optimal position parameter search is performed using a two-dimensional particle swarm algorithm; when the difference between the generated virtual low-resolution intensity image and the currently actually collected low-resolution intensity image is the smallest, it is considered that this position is the actual frequency domain position corresponding to the LED, and the position correction of the LED is completed; The method for collecting FPM with position error under low-power microscope in step S100 includes the following sub-steps: Step S101, using a programmable LED array with a specification of m×n as a light source; assuming that the interval between adjacent LEDs in the programmable LED array is d, the ideal position P _i,j of the LED in the i-th row and j-th column (1≤i≤m, 1≤j≤n) is expressed as: P _i,j = (x _i,j ,y _i,j ) = (i·d,j·d) Where x _i,j and y _i,j represent the two-dimensional coordinate component values of the horizontal plane where the LED is located in the airspace; Step S102, assuming that the light source is far enough from the sample, the illumination light wave emitted from the LED can be regarded as a monochromatic plane wave. Therefore, when a certain LED is used to illuminate the sample, the corresponding wave vector _{k i,j} is expressed as:

In the formula,

and

Respectively represent the wave vector component values of wave vector k _i,j in the x and y directions; λ represents the central wavelength value of LED radiation; h represents the vertical distance between the sample and the LED currently illuminating the sample; Step S103, normalize the intensity value of the incident light field of the sample to 1, then the spectrum emitted from the sample is expressed as: F{e(r)}＝F{o(r)exp(ik _i,j r)}＝O(kk _i,j ) Wherein, F{} represents the two-dimensional Fourier transform, e(r) represents the complex amplitude emitted from the sample, o(r) represents the complex amplitude modulation function of the sample space, that is, the complex amplitude function of the sample; o(k) represents the spectrum of the sample in the frequency domain; O(kk _i,j ) represents that the spectrum emitted by the sample is equivalent to the spectrum center of the sample being translated to the wave vector k _i,j ; Step S104, the complex amplitude emitted from the sample is low-pass filtered in the frequency domain by the pupil function P(k) of the low-power objective lens. The low-pass filtering process is expressed as: Gi _,j (k)＝O(kki _,j )P(k) Where G _i,j (k) represents the sample spectrum after passing through the low-power objective lens; Step S105: the sample spectrum is propagated to the image plane for imaging and is captured by the camera and converted into a low-resolution intensity image.

It is expressed as:

Where, g _i,j (r) represents the complex amplitude of the sample on the image plane, and F ^-1 {} represents the two-dimensional inverse Fourier transform; Step S106, lighting up the LEDs in the programmable LED array one by one, and executing steps S102 to S105 for each lit LED, thereby capturing a series of low-resolution intensity images of the sample at positions corresponding to different LEDs.

The reconstruction method in step S100 includes the following sub-steps: Step S111, spectrum initialization, that is, the result of interpolation of the low-resolution intensity image collected at the center position of the programmable LED array is used as the initial value of the intensity part of the reconstructed sample spectrum, and the zero phase is selected as the initial value of the relevant phase, and then combined with the pupil function of the low-power objective lens, the initial value of the reconstructed sample spectrum is solved as:

In the formula,

represents the initial value of the reconstructed sample spectrum; u and v represent the coordinate components in the two-dimensional frequency domain; R() represents the related interpolation process, I _center (r) represents the low-resolution intensity image collected when the LED is illuminated at the center of the LED array, and k _center represents the wave vector corresponding to the LED illumination sample at the center of the LED array;

represents the initial phase; Step S112, obtaining the sub-aperture low-resolution complex amplitude, that is, using the principle of translation effect to intercept the corresponding sub-aperture sample spectrum from the current sample spectrum, and then using two-dimensional inverse Fourier transform to obtain the sub-aperture low-resolution complex amplitude estimation value under the corresponding LED:

Where n represents the current number of iterations, 1≤n≤n _max , and n _max represents the set maximum number of iterations.

represents the sample spectrum of the sub-aperture corresponding to the updated position sequence (i, j); (u _{i, j} , vi _{, j} ) represents the translation effect of the position sequence (i, j) on the sample in the frequency domain, and P(u+u _{i, j} ,v+vi _{, j} ) represents the translation effect result;

Represents the estimated value of the complex amplitude of the sub-spectrum corresponding to the position sequence (i, j) at the nth iteration; Step S113, updating the sub-spectrum complex amplitude estimation value: The amplitude part of the sub-spectrum complex amplitude estimate is replaced by the low-resolution intensity image actually collected at the position sequence (i, j), and the phase part is kept unchanged. The update formula of the sub-spectrum complex amplitude estimate is:

In the formula,

Represents the updated sub-spectrum complex amplitude estimate of the position sequence (i, j) corresponding to the nth iteration; Step S114, using two-dimensional inverse Fourier transform to obtain the corresponding sub-aperture spectrum information for the sub-spectrum complex amplitude estimation value of the nth iteration corresponding position sequence (i, j):

In the formula,

Indicates the corresponding

Sub-aperture spectrum information, F{} represents two-dimensional Fourier transform; Afterwards, the overall sample spectrum is updated according to the sub-aperture positions corresponding to the position sequence (i, j):

In the formula,

represents the updated sample spectrum; Step S115, updating the sub-aperture spectrum information corresponding to all position sequences {(i,j)|1≤i≤m,1≤j≤n} according to steps S112 to S114, and obtaining the n-th iterative updated sample spectrum as the reconstruction result O ⁿ (u,v); Step S116, iterate the steps S112 to S115 until the reconstruction result converges or the loop ends, and obtain the final reconstruction result, that is, the reconstructed sample spectrum with position error, which is expressed as

The initial value of each iteration cycle is the sample spectrum result after the previous cycle update; Step S117, the reconstructed sample spectrum with position error is subjected to a two-dimensional inverse Fourier transform to obtain a reconstructed sample high-resolution complex amplitude with position error.

The sample high-resolution complex amplitude

It is divided into intensity part and phase part:

In the formula,

represents the reconstructed sample high-resolution complex amplitude with position error,

represents the intensity part of the reconstructed sample high-resolution complex amplitude with position error,

Represents the phase part of the reconstructed sample high-resolution complex amplitude with position error; The method for determining whether the reconstruction result converges in step S116 is: The loss judgment function corresponding to the nth iteration cycle is:

In the formula, when ^En reaches the minimum value or ( ^En - ^En-1 )/En ^-1 is less than the set threshold, the reconstruction result is considered to have converged; Step S200 includes the following sub-steps: Step S201, performing a coarse matching of the display area of the low-resolution intensity image and the high-resolution intensity image through a correlation matching operation to obtain a coarse matching area of the two images; Step S202, performing a binarization operation on the rough matching area, extracting the coordinates of the top, bottom, left and right four points used to describe the boundary of the area, and the area selected by the four points is identified as the area representing the same information; Step S203, taking the FPM reconstructed image resolution as a reference, interpolating the regions representing the same information in the high-resolution intensity image to obtain a registered high-resolution intensity image I ^hr ; Step S300 includes the following sub-steps: Step S301, using the registered high-resolution acquisition image I ^hr to replace the intensity part of the high-resolution intensity and phase complex amplitude with position error

The phase part remains unchanged, and the virtual spectrum ^Ovir is obtained after two-dimensional Fourier transform:

Step S302, following the FPM acquisition process in step S100, sub-aperture sampling is performed according to the ideal frequency domain position (u+u _i,j ,v+vi _,j ) corresponding to the position sequence (i,j):

In the formula,

represents the virtual sampling sub-aperture spectrum of the position sequence (i, j); Step S303: Perform a two-dimensional Fourier transform on the virtual sampling sub-aperture spectrum to obtain a corresponding virtual low-resolution intensity map.

Step S304, repeat steps S302 to S303 until all position sequence LEDs are sampled, generating a series of virtual low-resolution intensity images

Step S400 includes the following sub-steps: Step S401, introduce an independent frequency domain position error model (Δu _i,j , Δv _i,j ) for each position sequence (i,j). Then, the corresponding virtual low-resolution intensity image after the frequency domain position error model is introduced is expressed as:

Where Δu _i,j and Δv _i,j represent the position error components of the position sequence (i, j) in the two-dimensional frequency domain; Step S402, using the two-dimensional correlation function Corr2 to evaluate the difference between the virtual low-resolution intensity image corresponding to the position sequence (i, j) and the low-resolution intensity image actually acquired in step S100:

In the formula,

represents the mean of the virtual low-resolution intensity image corresponding to the position sequence (i, j),

represents the mean value of the low-resolution intensity image actually acquired in step S100 corresponding to the position sequence (i, j); Step S403, using a two-dimensional particle swarm algorithm (PSO) to search for the best frequency domain position parameters of the position sequence (i, j) in the frequency domain: (Δu _i,j ,Δv _i,j )=argmin(fitness) fitness＝1-Corr2for(i,j) Where fitness represents the fitness function describing the particle swarm, that is, the loss function. When the fitness function value of all particles converges to the minimum value, it means that the difference between the generated virtual low-resolution intensity image and the currently actually collected low-resolution intensity image is the smallest. At this time, (Δu _i,j ,Δv _i,j ) is the frequency domain correction parameter of the position sequence (i,j). Then, the frequency domain correction parameters are used to correct the frequency domain position parameters of the position sequence (i,j) to obtain the corrected optimal frequency domain position parameters. 2. The method for accurate correction of frequency domain light source position based on Fourier stack microscopy imaging according to claim 1, characterized in that step S404 further comprises: Step S404, using the corrected optimal frequency domain position parameters, perform FPM acquisition and reconstruction again according to the method of step S100, that is, obtain the corrected and reconstructed high-resolution intensity and phase complex amplitude with position error.

Images