Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Performance enhancement of coherent modulation imaging in the presence of missing data

Open Access Open Access

Abstract

Coherent diffraction imaging (CDI) has become a powerful imaging modality in synchrotron x-ray imaging and electron microscopy communities. In the far-field geometry, image quality of CDI depends strongly on the performance of detector; specifically, the dynamic range, pixel size, and the absence of missing data. Coherent modulation imaging (CMI), an innovative variant of CDI, improves the algorithmic convergence by inserting a modulator upstream of the detector. Here, we explore the potential of CMI in eliminating nonideal effects of detector by modifying the modulus constraint to extrapolate the missing part of diffraction pattern. Nine folds of extrapolation in area of diffraction pattern have been shown feasible in experiment; while sixteen folds in simulation. For image quality measured by Structural Similarity (SSIM), our method shows a maximum of 32% improvement over the traditional method. Our method provides a way to alleviate the effects of beamstop, gaps between modules, limited dynamic range, and limited detector size for CMI.

© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

1. Introduction

Coherent diffraction imaging (CDI) has been widely applied to image noncrystalline nanoparticles and biological samples [16]. It requires no lens but the diffraction pattern, iterative algorithm, and support information to recover the object wave [710]. And it can offer wavelength limited resolution [11,12]. However, loose support, missing diffraction data, and noise degrade reconstructed quality, severely [7,1317]. Also, for nanometer or atomic resolution imaging with large field of view, CDI is commonly conducted in far-field geometry; for which the detector needs to have high dynamic range, a large number of pixels, and small pixel size [11,18,19]. In practice, one would encounter nonideal situations as illustrated in Fig. 1(b)–(e), namely, missing data due to the beamstop, limited detection area, missing data due to the module gaps, or finite pixel size unfulfilling the Nyquist requirement. To explain it further, the employed beamstop blocks a part of central diffraction data so, the low-frequency diffraction data is partially lost as shown in Fig. 1(b). Also, using a detector with limited detection region cannot capture the high-frequency diffraction data distributed in the margin region which is shown in Fig. 1(c). Additionally, X-ray detector with large number of pixels is usually composed of several modules as Fig. 1(d). The data transmission circuits are generally installed in between the modules, thus, the diffraction data distributed in these regions cannot be detected. Also, physical pixel with large size as shown in Fig. 1(e) can have high dynamic range but it may not fulfill the Nyquist requirement. Various techniques such as structured illumination, lateral displacements of illumination, and subpixel translating of detector have been proposed to improve resolution and image quality by increasing the diversity of measurement and recovering the undetected information [12,2027]. These techniques can provide satisfying results when tens to thousands of measurements can be taken from a stable apparatus.

 figure: Fig. 1.

Fig. 1. Geometry of far-field imaging system. (a) Modulated diffraction pattern occupies a larger area comparing with unmodulated one due to the convolution effect and strong scattering nature of modulator. ${\otimes} $ is the convolution operator, $F\{{\cdot} \}$ is the Fourier transform operator, $AS\{{\cdot} \}$ is the angular-spectrum propagation operator, P is the probe, O is the object, M is the modulator, D is the detector, ${Z_1}$ is the distance between object and modulator, ${Z_2}$ is the distance between modulator and detector. (b) The beamstop blocks diffraction data distributed in central region. (c) The diffraction data in the gap between modules cannot be detected. (d) Data distributed beyond detector is lost. (e) The size of physical pixels is larger than the virtual pixels that obeys the Nyquist requirement.

Download Full Size | PDF

Coherent modulation imaging (CMI), a single-shot imaging technique, breaks the direct Fourier transform relation between the exit wave and diffraction data by inserting a modulator between the object and detector. Optically, the unknown exit wave gets modulated with a known modulator. At the far-field diffraction plane, the diffracted object wave gets convoluted with the Fourier spectrum of the modulation function, resulting in an expanded diffraction pattern as illustrated in Fig. 1(a). Encoding the sample information over a greater number of detector pixels enhances the robustness against imperfections of detector [28,29].

In this paper, we exploit the capability of single-shot imaging of CMI and extend the current CMI algorithm to overcome the problem of missing data. This paper is organized as follows. Foundations and algorithm of our method are first explained in Section 2. Results of simulations and proof-of-concept experiments are presented in Section 3 and Section 4. In Section 5 the performance of our algorithm is discussed. The overall summary is given in Section 6.

2. Methods

2.1 Foundations

We limited our discussion to setup shown in Fig. 1(a) where the modulator was placed downstream of the object, closely. And the detector was located in the far field diffraction region. For such geometry, the information of object wave is intertwined with the modulator function via convolution and distributed over a larger region than the case without modulator. Under such arrangement, it is feasible to accomplish reconstruction when only a portion of diffraction data is available. Additionally, the modulus of the Fourier spectrum of weak scattering sample concentrate on central region. And the modulator is strong scattering which leads a large diffraction pattern with specific distribution. Under the geometry mentioned above, the presence of weak sample will lead tiny modification on the diffraction pattern of modulator according to the property of convolution. With this notion, we can provide a reference for iterative algorithm to recover the missing part of diffraction pattern by making up the undetected part with the diffraction pattern of modulator.

To quantitatively evaluate the results in our simulations and experiments, we adopt the following metrics.

For simulation, we evaluate the reconstructed wavefront via invariant normalized root-mean-square (RMS) error metric [30], which is stated as follows

