Three-dimensional optical imaging in layered media

The present paper deals with the reconstruction of threedimensional objects from the scattered far-field. The config uration under study is typically the one used in the Optical Diffraction To m graphy (ODT), in which the sample is illuminated with various angle s of incidence and the scattered field is measured for each illumination. Th e retrieval of the sample from the scattered field is accomplished numeri cally by solving the inverse scattering problem. We present herein a fast method for solving the inverse scattering problem based on the Coup led Dipole Method (CDM) and applied it for complex background configura tion such as buried objects in a layered medium. Numerical experi ments are reported and robustness against the presence of noise in the data is analyzed. © 2006 Optical Society of America OCIS codes: (180.6900) Three-dimensional microscopy; (100.6640) Super resolution; (100.6890) Three-dimensional image processing References and links 1. A. Chomik, A. Dieterlen, C. Xu, O. Haeberlé, J. J. Meyer and S . Jacquey, “Quantification in optical sectioning microscopy: a comparison of some deconvolution algorithms in vi ew of 3D image segmentation," J. Opt. 28, 225 (1997). 2. J. O. Tegenfeldt, O. Bakajin, C.-F Chou, S. S. Chan, R. Aust in, W. Fann, L. Liou, E. Chan, T. Duke, E. C. Cox, “Near-field Scanner for Moving Molecules," Phys. Rev. Lett. 86, 1378 (2001). 3. L. A. Ghebern, J. Hwang, M. Edidin, “Design and optimizatio n of a near field scanning optical microscope for imaging biological samples in liquid," Appl. Opt. 37, 3574 (1998). 4. A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging , Society of Industrial and Applied Mathematics, (2001). 5. V. Lauer, “New approach to optical diffraction tomography yielding a vector equation of diffraction tomography and a novel tomographic microscope," J. Microsc. 205, 165 (2002). 6. P. S. Carney and J. C. Schotland, “Three-Dimensional Total Internal Reflection Microscopy," Opt. Lett. 26, 1072 (2001). 7. P. S. Carney and J. C. Schotland, “Theory of total-interna l-reflection tomography," J. Opt. Soc. Am. A 20, 542 (2003). 8. K. Belkebir, P. C. Chaumet, A. Sentenac, “Influence of multip le scattering on three-dimensional imaging with optical diffraction tomography," J. Opt. Soc. Am. A. 23, 586 (2006). 9. K. Belkebir, P. C. Chaumet, A. Sentenac, “Superresolution in total-internal reflection tomography," J. Opt. Soc. Am. A. 22, 1889 (2005). 10. P. C. Chaumet, A. Sentenac, and A. Rahmani, “Coupled dipole method for scatterers with large permittivity," Phys. Rev. E70, 036606 (2004). 11. A. Rahmani, P. C. Chaumet, and F. de Fornel, “Environment-in duced modification of spontaneous emission: Single-molecule near-field probe," Phys. Rev A 63, 023819 (2001). 12. P. C. Chaumet, K. Belkebir, A. Sentenac, “Three-dimension al sub-wavelength optical imaging using the coupled dipole Method," Phys. Rev. B, 69, 245405 (2004). #10347 $15.00 USD Received 17 January 2006; revised 30 March 2006; accepted 1 April 2006 (C) 2006 OSA 17 April 2006 / Vol. 14, No. 8 / OPTICS EXPRESS 3415


