High-precision coded aperture imaging method

Abstract

The application discloses a high-precision coded aperture imaging method and belongs to the technical field of lensless imaging. The method uses image spectrum information of one image to cross-iterate optical transfer function and spectrum information of an object, repeatedly decodes and compares coded intensity images, captures and superimposes detail information of each iteration to recover an image, finally converges to obtain high-precision optical transfer function, and uses the high-precision optical transfer function and image spectrum obtained by Fourier transform of an image obtained by a CCD to obtain high-precision recovered object information. The method has good high-frequency detail information recovery effect of an image, fast convergence speed and high reconstruction efficiency.

Description

TECHNICAL FIELD

[0001] The application relates to a high-precision coded aperture imaging method and belongs to the technical field of lensless imaging. BACKGROUND

[0002] As a lensless imaging technology, coded aperture imaging does not need to use a lens, but encodes a light beam through a specially designed coded board, and belongs to the category of lensless imaging technology. Since the 1960s, coded aperture imaging has been widely used in short-wavelength (x-ray, beta-ray, etc.) imaging, medical and astronomical fields due to its high light transmittance and signal-to-noise ratio. However, the traditional coded aperture imaging has low resolution accuracy and relatively poor imaging effect, which restricts its development and application in the field of lensless imaging.

[0003] In the traditional coded aperture imaging, the Wiener filtering method is usually used to reconstruct the information of an object from a blurred image collected, and the accuracy of the reconstruction method depends largely on the intensity point spread function of the coded aperture imaging device. However, the intensity point spread function is not easy to accurately obtain in actual experiments, which is the main reason for limiting the high-precision imaging of coded aperture imaging. At present, the methods for obtaining the intensity point spread function include direct projection method, pinhole photography method and point spread function formula derivation considering scalar diffraction effect. Considering the actual measurement, the pinhole photography method is commonly used in these methods. In actual measurement, a micro-hole is placed on the object plane, and the blurred spot obtained on the imaging device by the illumination beam through the coded aperture imaging device can be regarded as the intensity point spread function image. However, due to the limitation of the size of the micro-hole (usually microns) and the limited light flux, the collected intensity point spread function image usually has low accuracy, which has limitations for the recovery of image details when used for image reconstruction, and the effect of the reconstructed image is not ideal.

[0004] In order to improve the problem of low image reconstruction accuracy caused by inaccurate intensity point spread function measurement, the existing method proposes a method for improving the accuracy of the reconstructed image by collecting multiple coded board patterns to perform multiple iterations of object spectral information, that is, first randomly guess the spectral information of the object, in the light path of the coded aperture imaging system, by changing the pattern of the coded board each time (for example, a certain small angle (such as 7 degrees) can be rotated in each experiment) to collect the image obtained each time, the Fourier transform is performed to obtain the image spectrum information, the intensity point spread function of the coded aperture imaging light path of this time, and the optical transfer function after the Fourier transform. Multiply the optical transfer function obtained this time and the spectral information of the object obtained last time (take the spectral information of the object randomly guessed at the beginning as the spectral information of the object obtained for the 0th time) to obtain the image spectrum information, compare it with the spectral information of the blurred image collected by the experiment, perform Wiener filtering deconvolution on the residual image spectrum, and correct the spectral information of the object obtained last time to obtain the spectral information of the object this time. According to the above steps, multiple iterations are performed until the spectral information of the object almost does not need to be corrected (that is, it meets certain accuracy requirements), and then the final spectral information of the object can be obtained, and the spectral information of the image can be obtained by multiplying the optical transfer function of the coded aperture imaging system. After inverse Fourier transform, the image plane information can be obtained. However, the accuracy of the reconstructed image of the method depends largely on the intensity point spread function under the coded board pattern each time, and if the intensity point spread function collected in multiple iterations is inaccurate, the error of the recovered object will be large, resulting in that the recovered image is still not clear enough. That is, the method cannot guarantee the accuracy of the transfer function by changing the pattern of the coded board multiple times, and thus the imaging accuracy cannot be guaranteed.

