Improvement of optical image by measurement reduction technique at parametric multiplexing

TL;DR

This paper demonstrates how measurement reduction techniques can enhance the quality of optical images obtained through parametric multiplexing, leveraging quantum correlations and prior object information.

Contribution

It introduces a novel approach combining measurement reduction with parametric multiplexing to improve optical image quality using quantum correlations and object sparsity.

Findings

Enhanced image quality through measurement reduction.

Utilization of quantum correlations in multiplexed images.

Improved image reconstruction by incorporating object sparsity.

Abstract

In the process of parametric optical image amplification, images are formed at new frequencies in addition to the amplified original image. We show that the parametric multiplexing of optical images can be used to produce an image with improved quality. As an example, we study the parametric amplification of an optical image at low-frequency pumping in which multiplexed optical images turn out to be quantum-correlated. Additional improvement is made possible by using the information about the object that is available to the researcher, in particular, about sparsity of its image. To take the available information into account, we apply the measurement reduction technique.

Figures

Equations

$$A_j^{\text{in}}({\mathbf{r}}_1)=\frac{1}{{\lambda }_{j}f}\underset{-\infty }{\overset{\infty }{\int }}{A}_{j0}(\mathbf{\rho })\exp \left(-i\frac{2\pi }{{\lambda }_{j}f}\mathbf{\rho }{\mathbf{r}}_1\right)d\mathbf{\rho },$$

$$A_j(\mathbf{r})=\frac{1}{{\lambda }_{j}f}\underset{-\infty }{\overset{\infty }{\int }}\left[A_j^{\text{out}}({\mathbf{r}}_2)P({\mathbf{r}}_2)+(1-P^2({\mathbf{r}}_2){)}^{1/2}\hat{v}({\mathbf{r}}_2)\right]\exp \left(-i\frac{2\pi }{{\lambda }_{j}f}\mathbf{r}{\mathbf{r}}_2\right)d{\mathbf{r}}_2.$$

$$A_j(\mathbf{r})=\frac{(2\pi {)}^2}{{\lambda }_{j}f}\underset{-\infty }{\overset{\infty }{\int }}{a}_j^{\text{out}}(\mathbf{q})P\left(\mathbf{q}-\frac{k_j}{f}\mathbf{r}\right)d\mathbf{q},$$

$$a_j(\mathbf{q})=\frac{1}{(2\pi {)}^2}\underset{-\infty }{\overset{\infty }{\int }}{A}_j(\mathbf{r})e^{-i\mathbf{qr}}d\mathbf{r},$$

$$[A(\mathbf{r},z),A^{†}({\mathbf{r}}^{\prime },z)]=\delta (\mathbf{r}-{\mathbf{r}}^{\prime }),[A(\mathbf{r},z),A({\mathbf{r}}^{\prime },z)]=0,$$

$$\begin{array}[]{l}\omega_{p}=\omega_{1}+\omega_{2},\\ \omega_{p}+\omega_{1}=\omega_{3}.\end{array}$$