Introduction
The classical optical microscope have a resolution above the well known Rayleigh criterion and give a two dimensional image of three dimensional objects .Using the optical sectioning technique a three dimensional image can be obtained and the resolution can be improved with a deconvolution technique but this is time consuming.[1] In the last decades, intensive developments of methods have been carried out to image biological samples using electromagnetic probes [2,3].The resolution is in some cases below the Rayleigh criterion, i.e., better than a half of the wavelength of the exciting field.Microscopes operating in the near field, for example the Scanning Near-Field Optical Microscope (SNOM) [3], have the disadvantage of approaching the probe close to the sample.Thus, the interaction between the probe and the sample -which is not easy to model, since it depends on the shape of the probe and on the constitutive materials of both sample and the probe-blurs the image of the sample.Furthermore, the probe scanning on the top of the surface, it is not obvious to extract information related to objects buried in a substrate.In addition, for the case of a specimen deposited on a substrate, moving the probe along the surface may damage the sample which aimed to be imaged.
In the present paper, we consider an optical imaging system based on the Optical Diffraction Tomography (ODT) technique, which circumvents aforementioned disadvantages, with however an inferior resolution to the one achieved with near field microscopes.The basic idea underlying the ODT technique is firstly lighting the sample with various illuminations and secondly retrieving the object under test from the scattered field supposed known (modulus and phase).This second part requires a numerical procedure for solving the inverse scattering problem.Under the Born approximation the inverse scattering problem is linear and an image of the object under test may be performed through a simple inverse Fourier transform [4,5] or thanks to the singular value decomposition [6,7].
We have recently developed an algorithm that permits to retrieve accurately the threedimensional relative permittivity of scattering objects present in homogeneous background [8] or above a dielectric substrate [9] with less than 100 points of observations per illumination.The main bottleneck of this method is its greed in terms of time computation.Typically, several hours are needed with a rather simple configuration, e. g., objects of characteristic dimension of λ /4 present in a homogeneous background medium, and inversion performed with an investigating domain of volume size 8λ 3 (λ being the wavelength of the incident field) [8].In the biological application when studying many samples, this amount of computation time is overthrowing.In the present paper, we present a method which allows to locate the unknown objects in complex background environment such as layered media, and to differentiate the constitutive materials of the unknown objects in terms of absorbing or transparent materials.However, this method is not quantitative since the refractive indexes of the objects under test are not accurately retrieved.On the other hand the inversion is performed within only few minutes for a configuration involving layered media while the full solution [8] is obtained within several hours when the background is homogeneous.

Forward scattering problem
We use the coupled dipole method to compute the scattering of light by arbitrary objects.As this method has been presented in previous article, we only recall the main steps [10].The objects under study are discretized into N dipolar subunits, and the field at each subunit satisfies the following self consistent equation: where E inc (r i ) is the incident field, and α(r j ) the polarizability of the jth subunit which meet the Claussius-Mossotti relation: ε(r j ) is the relative permittivity of the subunit j, and d the size of the subunit.S is a tensor which correspond to the linear response of a dipole in the system of reference, i.e., homogeneous space [10], a substrate, or a multilayered system [11].Notice that when i = j in Eq. ( 1) the contact term is take into account through the Clausius-Mossotti relation [10].This is the weak form of the CDM which presents enough accurate for our aims.The dipole moment of the subunit i is written as p(r i ) = α(r i )E(r i ), hence Eq. ( 1) can be written under this symbolic form: where A is a matrix (3N × 3N) which contains all the tensors S. E, E inc and p are vector (3N) which contain all the local field, the incident field, and the dipole moment, respectively.The field scattered by the objects at an arbitrary position r reads as If we have a set of M points of observation, one can use this symbolic form: where B is a matrix (3M ×3N) and f a vector (3M) which contains all the diffracted field.Notice that the matrices A and B do not depend of the incident field and of the nature of the object.

Inversion algorithm
The object is assumed to be confined in a bounded box Ω (test domain or an investigating domain) and illuminated successively by For each excitation l, the scattered field f l is measured on a surface Γ at M points.The inverse scattering problem is now stated as determining properties of unknown objects present in the investigating domain Ω from f l .The approach used for solving this inverse scattering problem is based on a previous work of authors [12] where a simple configuration (homogeneous background medium) is involved.In the present paper we extend the method in order to handle layered media configurations.The basic ideas underlying the inversion algorithm remain the same as in Ref. [12].In this method, a sequence p l,n is built up according to the following recursive relation: where d l,n is an updating direction, β l,n a scalar number determined at each iteration step by minimizing the cost functional F n that represents the discrepancy between the data (measurements) and the scattered field corresponding to the best available estimate of the object p l,n .The cost functional F n is defined thus: Q Γ is the deduced norm from the inner product of two vectors < R, Q > Γ defined on Γ.This inner product reads as < R, Q > Γ = ∑ r k ∈Γ R ⋆ (r k ).Q(r k ) where R ⋆ denotes the complex conjugate of R. Now using p l,n and Eq. ( 6), Eq. ( 7) leads to a polynomial expression of the weighting coefficients β l,n .Then the cost functional is minimized with respect to L scalar coefficients β l,n .As updating direction d l,n the authors took the Polak-Ribière conjugate gradient direction where < ., .> Ω is the same inner product as defined previously but acting on vectors defined on Ω.The vector function g l;p is the gradient of the cost functional F with respect to p l evaluated for the (n − 1) th quantities.This gradient reads as: is the transpose complex conjugate matrix of the matrix B. Once the sources p l are reconstructed, one can determine the fields E l inside Ω using the Eq. ( 3).The polarizability α at the position r j is then given by where * denotes the complex conjugate.The permittivity ε distribution is determined easily using Eq. ( 2).

