## Abstract

In this paper, we make a proposal to obtain the Hilbert-transform for each entry of the projection data leaving the slice of a thin phase object. These modified projections are stacked in such a way that they form a modified sinogram called Hilbert-sinogram. We prove that the inverse Radon-transform of this sinogram is the directional Hilbert-transform of the slice function, and the reconstructed image is the directional edge enhancement of the distribution function on the slice. The Hilbert-transform is implemented by a 4*f* optical Fourier-transform correlator and a spatial filter consisting of a phase step of *π* radians. One important feature of this proposal is to perform a turn of 180° in the spatial filter at a certain value of the projection angle within the range $\left[0\xb0,\text{}360\xb0\right]$. The desired direction of enhancement can be chosen by the proper selection of such turning angle. We present both the mathematical modeling and numerical results.

©2011 Optical Society of America

## 1. Introduction

Optical Tomography (OT) is a form of computed tomography (CT) with the purpose of obtaining digital spatial distributions of some physical property lying within a given object by reconstructing images that is produced by light transmitted and scattered through an object [1]. When the probe beam travels through the object interior with no deviation (refractionless), the beam can be considered as composed of parallel rays. Since such a beam composed of only parallel rays, where each one of them remains parallel to the other, it is called parallel ray tomography. At the output, the probe is modified in a distribution known as projection, which is produced with a direction fixed to the object (projection angle,) and each projection depends on this angle. An important feature of this particular case consists in a description in which the use of Fourier methods is very useful. In this context, a parallel projection can be conveniently described by the Radon transform (RT), as it has been done in several fields, including computer aided x-ray tomography, radar imagery, geophysics and optics [2].

OT is used principally as a form of research in medical imaging, so it works best on soft tissues; for example, in imaging of breast and brain tissue [1]. However, those objects that modify only the phase in the probe beam are also of interest. These objects are known as phase-only objects and, specifically, when the optical path variations become much smaller than a wavelength, they are known as thin phase objects ($<\pi /3$); e.g., liquid flows, biological samples, thin films, etc. In OT of phase objects, the physical property that is used is the index of refraction. There has been a lot of research of phase objects using interferometric techniques for phase retrieval and in order to obtain a convenient algorithm for reconstruction. This method is known as interferometric tomography. For example, the reconstruction of flow fields in a simple structure with significant noise present is used in the measurement of asymmetric temperature fields [3]. In this proposal, there are four interferograms that are captured and processed with the use of automatic interferogram-processing software allowing us to obtain and reconstruct the projection data. We also use interferometric tomography to reconstruct the asymmetric three-dimensional temperature field generated by radiators and electronic chips [4]. In this situation, the fringe tracing method is used to process the interferograms, and considering the significant refraction present, which has been mentioned by other authors in order to reconstruct the index of refraction fields [5] or for the measurement of concentration profiles in the boundary layer formed at the cathode of an electrolytic cell containing ZnCI_{2} [6]. In a special case, an interfero-sinogram is formed in a two aperture common-path interferometer and afterwards, a phase-shifting is applied in order to obtain a sinogram and reconstruct its form by means of a back-projection algorithm [7]. In this proposal, oils, acetates, and glasses are used as a phase objects. Most of the techniques mentioned above obtain information from the interferograms and these have been oriented to get the distribution of slice functions. However, in many practical situations the phase variations are too large and the conventional interferometry methods are not practical. In order to overcome or minimize these deficiencies, other proposals to obtain some operation on the slice function have been made, such as obtaining the angular derivative in the slice function by using speckle interferometry [8], In this article, it is demonstrated that the angular derivative of projection data is proportional to the projection of angular derivative of slice function, which gives an edge enhancement on the reconstructed image. Another type of edge enhancement is obtained with the fractional derivative of the projection [9], since it is shown that it is proportional to directional derivative of the slice function. This technique could be experimentally implemented by using a $4{f}_{L}$ system of double Fourier-transform.

In the present manuscript, we discuss a new method to obtaining edge enhancement on the slices of thin phase objects by calculating the HT optical field, under the refractionless approximation. Therefore this proposal is made under the optical phase tomography of the parallel projections. We prove that the Hilbert transform (HT) of the projection data can be described analytically and that the inverse RT can also be established formally. Furthermore, we show that the reconstructed image is the directional HT of the slice distribution and that this image shows a directional edge enhancement. The optical HT is implemented by using a system of double Fourier-transform and a phase-step of *π* radians as a spatial filter in the frequency space, such as it has been proposed in several investigations [10–13].

