Phase recovery method based on optically limited imaging system and analytic constraint

Abstract

The invention provides a phase recovery method based on an optically limited imaging system. The phase recovery method comprises the following steps: acquiring first and second diffraction plane intensities; according to whether the intensity of the second diffraction plane is cut off or not, the phase recovery imaging method of the first type or the second type is executed, phase recovery imaging comprises multiple iterations, constraint is applied to a diffraction plane light field through a feedback mechanism in each iteration, and analytical constraint is applied. Wherein for the first type of phase recovery imaging, a spatial spectrum limited constraint and an analytic constraint of a first diffraction plane are applied; for a second type of phase recovery imaging, a limited illumination / limited field of view constraint and a second diffraction plane analytic constraint are applied. According to the phase recovery method, the defect that uniqueness of an existing iterative phase recovery method depends on object priori and multiple exposure is overcome by utilizing the constraint of an optically limited imaging system, and meanwhile, the problems that a traditional iterative projection algorithm is poor in reconstruction precision, too low in convergence speed and the like are solved by utilizing analytic constraint.

Description

Technical Field

[0001] The present invention relates to the field of phase recovery imaging, and more specifically to a phase recovery method based on an optically restricted imaging system and analytical constraints. Background Art

[0002] The imaging process of the camera projects a 3D scene onto a 2D image to capture the intensity distribution of the reflected light. While this approach is commonly used to record objects, it is difficult to effectively represent 3D scenes or uniformly absorbing objects. In physics and mathematics, "phase" plays a vital role as a real-valued scalar that determines the relative position of each point in a periodic signal. Due to the wave nature of light, the optical delay that occurs when light is projected onto a scattering medium can be defined as the phase change of the periodic light wave during the propagation process, which contains important optical path information and thus provides a means to quantify the depth information of the imaged object. By acquiring both amplitude and phase simultaneously, all information about the light field of the scatterer can be reconstructed.

[0003] Phase imaging technology is a powerful label-free imaging method that can recover phase information from intensity information detection. Interference method is a classic phase imaging method, mainly including digital holographic microscopy, phase contrast microscopy, differential interference microscopy, optical diffraction tomography (ODT) and other technologies. This type of method can visualize phase information by interfering with a known reference beam or two polarized beams. However, in order to achieve quantitative recovery of the phase information of the object, this type of method usually needs to be combined with phase shift technology or tomography technology to collect multiple frames of interference patterns, which makes this type of technology have the disadvantage of long measurement time and is not conducive to the application requirements of real-time imaging. Iterative phase recovery technology is a non-interference method, mainly including coherent diffraction imaging (CDI), self-interference holographic imaging, ptychography, Fourier ptychographic microscopy (FPM) and other technologies, providing another efficient imaging method for quantitative phase imaging technology. As a lensless imaging method, CDI records the diffraction image of the object, uses an iterative projection algorithm and imposes constraints such as the size of the object to inversely solve the phase of the object. Similarly, self-interference holographic imaging records the diffraction image after the interference of the object light and the illumination light, uses an iterative algorithm to obtain the phase and solve the twin image problem. These technologies have the ability to image in a single frame and can obtain good temporal resolution. However, the uniqueness of the solution of this type of method depends on oversampling, and its space-bandwidth product (i.e., information capacity) is inversely proportional to the sampling rate, so oversampling reduces the information capacity of the imaging system. In addition, the speed and robustness of phase recovery are heavily dependent on the convergence of the algorithm, making the accurate spatial size of the object and other prior constraints very important. However, the accurate object size is difficult to obtain in the application of transparent biological imaging. Even if the accurate object size is obtained, the convergence effect will also obtain different results with different initializations, and the success rate cannot be effectively guaranteed. Therefore, single recovery cannot be achieved.

[0004] In order to improve the robustness of phase recovery, speed up the calculation, and solve the problem of single-shot recovery, Rodenburg et al. proposed stacked diffraction imaging in 2004. Multiple diffraction images were obtained by scanning different overlapping local illumination areas of the sample, and the redundant information of the overlapping area was used for iteration, and a large field of view image was obtained at the same time. In 2005, Pedrini et al. proposed self-interference holographic imaging for multi-transmission distance measurement. Iterative phase recovery was performed by measuring the diffraction patterns of multiple transmission distances. In 2013, Guoan Zheng et al. proposed FPM technology based on a microscope system. Multiple microscopic images were obtained by overlapping multi-angle illumination, and the redundant information of the overlapping area was used for iteration. At the same time, a low-power objective lens was used to obtain a large field of view image, and a synthetic aperture was used to improve the spatial resolution. In 2016, Horstmeyer et al. proposed ODT technology without interference. Multiple images were obtained by multi-angle illumination. In the process of iterating the redundant information of the overlapping area, diffraction tomography models such as the Born approximation were introduced to achieve three-dimensional imaging of thick objects. These technologies are widely used in living cell imaging.

[0005] However, all of the above technologies require redundant measurements of multiple frames, thus reducing the temporal resolution of imaging. Although the existing non-interference phase recovery technology reduces the system complexity of interferometric phase imaging, it is an ill-posed mathematical problem due to the existence of fuzzy solutions to the phase problem. It requires additional redundant information such as oversampling, object size priors, or stacking and aperture synthesis to achieve good algorithm convergence. The need for redundant information in phase recovery reduces the efficiency of acquiring imaging image information, requiring phase imaging to collect multiple images to meet the uniqueness requirement.

[0006] In fact, all existing optical imaging systems are limited imaging systems. Due to the limited size of the imaging lens, the optical imaging system is actually a low-pass filter that filters out the high-frequency components of the object and only allows low-frequency components within a certain range to pass through the imaging system. The camera system and the microscope system are the two most typical optical limited imaging systems, and their structures are as follows: Figure 1 , Figure 2 shown.

[0007] The limitations of the camera system lie in the limited size of the imaging lens, the aperture size, and the field of view angle. On the Fourier plane, the limited size of the imaging lens and aperture filters out the high-frequency components of the object and determines the pupil function of the imaging system. On the real image plane, the limited size of the field of view angle limits the size of the real space that the image sensor can cover. Suppose the focal length of the camera is f, the half-height width of the sensor is h, the field of view angle is θ, the object moment is u, and the field of view height is FOV, then there is the following relationship:

[0008]

[0009]

[0010] The limitation of the microscope system lies in the limited pupil size, numerical aperture and field of view of the microscope objective. In the Fourier plane, the limited numerical aperture filters out the high-frequency components of the object and determines the pupil function of the imaging system. For the coherent light illumination system, let the cutoff frequency be f c , the numerical aperture of the objective is NA, and the illumination wavelength is λ, then there is the following relationship:

[0011]

[0012] On the real image plane, the field number FN of the objective lens and the magnification M jointly determine the field of view of the microscope system. Assume that the maximum field of view diameter of the microscope system is d FOV , then there is the following relationship:

[0013]

[0014] In order to collect the intensity distribution of the image plane and the intensity distribution of the Fourier plane of the imaging target, such as Figure 3 An optically confined imaging system is shown, which includes at least the following components: a light source 101, a first imaging lens 102, a beam splitter 103, a second imaging lens 104, a second two-dimensional array detector 105 (to obtain Fourier surface intensity, i.e., second diffraction plane intensity), and a first two-dimensional array detector 106 (to obtain image plane intensity, i.e., first diffraction plane intensity) that receives light from the first imaging lens 102 and reflected by the beam splitter 103.

[0015] The first imaging lens 102 may be a photographic system or a microscope system. The specific acquisition process is as follows: when acquiring the image plane intensity, the incident light sequentially passes through the first imaging lens 102, the beam splitter 103, and the first two-dimensional array detector 106; when acquiring the Fourier surface intensity, the incident light sequentially passes through the first imaging lens 102, the beam splitter 103, the second imaging lens 104, and the second two-dimensional array detector 105.

[0016] In addition, in recent years, many phase recovery technologies have utilized analytical constraints, among which the most representative are the off-axis holographic systems that utilize the Kramer-Krong (KK) relationship and the synthetic aperture technology that also utilizes the KK relationship. The signal is separated into real and imaginary parts by operating the real image plane light field through a nonlinear logarithmic function, and then the phase information is directly calculated by taking the logarithm of the real image plane intensity through the KK relationship, and the spatial bandwidth product of the traditional off-axis system is increased by 3 to 4 times. However, the above-mentioned technology also has shortcomings. For example, the former uses a complex interference system, which requires the coherence of the reference light and the object light, thereby increasing the complexity of the system. Secondly, although the spatial bandwidth product is increased, it is still lower than the spatial bandwidth product of the present invention. Similarly, although the synthetic aperture technology represented by the KK relationship can make up for the defect of insufficient spatial bandwidth product through multiple exposures, multiple exposures lead to a decrease in its temporal resolution.