Numerical results
In this section are presented numerical experiments to illustrate the performance of the imaging method described previously.The section is subdivided in three parts.The first one is devoted to the simple case involving a homogeneous background.The second part deals with objects embedded or deposited on a substrate.Finally, the last part treats the complicated case of targets buried in a layered medium.In all cases, results of reconstruction from corrupted data are presented.Thus, the robustness of the inversion scheme against the presence of noise in data is emphasized.The synthetic data used for inversion were obtained thanks to a forward solver with a mesh size of λ /20 while inversions were performed with a different mesh of size λ /10.All presented results were obtained without any post-treatment.

Simple configuration: case of homogeneous background medium
We start with the study of two objects present in a homogeneous background medium .The two objects are cubes of side a = λ /4 and separated by a distance of c = λ /3 (λ being the wavelength of vacuum).The relative permittivity of the cube located at x ≈ −0.3λ is ε = 2.25 while the relative permittivity of the other one located at x ≈ 0.3λ is ε = 2.25 + i0.5.The illumination of the samples is as described in Fig. 1, i.e., 16 plane waves in the two perpendicular planes (x, z) and (y, z).The electrical field remains in the incidence plane (Fig. 1).Let us denote by k inc the wavevector of the incident field and k d the wavevector of the diffracted field.The investigated domain Ω is a large cube of volume 2λ × 2λ × 2λ .Figures 2(a-d) present results of the inversion described in the previous section.The chosen representation of this configuration is used in the entire paper.For each set of four figures, first row corresponds to the reconstructed real part of the relative permittivity, in the plane (x, z) for the left image and in the plane (x, y) for the right image.The second row is as for the first row but for imaginary part of the relative permittivity instead of the real part.The full line curves represent the boundaries of the actual objects.
The convergence was achieved after 20 iteration steps, we did not observe any marked changes when pursuing the iterative process.Henceforth, all numerical experiments reported in the present paper correspond to the 20th iterate solution.The needed CPU time to complete the inversion was 5 mn using a standard commercial Personal Computer with an internal clock frequency of 3 GHz.Notice that p is obtained only in 80 seconds by minimizing the cost functional F defined in Eq. (7).Therefore, the main computation time is spent to obtained the internal field via Eq.(3).Hence the main advantage of this method is that Eq. ( 3) is used only once.In the method presented in Refs.[8,9] the linear system represented by Eq. ( 3) is solved at each iteration step for each angle of incidence, this explains its amount of the computation time about 10 hours.Figures 2(a) and 2(b) show that the resolution for the real part of the permittivity is better along the x axis than along the z axis.Using Ewald's sphere, the projection of k d − k inc along the z axis yields to an interval of [−k inc ; k inc ] while the projection along the x axis leads to an interval twice larger.This explains the better resolution along x axis.However, the reconstructed objects are not well separated along the x axis.Regarding the reconstruction of the imaginary part of relative permittivity, Figs.2(c) and 2(d), the single absorbing object is perfectly located (the left cube), yet it is shifted downwards.This may be due to the illumination done only in the direction of the z axis positive.Using a symmetric illumination in both directions (z positive and negative), the reconstructed objects would be perfectly centered on the actual objects.Note that a slight absorbing part appears on the right-hand cube, due to a weak coupling between both cubes.Now, we investigate the robustness of the algorithm of reconstruction against a presence of noise in data.In view of approaching the experimental conditions, we corrupt the scattered far-field data, f l=1,•••,L , by an additive uncorrelated noise on each component of the electric field at each position of observation, v stands for the components along x, y, or φ is a random number taken for each component of the positions of observation and incident angles with uniform probability density in [0, 2π]; u is a real number smaller than unity that monitors the noise level.Figures 2(e-h) show the reconstruction in the presence of noise (u = 10%) .It appears evident that the obtained results are similar to those shown in Figs.2(a-d).
We examine now the ability of the inverse algorithm to retrieve a complex shaped object.Therefore, we consider an inhomogeneous "U-shaped" object that is constituted by a dielectric bar and two absorbing cubes.The bar is of width λ and of section (λ /4 × λ /4).The relative permittivity of this bar is ε = 2.25.The two cubes are of side λ /4 and of relative permittivity ε = 2.25 + 0.5i.The illumination as well as the observation is unchanged.Figure 3 presents results of the inversion.It is clearly shown that the "U-shaped" object is accurately retrieved even from corrupted data with a value of u as high as 30%.Thus, the reconstruction method presented here is very robust against a presence of the noise in data.In addition, this fast method can provide us with a 3-D cartography showing objects that are absorbing or not.In the following sections, we will attempt to show that this method has the major advantage of being applicable to extremely complex situations, which would be tedious to carry out with a method as described in Ref. [9].