$$\left\{\begin{array}[]{l}\displaystyle\frac{\partial A_{1}}{\partial z}-\frac{i}{2k_{1}}\triangle_{\perp}A_{1}=i\beta A_{2}^{\dagger}+i\gamma A_{3},\\ \displaystyle\frac{\partial A_{2}}{\partial z}-\frac{i}{2k_{2}}\triangle_{\perp}A_{2}=i\beta A_{1}^{\dagger},\\ \displaystyle\frac{\partial A_{3}}{\partial z}-\frac{i}{2k_{3}}\triangle_{\perp}A_{3}=i\gamma A_{1}.\end{array}\right.$$

$$A_j(\mathbf{r},z)=\underset{-\infty }{\overset{\infty }{\int }}{a}_j(\mathbf{q},z)e^{i\mathbf{qr}}d\mathbf{q},$$

$$\left\{\begin{array}[]{l}\displaystyle\frac{da_{1}}{dz}=-i\mu_{1}a_{1}+i\beta a_{2}^{\dagger}+i\gamma a_{3},\\ \displaystyle\frac{da_{2}^{{\dagger}}}{dz}=i\mu_{2}a_{2}{{}^{\dagger}}-i\beta a_{1},\\ \displaystyle\frac{da_{3}}{dz}=-i\mu_{3}a_{3}+i\gamma a_{1},\end{array}\right.$$

$$\mathbf{a}=Q{\mathbf{a}}_0,$$

$$Q=\left(\begin{array}{ccc}{Q}_{11}& Q_{12}& Q_{13}\\ Q_{21}& Q_{22}& Q_{23}\\ Q_{31}& Q_{32}& Q_{33}\end{array}\right).$$

$$\begin{array}[]{l}\left\langle\hat{N}_{3}(\mathbf{r})\right\rangle=S_{p}|Q_{33}(k_{3}r/f)|^{2}\left\langle\hat{N}_{30}(-\mathbf{r})\right\rangle,\\ \left\langle\hat{N}_{2}(\mathbf{r})\right\rangle=(\lambda_{3}/\lambda_{2})^{2}S_{p}|Q_{23}(k_{2}r/f)|^{2}\left\langle\hat{N}_{30}(-(\lambda_{3}/\lambda_{2})\mathbf{r})\right\rangle,\\ \left\langle\hat{N}_{1}(\mathbf{r})\right\rangle=(\lambda_{3}/\lambda_{1})^{2}S_{p}|Q_{13}(k_{1}r/f)|^{2}\left\langle\hat{N}_{30}(-(\lambda_{3}/\lambda_{1})\mathbf{r})\right\rangle\end{array}$$

$$\frac{1}{l_{j1}}+\frac{1}{l_{j2}}=\frac{1}{f_j}.$$

$$⟨{\hat{N}}_j^{tr}(\mathbf{r})⟩={\left(\frac{l_{j1}}{l_{j2}}\right)}^2⟨{\hat{N}}_j\left(-\frac{l_{j1}}{l_{j2}}\mathbf{r}\right)⟩={\left(\frac{{\lambda }_3{l}_{j1}}{{\lambda }_j{l}_{j2}}\right)}^2{S}_p{\left|Q_{j3}\left(\frac{k_{2}r{l}_{j1}}{f{l}_{j2}}\right)\right|}^2⟨{\hat{N}}_{30}\left(\frac{{\lambda }_3{l}_{j1}}{{\lambda }_j{l}_{j2}}\mathbf{r}\right)⟩.$$

$$\frac{l_{j1}}{l_{j2}}=\frac{{\lambda }_j}{{\lambda }_3}.$$

$$⟨{\hat{N}}_j^{(tr)}(\mathbf{r})⟩=S_p{\left|Q_{j3}\left(\frac{k_{3}r}{f}\right)\right|}^2⟨{\hat{N}}_{30}(\mathbf{r})⟩.$$

$$\begin{array}[]{l}\sigma^{2}_{3}=S_{p}\left[1+2|Q_{32}(k_{3}r)/f)|^{2}\right]|Q_{33}(k_{3}r/f)|^{2}\left\langle\hat{N}_{30}(-\mathbf{r})\right\rangle,\\ \sigma^{2}_{2}=S_{p}\left[1+2\frac{S_{a}S_{p}}{(f\lambda_{3})^{2}}|Q_{21}(k_{3}r/f)|^{2}\right]|Q_{23}(k_{3}r/f)|^{2}\left\langle\hat{N}_{30}(\mathbf{r})\right\rangle,\\ \sigma^{2}_{1}=S_{p}\left[1+2|Q_{12}(k_{3}r/f)|^{2}\right]|Q_{13}(k_{3}r/f)|^{2}\left\langle\hat{N}_{30}(\mathbf{r})\right\rangle.\end{array}$$

$$C_{jl}(\mathbf{r})=⟨{\hat{N}}_j(\mathbf{r}){\hat{N}}_l(\mathbf{r})⟩-⟨{\hat{N}}_j(\mathbf{r})⟩⟨{\hat{N}}_l(\mathbf{r})⟩$$

$$C_{31}(\mathbf{r})=C_{13}(\mathbf{r})=\frac{S_a{S}_p^2}{(f{\lambda }_3{)}^2}\left\{|Q_{12}(k_{3}r/f)Q_{33}(k_{3}r/f){|}^2+|Q_{13}(k_{3}r/f)Q_{32}(k_{3}r/f){|}^2\right\}⟨{\hat{N}}_{30}(\mathbf{r})⟩,$$

$$C_{32}(\mathbf{r})=C_{23}(\mathbf{r})=\frac{S_a{S}_p^2}{(f{\lambda }_3{)}^2}|Q_{23}(k_{3}r/f){|}^2\left\{|Q_{32}(k_{3}r/f){|}^2+|Q_{33}(k_{3}r/f){|}^2\right\}⟨{\hat{N}}_{30}(\mathbf{r})⟩,$$

$$C_{12}(\mathbf{r})=C_{21}(\mathbf{r})=\frac{S_a{S}_p^2}{(f{\lambda }_3{)}^2}|Q_{23}(k_{3}r)/f){|}^2\left\{|Q_{12}(k_{3}r/f){|}^2+|Q_{13}(k_{3}r)/f){|}^2\right\}⟨{\hat{N}}_{30}(\mathbf{r})⟩,$$

$$\begin{array}[]{l}Q_{12}(0)=-Q_{21}(0)=i(\beta/\Gamma)\sinh\Gamma z,\\ Q_{13}(0)=Q_{31}(0)=i(\gamma/\Gamma)\sinh\Gamma z,\\ Q_{23}(0)=-Q_{32}(0)=(\beta\gamma/\Gamma^{2})(\cosh\Gamma z-1),\\ Q_{33}(0)=1-(\gamma/\Gamma)^{2}(\cosh\Gamma z-1),\end{array}$$

$$\xi =Ag+\nu ,$$

$$A=\left(\begin{array}{c}{B}_1{C}_1\\ B_2{C}_2\\ B_3{C}_3\end{array}\right).$$