[0017] Therefore, it is necessary to propose a phase retrieval method and algorithm based on optically constrained imaging systems and analytical constraints, so as to not only solve the defects of the existing iterative phase retrieval method that the uniqueness depends on object priors and multiple exposures, but also solve the problems of poor reconstruction accuracy and slow convergence speed of traditional iterative projection algorithms. Summary of the invention

[0018] The purpose of the present invention is to provide a phase recovery method based on an optically constrained imaging system and analytical constraints, so as to not only solve the defects of the existing iterative phase recovery method that the uniqueness depends on object priors and multiple exposures, but also solve the problems of poor reconstruction accuracy and slow convergence speed of traditional iterative projection algorithms.

[0019] In order to achieve the above object, the present invention provides a phase recovery method based on an optically constrained imaging system and analytical constraints, comprising:

[0020] S1: collecting a first diffraction plane intensity and a second diffraction plane intensity of an imaging target through the optically limited imaging system;

[0021] S2: if the intensity of the second diffraction plane is truncated, the first type of phase recovery imaging method is executed, otherwise, the second type of phase recovery imaging method is executed;

[0022] The first type of phase recovery imaging method specifically includes:

[0023] S21: Based on the collected first diffraction plane intensity I(x n ,y m ) and the second diffraction plane intensity The amplitude of the light field in the first diffraction plane collected is obtained by square root calculation |f(x n ,y m)| and the amplitude of the light field collected in the second diffraction plane

[0024] S22: Obtain the initial phase of the light field in the first diffraction plane through the Hilbert transform characteristics of the light field analysis signal

[0025] S23: The initial phase of the first diffraction plane light field The amplitude of the light field in the first diffraction plane collected |f(x n ,y m )| multiply to get the first diffraction plane light field f1(x n ,y m ); At the same time, for the second diffraction plane light field of the first iteration Initialize;

[0026] S24: Input the first diffraction plane light field f for the i-th iteration i (x n ,y m ) performs the diffraction transformation calculation between the first diffraction plane and the second diffraction plane to obtain the output second diffraction plane light field of the i-th iteration

[0027] S25: Using the feedback mechanism to output the second diffraction plane light field of the i-th iteration Apply the amplitude constraint and spatial spectrum restriction constraint of the second diffraction plane light field to obtain the input second diffraction plane light field of the i+1th iteration

[0028] S26: Input second diffraction plane light field for the i+1th iteration Perform the diffraction transformation calculation between the second diffraction plane and the first diffraction plane to obtain the output first diffraction plane light field f′ of the i+1th iteration i+1 (x n ,y m );Then the output first diffraction plane light field f′ of the i+1th iteration i+1 (x n ,y m ) is projected onto the amplitude of the light field of the first diffraction plane collected |f(x n ,y m )|, and obtain the first diffraction plane light field of the input to be corrected in the i+1th iteration

[0029] S27: Input the first diffraction plane light field to be corrected for the i+1th iteration Applying analytical constraints, we can obtain the input first diffraction plane light field f of the i+1th iteration i+1 (xn ,y m ), including:

[0030] S271: Input the first diffraction plane light field to be corrected for the i+1th iteration Perform Shannon interpolation to obtain the first diffraction plane light field to be corrected after Shannon interpolation

[0031] S272: Determine whether the number of Shannon interpolation k reaches the total number of Shannon interpolation N′. If so, input the corrected light field after Shannon interpolation into the first diffraction plane light field. As the input of the i+1th iteration, the first diffraction plane light field f i+1 (x n ,y m ); otherwise, the Shannon interpolated light field to be corrected is input into the first diffraction plane light field As the first diffraction plane light field to be corrected input for the i+1th iteration And return to step S271;

[0032] S28: add 1 to the current iteration number i as the current iteration number i of the new round, and return to step S24 to perform the next round of iteration; until the iteration end condition is met, the iteration ends at this time, and the input first diffraction plane light field f of the i+1th iteration is i+1 (x n ,y m ) as the reconstructed first diffraction plane light field;

[0033] The second type of phase recovery imaging method specifically includes:

[0034] S21': Based on the collected first diffraction plane intensity I(x n ,y m ) and the second diffraction plane intensity The amplitude of the light field in the first diffraction plane collected is obtained by square root calculation |f(x n ,y m )| and the amplitude of the light field collected in the second diffraction plane

[0035] S22': The initial phase of the light field in the first diffraction plane is obtained by the Hilbert transform characteristics of the light field analysis signal

[0036] S23': The initial phase of the first diffraction plane light field The amplitude of the light field collected in the second diffraction plane Multiply to get the input second diffraction plane light field of the first iteration At the same time, the first diffraction plane light field f1(x n ,y m ) to initialize;

[0037] S24': Input second diffraction plane light field for the i-th iteration Perform the diffraction transformation calculation between the second diffraction plane and the first diffraction plane to obtain the output first diffraction plane light field f of the i-th iteration i ′(x n ,y m );

[0038] S25': Using the feedback mechanism to output the first diffraction plane light field f' of the i-th iteration i (x n ,y m ) applies the amplitude constraint of the first diffraction plane light field and the constraint of limited illumination or limited field of view to obtain the input first diffraction plane light field f of the i+1th iteration i+1 (x n ,y m );

[0039] S26': Input first diffraction plane light field f for the i+1th iteration i+1 (x n ,y m ) performs the diffraction transformation calculation between the first diffraction plane and the second diffraction plane to obtain the output second diffraction plane light field of the i+1th iteration Then the output second diffraction plane light field of the i+1th iteration is The amplitude of the light field projected onto the second diffraction plane collected Get the second diffraction plane light field to be corrected input for the i+1th iteration

[0040] S27': Input the second diffraction plane light field to be corrected for the i+1th iteration Apply an analytical constraint, which is used to make the amplitude of the collected second diffraction plane light field contain phase information, and obtain the input second diffraction plane light field of the i+1th iteration Specifically include:

[0041] S271': Input the second diffraction plane light field to be corrected for the i+1th iteration Perform Shannon interpolation to obtain the second diffraction plane light field to be corrected after Shannon interpolation

[0042] S272': Determine whether the number of Shannon interpolation k reaches the total number of Shannon interpolation N'. If so, input the Shannon interpolated value to be corrected into the second diffraction plane light field. As the input second diffraction plane light field of the i+1th iteration Otherwise, the Shannon interpolated correction is input into the second diffraction plane light field As the second diffraction plane light field to be corrected input for the i+1th iteration And return to step S271;

[0043] S28': add 1 to the current iteration number i as the current iteration number i of the new round, and return to step S24' to perform the next round of iteration; until the iteration end condition is met, the iteration ends at this time, and the input first diffraction plane light field f of the i+1th iteration is i+1 (x n ,y m ) as the reconstructed first diffraction plane light field.

[0044] The first diffraction plane is an image plane, and the second diffraction plane is a Fourier plane;

[0045] In the step S1, the second diffraction plane intensity distribution is a Fourier surface space spectrum, and whether the imaging target is a strong scattering object or a weak scattering object is determined according to whether the Fourier surface space spectrum is truncated;

[0046] In step S2, the first type of phase recovery imaging method is a phase recovery imaging method for strongly scattering objects, and the second type of phase recovery imaging method is a phase recovery method for weakly scattering objects.

[0047] In step S24, step S26, step S24' and step S26', the diffraction transformation between the first diffraction plane and the second diffraction plane is calculated as a Fraunhofer diffraction transformation, and the diffraction transformation between the second diffraction plane and the first diffraction plane is calculated as an inverse Fraunhofer diffraction transformation.

[0048] Light scattered by the imaging target in an optically constrained imaging system causes a hard truncation of the Fourier plane spatial spectrum or the imaging field of view on the image plane due to the limited numerical aperture of the optical system; the limited numerical aperture of the optical system includes the size limitation of the light source, the aperture limitation of the objective lens, the pupil, numerical aperture and imaging field angle of the imaging lens, the size limitation of the detector's photosensitive area, the size limitation of the rectangular or circular aperture, and the probe size limitation of the fiber optic detector;

