Reference-less measurement of the transmission matrix of a highly scattering material using a DMD and phase retrieval techniques

Angélique Drémeau; Antoine Liutkus; David Martina; Ori Katz; Christophe Schülke; Florent Krzakala; Sylvain Gigan; Laurent Daudet

doi:10.1364/OE.23.011898

1. Introduction

Wave propagation in complex media is a fundamental problem in physics, be it in acoustics, optics, or electromagnetism [1]. In optics, it is particularly relevant for imaging applications. Indeed, when light passes through a multiply scattering medium, such as a biological tissue or a layer of paint, ballistic light is rapidly attenuated, preventing conventional imaging techniques, and random scattering events generate a so-called speckle pattern that is usually considered useless for imaging. Recently, wavefront shaping using spatial light modulators (SLM) has emerged as a unique tool to manipulate multiply scattered coherent light, for focusing or imaging in scattering media [2]. In essence, these methods use the linearity and time-reversal symmetry of the wave propagation, whatever the complexity of the medium, to control the output speckle field, by manipulating the light beam impinging on the scattering sample. Different wavefront shaping approaches rely on digital phase-conjugation [3, 4] or iterative algorithms [5], but it is also possible to measure the so-called transmission matrix (TM) of the medium [6], which fully describes light propagation through the linear medium, from the modulator device to the detector. This approach has been particularly efficient for focusing, imaging [7, 8] and for studying the transmission modes of the medium [9]. These methods are not only valid for scattering material but can also be applied to other complex transmission system, most notably multimode fibers, turning them into minimal footprint endoscopes [10–13].

A major limitation of most of these techniques for imaging is their speed. Indeed, the wavefront shaping process must be faster than the stability time of the medium, which can be of only a few milliseconds in biological tissues. Yet, most of the works reported so far have relied on phase modulators which are usually slow (few tens of Hertz for liquid crystal modulators). Micro Electro-Mechanical Systems (MEMS) modulators are much faster, but are usually not phase-only. As a promising alternative for wave shaping in complex media, Digital Micromirror Device (DMD) technology [14] offers binary amplitude modulators (i.e., “ON” or “OFF”) operating at > 20kHz, with high pixel counts (10⁶) and low pitch (around 10 microns), all this at low cost. These binary amplitude modulators have been used as phase modulators, using appropriate diffraction and filtering, e.g. by Lee-type amplitude holography [15, 16], as shown on Fig. 1(b). While phase control is more effective for wavefront shaping than amplitude control, some works reported on using DMD as genuine binary amplitude modulators for wavefront shaping through opaque scattering media, albeit usually yielding lower overall efficiency than phase modulators for focusing or mode matching [17–19]. The DMD configuration can also be optimized using genetic algorithms [20] to maximize the intensity enhancement.

Fig. 1 Different experimental approaches for measuring the complex-valued transmission matrix of a scattering medium with a binary DMD amplitude modulator: (a) using a reference arm for retrieving the phase of the output field by off-axis or phase-shifting holography; (b) using the DMD as a spatial phase modulator by displaying amplitude holograms, and using the unmodulated parts of the field as a phase-stable reference; (c) presented approach where only intensity values are measured.

Download Full Size | PDF

For the measurement of a TM, an additional issue lies in accessing the amplitude and phase of the output field, that in optics usually requires a holographic measurement, i.e., a reference beam, as shown on Fig. 1(a). This reference beam can either be co-propagating in the medium [7, 21], or use an external reference arm [8, 22]. The phase and amplitude of the measured field can then be extracted by simple linear combinations of interference patterns with a phase-shifted or off-axis reference. This however poses the unavoidable experimental problem of the interferometric stability of the reference arm.

In this work we report on the full measurement of the complex TM of a multiply scattering medium, using a DMD binary amplitude modulator as an SLM, with no reference on the detection side, as shown on Fig. 1(c). This approach combines the high-speed and high pixel counts allowed by DMD devices, with the simplicity and robustness of a reference-less optical setup. However, it involves advanced signal processing algorithms for phase retrieval, run on a sufficiently large number of input-output calibration measurements. Here, we use a Bayesian phase retrieval algorithm [23] for the estimation of a TM, based on actual noisy experimental measurements. This is done one output pixel at a time, a process that can easily be parallelized for speed. Although a number of phase estimation techniques could also be used, our approach provides good results at moderate computational costs. This technique is validated on experimental data, with three criteria. First, for a given arbitrary input pattern, we compare the measured output and the predicted output computed with the measured TM. Then, we show that the distribution of the singular values of the measured TM varies according to random matrix theory. Finally, we demonstrate that single- or multi-point light focusing can be achieved. Interestingly, our Bayesian formulation allows the use of the same model for the estimation of the optimal DMD binary input pattern. In addition to being an interesting signal processing problem, this approach is particularly relevant for real-life applications of the TM approach, since it allows a simple, fast and robust implementation.