$${\Sigma }_{\nu }=\left(\begin{array}{ccc}{B}_1{\Sigma }_{11}(g)B_1^{\ast }& B_1{\Sigma }_{12}(g)B_2^{\ast }& B_1{\Sigma }_{13}(g)B_3^{\ast }\\ B_2{\Sigma }_{21}(g)B_1^{\ast }& B_2{\Sigma }_{22}(g)B_2^{\ast }& B_2{\Sigma }_{23}(g)B_3^{\ast }\\ B_3{\Sigma }_{31}(g)B_1^{\ast }& B_3{\Sigma }_{32}(g)B_2^{\ast }& B_3{\Sigma }_{33}(g)B_3^{\ast }\end{array}\right)+{\Sigma }_{{\nu }^{\prime }}(g).$$

$$R_{\ast }\stackrel{\text{def}}{=}U(A^{\ast }{\Sigma }_{\nu }^{-1}A{)}^{-}{A}^{\ast }{\Sigma }_{\nu }^{-1},$$

$${\Sigma }_{R_{\ast }\xi }=U(A^{\ast }{\Sigma }_{\nu }^{-1}A{)}^{-1}{U}^{\ast }.$$

$$\hat{u}={\Pi }_{{\Sigma }_{R_{\ast }\xi }}\left({\tilde{R}}_{{\Sigma }_{R_{\ast }\xi }}{\left({\xi }^T,{\hat{u}}^T\right)}^T\right)$$

$${\Pi }_{{\Sigma }_{R_{\ast }\xi }}(u)\stackrel{\text{def}}{=}\underset{v\in {\mathcal{U}}_{\text{pr}}}{\mathrm{argmin}}(v-u,{{\Sigma }_{R_{\ast }\xi }}^{-1}(v-u))$$

Full text

Introduction

As it is well-known, in traditional parametric amplification of an optical image with high-frequency pumping, an additional image appears at the so-called idle frequency (see, for example [1, 2]). In the case of processes of optical parametric amplification with low-frequency pumping that can be realized in coupled parametric processes [3, 4, 5] optical images are formed at more than two frequencies [6, 7, 8, 9, 10]. In other words, in coupled parametric interactions frequency multiplexing of optical images occurs. The quantum theory of such interactions is developed in [6, 7, 8, 10], where various parametric image amplification schemes are examined. Experimental studies on the optical image multiplexing are presented in the articles [9, 10].

In the works [6, 7, 8] optical image quality is characterized by the signal-to-noise ratio (SNR). It is established that SNR of the image at the main frequency decreases in spite of increasing the mean photon number. Meanwhile, SNR of the images at additional frequencies increases along with their mean photon number.

Recently, image multiplexing has been used in ghost image acquisition schemes [11, 12, 13, 14, 15, 16, 17]

It is shown in our works [11, 12, 13, 14] that using quantum correlations of ghost images and the reduction technique in image processing, we can improve the characteristics of the reconstructed optical image.

The purpose of this paper is to study the quantum correlations of optical images during their multiplexing in two coupled parametric processes and apply the method of measurement reduction for obtaining an optical image with improved quality. The analysis is based on the coupled parametric interactions that realize parametric image amplification with low-frequency pumping, and the scheme with far-away object is considered. Note that some issues of such a process are studied in [7].

This process is of particular interest when the wavelength of the original image is in the ultraviolet range and the pump radiation wavelength of the conventional three-frequency interaction falls within the absorption region of a nonlinear crystal. The process of parametric amplification with low-frequency pumping allows the use of visible radiation as pumping, while the images at additional frequencies will be in the near-infrared range.

The article structure is as follows. In section 1, we discuss the optical setup for image amplification and multiplexing and the parametric processes taking place therein. In section 2 a specific variant of the optical setup is considered and the photon number means, variances and covariances are derived. In section 3 the measurement model and the measurement reduction method is outlined, including the notions of a measuring transducer, an ideal measuring transducer and the conditions for the possibility of image reconstruction. The information about the object that is available to the researcher and that is employed in reduction is summarized in subsection 3.1. In subsection 3.2, the algorithm of image processing using reduction method that takes this information into account is described. Computer modeling results are given in section 4. Main results of the article are summarized in the conclusion.

1 Amplification and multiplexing of an optical image

The scheme of parametric optical image amplification is depicted in Fig. 1. The weak optical image that is to be amplified is located in the object plane $P_{1}$. This image is projected by the lens $L_{1}$ onto the input $P_{2}$ of an aperidic nonlinear photon crystal (ANPC), for example, LiNbO3, in which coupled parametric interactions, specifically, down- and up-conversion processes, occur simultaneously. The amplified image and the images generated at two new frequencies (see below) are projected by the lens $L_{2}$ from the output of the crystal (plane $P_{3}$) onto the image plane $P_{4}$. The lenses $L_{1}$ and $L_{2}$ have the same focal length $f$. The object and image planes, as well as the input and the output of the ANPC, are at the distance $f$. This is the scheme with the so-called far-away object, it is similar to the scheme considered in [2, 7].

We denote the field operators in the object and the image planes as $A_{j0}(\boldsymbol{\rho})$ and $A_{j}(\mathbf{r})$, respectively, and those in the input and the output planes of the NPC as $A_{j}^{\textnormal{in}}(\mathbf{r}_{1})$ and $A_{j}^{\textnormal{out}}(\mathbf{r}_{2})$, respectively. The index $j$ is associated with the wavelength $\lambda_{j}$ . The operators $A_{j0}$, $A_{j}^{\textnormal{in}}$ and $A_{j}^{\textnormal{out}}$, $A_{j}$ are related by the Fourier transformation performed by the lens $L_{1}$