Objects above a dielectric substrate
We consider the same objects as in Fig. 2 (a = λ /4 and separated by c = λ /3) of the same relative permittivity but now deposited on a flat interface separating the superstrate (air, ε = 1) and the substrate (glass, ε = 2.25).The total reflection angle is in this case θ c = 41.8 • .The investigating domain is of volume (2λ × 2λ × λ ) and located above the interface.
The computation time (5 mn 20 s) is similar to the one may obtained in the case of the homogeneous background medium (5 mn).
Figures 4(a-d) show the reconstruction obtained using only propagative illumination, i.e., |θ inc | < θ c .In this case, Figs.4(a) and 4(b) are similar to Figs. 2(a) and 2(b) in terms of resolution.This is due to the fact that Ewald's spheres are identical in both configurations.However, Figs. 4(c) and 4(d) show an accurate localization of the imaginary part of permittivity.We interpret this superior resolution as being the result of a coupling effect between objects and the substrate, with multiple scattering improving the resolution [9].Figures 4(e-h) show the result of reconstruction using only evanescent illumination, i.e., |θ inc | > θ c .Considering the real part of permittivity, both cubes are now perfectly resolved as shown in Figs.4(e) and 4(f).However, the objects do not lie on the surface.The result for the imaginary part is less accurate, because the absorbing cube seems to float above the surface; and, a significant portion of the imaginary part appears at the object on the right.The use of evanescent illumination thus yields to good resolution, due to the high spatial frequencies provided by the incident wave (enlarged Ewald's sphere), yet the object is not well located.
We now use both propagative and evanescent illumination, i.e., −80 • < θ inc < 80 • .Further- more, noise has been added to the diffracted field (u = 10%) in order to mimic experimental conditions.Figure 5(a-d) clearly show that adding propagative and evanescent waves improves the result: both the real and the imaginary parts show excellent adequation between the actual profile and the reconstructed one.The two cubes are well located and perfectly resolved.Furthermore, Figs.5(c) and 5(d) show only one truly absorbing cube, which is not the case in Figs 5(g) and 5(h) where only evanescent waves are used nor in Figs 4(g) and 4(h) where noiseless data are used.The effect of noise in the reconstruction, when using evanescent waves, Figs.5(e-h) (in particular, artifacts appear at the top of the investigating box in Fig. 5(e)), is much higher than the one observed when using both propagative and evanescent waves Fig. 5(ad).In fact, an illumination with evanescent wave contains high spatial frequencies, which are known to be sensitive to noise.Illuminating targets with both propagative and evanescent waves leads to combine the robustness and the accuracy of the reconstruction.

