## Abstract

In the paper the optical diffraction tomographic system for reconstruction of the internal refractive index distribution in optical fiber utilizing grating Mach-Zehnder interferometer configuration is explored. The setup applies afocal imaging. Conventional grating application gives, however, highly aberrated object beam producing incorrect refractive-index reconstructions. The grating inherent aberrations are characterized, its influence on both image projections and refractive index reconstructions is presented. To remove aberrations and enable tomographic reconstruction a novel digital holographic algorithm, correcting optical system imaging, is developed. The algorithm uses plane wave spectrum decomposition of optical field for solving diffraction problem between parallel and tilted planes and enabling correction of imaging system aberrations. The algorithm concept was successfully proved in simulations and the experiment.

© 2009 Optical Society of America

## 1. Introduction

Optical Diffraction Tomography (ODT) [1] is the method for characterization of 3D distribution of material parameters in optical semitransparent elements. These include refractive-index, absorption or birefringence distribution. 3D object information is reconstructed from a set of object scattered complex optical fields obtained from experiment with a rotated object and coherent incident field. For every object angular position projection images are captured by means of interferometry [2], digital holography [3] or photoelasticity [4]. To reconstruct an object internal structure the captured projection should well approximate object integrated phase and amplitude distributions. Various measurement techniques, especially coherent ones share problems of noise and high sensitivity to external factors (vibrations, temperature variations, etc.). In case of ODT the influence of external factors is even higher as the setup has to be stable during entire capturing process (usually over 100 images). Therefore we examine the applicability of grating interferometric setup, which allows to minimize the problems mentioned above.

In the paper we utilize the grating ODT (GODT) system presented in the Fig. 1(a). The setup bases on the classical Mach-Zehnder interferometer with high frequency diffraction gratings playing the roles of a beam splitter and beam recombiner, respectively [5]. The main difference between the classical Mach-Zehnder interferometric system and the proposed configuration is the off-axis propagation of the object beam in respect to the imaging system. It may cause some unexpected degradation of the quality of results and needs investigation. The GODT system utilizes collimated coherent illumination and afocal imaging optics. The system is characterized by compact design and decreased sensitivity to vibration. Detected fringe pattern is not influenced by vibrations present within illumination and imaging modules [6,7]. We are developing new design of massive glass or PMMA based integrated GODT measurement head (Fig. 1(b)) [8] providing fully insensitive instrumentation to vibration. Finally the GODT configuration will by applied in compact, portable and outdoor measurement system, e.g. for optical fiber splices quality assessment. To prove system concept the free-space configuration presented in Fig. 1(a) has been developed. Applied grating configuration, without presence of a measured object, gives achromatic interference fringes [9]. The source with broad spectrum could be used as an illumination to reduce coherent noise. Unfortunately an object presence introduces wavelength dependent optical path difference reducing fringe contrast.

In the presented experimental configuration the object angular orientation relative to the incident plane wave is sequentially changed in the range from 0° to 360°. For every angular position *ϕ* the complex optical field scattered by the object is calculated using the temporal phase shifting method [10]. In ODT two main algorithms are used for tomographic reconstructions: the filtered back propagation [11] and filtered back projection [12]. Both algorithms require small amount of diffraction generated by the characterized structure. The filter back projection algorithm bases on the assumption that a phase image of the object center well approximates object projection [13,14]. The filtered back propagation algorithm accounts for diffraction applying Born or Rytov approximation. Experimentally the best ODT reconstructions are obtained if prior to application of the reconstruction algorithm complex optical field is propagated to the object center numerically or optically. In the system with the applied optical magnification the imaging plane is located inside the object, where the observed diffraction effects in the object beam are minimal.