[TABLE]

and by the lens $L_{2}$

[TABLE]

$P(\mathbf{r}_{2})$ is the pupil frame function that accounts for the finite area $S_{a}$ of the pupil. Taking it into account is necessary for the correct analysis of the image amplification scheme [2] due to vacuum fluctuations outside of the pupil’s aperture, described by the operator $\hat{v}(\mathbf{r}_{2})$ (the second term in Eq. (2)). However, they do not contribute to any normal-ordered operational expressions associated with measurable values, so we omit this term below. Therefore, the expression (2) can be presented as

[TABLE]

where $a_{j}^{\textnormal{out}}(\mathbf{q})$ and $P(\mathbf{q})$ are the Fourier transforms of $A_{j}^{\textnormal{out}}(\mathbf{r})$ and $P(\mathbf{r})$, respectively:

[TABLE]

$\mathbf{q}$ is the transversal wave vector.

All operators under consideration $A_{j0}(\mathbf{r})$, $A_{j}(\mathbf{r})$, $A_{j}^{\textnormal{in}}(\mathbf{r})$ and $A_{j}^{\textnormal{out}}(\mathbf{r})$ obey the commutation relations

[TABLE]

where $z$ is the direction of wave propagation. Mean value $\langle\hat{N}(\mathbf{r},z)\rangle$ of the operator $\hat{N}(\mathbf{r},z)=A^{\dagger}(\mathbf{r},z)A(\mathbf{r},z)$ is the mean photon flux density in cross-section $z$, measured in photons per cm2.

To find the connection between operators $A_{j}^{\textnormal{in}}(\mathbf{r}_{1})$ and $A_{j}^{\textnormal{out}}(\mathbf{r}_{2})$, i. e. the connection between fields at the output and the input of the ANPC, it is necessary to consider the nonlinear processes in the ANPC. The processes under study are two coupled processes

[TABLE]

Here $\omega_{p}$ is the frequency of intense pump wave, and $\omega_{1}$, $\omega_{2}$ and $\omega_{3}$ are the frequencies of generated waves with $\omega_{p}$ and $\omega_{1}$ being the shared frequencies of two processes. The first down-conversion process in Eq. (6) represents parametric amplification during high-frequency pumping, and the second one is the up-conversion process. They can be implemented simultaneously in an aperiodical NPC.

In the undepleted pump plane wave approximation taking into account the diffraction phenomenon the processes in Eq. (6) can be described by the system of equations

[TABLE]

Here $\triangle_{\perp}=\triangle_{\perp}(x,y)$ is the transversal Laplacian, $A_{j}^{{\dagger}}=A_{j}^{{\dagger}}(\mathbf{r},z)$ and $A_{j}=A_{j}(\mathbf{r},z)$ are the creation and annihilation operators of photons with frequency $\omega_{j}$ ($j=1,2,3$) respectively, $\beta$ and $\gamma$ are real nonlinear coupling coefficients that are proportional to second order nonlinear susceptibility and the absolute value of pump wave amplitude [3]. Eqs. (7) are derived for a lossless ANPC and for interaction of monochromatic waves.

The system of Eqs. (7) is solved by applying Fourier transform

[TABLE]

after which the Eqs. (7) become

[TABLE]

where $a_{j}=a_{j}(\mathbf{q},z)$, $\mu_{j}=\frac{q^{2}}{2k_{j}}$.

In the matrix form the solution of Eqs. (9) has the form

[TABLE]

where $\mathbf{a}_{0}^{T}=(a_{10},a_{20}^{\dagger},a_{30})$ is determined by the values of the operators at ANPC input ($z=0$), the index $T$ denotes transposition. In the case considered below, the operators $\mathbf{a}_{10}(\mathbf{q})$, $\mathbf{a}_{20}(\mathbf{q})$ describe the vacuum state and the operator $\mathbf{a}_{30}(\mathbf{q})$ describes a coherent state.

The matrix $Q$ consists of transfer functions $Q_{nm}$:

[TABLE]

The functions $Q_{nn}=Q_{nn}(q,z)$ are the self-transfer functions because they describe the amplification at frequencies $\omega_{n}$ ($n=1,2,3$), while the cross-transfer functions $Q_{nm}=Q_{nm}(q,z)$ describe the conversion from frequency $\omega_{m}$ to $\omega_{n}$. Elements of the matrix (11) can be found in [7] and are not given here due to being cumbersome.

2 Formulation of the problem

In the previous section, the quantum theory of two coupled parametric processes is presented in relation to amplification and frequency conversion of an optical image, which can arrive at the nonlinear crystal at any of the frequencies $\omega_{1}$, $\omega_{2}$ or $\omega_{3}$. Here we turn to the case when the image with a mean photon number density $\left\langle\hat{N}_{30}(\mathbf{r})\right\rangle$ is fed to the crystal in a coherent state at the frequency $\omega_{3}$. In the framework of the monochromatic waves under consideration, results given below are valid if the image registration time is less than the characteristic time of image change, for example, the correlation time.

In the image plane $P_{4}$, the mean photon number over pixel area $S_{p}$ in the amplified image and the additional ones is given by the expressions

[TABLE]