2. Experimental setup

Our experimental setup, described in Fig. 2, uses a DMD-array from Texas Instrument (1920×1080 tilting micromirrors), driven by the DLP V-9500 VIS module (Vialux). The DMD is made of mirrors that can switch between two angular positions separated by 24°, thus reflecting each pixel either toward a beam dump (pixel “OFF”) or towards the focusing system (pixel “ON”).

Fig. 2 Experimental scheme: A 532 nm CW laser is expanded through a telescope in order to obtain an homogeneous beam. Through a rectangular mask, it illuminates the DMD which acts as binary amplitude spatial light modulator. The DMD reflects the light in two different directions corresponding to either “ON” (unit transmission) or “OFF” (the light is deviated towards a beam dump). The transmitted pattern is focused by a first lens L1 on the scattering medium – here a white paint layer –, acting as a thick multiply scattering medium. The transmitted speckle pattern is collected by a microscope objective and is observed through a polarizer P on a CCD camera.

Download Full Size | PDF

Under Matlab, an amplitude mask is computed and loaded on the DMD. The pattern corresponding to the ON pixels is focused on the surface of a thick scattering medium by means of a f = 100 mm lens L1 (thus the DMD pixels correspond roughly to incidence angles on the sample). The sample is a ∼ 100 microns thick layer of white paint, which is thick enough in order to considerably mix the light on the other side, producing a complex speckle interfering pattern. This speckle pattern is collected through a microscope objective (L2) and detected on a camera (AVT Pike F-100B). In order to measure the TM, we need to send a large series of input patterns (typically a few times the number of input pixels we wish to control), in a time over which the medium can be considered stationary. For this purpose, we use the “high speed” driver provided with the DMD in order to load all the to-be-projected random amplitude masks to the memory of the DMD driver module, and we trigger the display of each mask via a DAQ card (National Instruments, PCI-6221) and a waveform generator. In the same way, in order to be as fast as possible, we also only consider a subregion on the camera of size of 400×400 pixels. The overall acquisition rate is 31 images per second. To monitor the stability of the medium, we periodically measure the correlation of the speckle image corresponding to the same input mask. We therefore quantify the stability of the medium, which is better than 98% over the total measurement time (typically around 5 minutes).

3. Phase retrieval for optics through complex media

The experimental setup described in the previous section presents different challenges in order to exploit the TM framework: the absence of reference requires a solution exploiting only intensity measurements, and this must be robust to experimental noise. This paper shows that signal processing techniques for phase retrieval can be used, first for the estimation of the TM, and then for getting the optimal DMD configuration in a focalization task.

3.1. Calibration as a phase retrieval problem

If the focalization task - or in a more general view, any front-shaping task using spatial light modulators - can be easily expressed as a phase retrieval problem - knowing the TM, which SLM setup best explains the intensity observations or the desired output - , this is not the case for the estimation of the TM.

We formalize the latter as a calibration problem: given P incoming waves, assumed perfectly known, which model explains at best the observed outputs?

Formally, let x_μ ∈ {0,1}^N stand for the binary DMD inputs related to the μ-th acquisition, where N is the number of pixels (mirrors) used on the DMD. We assume that the partial observations of the sole moduli of the transmitted waves (the square root of the camera measured intensities), denoted by $y_{μ} \in ℝ_{+}^{M}$ , obey

\begin{matrix} y_{μ} = | D x_{μ} |, & \forall μ \in {1, ..., P}, \end{matrix}

where D is the unknown complex-valued transmission matrix characterizing the scattering material, and M is the number of observed pixels on the camera.

Then, adopting a matrix formulation and conjugating-transposing the system, we get

Y^{H} = | X^{H} D^{H} |,

where Y = [y₁,…,y_P], X = [x₁,…,x_P] and .^H denotes the conjugate-transpose of a matrix/vector. This reveals a “classic” phase retrieval problem: given the matrix of inputs X^H, each column of Y^H is used to estimate each complex-valued column of D^H.

3.2. Phase retrieval techniques