$$RMS = \sqrt {1 - \frac{{\sum {{|{{\mathrm{\Phi }^{\ast }}({{{\textbf r}_d}} )\mathrm{\Psi }({{{\textbf r}_d}} )} |}^2}}}{{\sum {{|{\mathrm{\Psi }({{{\textbf r}_d}} )} |}^2}{{|{\mathrm{\Phi }({{{\textbf r}_d}} )} |}^2}}}} $$

Symbol $\mathrm{\Psi }({{{\textbf r}_d}} )$ is the Fourier transform of exit wave reconstructed, $\mathrm{\Phi }({{{\textbf r}_d}} )$ denotes the Fourier transform of ground truth, ${\ast} $ is the complex conjugation.

Apart from RMS metric, we evaluate the recovered diffraction pattern with an error matrix which is stated as

$$\epsilon = {\sqrt {\frac{{{{\left( {|{\mathrm{\psi }({{{\textbf r}_d}} )} |- \sqrt {{\textrm{I}_\textrm{r}}({{{\textbf r}_d}} )} } \right)}^2}}}{{{\textrm{I}_\textrm{r}}({{{\textbf r}_d}} )+ \textrm{C}}}} _.}$$
where $\mathrm{\psi }({{{\textbf r}_d}} )$ is the recovered exit wave at the detector plane, ${\textrm{I}_\textrm{r}}({{{\textbf r}_d}} )$ is the real diffraction pattern, $\textrm{C}$ is a constant that prevent infinity after division.

To provide visual evaluation of reconstructed image, we also employ Structural Similarity (SSIM) [31] with definition of

$$SSIM = {\frac{{({2{\mathrm{\mu }_\mathrm{\psi }}{\mathrm{\mu }_\phi } + {\textrm{C}_1}} )({2{\mathrm{\sigma }_{\mathrm{\psi }\phi }} + {\textrm{C}_2}} )}}{{({\mathrm{\mu }_\mathrm{\psi }^2 + \mathrm{\mu }_\phi^2 + {\textrm{C}_1}} )({\mathrm{\sigma }_\mathrm{\psi }^2 + \mathrm{\sigma }_\phi^2 + {\textrm{C}_2}} )}}_.}$$
where ${\mathrm{\mu }_\mathrm{\psi }},{\mathrm{\mu }_\phi }$ are the mean of recovered and reference amplitude at the support plane, ${\mathrm{\sigma }_\mathrm{\psi }},{\; }{\mathrm{\sigma }_\phi }$ are the standard deviation of recovered and reference amplitude, ${\mathrm{\sigma }_{\mathrm{\psi }\phi }}$ is the covariance of recovered and reference amplitude, ${\textrm{C}_1},{\textrm{C}_2}$ are constants that prevent infinity after division.

In addition, we define two variables to quantify the degree of extrapolation from the perspective of area and energy.

First, ${\eta _A}$ measures the ratio of total area of calculation window with respect to the size of camera in the unit of pixel, and it is stated as

$${\eta _A} = {\frac{{\sum {\textrm{A}_\textrm{c}}({{{\textbf r}_d}} )}}{{\sum \textrm{A}({{{\textbf r}_{\boldsymbol d}}} )}}_.}$$

In this formula, ${\textrm{A}_\textrm{c}}({{{\textbf r}_{\boldsymbol d}}} )$ denotes the binary matrix of calculation window, $\textrm{A}({{{\textbf r}_{\boldsymbol d}}} )$ denotes the camera array which is a binary matrix.

Second, metric of extrapolation energy ratio ${\eta _E}$ measures the amount of energy of measured relative to total of diffraction patten, written as

$${\eta _E} = {\frac{{\sum {\textrm{I}_\textrm{r}}({{{\textbf r}_d}} )\textrm{A}({{{\textbf r}_d}} )}}{{\sum {\textrm{I}_\textrm{r}}({{{\textbf r}_d}} ){\textrm{A}_\textrm{c}}({{{\textbf r}_d}} )}}_.}$$

2.2 Algorithm