[0049] Among them, the finite numerical aperture of the optical system corresponding to the limited Fourier surface space spectrum includes the pupil and / or numerical aperture size limitation of the imaging lens; the limited Fourier surface space spectrum causes hard truncation of the Fourier surface space spectrum; the finite numerical aperture of the optical system corresponding to the limited illumination includes the size limitation of the light source, and the finite numerical aperture of the optical system corresponding to the limited field of view includes the pupil size limitation of the imaging lens and / or the limitation of the imaging field of view angle; limited illumination or limited field of view causes hard truncation of the imaging field of view of the image plane.

[0050] The step S22 specifically includes:

[0051] S221: The amplitude of the first diffraction plane light field collected |f(x n ,y m )| Perform Fourier transform to obtain the Fourier transform F of the real part of the light field to be restored R ;

[0052] S222: Establishing the Fourier transform F of the real part of the light field to be restored R The constraint equations for the real and imaginary parts of ;

[0053] S223: Based on the Fourier transform F of the real part of the light field to be restored R The constraint equations of the real and imaginary parts of the intersection point are solved to obtain the coordinates (R{F I},I{F I});

[0054] S224: According to the intersection coordinates (R{F I},I{F I}), calculate the initial phase

[0055] In step S222, the Fourier transform F of the real part of the light field to be restored R The constraint equations for the real and imaginary parts of are:

[0056]

[0057] Among them, I{F R} is the imaginary part of the Fourier transform of the real part of the light field to be restored, R{F I} is the real part of the Fourier transform of the imaginary part of the light field to be restored, I{F R} is the real part of the Fourier transform of the real part of the light field to be restored, I{F I} is the imaginary part of the Fourier transform of the imaginary part of the light field to be restored, F(f x ) is the intensity of the Fourier surface light field collected;

[0058] In step S224, the initial phase for:

[0059]

[0060] Among them, |f(x n ,y m )| is the amplitude of the collected image light field, R{F I}、I{F I} are the coordinate values ​​of the two dimensions of the intersection coordinates.

[0061] The analytical constraint in step S27 is obtained by the first diffraction plane light field to be corrected after Shannon interpolation in step S271. To achieve this, the first diffraction plane light field to be corrected is input after Shannon interpolation for:

[0062]

[0063] Where h(x0,y0) is the inverse Fourier transform of the finite aperture of the system;

[0064] The analytical constraint in step S27' is obtained by the second diffraction plane light field to be corrected after Shannon interpolation in step S271'. To achieve the second diffraction plane light field to be corrected after Shannon interpolation for:

[0065]

[0066] in, is the Fourier transform of a system with finite illumination or finite field of view.

[0067] The diffraction transformation between the first diffraction plane and the second diffraction plane is a Fresnel transformation, and the formula for calculating the diffraction transformation between the first diffraction plane and the second diffraction plane is:

[0068]

[0069] Among them, U zi (x′,y′) and U i (x, y) is the light field of the first diffraction plane and the second diffraction plane, z is the distance between the two planes, is the Fourier transform, k is the wave number, f x is the frequency in the x direction, f y is the frequency in the y direction, and λ is the wavelength.

[0070] The first diffraction plane and the second diffraction plane are two diffraction planes that are in one or more combinations of a near Fresnel zone, a Fresnel zone, and a far field.

[0071] In step S25, after applying the amplitude constraint and spatial spectrum constraint of the second diffraction plane light field, the input second diffraction plane light field of the i+1th iteration is for:

[0072]

[0073] in, is the input second diffraction plane light field of the i-th iteration; is the output second diffraction plane light field of the i-th iteration, is the negative feedback term of the second diffraction plane, β and β1 are iteration constants, and f max is the cutoff frequency of the calibrated pupil function;

[0074] In step S25', after applying the amplitude constraint and the finite illumination constraint of the first diffraction plane light field, the input first diffraction plane light field f of the i+1th iteration i+1 (x n ,y m )for:

[0075]

[0076] Among them, f i (x n ,y m ) is the input first diffraction plane light field of the i-th iteration, f′ i (x n ,y m ) is the output first diffraction plane light field of the i-th iteration, f di (x n ,y m ) is the negative feedback term of the first diffraction plane, β and β1 are iteration constants, and L is the cutoff radius of the calibrated illumination function or field distribution function.

[0077] Preferably, calibrating the pupil function specifically includes:

[0078] A point light source is placed in front of the optically limited imaging system, and a two-dimensional array detector is placed on the first diffraction plane to detect the image of the point light source. A point spread function is obtained from the spatial distribution of the intensity of the point light source image, and a pupil function of the optically limited imaging system is obtained by Fourier transforming the point spread function; or a point light source is placed in front of the optically limited imaging system, and a two-dimensional array detector is placed on the second diffraction plane to detect the spatial spectrum distribution of the point light source in the second diffraction plane, and the distribution is directly used as the calibrated pupil function of the optically limited imaging system;

[0079] The calibrating illumination function specifically includes: when the optically limited imaging system is subjected to limited illumination smaller than the field of view, removing the imaging target, placing a two-dimensional array detector on the first diffraction plane to detect the light intensity distribution function of the illumination light as the illumination function;

[0080] The field of view distribution function is calibrated, specifically including: under a limited field of view, irradiating the optically limited imaging system with parallel light larger than the field of view, placing a two-dimensional array detector on the first diffraction plane to detect the light intensity distribution function of the illumination light as the field of view distribution function.

[0081] Preferably, calibrating the pupil function further comprises: determining the noise intensity according to the intensity distribution of the pupil function, and determining the cutoff frequency of the calibrated pupil function according to the noise intensity, so that within the range of the cutoff frequency of the calibrated pupil function, the pupil function is a detection value, and outside the range of the cutoff frequency of the calibrated pupil function, the pupil function is 0;

[0082] When calibrating the illumination function or the field of view distribution function, the noise intensity is determined according to the signal intensity distribution of the illumination function or the field of view distribution function, and the cutoff radius of the illumination function or the field of view distribution function is set according to the noise intensity, so that the illumination function or the field of view distribution function signal within the cutoff radius is the detection value, and the illumination function or the field of view distribution function signal outside the cutoff radius is 0.

[0083] The phase recovery method based on the optically constrained imaging system and analytical constraints of the present invention not only utilizes the constraints of the optically constrained imaging system to solve the defects of the existing iterative phase recovery method that the uniqueness depends on the object prior and multiple exposures, but also utilizes the analytical constraints to solve the problems of poor reconstruction accuracy and slow convergence speed of the traditional iterative projection algorithm. It realizes high-precision phase reconstruction without losing the spatial bandwidth product of the coaxial optical system and without the need for multiple exposures, and has a 100% success rate.

[0084] Specifically, compared with single-shot exposure technologies such as CDI, the phase recovery method based on the optically restricted imaging system and analytical constraints of the present invention uses the limitations of the optical system itself as algorithm constraints, without the need for object priors. The constraints on the intensity of the collected real image plane and the Fourier plane eliminate the interference of trivial solutions to traditional iterative algorithms, solving the problem that they are prone to fall into local minima.

[0085] Compared with multiple exposure technologies such as FPM, the phase recovery method of the present invention only needs to collect intensity information of the real image plane and the Fourier plane to complete the recovery without redundant measurement.

[0086] Compared with the traditional iterative projection algorithm, the phase recovery method of the present invention can achieve high-precision reconstruction and has a 100% success rate.

[0087] Compared with solving the KK relationship constraint, the phase recovery method of the present invention does not require that the information of the Fourier surface exists only in the half plane. Therefore, compared with the off-axis holographic system, the system space-bandwidth product can be effectively improved, and compared with the synthetic aperture technology, a single exposure can be used for reconstruction. BRIEF DESCRIPTION OF THE DRAWINGS

[0088] Figure 1 This is a typical optical path diagram of a camera system.

[0089] Figure 2 This is a typical optical path diagram of a microscope system.

[0090] Figure 3 This is a typical optical path diagram of an optically constrained imaging system.