The problem of reconstructing a complex vector given only the magnitude of measurements is a non-convex optimization problem notoriously difficult to solve. Many algorithms have been devised in the literature to deal with this problem, especially in optics where most approaches aims at reconstructing a phase object from an intensity image in a Fourier plane [24, 25]. Here we are interested in a somewhat simpler and more general problem where the transform is completely general and the output points can be treated independently. Amongst recent contributions to this signal processing problem, we can mention [26,27] which approximate the phase recovery problem by relaxed problems enabling the use of standard optimization procedures, and the Bayesian procedures proposed in [23, 28], which circumvent the non-linearity of the modulus through the introduction of hidden variables and resort to variational approximations.

In this paper, we chose to use the so-called prVBEM Bayesian approach introduced in [23]. Our choice is motivated by two of its main advantages:

its generic framework, enabling its use in both applications we are here interesting in, namely the calibration of the TM and the focalization task,
its reasonable computational complexity, of crucial importance for a complete procedure including the calibration of the TM and its immediate use for focusing.

We will not detail here the derivation of the algorithm. Its particularizations to the two problems of interest, calibration of the TM and focalization, are described resp. in appendix 7.1 and 7.2. For a methodological justification of the approach, we refer the reader to the original paper [23].

4. Estimating the TM with intensity-only measurements and binary inputs

In this section, we focus on the calibration stage of the proposed approach. We assess the accuracy of the TM estimated through the phase retrieval technique [23] with regard to two different analyses: its prediction performance, namely its ability to predict the observations, and its nature, expected to satisfy random features.

4.1. Prediction performance

To assess the prediction performance of the estimated TM, we adopt a cross-validation-like experimental framework. The setup is as follows. We measure the M = 40000 camera pixels stemming from N = 900 DMD mirrors, 50% of them being turned on, the others off at each displayed pattern. The operation is repeated randomly P = 6000 times. Given this dataset, a row of the TM is then learned from p = αN calibration measurements, with α varying in {1,…,6}, and used in a second step to predict the P−p remaining measurements. This estimation is performed on 50 different rows of the TM.

The prVBEM algorithm presents of complexity of order $O (p^{2})$ . Run for 200 iterations, this leads to a computational time of order 0.6 s per estimated row of the TM, in Matlab, on a Mac-book Air with a 1.7Ghz i7 processor. This performance is good in comparison with other state-of-the-art algorithms (see e.g., [23]) and, keeping in mind that rows are independent, sufficient for making possible a complete procedure “calibration of the TM-focalization” in a reasonable time.

The prediction performance of the estimated TM is evaluated according to the mean-square error (MSE) and the normalized cross-correlation between the moduli of the P−p predicted measurements and the actual observed ones, such as, for the estimate $\hat{d}$ of the current (conjugate) row of the TM (see Eq. (2)),

M S E = \frac{‖ y_{P - p} - | X_{P - p}^{H} \hat{d} ‖ |_{2}^{2}}{‖ y_{P - p} | |_{2}^{2}},

C o r r = \frac{y_{P - p}^{H} | X_{P - p}^{H} \hat{d} |}{‖ y_{P - p}^{H} y_{P - p} | |_{2} ‖ {\hat{d}}^{H} X_{P - p} X_{P - p}^{H} \hat{d} | |_{2}},

where y_P₋_p stands for the current (real-valued) row of Y (see Eq. (2)) restricted to its last P−p elements, the first p ones being used for the estimation of

\hat{d}

.

These figures of merit give different insights into the estimations achieved by the algorithm: the MSE measures the distance (in an euclidian sense) to the actual observations, while the correlation evaluates their angular difference. Figure 3 show these quantities (resp. in (a) and (b)) averaged over the 50 rows of the TM considered for estimation (red lines), as well as the minimum and maximum (given by whiskers), for an increasing α = p/N.

Fig. 3 Prediction performance according to (a), the mean-square error (MSE), in log scale, and to (b), the normalized cross-correlation between observation predictions using the estimated TM, and actual measurements of the output moduli (square root of the camera intensity values), as a function of the number of calibration measurements (x-axis is α, such that p = αN calibration measurements are used).

Download Full Size | PDF

Both curves match admirably: increasing α leads to a decrease of the MSE and an increase of the correlation. Interestingly, we see that for α ≥ 3, that is, for at least 3 times more real measurements than complex unknowns, the estimation of the TM is accurate enough to predict the observations with an average correlation around 0.95 and an average MSE lower than 0.18 (Fig. 3(a) shows this value in log-scale).

4.2. Comparison of singular values to Random matrix theory

Interestingly, we can check that the measured TM presents some characteristics as predicted by random matrix theory. One practical way is to verify that the distribution of its normalized singular values obeys the Marčenko-Pastur law [29]. It should be noted that such apparently random signals are the hardest case for phase retrieval, where no specific structure can be taken into account.

