Compressive Schlieren Deflectometry

Schlieren deflectometry aims at characterizing transparent objects by studying the deflections undergone by refracted incident light rays at any point of their surface point [1]. Compared to other characterization techniques relying on interferometry, deflectometry is also insensitive to vibrations which makes it very attractive for industrial deployment (e.g., for quality control).

Consider a (thin) transparent object shined on one side with a beam of parallel light rays, as shown in Fig. 1(left). At each surface location ${\boldsymbol{p}}$, the refracted light is deviated in multiple directions characterized in a local coordinate system $({\boldsymbol{e}_{1}},{\boldsymbol{e}_{2}},{\boldsymbol{e}_{3}})$, with $\boldsymbol{e}_{3}$ parallel to incident light beam. Using the spherical coordinates $({\theta},{\varphi})$ in this system (see Fig. 1(left)), the resulting deflection spectrum $\tilde{s}_{\boldsymbol{p}}({\theta},{\varphi})\in{\mathbb{R}}_{+}$ represents the flux of light deviated in the direction $({\theta},{\varphi})$. The evolution of deflection spectra along the object surface encodes information about its local curvature. This is why deflectometry is often used for characterizing transparent object surfaces.

In this paper, $\tilde{s}_{\boldsymbol{b}}$ is conveniently represented by its projection in the $\Pi_{\boldsymbol{p}}=\boldsymbol{e}_{2}\boldsymbol{e}_{3}$ plane, i.e., according the projected function $s_{\boldsymbol{p}}(r(\theta),\varphi)=\tilde{s}_{\boldsymbol{p}}(\theta,\varphi)$ with $r(\theta)=\tan\theta$. Moreover, the object surface is assumed sufficiently smooth for being parametrized by a projection of $\boldsymbol{p}$ in the same plane (on a arbitrary fixed origin), so that $\boldsymbol{p}$ is basically a 2-D vector.

For most objects (e.g., with smooth surfaces) deflections at any location ${\boldsymbol{p}}$ occur in a limited range of angles. The deflection spectra therefore tend to be naturally sparse in plane $\Pi_{\boldsymbol{p}}$ or in some appropriate basis of this domain (e.g., wavelets). Fig. 1(right) shows an example of a discretized deflection spectrum $s_{\boldsymbol{p}}$ for one location of a plano convex lens obtained using the setup described in Sec. 2. The white spot in the image signifies that deflections occur at only a few angles (as governed by classical optics) and deflections elsewhere are negligible.

Measuring a deflection spectrum $s_{\boldsymbol{p}}$ for every ${\boldsymbol{p}}$ is not straightforward. The optical setup described in Sec. 2 measures them indirectly by optical comparison with a certain number of programmable modulation patterns. To maximize the amount of information collected in this process, we exploit the sparse nature of spectra in each plane $\Pi_{\boldsymbol{p}}$ and use the framework of spread spectrum¹¹1“Spread Spectrum” is not related to the studied deflection “spectrum” but it refers to the signal frequency spectrum. compressive sensing [2], described in Sec. 4. From a fine calibration of the system relative to its intrinsic noise, our approach is then experimentally validated using deflectometric measurements and the numerical spectrum reconstruction results are presented in Sec. 5. Further, from the possibility of performing signal processing in the compressed domain [3, 4], we also show in Sec. 6 that for localized deflection spectra summarized by a few optical parameters (e.g., location and width of the main peak) the compressed optical observations themselves can be processed. This is aimed at extracting the relevant information at a smaller number of measurements than those required for a full spectrum reconstruction.

Deflection spectra can be measured by the Schlieren deflectometer represented in Fig. 2. Its key components are (i) a Spatial Light Modulator (SLM), (ii) the Schlieren lens with focal length ${f}$, (iii) the Telecentric System (TS) and (iv) the Charged Coupled Device (CCD) camera collecting the light.

The object to be analyzed is placed in between the Schlieren lens and the telecentric system. It is shined on its left by a light source and, thanks to the telecentric system, only the parallel light rays emerging out from the object are collected by the CCD. Moreover, up to a flipping around the optical axis, each location ${\boldsymbol{p}}$ on the object, at a distance $\tau$ from the optical axis $O$ (dashed line), is probed by a corresponding CCD pixel also at a distance of $\tau$ from $O$. Each location $\boldsymbol{p}$ is thus in one-to-one correspondence with a CCD pixel and we will sometimes consider $\boldsymbol{p}$ as CCD pixel location.

