A new wavelet-based reconstruction algorithm for twin image removal in digital in-line holography

A new wavelet-based reconstruction algorithm for twin image removal in digital in-line holography

Optics and Lasers in Engineering 82 (2016) 159–172 Contents lists available at ScienceDirect Optics and Lasers in Engineering journal homepage: www...

4MB Sizes 129 Downloads 271 Views

Optics and Lasers in Engineering 82 (2016) 159–172

Contents lists available at ScienceDirect

Optics and Lasers in Engineering journal homepage: www.elsevier.com/locate/optlaseng

A new wavelet-based reconstruction algorithm for twin image removal in digital in-line holography Jamel Hattay a,c,n, Samir Belaid b, Taoufik Aguili a, Denis Lebrun c a

Communication System Laboratory Sys'Com, National Engineering School of Tunis, University Tunis El Manar, Tunisia University of Monastir, Tunisia c UMR 6614 CORIA, Normandy University, Avenue de l'Université, 76801 Saint-Etienne du Rouvray, France b

art ic l e i nf o

a b s t r a c t

Article history: Received 4 November 2015 Received in revised form 10 February 2016 Accepted 11 February 2016

Two original methods are proposed here for digital in-line hologram processing. Firstly, we propose an entropy-based method to retrieve the focus plane which is very useful for digital hologram reconstruction. Secondly, we introduce a new approach to remove the so-called twin images reconstructed by holograms. This is achieved owing to the Blind Source Separation (BSS) technique. The proposed method is made up of two steps: an Adaptive Quincunx Lifting Scheme (AQLS) and a statistical unmixing algorithm. The AQLS tool is based on wavelet packet transform, whose role is to maximize the sparseness of the input holograms. The unmixing algorithm uses the Independent Component Analysis (ICA) tool. Experimental results confirm the ability of convolutive blind source separation to discard the unwanted twin image from in-line digital holograms. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Twin image removal Digital holography Focus plane Convolutive blind source separation Wavelet transform

1. Introduction Unlike classical holography [1], digital holography [2] makes possible the numerical processing of digital holograms, such as the focusing at variable depths of 3D objects, tracking, localizing and estimating the sizes of particles in 3D volume [3]. The in-line configuration of digital holography uses a simple instrumental scheme, where the object and reference beams are the diffracted and non-diffracted parts of the same beam. This configuration has proved its effectiveness in the study of particles, aerosols, fluid-flow analysis [4], coherent diffractive imaging with X-ray holography [5], atomic imaging with complex γ -ray holography [6] and size recovering of nanoparticles with cavitation air bubbles [7]. However, both zero-order and twin images are serious obstacles to the retrieval of the reconstructed image with high quality. These parasitic zero-order and twin images are due respectively to the acquisition process, where different intensities of both waves and the object wavefront information are recorded, and to the complex conjugate of the object wave. This redundant information generates unwanted fringes that make difficult an accurate reconstruction of the image, and can even heavily disturb the recovering of its microphysical properties. n

Corresponding author. E-mail address: [email protected] (J. Hattay).

http://dx.doi.org/10.1016/j.optlaseng.2016.02.009 0143-8166/& 2016 Elsevier Ltd. All rights reserved.

In the last decade, many techniques were developed to solve the twin image removal problem. Usually, this problem was processed by both spatial carried-frequency approaches and numerical in-line geometry means [8]. In [9], a suitable algorithm was developed for simultaneous measurement of refractive index and thickness which enabled the removal of the unwanted twin image. This tool involves the dual-wavelength diffraction phase microscopy based on the spatial filter as the form of a common-path imaging interferometer. In [10,11], various advanced spatial carrier phase-shifting (SCPS) manipulations were described for the removing of the twin image on the reconstructed hologram. In other respects, many algorithms were proposed to track the unwanted fringes based on the principal-vector-directed fringetracking technique [12]. Using the Gaussian derivatives to estimate the fringes gradients, this technique involves some thresholding hysteresis for singular points segmentation. By means of some deconvolution techniques [14], the authors of [13] propose an iterative algorithm to discard the unwanted fringes from the reconstructed images. Based on wavelet transform theory, the authors of [15] proved the perfect suppression of twin image by angular spectrum diffraction at the ridge of Gabor wavelet transform (GWT). In our preview work, we presented a robust twinimage removal method based on blind source separation (BSS), combined with an efficient quincunx lifting scheme [16]. Convolutive blind source separation (BSS), also known as blind source deconvolution, consists in separating a set of source signals from a set of convolutely mixed signals without any information

160

J. Hattay et al. / Optics and Lasers in Engineering 82 (2016) 159–172

about the source signals and the mixing process. The inputs of BSS unmixing process are the mixed images received over Q sensors (Q observed images or Q observations), and the outputs are P estimated images (i.e. estimated source images). Several methods have been proposed to solve convolutive BSS, among which is our contribution based on sparse and multiscale representation of image mixtures [17,18]. In this paper, we present our enhanced BSS-based approach for twin image removal. As explained in [14], the mathematical model of digital Fresnel holography shows that the image in the reconstructed plane can be written as a convolution formalism of different physical functions involved in digital holography. In this context, these physical effects, that are the image of the object, the Fresnel kernel and the twin image, convolutely mixed, form the observed images in the reconstructed plane. The input of our deconvolution algorithm is the hologram fringe pattern, the real part and the imaginary part of the reconstructed hologram to obtain the original image of the object. The paper is organized as follows. In Section 2, a brief overview of hologram recording and reconstruction of in-line holograms is given. In Section 3, we present a new method for focus plane retrieving based on an entropy criterion. In Section 4, we emphasize the concept of sparse image representations as a preprocessing for the blind source separation task. In Section 5, we describe the proposed reconstruction algorithm and its underlying sparseness measure. In Section 6.2, we illustrate and validate the proposed method using various real and simulated reconstructed holograms. Finally, Section 7 contains the summary and concluding remarks.