## 2. Basic considerations

Let us consider a monochromatic plane wave $A(r)\phantom{\rule{.2em}{0ex}}\mathrm{exp}\phantom{\rule{.2em}{0ex}}\left[i\left(k\cdot r-\omega t\right)\right]$ which is traveling on direction given by wave vector $k=\left(2\pi /\lambda \right)\left(-\mathrm{sin}\phi ,\mathrm{cos}\phi ,0\right)$, where *λ* and *ω* are the wavelength and the frequency of light respectively, $i=\sqrt{-1}$ is the imaginary unit, **r** is a position vector, and *φ* is the azimuth rotation angle with respect to *x*-axis, as is illustrated in Fig. 1
. In rotated coordinates this plane wave can be rewritten as

*z*by traveling on the $(x,y)$ plane. Observing only one object slice at $z={z}_{0}$ (see Fig. 1b), these relations can be rewritten without

*z,*and in particular, the optical path can be expressed in a more convenient form, relating it with the RT [2] by means of

*φ*in the ${\stackrel{\vee}{f}}_{\phi}(p)$ indicate that the RT is applied on the slice function $f\left(x,y\right)$ when the projection angle

*φ*is kept fixed, while the projection coordinate

*p*is varied, so that under these conditions Eq. (3) are considered to be a sample of the RT, and they are known as a profile or a projection at

*φ*angle. As it is well known in tomographic theory, these projection data constitute the basis for the reconstruction of an object slice [1], and thus the index of refraction $f\left(x,y\right)$ can be reconstructed from a set of projections in the range $\left(0,\pi \right)$ by using the back-projection algorithm, which is the numerical implementation of the inverse RT [1].But, for a thin phase object and for only one slice at $z={z}_{0}$, Eq. (2) can be written as

*p*coordinate). A given object slice defines the function$f\left(x,y\right)$ at height $z={z}_{0}$ and the optical field leaving of the slice ${A}_{\phi}(p),$ given by Eq. (4). Lens

*L*

_{1}performs the FT of the wave leaving the object $\stackrel{~}{{A}_{\phi}}(w)=\Im \left\{{A}_{\phi}(p)\right\}$ located on the back focal plane (system’s frequency plane). Lens ${L}_{2}$ performs the inverse FT conventionally after a filter ${\stackrel{~}{h}}_{\phi}(w),$ which is placed on the frequency plane. The filter function consists of a glass plate which is covered on one half by a

*π*phase-shifting layer, which is the step phase, and it is mathematically approximated by the

*signum*function $\text{sign}(w)$, and for our convenience $-i\cdot \text{sign}(w)$. It is important to note that there is an ambiguity at $w=0$, however a adequate discussion about this ambiguity is out of the scope of this manuscript. Thus, on the image plane the transmittance is given by the convolution of input function and the system’s impulse response, ${h}_{\phi}(p)={\Im}^{-1}\left\{-i\text{sign}(w)\right\}=1/\pi p$. In others words, we have

The principal aim of this paper is first, to consider that the projection data ${\stackrel{\vee}{f}}_{\phi}(p)$ is obtained using Eq. (6). Afterwards, we analytically prove that when the back-projection algorithm is used, the HT of object slice on vertical direction is reconstructed, and this reconstruction results in an edge-enhancement along the vertical direction. Finally, we deduce a mathematical model to obtain the directional HT of the object slice and thus, the directional edge-enhancement of the object slice.

## 3. Theoretical analysis

We will now proceed to show that it is possible to obtain the HT of the object slice computing the HT of each projection data for all *φ* in the range $\left(0,2\pi \right)$. First, the projection or profile at angle *φ*, as described in Eq. (3), is denoted by

In the Fourier space, the same relation is given by