The proposed algorithm is modified on the conventional CMI algorithm. Proposed method finishes reconstruction via single measurement, also, attempts to enhance the performance by extrapolating the diffraction data beyond the detection region of detector. The flowchart of the performance enhanced CMI algorithm is shown in Fig. 2. And implementation details are stated as follows:

  • 0) Initiate the combined diffraction pattern as input. The combined diffraction pattern composes of the measured part and calculated part. The measured part is the diffraction data that camera captured. The rest of calculation window is filled with the calculated diffraction pattern. The calculated diffraction pattern is obtained through forward propagation of a known probe or energy normalized support when a flat probe is assumed.
  • 1) Apply the support constraint.
    $${\mathrm{\psi }^\textrm{k}}({{{\textbf r}_{\boldsymbol o}}} )= {\mathrm{\hat{\psi }}^{\textrm{k} - 1}}({{{\textbf r}_{\boldsymbol o}}} )\textrm{S}$$
    where ${\mathrm{\psi }^\textrm{k}}({{{\textbf r}_{\boldsymbol o}}} )$ is the wavefront after support constraint in the current iteration, ${\mathrm{\hat{\psi }}^{\textrm{k} - 1}}({{{\textbf r}_{\boldsymbol o}}} )$ denotes the object wave after the modified modulus constraint in the $\textrm{k} - 1$ iteration, $\textrm{S}$ denotes support function.
  • 2) Propagate the exit wave to the front side of modulator.
    $${\mathrm{\psi }^\textrm{k}}({{{\textbf r}_m}} )= AS\{ {\mathrm{\psi }^\textrm{k}}({{{\textbf r}_{\boldsymbol o}}} )\} $$

    Symbol ${\mathrm{\psi }^\textrm{k}}({{{\textbf r}_m}} )$ is the wavefront at the front side of the modulator, $AS\{{\cdot} \}$ is the angular-spectrum propagation method.

  • 3) Apply modulation to get the wavefront at the rear side of modulator.
    $${\mathrm{\psi }^\textrm{k}}({{{\textbf r}_M}} )= {\mathrm{\psi }^\textrm{k}}({{{\textbf r}_m}} )\textrm{M}({{{\textbf r}_m}} )$$
    where ${\mathrm{\psi }^\textrm{k}}({{{\textbf r}_M}} )$ is the wavefront at the rear side of the modulator, $\textrm{M}({{{\textbf r}_m}} )$ is the transmission function of modulator.
  • 4) Propagate the modulated wavefront to the detector plane with Fourier transform.
    $${\mathrm{\psi }^\textrm{k}}({{{\textbf r}_{\boldsymbol d}}} )= F\{{2{\mathrm{\psi }^\textrm{k}}({{{\textbf r}_M}} )- {\mathrm{\psi }^{\textrm{k} - 1}}({{{\textbf r}_M}} )} \}$$

    Here, ${\mathrm{\psi }^\textrm{k}}({{{\textbf r}_{\boldsymbol d}}} )$ is the wavefront at detector plane, ${\mathrm{\psi }^{\textrm{k} - 1}}({{{\textbf r}_M}} )$ is the wavefront at the rear side of modulator in the $\textrm{k} - 1$ iteration. And ${\mathrm{\psi }^{\textrm{k} - 1}}({{{\textbf r}_M}} )$ is set as ${\mathrm{\psi }^\textrm{k}}({{{\textbf r}_M}} )$ when $\textrm{k}$ equals to the integral multiple of 10 to escape from stagnation.

  • 5) Apply the modified modulus constraint.
    $${\mathrm{\hat{\psi }}^\textrm{k}}({{{\textbf r}_d}} )= \textrm{A}({{{\textbf r}_{\boldsymbol d}}} )\frac{{{\mathrm{\psi }^\textrm{k}}({{{\textbf r}_d}} )\sqrt {\textrm{I}({{{\textbf r}_d}} )} }}{{|{{\mathrm{\psi }^\textrm{k}}({{{\textbf r}_d}} )} |}} + \gamma ({{\textrm{A}_\textrm{c}}({{{\textbf r}_{\boldsymbol d}}} )- \textrm{A}({{{\textbf r}_{\boldsymbol d}}} )} )\frac{{{\mathrm{\psi }^\textrm{k}}({{{\textbf r}_d}} )\sqrt {\textrm{I}({{{\textbf r}_d}} )} }}{{|{{\mathrm{\psi }^\textrm{k}}({{{\textbf r}_d}} )} |}}$$
    where ${\mathrm{\hat{\psi }}^\textrm{k}}({{{\textbf r}_d}} )$ is the wavefront after the modified modulus constraint, $\textrm{I}({{{\textbf r}_d}} )$ is the combined diffraction pattern, $\gamma $ is a weighted coefficient that controls the calculated diffraction intensity outside the camera. Modified modulus constraint reinforces the retrieved diffraction pattern outside the aperture of camera to be consistent with the calculated diffraction pattern, for the first few iterations. The calculated part of diffraction pattern provides a reference for iterative algorithm to make up the part of missing data outside the detection region. For the part within the aperture of camera, the reconstructed diffraction amplitude is enforced be consistent with the measured diffraction data for whole loop to ensure the algorithm can search the accurate object wave.

    In our algorithm, the value of $\gamma $ and how it changes along with iterations are key factors. With too few iterations or not employing the calculated diffraction pattern, iterative algorithm usually yields unsatisfying result. In practice, we set $\mathrm{\gamma } = 1$ for the first 5–30 iterations to eliminate the potential solution space beyond the aperture of detector; then, we set $\mathrm{\gamma } = 0$ to ensure the consistency between the measured one and recovered one within the detection area.

  • 6) Propagate the revised wavefront back to the rear side of modulator with inverse Fourier transform.
    $${\mathrm{\hat{\psi }}^\textrm{k}}({{{\textbf r}_M}} )= {F^{ - 1}}\{{{{\mathrm{\hat{\psi }}}^\textrm{k}}({{{\textbf r}_d}} )} \}$$
  • 7) Undo modulation to update the wavefront at front side of modulator.
    $${\mathrm{\hat{\psi }}^\textrm{k}}({{{\textbf r}_m}} )= {\mathrm{\psi }^\textrm{k}}({{{\textbf r}_m}} )+ \frac{{\textrm{M}{{({{{\textbf r}_m}} )}^\ast }}}{{|{\textrm{M}({{{\textbf r}_m}} )} |_{\textrm{max}}^2}}({{{\mathrm{\hat{\psi }}}^\textrm{k}}({{{\textbf r}_M}} )\; \; - {\mathrm{\psi }^\textrm{k}}({{{\textbf r}_M}} )+ {\mathrm{\psi }^{\textrm{k} - 1}}({{{\textbf r}_M}} )} )$$
  • 8) Propagate the updated wavefront from modulator plane to support plane.
    $${\mathrm{\hat{\psi }}^\textrm{k}}({{{\textbf r}_{\boldsymbol o}}} )= A{S^{ - 1}}\{{{{\mathrm{\hat{\psi }}}^\textrm{k}}({{{\textbf r}_m}} )} \}$$
  • 9) Repeat steps 1-8 until the predetermined iteration number is reached.

 figure: Fig. 2.