Fig. 1. Recording system of digital in-line holograms ðξ; ηÞ: object plane. (x,y): sensor plane.

2. Digital in-line holography 2.1. Digital hologram recording In the detector and reconstruction planes, in-line holography uses the convolution formalism for mathematic model of recorded and reconstructed images. We note by 1  Tðξ; ηÞ the binary amplitude distribution of an opaque object illuminated by a monochromatic plane wave and located at a distance ze from a camera as shown in Fig. 1. As explained in [19], the particles are

Fig. 3. Mixing and unmixing system.

Fig. 4. Illustration of the in-line holography process as convolutive BSS problem: X1: intensity function of hologram, X2: real part of the reconstructed image and X3: imaginary part of the reconstructed image.

Fig. 2. QLS decomposition of digital hologram at J¼4.

J. Hattay et al. / Optics and Lasers in Engineering 82 (2016) 159–172

161

2.2. Digital hologram reconstruction The amplitude distribution in the reconstructed plane located at a distance zr from the camera can be calculated by a convolution operation: Rðx; yÞ ¼ ½I ze nnhzr ðx; yÞ

ð6Þ

In order to weaken the influence of the high spatial frequencies (generally undersampled), the function hzr ðx; yÞ is often windowed by a Gaussian function. In our laboratory, Rðx; yÞ is either computed by applying a wavelet transformation or a Fractional Fourier Transformation [21,22]. Here, we note that the diffraction kernel has not been filtered and the convolution operation described by Eq. (6) has been directly applied to the holograms (Fresnel transform). Using Eqs. (4) and (6) and assuming that the right focus plane is selected (i.e. ze ¼ zr ¼ z), the complex reconstructed image is simply rewritten as Rðx; yÞ ¼ 1  Tðx; yÞ  ½T nnh2z ðx; yÞ

ð7Þ

The simplifications leading to Eq. (7) are due to the fact that hz ðx; yÞnnh z ðx; yÞ ¼ δðx; yÞ and hz ðx; yÞnnhz ðx; yÞ ¼ h2z ðx; yÞ, where δðx; yÞ is the Dirac distribution [20] and the last term in the right hand side of this equation represents the twin image.

3. Autofocus plane using entropy criterion

Fig. 5. General block diagram of BSS-based reconstruction algorithm.

often approximated as planar opaque objects (2D) since this approximation is valid for small objects located far away from the sensor. From the Huygens–Fresnel integral, the complex amplitude Aze in the detector plane is calculated by 2-D convolution operation [20]. Aze ðx; yÞ ¼ ½ð1  TÞnhze ðx; yÞ

ð1Þ

where

  iλπz ðx2 þ y2 Þ 1 e hze ðx; yÞ ¼ e iλze

ð2Þ

From this complex amplitude, the recorded intensity in the sensor plane is given as follows. I ze ðx; y; ze Þ ¼ Aze ðx; yÞA ze ðx; yÞ

ð3Þ

By means of Eqs. (1) and (3), I ze ðx; yÞ can be expressed as follows: I ze ðx; y; ze Þ ¼ 1  ½T nnðhze þ h ze Þðx; yÞ

ð4Þ

Note that, in Eq. (4), there is the term ðj ½T nnhze ðx; yÞj Þ that has been omitted. This approximation is valid, provided that the Farfield domain is considered [21]. We obtain therefore:   2 π  2 2 sin x þy : ð5Þ I ze ðx; y; ze Þ ¼ 1  Tðx; yÞnn λze λze 2

In many practical cases, the distance between the object and the measurement plane is unknown. In order to correctly reconstruct the image, it is necessary to know this information. Some methods involve the mutual information to find the focus plane of the reconstructed holograms. The entropy of the reconstructed image is calculated in the optimal distance between the object and the detector plane usually noted by reconstruction distance zr [23]. According to the mutual information theory [24], the entropy has been especially used as a quantitative measurement for the dispersion of joint distribution and the degree of independence between the variables. In this context, the author of [25] shows that the entropy variation of the reconstructed hologram is symmetric about the focus plane. In this plane, the image of the object is perfectly reconstructed. Furthermore, as shown by the authors of [26], the variations of the imaginary part are theoretically minimum in the image plane. Consequently, when the best focus plane is selected, the real part's entropy is minimal and the entropy of the imaginary part is maximal. Several works use a sparse hologram representation to measure the focus plane [27]. The focus plane has been improved by canceling the autocorrelation of the input hologram owing to multiresolution analysis through the wavelet packet transform [28]. Motivated by the latter work, we propose an entropy-based method to estimate the focus plane. Firstly, in order to provide a sparse representation, the real/imaginary parts of the reconstructed image are decomposed by an extended wavelet transform also called Quincunx Lifting Scheme (QLS) [29]. Recall that, as shown in Fig. 2, the QLS decomposes the original hologram into approximation subband at level J noted by aJ and several detail subbands at resolution levels 1; …; J [16]. Secondly, the wavelet coefficients of the imaginary part, obtained by the latter decomposition, are modeled by a suitable Generalized Gaussian Distribution (GGD) to find the optimal entropy that corresponds to the focus plane. Thus, the weighted entropy of imaginary part of reconstructed hologram HJ is the weighted sum of the entropy HðaJ=2 Þ of the coarsest approximation and the entropies Hðdj=2 Þ of the detail sub-