The exemplary intensity image of an object beam generated by a multimode fiber serving as an object acquired in our setup is presented in Fig. 2. The image was taken with a special care in order to minimize visible diffraction effects. In a well aligned ODT system and with proper choice of immersion liquid refractive index the multimode fiber should not be visible within the object beam image. Additionally we acquired the image of the optical fiber illuminated with 0 order, while ±1 orders were blocked. There were no visible diffraction effects. We conclude that observed strong and nonsymmetrical diffraction in the object beam is an indication of high aberrations. The analysis of aberrations introduced by diffraction gratings illuminated by a non collimated beam has been given elsewhere [15]. It has been demonstrated that in the system where a grating is illuminated with a spherical beam wave front aberrations appear. Presented analysis was focused on converging-diverging beams (point sources). In our case (object dependent optical field) the aberration analysis based on a plane wave decomposition is more suited then the analysis based on sets of the object points.

Section 2 of the paper is devoted to the aberration analysis of the GODT system. The analysis utilizing plane wave harmonics illustrates the origin of aberrations and the nature of imaging errors. In Sec. 3 we introduce a digital holographic algorithm which removes the aberration errors mentioned above. The applicability of the developed holographic algorithm is proven by simulation and experiment.

## 2. Analysis of imaging system errors

In Fig. 3 the imaging of the object in the GODT system is presented. The optical system is focused inside the object center, this is required by ODT. The optical axis of imaging system and optical axis of the object beam are localized within beams referring to different grating diffraction orders (+1^{th} and 0^{th} respectively). Therefore one is unable to conjugate the object and the detector planes optically. The optical system conjugates detector and distorted image planes introducing severe aberrations.

To study the imaging aberrations introduced by the grating we perform a spatial frequency analysis. The imaging system aberrations and grating spatial limitations are neglected. Presented grating introduces nonlinear transform of a single object spatial harmonics *f _{xi}* from object plane into a single image harmonics in the distorted image plane

where Λ is grating period. The object spatial Fourier spectrum modification is presented in Fig. 4. The grating introduces a highly nonlinear frequency shift into the imaging system. It can be expected, therefore, that the image composed of shifted harmonics will be strongly aberrated.

To visualize the effect of errors on the observed image, we have simulated the GODT system presented in Fig. 1. Unfortunately no single modeling tool can accomplish that task. We have used combination of two methods: free space propagation method (FSPM) [16] and wave propagation method (WPM) [17]. The WPM is numerically a very efficient tool to simulate light propagation of an optical beam in an optical model with slow variation of the refractive index distribution. The method is used to propagate the light through the measured specimen (optical fiber) placed in the immersion liquid. In all other system regions free space propagation between parallel and titled planes is applied. The propagation between tilted planes can be simulated by two methods. In the first approach standard propagation algorithms are applied point by point - this approach is numerically not efficient. Second approach applies nonlinear modification of the plane wave spectrum [18,19]. The method uses a spectrum approximation, therefore it can introduce some errors. To present the imaging errors we have applied the point by point approach.

As it was mentioned above, the image created by the grating is strongly aberrated. The question is: where the image plane shall be located. We selected the image position of on-axis object point given by low aperture object rays (Fig. 3). In Fig. 5 the object projection image of an optical multimode fiber is compared with its phase projection. The phase image is distorted and magnified what leads to the most disturbing error. The phase variation error could be disregarded in the measurement setup. Observed variations can have the same order of magnitude as the measurement noise. Therefore after the measurement of an optical structure with tomograph in off-axis configuration, one may come to the conclusion that measured image projection is a correct one. Unfortunately in this case, the tomographic reconstruction will by considerably corrupted (Fig. 6).

The nature of the imaging error can be deduced from Eq. (1). Lets compare Fourier spectra of optical fields in planes X and X’. The Fourier spectrum in plane X’ is a nonlinearly modified copy of Fourier spectrum from plane X. There are three modifying terms in Eq. (1): first term represents linear modification, second nonlinear one and the third - the frequency shift introduced by the grating.

## 3 Digital holographic reconstruction of the projection image

#### 3.1. Algorithm Concept