[0005] In addition, the resolution accuracy of the imaging receiving device (usually a CCD or CMOS camera) is in the micrometer level, and the resolution accuracy of the coded aperture imaging is usually at least one order of magnitude lower than that, so the low resolution accuracy of the coded aperture imaging itself becomes an important factor limiting its wide application. SUMMARY

[0006] In order to solve the problem of poor imaging effect caused by the difficulty in accurately measuring the intensity point spread function in the current coded aperture imaging, the application provides a high-precision coded aperture imaging method. The method does not need to change the pattern of the coded board, so that the optical transfer function of the coded aperture imaging system does not change, thereby guaranteeing the accuracy of the optical transfer function of the system. The application cross-iterates the optical transfer function and the spectral information of the object by using the image spectrum information of an image, finally converges to obtain a high-precision optical transfer function, and obtains high-precision recovered object information by using the high-precision optical transfer function and the image spectrum obtained by Fourier transform of the image acquired by the CCD, thereby improving the clarity of the recovered image.

[0007] A high-precision coded aperture imaging method, the method is based on a coded aperture imaging system, the method comprises:

[0008] Step S1: obtaining a blurred image Image of an object to be imaged and an intensity point spread function g(x, y);

[0009] Step S2: according to the blurred image Image of the object to be imaged and the intensity point spread function g(x, y) obtained, calculating the optical transfer function OTF(0) in the frequency domain and the spectrum f_Image of the image;

[0010] Step S3: according to the optical transfer function OTF(0) in the frequency domain and the spectrum f_Image of the image, obtaining the spectrum f_Object(0) of the object to be imaged by Wiener filtering method;

[0011] Step S4: taking the optical transfer function OTF(0) in the frequency domain as the initial value of the optical transfer function OTF distribution, and taking f_Object(0) as the initial value of the spectrum of the object to be imaged, cross-iterating the spectrum of the optical transfer function and the object to be imaged, so that the spectrum information f_Object of the object to be imaged and the optical transfer function OTF converge;

[0012] Step S5: using the iteratively updated spectrum information f_Object(i) of the object to be imaged and the optical transfer function OTF(i) to obtain the updated spectrum f_Image(i) of the image;

[0013] Step S6: performing inverse Fourier transform on the updated spectrum f_Image(i) of the image to obtain the image Image(i), and calculating the image residual ΔI(i) between the image Image(i) and the blurred image Image of the object to be imaged obtained in step S1;

[0014] Step S7: setting a critical condition ε0 according to the imaging accuracy, when the image residual ΔI(i) is less than the critical condition, calculating the accurate object Object according to the corresponding spectrum information f_Object(i) of the object to be imaged and the optical transfer function OTF(i).

[0015] Optionally, the coded aperture imaging system comprises a light source, an aperture, a coding plate and an imaging device arranged in sequence; wherein the coding plate is an arbitrary shape coding plate; the diameter of the aperture is in the order of microns.

[0016] Optionally, the step S1 comprises:

[0017] Step a1, turning on the light source, so that the light source passes through the aperture and the coding plate, and the imaging device obtains the corresponding image, encodes the obtained image as the intensity point spread function g(x, y);

[0018] Step a2: placing the object to be imaged at the stop instead of the stop so that the center of the object to be imaged is on the optical axis; turning on the light source, the imaging device obtaining a blurred image Image of the object to be imaged, and obtaining an encoded intensity image Image(x, y) of the object to be imaged according to the blurred image Image.

[0019] Optionally, the step S2 comprises:

[0020] obtaining the optical transfer function OTF(0) and the spectrum f_Image of the image in the frequency domain by Fourier transform

[0021]

[0022] where F{} is a Fourier transform operator.

[0023] Optionally, the step S3 comprises:

[0024] calculating the spectrum f_Object(0) of the object to be imaged according to the following formula:

[0025]

[0026] where ξ is a set minimum positive number; OTF(0) * is the conjugate complex of the initial value of the optical transfer function OTF(0).

[0027] Optionally, the step S4 comprises:

[0028] Step b1: taking OTF(0) as the initial value of the optical transfer function OTF distribution and f_Object(0) as the initial value of the object to be imaged distribution in the frequency domain;