From classical optics, a light ray that is incident on location ${\boldsymbol{p}}$ at an angle ${\theta}_{\boldsymbol{p}}$ originates from the light source at a distance of $\Delta x={f}\tan{\theta}_{\boldsymbol{p}}$ from the optical axis. Likewise, the light rays originating at different locations on the source have different incident angles at ${\boldsymbol{p}}$. Since we can always virtually invert the light propagation in the system, everything works as if the object was shined on its right by a beam of parallel light rays. Therefore, up to a global scaling by $f$, the SLM plane is actually the local plane $\Pi_{\boldsymbol{p}}$ of the deflection spectrum occurring at $\boldsymbol{p}$. Modulating the SLM corresponds to modulating $s_{\boldsymbol{p}}$, while the light collected in CCD pixel $\boldsymbol{p}$ is just an inner product of $s_{\boldsymbol{p}}$ with the modulation.

If we generate $M$ such modulations $\phi_{i}\in\mathbb{R}^{N}$ with $1\leq i\leq M$ in the SLM of $N$ pixels, considering the discrete nature of the CCD camera (having $N_{\rm C}$ pixels), the discretized deflection spectra are observed through

$$\boldsymbol{y}_{k}=\boldsymbol{\Phi}\boldsymbol{s}_{k}+\boldsymbol{n},\quad 1\leq k\leq N_{\rm C},$$

(1)

where $\boldsymbol{\Phi}^{T}=(\boldsymbol{\phi}_{1},\cdots,\boldsymbol{\phi}_{M})\in\mathbb{R}^{N\times M}$ is the sensing matrix, $k$ is a CCD pixel index, $\boldsymbol{s}_{k}\in\mathbb{R}^{N}$ is the discretized spectrum at the $k^{\rm th}$ pixel/object location, and $\boldsymbol{n}$ models the measurement noise (assumed Gaussian). Notice that the SLM and the CCD 2-D grids are represented as 1-D spaces in order to ease the notation, so that $\boldsymbol{\Phi}$ is then a sensing 2-D matrix acting on 1-D vectors.

The quest now is to optimize the design of $\boldsymbol{\Phi}$ in order to maximize the spectrum information captured in each $\boldsymbol{y}_{k}$, i.e., using Compressed Sensing theory.

A well known example of inner products based imaging system is the Rice University’s single pixel camera [5]. This camera uses a single photosensor to capture the inner products between the scene to be imaged and random binary modulation patterns, realised using a micro-mirror array. In our deflectometric system, each CCD pixel behaves like a single pixel camera, which poses a great computational challenge for reconstruction. We can notice also that a binary SLM modulation for compressive imaging is also advocated in [6] for performing random convolutions of images. However, this method modulates the Fourier transform of the signal of interest and not its spatial domain representation.

The previous use of Schlieren Deflectometer is in the Phase Shifting Schlieren (PSS) system. PSS quantitively measures the deviation angles from the deflectometric measurements, by assuming single deflection for each object location. It overcomes several calibration issues of the traditional deflectometry by using multiline phase shifted spatial filters in the SLM [7]. In this case, the CCD pixel values $\boldsymbol{y}_{k}$ directly encode the deviation angle of the corresponding object location $k$. By making several measurements with phase shifted sinusoidal modulations as $\boldsymbol{\Phi}$, the deviation angle is then numerically decoded from the $\boldsymbol{y}_{k}$, using $n$-step algorithms [8]. PSS has been successfully used in tomographic applications such as refractive index map reconstruction [9, 10]. However, the phase shifting method is unable to recover the full deflection spectra or to estimate several main deflection angles per point. Moreover, the use of non-binary modulation patterns brings in the problem of non-linearity in the SLM response. This is one of the reasons why we advocate using binary modulation patterns as described in Sec. 5.

Compressive Sensing (CS) shows that any signal $\boldsymbol{x}=\boldsymbol{\Psi}\boldsymbol{\alpha}\in\mathbb{C}^{N}$, having a sparse representation in an orthonormal sparsity basis $\boldsymbol{\Psi}\in\mathbb{C}^{N\times N}$, i.e., $\|\boldsymbol{\alpha}\|_{0}:=\#\{j:\alpha_{j}\neq 0\}\leq K\ll N$, can be tractably recovered from a few corrupted linear measurements of the form