As expected, the output optical images are inverted relative to the initial image. It is important to note that the scales of output images at different frequencies are different, and the change in spatial scale is determined by the coefficient $\lambda_{3}/\lambda_{j}$, where $\lambda_{j}$ is the image wavelength. It should be noted that this fact was not taken into account in [7].

Different image scales somewhat complicate the reduction algorithm, as this means that sizes of pixels are different in different arms. One can proceed further in several ways.

• •

In the general approach without the assumption that the image is piecewise constant (see sec. 3), dealing with different pixel sizes in different arms is avoided, since in the infinite-dimensional case pixel size does not affect image representation.

• •

The images can be rescaled during processing if sensor point spread functions allow to do this both accurately and without loss of data.

• •

Finally, one can bring the images to the same scale using additional lenses. This approach is considered below.

It should be noted that these approaches provided the same results when they are valid.

In the image plane $P_{4}$, images with different frequencies are directionally separated (for example, using a prism). A lens with the focal distance $f_{j}$ ($j=1,2$) is placed on the path of radiation with frequency $\omega_{j}$ in order to bring the image to the spatial scale that coincides with the scale of the image at the frequency $\omega_{3}$. The lens is located at a distance $l_{j1}$ from the plane $P_{4}$ and a distance $l_{j2}$ to the measurement plane. The specified distance must satisfy the lens formula

[TABLE]

In this case, the distribution of the image mean photon number in the optically conjugate plane is given by ($j=1,2$):

[TABLE]

According to Eq. (14), to match the image scales at the wavelengths $\lambda_{j}$ and $\lambda_{3}$ the ratio

[TABLE]

must be satisfied. Under this condition, we have

[TABLE]

The photon number variances $\sigma^{2}_{j}=\left\langle\hat{N^{2}_{j}}(\mathbf{r})\right\rangle-\left\langle\hat{N}_{j}(\mathbf{r})\right\rangle^{2}$ are given by the following formulas

[TABLE]

Finally, the mutual correlations of the image fluctuations (covariances) between different frequencies

[TABLE]

have the forms

[TABLE]

[TABLE]

[TABLE]

To estimate the above moments, we use their values at $\mathbf{r}=0$. In this case, the matrix elements (11) have a simple analytical form:

[TABLE]

where $\Gamma=(\beta^{2}-\gamma^{2})^{1/2}$. $|Q_{33}(0)|^{2}$ as a function of $\gamma/\beta$ and $\beta z$ is shown in Fig. 2.

It follows from Eq. (22) that the amplification of image takes place if $\beta>\gamma$.

In the absence of a down-conversion process ($\beta=0$) and, therefore, without generation of additional frequencies, there is no amplification. In the coupled parametric process under the condition $\beta>\gamma$, the original image initially decays and at the interaction length $z_{0}$, $\cosh z_{0}=(\beta/\gamma)^{2}$, the mean photon number $\langle\hat{N}_{30}(\mathbf{r})\rangle$ of the image becomes zero. The image gain process begins after the interaction length $z_{m}$, $\cosh z_{m}=2(\beta/\gamma)^{2}-1$ (see also Fig. 2). As for the photon numbers $\langle\hat{N}_{1}(\mathbf{r})\rangle$, $\langle\hat{N}_{2}(\mathbf{r})\rangle$ at frequencies $\omega_{1}$ and $\omega_{2}$, they monotonously grow as the interaction length increases according to expressions (12), (22).

3 Processing of acquired images

The output of sensors in the $i$-th arm, denoted as $\xi^{(i)}(\mathbf{r})$, can be considered as the output of a measuring transducer (MT) for input signal $g(\mathbf{r})\sim\langle\hat{N}_{30}(-\mathbf{r})\rangle$.

We will consider piecewise constant images, i. e. transparency of the research object is constant within each pixel. Areas of constant transparency and constant brightness corresponding to pixels are considered to be ordered in an arbitrary but fixed way. Due to that it is sufficient for us to consider a finite number of values of $\mathbf{r}$. Thus, $g$ as the vector of transparencies is an element of finite-dimensional Euclidean space $\mathcal{F}$.

This assumption is made for simplicity (in order to avoid working with infinite-dimensional spaces) and is not crucial to the reduction method. For examples of reduction in the infinite-dimensional case see, e. g., [18, ch. 10] and [19]. If the assumption is invalid, the estimate of the algorithm below estimates linear combinations of $\langle\hat{N}_{30}(-\mathbf{r})\rangle$ around the values of $\mathbf{r}$ determined by locations of the sensors with weights that are dependent on the MT and the ideal MT $U$ specified by the researcher (see below). For example, if sensors have uniform light sensivity throughout their area, the ideal MT corresponds to unit-sized sensors of unit size and the images were acquired by sensors that are several times larger, $g$ is averaged over unit-sized areas centered at sensor locations.

Let us formulate the measurement model as

[TABLE]

where $g$ is an unknown vector that describes the transparency distribution of the object, $\nu$ is the measurement error with zero expectation, $\operatorname{\mathbb{E}}\nu=0$, which means absence of systematic measurement error, and covariance matrix $\Sigma_{\nu}=\operatorname{\mathbb{E}}\nu\nu^{*}$. The dimension of vector $g$ is the number of pixels in the object image, while the dimension of $\xi$ is the total number of sensors. The condition of systematic measurement error absense $\operatorname{\mathbb{E}}\nu=0$ means, in particular, that the expectation of the component of measurement results caused by sensor dark noises is subtracted from the measurement results, similar to [20, 14] for ghost images.