[0029] Step b2: the image spectrum residual Δf_Image(i) after the i-th iteration is represented as:

[0030] Δf_Image(i) = f_Image - f_Object(i-1) * OTF(i-1)

[0031] where i = 1, 2, …; f_Image represents the imaging spectrum information, f_Object(i-1) represents the object spectrum information, and OTF(i-1) represents the optical transfer function (i-1)th step b3: alternately updating the object spectrum information f_Object(i) and the OTF(i) information so that the object spectrum information f_Object and the optical transfer function OTF converge;

[0032]

[0033] Wherein, the spectrum information f_Object of the object and the optical transfer function OTF convergence condition is that the MSE of the image residual obtained by the two is less than the set accuracy requirement.

[0034] Optionally, the step S5 comprises:

[0035] After the i-th iteration, the updated image spectrum f_Image(i) is obtained by using the updated object spectrum f_Object(i) and the optical transfer function OTF(i) as follows:

[0036] f_Image(i) = f_Object(i) * OTF(i).

[0037] Optionally, the step S6 comprises:

[0038] The inverse Fourier transform is performed on the updated image spectrum f_Image(i) to obtain the image Image(i)

[0039] Image(i) = |F -1 {f_Image(i)}|

[0040] Wherein, F -1 {} is the inverse Fourier transform operator;

[0041] The image residual ΔI(i) between the image Image(i) and the actually collected blurred intensity image Image(x, y) is calculated, and the MSE is calculated:

[0042] ΔI(i) = Image(i) - Image

[0043]

[0044] Wherein, M represents the number of pixel columns of the image, and N represents the number of pixel rows of the image.

[0045] Optionally, the critical condition ε0 in the step S7 represents the difference degree of the image surface information obtained by multiplying the optical transfer functions and the object spectrum information obtained by using adjacent two iterations and the image surface information actually captured by the CCD;

[0046] When the image residual ΔI(i) is not less than the critical condition ε0, the cross iteration is continued until the image residual ΔI(i) is less than the critical condition ε0.

[0047] Optionally, the imaging device adopts a CCD camera.

[0048] The present application has the following beneficial effects:

[0049] The method of the present application uses the image frequency spectrum information to cross-iterate the optical transfer function and the spectrum information of the object, repeatedly decodes and compares the coded intensity image, captures the detail information of each iteration, superimposes the image recovery, finally converges to obtain the high-precision optical transfer function, and uses the high-precision optical transfer function and the image frequency spectrum obtained by the Fourier transform of the image obtained by the CCD to obtain the high-precision recovered object information, the image high-frequency detail information recovery effect is good, the convergence speed is fast, and the reconstruction efficiency is high. BRIEF DESCRIPTION OF DRAWINGS

[0050] In order to more clearly illustrate the technical solutions in the embodiments of the present application, the drawings needed to be used in the embodiment description will be briefly introduced. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can be obtained by those skilled in the art without creative labor.

[0051] Figure 1 is the system structure diagram for coding aperture imaging reconstruction in the present application.

[0052] Figure 2 is the total flow chart of the reconstruction method of the high-precision coding aperture imaging system disclosed in the present application.

[0053] Figure 3 is the simulation comparison diagram of the optical transfer function recovered by the method of the present application, the actually obtained optical transfer function and the theoretically obtained optical transfer function.

[0054] Figure 4 is the iteration error curve diagram of the root mean square error (MSE) between the image obtained by iteration and updating by the method of the present application and the actually collected image.

[0055] Figure 5 is the comparison diagram of the object recovered by the method of the present application and the object itself. DETAILED DESCRIPTION

[0056] In order to make the purpose, technical solutions and advantages of the present application more clear, the embodiments of the present application will be further described in detail below with reference to the drawings.

[0057] Embodiment one:

[0058] The present embodiment provides a high-precision coding aperture imaging method, the method comprising:

[0059] Step S1: obtaining the blurred image Image of the object to be imaged and the intensity point spread function g(x, y);

[0060] Step S2: According to the obtained blurred image Image of the object to be imaged and the intensity point spread function g(x, y), the optical transfer function OTF(0) in the frequency domain and the spectrum f_Image of the image are calculated;