[0091] Figure 4-Figure 6 is a flow chart of a phase recovery method based on an optically constrained imaging system and analytical constraints according to an embodiment of the present invention, wherein Figure 4 The overall flow chart is shown. Figure 5 shows the phase retrieval imaging method for strongly scattering objects, Figure 6 A phase retrieval imaging method for weakly scattering objects is shown.

[0092] Figure 7 Graph of the truncation effect for an optically confined system.

[0093] Figure 8 When the object is a strong phase object, the real image plane intensity image and the Fourier plane intensity image are collected, which shows the truncation effect of the pupil function.

[0094] Fig. 9 When the object is a weak phase object, the real image plane intensity image and the Fourier plane intensity image are collected, which shows the truncation effect of the limited field of view.

[0095] Fig.10 The initialization results using analytical constraints and the comparison with other methods.

[0096] Fig.11 It is a light field complex amplitude distribution diagram and a phase image reconstructed by the phase recovery method based on the optically restricted imaging system and analytical constraints of the present invention.

[0097] Fig.12 It is a comparison diagram of the phase recovery method based on the optically restricted imaging system and analytical constraints of the present invention and the traditional projection iteration algorithm. DETAILED DESCRIPTION

[0098] The following describes the implementation of the present invention through specific embodiments, and those skilled in the art can easily understand other advantages and effects of the present invention from the contents disclosed in this specification. The present invention can also be implemented or applied through other different specific implementations, and the details in this specification can also be modified or changed in various ways based on different viewpoints and applications without departing from the spirit of the present invention.

[0099] It should be noted that the following embodiments and features therein may be combined with each other if there is no conflict.

[0100] like Figure 4-Figure 6 Shown is a flow chart of a phase recovery method based on an optically constrained imaging system and analytical constraints according to the first embodiment of the present application, which can be applicable to any band of visible light and to an optically constrained imaging system capable of collecting the image plane intensity distribution and Fourier plane intensity distribution of the imaging target.

[0101] Among them, existing optical imaging systems are all limited imaging systems. Figure 3 As shown, the optically constrained imaging system capable of collecting the image plane intensity distribution and the Fourier plane intensity distribution of the imaging target includes at least the following components: a light source 101, a first imaging lens 102, a beam splitter 103, a second imaging lens 104, a second two-dimensional array detector 105 (for obtaining the Fourier plane intensity, i.e., the second diffraction plane intensity), and a first two-dimensional array detector 106 (for obtaining the image plane intensity, i.e., the first diffraction plane intensity) which receives light from the first imaging lens 102 and is reflected by the beam splitter 103.

[0102] like Figure 4-Figure 6 As shown, according to the first embodiment of the present invention, the phase recovery method based on the optically constrained imaging system and analytical constraints specifically includes:

[0103] Step S1: Collect the image plane intensity I(x n ,y m ) and Fourier surface intensity Therefore, the Fourier surface intensity can be used as the Fourier surface spatial spectrum, and whether the imaging target is a strong scattering object or a weak scattering object can be determined according to whether the Fourier surface spatial spectrum of the imaging target is truncated.

[0104] That is, the definition of strong scattering objects and weak scattering objects is based on whether the Fourier surface space spectrum is truncated by the pupil function. The Fourier surface space spectrum of a strong scattering object is truncated. The specific truncation effect is as follows: Figure 7 The Fourier space spectrum of a weakly scattering object is not truncated.

[0105] In this embodiment, for a strong scattering object, the image plane intensity I(x n ,y m ) and Fourier surface intensity like Figure 8 As shown in the figure, for weakly scattering objects, the image plane intensity I(x n ,y m ) and Fourier surface intensity like Fig. 9 shown.

[0106] In an optically constrained imaging system, the light scattered by the imaging target is hard-truncated in the Fourier plane spatial spectrum or the imaging field of view on the image plane due to the limited numerical aperture of the optical system.

[0107] When the Fourier surface spatial spectrum is further limited due to the numerical aperture of the optical imaging system, the finite numerical aperture of the optical system corresponding to the limited Fourier surface spatial spectrum includes the pupil and / or numerical aperture size limitation of the imaging lens, and the limited Fourier surface spatial spectrum causes hard truncation of the Fourier surface spatial spectrum. The collected data is suitable for the phase recovery imaging method of strongly scattering objects.

[0108] When the imaging field of the image plane is more limited due to limited illumination (or limited field of view), the limited numerical aperture of the optical system corresponding to the limited illumination includes the size, direction and shape restrictions of the light source, and the limited numerical aperture of the optical system corresponding to the limited field of view includes the pupil size restriction of the imaging lens and / or the imaging field of view angle restriction. Limited illumination or limited field of view causes hard truncation of the imaging field of the image plane, and the collected data is suitable for the phase recovery imaging method of weakly scattering objects.

[0109] In summary, the limited numerical aperture of the optical system includes, but is not limited to, the size, direction and shape limitations of the light source, the aperture limitations of the objective lens (if there is an objective lens), the pupil, numerical aperture and imaging field angle limitations of the imaging lens, the size limitations of the detector's photosensitive area, the size limitations of the rectangular or circular aperture (if there is an aperture), the probe size limitations of the fiber optic detector (if there is a fiber optic detector), and so on.

[0110] Step S2: Determine whether the Fourier plane spatial spectrum of the imaging target is truncated. If so, it indicates that it is a strong scattering object, and the phase recovery imaging method for strong scattering objects is executed accordingly. Otherwise, it indicates that it is a weak scattering object, and the phase recovery imaging method for weak scattering objects is executed.

[0111] like Figure 5 As shown, the phase recovery imaging method of a strongly scattering object specifically includes:

[0112] Step S21: Based on the acquired image plane intensity I(x n ,y m ) and Fourier surface intensity (like Figure 8 As shown), the amplitude of the collected image light field |f(x n ,y m )| and the amplitude of the Fourier surface light field collected

[0113] Step S22: Obtain the initial phase of the image plane light field through the Hilbert transform characteristics of the light field analysis signal

[0114] The step S22 specifically includes:

[0115] Step S221: The amplitude of the collected image plane light field |f(x n ,y m )| Perform Fourier transform to obtain the Fourier transform F of the real part of the light field to be restored R ;

[0116] That is, the Fourier transform F of the real part of the light field to be restored R for:

[0117]

[0118] Step S222: Establish the Fourier transform F of the real part of the light field to be restored R The constraint equations for the real and imaginary parts of ;

[0119] Among them, the Fourier transform F of the real part of the light field to be restored is taken R The real part is R{F R}, take the Fourier transform F of the real part of the light field to be restored R The imaginary part is I{F R}, the Fourier transform F of the real part of the light field to be restored R The constraint equations for the real and imaginary parts of are:

[0120]

[0121] Among them, I{F R} is the imaginary part of the Fourier transform of the real part of the light field to be restored, R{F I} is the real part of the Fourier transform of the imaginary part of the light field to be restored, I{F R} is the real part of the Fourier transform of the real part of the light field to be restored, I{F I} is the imaginary part of the Fourier transform of the imaginary part of the light field to be restored, F(f x ) is the collected Fourier surface light field intensity.

[0122] Because I{F R}、I{F R}、F(f x )、F(-f x ) It is known that the Fourier transform F of the real part of the light field to be restored is R The constraint equations for the real and imaginary parts of can be used as two circle equations, whose centers are (I{F R},-R{F R})、(-I{F R},-R{F R}), the corresponding circle radii are

[0123] Step S223: Based on the Fourier transform F of the real part of the light field to be restored R The constraint equations of the real and imaginary parts of the intersection point are solved to obtain the coordinates (R{F I},I{F I}).

[0124] Among them, the two circle equations are combined to solve the intersection coordinates (R{F I},I{F I}).

[0125] Step S224: According to the intersection coordinates (R{F I},I{F I}), calculate the initial phase

[0126] Among them, the initial phase for:

[0127]

[0128] Among them, |f(x n ,y m )| is the amplitude of the collected image light field, R{F I}、I{F I} are the coordinate values ​​of the two dimensions of the intersection coordinates.

[0129] Step S23: Set the initial phase of the image plane light field The amplitude of the image light field collected |f(x n ,y m )| multiply to get the input image light field f1(x n ,y m ); At the same time, the input Fourier surface light field of the first iteration Initialize;

[0130] At this time, the current iteration number i=1; that is, the input image plane light field f1(xn ,y m )for:

[0131]

[0132] Among them, |f(x n ,y m )| is the amplitude of the collected image light field, is the initial phase of the image light field.