The projection image given in our system is highly distorted what prevents correct reconstruction of the object refractive index distribution. To recover proper projection image, it is necessary to apply a special digital holographic reconstruction algorithm. The problem of tilted plane imaging was also recognized in digital holography [21,22]. However none of elaborated methods could be applied in our situation. Our system is highly off-axis and complex one, the grating is located on the optical path changing propagation direction. The goal of the developed algorithm is to convert the image observed by the camera to the object projection image. If the detector is conjugated with distorted image plane given by the grating, the reconstruction algorithm should remove distortion from the imaging process. Generally, the detector can be conjugated with any other plane. Our algorithm propagates the beam between tilted planes including change of propagation direction introduced by the grating.

The holographic reconstruction algorithm visualized in Fig. 7 consists of a few steps. At first the captured complex optical field at plane X_{h}’ is propagated to the grating plane. At the grating plane X_{g} the direction of propagation is changed into highly off-axis one. Finally such a field is propagated to an object centre (plane X_{O}). The digital holographic algorithm could be composed of standard free space propagation algorithms and could perform the above mentioned job, unfortunately in an highly inefficient manner. The sampling have to be very dense and the size of processed data will be very large (both on-axis and off-axis fields are processed by the algorithm). Additionally splitting algorithm into unnecessary steps highly enlarges algorithm computational load. Therefore we have elaborated a special holographic reconstruction algorithm.

Let us consider diffraction in a linear, homogenous and isotropic medium where a plane wave spectrum [23] is related to optical field *u*(**x**) by equation

where corresponding vectors are given by **x**=*xê** _{x}*+

*yê*

*,*

_{y}**f**

*=*

_{t}*f*

_{x}ê*+*

_{x}*f*

_{y}ê*. The*

_{y}*ê*

*are the unit vectors and dot denotes the vector scalar product. The spatial frequencies ft are connected with wave vector*

_{x,y,z}**k**=2

*π*(

*ê*

*+*

_{x}f_{x}*ê*

*+*

_{y}f_{y}*ê*

*(*

_{z}*λ*

^{-2}-

**f**·

_{t}**f**)

_{t}^{1/2}). The evanescent waves are not considered within this paper. The above formula allows to decompose optical field into a collection of plane waves characterized by wave vectors and complex amplitudes. Such plane waves can be defined in an arbitrarily tilted coordinate system.

Within this section, we elaborate a digital holographic reconstruction algorithm utilizing the plane wave spectrum decomposition The method allows to develop one special numerical tool consisting of multiple free space propagation steps including diffraction between parallel and tilted planes. Both diffraction algorithms consist in modification of plane wave spectrum. Our digital hologram reconstruction algorithm can be split into two major steps: step 1 - conversion of input optical field given at the holographic image plane X_{h}’ into a plane wave spectrum decomposition defined in the tilted object projection image plane X_{O}’ using propagation algorithm between parallel planes; step 2 - conversion of the plane wave spectrum decomposition defined in the tilted object projection image plane X_{O}’ into optical field in the object projection image plane X_{O} using propagation algorithm between two tilted planes.

#### 3.2. Algorithm step 1 - propagation between parallel planes

The plane wave spectrum at the plane XO’ is related to plane wave spectrum at the plane X_{h}’ by

The formula is simple but numerically inefficient. There are highly on-axis and off-axis fields distributions. In Fig. 8 the spatial shift between input and output fields is visualized. The spatial shift requires both large dimensions of a processed matrix and dense sampling. For an exemplary practical situation: Λ=1 µm, z_{2}=5 mm, the required sampling Δ_{x}≈0.3 µm, the spatial shift 3.14 mm, we need *N _{x}*=13000 sampling points.

To obtain a robust algorithm we assume that our optical field is a paraxial one and the input optical frequency spectrum is band limited to frequency **f**
_{m}=[-*f _{m}*,

*f*]. In diffraction tomography slow variations of refractive index distribution are characterized only. Additionally, during the tomographic reconstruction the low pass filtering operations are applied. For computational simplicity it is convenient to design algorithm (Eq. (3)) working with the shifted frequency