$$\boldsymbol{y}=\boldsymbol{\Gamma}_{\Omega}^{\ast}\boldsymbol{\Psi}\boldsymbol{\alpha}\;+\;\boldsymbol{n},$$

(2)

where ^∗ denotes the conjugate transpose, $\boldsymbol{\Gamma}\in\mathbb{C}^{N\times N}$ is an orthonormal sensing basis, $\boldsymbol{\Gamma}_{\Omega}$ is the submatrix formed by restricting the columns of $\boldsymbol{\Gamma}$ to those belonging to $\Omega\subset[N]:=\{1,\,\cdots,N\}$ and $\boldsymbol{n}$ is a Gaussian noise vector.

In particular, if $\Omega$ has size $M$ and is drawn uniformly at random in $[N]$, and if $\boldsymbol{\Gamma}$ and $\boldsymbol{\Psi}$ are incoherent, meaning that the coherence $\mu:=\sqrt{N}\max_{1\leq i,j\leq N}|\langle\boldsymbol{\Gamma}_{j},\boldsymbol{\psi}_{i}\rangle|$ is very close to 1, then,

$$M=O\big{(}\mu^{2}K\log^{4}(N)\big{)}$$

measurements are enough to recover a good estimate of $\boldsymbol{x}$ [11]. This reconstruction is performed by solving the following convex optimization

$$\widehat{\boldsymbol{\alpha}}:=\operatorname*{\arg\min}_{\widetilde{\boldsymbol{\alpha}}\in\mathbb{C}^{N}}\;\|\widetilde{\boldsymbol{\alpha}}\|_{1}\;\text{subject to}\;\|\boldsymbol{y}-\mathbf{\Phi}\widetilde{\boldsymbol{\alpha}}\|_{2}\leq\epsilon,$$

(3)

where $\boldsymbol{\Phi}=\boldsymbol{\Gamma}^{*}_{\Omega}\boldsymbol{\Psi}$, and $\epsilon$ is a bound on $\|\boldsymbol{n}\|_{2}\leq\epsilon$. Under the previous conditions, the theory guantees that [12]

$$\|\boldsymbol{\alpha}-\widehat{\boldsymbol{\alpha}}\|_{2}=O\big{(}\tfrac{\|\boldsymbol{\alpha}-\boldsymbol{\alpha}_{K}\|_{1}}{\sqrt{K}}\;+\;\epsilon\big{)}$$

(4)

holds with a probability at least $1-N^{-\gamma\log^{3}(N)}$, where $\boldsymbol{\alpha}_{K}$ is the best $K$-term approximation of the vector $\boldsymbol{\alpha}$.

As shown above, the number of the measurements scales quadratically as $\mu$ and hence it is desirable to have the two bases as incoherent as possible to make $\mu$ close to $1$. In order to circumvent this problem, [2] cleverly makes the sensing and sparsity bases incoherent by a simple pre-modulation, or spread spectrum, of the data vector $\boldsymbol{x}$. This technique is optimal for universal sensing basis, i.e., when all the entries of $\boldsymbol{\Phi}$ have the same complex amplitude, which is the case for Fourier and Hadamard bases.

Mathematically, with random modulation, the sensing matrix becomes $\boldsymbol{\Phi}=\boldsymbol{\Gamma}_{\Omega}^{\ast}\boldsymbol{M}\boldsymbol{\Psi}$, where $\boldsymbol{M}=\operatorname{diag}(\boldsymbol{m})$ is a diagonal matrix whose diagonal is a random vector $\boldsymbol{m}$ such that $|m_{i}|=1$, e.g., a Steinhaus or Rademacher sequence. In this case, we need

$$M\geq C_{\rho}\,K\log^{5}(N)$$

measurements in order to recover a solution $\boldsymbol{\alpha}^{\star}$ of (3) satisfying (4) with a probability at least $1-O(N^{-\rho})$, for some $0<\rho<\log^{3}(N)$. Noticeably, the coherence has disappeared from the condition implying that with spread spectrum and universal sensing basis, the recovery guarantee is universal, irrespective of the sparsity basis. We see in next section how to exploit the spread spectrum method in our optical setup.

We propose to combine the Schlieren deflectometer of Sec. 2 with the spread spectrum CS principles described above.

Since the SLM accepts only positive values, once the random set $\Omega\subset[N]$ has been drawn, the corresponding sensing matrix $\boldsymbol{\Phi}=\boldsymbol{H}^{T}_{\Omega}\boldsymbol{M}$ must be suitably biased and scaled as