*φ*, $\stackrel{~}{f}\left(\mu ,\nu \right)={\Im}_{2D}\left\{f\left(x,y\right)\right\}$ is the FT of the slice function $f\left(x,y\right)$ with ${\Im}_{2D}\{\cdots \}$ denoting the bi-dimensional FT operator, and finally it is easy to prove that $\delta \left({w}_{\perp}\right)={\Im}_{2D}\left\{\delta \left(p\right)\right\}$, where $w=\mu \mathrm{cos}\phi +\nu \mathrm{sin}\phi $, and ${w}_{\perp}=-\mu \mathrm{sin}\phi +\nu \mathrm{cos}\phi $ are the conjugate variables of

*p*and ${p}_{\perp}$, both of them result from a coordinate rotation of

*μ*and

*ν*(the spatial frequencies) at the

*φ*angle. Equation (8) is known as the Fourier slice theorem, which means that the projection spectrum is equal to a sample of the slice function bi-dimensional spectrum along of the line ${w}_{\perp}=0$.

#### 3.1. HT of projections: Proof

We propose to calculate the HT of the projection ${\stackrel{\vee}{f}}_{\phi}(p)$ for every *φ* projection angle by using Eq. (6). The *φ* angle is carried out by rotating the object around the *z* axis and by considering that the reference system turns together the object, while the optical system is kept fixed; i. e., $\left(x,y\right)$ turn and $\left(p,{p}_{\perp}\right)$ are keep fixed. In this way, a modified sinogram can be generated, which is known as the Hilbert-sinogram. The HT operation is rewritten as:

In order to prove that Eq. (9) is a valid projection; both the symmetry and the zero-moment properties of the RT must be satisfied. Thus, ${\stackrel{\cap}{\stackrel{\vee}{f}}}_{\phi +\pi}(p)=\stackrel{\cap}{\stackrel{\vee}{{f}_{\phi}}}(-p)$and ${\int}_{-\infty}^{\infty}dp\stackrel{\cap}{\stackrel{\vee}{{f}_{\phi}}}\left(p\right)}=const.$, must be verified for each *φ*. The symmetry property states that

*φ*∈ (0,

*π*) and negative when $\phi \in \left(\pi ,2\pi \right)$. Due to this, the sign function is introduced as a factor in Eq. (9)

#### 3.2. Vertical edge-enhancement: particular case

Without loss of generality, after applying the convolution property of the FT, Eq. (12) can be written as

*φ*)sgn(

*w*), and $\stackrel{~}{\stackrel{\vee}{{f}_{\phi}}}\left(w\right)$ is the FT of the projection ${\stackrel{\vee}{f}}_{\phi}\left(p\right)$. Thereafter, by applying the Fourier slice theorem to Eq. (16) and changing to coordinates $\left(\mu ,\nu \right)$, we obtain the following relation

The Fourier inverse of this triple product results in a triple convolution, namely,

*f*(

*x*,

*y*), thus

Using Eq. (16), the filter on the Fourier space can be written as

#### 3.2.1. Numerical simulation 1

Figure 3
shows a numerical simulation of the vertical edge-enhancement in optical tomography, as has been studied in before. The images are presented in 8-bit gray levels. Column 2a of Fig. 3 shows an object’s slice$f\left(x,y\right)$ with $200\times 200$ pixels consisting of a rectangle of unit height (upper row), while in the lower row there is the corresponding function’s profile along a vertical raster line crossing by the rectangle’s center. Column 2b shows the sinogram of $f\left(x,y\right)$ as resulting from an algorithm implementing Eq. (3). Coordinate *p* has 200 data entries and the angular step used is 1.8°, so the projection number is 200 for $\phi \in \left(0,2\pi \right)$. This means that the sinogram has $200\times 200$ pixels. The upper row of column 2c shows the convolution of the RT of $f\left(x,y\right)$ and the system’s impulse response $1/\left(\pi p\right)$ (Hilbert-sinogram) with the corresponding sign compensation, as is indicated in Eq. (12). Note that the sign changes around $\phi =\pi $, in order to invert the sign filter. Under a filtered back-projection algorithm, we obtain the tomographic reconstruction, where the Hilbert-sinogram shown in Fig. 3c is used as an entrance sinogram, the upper rows of columns 2d-e show the results of both in a gray-levels plot and in a 3-D plot, respectively. The lower rows of the same columns show representative raster lines of the reconstructions. The line in (d) is a row of the reconstructed image as shown in the column (d) (is a row close to the central one), which appears to consists of noise around a constant. The line in (e) is a column close to the central one. This profile shows enhancements along the vertical direction.