[0061] Step S3: According to the optical transfer function OTF(0) in the frequency domain and the spectrum f_Image of the image, the spectrum f_Object(0) of the object to be imaged is obtained by Wiener filtering method;

[0062] Step S4: The optical transfer function OTF(0) in the frequency domain is taken as the initial value of the optical transfer function OTF distribution, and the spectrum f_Object(0) of the object to be imaged is taken as the initial value of the object to be imaged in the frequency domain. Cross iteration is performed on the spectrum of the optical transfer function and the object to be imaged, so that the spectrum information f_Object of the object to be imaged and the optical transfer function OTF converge;

[0063] Step S5: The updated spectrum f_Image(i) of the image is obtained by using the updated spectrum information f_Object(i) of the object to be imaged and the optical transfer function OTF(i);

[0064] Step S6: The inverse Fourier transform is performed on the updated spectrum f_Image(i) of the image to obtain the image Image(i), and the image residual ΔI(i) between the image Image(i) and the blurred image Image of the object to be imaged obtained in step S1 is calculated;

[0065] Step S7: According to the imaging accuracy, the critical condition is set, when the image residual ΔI(i) is less than the critical condition ε0, the accurate object Object is calculated according to the corresponding spectrum information f_Object(i) of the object to be imaged and the optical transfer function OTF(i).

[0066] Embodiment two:

[0067] The embodiment provides a high-precision coded aperture imaging method, which is realized based on a high-precision coded aperture imaging system, as shown in the figure, the system comprises a light source, an aperture, a coded board and a CCD camera arranged in sequence, wherein the shape of the coded board is not limited and can be a round hole board, a square hole board, a single ring hole board, a multi-ring hole board and the like. When imaging, the object to be imaged is placed at the position of the aperture. Figure 1

[0068] As shown in the figure, the high-precision coded aperture imaging method comprises the steps of obtaining the initial value of the object to be imaged and the optical transfer function, and iteratively updating the object to be imaged and the optical transfer function; specifically, comprising: Figure 2

[0069] (1) the step of obtaining the initial value of the object to be imaged and the optical transfer function comprises: ​​

[0070] Step a1, obtaining intensity point spread function g(x, y): light source, diaphragm, encoding board and CCD camera are placed in turn, wherein the diaphragm is placed on the object plane and coaxial with the optical axis, and the diameter of the diaphragm is micron level; turn on the light source, so that the light source passes through the diaphragm and the encoding board, and the CCD camera collects the intensity point spread function g(x, y) after the illumination light passes through the diaphragm and the encoding board;

[0071] The CCD camera collects an image, and the corresponding intensity point spread function g(x, y) is obtained by encoding the image.

[0072] Step a2, obtaining the initial value of the distribution of the object to be imaged: placing the object to be measured on the object plane, replacing the diaphragm with the object to be imaged at the position of the diaphragm, so that the center of the object to be imaged is on the optical axis, and the CCD camera imaging device collects the encoded intensity image Image(x, y) after the illumination light passes through the imaging system device;

[0073] The CCD camera collects a blurred image of the object to be imaged, and the encoded intensity image Image(x, y) of the object to be imaged is obtained according to the blurred image.

[0074] Step a3: obtaining the optical transfer function OTF(0) and the spectrum f_Image in the frequency domain by Fourier transform

[0075]

[0076] Where F{} is the Fourier transform operator.

[0077] Step a4: obtaining the spectrum f_Object(0) of the object to be imaged by Wiener filtering method

[0078]

[0079] Where ξ is a small positive number set; OTF(0) * is the conjugate complex of the initial value of the optical transfer function OTF(0).

[0080] (2) The object to be imaged and the optical transfer function iterative updating step include:

[0081] Step b1: taking OTF(0) as the initial value of the optical transfer function OTF distribution, and f_Object(0) as the initial value of the distribution of the object to be imaged in the frequency domain;

[0082] Step b2: the image spectrum residual Δf_Image(i) after the i-th (i=1, 2, …) iteration is represented as

[0083] Δf_Image(i) = f_Image - f_Object(i - 1) * OTF(i - 1)