[0133] Step S24: The input image light field f of the i-th iteration i (x n ,y m ) performs Fraunhofer diffraction transform (i.e. standard two-dimensional Fourier transform) to obtain the output Fourier surface light field of the i-th iteration

[0134] That is, the output Fourier surface light field of the i-th iteration is for:

[0135]

[0136] Among them, f i (x n ,y m ) is the input image light field of the i-th iteration, is the Fourier transform.

[0137] Step S25: Using the feedback mechanism to output the Fourier surface light field of the i-th iteration Apply the amplitude constraint and spatial spectrum restriction constraint of the Fourier surface light field to obtain the input Fourier surface light field of the i+1th iteration

[0138] After applying the amplitude constraint and spatial spectrum constraint of the Fourier surface light field, the input Fourier surface light field of the i+1th iteration is for:

[0139]

[0140] in, is the input Fourier surface light field of the i-th iteration, which comes from the previous iteration. The input Fourier surface light field of the i-th iteration when i=1 comes from the initialization above; is the output Fourier surface light field of the i-th iteration, is the negative feedback term of the Fourier surface, β and β1 are iteration constants ranging from 0.5 to 1, and f max is the cutoff frequency of the calibrated pupil function. The calibration method of the pupil function is detailed below.

[0141] It represents the negative feedback of the finite Fourier surface intensity signal so that the signal outside the cutoff frequency tends to 0 during the iteration process; Represents the negative feedback of the intensity of the collected Fourier surface light field so that it tends to the amplitude of the collected Fourier surface light field during the iteration process

[0142] Among them, the Fourier surface negative feedback term Defined as:

[0143]

[0144] in, is the amplitude of the collected Fourier surface light field, is the output Fourier surface light field of the i-th iteration, is the amplitude of the output Fourier surface light field of the i-th iteration, is the input Fourier surface light field of the i-th iteration, is the amplitude of the input Fourier surface light field at the i-th iteration.

[0145] Step S26: Input Fourier surface light field of the i+1th iteration Perform the reverse Fraunhofer diffraction transformation to obtain the output image light field f′ of the i+1th iteration i+1 (x n ,y m );Then the output image light field f′ of the i+1th iteration i+1 (x n ,y m ) is projected onto the collected image plane light field |f(x n ,y m )|, and obtain the input image plane light field to be corrected for the i+1th iteration

[0146] That is to say, the output image light field f′ of the i+1th iteration i+1 (x n ,y m )for:

[0147]

[0148] The input image light field to be corrected at the i+1th iteration for:

[0149]

[0150] Among them, |f(x n ,y m)| is the amplitude of the collected image light field.

[0151] Step S27: The input image plane light field to be corrected in the i+1th iteration Apply analytical constraints, which are used to make the amplitude of the collected image plane light field contain phase information, and obtain the input image plane light field f of the i+1th iteration i+1 (x n ,y m ).

[0152] In this embodiment, the analytical constraint adopts Shannon interpolation theorem, so the step S27 specifically includes:

[0153] Step S271: The input image plane light field to be corrected in the i+1th iteration Perform Shannon interpolation to obtain the input image plane light field to be corrected after Shannon interpolation

[0154] Therefore, the analytical constraint in step S27 is obtained by the first diffraction plane light field to be corrected after Shannon interpolation in step S271. To achieve, the input image plane light field to be corrected after Shannon interpolation for:

[0155]

[0156] in, is the input image light field to be corrected after Shannon interpolation, |f(x n ,y m )| is the amplitude of the collected image light field, is the input image light field to be corrected at the i+1th iteration, h is the Shannon interpolation function, x o ,y o is the integration variable.

[0157] Step S272: Determine whether the number of Shannon interpolation k reaches the total number of Shannon interpolation N′. If so, the input image plane light field to be corrected after Shannon interpolation is As the input image light field f of the i+1th iteration i+1 (x n ,y m ); otherwise, the input image plane light field to be corrected after Shannon interpolation As the input image light field to be corrected for the i+1th iteration And return to step S271 to continue Shannon interpolation.

[0158] The total number of Shannon interpolation times N′ is usually 10.

[0159] Step S28: add 1 to the current iteration number i as the current iteration number i of the new round, and return to step S24 to perform the next round of iteration; until the iteration end condition is met, the iteration ends at this time, and the input image plane light field f of the i+1th iteration is i+1 (x n ,y m ) as the reconstructed image plane light field. Thus, phase recovery is completed, such as Fig.11 shown.

[0160] In this embodiment, the iteration end condition includes: error function ΣΣ||U i ′|-|U||<ε, and the number of iterations i> the iteration number threshold N, at least one of the two is true. Where ε is the noise level of the two-dimensional array detector, and the iteration number threshold N is usually 3000.

[0161] That is to say, if neither of them holds, it is necessary to continue the next round of iteration; if one or both of them hold, the iteration ends, and the input image light field f i+1 (x n ,y m ) is to reconstruct the image plane light field, that is, to complete the phase recovery.

[0162] like Figure 6 As shown, the phase recovery imaging method of a weakly scattering object specifically includes:

[0163] Step S21': Based on the acquired image plane intensity I(x n ,y m ) and Fourier surface intensity (like Fig. 9 As shown), the amplitude of the collected image light field |f(x n ,y m )| and the amplitude of the Fourier surface light field collected

[0164] Step S22': Obtain the initial phase of the image plane light field through the Hilbert transform characteristics of the light field analysis signal

[0165] The step S22' specifically includes:

[0166] Step S221': The amplitude of the collected image plane light field |f(x n ,y m )| Perform Fourier transform to obtain the Fourier transform F of the real part of the light field to be restored R ;

[0167] That is, the Fourier transform F of the real part of the light field to be restored R for:

[0168]

[0169] Step S222': Establish the Fourier transform F of the real part of the light field to be restored R The constraint equations for the real and imaginary parts of ;

[0170] Among them, the Fourier transform F of the real part of the light field to be restored is taken R The real part is R{F R}, take the Fourier transform F of the real part of the light field to be restored R The imaginary part is I{F R}, the Fourier transform F of the real part of the light field to be restored R The constraint equations for the real and imaginary parts of are:

[0171]

[0172] Among them, I{F R} is the imaginary part of the Fourier transform of the real part of the light field to be restored, R{F I} is the real part of the Fourier transform of the imaginary part of the light field to be restored, I{F R} is the real part of the Fourier transform of the real part of the light field to be restored, I{F I} is the imaginary part of the Fourier transform of the imaginary part of the light field to be restored, F(f x ) is the collected Fourier surface light field intensity.

[0173] Because I{F R}、I{F R}、F(f x )、F(-f x ) It is known that the Fourier transform F of the real part of the light field to be restored is R The constraint equations for the real and imaginary parts of can be used as two circle equations, whose centers are (I{F R},-R{F R})、(-I{F R},-R{F R}), the corresponding circle radii are

[0174] Step S223': According to the Fourier transform F of the real part of the light field to be restored R The constraint equations of the real and imaginary parts are solved to obtain the intersection coordinates (R{F I},I{F I});

[0175] Step S223': According to the intersection coordinates (R{F I},I{F I}) Calculate the initial phase

[0176] Initial Phase for:

[0177]

[0178] Among them, |f(x n ,y m )| is the amplitude of the image light field.

[0179] Therefore, the value of the initial phase obtained in step S22' is exactly the same as the value of the initial phase obtained in step S22 above.

[0180] Step S23': Set the initial phase of the image plane light field The amplitude of the Fourier surface light field Multiply to get the input Fourier surface light field of the first iteration At the same time, the input image light field f1(x n ,y m ) to initialize;

[0181] At this time, the number of iterations i = 1: that is, the input Fourier surface light field of the first iteration for:

[0182]

[0183] Among them, |f(x n ,y m )| is the amplitude of the collected Fourier surface light field, is the initial phase of the image light field.

[0184] Step S24': Input Fourier surface light field of the i-th iteration Perform the reverse Fraunhofer diffraction transformation to obtain the output image light field f of the i-th iteration i ′(x n ,y m );

[0185] That is to say, the output image light field f′ of the i-th iteration i (x n ,y m )for:

[0186]

[0187] in, is the input Fourier surface light field of the i-th iteration, is the inverse Fourier transform.