In order to reduce the influence of specifics of our experimental setting, we perform the following operations, as in [30]:

We normalize over the rows and columns, to attenuate the illumination artifacts: residual illumination “by default” on each pixel of the camera for the rows, and inhomogeneous contribution of each DMD mirror on the entire set of camera pixels for the columns.
Because of the size of the speckle grains, two neighboring DMD mirrors may affect the material in the same way, as well, two pixels of the camera will be potentially correlated. To avoid this effect, we subsample the rows and columns of the matrix.

To draw the empirical spectral density, we then consider the following setup. We subsample the columns of the matrix up to N = 200 and leave the number of rows varying, more precisely M = γN, with γ ∈ {1,…,6}. These sub-matrices thus constitute partitions of the estimated matrix, randomly picked 100 times to average the resulting densities. Figure 4 compares the experimental curves to the theoretical ones drawn according to the Marčenko-Pastur law. We see that the experiments qualitatively follow the predictions. We remark however that the larger γ is, the more chances we have to consider the contributions of neighboring correlated pixels. This partly explains the increasing gap between both curves.

Fig. 4 Density of the normalized singular values for different γ = M/N. Stamped line: experimental results, continuous line: Marčenko-Pastur law.

Download Full Size | PDF

5. Focusing with the DMD

Knowing the TM gives a powerful and flexible tool to control light within the scattering medium [30]. In particular, it can be used to compute which DMD input has to be set, in order to display a given arbitrary pattern at the receiver end. In this section, we demonstrate the special case of focusing light with maximum intensity on a desired pattern (a chosen sparse subset of the output pixels), with the TM measured experimentally as in the section 4. It should be emphasized that we keep the same experimental setup, with the binary DMD as input device. Here, simple inversion methods such as [30] cannot be used, as these require phase-modulated input and amplitude and phase detection.

We propose here to resort to the Bayesian variational approach [23], in a similar way as for the calibration but adapted here to the binary nature of the DMD inputs. The particularization of the algorithm to this case is given in appendix 7.2.

We assess the performance of the proposed focusing approach through different experiments. The general setting is as follows. The DMD inputs, here taken of dimension N = 1600, are estimated from the desired outputs (see appendix 7.2), focusing on 1 to 4 target points. The procedure exploits the TM, reduced to its rows of interest and previously measured as discussed in section 4.

Figure 5 shows an example of the observed output field, corresponding to the estimated DMD configuration, optimized to focus on 3 points. To quantitatively evaluate the focusing performance, we measure the intensity enhancement factor, as:

η ≜ \frac{I_{foc}}{I_{back}},

where I_foc is the intensity inside the target area after spatial binary amplitude modulation is performed, I_back is the average background intensity. This value is measured for 100 trials, as a function of the number of calibration measurements used to learn the TM.

Fig. 5 Illustration of light focusing on 3 points. The circles mark the positions of the targets.

Download Full Size | PDF

Two different setups are then considered: the single-point focusing case and the multi-target case.

5.1. Focusing on a single point

Figure 6 compares the enhancement factors achieved by two different focusing methods, namely a simple binary phase-conjugation - performing $\hat{x} = [ℜ (D^{H} y) > 0]$ - and the proposed method, in the case where only one target point is focused. Results are presented under a “box” format, where:

the middle segment stands for the average enhancement $\bar{η}$ over the 100 trials,
the upper and lower bounds of the rectangle define the interval $[\begin{matrix} \bar{η} - σ_{η} & \bar{η} + σ_{η} \end{matrix}]$ (where σ_η is the experimentally computed standard deviation), in which lies, under the Gaussian assumption, 68 % of the trials,
the whiskers represent the minimum and maximum values observed over the entire set of trials.

Fig. 6 Single target experiment. Enhancement factor as a function of the number of measurements used to learn the TM (x-axis is α, such that p = αN calibration measurements are used). For the same estimation of the TM, 2 focusing techniques are compared: binary phase conjugation (blue boxes), and the Bayesian technique [23] (red boxes).

Download Full Size | PDF

For each experiment point α ∈ {2,…, 6}, such that p = αN calibration measurements are used to compute the TM, we display the boxes related to the binary phase-conjugation method (blue boxes), and the Bayesian technique (red boxes) [23] particularized in appendix 7.2.

As a first observation, we can see that the general dependency with regard to α noticeably resonates with the curve of the prVBEM algorithm in Fig. 3(b): there is a clear gap between the performance achieved for α = 2 and for α = 3, while, for α ≥ 3, the intensity enhancement keeps increasing but less significantly.