[0084] Step b3: alternately iteratively update the object's spectral information f_Object(i) and OTF(i) information, so that the object's spectral information f_Object and the optical transfer function OTF converge.

[0085]

[0086] Wherein, the object's spectral information f_Object and the optical transfer function OTF convergence condition is the MSE root mean square error of the image residual obtained by the two is less than the set accuracy requirement.

[0087] Step b4: after the i-th iteration, the updated object spectrum f_Object(i) and the optical transfer function OTF(i) are used to obtain the updated image spectrum f_Image(i)

[0088] f_Image(i) = f_Object(i) * OTF(i)

[0089] Step b5: inverse Fourier transform the updated image spectrum f_Image(i) to obtain the image Image(i)

[0090] Image(i) = |F -1 {f_Image(i)}|

[0091] Wherein F -1 {} is the inverse Fourier transform operator.

[0092] Step b6: compare the image residual ΔI(i) between the image Image(i) and the actually collected blurred intensity image Image(x, y), and calculate the root mean square error MSE.

[0093] ΔI(i) = Image(i) - Image

[0094]

[0095] Wherein, M represents the number of pixel columns of the image, and N represents the number of pixel rows of the image.

[0096] Step b7: set the critical condition ε0, when the MSE calculated by step b6 is greater than the set critical condition ε0, then return to step b2 to reiterate the loop, until the MSE calculated by step b6 is less than the set critical condition ε0.

[0097] When the MSE calculated by step b6 is less than the set critical condition ε0, it is considered that the accurate object spectrum f_Object(i) and optical transfer function OTF(i) have been obtained. The accurate object Object after the accuracy has been improved by the coded aperture iterative loop can be obtained by the inverse Fourier transform.

[0098]

[0099] It should be noted that the critical condition ε0 represents the difference between the image plane information obtained by multiplying the optical transfer function obtained from two consecutive iterations and the object's spectral information, and then performing an inverse Fourier transform, and the image plane information actually captured by the CCD. When the difference meets a certain accuracy requirement, it indicates that the required high-precision imaging has been obtained. In practical applications, this value can be determined based on the required accuracy of object information recovery; for example, it can be set to 10. -7 .

[0100] To verify the effectiveness of the method proposed in this application, the optical transfer function recovered using this method was compared with the actual optical transfer function and the theoretically obtained optical transfer function through simulation. Figure 3 As shown, Figure 3 The right side shows an enlarged OTF plot of the optical transfer function taken along y=257. It can be seen that the optical transfer function recovered using the method of this application is infinitely close to its theoretical value. Therefore, the method of this application has high accuracy. The iteration error curve is shown below. Figure 4 As shown, the root mean square error between the image obtained in each iteration and the actual acquired image decreases continuously with the number of iterations and gradually converges, indicating that the method of this application makes the recovered image more accurate.

[0101] The high-precision optical transfer function recovered based on the method of this application and the comparison image of the recovered object and the object itself are shown below. Figure 5 As shown in the figure, the object recovery effect is good, indicating that the method of this application has high precision.

[0102] Some steps in the embodiments of the present invention can be implemented using software, and the corresponding software program can be stored in a readable storage medium, such as an optical disc or a hard disk.