#### 3.3. Directional edge-enhancement: general case

To describe the reconstruction process used in order to obtain edge enhancement along an arbitrary direction, we first define the HT along a given direction *α* as

*α*according with Eq. (22). The FT of Eq. (23) can thus be written as

What is in rectangular parenthesis is the FT of Eq. (22) and it is possible to prove that $-\mathrm{sgn}(\mu \mathrm{cos}\alpha +\nu \mathrm{sin}\alpha )=\Im \left\{\frac{\delta (-x\mathrm{sin}\alpha +y\mathrm{cos}\alpha )}{\pi (x\mathrm{cos}\alpha +y\mathrm{sin}\alpha )}\right\}$. In turn, Eq. (24) is the FT of the RT of the HT along direction indicated by *α* of the slice function $\Im \left\{{\Re}_{\phi}\left\{{\mathbf{H}}^{\alpha}\left\{f\left(x,y\right)\right\}\right\}\right\}$, and it can be rewritten by using rotated coordinates $\left(w,{w}_{\perp}\right)$ for example

This filter basically consists of the signum function and a change of sign given by the direction *α* and the *φ* projection angle through the cosine function. Therefore, the experimental implementation is viable for any given direction of enhancement by using the same scheme depicted in Fig. 2. For the case of $\alpha =\pi /2$, the situation described by Eq. (21) is achieved. Note that, for the case of $\alpha =0$, the enhancement direction becomes horizontal and the filter can be given by

By applying the inverse FT to Eq. (25), the following relation is obtained:

This expression shows the equivalence between the RT and the HT of $f(x,y)$ along the direction *α* and the HT of the RT of $f(x,y)$ with a sign factor. In other words,

As was noted before, in the case in which $\alpha =\pi /2$ the situation described in Eq. (20) is achieved. It is important to note that what is inside the brackets of the Radon operator on the left side of the equality ${\mathbf{H}}^{\alpha}\left\{f\left(x,y\right)\right\}$ is the slice function to be reconstructed.

#### 3.3.1. Numerical simulation 2

An experimental implementation of this HT operation at *α* angle can be done with the same optical system of Fig. 2 and the same phase-step filter, but with the difference of choosing its initial position according to the angle *α*. Such a position is the angle value at which the turning of the filter has to be done in order to fulfill the RT symmetry properties.

Column (a) of Fig. 4 shows two slices: a phantom (upper row) and a uniform ring (lower row). Column 3b presents respective sinograms and Column 3c the corresponding Hilbert-sinograms obtained by using Eq. (24) with $\alpha =0\xb0$. The sinogram and Hilbert-sinogram have 200 projections for $\phi \in \left(0,2\pi \right)$. The images in column 3d are their respective reconstructions, showing edge enhancement along the horizontal direction, as expected. Columns (e) and (f) show Hilbert-sinograms and reconstructions for $\alpha =45\xb0$, while columns (g) and (h) present similar images for the cases in which$\alpha =90\xb0.$

The expected direction of enhancement is to be seen. The digital image parameters are the same as those described in section of the simulation 1. In the Hilbert-sinograms, which are depicted in Figs. 4c, 4e, and 4d, the change of sign in the same projections is performed according to the factor $\mathrm{sgn}\left[\mathrm{cos}\left(\phi -\alpha \right)\right]$ of the relationship stated in Eq. (27). These changes of sign can also be associated with the fact of the inversion of the sign filter on the Fourier space of an $4{f}_{L}$ image selected according to choice the enhancement direction angle and the projection angle.

## 4. Conclusions

A theoretical analysis has been developed in order to obtain a directional edge enhancement of tomographic images using a refractiveless approximation (parallel projections) for thin phase objects. In this proposal, the Hilbert transform of each parallel projection from a given thin phase object slice is obtained first. This was carried out by spatial filtering, with a filter consisting of a phase step of *π* radians. In general, the filter has to be rotated 180° once around the optical axis at a certain projecting angle *α* (in practice, when $\alpha =90\xb0$ and one works within the range ${0}^{\xb0}\le \phi \le {180}^{\xb0}$. There is no need of a physical rotation because the usual back-projection algorithm virtually makes the rotation for ${180}^{\xb0}<\phi <{360}^{\xb0}$.) An important parameter is the projection angle value at which the filter has to be turned in the Fourier plane. Depending on such a value, the enhancement direction results at the same angle *α*. Thus, it is possible to choose the enhancement direction only selecting the projection angle value at which the turn is made. Numerical results for two different slices were presented, using three different values of *α*. An experimental implementation of this proposal seems thus to be possible and work on it is in progress.