Interestingly, the prVBEM algorithm seems to outperform binary phase-conjugation, with regard to the mean and maximum values measured, but not in a statistically significant manner. Focusing on the most favorable case considered here, namely with α = 6, the best intensity enhancement factor lies around 140, to be compared with the ideal expected enhancement given by $1 + \frac{1}{π} (\frac{N}{2} - 1) ≃ 255$ , see [17].

5.2. Focusing on multiple points

For this second setup, we are interested in the performance of the prVBEM approach in a context of multiple target points. Additionally to the intensity enhancement, we consider here the missed detection rate, defined as the number of trials (expressed in percentage) failing to focus on at least one of the multiple target points, i.e., the number of trials for which at least one of the T largest intensity peaks in the output image does not match any of the T targets.

Figure 7 represents these two figures of merit under diagram formats. They present an interesting general symmetry: increasing the number of targets or decreasing the number of calibration points leads to an increase of the missed detections and a decrease of the enhancement factor. The missed detection rate seems however more sensitive to the number of calibration points used to learn the TM: for α = 2 and 2 target points, the algorithm fails with a rate approaching 40%, while for α = 3 and the same number of targets, we keep a reasonable performance (around 10%). In a more general view, these figures greatly highlight the deep relation between the quality of the calibration and the focusing performance.

Fig. 7 Multiple target experiment. (a) Average enhancement factor as a function of the number of measurements used to learn the TM (x-axis is α, such that p = αN calibration measurements are used), and the number of target points (y-axis). (b) Missed detection rate (same axis as in (a)).

Download Full Size | PDF

6. Conclusion

This paper shows that the full complex-valued transmission matrix of a strongly scattering material can be estimated, up to a global phase factor on each of its rows, with a simple experimental setup involving only real-valued inputs and outputs. In our experiment, the inputs are amplitude modulations on a binary DMD, and the output is the field intensity measured on a CCD camera, that gathers a significant amount of measurement noise. Note that no reference arm is used, that would allow interferometric measurements, but that would make the experimental setup more complex and considerably more unstable.

We here resort to Bayesian phase retrieval techniques, and we have shown that, amongst such techniques, a recently proposed variational approach (the so-called prVBEM algorithm [23]) allows a precise estimation of the transmission matrix, tractable in computational complexity and scalable for large-size signals, provided that we have a sufficiently large number of input-output calibration signals. Experimental results validate this concept, both in terms of output prediction, distribution of singular values, and in an application of light focusing onto a number of target points in the output plane. It should be emphasized that this estimation of the transmission matrix opens many applications beyond light focusing, may it be for imaging through the scattering material [30,31], or for obtaining information about the scattering material itself.

7. Appendix: particularization of the prVBEM algorithm

In this appendix, we outline the main lines of the prVBEM algorithm particularized to its use for calibration of the TM and focalization. We refer the reader to the original paper [23] for a methodological presentation of the procedure.

7.1. prVBEM for the calibration of the TM

To circumvent the non-linearity induced by the absence of phase observations, the procedure introduces new variables modeling, on the one hand, the missing phases of the observations, and on the other hand, some acquisition noise. Thus, recalling that we resort to a conjugate-transposition of the matrix system (see Eq. (2)), each absolute-valued measurement y_μ, μ ∈ {1…P}, of any row y of Y, is expressed as

y_{μ} = e^{j θ_{μ}} (\sum_{i = 1}^{N} x_{μ i} d_{i}^{*} + ω_{μ}),

where θ_μ ∈ [0,2π) stands for its missing conjugate phase, x_μi is the ith element of the μth row in X,

d_{i}^{*}

corresponds to the ith conjugate element in the current estimated row d of D and ω_μ is an additive noise, assumed centered isotropic Gaussian (denoted

C N

in the following) with variance σ². We moreover suppose that the probability distributions for the entries of the matrix and for the missing phases are:

\begin{matrix} p (d) = \prod_{i = 1}^{N} p (d_{i}) & with & p (d_{i}) = C N (0, σ_{d}^{2}), \end{matrix}

\begin{matrix} and & p (θ) = \prod_{μ = 1}^{P} p (θ_{μ}) & with & p (θ_{μ}) = \frac{1}{2 π} . \end{matrix}

Under these assumptions, the absence of phases in the observations is naturally taken into account in the model since marginalizing on θ_μ leads to a distribution on y_μ which only depends on the moduli of y_μ and $\sum_{i = 1}^{N} x_{μ i} d_{i}^{*}$ .

Within model (6)–(8), the recovery of the complex row d of the TM can be expressed as the solution of the following marginalized Maximum A Posteriori (MAP) estimation problem

\hat{d} = \arg \max p (d | y),