[0103] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A high-precision coded-aperture imaging method, characterized by, The method is implemented based on an encoding aperture imaging system, and the method comprises the following steps: Step S1: obtaining a blurred image Image of an object to be imaged and an intensity point spread function g(x, y); Step S2: calculating an optical transfer function OTF(0) in a frequency domain and a spectrum f_Image of an image according to the obtained blurred image Image of the object to be imaged and the intensity point spread function g(x, y); Step S3: obtaining a spectrum f_Object(0) of the object to be imaged by a Wiener filtering method according to the optical transfer function OTF(0) in the frequency domain and the spectrum f_Image of the image; Step S4: taking the optical transfer function OTF(0) as an initial value of an optical transfer function OTF distribution and taking the spectrum f_Object(0) of the object to be imaged as an initial value of a distribution in the frequency domain, and performing cross iteration on the spectrum of the optical transfer function and the object to be imaged, so that the spectrum information f_Object of the object to be imaged and the optical transfer function OTF converge; Step S5: obtaining an updated spectrum f_Image(i) of the image by using the updated spectrum information f_Object(i) of the object to be imaged and the optical transfer function OTF(i) after iteration; Step S6: performing inverse Fourier transform on the updated spectrum f_Image(i) of the image to obtain an image Image(i), and calculating an image residual error ΔI(i) between the image Image(i) and the blurred image Image of the object to be imaged obtained in step S1; Step S7: setting a critical condition ε0 according to imaging accuracy, and when the image residual error ΔI(i) is less than the critical condition ε0, calculating an accurate object Object according to the corresponding spectrum information f_Object(i) of the object to be imaged and the optical transfer function OTF(i). The step S4 comprises: Step b1: taking OTF(0) as an initial value of an optical transfer function OTF distribution and taking f_Object(0) as an initial value of a distribution of an object to be imaged in a frequency domain; Step b2: the image spectrum residual error Δf_Image(i) after the i th iteration is represented as: Δf_Image(i) = f_Image-f_Object(i-1)*OTF(i-1) Wherein, i = 1, 2, …; f_Image represents imaging spectrum information, f_Object(i-1) represents object spectrum information, and OTF(i-1) represents the (i-1) th iteration result of an optical transfer function; Step b3: alternately updating the spectrum information f_Object(i) of the object and the OTF(i) information, so that the spectrum information f_Object of the object and the optical transfer function OTF converge; Wherein, the convergence condition of the spectrum information f_Object of the object and the optical transfer function OTF is that the MSE (Mean Square Error) of the image residual error obtained by the two is less than the set accuracy requirement.

2. The method of claim 1, wherein, The encoding aperture imaging system comprises a light source, a diaphragm, an encoding plate and an imaging device arranged in sequence; wherein the encoding plate is an encoding plate of any shape; the diameter of the diaphragm is in the order of microns.

3. The method of claim 2, wherein, The step S1 comprises: Step a1, turn on the light source so that the light source passes through the diaphragm and the encoding plate, the imaging device obtains the corresponding image, and the obtained image is encoded as an intensity point spread function g(x, y); Step a2, place the object to be imaged at the diaphragm to replace the diaphragm, so that the center of the object to be imaged is on the optical axis; turn on the light source, and the imaging device obtains a blurred image Image of the object to be imaged, and obtains an encoded intensity image Image(x, y) of the object to be imaged according to the blurred image Image.

4. The method of claim 3, wherein, The step S2 comprises: Obtaining the optical transfer function OTF(0) and the spectrum f_Image of the image in the frequency domain by Fourier transform Where F{} is the Fourier transform operator.

5. The method of claim 4, wherein, The step S3 comprises: The spectrum f_Object(0) of the object to be imaged is calculated according to the following formula: where ξ is a set small positive number; OTF(0) * is the conjugate complex of the OTF(0) optical transfer function initial value.

6. The method of claim 5, wherein, The step S5 comprises: After the i-th iteration, the updated object spectrum f_Object(i) and the optical transfer function OTF(i) are used to obtain the updated spectrum f_Image(i) of the image as follows: f_Image(i) = f_Object(i) * OTF(i).

7. The method of claim 6, wherein, The step S6 comprises: The inverse Fourier transform of the updated spectrum f_Image(i) of the image is performed to obtain the image Image(i) Image(i) = |F -1 {f_Image(i)}| where F -1 is the inverse Fourier transform operator; The image residual ΔI(i) between the image Image(i) and the actually captured blurred intensity image Image(x, y) is calculated, and the mean square error MSE is calculated: ΔI(i) = Image(i) - Image Where M represents the number of pixel columns of the image, and N represents the number of pixel rows of the image.

8. The method of claim 7, wherein, The critical condition ε0 in the step S7 represents the difference degree between the image surface information obtained by multiplying the optical transfer function and the object spectrum information obtained by using two adjacent iterations and the image surface information actually captured by the CCD; When the image residual ΔI(i) is not less than the critical condition ε0, return to step b2 to continue cross iteration until the image residual ΔI(i) is less than the critical condition ε0.

9. The method of claim 8, wherein, The imaging device adopts a CCD camera.

Images