Fig. 2. The illustration of the performance enhanced CMI. (a) The combined diffraction pattern is the input of our algorithm. (b) Flowchart of proposed algorithm.

Download Full Size | PDF

Our code is written in MATLAB using parallel computing toolbox.

3. Simulations

In this section, simulations were performed to validate the feasibility of proposed method. A $1mm$ pinhole was illuminated by a flat $632.8nm$ laser beam with a total of $2 \times {10^9}$ photons. The complex-valued object composed of a positive USAF resolution target image in amplitude and a cat image in phase is shown in Fig. 3(a) and 3(b). Distance among object, modulator and detector were set as ${Z_1} = 30mm$ and ${Z_2} = 30mm$, respectively. Phase modulator was a gaussian filtered random matrix with range of $0\sim 2\mathrm{\pi }$. We simulated a camera with 14bits in dynamic range and $5.04\mu m$ in pixel size. A $1024 \times 1024$ diffraction pattern was first generated and then its central part was cropped out for use.

 figure: Fig. 3.

Fig. 3. (a-b) Ground truth of amplitude and phase. (c) Simulated diffraction pattern and the detection area are marked by white dotted squares. (d) Simulated diffraction pattern without object.

Download Full Size | PDF

In Fig. 4(a1), (a2), (d1) and (d2), the quality of both conventional CMI algorithm and our method is comparable when ${\eta _A} = 4$. And the result of conventional algorithm deteriorates rapidly with the increasing area of missing data, while our method can still maintain satisfying quality. In Fig. 4(c1), (c2), (f1) and (f2); conventional algorithm fails to search optimal result due to the limited size of calculation window and the lack of diffraction data when ${\eta _A} = 16$. By our method, the calculated part of combined diffraction pattern provides a reference for iterative algorithm to complete the missing part of diffraction pattern, so our method can yield satisfying result when a large portion of diffraction data is lost.

 figure: Fig. 4.

Fig. 4. Crop from center $1m{m^2}$ of reconstructed amplitude and phase when ${\eta _A} = 4,9,16$. (a-c) Reconstructed object by conventional method. (d-f) Results of our method.

Download Full Size | PDF

The linecut plots of Group −2 and $\textrm{RMS}$ curves over 200 iterations are illustrated in Fig. 5 to evaluate the performance of our method. In Fig. 5(a), linecut plot of our method fits the ground truth better than the traditional method because the size of calculation window is an important factor for reconstruction. With conventional algorithm, the calculation window depends on the size of detection area; thus, the reconstruction quality will deteriorate rapidly if the integrity of diffraction pattern cannot be guaranteed. On the contrary, a larger calculation window and the modified modulus constraint enables the proposed algorithm to make up the missing diffraction data and improve convergency as shown in Fig. 5(b). And the gap between two methods becomes evident along with the increasing area of missing diffraction data.

 figure: Fig. 5.

Fig. 5. Analysis for reconstruction results and convergency when ${\eta _A} = 4,9,16$. (a) The linecut plots of recovered amplitude by different methods. (b) Evolution of $\textrm{RMS}$ over 200 iterations with different methods.

Download Full Size | PDF

Figure 6 provides further evaluation for our method from the perspective of diffraction pattern when ${\eta _A} = 16$. Figure 6(a) and 6(b) show a visual comparison of recovered and measured diffraction pattern. Speckle distribution of both retrieved and simulated diffraction pattern shows a good agreement. The matrix $\epsilon $ is shown in Fig. 6(c) and the white dotted square marks the center detector region of $ 256 \times 256$ pixels. The ring-average of $\epsilon $ is shown in Fig. 6(d). Here, we define the criterion of successful reconstruction as $\epsilon - \textrm{ring}\; \textrm{average} < 0.07$. The dotted red and green line indicate the detector region and extrapolated part of diffraction patten.

 figure: Fig. 6.

Fig. 6. Analysis for recovered diffraction pattern when ${\eta _A} = 16$. (a-b) Measured and recovered diffraction pattern, respectively. (c)The matrix of $\epsilon $. The red dotted line and white square indicate the extent of input diffraction pattern. Green circle marks the region extrapolated, successfully. (d) Ring average of (c). The red and green line correspond to the red and green circle of (c).

Download Full Size | PDF

4. Experiment

We used a stabilized HeNe laser to generate highly coherent laser beam with center wavelength of $632.8nm$, and a flux of $2 \times {10^9}$ photons. The beam was expanded and collimated with two lenses. For sample, a positive USAF resolution target was sticked on a $1.1mm$ pinhole. The sample was $30.68mm$ away from the modulator. Modulator employed was a thin phase plate with divergence angle of ${3.5^ \circ }$. Far-field diffraction data was formed on the detector by a lens with $30mm$ focal length. The bit depth and pixel size of employed detector (CS2100) were 16bits and $5.04\mu m$, respectively. The photo of our built setup is shown in Fig. 7.

 figure: Fig. 7.

Fig. 7. Apparatus of experiment.

Download Full Size | PDF