162

J. Hattay et al. / Optics and Lasers in Engineering 82 (2016) 159–172

Fig. 6. Intensity functions of the simulated holograms.

Fig. 7. Focus plane calculation of simulated small particle (estimated ze ¼ 120 mm) and big particle (estimated ze ¼ 200 mm) holograms.

images at every resolution levels j ¼ 1; …; J: HJ ¼

1 2

J

HðaJ Þ þ

J X 1 j¼12

j

Hðdj=2 Þ

ð8Þ

Recall that the entropy of a discrete random variable X taking values in ω with probability distribution p(X) is given by [30] X HðXÞ ¼  pðxÞ log pðxÞ xAω

As explained in [31], the details subbands are seen as a continuous random variable ξ whose realizations are the coefficients of each subband. The distribution of ξ, denoted by f, is modelled as a zero-mean Generalized Gaussian Density (GGD) function expressed as follows: 8 ξ A R;

βj=2 βj=2 e  ðj ξ j =αj=2 Þ f ðξÞ ¼ 1 2αj=2 Γ ðβj=2 Þ R þ1

ð9Þ

e dt [32]. Either the moments method where Γ ðzÞ ¼ 0 t [33] or the Maximum Likelihood (ML) method [34] can be used to estimate the scale parameter αj=2 and the shape parameter βj=2 . Using the approximation proposed in [35], the entropy is z1 t

equivalent to the differential entropy given by the following equation:   ~ αj=2 ; β Þ ¼  E½log f ðξÞ  Hð j=2 By returning to Eq. (9), the approximated differential entropy becomes 1 ! 2αj=2 Γ ðβ j=2 Þ 1 ~ αj=2 ; β Þ ¼ log Hð ð10Þ þ βj=2 j=2

βj=2

The global entropy of the imaginary part HJ given in Eq. (8) could be approximated as follows: HJ C

1 2

J

HaðJÞ=2 þ

J X 1 ~ ðαj=2 ; β j=2 Þ H j dðjÞ=2 2 j¼1

ð11Þ

From Eq. (11), we deduce that the QLS operators pj=2 defined in [36] that minimize the global entropy HJ also minimizes ~ d ðαj=2 ; β Þ and then maximizes the maximum likelihood H j=2 ðjÞ=2 Lðpj=2 ; αj=2 ; β j=2 Þ ¼

K j=2 X

  log f ðξ~ j=2 ðkÞ  ξj=2 ðkÞ pj=2

ð12Þ

k¼1

where ξ~ j=2 ð1Þ; …; ξ~ j=2 ðK j=2 Þ and ξj=2 ð1Þ; …; ξj=2 ðK j=2 Þ are realizations and ξ . of ξ~ j=2

j=2

J. Hattay et al. / Optics and Lasers in Engineering 82 (2016) 159–172

163

Fig. 8. BSS algorithm applied to the simulated hologram of small particle reconstructed at ze ¼ 120 mm.

later in Section 6.2. They confirm that the use of the proposed entropy model is a useful technique to measure the focus plane regardless of the effect of instruments errors and additive noise in the hologram recording phase.

Finally, the ML estimator of pj=2 can be expressed as X K j=2

Lðpj=2 ; αj=2 ; βj=2 Þ ¼

k¼1

1

log ðβ j=2 Þ  log ð2αj=2 Γ ðβj=2 Þ "

 exp β j=2 log

j ξ~ j=2 ðkÞ  ξj=2 ðkÞpj=2 j

αj=2

!# ð13Þ

and minimizes the following ℓβj=2 criterion that will be mainly used to optimize the entropy of imaginary part of the reconstructed hologram: ℓ

βj=2

ðpj=2 Þ ¼

K j=2 X

j ξ~

ξ

j=2 ðkÞ  j=2 ðkÞpj=2 j

βj=2

ð14Þ

k¼1

In our case, the Gauss–Newton tool is used as an iterative algorithm to minimize this criterion [37]. Next, we define an interval of z supposed to contain the axial coordinate of the focus plane. For each value of this interval and using a very small pitch, we calculate successively the optimal entropy of the imaginary part of the reconstructed image at this distance. As a result, we obtain a set of entropy values versus the reconstruction distance z where the maximum of this set corresponds to the axial coordinate of optimal focus plane. The performance of this analysis, approved by different simulated hologram where the axial coordinate z is known, will be presented