The matrix $A$ describes image amplification, multiplexing and acquisition: the matrix element $A_{ij}$ is equal to the mean output of $i$-th sensor for unit transparency of $j$-th element of the illuminated object and zero transparency of other object elements (i. e. whose indices differ from $j$). Due to the measuring setup with three arms, it is a block matrix and consists of three blocks describing different arms:

[TABLE]

Under the conditions used to derive the mean numbers of photons, their variances and mutual correlations the matrices $C_{1}$–$C_{3}$ are identity ones multiplied by the factor before $\langle\hat{N}_{30}(-\mathbf{r})\rangle$ in expression (16) for the mean numbers of photons. The matrices $B_{1}$–$B_{3}$ model the sensors. Specifically, the matrix element $(B_{i})_{pk}$ is equal to the output of the sensor in $i$-th arm at $p$-th position for unit brightness of $k$-th pixel of the image formed in that arm and zero brightness of other pixels. For the same reason, the noise covariance matrix has block form as well:

[TABLE]

Here the element with indices $k$, $k^{\prime}$ of the block $\Sigma_{ij}$ is equal to either photocount variance in Eq. (17) if $i=j$ or covariance of photocounts in Eq. (19)–(21) if $i\neq j$ for the same pixel ordering as in the matrix $\mathbf{A}$. Hence, the dependence of the matrix (25) on $g$ is caused by the dependence of variances and covariances on $\langle\hat{N}_{30}(-\mathbf{r})\rangle$. The term $\Sigma_{\nu^{\prime}}$ is the covariance matrix of the noise component $\nu^{\prime}$ that is unrelated to parametric amplification and multiplexing, e. g. thermal noise in circuits, detection of outside photons, non-unit quantum efficiency of the sensors and their dark noise.

The objective of an image processing algorithm is to output the most accurate estimate of the signal $Ug$ from the measurement result $\xi$, where the matrix $U$ describes a measuring device that is ideal (for the researcher). Hence, $Ug$ is the feature of the original image $g$ that is of interest to the researcher. We consider the case when the researcher is interested in reconstruction of the object image itself, and imaging does not distort the object, therefore, $U=I/n$, where $n$ is the average number of photons per pixel of the illuminated object. One way of achieving this is the measurement reduction method described in [18], see also [21, 22, 23, 24]. If the estimation process is described by a linear operator $R$ ($R\xi$ is the result of processing the measurement $\xi$), the corresponding mean squared error (MSE) in the worst case of $g$, $h(R,U)=\sup\limits_{g\in\mathcal{F}}\operatorname{\mathbb{E}}\lVert R\xi-Ug\rVert^{2}$, as shown in [18], is minimal for $R$ that is equal to the linear unbiased reduction operator

[TABLE]

where - denotes pseudoinverse. $h(R_{*},U)=\operatorname{tr}U(A^{*}\Sigma_{\nu}^{-1}A)^{-1}U^{*}$, and the covariance matrix of the linear reduction estimate $R_{*}\xi$ is

[TABLE]

Estimation is possible (MSE is finite) if the condition $U(I-A^{-}A)=0$ holds, where, as noted above, $A$ characterizes the real measuring device, while $U$ characterizes an ideal one with the point spread function required by the researcher, and, therefore, the desired resolution, if this condition if fulfilled. This condition essentially means that only the features of the object that are measured by the real measuring device (that is, affect its output) can be estimated. Unlike fluorescence-based superresolution techniques, see e. g. [25], the proposed technique does not require attaching fluorescent molecules to the object. However, in addition to the above condition, the error of the obtained estimate can be too large to distinguish the signal from the noise in practice, and usually the better the desired resolution of the ideal measuring device compared to the resolution of the real one, the larger MSE of the obtained estimate. Nevertheless, by choosing $U$ one can select an acceptable (to him) compromise between obtained resolution and noise magnitude, that is, to estimate with acceptable resolution and with an tolerable noise level.

In the case under consideration, as seen from Eq. (24), the diagonal elements of $C_{1}$–$C_{3}$ are nonzero. Therefore, each block of $C_{j}$ is non-degenerate, so for non-degenerate $B_{j}$ the reduction error takes only finite values.

The measurement reduction technique for the case when it is known that $u=Ug\in\mathcal{U}_{\textnormal{pr}}\subset\mathcal{U}$, where $\mathcal{U}_{\textnormal{pr}}$ is convex and closed, in other words, when the feature of interest of the object is known to satisfy certain given constraints, was considered in [24, 26]. The linear estimate (26) is refined using this information by solving the equation

[TABLE]

for $\hat{u}$, where $\tilde{R}_{\Sigma_{R_{*}\xi}}$ is the measurement reduction operator for a MT $\left(A^{T},U^{T}\right)^{T}$ and noise with covariance matrix $\begin{pmatrix}\Sigma_{\nu}&0\\ 0&\Sigma_{R_{*}\xi}\end{pmatrix}$, and the operator