_{m}**f**

_{ts}=(

*f*=

_{xs}*f*-Λ

_{x}^{-1},

*f*):

_{y}The object spatial spectrum is now centered in frequency domain. However the third factor of Eq. (4) has highly inclined phase at the object frequencies, necessary for an object spatial shift. After removal of the spatial frequency tilt we get the final formula for PWS at the plane X_{O’}:

where the phase accommodation kernel is $H\left({\mathbf{f}}_{\mathrm{ts}}\right)=\mathrm{exp}\left\{ik{z}_{1}{(1-{\lambda}^{2}{\mathbf{f}}_{\mathrm{ts}}\xb7{\mathbf{f}}_{\mathrm{ts}})}^{1\u20442}\right\}\mathrm{exp}\left\{ik{z}_{2}{{(1-{\lambda}^{2}\left({f}_{\mathrm{xs}}+{\mathbf{\Lambda}}^{-1}\right)}^{2}+\lambda}^{2}{{f}_{y}^{2})}^{1\u20442}\right\}\mathrm{exp}\left\{\genfrac{}{}{0.1ex}{}{2\pi i\lambda {z}_{2}{f}_{\mathrm{xs}}}{\mathbf{\Lambda}{\left(1-{\lambda}^{2}{\mathbf{\Lambda}}^{-2}\right)}^{1\u20442}}\right\}\xb7.$ With this formula spatial shift and high frequency phase modulation are removed. Therefore computed optical field is centered in both spatial and frequency domains. To get required number of samples we have performed the spatial frequency position analysis [24]:

The algorithm configuration requires that spatial frequency position can not exceed dimensions of computational window, i.e., *x _{lf}*≤

*N*Δ

_{x}*/2. This gives required number of samples:*

_{x}Considering parameters of practical tomographic imaging system: sample size Δ=1 µm, wavelength *λ*=0.532 µm, Λ=1 µm, *z _{1}*=6.459 mm,

*z*=-5 mm and

_{2}*f*=326.4 mm

_{m}^{-1}(such a bandwidth corresponds 10° aperture angle), we get necessary computational window dimension N

_{x}>1594. By using computation window of that size, we ensure that the aliasing in the phase accommodation kernel occurs at higher frequency than fm and the entire optical field bandwidth is processed correctly.

#### 3.3. Algorithm step 2 - propagation between tilted planes

The plane wave vector can be defined in an arbitrarily rotated coordinate system. The arbitrary wave vector **k**′=*k*′_{x}ê* _{x}*+

*k*′

*+*

_{y}ê*k*′

_{z}ê*defined in plane (x’, y, z’) is transferred into*

_{z}**k**=

**Tk**′ defined in (x, y, z), where

is the matrix of coordinate rotation according to axis y by angle α. By application of matrix **T**, the full collection of the plane waves defining an arbitrary optical field can be characterized in a new rotated coordinate system. In this way we can perform mapping of plane wave spectra between planes X_{O’} and X_{O} by

This equation expresses an optical field propagation between two tilted planes as nonlinear mapping of frequency spectra. The approximation methods have to be applied in reality. Nonlinearly distributed nodes of frequency spectrum have to be mapped into uniformly spaced frequency grid. It is required for application of the Discrete Fourier Transform. The best results were obtained with cubic interpolation [25].

Nonlinear spectrum mapping with Eq. (9) causes a shift of *U ^{X}*′(

**f**

*) by the carrier Λ*

_{t}^{-1}. Such a frequency shift unnecessarily enlarges computer power consumption of algorithm and requires dense sampling. For this reason we propose application of a modified mapping with shifted frequency coordinates:

Nonlinear mapping of plane wave spectrum with this equation works within low frequency bandwidth only and sampling below Λ^{-1} can be applied.