To provide a reference, we recovered the exit wave with full diffraction pattern; and the amplitude shown in Fig. 8(a) has been divided by the known probe. In Fig. 8(b) and 8(c), the speckle distribution varies little, regardless of the absence of sample which is consistent with our theory. The dotted squares represent the extent of the detector. Obvious refinement concerning quality and resolution is presented in Fig. 9. The robustness of our method begins to show with the decreasing area of the aperture of detector. Specifically, a clear deviation between the proposed algorithm and conventional algorithm in terms of quality shows when extrapolation ratio ${\eta _A} = 4({\eta _E} = 72{\%}$). Group 6 element 2($71lp/mm$) can be resolved by our method, in contrast, we can hardly resolve Group 6 element 1($64lp/mm$) by conventional method due to the discontinuity. The distinction between two methods becomes evident when ${\eta _A} = 9({{\eta_E} = 45\%} )$, we can only recognize Group 4 element 5 $({28lp/mm} )$, but we can still recognize Group 6 element 1($64lp/mm$) from the resolution target reconstructed by our method.

 figure: Fig. 8.

Fig. 8. Experiment results. (a) Crop center $1.1m{m^2}$ from reconstructed amplitude using full diffraction pattern. And it was divided by the known probe. (b) The acutal diffraction pattern. White dotted squares indicate the extent of detector which corresponds to ${\eta _A} = 2.25,\; 4,\; 6.24,\; 9$, respectively. (c) Diffraction pattern in absence of sample.

Download Full Size | PDF

 figure: Fig. 9.

Fig. 9. Visual comparation of recovered resolution target. (a-d) Recovered amplitude of center $1.1m{m^2}$ by conventional CMI after the division of known probe. (e-h) Results of the same region by our method.

Download Full Size | PDF

Despite the visual comparation, we also calculated the $\textrm{SSIM}$ to provide a quantitative comparation, which is shown in Fig. 10. And the resolution target shown in Fig.8(a) was used as reference. Results of our method show an observable improvement over conventional method which is consistent with the visual comparation. Specifically, results of our method have maximum $32\%$ of improvement than the conventional method in terms of $\textrm{SSIM}$.

 figure: Fig. 10.

Fig. 10. Comparation of different methods by $\textrm{SSIM}$.

Download Full Size | PDF

In addition to the analysis of recovered object, evaluation for diffraction pattern is presented in Fig. 11. There is a clear consistency shown in Fig. 11(a) and 11(b) between the real diffraction pattern and recovered one within the region ($340 \times 340$ pixels) marked by white square. But, recovered diffraction pattern decays slightly slower than the actual one which may attribute to noise. The $\epsilon $ matrix in Fig. 11(c) confirms that our algorithm can not only retrieve the central part of diffracted intensity but also extrapolate the part of area surrounded by green dotted circle. Figure 11(d) shows the ring average of Fig. 11(c), and the dotted red line represents the extent of input diffracted intensity. Dotted green line represents the region where $\mathrm{\epsilon } - \textrm{ring}\; \textrm{average} < 0.07$.

 figure: Fig. 11.

Fig. 11. Analysis for recovered diffraction pattern when $ {\eta _A} = 9$. (a-b) Real and recovered diffraction pattern. Dotted square is the detector region. (c) The $\epsilon $ matrix. And the while square indicates the extent of detector. Green circle is the part extrapolated successfully. (d) Ring average of (c). And dotted red and green lines represent the size of detector and region extrapolated, respectively.

Download Full Size | PDF

5. Discussion

In this section, we discuss several key factors for our method, by simulation and experiment.

The modulation effect is critical for our method. So, we conducted simulations for continues-type and binary-type modulator to discuss the performance of proposed method when employing different types of modulators. We set the criterion for successful reconstruction as $\textrm{RMS} < 0.15$. The extrapolation ratio ${\eta _A}$ is related to modulator, so, we choose the minimum extrapolation energy ratio ${\eta _E}$ as our metric. Continues-type modulator was generated by applying a low-pass filter and gaussian filter to a random matrix. The low-pass filter was used to control the spatial frequency of modulator which was achieved by multiplying a circular mask with predetermined radius in Fourier domain. For binary-type modulator, we generated random matrix according to the predefined grain size. Here, the grain size defines the minimum scattering unit which controls the spatial frequency of phase. And the range of phase is another variable of modulator which defines the phase delay at each point. By changing these parameters, quantitative analysis for three major parameters phase frequency, phase range and modulator type, is presented in Fig. 12.

 figure: Fig. 12.

Fig. 12. Quantitative evaluation of the parameters of modulator. (a) Plot of the minimum extrapolation energy versus radius of low-pass filter for continues-type modulator. (b) Plot of the minimum extrapolation energy versus grain size for binary-type modulator.

Download Full Size | PDF

Simulation results indicate that modulator with high spatial frequency is preferred for both continues and binary type modulator. The best performance can be obtained when using continues-type modulator with phase range of $0\sim 2\mathrm{\pi \;\ }$ and high spatial frequency or binary-type modulator with phase range of $0\sim \mathrm{\pi \;\ }$ and high spatial frequency. In simulation, the binary-type modulator can have better performance compared with continuous-type modulator when grain size smaller than 4, however, it is difficult to manufacture such binary-type modulator. In practice, continues-type modulator can have better performance considering the real manufacture craft.

Figure 13 shows the results of our method when central diffraction data is missing. These results prove that the information of object and modulator is intertwined which is a consistency between theory and experiment. In Fig. 13(c), satisfying result can be obtained by our method when nearly half of energy of diffraction pattern is missing. Such a promising feature provides a flexibility for the design of detector and joint measurement.

 figure: Fig. 13.

Fig. 13. Experimental results of our method when central diffraction data is missing. (a) Real diffraction pattern. And the white squares represent the extent of central missing data. (b) The cropped amplitude recovered with full diffraction pattern after division of the know probe. (c-f) The cropped amplitude recovered when central diffraction data is missing.