$$\overline{\boldsymbol{\Phi}}=\tfrac{1}{2}(\boldsymbol{\Phi}+\mathbf{1}_{N}\mathbf{1}_{N}^{T})\in\{0,1\}^{M\times N},$$

(5)

where $\mathbf{1}_{N}\in\{1\}^{N}$ is the vector of ones. The choice of Hadamard makes the biased sensing binary. This has interesting advantages in the optical setup since it avoids any non-linear response of the SLM for non-binary values.

Moreover, thanks to an extra measurement of each ${\boldsymbol{s}}_{k}$, i.e., the value $\overline{z}_{k}=\langle\boldsymbol{1}_{N},\boldsymbol{s}_{k}\rangle$ obtained by an entirely transparent SLM, any measurement vector $\overline{\boldsymbol{y}}_{k}=\overline{\boldsymbol{\Phi}}\boldsymbol{s}_{k}$ can be debiased by

$$\boldsymbol{y}_{k}:=\mathbf{\Phi}{\boldsymbol{s}}_{k}=2\overline{\boldsymbol{y}}_{k}-\overline{z}_{k}\,\mathbf{1}_{M}.$$

(6)

Up to a correct evaluation of the noise, which corrupts $\overline{\boldsymbol{y}}_{k}$, ${\boldsymbol{z}}_{k}$ and the debiasing process (6), we can therefore reconstruct ${\boldsymbol{s}}_{k}$ from the unbiased $\boldsymbol{y}_{k}$.

From actual measurements in the absence of test object, we obtain, on an arbitrary CCD pixel, $\boldsymbol{y}^{\rm no}=\mathbf{\Phi}({\boldsymbol{s}}^{\rm no}+\boldsymbol{n}_{{\boldsymbol{s}}})+\boldsymbol{n}_{\boldsymbol{y}}$, where $\boldsymbol{n}_{{\boldsymbol{s}}}$ and $\boldsymbol{n}_{\boldsymbol{y}}$ are the unknown signal and observation noises. After a small calibration of the SLM origin, and up to a small optimization of the disk height in ${\boldsymbol{s}}^{\rm no}$, we can therefore compute a bound on the noise power as $\epsilon=\|\mathbf{\Phi}\boldsymbol{n}_{{\boldsymbol{s}}}+\boldsymbol{n}_{\boldsymbol{y}}\|=\|\boldsymbol{y}^{\rm no}-\mathbf{\Phi}{\boldsymbol{s}}^{\rm no}\|$. We can either obtain this value for every $M$ or estimate it for $M=N$ only and scale the result as $\epsilon(N)=\sqrt{M+2\sqrt{M}}\ \epsilon(M)/\sqrt{N}$ for $M<N$.

For evaluating compressive reconstruction performance, (3) was solved with $M=N$ measurements to obtain the reference reconstruction $\tilde{{\boldsymbol{s}}}_{k}$. Reconstructions for $M<N$ were compared with $\tilde{{\boldsymbol{s}}}_{k}$ using the (output) Signal-to-Noise Ratio $\text{oSNR}:=20\log_{10}(\|\tilde{{\boldsymbol{s}}}_{k}\|_{2}/\|\tilde{{\boldsymbol{s}}}_{k}-\hat{{\boldsymbol{s}}}_{k}\|_{2})$.

Fig. 3(left) shows a deflection spectrum reconstructed using the full set of measurements for the lens with $60D$ optical power and Fig. 3(right) shows the reconstruction with $400$ measurements ($M/N\approx 10\%$). Note that the spectrum is well localized, corroborating our initial observation.

Fig. 4 shows the plot of oSNR versus the number of measurements $M/N$ (in %), averaged over the trials and locations. The performance curves for both the lenses almost overlap, and the oSNR improves as $M/N$ increases. Though the absolute values of oSNR seem low, its significance has be understood in the light of the input SNR, which is approximately computed as $\text{iSNR}\ :=20\log_{10}(\|\mathbf{\Phi}{\boldsymbol{s}}^{\rm no}\|_{2}/\|\boldsymbol{y}^{\rm no}-\mathbf{\Phi}{\boldsymbol{s}}^{\rm no}\|_{2})\simeq 4.34$ dB. The horizontal dotted line on the plot indicates the iSNR for our experiments, and it is clear that the reconstruction procedure improves the oSNR, beyond the iSNR, thereby demonstrating the ability of CS reconstruction of deflection spectra in low input SNR regime.