\begin{matrix} with & p (d | y) = \int_{θ} p (d, θ | y) \end{matrix} .

Because of the marginalization on the hidden variables θ, the direct computation of p(d|y) is however intractable in general. The prVBEM algorithm is based on the computation of a particular (the so-called “mean-field”) approximation $\hat{q} (d, θ) = \prod_{i} q (d_{i}) \prod_{μ} q (θ_{μ})$ of the posterior joint distribution p(d,θ|y). The procedure is then iterative, updating the factors q(d_i) as

q (d_{i}) = C N (m_{i}, Σ_{i}),

with

\begin{matrix} m_{i} = \frac{σ_{d}^{2}}{σ_{n}^{2} + σ_{d}^{2} x_{i}^{H} x_{i}} {〈 r_{i} 〉}^{H} x_{i}, & Σ_{i} = \frac{σ_{n}^{2} σ_{d}^{2}}{σ_{n}^{2} + σ_{d}^{2} x_{i}^{H} x_{i}}, & 〈 r_{i} 〉 = \bar{y} - \sum_{k \neq i} m_{k}^{*} x_{k}, \end{matrix}

and

\bar{y} = {[y_{μ} e^{(j \arg (y_{μ}^{*} 〈 z_{μ} 〉))} \frac{I_{1} (\frac{2}{σ_{n}^{2}} | y_{μ}^{*} 〈 z_{μ} 〉 |)}{I_{0} (\frac{2}{σ_{n}^{2}} | y_{μ}^{*} 〈 z_{μ} 〉 |)}]}_{μ = {1 \dots P}}

〈 z_{μ} 〉 = \sum_{i = 1}^{N} m_{i}^{*} x_{μ i} .

In the equations above, I₀ (resp. I₁) stands for the modified Bessel function of the first kind for order 0 (resp. 1).

Coming back to problem (9), an approximation of p(d|y) simply follows from

p (d | y) = \int_{θ} p (d, θ | y),

≃ \int_{θ} q (d, θ),

≃ \int_{θ} \prod_{i} q (d_{i}) \prod_{μ} q (θ_{μ}),

≃ \prod_{i} q (d_{i}),

resulting in an element-wise estimation of row

\hat{d}

, such as

{\hat{d}}_{i} = m_{i}

.

7.2. prVBEM for focusing

The problem of focusing is expressed as an inverse problem, where, knowing the TM D and the observation y, we look for the DMD input x such as described in (1). Adopting a similar modeling as in previous section, we then assume, for all elements y_μ with μ ∈ {1,…,M},

y_{μ} = e^{j θ_{μ}} (\sum_{i = 1}^{N} d_{μ i} x_{i} + ω_{μ}),

where θ_μ ∈ [0,2π) stands for the missing conjugate phase, d_μi is the μth element of the ith column in D, x_i ∈ {0,1} corresponds to the state of the ith DMD pixel and ω_μ is an additive noise, assumed centered isotropic Gaussian of variance σ². As previously, we suppose that the elements θ_μ are independently and uniformly distributed in the interval [0,2π), however, in order to accommodate for binary inputs, we consider here a Bernoulli model for x:

\begin{matrix} p (x) = \prod_{i = 1}^{N} p (x_{i}) & with & p (x_{i}) = Ber (p_{i}) {\begin{array}{l} p_{i} \\ 1 - p_{i} \end{array} & \begin{matrix} if x_{i} = 1, \\ if x_{i} = 0, \end{matrix} \end{matrix}

where we set p_i to 0.5, noticing that asymptotically half of the DMD pixels are expected to be “ON” [17].

Then, within model (19)–(20), the recovery of the DMD inputs is expressed as the marginalized MAP estimation problem

\hat{x} = \underset{x}{\arg \max} p (x | y),

\begin{matrix} with & p (x | y) = \int_{θ} p (x, θ | y) \end{matrix} .

The procedure is then similar to the one described previously for the calibration task. The prVBEM algorithm iteratively updates the factors of the mean-field approximation q(x,θ) = ∏_i q(xi)∏_μ q(θ_μ) of p(x,θ|y), which leads, in the particular case of focusing, to the expression

q (x_{i}) = p (x_{i}) \exp (x_{i} \frac{2 ℜ (d_{i}^{H} 〈 r_{i} 〉) - d_{i}^{H} d_{i}}{σ^{2}}),

where

〈 r_{i} 〉 = \bar{y} - \sum_{k \neq i} q (x_{k} = 1) d_{k},

\bar{y} = {[y_{μ} e^{(j \arg (y_{μ}^{*} 〈 z_{μ} 〉))} \frac{I_{1} (\frac{2}{σ^{2}} | y_{μ}^{*} 〈 z_{μ} 〉 |)}{I_{0} (\frac{2}{σ^{2}} | y_{μ}^{*} 〈 z_{μ} 〉 |)}]}_{μ = {1 \dots M}}

〈 z_{μ} 〉 = \sum_{i} q (x_{i} = 1) d_{μ i},

and I₀ (resp. I₁) stands for the modified Bessel function of the first kind for order 0 (resp. 1).

Then, an approximation of p(x|y) is given by

p (x | y) = \int_{θ} p (x, θ | y),

≃ \int_{θ} \prod_{i} q (x_{i}) \prod_{μ} q (θ_{μ}),

= \prod_{i} q (x_{i}) .

Using this approximation, the problem is easy to solve by a simple thresholding operation, i.e., ${\hat{x}}_{i} = 1$ if q(x_i = 1) > 0.5 and ${\hat{x}}_{i} = 0$ otherwise.

Acknowledgments

AD is currently working at ENSTA Bretagne, LabSTICC (UMR 6285), 2 rue François Verny, F-29200 Brest, France. OK acknowledges the support of the Marie Curie Intra-European Fellowship for career development (IEF). LD acknowledges a joint research position with the Institut Universitaire de France. This work has been supported in part by the CSI:PSL grant and LABEX WIFI (Laboratory of Excellence within the French Program “Investments for the Future”) under references ANR-10-LABX-24 and ANR-10-IDEX-0001-02 PSL*, and by the ERC under the European Union’s 7th Framework Programme Grant Agreements 307087-SPARCS and 278025-COMEDIA.

References and links

1. P. Sebbah, Waves and Imaging Through Complex Media (Springer, 2001). [CrossRef]

2. A. P. Mosk, A. Lagendijk, G. Lerosey, and M. Fink, “Controlling waves in space and time for imaging and focusing in complex media,” Nature Photonics 6, 283–292 (2012). [CrossRef]

3. M. Cui and C. Yang, “Implementation of a digital optical phase conjugation system and its application to study the robustness of turbidity suppression by phase conjugation,” Opt. Express 18, 3444–3455 (2010). [CrossRef] [PubMed]

4. I. N. Papadopoulos, S. Farahi, C. Moser, and D. Psaltis, “Focusing and scanning light through a multimode optical fiber using digital phase conjugation,” Opt. Express 20, 10583 (2012). [CrossRef] [PubMed]

5. I. M. Vellekoop and A. P. Mosk, “Focusing coherent light through opaque strongly scattering media,” Opt. Express 32, 2309 (2007).

6. S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, “Measuring the transmission matrix in optics: An approach to the study and control of light propagation in disordered media,” Phys. Rev. Lett. 104, 100601 (2010). [CrossRef] [PubMed]

7. S. Popoff, G. Lerosey, M. Fink, A. C. Boccara, and S. Gigan, “Image transmission through an opaque material,” Nature Commun. 1, 81 (2010). [CrossRef]

8. Y. Choi, T. D. Yang, C. Fang-Yen, P. Kang, K. J. Lee, R. R. Dasari, M. S. Feld, and W. Choi, “Overcoming the diffraction limit using multiple light scattering in a highly disordered medium,” Phys. Rev. Lett. 107, 023902 (2011). [CrossRef] [PubMed]

9. M. Kim, Y. Choi, C. Yoon, W. Choi, J. Kim, Q.-H. Park, and W. Choi, “Maximal energy transport through disordered media with the implementation of transmission eigenchannels,” Nature Photonics 6, 581–585 (2012). [CrossRef]

10. Y. Choi, C. Yoon, M. Kim, T. D. Yang, C. Fang-Yen, R. R. Dasari, K. J. Lee, and W. Choi, “Scanner-free and wide-field endoscopic imaging by using a single multimode optical fiber,” Phys. Rev. Lett. 109, 023902 (2012). [CrossRef]

11. I. N. Papadopoulos, S. Farahi, C. Moser, and D. Psaltis, “High-resolution, lensless endoscope based on digital scanning through a multimode optical fiber,” Biomed. Opt. Express 4, 260–270 (2013). [CrossRef] [PubMed]

12. S. Bianchi and R. Di Leonardo, “A multi-mode fiber probe for holographic micromanipulation and microscopy,” Lab on a Chip 12, 635 (2012). [CrossRef]

13. T. Čižmár and K. Dholakia, “Shaping the light transmission through a multimode optical fibre: complex transformation analysis and applications in biophotonics,” Opt. Express 19, 18871 (2011). [CrossRef] [PubMed]

14. J. B. Sampsell, “Dmd display system,” (1995). US Patent 5,452,024.

15. D. B. Conkey, A. M. Caravaca-Aguirre, and R. Piestun, “High-speed scattering medium characterization with application to focusing light through turbid media,” Opt. Express 20, 1733–1740 (2012). [CrossRef] [PubMed]

16. S. A. Goorden, J. Bertolotti, and A. P. Mosk, “Superpixel-based spatial amplitude and phase modulation using a digital micromirror device,” arXiv:1405.3893 [physics] (2014).

17. D. Akbulut, T. J. Huisman, E. G. van Putten, W. L. Vos, and A. P. Mosk, “Focusing light through random photonic media by binary amplitude modulation,” Opt. Express 19, 4017–4029 (2011). [CrossRef] [PubMed]

18. D. Kim, W. Choi, M. Kim, J. Moon, K. Seo, S. Ju, and W. Choi, “Implementing transmission eigenchannels of disordered media by a binary-control digital micromirror device,” Opt. Commun. 330, 35–39 (2014). [CrossRef]

19. J. W. Tay, J. Liang, and L. V. Wang, “Amplitude-masked photoacoustic wavefront shaping and application in flowmetry,” Opt. Lett. 39, 5499–5502 (2014). [CrossRef] [PubMed]

20. X. Zhang and P. Kner, “Binary wavefront optimization using a genetic algorithm,” Journal of Optics 16, 125704 (2014). [CrossRef]

21. T. Chaigne, O. Katz, A. C. Boccara, M. Fink, E. Bossy, and S. Gigan, “Controlling light in scattering media non-invasively using the photoacoustic transmission matrix,” Nature Photonics 8, 58–64 (2014). [CrossRef]

22. D. Akbulut, T. Strudley, J. Bertolotti, T. Zehender, E. Bakkers, A. Lagendijk, W. Vos, O. Muskens, and A. Mosk, “Measurements on the optical transmission matrices of strongly scattering nanowire layers,” in Proc. Conf. on Lasers and Electro-Optics Europe and International Quantum Electronics Conf. (CLEO EUROPE/IQEC), (2013).

23. A. Drémeau and F. Krzakala, “Phase recovery from a bayesian point of view: the variational approach,” in Proceedings of IEEE Trans. Acoust. Speech Signal Process. (2015).

24. R. Gerchberg and W. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik 35, 237–246 (1972).

25. J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21, 2758–2769 (1982). [CrossRef] [PubMed]

26. E. Candès, T. Strohmer, and V. Voroninski, “Phaselift : exact and stable signal recovery from magnitude measurements via convex programming,” Commun. in Pure and Applied Mathematics 66, 1241–1274 (2013). [CrossRef]

27. I. Waldspurger, A. d’Aspremont, and S. Mallat, “Phase recovery, maxcut and complex semidefinite programming,” Mathematical Programming Series A-Springer (2013).

28. P. Schniter and S. Rangan, “Compressive phase retrieval via generalized approximate message passing,” in Proceedings of Communication, Control, and Computing (Allerton) (2012).

29. V. A. Marčenko and L. A. Pastur, “Distribution of eigenvalues for some sets of random matrices,” Mathematics of the USSR-Sbornik 1, 457–483 (1967). [CrossRef]

30. S. M. Popoff, G. Lerosey, M. Fink, A. C. Boccara, and S. Gigan, “Controlling light through optical disordered media: transmission matrix approach,” New Journal of Physics 13, 123021 (2011). [CrossRef]

31. A. Liutkus, D. Martina, S. Popoff, G. Chardon, O. Katz, G. Lerosey, S. Gigan, L. Daudet, and I. Carron, “Imaging with nature: Compressive imaging using a multiply scattering medium,” Sci. Rep.4 (2014). [CrossRef] [PubMed]

Reference-less measurement of the transmission matrix of a highly scattering material using a DMD and phase retrieval techniques

Abstract

1. Introduction

2. Experimental setup

3. Phase retrieval for optics through complex media

3.1. Calibration as a phase retrieval problem

3.2. Phase retrieval techniques

4. Estimating the TM with intensity-only measurements and binary inputs

4.1. Prediction performance

4.2. Comparison of singular values to Random matrix theory

5. Focusing with the DMD

5.1. Focusing on a single point

5.2. Focusing on multiple points

6. Conclusion

7. Appendix: particularization of the prVBEM algorithm

7.1. prVBEM for the calibration of the TM

7.2. prVBEM for focusing

Acknowledgments

References and links

Cited By

Figures (7)

Equations (29)

Optics Express