Download Full Size | PDF

In experiment, we observed a tradeoff between the resolution and extrapolation ratio. A larger ${Z_1}$ allows for a larger extrapolation ratio to be used. Figure 14 shows the reconstructions versus ${Z_1}$, the distance between the sample and modulator. The degraded resolution as shown by comparing Fig. 14(f) and 14(a) is believed to be caused by the angular-spectrum propagator.

 figure: Fig. 14.

Fig. 14. Center $1.1m{m^2}$ region of recovered amplitude under different distance. (a-e) Retrieved amplitude when distance between sample and modulator ${Z_1} = 11.66mm$. (f-j) Recovered amplitude when ${Z_1} = 30.68mm$.

Download Full Size | PDF

6. Conclusion

To conclude, we develop a new algorithm for CMI to address missing data issues in experiment. In detail, we investigate the combined diffraction pattern and modified the modulus constraint for the current CMI algorithm. Both simulation and experiment validate the effectiveness of our proposed method.

The completeness of diffraction data is essential for obtaining reconstruction of good quality in CDI. Although the demand for high dynamic range of detector can be partly relaxed in CMI, the robustness of current CMI algorithm against missing diffraction data is still limited. By our proposed method, extrapolation ratio ${\eta _A} = 16$ (${\eta _E} = 27{\%}$) and ${\eta _A} = 9({{\eta_E} = 45{\%}} )$ are realized in numerical simulations and optical experiments. In addition, our method has maximum $32\%$ of improvement in $\textrm{SSIM}$ over the traditional method. Moreover, the proposed method can also work when only peripheral diffraction data are available which provide flexibility in the design of detector.

Our method has broad applications in X-ray and electron microscopy communities. Image quality and resolution can be improved by the proposed method while reducing the effect of beamstop and modules gaps in synchrotron radiation and free-electron laser imaging systems. Also, performance enhanced CMI can be combined with electron energy loss spectroscopy (EELS) to perform a joint measurement. Additionally, our method can enhance the frame rate of imaging system algorithmically by reducing the data transmitted each frame.

There are rooms for further improvement. Optimal design of the modulator so that the resultant diffraction pattern fit with the geometry of detector could lead to better performance. Such a possibility is currently under investigation.

Funding

National Natural Science Foundation of China (11775105, 12074167); Shenzhen Science and Technology Innovation Program (KQTD20170810110313773); Centers for Mechanical Engineering Research and Education at MIT and SUSTech (6941806).

Disclosures

The authors declare no conflicts of interest. This work is original and has not been published elsewhere.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. M. Nakasako, A. Kobayashi, Y. Takayama, K. Asakura, M. Oide, K. Okajima, T. Oroguchi, and M. Yamamoto, “Methods and application of coherent X-ray diffraction imaging of noncrystalline particles,” Biophys Rev 12(2), 541–567 (2020). [CrossRef]  

2. J. Miao, P. Charalambous, J. Kirz, and D. Sayre, “Extending the methodology of X-ray crystallography to allow imaging of micrometre-sized non-crystalline specimens,” Nature 400(6742), 342–344 (1999). [CrossRef]  

3. J. Miao, C.-C. Chen, C. Song, Y. Nishino, Y. Kohmura, T. Ishikawa, D. Ramunno-Johnson, T.-K. Lee, and S. H. Risbud, “Three-Dimensional GaN − Ga 2 O 3 Core Shell Structure Revealed by X-Ray Diffraction Microscopy,” Phys. Rev. Lett. 97(21), 215503 (2006). [CrossRef]  

4. J. Miao, K. O. Hodgson, T. Ishikawa, C. A. Larabell, M. A. LeGros, and Y. Nishino, “Imaging Whole Escherichia coli Bacteria by Using Single-Particle X-Ray Diffraction,” Proc. Natl. Acad. Sci. U. S. A. 100(1), 110–112 (2003). [CrossRef]  

5. D. Shapiro, P. Thibault, T. Beetz, V. Elser, M. Howells, C. Jacobsen, J. Kirz, E. Lima, H. Miao, A. M. Neiman, and D. Sayre, “Biological imaging by soft x-ray diffraction microscopy,” Proc. Natl. Acad. Sci. U.S.A. 102(43), 15343–15346 (2005). [CrossRef]  

6. C. Song, H. Jiang, A. Mancuso, B. Amirbekian, L. Peng, R. Sun, S. S. Shah, Z. H. Zhou, T. Ishikawa, and J. Miao, “Quantitative Imaging of Single, Unstained Viruses with Coherent X Rays,” Phys. Rev. Lett. 101(15), 158101 (2008). [CrossRef]  

7. R. A. Dilanian, G. J. Williams, L. W. Whitehead, D. J. Vine, A. G. Peele, E. Balaur, I. McNulty, H. M. Quiney, and K. A. Nugent, “Coherent diffractive imaging: a new statistically regularized amplitude constraint,” New J. Phys. 12(9), 093042 (2010). [CrossRef]  

8. R. W. Gerchberg and W. O. Saxton, “A Practical Algorithm for the Determination of Phase from Image and Diffraction Plane Pictures,” 35, 6 (1972).

9. P. Kocsis, I. Shevkunov, V. Katkovnik, H. Rekola, and K. Egiazarian, “Single-shot pixel super-resolution phase imaging by wavefront separation approach,” Opt. Express 29(26), 43662 (2021). [CrossRef]  