4. Blind Source Separation (BSS) for twin image removal 4.1. Convolutive formalism of in-line digital holography In the Blind Source Separation BSS process, both the sources and the mixing process are unknown and only the recording mixtures are available. The goal is to recover the unknown sources from the observed mixtures. As illustrated in Fig. 3, the general model considers P source signals sp(n), p : 1; …; P and Q observations’ signals xq(n), q : 1; …; Q generated with the mixing system M. The linear and instantaneous model supposes that the coefficients M q;p are simple scalar values. The equation of the mixing process can be given as follows: xq ðnÞ ¼

P X

M q;p sp ðnÞ þ θðnÞ

ð15Þ

p¼1

where the additive noise at position n is denoted by θðnÞ. However, in real situations, such as in image processing, another difficulty

164

J. Hattay et al. / Optics and Lasers in Engineering 82 (2016) 159–172

Fig. 9. BSS algorithm applied to the simulated hologram of big particle reconstructed at ze ¼200 mm.

Fig. 10. Profile comparison applied to simulated hologram of small particle hologram.

cannot be omitted. In this case, all source pixels are weighted and delayed and the observations recorded by the sensors are seen as a sum with multiple delays of propagation. This mixing process is referred to as the convolutive mixing process [38].

The 2D convolutive model assumes that the P images are denoted by sðn1 ; n2 Þ ¼ ðs1 ðn1 ; n2 Þ; ‥; sP ðn1 ; n2 ÞÞT where the doubleindex notation ðn1 ; n2 Þ designates a pixel position in the whole image. M q;p is replaced, in this situation, by an L  L matrix which

J. Hattay et al. / Optics and Lasers in Engineering 82 (2016) 159–172

165

Fig. 11. Intensity distributions of the experimental holograms recorded in various experiments.

designates the convolutive mask. The convolutive mixing model, without noise, uses the Q A N n observed images xðn1 ; n2 Þ ¼ ðx1 ðn1 ; n2 Þ; ‥; xQ ðn1 ; n2 ÞÞT as follows: xq ðn1 ; n2 Þ ¼

P L L X X X

M q;p ðk1 ; k2 Þsp ðn1  k1 ; n2  k2 Þ

p ¼ 1 k1 ¼ 1 k2 ¼ 1

Q L L X X X

Hologram

Description

Cavitation bubble

The unmixing process consists in inverting operation described above as indicated in Fig. 3. It involves an L  L-matrix denoted by W p;q . The unmixing process is given by the following equation: yp ðn1 ; n2 Þ ¼

Table 1 Recording features of tested holograms (d is a bubble diameter, ze is a recording distance, λ is a wavelength and p is the pixel size).

W p;q ðk1 ; k2 Þxq ðn1  k1 ; n2  k2 Þ

Bubble recorded in a cavitation tunnel Real large particle Big bubble recorded in experimental project Calibration target Calibration target of opaque disks and bars Irregular object Unknown object carried in a cavitation tunnel

Recording parameters d ¼120 μm ze ¼90 mm λ ¼635 nm p ¼ 6.7 μm d ¼517 μm ze ¼ 117 mm λ ¼488 nm, p¼ 5.5 μm ze ¼117 mm λ ¼488 nm p ¼5.5 μm ze ¼67 mm λ ¼635 nm p ¼6.7 μm

q ¼ 1 k1 ¼ 1 k2 ¼ 1