[TABLE]

describes orthogonal projection onto $\mathcal{U}_{\textnormal{pr}}$ by minimizing the Mahalanobis distance $\lVert\Sigma_{R_{*}\xi}^{-1/2}\cdot\rVert$ associated with the covariance matrix $\Sigma_{R_{*}\xi}$ (27) of the linear reduction estimate $R_{*}\xi$. The earlier version of reduction technique proposed in [13] and in [24] for similar information used minimization of the ‘‘ordinary’’ Euclidean distance instead of Mahalanobis distance. In [26], the advantages of minimizing Mahalanobis distance instead of Euclidean distance during projection are shown. For such prior information, the covariance matrix (27) of linear reduction estimate error is an upper bound on the covariance matrix of the obtained estimate $\hat{u}$.

3.1 Prior object information

It is obvious that a priori the transparency distribution of the object takes values in $[0,1]$, hence $g\in[0,n]^{\dim\mathcal{F}}$, $Ug\in[0,1]^{\dim\mathcal{F}}$, where the average number $n$ of photons per pixel of the illuminated object is assumed to be known.

It is assumed that the transparency distribution of the object is not ‘‘entirely’’ arbitrary: transparencies of neighboring pixels usually do not differ much. As a result, the image is sparse (many of its components are zero) in a given basis, similarly to compressed sensing [27, 28, 29, 30]. The hypothesis ‘‘$i$-th component in the given basis of the estimate $\hat{u}$ is zero’’ is treated as a statistical hypothesis that is tested using the measurement data against the alternative that it is nonzero. Its testing is controlled by choosing the significance level (the probability of rejecting the hypothesis when it is true) of the rejection criterion or a parameter $\tau$ of the criterion that monotonously depends on it.

The researcher also knows the matrix $A$ (24) that describes image acquisition conditions and, up to the vector $g$, the matrix $\Sigma_{\nu}$ (25) that describes the magnitudes of measurement errors. Note that the worst case of $g$ is realized if all pixels are equally transparent.

3.2 Reduction algorithm

The proposed algorithm of multiplexed GI processing using measurement reduction technique that is based on the indicated prior information has the following form.

Calculation of linear reduction estimate $R_{*}\xi$ (26) based on the acquired images. For calculation of the covariance matrix (25), the worst case, that all pixels have the same brightness, is assumed. 2. 2.

Refinement of the estimate $R_{*}\xi$ using the information $\mathcal{U}_{\textnormal{pr}}=[0,1]^{\dim\mathcal{F}}$ by the method (28) by fixed-point iteration, i. e. by consecutive application of the mapping (28) with $\Pi_{\mathbf{\Sigma}_{R_{*}\xi}}(R_{*}\xi)$ as the initial approximation. We denote the obtained estimate by $\hat{u}$. 3. 3.

Application of the sparsity-inducing transformation $T$ to $\hat{u}$. ‘‘Sparsity-inducing’’ means that the researcher expects the chosen transform of the true transparency distribution of the object to be sparse. 4. 4.

Calculation of the worst-case (in $g$) variances $\sigma_{T\hat{u}}^{2}=(\sigma_{(T\hat{u})_{1}}^{2},\dots,\sigma_{(T\hat{u})_{\dim\mathcal{F}}}^{2})$ of the components of $T\hat{u}$ (the diagonal matrix elements of $T\mathbf{\Sigma}_{R_{*}\xi}T^{*}$) and calculation of $T\hat{u}_{\textnormal{thr}}$ in the following way: $(T\hat{u}_{\textnormal{thr}})_{i}\operatorname{\overset{\text{def}}{=}}0$ if $|(T\hat{u})_{i}|<\tau\sigma_{(T\hat{u})_{i}}$, otherwise $(T\hat{u}_{\textnormal{thr}})_{i}\operatorname{\overset{\text{def}}{=}}(T\hat{u})_{i}$. 5. 5.

Inverse transformation $T^{-1}$ of $T\hat{u}_{\textnormal{thr}}$ (if $T$ is a unitary transformation, then $T^{-1}=T^{*}$), i. e. calculation of $\hat{u}_{\textnormal{thr}}\operatorname{\overset{\text{def}}{=}}T^{-1}T\hat{u}_{\textnormal{thr}}$. 6. 6.

Calculation of the projection $\Pi_{\mathbf{\Sigma}_{R_{*}\xi}}(\hat{u}_{\textnormal{thr}})$ that is considered to be the result of processing.

The algorithm parameter $\tau\geq 0$ reflects a compromise between noise suppression (the larger the value of $\tau$, the greater the noise suppression) and distortion of images whose components are close to [math]. As mentioned above, the step 4 can be considered as testing statistical hypotheses $(TUf)_{i}=0$ (for the alternative $(TUf)_{i}\neq 0$) for all $i$. In this paper the criterion used in step 4 is based on Chebyshev’s inequality: if $(TUf)_{i}=0$, then $\Pr\left(|(T\hat{u})_{i}|\geq\tau\sigma_{(T\hat{u})_{i}}\right)\leq\tau^{-2}$ (hence, the significance level is at least $\tau^{-2}$). Step 4 can be also interpreted as replacement of the original matrix $U$ with one whose kernel contains the estimate components after the specified transform that are affected by noise of the specified magnitude or more.