10. J. R. Fienup, “Phase retrieval with continuous version of hybrid input-output,” in Frontiers in Optics (OSA, 2003), p. ThI3.

11. F. Pfeiffer, “X-ray ptychography,” Nat. Photonics 12(1), 9–17 (2018). [CrossRef]  

12. M. Li, L. Bian, G. Zheng, A. Maiden, Y. Liu, Y. Li, J. Suo, Q. Dai, and J. Zhang, “Single-pixel ptychography,” Opt. Lett. 46(7), 1624 (2021). [CrossRef]  

13. X. Huang, J. Nelson, J. Steinbrener, J. Kirz, J. J. Turner, and C. Jacobsen, “Incorrect support and missing center tolerances of phasing algorithms,” Opt. Express 18(25), 26441 (2010). [CrossRef]  

14. P. Thibault and M. Guizar-Sicairos, “Maximum-likelihood refinement for coherent diffractive imaging,” New J. Phys. 14(6), 063004 (2012). [CrossRef]  

15. A. Barty, J. Küpper, and H. N. Chapman, “Molecular Imaging Using X-Ray Free-Electron Lasers,” Annu. Rev. Phys. Chem. 64(1), 415–435 (2013). [CrossRef]  

16. H. N. Chapman, A. Barty, S. Marchesini, A. Noy, S. P. Hau-Riege, C. Cui, M. R. Howells, R. Rosen, H. He, J. C. H. Spence, U. Weierstall, T. Beetz, C. Jacobsen, and D. Shapiro, “High-resolution ab initio three-dimensional x-ray diffraction microscopy,” J. Opt. Soc. Am. A 23(5), 1179 (2006). [CrossRef]  

17. K. Giewekemeyer, P. Thibault, S. Kalbfleisch, A. Beerlink, C. M. Kewish, M. Dierolf, F. Pfeiffer, and T. Salditt, “Quantitative biological imaging by ptychographic x-ray diffraction microscopy,” Proc. Natl. Acad. Sci. U.S.A. 107(2), 529–534 (2010). [CrossRef]  

18. M. F. Duarte, M. A. Davenport, D. Takhar, J. N. Laska, T. Sun, K. F. Kelly, and R. G. Baraniuk, “Single-pixel imaging via compressive sampling,” IEEE Signal Process. Mag. 25(2), 83–91 (2008). [CrossRef]  

19. M. P. Edgar, G. M. Gibson, and M. J. Padgett, “Principles and prospects for single-pixel imaging,” Nature Photon 13(1), 13–20 (2019). [CrossRef]  

20. A. M. Maiden, M. J. Humphry, F. Zhang, and J. M. Rodenburg, “Superresolution imaging via ptychography,” J. Opt. Soc. Am. A 28(4), 604 (2011). [CrossRef]  

21. M. Stockmar, P. Cloetens, I. Zanette, B. Enders, M. Dierolf, F. Pfeiffer, and P. Thibault, “Near-field ptychography: phase retrieval for inline holography using a structured illumination,” Sci. Rep. 3(1), 1927 (2013). [CrossRef]  

22. J. Zhang, J. Sun, Q. Chen, J. Li, and C. Zuo, “Adaptive pixel-super-resolved lensfree in-line digital holography for wide-field on-chip microscopy,” Sci. Rep. 7(1), 11777 (2017). [CrossRef]  

23. W. Luo, Y. Zhang, A. Feizi, Z. Göröcs, and A. Ozcan, “Pixel super-resolution using wavelength scanning,” Light Sci Appl 5(4), e16060 (2016). [CrossRef]  

24. W. Xu, H. Lin, H. Wang, and F. Zhang, “Super-resolution near-field ptychography,” Opt. Express 28(4), 5164 (2020). [CrossRef]  

25. L. Bian, J. Suo, Q. Dai, and F. Chen, “Experimental comparison of single-pixel imaging algorithms,” J. Opt. Soc. Am. A 35(1), 78 (2018). [CrossRef]  

26. M. Aßmann and M. Bayer, “Compressive adaptive computational ghost imaging,” Sci. Rep. 3(1), 1545 (2013). [CrossRef]  

27. Q. Yang, L. Cao, H. Zhang, H. Zhang, and G. Jin, “Method of lateral image reconstruction in structured illumination microscopy with super resolution,” J. Innov. Opt. Health Sci. 09(03), 1630002 (2016). [CrossRef]  

28. F. Zhang, B. Chen, G. R. Morrison, J. Vila-Comamala, M. Guizar-Sicairos, and I. K. Robinson, “Phase retrieval by coherent modulation imaging,” Nat Commun 7(1), 13367 (2016). [CrossRef]  

29. J. Zhao, W. Xu, J. Yi, B. Wang, and F. Zhang, “Extended coherent modulation imaging for single-shot object retrieval free from illumination artifacts,” Ultramicroscopy 240, 113591 (2022). [CrossRef]  

30. J. R. Fienup, “Invariant error metrics for image reconstruction,” Appl. Opt. 36(32), 8352 (1997). [CrossRef]  

31. Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image Quality Assessment: From Error Visibility to Structural Similarity,” IEEE Trans. on Image Process. 13(4), 600–612 (2004). [CrossRef]  

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (14)