Finally, we can relate the complex optical fields in plane X_{O’} and projection image plane X_{O} by formula:

where $J\left({\mathbf{f}}_{t}\right)=\mathrm{cos}\alpha +\genfrac{}{}{0.1ex}{}{{f}_{x}}{{({\lambda}^{-2}-{\mathbf{f}}_{t}\xb7{\mathbf{f}}_{t})}^{1\u20442}}\mathrm{sin}\alpha $ and the plane wave spectrum is defined by Eq. (10). For paraxial field the Jakobian is negligible.

In Fig. 9 the phase projection image given by the outlined above image correction algorithm is compared with a phase projection of the object. One can notice that both distributions are very similar. The observed discrepancy is caused by diffraction of plane wave illuminating an optical fiber (step-wise structure).

#### 3.4. Experimental verification of the algorithm

To verify the analysis described above, we have implemented the digital holographic algorithm based on Eq. (11) in experimental set-up presented in Fig. 1. The exemplary input fringe pattern is shown in Fig. 10(a). For every angular object orientation five interferograms with i(π/2) phase shift were taken and the phase of the object beam was reconstructed using five frame phase shifting algorithm (Fig. 10(b)). The aberrations can be observed in both figures. Next, the complex optical field was processed using the developed holographic algorithm and the reconstruction in object plane was obtained. The phase was unwraped and phase tilt was removed giving the required phase projection image for tomographic reconstruction. The final phase projection image is presented in Fig.10(c). It is seen that the aberrations are significantly minimized. Next, the ODT data capturing was continued and for equally spaced 90 angular object positions the phase projection images were obtained. The result of the tomographic reconstruction is presented in the Fig. 10(d).

## 4. Conclusions

In the paper the optical diffraction tomographic system utilizing two gratings in the Mach-Zehnder interferometer configuration is explored. Within the system an object beam changes direction of propagation and its image given by optical system is highly aberrated. The presented plane wave spectrum aberration analysis reveals that the imaged object optical field is both linearly and nonlinearly disturbed. The tomographic reconstruction performed on such a beam will be spatially extended, reconstructed coefficients will be reduced and nonlinearly modified. In order to remove system aberrations and to enable accurate tomographic reconstruction, the spatial digital holographic algorithm correcting optical imaging system has been developed. The algorithm uses plane wave spectrum decomposition of optical field solving diffraction between parallel and tilted planes. The elaborated numerical tool applies the system dependent plane wave spectrum phase accommodation and changes nonlinear spacing between frequency nodes into a linear one. The special attention is devoted to the reduction of algorithm computer power requirements. Both simulations and experimental results are presented. Aberration and its correction are visualized showing both corrupted and corrected projections and reconstructions.

## Acknowledgements

We would like to acknowledge the support of the Ministry of Science and Higher Education within the projects N505 008 31/1374, N505 4285 33 and the statutory founds.

## References and links

**1. **M. Born and E. Wolf, *Principles of Optics*7th (expanded) edition, (Cambridge University Press, 1999).

**2. **W. L. Howes and D. R. Buchele, “Optical Interferometry of Inhomogeneous Gasses,” J. Opt. Soc. Am. **56**, 1517–1528 (1966). [CrossRef]

**3. **L. P. Yaroslavskii and N. Merzlyakov, *Methods of digital holography* (Consultants Bureau, New York, 1980).

**4. **S. K. Mangal and K. Ramesh, “Determination of characteristic parameters in integrated photoelasticity by phase-shifting technique,” Opt. Laser Eng. **31**, 263–278 (1999). [CrossRef]

**5. **W. Górski and M. Kujawińska, “Three-dimensional reconstruction of refractive index inhomogeneities in optical phase elements,” Opt. Lasers Eng. **38**, 373–385 (2002). [CrossRef]

**6. **E.N. Leith and L. Shentu, “Tomographic reconstruction of objects by grating interferometer,” Appl. Opt. **25**, 907–913 (1986). [CrossRef]