## Acknowledgment

This work was partially supported by PROMEP under grant PROMEP /103.5/09/4544 and one of the authors (AMP) appreciates CONACYT’s under grant scholarship number 160260/160260

## References and links

**1. **G. S. Abdoulaev and A. H. Hielscher, “Three-dimensional optical tomography with the equation of radiative transfer,” J. Electron. Imaging **12**(4), 594–601 (2003). [CrossRef]

**2. **S. R. Deans, *The Radon Transform and Some of its Applications* (Wiley, New York. 1983).

**3. **D. Yan and S. S. Cha, “Computational and interferometric system for real-time limited-view tomography of flow fields,” Appl. Opt. **37**(7), 1159–1164 (1998). [CrossRef]

**4. **D. Wu and A. He, “Measurement of three-dimensional temperature fields with interferometric tomography,” Appl. Opt. **38**(16), 3468–3473 (1999). [CrossRef]

**5. **I. H. Lira and C. M. Vest, “Refraction correction in holographic interferometry and tomography of transparent objects,” Appl. Opt. **26**(18), 3919–3928 (1987). [CrossRef] [PubMed]

**6. **S. Cha and C. M. Vest, “Tomographic reconstruction of strongly refracting fields and its application to interferometric measurement of boundary layers,” Appl. Opt. **20**(16), 2787–2794 (1981). [CrossRef] [PubMed]

**7. **C. Meneses-Fabian, G. Rodriguez-Zurita, and V. Arrizón, “Optical tomography of transparent objects with phase-shifting interferometry and stepwise-shifted Ronchi ruling,” J. Opt. Soc. Am. A **23**(2), 298–305 (2006). [CrossRef]

**8. **C. Meneses-Fabian, G. Rodriguez-Zurita, R. Rodriguez-Vera, and J. F. Vazquez-Castillo, “Optical tomography with parallel projection differences and Electronic Speckle Pattern Interferometry,” Opt. Commun. **228**(4-6), 201–210 (2003). [CrossRef]

**9. **G. Rodríguez-Zurita, C. Meneses-Fabián, J.-S. Pérez-Huerta, and J.-F. Vázquez-Castillo, “*“Tomographic directional derivative of phase objects slices using 1-D derivative spatial filtering of fractional order ½*,” ICO20,” Proc. SPIE **6027**, 410–416 (2006).

**10. **K. Sendhil, C. Vijayan, and M. P. Kothiyal, “Spatial phase filtering with a porphyrin derivative as phase filter in an optical image processor,” Opt. Commun. **251**(4-6), 292–298 (2005). [CrossRef]

**11. **J. A. Davis, D. E. McNamara, and D. M. Cottrell, “Analysis of the fractional hilbert transform,” Appl. Opt. **37**(29), 6911–6913 (1998). [CrossRef]

**12. **J. A. Davis, D. E. McNamara, D. M. Cottrell, and J. Campos, “Image processing with the radial Hilbert transform: theory and experiments,” Opt. Lett. **25**(2), 99–101 (2000). [CrossRef]

**13. **J. A. Davis and M. D. Nowak, “Selective edge enhancement of images with an acousto-optic light modulator,” Appl. Opt. **41**(23), 4835–4839 (2002). [CrossRef] [PubMed]

**14. **F. Zernike, “How I discovered phase contrast,” Science **121**(3141), 345–349 (1955). [CrossRef] [PubMed]

**15. **J. Glückstad, P. C. Mogensen and R. L. Eriksen, “The generalized phase contrast method and its applications,” DOPS-NYT **1–2001**, 49–54.

**16. **S. Fürhapter, A. Jesacher, S. Bernet, and M. Ritsch-Marte, “Spiral phase contrast imaging in microscopy,” Opt. Express **13**(3), 689–694 (2005). [CrossRef] [PubMed]