[0188] Step S25': Using the feedback mechanism to calculate the output image light field f of the i-th iteration i ′(xn ,y m ) imposes the amplitude constraint and limited illumination (or limited field of view) constraint on the image plane light field to obtain the input image plane light field f of the i+1th iteration i+1 (x n ,y m );

[0189] After applying the amplitude constraint and finite illumination constraint of the image plane light field, the input image plane light field f of the i+1th iteration is i+1 (x n ,y m )for:

[0190]

[0191] Among them, f i (x n ,y m ) is the input image plane light field of the i-th iteration, which comes from the previous iteration. The input image plane light field of the i-th iteration when i=1 comes from the initialization above. i ′(x n ,y m ) is the output image light field of the i-th iteration, f di (x n ,y m ) is the image plane negative feedback term, β and β1 are iteration constants, ranging from 0.5 to 1, and L is the cutoff radius of the calibrated illumination function or field of view distribution function. The calibration method is detailed below.

[0192] Therefore, f i (x n ,y m )-βf′ i (x n ,y m ) represents the negative feedback of limited illumination (or limited field of view) so that its signal outside the illumination tends to 0 during the iteration process; f′ i (x n ,y m )+β1f di (x n ,y m ) represents the negative feedback of the collected image plane light field amplitude so that it tends to the collected image plane light field amplitude |f(x n ,y m )|.

[0193] Among them, the image plane negative feedback term f di (x n ,y m ) is defined as:

[0194]

[0195] Among them, |f(x n ,y m )| is the amplitude of the collected image light field, f i ′(x n ,y m ) is the output image light field of the i-th iteration, |f i ′(x n ,y m )| is the amplitude of the output image light field of the i-th iteration, f i (x n ,y m ) is the input image light field of the i-th iteration, |f i (x n ,y m )| is the amplitude of the input image light field of the i-th iteration.

[0196] Step S26': the input image light field f of the i+1th iteration i+1 (x n ,y m ) performs Fraunhofer diffraction transformation to obtain the output Fourier surface light field of the i+1th iteration Then the output Fourier surface light field of the i+1th iteration is The amplitude of the light field projected onto the collected Fourier surface Get the input Fourier surface light field to be corrected for the i+1th iteration

[0197] That is, the output Fourier surface light field of the i+1th iteration is for:

[0198]

[0199] Input Fourier surface light field to be corrected at the i+1th iteration for:

[0200]

[0201] in, is the amplitude of the collected Fourier surface light field.

[0202] Step S27': for the input image plane light field to be corrected in the i+1th iteration Apply analytical constraints, which are used to make the amplitude of the collected Fourier surface light field contain phase information, and obtain the input Fourier surface light field of the i+1th iteration

[0203] In this embodiment, the analytical constraint adopts Shannon's interpolation theorem, so the step S27' specifically includes:

[0204] Step S271': the input Fourier surface light field to be corrected in the i+1th iteration Perform Shannon interpolation to obtain the Fourier surface light field to be corrected after Shannon interpolation

[0205] Therefore, the analytical constraint in step S27' is obtained by the second diffraction plane light field to be corrected after Shannon interpolation in step S271'. To achieve, the Fourier surface light field to be corrected after Shannon interpolation for:

[0206]

[0207] in, is the input Fourier surface light field to be corrected after Shannon interpolation, is the amplitude of the collected Fourier surface light field, is the input Fourier surface light field to be corrected at the i+1th iteration, h is the Shannon interpolation function, is the integration variable.

[0208] Step S272': determine whether the number of Shannon interpolation k reaches the total number of Shannon interpolation N'. If so, input the Fourier surface light field to be corrected after Shannon interpolation As the input Fourier surface light field of the i+1th iteration Otherwise, the Shannon interpolated correction is input into the Fourier surface light field As the input Fourier surface light field to be corrected for the i+1th iteration And return to step S271 to continue Shannon interpolation.

[0209] The total number of Shannon interpolation times N′ is usually 10.

[0210] Step S28': add 1 to the current iteration number i as the current iteration number i of the new round, and return to step S24' to perform the next round of iteration; until the iteration end condition is met, the iteration ends at this time, and the input image plane light field f of the i+1th iteration is i+1 (x n ,y m ) is used as the reconstructed image light field, thus completing the phase recovery.

[0211] In this embodiment, the iteration end condition includes: error function ∑∑||U′ i|-|U||<ε, and the number of iterations i> the iteration number threshold N, at least one of the two is true. Among them, ε is the noise level of the two-dimensional array detector, and the iteration number threshold N is usually 3000.

[0212] That is to say, if neither of them holds, it is necessary to continue the next round of iteration; if one or both of them hold, the iteration ends, and the input image light field f i+1 (x n ,y m ) is to reconstruct the image plane light field, that is, to complete the phase recovery.

[0213] Among them, calibrating the pupil function specifically includes:

[0214] A point light source is placed in front of the optically constrained imaging system, and a two-dimensional array detector is placed on the image plane to detect the image of the point light source. The point spread function is obtained from the spatial distribution of the intensity of the point light source image, and the pupil function of the optically constrained imaging system is obtained by Fourier transform. Alternatively, a point light source is placed in front of the optically constrained imaging system, and a two-dimensional array detector is placed on the Fourier plane to detect the Fourier plane spatial spectrum distribution of the point light source, and directly used as the calibrated pupil function of the optically constrained imaging system.

[0215] In particular, when the imaging system is a microscope system, fluorescent beads are placed as point light sources within the imaging depth of field of the front focal plane of the imaging system, a two-dimensional array detector is placed on the image plane to detect the image of the point light source, and a point spread function is obtained from the intensity spatial distribution of the point light source image, which is Fourier transformed to obtain the pupil function of the optically confined imaging system; or fluorescent beads are placed as point light sources within the imaging depth of field of the front focal plane of the imaging system, a two-dimensional array detector is placed on the Fourier plane to detect the Fourier plane spectrum distribution of the point light source, and directly used as the pupil function of the optically confined imaging system.

[0216] The method further comprises: determining the noise intensity of the pupil function according to the intensity distribution of the pupil function, and determining the cutoff frequency f of the calibrated pupil function according to the noise intensity. R , so that the cutoff frequency f of the calibrated pupil function R In the range of the pupil function, the detection value is at the cutoff frequency f R Outside the range of , the pupil function is 0.

[0217] At the cutoff frequency f of the calibrated pupil function R In the range, the detection intensity of the pupil function is greater than the noise intensity, so when the spatial frequency of the Fourier surface is less than the cutoff frequency f of the limited system R When f x 2 +f y 2 <f R2 When , the pupil function is the detection value. If the spatial frequency of the Fourier surface is greater than the cutoff frequency of the restricted system, the pupil function is 0.

[0218] It should be noted that when obtaining the calibrated pupil function, the detected object, ie, the imaging target, needs to be removed, because the pupil function is irrelevant to the imaging target.

[0219] Among them, calibrating the illumination function specifically includes: when the optically limited imaging system is subjected to limited illumination smaller than the field of view (that is, the size of the illumination light cannot completely cover the pupil of the imaging lens), removing the imaging target and placing a two-dimensional array detector on the image plane to detect the intensity distribution function of the illumination light as the illumination function.

[0220] The field of view distribution function is calibrated, specifically including: under a limited field of view (i.e., the pupil size of the imaging lens is limited, and the field of view angle of the imaging system is limited, resulting in a limited field of view), irradiating the optically limited imaging system with parallel light larger than the field of view, and placing a two-dimensional array detector on the image plane to detect the light intensity distribution function of the illumination light as the field of view distribution function.

[0221] When calibrating the illumination function or the field of view distribution function, the noise intensity is determined according to the signal intensity distribution of the illumination function or the field of view distribution function, and the cutoff radius L of the illumination function or the field of view distribution function is set according to the noise intensity, so that the illumination function or the field of view distribution function signal within the cutoff radius is the detection value, and the illumination function or the field of view distribution function signal outside the cutoff radius is 0.

[0222] Within the cut-off radius, the detection intensity is greater than the noise intensity. Therefore, when the spatial coordinate range is less than the cut-off radius L, that is, x 2 +y 2 <L 2 , the illumination function is the detection value; otherwise, outside the cutoff radius, the detection intensity is less than the noise intensity, so the illumination function is 0.