**7. **L. Sałbut and M. Kujawińska, “The optical measurement station for complex testing of microelements,” Opt. Laser Eng. **36**, 225–240 (2001). [CrossRef]

**8. **R. Krajewski, M. Kujawińska, B. Volckaerts, and H. Thienpont “Low-cost microinterferomtric tomograpy system for 3D refraction index distribution Measurements in optical fiber splices, Proc. SPIE **5855**, 17th OFS conference Brugge, 347–351 (2005). [CrossRef]

**9. **R. Czarnek, “High sensitivity moiré interferometry with compact achromatic interferometer,” Opt. Laser Eng. **13**, 93–101 (1990).

**10. **I. Yamaguchi and T. Zhang, “Phase-shifting digital holography,” Opt. Lett. , **22**, 1268–1270 (1997). [CrossRef]

**11. **A.J. Devaney, “A filtered backpropagation algorithm for diffraction tomography,” Ultrasonic imaging **4**, 336–350 (1982). [CrossRef]

**12. **J. Hseigh, *Computed Tomograph*, (SPIE Press, Washington, 2003).

**13. **T.C. Wedberg, J.J. Stamnes, and W. Singer, “Comparison of filtered backpropagation and backprojection algorithms for quantitative tomography,” Appl. Opt. **34**, 6575–6581 (1995). [CrossRef]

**14. **T. Kozacki, M. Kujawinska, and P. Kniazewski, “Investigating the limitation of optical scalar field tomography,” Opto-Electron. Rev. **15**, 102–109 (2007). [CrossRef]

**15. **P. Ferraro, S. De Nicola, A. Finizio, and G. Pierattini, “Reflective grating interferometer in a noncollimated configuration,” Appl. Opt. **39**, 2116–2121 (2000). [CrossRef]

**16. **F. Shen and A. Wang, “Fast-Fourier transform based numerical integration method for the Rayleigh-Sommerfeld diffraction formula,” Appl. Opt. **45**, 1102–1110 (2006). [CrossRef]

**17. **K.-H. Brenner and W. Singer, “Light propagation through microlenses: a new simulation method,” Appl. Opt. **32**, 4984–4988 (1993). [CrossRef]

**18. **N. Delen and B. Hooker, “Free-space beam propagation between arbitrarily oriented planes based on full diffraction theory: a fast Fourier transform approach,” J. Opt. Soc. Am. A **15**, 857–867 (1998). [CrossRef]

**19. **K. Matsushima, H. Schimmel, and F. Wyrowski, “Fast calculation method for optical diffraction on tilted planes by use of the angular spectrum of plane waves,” J. Opt. Soc. Am. A **20**, 1755–1762 (2003). [CrossRef]

**20. **C. J. Dasch, “One-dimensional tomography: a comparison of Abel, onion-peeling, and filtered backprojection methods,” Appl. Opt. **31**, 1146–1152 (1992). [CrossRef]

**21. **D. Lebrun, A. Benkouider, S. Coëtmellec, and M. Malek, “Particle field digital holographic reconstruction in arbitrary tilted planes,” Opt. Express **11**, 224–229 (2003). [CrossRef]

**22. **L. Yu, Y. An, and L. Cai, “Numerical reconstruction of digital holograms with variable viewing angles,” Opt. Express **10**, 1250–1257 (2002). [PubMed]

**23. **J. J. Stamnes, *Waves in Focal Regions* (Hilger, Bristol1986).

**24. **T. Kozacki, “Numerical errors of diffraction computing using plane wave spectrum decomposition,” Opt. Commun. **281**, 4219–4223 (2008). [CrossRef]

**25. **T.M. Lehmann, C. Gonner, and K. Spitzer, “Survey: interpolation methods in medical image processing,” IEEE Trans. Med. Imaging **18**, 1049–1075 (1999). [CrossRef]