where y1 ðn1 ; n2 Þ; …; yP ðn1 ; n2 Þ denote the estimated image sources. In this paper, the convolutive model of BSS is applied to solve the twin image problem in the in-line holography configuration. As described in Fig. 4, the number of image sources and observations is reduced to three: the original Hologram and both the real and imaginary parts of the reconstructed image at the optimal reconstruction distance zr ¼ zopt . According to Eq. (7), the reconstructed image is not perfect because it is surrounded by an additive complex noise due to the unwanted fringes corresponding to the term T nnh2z ðx; yÞ. By resorting to a BSS-based unmixing algorithm, we show how the source function 1  Tðx; yÞ, which represents the improved reconstructed image, is extracted from the three observed images indicated above (i.e. the hologram, the real and imaginary parts of

the reconstructed image given by Eq. (7). The second and third sources are the unwanted twin image function and the Fresnel kernel functions. 4.2. Adaptive Quincunx Lifting Scheme The Spatial Lifting Scheme is considered as an extension of classic wavelet transforms (LS) [39]. LS is made up of three steps. To a better understanding of this scheme, we start by explaining it in the case of a one dimensional signal. The splitting step decomposes, spatially, the input signal into two sample sets S 1 and S 2 . The second step is the prediction step where the S 1 samples

166

J. Hattay et al. / Optics and Lasers in Engineering 82 (2016) 159–172

Fig. 12. Entropy calculation of the optimal axial coordinate of the focus plane. (a) Focus Plane Calculation of Cavitation Bubble (small bubble) (Estimated ze ¼ 90mm) and Big Bubble (Estimated ze ¼ 117mm) holograms.(b) Focus Plane Calculation of Calibration Target (Estimated ze ¼ 117mm) and Irregular object (Estimated ze ¼ 67mm) holograms.

are predicted from the S 2 samples. Finally, the update step consists in smoothing the S 2 samples using the residual of prediction error. This process is repeated recursively on the resulting updated coefficients to provide a representation of the input signal at several resolution levels. The residual prediction error at each level is known as the details subband and the updated set is considered as the approximations subband. To reconstruct the original signal, inverse LS reverses the above-described steps using the same update and prediction operators. When the signal is two dimensional, such as in image processing domain, the LS is applied separately on the rows then on the columns. Recently, the nonseparable lifting structures are combined with quincunx sampling of input image to obtain the Quincunx Lifting Scheme (QLS) version of the conventional LS [36,31]. In what follows, we explain the extension of the three LS steps, described above, for the QLS concept. Let us denote by x the input hologram (rectangularly sampled image). The quincunx sampling provides the following two polyphase components, instead of the S 1 and S 2 sets of classic LS. ( x1=2 ðn1 ; n2 Þ ¼ xðn1  n2 ; n1 þ n2 Þ ð16Þ x~ 1=2 ðn1 ; n2 Þ ¼ xðn1  n2 þ 1; n1 þ n2 Þ The prediction step consists in estimating x~ 1=2 samples from those of x1=2 . The prediction residual coefficients forming vector d1 are used to update the x1=2 samples into an approximation a1. Passing from resolution 1/2 to the next resolution can, therefore, be

expressed as follows: ( d1 ¼ x~ 1=2  ⌊x1=2 p1=2 ⌉ a1 ¼ x1=2 þ ⌊d1=2 u1=2 ⌉

ð17Þ

where

 x1=2 ðn1 ; n2 Þ is the vector containing the four neighboring samples of x1=2 ðn1 ; n2 Þ;

 d1 ðn1 ; n2 Þ is the reference detail vector and  p1=2 and u1=2 are respectively the prediction and the update vectors. More generally, for a multiresolution analysis, aj=2 will represent the approximation of the hologram at resolution j=2, jA Nn . The next step consists in splitting the input signal in order to generate two polyphase components xj=2 and x~ j=2 as follows: ( xj=2 ðn1 ; n2 Þ ¼ aj=2 ðn1  n2 ; n1 þ n2 Þ ð18Þ x~ j=2 ðn1 ; n2 Þ ¼ aj=2 ðn1  n2 þ 1; n1 þ n2 Þ Similarly, the prediction step consists in approximating the x~ j=2 samples from those of xj=2 ones. Finally, the residual error prediction coefficients forming vector dj=2 are used to update the xj=2 samples into the approximation aðj þ 1Þ=2 at the next stage. The described scheme is iterated J times in order to provide a multiresolution representation. The transition between resolution j=2

J. Hattay et al. / Optics and Lasers in Engineering 82 (2016) 159–172

167

Fig. 13. BSS algorithm applied to the experimental hologram of cavitation bubble reconstructed at ze ¼90 mm.

and the next resolution, i.e., ðjþ 1Þ=2, can be expressed as follows: (

dðj þ 1Þ=2 ¼ x~ j=2  ⌊xj=2 pj=2 ⌉ aðj þ 1Þ=2 ¼ xj=2 þ ⌊dj=2 uj=2 ⌉

ð19Þ

where

 xj=2 is the vector containing the four neighboring samples of xj=2 ;

 dj=2 is the residual prediction vector and  pj=2 and uj=2 are respectively the prediction and the update coefficient vectors. Accordingly, the final multiresolution representation of the hologram after a recursive decomposition over J stages is equivalent to the sum of the last approximation aJ=2 and all details dj=2 , 1 r jr J. Besides, the performances of the later decomposition (i.e. QLS) can be enhanced by the use of adaptive predictors [29]. The Adaptive QLS (AQLS) considers a block-based adaptation procedure which is coupled with a classified prediction approach. In the context of holography applications, the adaptive prediction approach uses several pairs of predictors and update operators. The proposed adaptation must depend on the objects that form the region of interest in the input hologram.

Note that, the AQLS is investigated and adapted for digital holography [16]. It was found as a powerful tool for removing the unwanted conjugate image (twin image) derived from the complex conjugate of object wave using an efficient method based on blind source separation approach [40]. 4.3. Frequency domain ICA Independent Component Analysis (ICA) is considered as an efficient BSS algorithm [41]. It provides very good results when the assumption of the statistical independence is verified. However, the Convolutive BSS, where the input images are mixed through a convolutive formalism, requires an extension of the classic ICA [42]. In this context, many algorithms were developed to solve this problem in time domain [43] or in transform domain (i.e. frequency domain) by means of the FFT transform [44]. The main contribution of this paper consists in using the wavelet packet transform (i.e. adaptive quincunx lifting scheme AQLS described in Section 4) to obtain the frequency complex coefficient instead of the classic FFT usually used in other works. The separation task is achieved by the frequency ICA algorithm dedicated to the transformed complex coefficients [45] which uses, as input data, the wavelet coefficients obtained by the AQLS tool. Notice that the frequency ICA algorithm requires a perfect decorrelation of the

168

J. Hattay et al. / Optics and Lasers in Engineering 82 (2016) 159–172

Fig. 14. BSS algorithm applied to the experimental hologram of big bubble reconstructed at ze ¼ 117 mm.

input images [46]. For this reason the wavelet packet transform, based on AQLS, is seen as a powerful transform to decorrelation purpose thanks to its fundamental statistical properties.

5. Proposed BSS-based hologram reconstruction The procedure starts by defining an interval I supposed to contain the axial coordinate of the best focus plane (i.e. zr ¼ze). Then, the reconstructed image is calculated through the Fresnel method, described in Section 2, in order to obtain the real and imaginary parts. The imaginary part of the reconstructed image is processed by the wavelet transform tool based on Adpative Quincunx Lifting Scheme (AQLS) exposed in Section 4.2. Applying this tool at a given level J results in a WP-tree whose nodes correspond to the decomposed holograms at several resolution levels as shown in Fig. 2. Then, we locate the node corresponding to the sparsest component in the WP-tree where the entropy is the decision criterion for sparseness. The latter tasks are repeated for all z distances of the interval I. For each value of zr, the entropy of the imaginary part, which minimizes the criterion ℓβ of Eq. (14), is calculated through the model given by the Eq. (8) in Section 3. The scan of this interval leads to the variation of the entropy of imaginary part versus the reconstruction distance zr. As explained

above, the maximum entropy value corresponds to the axial coordinate of the best focus plane ze. When the distance ze is estimated, the process of Fresnel reconstruction and AQLS transform are achieved for this focusing plane. Notice that, the hologram is also transformed by the same tool (i.e. AQLS). As indicated by Fig. 5, the next step consists in denoising the input hologram and its two reconstructed parts (real and imaginary) by thresholding the AQLS-details coefficients in order to estimate the standard deviation and variance of unwanted noise due to the recording environment and instruments. In the present study, the Stein Threshold, described in [47], is applied. After this step, these three transformed and denoised images, by AQLS, are used as an input of the unmixing algorithm based on frequency Independent Component Analysis (ICA) as detailed in Section 4.3. An inverse AQLS transform is carried out to retrieve the estimated and improved reconstructed image in the best focus plane, the unwanted twin image as the second source image in the convolutive formalism. The Fresnel kernel is supposed to be the third unknown source image convolutely mixed with the first source image (i.e. hologram) in the holography recording step. Algorithm 1 summarizes the overall proposed reconstruction process of digital hologram to remove the unwanted twin image. Note that the MATLAB codes can be obtained from the authors upon request.

J. Hattay et al. / Optics and Lasers in Engineering 82 (2016) 159–172

169

Algorithm 1. BSS-based reconstruction algorithm.

6. Experimental results 6.1. Simulated holograms We start by showing the behavior of the proposed reconstruction algorithm in case of simulated holograms. Fig. 6 shows two holograms that correspond respectively to small and big particles. Fig. 7 shows the entropy measurement of the focus plane applied to these two simulated holograms. As indicated by this figure, the axial coordinates ze are well estimated for both small and big particles.

The performance of our proposed algorithm has been tested also on the simulated holograms for twin image suppression. Figs. 8 and 9 indicate the results obtained by BSS-based algorithm applied respectively to simulated small and big particle holograms. It is clear in Figs. 8(d) and 9(d) that the twin image removal is successfully achieved where the particle is shown without any artifact due to the twin image. Fig. 10 shows the perfect reconstruction of simulated small particle in terms of radial intensity distribution. Particle profile is reconstructed without any artifacts of twin image.

170

J. Hattay et al. / Optics and Lasers in Engineering 82 (2016) 159–172

Fig. 15. BSS algorithm applied to the experimental hologram of calibration target reconstructed at ze ¼117 mm.

6.2. Experimental holograms Fig. 11 shows four different Regions Of Interest (ROI) cropped in digital holograms used in our experiments. Table 1 lists all recording features of these holograms. Notice that, in order to test the proposed algorithm with different objects, these holograms contain different shapes (bubble, disks, bars and irregular object). Note that the performances of the proposed approach are not compared with any existing BSS algorithms because these algorithms are not applicable in this way. It is worth mentioning that we have tested the implemented algorithms in NMFLAB platform [42], dedicated to blind source separation, but the separation task failed to discard the twin image. Firstly, we start by testing our entropy criterion to calculate the focus plane of tested holograms. Fig. 12 shows that the recording distances ze are well recovered for all tested holograms. As indicated by this figure, the optimal distance is retrieved as the maximum of entropy of imaginary part of reconstructed hologram because the signal variations are minimum in the focus plane (see [26]). It is clear that the proposed entropy-based method is suitable to these experimental holograms and efficient to any object shapes (bubble, bars, disks and irregular object). In order to test the reconstruction algorithm in a real situation, different experimental holograms have been processed. Figs. 13(d),

14(d), 15(d) and 16(d) show that the twin images are removed and the improved reconstructed images are retrieved without unwanted artifacts. These manipulations have two aims. Firstly, they show the behavior of the proposed BSS approach in several real situations with different physical parameters (particle diameter, recording distance, wavelength, etc.) and different shapes (bubble, disks, bars and irregular object). Secondly, they illustrate that our separation approach is powerful regardless of the impact of double artifact between two neighboring objects shown in Fig. 15. In this figure, the disks and bars are very close and the artifacts are overlapped. Note also that, in Fig. 14(d), the improved reconstructed image is not an opaque disk because in reality the experimentation involves a large transparent bubble that cannot be assumed as an opaque disk through the approximation demonstrated in [19]. This observation suggests that there is a need to use a more complete model that takes into account the refraction in the bubble . Fig. 17 shows a comparison between the radial intensity distribution of hologram cavitation bubble reconstructed by using the convolution of Eq. (6) and the estimated cavitation bubble function recovered by applying the proposed algorithm. The aim behind this figure is to objectively evaluate the proposed twin image removal algorithm. As we can see here, the profile of the improved image is reconstructed without any fluctuations whereas the real

J. Hattay et al. / Optics and Lasers in Engineering 82 (2016) 159–172

171

Fig. 16. BSS algorithm applied to the experimental hologram of irregular object reconstructed at ze ¼ 67 mm.

Fig. 17. Profile comparison between the real part of reconstructed hologram of cavitation bubble and the improved reconstructed bubble by BSS.

part of the reconstructed image is corrupted by the fluctuations coming from the unwanted twin image.

7. Summary and concluding remarks In digital holography, the in-line architecture offers significant advantages over its off-axis counterpart, apart the main problem of twin image. In this paper, we have succeeded to elaborate a

focus plane estimation method based on the entropy calculation of the imaginary part of the reconstructed image. This method has been successfully applied to simulated and experimental situations to estimate the focus plane of different kinds of objects (bubbles, calibration target, irregular objects). To complete the proposed reconstruction algorithm, we have proposed a novel method for suppression of the twin image in the reconstructed image that is based on the convolutive blind source separation technique.

172

J. Hattay et al. / Optics and Lasers in Engineering 82 (2016) 159–172

The proposed algorithm has been tested in both simulated and experimental holograms. The tested holograms contain several objects with different physical properties such as the object shapes (bubbles, disks, bars and irregular form), sizes (small and big particles) and recording parameters. The obtained results show that the data sparseness, by means of the AQLS transform, is powerful to find the optimal axial coordinate of the focus plane in the simulated and real cases described above. Then the AQLS tool is followed by the extended complex ICA algorithm. The experimental results indicate the perfect suppression of the unwanted twin image.

References [1] [2] [3] [4]

[5]

[6] [7]

[8] [9]

[10]

[11] [12] [13]

[14] [15]

[16]

[17]

[18]

[19] [20]

Gabor D. A new microscopic principle. Nature 1948;161:777–8. Schnars U, Juptner WPO. Digital holography. Berlin: Springer; 2004. Kreis TM. Handbook of holographic interferometry. Colorado: Wiley; 2005. Dubois Frank, Callens Natacha, Yourassowsky Catherine, Hoyos Mauricio, Kurowski Pascal, Monnom Olivier. Digital holographic microscopy with reduced spatial coherence for three-dimensional particle flow analysis. Appl Opt 2006;45(5):864–71. De Caro L, Giannini C, Pelliccia D, Mocuta C, Metzger TH, Guagliardi A, et al. Inline holography and coherent diffractive imaging with x-ray waveguides. Phys Rev B 2008;77 081408 (R). Korecki P, Materlik G, Korecki J. Complex γ-ray hologram: solution to twin images problem in atomic resolution imaging. Phys Rev Lett 2001;86:1534. Coetmellec S, Pejchang D, Allano D, Grehan G, Lebrun D, Brunel M, et al., Digital in-line holography in a droplet with cavitation air bubbles. J Eur Opt Soc - Rapid Publ. 2014;9. http://dx.doi.org/10.2971/Jeos.2014.14056. Stoykova Elena, Kang Hoonjong, Park Jiyung. Twin-image problem in digital holography—a survey. Chin Opt Lett 2014;12(6):060013 1671-7694 COL. Reza Jafarfard Mohammad, Moon Sucbei, Tayebi Behnam, Young Kim Dug. Dual-wavelength diffraction phase microscopy for simultaneous measurement of refractive index and thickness. Opt Lett 2014;39(10):2908–11. Wu Jian, Feng Lu Ming, Chao Dong Yan, Zheng Ming, Huang Meng, Nan Wu Yi. Zero-order noise suppression with various space-shifting manipulations of reconstructed images in digital holography. Appl Opt 2011;50(34):H56–61. Xu Jiancheng, Xu Qiao, Peng Hansheng. Spatial carrier phase-shifting algorithm based on least-squares iteration. Appl Opt 2008;47(29):5446–53. Zhang Zhihui, Guo Hongwei. Principal-vector-directed fringe-tracking technique. Appl Opt 2014;53(31):7381–93. Denis L, Fournier C, Fournel T, Ducottet C. Numerical suppression of the twinimage in in-line holography of a volume of micro-objects. Meas Sci Technol 2008;19:074004. Picart P, Leval J. General theoretical formulation of image formation in digital Fresnel holography. J Opt Soc Am A Opt Image Sci Vis 2009;26(2):244. Weng J, Zhong J, Hu C. Digital reconstruction based on angular spectrum diffraction with the ridge of wavelet transform in holographic phase-contrast microscopy. Opt Express 2008;16(26):21971–81. Hattay Jamel, Belaid Samir, Lebrun Denis, Naanaa Wady. Digital in-line particle holography: twin-image suppression using sparse blind source separation. J Signal, Image Video Process 2015;9(8), ISSN 1863-1703, Springer. Hattay Jamel, Belaid Samir, Naanaa Wady. Non-negative matrix factorisation for blind source separation in wavelet transform domain. IET Signal Process 2015;9(2):111–9. http://dx.doi.org/10.1049/iet-spr.2013.0409 Print ISSN 17519675, Online ISSN 1751-9683. Hattay Jamel, Belaid Samir, Naanaa Wady. Multiresolution convolutive blind source separation using adaptive lifting scheme. In: IEEE 20th international conference on electronics, circuits, and systems (ICECS), 8–11 December 2013, Abu Dhabi UAE; 2013. p. 273–6. http://dx.doi.org/10.1109/ICECS.2013.6815407. Slimani F, Grhan G, Gouesbet G, Allano D. Near-field Lorenz–Mie theory and its application to microholography. Appl Opt 1984;23(22):4140–8. Onural L. Diffraction from a wavelet point of view. Opt Lett 1993;18:846–8.

[21] Belaid S, Lebrun D, Ozkul C. Application of two dimensional wavelet transform to hologram analysis: visualization of glass fibers in a turbulent flame. Opt Eng 1997;36:1947–51. [22] Coetmellec S, Lebrun D, Ozkul C. Characterization of diffraction patterns directly from in-line holograms with the fractional Fourier transform. Appl Opt 2002;41:312–9. [23] Gillespie John, King RA. The use of self-entropy as a focus measure in digital holography. Pattern Recognit Lett 1989; 9 [North-Holland]. [24] Virtanen Ilkka. Entropy correlation coefficient, a measure of statistical dependence for categorized data. In: Paper presented at the 15th European meeting of statisticians, Palermo, Italy; 1982. [25] Thum C. Entropy of an image with application to focussing. Opt Acta 1984;31:203. [26] Pan Gang, Meng Hui. Digital holography of particle fields: reconstruction by use of complex amplitude. Appl Opt 2003;42(5):827–33. http://dx.doi.org/ 10.1364/AO.42.000827. [27] Liebling Michael, Unser Michael. Autofocus for digital Fresnel holograms by use of a Fresnelet-sparsity criterion. J Opt Soc Am A 2004;21(December (12)):2424–30. [28] Widjaja J, Jutamulia S. Use of wavelet analysis for improving autofocusing capability. Opt Commun 1998;151:12–4. [29] Hattay J, Benazza-Benyahia A, Pesquet J-C. Adaptive lifting schemes using variable size block segmentation. In: Proceedings of the ACIVS Conference, Brussels, Belgium, August 31–September 3, 2004. [30] Rao CR, Nayak TK. Cross entropy, dissimilarity measures, and characterizations of quadratic entropy. IEEE Trans Inf Theory 1985;IT-31(5):589–93. [31] Amel Benazza-Benyahia, Jean-Christophe Pesquet, Jamel Hattay, Hela Masmoudi. Block-based adaptive vector lifting schemes for multichannel image coding. J Image Video Process 2007(1) Eurasip 2007. [32] Antonini M, Barlaud M, Mathieu P, Daubechies I. Image coding using wavelet transform. IEEE Trans Image Process 1992;1(2):205–20. [33] Sharifi K, Garcia L. Estimation of shape parameter for generalized Gaussian distribution in subband decompositions of video. IEEE Trans Circuits Syst Video Technol 1995;5(1):52–6. [34] Do MN, Vetterli M. Wavelet-based texture retrieval using generalized Gaussian density and Kullback–Leibler distance. IEEE Trans Image Process 2002;11 (2):146–58. [35] Gish H, Pierce JN. Asymptotically efficient quantizing. IEEE Trans Inf Theory 1968;14(September):676–83. [36] Hattay J, Benazza-Benyahia A, Pesquet J-C. Adaptive lifting for multicomponent image coding through quadtree partitioning. In: Proceedings of the IEEE international conference on acoustics, speech, signal processing, Philadelphia, USA; March 2005. [37] Bjorck Ake. Numerical methods for least squares problems. Philadelphie, SIAM, ISBN 978-0-89871-360-2, LCCN 96003908; 1996. [38] Vaya C, Rieta JJ, Sanchez C, Moratal D. Performance study of convolutive BSS algorithms applied to the electrocardiogram of atrial fibrillation. In: ICA'06; 2006. p. 495–502. [39] Sweldens W. The lifting scheme: a new philosophy in biorthogonal wavelet constructions. SPIE Wavelet Appl Signal Image Process III 1995;2569:68–79. [40] Hattay Jamel, Belaid Samir, Naanaa Wady. Geometric blind source separation using adaptive lifting scheme. In: 17th IEEE SPA conference SPA'2013, Poland; September 2013. [41] Comon P. Independent component analysis: a new concept. Signal Process 1994;36:287–314. [42] Amari SI, Cichocki A, Yang HH. A new learning algorithm for blind source separation. In: Advances in Neural Information Processing Systems 8, MIT Press, 1996, p. 757–763. [43] Hyvarinen A, Karhunen J, Oja E. Independent component analysis, Wiley, June 2001, p. 504, ISBN: 978-0-471-40540-5. [44] Torkkola K. Blind separation of convolved sources based on information maximization. In: Proceedings of IEEE workshop on neural networks and signal processing (NNSP'96), Kyoto, Japan. [45] Bingham E, Hyvarinen A. A fast fixed-point algorithm for independent component analysis of complex-valued signals. Int J Neural Syst 2000;10(1):1–8. [46] Bell A, Sejnowski T. An information–maximization approach to blind separation and blind deconvolution. Neural Comput 1995;7:1129–59 [1996, 423–32]. [47] Donoho D. Denoising by soft thresholding. IEEE Trans Inf Theory 1995;41(May (3)):613–27.