[0223] In other embodiments of the present invention, the phase recovery method based on the optically constrained imaging system and analytical constraints of the present invention is not limited to the diffraction transformation between the real image plane and the Fourier plane and its iteration. The diffraction transformation between the real image plane and the Fourier plane can also be replaced by a diffraction transformation between two diffraction planes in one or more combinations of the near Fresnel zone, the Fresnel zone, and the far field.

[0224] For example, in the second embodiment of the present invention, the Fraunhofer diffraction transform (i.e., the standard two-dimensional Fourier transform) between the real image plane and the Fourier plane is replaced by the diffraction transform between two diffraction planes in the Fresnel zone. Therefore, the other steps of this embodiment are the same as the steps of the phase recovery method based on the optically restricted imaging system and analytical constraints of the first embodiment of the present application. The only difference is that in all the steps, the image plane is replaced by the first diffraction plane in the Fresnel zone, the Fourier plane is replaced by the second diffraction plane in the Fresnel zone, the phase recovery imaging method for strongly scattering objects is replaced by the first type of phase recovery imaging method, and the phase recovery method for weakly scattering objects is replaced by the second type of phase recovery imaging method.

[0225] Correspondingly, in step S24, step S26, step S24', and step S26', the Fraunhofer diffraction transform is replaced by a diffraction transform calculation between the first diffraction plane and the second diffraction plane in the Fresnel zone (i.e., Fresnel transform), and the reverse Fraunhofer diffraction transform is replaced by a diffraction transform calculation between the second diffraction plane and the first diffraction plane in the Fresnel zone.

[0226] Therefore, the formula for calculating the diffraction transformation between the first diffraction plane and the second diffraction plane is as follows:

[0227]

[0228] Among them, U zi (x′,y′) and U i (x, y) is the light field of the first diffraction plane and the second diffraction plane, z is the distance between the two planes, is the Fourier transform, k is the wave number, f x is the frequency in the x direction, f y is the frequency in the y direction, and λ is the wavelength.

[0229] It should be noted that in any embodiment of the present invention (such as the first diffraction plane and the second diffraction plane are in one or more combinations of the near Fresnel zone, the Fresnel zone, and the far field), the truncation of the first diffraction plane in the present invention corresponds to limited illumination or limited field of view, and the truncation of the second diffraction plane corresponds to limited spatial spectrum.

[0230] In the present invention, the results of initialization using analytical constraints and comparison with other methods are Fig.10 As shown by Fig.10 It can be seen that the initialization error is lower when using analytical constraints for initialization. The comparison diagram of the phase recovery method based on optically constrained imaging system and analytical constraints of the present invention and the traditional projection iteration algorithm is shown in FIG. Fig.12 As shown by Fig.12 It can be seen that the phase recovery method of the present invention has a lower reconstruction error.

[0231] The above is only a preferred embodiment of the present invention, and is not intended to limit the scope of the present invention. The above embodiments of the present invention can also be modified in various ways. All simple, equivalent changes and modifications made according to the claims and the description of the present invention fall within the scope of protection of the claims of the present invention. The contents not described in detail in the present invention are all conventional technical contents.

Claims

1. A phase retrieval method based on an optically constrained imaging system and analytical constraints, characterized in that: include: Step S1: collecting a first diffraction plane intensity and a second diffraction plane intensity of an imaging target through the optically limited imaging system; Step S2: if the second diffraction plane intensity is truncated, then the first type of phase recovery imaging method is executed; otherwise, the second type of phase recovery imaging method is executed; The first type of phase recovery imaging method specifically includes: Step S21: Based on the collected first diffraction plane intensity I(x n ,y m ) and the second diffraction plane intensity The amplitude of the light field in the first diffraction plane collected is obtained by square root calculation |f(x n ,y m )| and the amplitude of the light field collected in the second diffraction plane Step S22: Obtaining the initial phase of the first diffraction plane light field through the Hilbert transform characteristics of the light field analysis signal Step S23: The initial phase of the first diffraction plane light field The amplitude of the light field in the first diffraction plane collected |f(x n ,y m )| multiply to get the first diffraction plane light field f1(x n ,y m ); At the same time, for the second diffraction plane light field of the first iteration Initialize; Step S24: Input the first diffraction plane light field f of the i-th iteration i (x n ,y m ) performs the diffraction transformation calculation between the first diffraction plane and the second diffraction plane to obtain the output second diffraction plane light field of the i-th iteration Step S25: Using the feedback mechanism to output the second diffraction plane light field of the i-th iteration Apply the amplitude constraint and spatial spectrum restriction constraint of the second diffraction plane light field to obtain the input second diffraction plane light field of the i+1th iteration Step S26: Input the second diffraction plane light field of the (i+1)th iteration Perform the diffraction transformation calculation between the second diffraction plane and the first diffraction plane to obtain the output first diffraction plane light field f of the i+1th iteration i ′ +1 (x n ,y m ); then the output first diffraction plane light field f of the i+1th iteration is i ′ +1 (x n ,y m ) is projected onto the amplitude of the light field of the first diffraction plane collected |f(x n ,y m )|, and obtain the first diffraction plane light field of the input to be corrected in the i+1th iteration Step S27: input the first diffraction plane light field to be corrected in the (i+1)th iteration Applying analytical constraints, we can obtain the input first diffraction plane light field f of the i+1th iteration i+1 (x n ,y m ), including: Step S271: input the first diffraction plane light field to be corrected in the (i+1)th iteration Perform Shannon interpolation to obtain the first diffraction plane light field to be corrected after Shannon interpolation Step S272: Determine whether the number of Shannon interpolation k reaches the total number of Shannon interpolation N ′ If so, then the Shannon interpolated input to be corrected is input into the first diffraction plane light field As the input of the i+1th iteration, the first diffraction plane light field f i+1 (x n ,y m ); otherwise, the Shannon interpolated light field to be corrected is input into the first diffraction plane light field As the first diffraction plane light field to be corrected input for the i+1th iteration And return to step S271; Step S28: add 1 to the current iteration number i as the current iteration number i of the new round, and return to step S24 to perform the next round of iteration; until the iteration end condition is met, the iteration ends at this time, and the input first diffraction plane light field f of the i+1th iteration is i+1 (x n ,y m ) as the reconstructed first diffraction plane light field; The second type of phase recovery imaging method specifically includes: Step S21': Based on the collected first diffraction plane intensity I(x n ,y m ) and the second diffraction plane intensity The amplitude of the light field in the first diffraction plane collected is obtained by square root calculation |f(x n ,y m )| and the amplitude of the light field collected in the second diffraction plane Step S22': Obtain the initial phase of the first diffraction plane light field through the Hilbert transform characteristics of the light field analysis signal Step S23': Set the initial phase of the first diffraction plane light field The amplitude of the light field collected in the second diffraction plane Multiply to get the input second diffraction plane light field of the first iteration At the same time, the first diffraction plane light field f1(x n ,y m ) to initialize; Step S24': Input the second diffraction plane light field of the i-th iteration Perform the diffraction transformation calculation between the second diffraction plane and the first diffraction plane to obtain the output first diffraction plane light field f of the i-th iteration i ′ (x n ,y m ); Step S25': Using the feedback mechanism to output the first diffraction plane light field f of the i-th iteration i ′ (x n ,y m ) applies the amplitude constraint of the first diffraction plane light field and the constraint of limited illumination or limited field of view to obtain the input first diffraction plane light field f of the i+1th iteration i+1 (x n ,y m ); Step S26': Input the first diffraction plane light field f of the i+1th iteration i+1 (x n ,y m ) performs the diffraction transformation calculation between the first diffraction plane and the second diffraction plane to obtain the output second diffraction plane light field of the i+1th iteration Then the output second diffraction plane light field of the i+1th iteration is The amplitude of the light field projected onto the second diffraction plane collected Get the second diffraction plane light field to be corrected input for the i+1th iteration Step S27': input the second diffraction plane light field to be corrected for the i+1th iteration Apply analytical constraints to obtain the input second diffraction plane light field of the i+1th iteration Specifically include: Step S271': input the second diffraction plane light field to be corrected in the i+1th iteration Perform Shannon interpolation to obtain the second diffraction plane light field to be corrected after Shannon interpolation Step S272': Determine whether the number of Shannon interpolation k reaches the total number of Shannon interpolation N ′ If so, then the Shannon interpolated input to be corrected is input into the second diffraction plane light field As the input second diffraction plane light field of the i+1th iteration Otherwise, the Shannon interpolated correction is input into the second diffraction plane light field As the second diffraction plane light field to be corrected input for the i+1th iteration And return to step S271; Step S28': add 1 to the current iteration number i as the current iteration number i of the new round, and return to step S24' to perform the next round of iteration; until the iteration end condition is met, the iteration ends at this time, and the input first diffraction plane light field f of the i+1th iteration is i+1 (x n ,y m ) as the reconstructed first diffraction plane light field.