Objects buried in the dielectric substrate
We consider herein a configuration of objects buried in a glass-substrate (ε = 2.25).In Figs. 6  and 7 the objects are the same as in Section 3.2.1 : cubes with a = λ /4 separated by c = λ /3 with a relative permittivity ε l = 2.25 + 0.5i for the cube located at (−0.3λ , 0, −0.65λ ) and ε r = 2.25 for the cube located at (0.3λ , 0, −0.65λ ).Illumination remains at −80 • < θ inc < 80 • and the dimension of the investigating domain is 2λ × 2λ × 1.5λ and includes the air-glass interface.This configuration would be more difficult to study with an optical microscope.It is typically the one that would be used for mines detection if the scattered fields were measured closed to the interface.This is not a limitation for our method, it suffices to replace the far-field tensor into the near-field one.
Figures 6(a-d) show the reconstruction obtained without noise when the observation points are located only above the substrate as in all the previous figures.The objects are clearly well located below the interface and the cubes are separated.However, the objects are not resolved as clearly as in Figs.5(a-d).This may be due to interactions between the objects and substrate's surface.The imaginary part of the relative permittivity is extremely well located on the left-hand cube, although it appears above the surface like an absorbing object.To improve the quality of reconstruction, observations were carried out below and above the surface.The illumination remains unchanged.We observe in Figs.7(a-d) a good localization along the z axis for the real part of permittivity, and a good lateral separation.In addition, the imaginary part is now perfectly located on the actual cube.We explain this resolution by looking at the Ewald sphere, which is slightly enlarged with respect to the x-axis and much more enlarged in the z-axis direction.In the case of the scattered field corrupted with noise (u = 10%), Figs.6(e-h) and Figs.7(e-h), show that the reconstruction is almost not altered by the presence of noise in the data.This is particularly true when observation points are located above and below the surface (Figs.7(e-h)).
In all previous examples we have restricted our study to objects situated in a single (x, z) plane.One may wonder what would be the reconstruction if the objects are placed somewhere else.This is investigated in Fig. 8 where the observation points are located above and below the surface, the illumination being the same as previously.Figure 8 presents the result of the reconstruction of four targets located in different (x, z) planes.All objects are well retrieved.In fact, objects can be distributed anywhere, we have only chosen objects in particular plane for sake of simplicity.In addition, the computational time is identical for all cases presented in Figs.6-8.Indeed, the computational time is completely independent of the shape and the distribution of the objects inside the investigating domain.It depends only of the size of the investigating domain.

Complex configuration: Case of layered medium
The last configuration studied in this paper concerns a relatively complex configuration.This configuration involves a multilayered system with objects of different permittivities.The chosen geometry is depicted in Fig. 9(a).The illumination is the same as in Section 3.2.2,−80 • < θ inc < 80 • , and the observation points are located only above the surface.
Note that the computational time in this rather complex case is only about 16 mn, which is mainly due to the construction of the multi-layer tensor, i.e. matrix A. Once the matrix A is built the needed computational time for solving the inverse scattering problem remains almost unchanged.Saving the matrix A would be preferred for a repeated imaging objects present in an investigated domain of fixed size.
Figures 9(b-e) show the reconstruction obtained with the geometry described in Fig. 9(a).Note that each layer has its own color scale.In the case of noiseless data, the map of the real part of the relative permittivity, Fig. 9(b), shows that the objects are perfectly retrieved, except the one slightly above the substrate.Figure 9(c), represents the imaginary part of the relative permittivity.It is shown that a single absorbing cube is present in the layered medium.The case where the scattered field is corrupted with noise (u = 10%) is presented in Figs.9(d) and 9(e).The main effect of the noise is to perturb the map of the imaginary part of the relative permittivity.We noticed that when two cubes are located in different layers but one on top of the other, coupling occurs, thus hindering reconstruction.This type of coupling, between different objects present in different layers, deserves to be investigated more in depth.

Conclusion
We have simulated an experiment of optical diffraction tomography.The method that we proposed is a full vectorial inversion scheme.It permits to localize the objects and to discriminate absorbing objects and transparent ones.The objects can be in homogeneous space or put upon a flat substrate or buried in a substrate.Is is also possible to handle a more complex configuration with many objects in a layered medium.The main advantage of our method is the rapidity.The requested computational time is only few minutes.This low computational time can be useful in the biological applications, if one wants to localize the objects in stratified media.Notice that if one knows the value of the relative permittivity of the objects, then a post treatment can be done which will enhance the resolution.