Fig. 1.
Fig. 1. Geometry of far-field imaging system. (a) Modulated diffraction pattern occupies a larger area comparing with unmodulated one due to the convolution effect and strong scattering nature of modulator. ${\otimes} $ is the convolution operator, $F\{{\cdot} \}$ is the Fourier transform operator, $AS\{{\cdot} \}$ is the angular-spectrum propagation operator, P is the probe, O is the object, M is the modulator, D is the detector, ${Z_1}$ is the distance between object and modulator, ${Z_2}$ is the distance between modulator and detector. (b) The beamstop blocks diffraction data distributed in central region. (c) The diffraction data in the gap between modules cannot be detected. (d) Data distributed beyond detector is lost. (e) The size of physical pixels is larger than the virtual pixels that obeys the Nyquist requirement.
Fig. 2.
Fig. 2. The illustration of the performance enhanced CMI. (a) The combined diffraction pattern is the input of our algorithm. (b) Flowchart of proposed algorithm.
Fig. 3.
Fig. 3. (a-b) Ground truth of amplitude and phase. (c) Simulated diffraction pattern and the detection area are marked by white dotted squares. (d) Simulated diffraction pattern without object.
Fig. 4.
Fig. 4. Crop from center $1m{m^2}$ of reconstructed amplitude and phase when ${\eta _A} = 4,9,16$. (a-c) Reconstructed object by conventional method. (d-f) Results of our method.
Fig. 5.
Fig. 5. Analysis for reconstruction results and convergency when ${\eta _A} = 4,9,16$. (a) The linecut plots of recovered amplitude by different methods. (b) Evolution of $\textrm{RMS}$ over 200 iterations with different methods.
Fig. 6.
Fig. 6. Analysis for recovered diffraction pattern when ${\eta _A} = 16$. (a-b) Measured and recovered diffraction pattern, respectively. (c)The matrix of $\epsilon $. The red dotted line and white square indicate the extent of input diffraction pattern. Green circle marks the region extrapolated, successfully. (d) Ring average of (c). The red and green line correspond to the red and green circle of (c).
Fig. 7.
Fig. 7. Apparatus of experiment.
Fig. 8.
Fig. 8. Experiment results. (a) Crop center $1.1m{m^2}$ from reconstructed amplitude using full diffraction pattern. And it was divided by the known probe. (b) The acutal diffraction pattern. White dotted squares indicate the extent of detector which corresponds to ${\eta _A} = 2.25,\; 4,\; 6.24,\; 9$, respectively. (c) Diffraction pattern in absence of sample.
Fig. 9.
Fig. 9. Visual comparation of recovered resolution target. (a-d) Recovered amplitude of center $1.1m{m^2}$ by conventional CMI after the division of known probe. (e-h) Results of the same region by our method.
Fig. 10.
Fig. 10. Comparation of different methods by $\textrm{SSIM}$.
Fig. 11.
Fig. 11. Analysis for recovered diffraction pattern when $ {\eta _A} = 9$. (a-b) Real and recovered diffraction pattern. Dotted square is the detector region. (c) The $\epsilon $ matrix. And the while square indicates the extent of detector. Green circle is the part extrapolated successfully. (d) Ring average of (c). And dotted red and green lines represent the size of detector and region extrapolated, respectively.
Fig. 12.
Fig. 12. Quantitative evaluation of the parameters of modulator. (a) Plot of the minimum extrapolation energy versus radius of low-pass filter for continues-type modulator. (b) Plot of the minimum extrapolation energy versus grain size for binary-type modulator.
Fig. 13.
Fig. 13. Experimental results of our method when central diffraction data is missing. (a) Real diffraction pattern. And the white squares represent the extent of central missing data. (b) The cropped amplitude recovered with full diffraction pattern after division of the know probe. (c-f) The cropped amplitude recovered when central diffraction data is missing.
Fig. 14.
Fig. 14. Center $1.1m{m^2}$ region of recovered amplitude under different distance. (a-e) Retrieved amplitude when distance between sample and modulator ${Z_1} = 11.66mm$. (f-j) Recovered amplitude when ${Z_1} = 30.68mm$.

Equations (13)

Equations on this page are rendered with MathJax. Learn more.

R M S = 1 | Φ ( r d ) Ψ ( r d ) | 2 | Ψ ( r d ) | 2 | Φ ( r d ) | 2
ϵ = ( | ψ ( r d ) | I r ( r d ) ) 2 I r ( r d ) + C .
S S I M = ( 2 μ ψ μ ϕ + C 1 ) ( 2 σ ψ ϕ + C 2 ) ( μ ψ 2 + μ ϕ 2 + C 1 ) ( σ ψ 2 + σ ϕ 2 + C 2 ) .
η A = A c ( r d ) A ( r d ) .
η E = I r ( r d ) A ( r d ) I r ( r d ) A c ( r d ) .
ψ k ( r o ) = ψ ^ k 1 ( r o ) S
ψ k ( r m ) = A S { ψ k ( r o ) }
ψ k ( r M ) = ψ k ( r m ) M ( r m )
ψ k ( r d ) = F { 2 ψ k ( r M ) ψ k 1 ( r M ) }
ψ ^ k ( r d ) = A ( r d ) ψ k ( r d ) I ( r d ) | ψ k ( r d ) | + γ ( A c ( r d ) A ( r d ) ) ψ k ( r d ) I ( r d ) | ψ k ( r d ) |
ψ ^ k ( r M ) = F 1 { ψ ^ k ( r d ) }
ψ ^ k ( r m ) = ψ k ( r m ) + M ( r m ) | M ( r m ) | max 2 ( ψ ^ k ( r M ) ψ k ( r M ) + ψ k 1 ( r M ) )
ψ ^ k ( r o ) = A S 1 { ψ ^ k ( r m ) }
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.