2. The phase retrieval method based on an optically constrained imaging system and analytical constraints according to claim 1, characterized in that: The first diffraction plane is an image plane, and the second diffraction plane is a Fourier plane; In the step S1, the second diffraction plane intensity distribution is a Fourier surface space spectrum, and whether the imaging target is a strong scattering object or a weak scattering object is determined according to whether the Fourier surface space spectrum is truncated; In step S2, the first type of phase recovery imaging method is a phase recovery imaging method for strongly scattering objects, and the second type of phase recovery imaging method is a phase recovery method for weakly scattering objects.

3. The phase retrieval method based on an optically constrained imaging system and analytical constraints according to claim 2, characterized in that: In step S24, step S26, step S24' and step S26', the diffraction transformation between the first diffraction plane and the second diffraction plane is calculated as a Fraunhofer diffraction transformation, and the diffraction transformation between the second diffraction plane and the first diffraction plane is calculated as an inverse Fraunhofer diffraction transformation.

4. The phase recovery method based on optically constrained imaging system and analytical constraints according to claim 2, characterized in that: In an optically constrained imaging system, the light scattered by the imaging target causes a hard truncation of the Fourier plane spatial spectrum or the imaging field of view on the image plane due to the limited numerical aperture of the optical system. The limited numerical aperture of the optical system includes the size limitation of the light source, the aperture limitation of the objective lens, the pupil, numerical aperture and imaging field angle of the imaging lens, the size limitation of the detector's photosensitive area, the size limitation of the rectangular or circular aperture, and the probe size limitation of the fiber optic detector; The limited numerical aperture of the optical system corresponding to the limited Fourier surface spatial spectrum includes the pupil and / or numerical aperture size limitation of the imaging lens; the limited Fourier surface spatial spectrum causes the Fourier surface spatial spectrum to be hard truncated; Limited illumination corresponds to the size limitation of the illumination source on the object, and limited field of view corresponds to the limited numerical aperture of the optical system, including the pupil size limitation of the imaging lens and / or the limitation of the imaging field of view angle; limited illumination or limited field of view leads to hard truncation of the imaging field of view on the image plane.

5. The phase retrieval method based on an optically constrained imaging system and analytical constraints according to claim 1, characterized in that: The step S22 specifically includes: Step S221: The amplitude of the first diffraction plane light field collected |f(x n ,y m )| Perform Fourier transform to obtain the Fourier transform F of the real part of the light field to be restored R ; Step S222: Establish the Fourier transform F of the real part of the light field to be restored R Constraint equations for the real and imaginary parts of ; Step S223: Based on the Fourier transform F of the real part of the light field to be restored R The constraint equations of the real and imaginary parts of the intersection point are solved to obtain the coordinates (R{F I },I{F I }); Step S224: According to the intersection coordinates (R{F I },I{F I }), calculate the initial phase Wherein, in the step S222, the Fourier transform F of the real part of the light field to be restored R The constraint equations for the real and imaginary parts of are: Among them, I{F R } is the imaginary part of the Fourier transform of the real part of the light field to be restored, R{F I } is the real part of the Fourier transform of the imaginary part of the light field to be restored, I{F R } is the real part of the Fourier transform of the real part of the light field to be restored, I{F I } is the imaginary part of the Fourier transform of the imaginary part of the light field to be restored, F(f x ) is the intensity of the Fourier surface light field collected; In step S224, the initial phase for: Among them, |f(x n ,y m )| is the amplitude of the collected image light field, R{F I }、I{F I } are the coordinate values ​​of the two dimensions of the intersection coordinates.

6. The phase retrieval method based on optically constrained imaging system and analytical constraints according to claim 1, characterized in that: The analytical constraint in step S27 is obtained by the first diffraction plane light field to be corrected after Shannon interpolation in step S271. To achieve this, the first diffraction plane light field to be corrected is input after Shannon interpolation for: Where h(x0,y0) is the inverse Fourier transform of the finite aperture of the system; The analytical constraint in step S27' is obtained by the second diffraction plane light field to be corrected after Shannon interpolation in step S271'. To achieve the second diffraction plane light field to be corrected after Shannon interpolation for: in, is the Fourier transform of a system with finite illumination or finite field of view.

7. The phase recovery method based on optically constrained imaging system and analytical constraints according to claim 1, characterized in that: The diffraction transformation between the first diffraction plane and the second diffraction plane is a Fresnel transformation, and the formula for calculating the diffraction transformation between the first diffraction plane and the second diffraction plane is: Among them, U zi (x ′ ,y ′ ) and U i (x, y) is the light field of the first diffraction plane and the second diffraction plane, z is the distance between the two planes, is the Fourier transform, k is the wave number, f x is the frequency in the x direction, f y is the frequency in the y direction, and λ is the wavelength.

8. The phase recovery method based on an optically confined imaging system according to claim 1, characterized in that: The first diffraction plane and the second diffraction plane are two diffraction planes that are in one or more combinations of a near Fresnel zone, a Fresnel zone, and a far field.

9. The phase recovery method based on optically constrained imaging system and analytical constraints according to claim 1, characterized in that: In step S25, after the amplitude constraint and spatial spectrum constraint of the second diffraction plane light field, the second diffraction plane light field of the i+1th iteration for: in, is the input second diffraction plane light field of the i-th iteration; is the output second diffraction plane light field of the i-th iteration, is the negative feedback term of the second diffraction plane, β and β1 are iteration constants, and f max is the cutoff frequency of the calibrated pupil function; In step S25', after applying the amplitude constraint and the finite illumination constraint of the first diffraction plane light field, the input first diffraction plane light field f of the i+1th iteration i+1 (x n ,y m )for: Among them, f i (x n ,y m ) is the input first diffraction plane light field of the i-th iteration, f i ′ (x n ,y m ) is the output first diffraction plane light field of the i-th iteration, f di (x n ,y m ) is the negative feedback term of the first diffraction plane, β and β1 are iteration constants, and L is the cutoff radius of the calibrated illumination function or field distribution function.

10. The phase recovery method based on optically constrained imaging system and analytical constraints according to claim 9, characterized in that: In step S25, calibrating the pupil function specifically includes: A point light source is placed in front of the optically limited imaging system, and a two-dimensional array detector is placed on the first diffraction plane to detect the image of the point light source. A point spread function is obtained from the spatial distribution of the intensity of the point light source image, and a pupil function of the optically limited imaging system is obtained by Fourier transforming the point spread function; or a point light source is placed in front of the optically limited imaging system, and a two-dimensional array detector is placed on the second diffraction plane to detect the spatial spectrum distribution of the point light source in the second diffraction plane, and the distribution is directly used as the calibrated pupil function of the optically limited imaging system; In the step S25', calibrating the illumination function specifically includes: when the optically limited imaging system is subjected to limited illumination smaller than the field of view, removing the imaging target, placing a two-dimensional array detector on the first diffraction plane to detect the light intensity distribution function of the illumination light as the illumination function; The field of view distribution function is calibrated, specifically comprising: under a limited field of view, irradiating the optically limited imaging system with parallel light larger than the field of view, placing a two-dimensional array detector on a first diffraction plane to detect a light intensity distribution function of the illumination light as the field of view distribution function; Calibrating the pupil function also includes: determining the noise intensity according to the intensity distribution of the pupil function, and determining the cutoff frequency of the calibrated pupil function according to the noise intensity, so that within the range of the cutoff frequency of the calibrated pupil function, the pupil function is a detection value, and outside the range of the cutoff frequency of the calibrated pupil function, the pupil function is 0; When calibrating the illumination function or the field of view distribution function, the noise intensity is determined according to the signal intensity distribution of the illumination function or the field of view distribution function, and the cutoff radius of the illumination function or the field of view distribution function is set according to the noise intensity, so that the illumination function or the field of view distribution function signal within the cutoff radius is the detection value, and the illumination function or the field of view distribution function signal outside the cutoff radius is 0.

Images