Fig. 1 .
Fig.1.Sketch of the illumination and detection configuration.The observation points are in the far field zone regularly placed on a half sphere.The illumination corresponds to 16 plane wave propagating towards to the positive values of z in both planes (x, z) and (y, z) such that the electric field is in the incidence plane.The angle between the incident wave vector and the z axis ranges over −80 • to 80 • .The cube represents the objects under test and the background can be either homogeneous or layered medium.

Fig. 2 .
Fig. 2. Map of the reconstructed relative permittivity using our inversion scheme.The size of the investigating domain is 2λ × 2λ × 2λ .Objects under test are cubes (boundaries of these cubes are plotted in black) of side a = λ /4 and separated by c = λ /3.The actual relative permittivity of the cube is ε l = 2.25 + 0.5i and ε r = 2.25 for the left and right cube, respectively.The four left figures are obtained from noiseless data, while the right figures are obtained from a corrupted scattered field with u = 10%.a) and e) shows the real part of the relative permittivity in the plane y = 0. b) and f) shows the real part of the permittivity in the plane z = 0. c) and g) shows the imaginary part of the permittivity in the plane y = 0. d) and h) shows the imaginary part of the relative permittivity in the plane z = 0.

Fig. 3 .
Fig. 3. Same as Fig. 2 but with an inhomogeneous U-shaped object.This object is constituted by a dielectric bar (ε = 2.25) and two absorbing cubes located at the extremities (ε = 2.25 + 0.5i).Maps a), c), e) and g) are plotted in the plane y = 0; b) and f) in the plane z = 0; and d) and h) in the plane z = −0.4λ .The dashed lines in the figures a), c), e) and g) represent, the cross sections maps of b), f), d) and h) in the (x, y) plane.The figures on the right are organized as the figures on the left but with a scattered field corrupted with noise (u = 30%).

Fig. 4 .
Fig.4.Map of the reconstructed relative permittivity using the inversion algorithm.The size of the investigating domain is 2λ ×2λ ×λ .We have a = λ /4, c = λ /3, and ε l = 2.25+0.5i,ε r = 2.25, and ε s = 2.25.The four left figures are obtained with only propagative waves, while the figures on the right are obtained with solely evanescent waves.a) and e) show the real part of the permittivity in the plane y = 0. b) and f) show the real part of the permittivity in the plane z ≈ λ /7 (dashed line).c) and g) show the imaginary part of the permittivity in the plane y = 0. d) and h) show the imaginary part of the permittivity in the plane z = λ /7 (dashed line).

Fig. 5 .
Fig. 5. Map of the relative permittivity given by our inversion scheme when the scattered field is corrupted with noise (u = 10%).The size of the investigating domain is 2λ × 2λ × λ .We have a = λ /4, c = λ /3, and ε l = 2.25 + 0.5i, ε r = 2.25, and ε s = 2.25.The four figures on the right are obtained with only evanescent waves, while the figures on the left are obtained with both evanescent and propagative waves.a) and e) show the real part of the permittivity in the plane y = 0. b) and f) show the real part of the permittivity in the plane z = λ /7 (dashed line).c) and g) show the imaginary part of the permittivity in the plane y = 0. d) and h) show the imaginary part of the permittivity in the plane z = λ /7.

Fig. 6 .Fig. 7 .Fig. 8 .
Fig.6.Map of the relative permittivity when the objects are embedded in the substrate.The size of the investigating domain is 2λ × 2λ × 1.5λ .We have a = λ /4, c = λ /3, and ε l = 1.5 + 0.5i, ε r = 1.5, and ε s = 2.25.Left figures (a-d) are obtained without noise and figures on the right (e-h) are obtained with a scattered field corrupted with noise (u = 10%).a) and e) show the real part of the permittivity in the plane y = 0. b) and f) show the real part of the permittivity in the (x, y) plane located at z = −0.65λ(dashed line).c) and g) show the imaginary part of the permittivity in the plane y = 0. d) and h) show the imaginary part of the permittivity in the plane (x, y) located at z = −0.65λ(dashed line).

Fig. 9 .
Fig. 9. a) Sketch of the studied configuration.The dimension of the investigating domain is 2λ × 2λ × 2.2λ and a = λ /4.b) Real part of the relative permittivity in the (x, z) plane.Notice that each layer has its own color scale.c) Imaginary part of the relative permittivity in the (x, z) plane.d) and e) same as b) and c), respectively, but with noisy data (u = 10%).