Suppose we are interested in summarizing the deflection spectra by a few parameters that characterize them at each location. For example, we could summarize a spectrum by the location of the bright feature it contains. To address this, we propose the following robust method to localize the feature.

Assuming that the feature is circular and its radius $\rho$ can be known a priori, the feature can be localized using a matched filter. To this end, we employ a template $\mathbf{G}_{\rho}$, which contains a two-dimensional Gaussian of width $\rho$ at its center. To localize the feature in a spectrum $\tilde{{\boldsymbol{s}}}_{k}$ we simply translate $\mathbf{G}_{\rho}$ over the support of the spectrum and find the best translation $\hat{\tau}\in{\mathbb{R}}^{2}$ that maximizes the correlation between the spectrum and the translated template $\mathcal{T}_{{\tau}}\mathbf{G}_{\rho}$. By letting ${\boldsymbol{g}}^{\rho}_{\tau}$ to denote the vectorized form of $\mathcal{T}_{{\tau}}\mathbf{G}_{\rho}$, we solve

$$\hat{\tau}_{k}=\arg\max_{\tau}\,|\langle\tilde{{\boldsymbol{s}}}_{k},{\boldsymbol{g}}^{\rho}_{\tau}\rangle|.$$

(7)

With an abused terminology we designate the translation parameters $\hat{\tau}$ as the centroids of the spectral features. These centroids also provide local information about the mean deflection angles at various locations.

In order to characterize an object using the centroids defined in (7), we have to reconstruct the spectra at all $N_{\rm C}$ locations individually and then compute the centroids, which requires huge computational time. It essential therefore to explore the possibility of computing the centroids right from the measurements $\boldsymbol{y}_{k}$, without even reconstructing the spectra.

Davenport et al. have proposed the idea of performing common signal processing tasks such as detection and classification using the compressive samples obtained using random measurements [3, 4]. In the same spirit, we compute the centroids directly from the measurements $\boldsymbol{y}_{k}$ defined in (6). To accomplish this compressive centroid estimation, we carry the definition of the centroid of the spectrum (7) into the measurement domain and solve

$$\widetilde{\tau}_{k}=\arg\max_{\tau}\,|\langle\mathbf{\Phi}^{T}\boldsymbol{y}_{k},{\boldsymbol{g}}^{\rho}_{\tau}\rangle|.$$

(8)

Notice that solving (8) is the same as solving (7), but instead of using $\tilde{{\boldsymbol{s}}}_{k}$, obtained by solving (3), we simply compute $\mathbf{\Phi}^{T}\boldsymbol{y}_{k}$ and use it for centroid estimation.

For experimental evaluation, we retained the configuration described in Sec. 5. For each of the 5 CCD locations $k$, the centroid $\widetilde{\tau}_{k}$ was computed by solving (8). A “ground truth” centroid $\hat{\tau}_{k}$ was also found by solving (7), with ${\boldsymbol{s}}_{k}$ reconstructed (solving (3)) using full $M=4096$ measurements.

Fig. 5 shows the centroid computation error $\|\hat{\tau}-\tilde{\tau}\|_{2}$ as a function of the number of measurements $M/N$. Each data point is an average over 50 independent trials for each value of $M$ and over all the locations. The horizontal dotted line indicates a unit pixel error and it can be seen that the compressive centroid estimation achieves sub-pixel accuracy, even with the number of measurements as low as $2.4\%$ ($50$ measurements) for the $60$D lens, and about $3.7\%$ for the $10.03$D lens. This demonstrates the ability of the measurements to capture sufficient information about the feature of the spectra.

This paper presents a novel approach to perform schlieren deflectometry using compressive sensing principles. In contrast to the existing system, the proposed approach enables the reconstruction of deflection spectra, instead of a single deflection angle at each location. Although it involves computationally intensive procedure, spectra reconstruction leads to better characterization of the studied object. The empirical results also demonstrate the capability of our method to perform meaningful inferences about the deflection spectra using only the compressive measurements, without involving any reconstruction.

Several challenges remain to be addressed in the future work. It is of foremost importance to fully understand the noise model to tune the reconstruction algorithm. Also, incorporating further priors such as non-negativity helps the reconstruction. We also intend to develop approaches to simultaneously reconstruct spectra, corresponding to neighbouring CCD pixels, by exploiting the correlation amongst them.