4 Computer modeling results

The results of image processing according to the described algorithm are shown in Figs. 3–7. The computer modeling was carried out for wave lengths $(\lambda_{1})^{-1}=$1.2\text{,}{\mathrm{\SIUnitSymbolMicro m}}^{-1}, $(\lambda_{2})^{-1}=$0.8\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-1}, $(\lambda_{3})^{-1}=$3.2\text{,}{\mathrm{\SIUnitSymbolMicro m}}^{-1}, aperture area $S_{a}=$25\text{\,}{\mathrm{cm}}^{2}, pixel area $S_{p}=$100\text{,}{\mathrm{\SIUnitSymbolMicro m}}^{2}, focal distance $f=$10\text{\,}\mathrm{cm} and the value of crystal parameter $\epsilon=\gamma/\beta$ and the dimensionless crystal length $\beta z$ indicated in figure captions, with $\epsilon$ ranging from $0.4$ to $0.8$ and $\beta z$ ranging from $1$ to $5$. The sensors in arms are identical ones that are three times as large as an element of the object image. Therefore, image processing via measurement reduction increases resolution in addition to noise suppression. It should be noted, however, that the objectives of superresolution and reconstruction of the image with a small number of photons are generally at odds with each other: relaxing resolution requirements allows to reconstruct the image using less photons, as less components of the image have to be recovered.

One can see that additional information about sparsity enables higher noise suppression without compromising reducing obtained resolution too much. Increase of $\tau$ leads to better noise suppression (cf., e. g., Figs. 3g and 3h), but also worse distortions caused by discarding ‘‘significant’’ image components as well (cf., e. g., Figs. 4e and 4f). Too large values of $\tau$ cause degradation of image fidelity due to distortion outweighing improved noise suppression, as small-scale image details are suppressed as well. Therefore, one should choose the maximal value of $\tau$ that preserves the details of interest. To do that, one can model processing of a test image that contains the required details and choose the largest value of $\tau$ that preserves them, or visually compare the reduction results for different $\tau$ and select its final value by ternary search.

The transform whose result for the transparency distribution of the object is sparse that is usually employed in image processing by the means of compressed sensing is discrete cosine transform (DCT) [28, 29, 30]. In [31], several transforms (identity transform, discrete wavelet transform and DCT) were reviewed and the advantages of DCT were shown. However, in [14] it was shown that Haar transform may be preferable in the case of a transparency distribution that contains areas of weakly changing transparency with sharp borders if these areas are large compared to the resolution of the ideal MT and the location of the borders is important to the researcher, since those features align well with the vectors of the Haar transform basis. As the object images in Figs. 3–7 are of this type, the sparsity-inducing transform used in these figures is the Haar transform.

In all these figures results of processing of all images are compared to results of processing of a single image, namely, the one with the best signal-to-noise ratio. Since processing only a single image means that some of the obtained measurements are discarded, this results in worse estimate quality, even if the sum of acquired images is not noticeably different from the best image. For example, in Fig. 5 one can see that the same values of $\tau$ cause more distortion when processing a single image. In Fig. 5i only a single component, corresponging to uniform transparency, remains in the result of processing a single image, while the same value of $\tau$ produces an acceptable, although suboptimal, image shown in Fig. 5g. Furthermore, in the general case of additional noise multiplexing provides the means for further noise suppression if noise photons in different arms are detected independently.

As far as the parameters of the ANPC are concerned, larger crystal length $\beta z$ leads to more amplification and more photons. The improvement of the number of photons is especially large when $\beta z$ is increased from $1$ to $2$. The results do not depend on the value of the parameter $\epsilon$ as much.

Conclusion

The problem of recovering images acquired in photon-sparse conditions can be solved by taking advantage of the additional information available to the researcher about the measurement process and about the object. Alternatively one can make the detection conditions worse (e. g. use sensors with less resolution) while preserving the same estimation quality. In this work, the additional information about the object is the information that the object transparency distribution is not arbitrary, namely, transparencies of close pixels tend to be close. This information is formalized as sparsity of the result of a given transform of the transparency distribution, similar to compressed sensing.

In compressed sensing, as a rule, the measurement error is modeled as an arbitrary vector with bounded norm. In the proposed method it is modeled as a random vector, and selection of the estimate components which are considered to be zero is based on the statistical properties of the estimate components, namely, their variances. The use of covariances of the estimate components in addition to their variances is a subject of further research. Another subject of further research are the opportunities provided by multiplexing for analysis of the measurement data, e. g., for verifying the reliability of both the measurement model [32, 33] and the results of reduction [34].

We consider that computer modeling based on the developed algorithm showed high efficiency of the developed reduction technique for parametric amplification of images and frequency conversion in the sense of improvement of both their quality and their noise immunity.

Finally, we emphasize once again that we are dealing with parametric interaction with low-frequency pumping. The carrier wavelength of the original image may be in the ultraviolet range, while the pump frequency may belong to the visible range. Further applications of the measurement reduction technique to processing of quantum images are being developed. This work was supported by Russian Foundation for Basic Research, grant 18-01-000598 A.

Acknowledgments

The authors acknowledge discussions with Prof. A. V. Belinsky.