Fast calculation of spherical computer generated hologram using spherical wave spectrum method

Boaz Jessie Jackin; Toyohiko Yatagai

doi:10.1364/OE.21.000935

1. Introduction

Holograms are capable of recording and reproducing all three dimensional information like motion parallax, accomodation, occlusion, convergence and so on. For three-dimensional information storage and retrieval(of an object) to be complete, recording and reconstruction is to be done on all sides(directions) of an object, i.e, for 360° on both azimuth and polar directions. This is achieved by recording and reconstructing on a spherical surface surrounding the object, which is called as spherical holography. However available optical techniques and numerical methods have restricted holography to be done on plane and cylindrical surfaces only. The restriction due to optical techniques arises from achieving coherent illumination of object and reference wave on a spherical surface. One solution is to model the object on the computer and numerically generate the hologram. But still the following issues prevail, i) The unavailability of fast calculation methods or availability of only approximated fast calculation methods for spherical CGH and ii) the need to illuminate the whole spherical surface using a single reference wave front for optical reconstruction. However the availability of optic fibers and developments in multidimensional lasers [1, 2] that include microball lasers or lasing microspheres have provided hope in realizing spherical holography. Motivated by these appealing facts it was decided to make some developments to numerical methods for spherical holography. Accrodingly, this research work is an attempt to devise a new fast calculation method for computer generation of spherical holograms.

Wave propagation in the paraxial regime is defined by the Rayleigh Sommerfeld diffraction integral [3]. This formula enables us to calculate the optical field(complex amplitude) at an arbitrary point or surface(hologram) knowing the radiating source point or surface(object). This method also known as the direct integration method and is slow due to large sampling points in object and hologram. To speed up computation the Rayleigh Sommerfeld integral is expressed as convolution integral or as angular spectrum of plane waves and then evaluated using fast Fourier transform(FFT). This speeds up the calculation considerably, however demands the object surface and hologram surface to satisfy shift invariance. Accordingly, to generate a plane hologram using FFT, the object should also be a parallel plane surface or set of parallel plane surfaces. If the hologram surface is non planar or not parallel to the object surface then the direct integration method has to be used, which consumes O(N⁴) or O(N³) operations for N sampling points. But later on smart ways to use FFT for calculations even when the hologram surface is non planar or parallel were devised. Tommasi et al., [4] developed a fast computation method where the hologram is a tilted plane with respect to the object surface. Sakamoto et al.[5] took advantage of the shift invariance in rotation and developed a method to compute cylindrical hologram for a plane object using FFT. Later the method was extended to a volume object by slicing it into planar segments [6]. Yamaguchi et al.[7] proposed a fast computation method for cylindrical holograms by segmenting the object and hologram surface into elemental plane and parallel surfaces. Sando et al.[8] came up with a brilliant idea by choosing the object to be a concentric cylindrical surface with the hologram surface which makes the system completely shift invariant. They derived the point spread function(PSF) for this configuration and developed the convolution integral to calculate wave progagation from one cylindrical surface(object) to other(hologram).The convolution method was evaluated using FFT, making it a fast computation method. Later Jackin et al.[9], devised a new and fast computation method for the same problem by deriving the Transfer Function and defining wave propagation in spectral domain. The advantage of this method is, it uses one less FFT operation for its computation compared to the convolution method. Then the method was extended to a volume object by slicing in to cylindrical surfaces [10].

Fast computation solutions for spherical computer generated hologram emplyoing PSF(convolution method) was proposed by Tachiki et al.[11]. Though the object was assumed to be a concentric spherical surface with the hologram surface, shift invariance does not exist due to unequal sampling points on a spherial grid(ie., the grid points are more crowded at the poles). To facilitate fast calculation using FFT, an approximation to the convolution integaral was proposed, which forced the PSF to be spatially invariant. This research work is an attempt to achieve the same but by deriving the transfer function and defining wave propagation in spectral domain. In other words, it is to find a spectral propagation formula (for spherical surfaces in spherical coordinates) analogous to the angular spectrum formula(for plane surfaces in cartesian coordinates). The expected advantages compared to the earlier method are faster calculation and no approximations. Section 2 explains the theoretical development of the calculation method. The numerical procedure for evaluating the formula is described in section 3. Simulation results to verify the proposed method are presented in Section 4. Section 5 summarises the work with possible applications and future improvements to this work.

2. Theoretical background

We will describe points on the sphere by their latitude θ ∈ [−π/2, π/2] and longitude ϕ ∈ [0, 2π] as shown in Fig. 1. Throughout this work we will be dealing with square-integrable functions that span a space L²(S) on the unit sphere S. Fast computation algorithms take advantage of the cyclic and peroidic properties of the transformation kernel for fast calculations. The cyclic and periodic properties exists only if the system is shift invariant between the transformation planes for that operation. Which means that, the object and hologram surface should be shift invariant in order to devise a fast computation scheme. To achieve this the object and hologram surface are chosen to be concentric spherical surfaces, so that they can remain shift invariant in rotation along ϕ and θ directions. But these surfaces are defined on a spherical grid, where the sampling points are more dense at the poles than at the equator. Hence shift invariance is not satisfied. However in the hologram generation process, the hologram and object are band-limited functions on L²(S) space. The band-limited functions on L²(S) space has a very useful and important property that any rotated version of a band-limited function is also a band-limited function with the same bandwidth[12, 13]. Thus, they are refered to as having uniform resolution, at all points on the sphere, meaning that they are shift-invariant. The triangular truncation and Gauss-legendre quardature method that occurs in the transformation operation is responsible for this property(explained in the next section). Hence the system shown in Fig. 1 does have shift invariance relationship between object and hologram surfaces and approves the possiblity of fast computation formula. With this assurance we proceed to develop the fast calculation method for computer generated spherical holography starting from the basic equations of electromagnetism.

Fig. 1 Geometry of the problem.

Download Full Size | PDF

An electromagnetic field is defined by Maxwell’s equations and its propagation by the Helmholtz wave equation. Hence for any particular system the complex amplitude of a propagating wave at any instance of time and any where in space, can be found by solving the wave equation, applying its constraints and conditions. Accordingly for the system shown in Fig. 1 the solution can be derived starting from the wave equation as follows, The time dependent vector wave equation u(r,θ, ϕ) is expressed as

\nabla^{2} u - ε μ \frac{\partial^{2} u}{\partial t^{2}} = 0

where r is the radius of the spherical surface of interest and θ, ϕ represent the azimuthal and polar angles in the surface. The Laplacian operator ∇² in spherical co-ordinates is defined as

\nabla^{2} = \frac{1}{r^{2}} \frac{\partial}{\partial r} (r^{2} \frac{\partial}{\partial r}) + \frac{1}{r^{2} sin θ} \frac{\partial}{\partial θ} (sin θ \frac{\partial}{\partial θ}) + \frac{1}{r^{2} {sin}^{2} θ} \frac{\partial^{2}}{\partial ϕ^{2}}

When the wave is non polarised and the medium of propagation is homogeneous then the vector nature of the function “u” can be discarded. The scalar wave equation given in spherical coordinates becomes

\frac{1}{r^{2}} \frac{\partial}{\partial r} (r^{2} \frac{\partial u}{\partial r}) + \frac{1}{r^{2} sin θ} \frac{\partial}{\partial θ} (sin θ \frac{\partial u}{\partial θ}) + \frac{1}{r^{2} {sin}^{2} θ} \frac{\partial^{2} u}{\partial ϕ^{2}} - \frac{1}{c^{2}} \frac{\partial^{2} u}{\partial t^{2}} = 0

The solution of Eq. (3) can be found by separation of variables [14, 15], which can be expressed as shown below

u (r, θ, ϕ, t) = R (r) Θ (θ) Φ (ϕ) T (t)

The separable variables obey the following four ordinary differential equations.

\frac{d^{2} Φ}{d ϕ^{2}} + m^{2} Φ = 0

\frac{1}{sin θ} \frac{d}{d θ} (sin θ \frac{d Θ}{d θ}) + [n (n + 1) - \frac{m^{2}}{{sin}^{2} θ}] Θ = 0

\frac{1}{r^{2}} \frac{d}{d r} (r^{2} \frac{d R}{d r}) + k^{2} R - \frac{n (n + 1)}{r^{2}} R = 0

\frac{1}{c^{2}} \frac{d^{2} T}{d t^{2}} + k^{2} T = 0

The solution to all the above equations are derived in by Arfken [15]. Only the final results are used in this research work. The solution to azimuthal Eq. (5) is

Φ (ϕ) = Φ_{1} e^{im ϕ} + Φ_{2} e^{- im ϕ}

where m must be an integer so that there is continuity and periodicity in Φ(ϕ). Φ₁ and Φ₂ are constants.

The solution to polar Eq. (6) is

Θ (θ) = Θ_{1} P_{n}^{m} (cos θ) + Θ_{2} Q_{n}^{m} (cos θ)

where,

P_{n}^{m}

and

Q_{n}^{m}

are the associated Legendre polynomials of first and second kind respectively and Θ₁ and Θ₂ are constants.

Q_{n}^{m}

is not finite at the poles where cos(θ) = ±1, so this solution is discarded(Θ₂ = 0).

For the radial differential equation Eq. (7) the solutions are

R (r) = R_{1} h_{n}^{(1)} (k r) + R_{2} h_{n}^{(2)} (k r)

where,

h_{n}^{(1)}

and

h_{n}^{(2)}

are the Spherical Hankel functions of the first and second kind respectively. Since we are interested only in outgoing wave, we can neglect the second term (R₂ = 0)

The solution to Eq.(8) is

T (ω) = T_{1} e^{i ω t} + T_{2} e^{- i ω t}

Again here we assume T₂ = 0 due to the convention of time

The angle functions in Eq. (9) and Eq. (10) are conveniently combined into a single function called as spherical harmonic $Y_{n}^{m}$ [14, 15] defined by

Y_{n}^{m} (θ, ϕ) \equiv {(- 1)}^{m} \sqrt{\frac{(2 n + 1) (n - m)!}{4 π (n + m)!}} P_{n}^{m} cos (θ) e^{im ϕ}

where the quantity

{\bar{P}}_{n}^{m} = \sqrt{\frac{(2 n + 1) (n - m)!}{4 π (n + m)!}} P_{n}^{m} (cos θ)

is known as the orthonormalized associated Legendre polynomial. The term (−1)^m is called as the Condon-Shortley phase. Hence the spherical harmonics can also be represented in a short form as shown below(neglecting the Condon-Shortley phase)

Y_{n}^{m} (θ, ϕ) = {\bar{P}}_{n}^{m} (cos θ) e^{im ϕ}

Combining all the above equations, the travelling wave solution to Eq. (3) can be represented as

u (r, θ, ϕ, ω) = \sum_{n = 0}^{\infty} \sum_{m = - n}^{n} A_{m n} (ω) h_{n} (k r) Y_{n}^{m} (θ, ϕ)

The radiated field is completely defined when the coefficient A_mn is determined. This is achieved by using the orthonormal property of the spherical harmonics. Assume that the wave-field u(r,θ, ϕ) is known on a sphere of radius r = a. We also drop the time variable(which is not important) for simplicity.Now multiplying each side of Eq. (16) (evaluated at r = a) by $Y_{n}^{m} {(θ, ϕ)}^{*}$ and integrating over the sphere, gives

A_{m n} = \frac{1}{h_{n} (k a)} \int_{\frac{- π}{2}}^{\frac{π}{2}} \int_{0}^{2 π} u (a, θ, ϕ) Y_{n}^{m} {(θ, ϕ)}^{*} sin (θ) d θ d ϕ

where, dΩ = sin(θ)dθdϕ′, is the solid angle for integrating on a sphere. Inserting the expression for A_mn back into Eq. (16), we get,

u (r, θ, ϕ) = \sum_{n = 0}^{\infty} \sum_{m = - n}^{n} Y_{n}^{m} (θ, ϕ) ([\int_{\frac{- π}{2}}^{\frac{π}{2}} \int_{0}^{2 π} u (a, θ^{'}, ϕ^{'}) Y_{n}^{m} {(θ^{'}, ϕ^{'})}^{*} d Ω^{'}] \frac{h_{n} (k r)}{h_{n} (k a)})

Hence the wavefield at any spherical surface u(r,θ, ϕ) can be calculated knowing the wave-field at u(a,θ′, ϕ′).

For easy understanding and interpretation of Eq. (18), it is compared with the well known angular spectrum of plane waves equation [3] which is given by

u (x, y, z) = \frac{1}{4 π^{2}} \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} e^{i (k_{x} x + k_{y} y)} d k_{x} d k_{y} ([\int_{- \infty}^{\infty} \int_{- \infty}^{\infty} u (x^{'}, y^{'}, 0) e^{- i (k_{x} x^{'} + k_{y} y^{'})} d x^{'} d y^{'}] e^{i k_{z} z})

It is well known that in Eq. (19) the quantity within square brackets represents the source wave-field decomposed into spectrum of plane waves in (k_x, k_y). The propagation of the spectrum is defined by the quantity e^ik_zz which is known as the transfer function. The propagated spectrum is recomposed into the wavefield at destination by the integral over k_x, k_y. Thus this equation defines wave propagation from one plane surface u(x′, y′,0) to another surface u(x,y,z). When comparing Eq. (18) with Eq. (19), the following can be deduced

Wave field in (θ, ϕ) at a radiating spherical surface of radius r = a is decomposed into its wave spectra in (m,n) defined by $U_{m n} (a) = \int u (a, θ, ϕ) Y_{n}^{m} {(θ, ϕ)}^{*} d Ω$
The decomposed wave components(spectra) are expressed by Eq. (13) which is composed of a travelling wave component in ϕ and defined by e^imϕ and a standing wave component in θ given by $P_{n}^{m} (cos θ)$
The decomposed wave components in this system can be named as spherical wave components in analogous to plane wave components for planar system. Similarly Eq. (20) can be termed as Spherical wave spectrum as analogus to Angular spectrum of plane waves.
The wavenumbers k_x and k_y are imitated by the quantities m/a and n/a. Hence we can refer to the spherical wave spectrum as a k-space spectrum due to this anlogy.
Equation (20) can be viewed as a forward Fourier transform using $Y_{n}^{m} (θ, ϕ)$ as the basis function. In other words, the spectral space for spherical system is spanned by Spherical harmonic functions( $Y_{m}^{n} (θ, ϕ)$ ). This is also termed as Spherical Harmonic Transform (SHT)
The propagation of spherical wave spectrum from one spherical surface of radius a to another of radius r is given by $U_{m n} (r) = \frac{h_{n} (k r)}{h_{n} (k a)} U_{m n} (a)$
Hence the quantity h_n(kr)/h_n(ka) can be refered to the Transfer function(TF) for spherical system, as opposed to the quantity e^ik_zz for planar system.
The inverse spherical harmonic transform(ISHT) that recomposes the wave field back from the spectrum is given by $u (r, θ, ϕ) = \sum_{n = 0}^{\infty} \sum_{m = - n}^{n} U_{n m} (r) Y_{n}^{m} (θ, ϕ)$

The transfer function is crucial since it completely defines propagation and hence it is worth discussing some of its properties. During propagation amplitude and phase of the spectral components change with distance as defined by the transfer function. What is most important is the rate of change of phase of the transfer function which determines the sampling requirements. Accordingly the plot shown in Fig. 2 reveals that the phase change increases with increasing orders ’n’ of the transfer function (The plot was generated for wavelength 100μm, radius 10 and 0.5 cm, and upto 256 orders). Hence sampling requirements will be satisfied if highest order ’n’ of the transfer function is sampled according to Nyquist criteria. It can also be understood that, the rate of phase change becomes higher as the distance of propagation increases. This will demand very large sampling and also increase the numerical errors. It is worthy to note here that the spherical Hankel functions are asymptotic in nature. In the far field the spherical Hankel functions can be approximated by their asymptotic expressions which is as shown in Eq. 23.

h_{n}^{(1)} (x) = {(- i)}^{n + 1} \frac{e^{i x}}{i x}

This could yield a formula analogous to the farfield Fresnel diffraction formula. However it requires systematic development of theory with proper analysis of approximations and are not discussed in this paper. But this can be considered as a future work to the proposed method.

Fig. 2 Plot of phase(in radians) of the transfer function for increasing order (n).

Download Full Size | PDF

Thus the devised formula which is analogus to the angular spectrum of plane waves(AS) formula defines wave propagation between spherical surfaces. The numerical procedure to implement the devised formula is discussed in the next section.

3. Numerical computation

The numerical computation of AS method heavily depends on the FFT operations for which a lot of tools and methods are available. But the numerical computation of the proposed method heavily depends on the SHT operations. Fast computation was guaranteed from the theory and from the geometry of the system. Now a numerical procedure that takes advantage of this is required. Fortunately lot of fast computation numerical methods have been reported for SHT and this research work could make use of it. Since FFT and numerical computation of AS method are well understood, this section intends to introduce SHT and numerical computation of proposed method in close anology to the former.

Continuing with the comparision from previous section, the numerical computation of wave propagation for planar and spherical systems according to Eq. (18) and Eq. (19) can be represented respectively as

u (r, θ, ϕ) = ISHT [SHT (u (a, θ, ϕ)) \times T F_{s}]

u (x, y, z) = IFFT [FFT (u (x, y, 0)) \times T F_{c}]

where, FFT[...] and IFFT[...] denote the foward and inverse fast Fourier transform operations, while SHT[...] and ISHT[...] denote the forward and inverse spherical harmonic transform operations. TF_s and TF_c are the transfer functions for spherical and planar wave propagations as defined by Eq. (18) and Eq. (19) respectively.Comparison reveals that the computation method (for spherical system) is analogous to the angular spetrum of plane waves method (for plane system) except the fact that, the Fourier Transform is replaced by the Spherical harmonic transform. Hence evaluation of SHT becomes the key, the others being only basic mathematical evalutations. Spherical harmonic transform have been studied extensively and fast computation algorithms and optimisation methods have been proposed. A method that requires only O(N²logN) operations for N grid points as opposed to the standard N³ operations was proposed by Chien et al.[12]. They imposed truncations on the spectral components and used fast multipole method and fast Fourier transform for evaluation. Their method is called as “spectral truncation method” and the errors due to truncation were well within acceptable limits. Later on Healey et al.[16] could achieve the same using O(N(logN))² operations. They took advantage of the recursive properties of associated Legendre polynomial for fast computation. This method was found to be more efficient and hence was choosen for numerical evaluation the Spherical harmonic transforms in this work. A brief outline of the numerical evaluation is presented. For more details please refer to Healey [16] and Driscoll [13].

A 2D DFT is computed by separating the transformation kernel into its variables, and evaluating it as column transform followed by a row tranform. Similarly the spherical harmonic transform kernel $Y_{n}^{m}$ given by Eq. (15) is also variable separable and can be separated into ϕ component and θ component. Accordingly the transformation given by Eq. (20) can be represented as shown below

U_{m n} (r) = \int_{- π / 2}^{π / 2} (\int_{- π}^{π} u (r, θ, ϕ) e^{im ϕ} d ϕ) {\bar{P}}_{n}^{m} (cos θ) d θ

First, the quantity within the round brackets alone is to be computed which is nothing but a Fourier transform operation. The Fourier coefficients U^m(θ) are evaluated for m = −N,...,N as shown below

U_{m} (θ) = \int_{- π}^{π} u (r, θ, ϕ) e^{i m ϕ} d θ

= \frac{1}{I} \sum_{i = 1}^{I} u (r, θ, ϕ_{i}) e^{im ϕ_{i}}

where ϕ_i = 2πi/I for i = 1,...,I. The equispaced longitudes ϕ_i enables the use of fast Fourier transform.

Second, the Legendre transform of the Fourier coefficients U_m(θ) is to be evaluated for |m| ≤ n ≤ N. This is done using the Gaussian-Legendre quadrature as shown below,

U_{n m} = \int_{- π / 2}^{π / 2} U_{m} (θ_{j}) {\bar{P}}_{n}^{m} (cos θ) sin θ d θ

= \sum_{j = | m |}^{N} U_{m} (θ_{j}) {\bar{P}}_{n}^{m} (cos θ_{j}) w_{j}

where θ_j and w_j are respectively the Gauss nodes and weights and are calculated using the Fourier-Newton method as described by Swarztrauber [17]. The Gauss-Legendre quadrature replaces the integral by the sum. The fact that the summation runs only from |m| to N is refered to as the triangular truncation. The use of Gauss-Legendre quadrature method redistributes θ into Gaussian nodes θ_j. This along with the triangular truncation are responsible for uniform resolution on the latitudinal points. Again the triangular truncation along with the recurrance property of Legendre polynomial helps to achieve fast computation.

In the similar way the inverse spherical harmonic transform can be represented as

u (θ, ϕ) = \sum_{m = - N}^{N} (\sum_{n = | m |}^{N} U_{n m} P_{n}^{m} (cos θ)) e^{im ϕ}

The inverse transform also follows the same procedure and is computed in two steps but in the reverse order (ie, legendre transform first and Fourier transform next) as shown below

U_{m} (θ) = \sum_{n = | m |}^{N} U_{n m} {\bar{P}}_{n}^{m} (cos θ)

u (θ, ϕ) = \sum_{m = - N}^{N} U_{m} (θ) e^{im ϕ}

Thus using this numerical procedure fast computation of wave propagation in spherical computer generated holograms is achieved, which uses only O(N(logN)² operations for the SHT computation.

There are a lot of software tools available on the internet to do the SHT operation. Most of them are tuned and dedicated for a geophysical process which requires only real SHT but holography requires complex SHT. However the package SHTools by Mark Wieczorek [18] could do a complex SHT operation in Fortran language. Though not tested by us this is the best one available that suits the work related to this paper. Hence we recommend using this if it is required to quickly reproduce the work in this paper. The next section describes the testing and verification of this numerical procedure using simulations.

4. Simulation results

The system considered for simulation experiments is shown in Fig. 1. The object(O(a,θ, ϕ)) is a spherical surface of radius 1cm and the hologram(H(r,θ, ϕ)) is another concentric spherical surface of radius 10cm. The reference is considered to be a virtual source emitting spherical waves from the center, ie., the wavefield due to reference has same phase and amplitude throughout the hologram plane. This is similar to using a plane reference wave with normal incidence in plane holography.

It is first required to verify the correctness of the method by expecting it to reproduce the already known diffraction results. Hence the proposed method is subjected to generate already reported diffraction patterns. For this the diffraction pattern reported by Tachiki et al [11] for spherical surfaces is used as reference. Accordingly, a simple object was chosen which is a spherial surface with two irradiating points at ϕ = −π/2 and ϕ = π/2, as shown in Fig. 3. The object and hologram are composed of 256 pixels in the longitudinal (north-south) direction and 512 pixels in the latitudinal(east-west) direction. The wavelength was chosen to be λ = 100μm, in order to reduce sampling requirements and visiblity of fringes. The procedure for numerical generation of hologram using the proposed method is expressed in an abstract form as shown below.

AmplitudeHologram = {| (ISHT [SHT (Object) \times T F]) + ISHT [SHT (Reference) \times T F] |}^{2}

Then a hologram for the same object was simulated using the the well known direct integration formula defined in Eq. (34).

H (r, θ, ϕ) = \iint \frac{O (θ^{'}, ϕ^{'}) exp (i k L)}{L} d x d y

where,

L = \sqrt{r^{2} + a^{2} - 2 r a [sin (θ) sin (θ^{'}) + cos (θ) cos (θ^{'}) cos (ϕ - ϕ^{'})]}

Hologram = {| H_{object} (r, θ, ϕ) + H_{reference} (r, θ, ϕ) |}^{2}

The simulation results are shown in Fig. 4. The pattern generated by the proposed method matches with the one generated by direct integration method. However the distribution of brightness and contrast across the pattern is constant for the direct integration method while it decreases gradually from the center for the proposed method. This inconsistency can be explained as follows. The direct integration formula Eq. (34) is the Rayleigh-Sommerfeld diffraction formula [11] without the obliquity factor. The obliquity factor is the cosine of the angle between the normal of the radiating surface to the direction of the observation point. This is responsible for the distribution of light intensity based on the angle (i.e,more bright at the center and gradually decreases outwards and no radiation backwards). However the spectral method which is the solution to the boundary value problem of the wave equation, incorporates the obliquity factor and hence the brightness and contrast varies radially. More over the obliquity factor does not alter the phase of the traveling wave which inturn does not affect interference pattern and hence guarantees a fair comparison.

Fig. 3 Object.

Download Full Size | PDF

Fig. 4 Computed hologram(intensity) using a)proposed method and b)direct integration.

Download Full Size | PDF

Then the proposed method is tested to see whether it obeys the fundamental laws of diffraction and interference. In other words, it is to verify that it has the same qualities as angular spectrum method. For this two simulation experiments for qualitative analysis was performed. First it was intended to analyze the change in interference pattern with change in wavelength which is a fundamental law. Accordingly for the object shown in Fig. 3, the hologram was computed for wavelengths varying as a)150μm, b) 200 μm, c) 250 μm, d) 300 μm, e)350μm and f)400μm respectively. The corresponding holograms generated are shown in Fig. 5. As expected the fringe density decreases with increase in wavelength. Second another fundamental law which is the change in interference pattern with the change in distance between the coherent sources was verified. This also corresponds to the youngs double slit experiment. Accordingly, the hologram is computed for varying position of the pair of the point sources on the spherical surface. This setup is similar to the young double slit experiment. The positions of point sources were set to be a)(ϕ = −π/6, ϕ =π/6), b)(ϕ = −π/8, ϕ =π/8), c) (ϕ = −π/16, ϕ =π/16) and d) (ϕ = −π/32, ϕ = π/32) as shown in Fig. 6(a)–(d) respectively. The corresponding computed hologram pattern are shown in Fig. 6(e)–(h). As expected the fringe density decreases with decrease in distance between the point sources. Since the diffraction formula agrees well with the fundamental laws of diffraction it is confirmed that it behaves like the any other diffraction formula and can be used to simulate wave propagation. Moreover the above mentioned results also reveals that the wavefield on the spherical surface was computed correctly and as expected.

Fig. 5 Computed hologram(intensity) for wavelength a) 150μm, b) 200μm, c) 250μm, d) 300μm, e) 350μm, f) 400μm.

Download Full Size | PDF

Fig. 6 Object with point sources at a)(θ = −π/6, θ = π/6), b)(θ = −π/8, θ = −π/8), c) (θ = −π/16, θ = π/16), d) (θ = −π/32, θ = π/32) and their corresponding hologram pattern(intensity).

Download Full Size | PDF

Since the theory is developed in context to computer generated holography, it is mandatory to verify its applicability to the same. For this it was decided to perform spherical hologram generation and then reconstruction from the hologram on the computer using the proposed method. The spherical object for which the hologram is to be made is shown in Fig. 7. The object and hologram were composed of 256 × 512 pixels. The wavelength for simulation was chosen to be λ = 30 μm. The generated hologram is shown in Fig. 8.

Fig. 7 Object.

Download Full Size | PDF

Fig. 8 Hologram(intensity).

Download Full Size | PDF

From this hologram the object was reconstructed back onto the original spherical surface using Eq. 18. This means that we are attempting to reconstruct a real image of the hologram and not the virtual one. For this the conjugate of the reference beam has to be used. Accordingly the numerical procedure of the for reconstruction is expressed in an abstract form as shown below.

Reconstruction = {| ISHT [SHT (Hologram \times Conjugate [Reference]) \times T F] |}^{2}

The reconstruction result is shown in Fig. 9. The reconstruction matches exactly with the object chosen. As mentioned earlier, the object and hologram are square integrable band limited functions on a closed surface. Hence a rotated(shifted in theta or phi) version of a the object or hologram will produce a rotated version of the reconstruction. The wave propagation calculation requires O(N(logN))² operations for N sampling points and hence it is a fast computation formula. The calculations were executed using a scripted language-python in a Dell Precision T7400 machine with 12 GB of RAM memory. A comparision of calculation time for the direct integration, convolution and spectral methods is shown in Table 1.

Fig. 9 Reconstruction.

Download Full Size | PDF

Table 1. Comparision of calculation speed

View Table

5. Concluding remarks

The solution to Helmholtz wave equation in spherical coordinates is derived using variable separable method. The spherical wave spectrum and transfer function were defined from the solution. A formula for computing the wave propagation from irradiating spherical surfaces is devised. A fast computation method that evaluates the wave propagation formula in O(N(logN))² operations was suggested. The proposed method was tested for correctness using simulated experiments. Generation and reconstruction of a spherical hologram for a spherical object was also successful. Hence a new and fast computation method for computer genetated spherical holograms is realized.

This method can be extended by considering it as a farfield diffraction problem by utilizing its asymptotic form as shown in Eq. (23). This will be analogous to the Fresnel diffraction formula which will take this development more close to practical applications. More over other similar problems (in acoustics and geophysics) have shown that far field diffraction pattern is directional and not uniformly distributed over a spherical surface. Hence we do not need a spherical detector (which is not available now) or scan over the whole spherical surface in order to holographically image a spherical volume. This interesting fact takes it even more closer to realizing practical applications. The availability of optic fibers will also be a great aid in this process. When considering holographic display, printing a spherical hologram on a flexible film and reconstructing using a laser will be a challenging task. However with the availability of high precision lithographic machines and optic fibers this is not an impossible task. It is also worth refereing to Horvath [19] for more interesting facts on the importance of this research work to the future of lasers. Hence the days when computer generated spherical holography will have an application is not so far. We believe this research will be a small step in this regard and will aid the development of computer generated spherical holography.

References and links

1. C. Gmachl, F. Capasso, E. E. Narimanov, J. U. Nockel, A. D. Stone, J. Faist, D. L. Sivco, and A. Y. Cho, “High power directional emission from microlasers with chaotic resonators,” Science 280, 1544–1545 (1998). [CrossRef]

2. N. M. Lawandy and R. M. Balachandran “Random laser?,” Nature 373, 204, (1995). [CrossRef]

3. J. W. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1968).

4. T. Tommasi and B. Bianco, “Computer-generated holograms of tilted planes by a spatial frequency approach,” J. Opt. Soc. Am. A 10, 299–305 (1993), http://www.opticsinfobase.org/josaa/abstract.cfm?URI=josaa-10-2-299. [CrossRef]

5. Y. Sakamoto and M. Tobise, “Computer generated cylindrical hologram,” in Practical Holography XIX: Materials and Applications Tung H. Jeong and Hans I. Bjelkhagen, eds., Proc. SPIE 5742, 267–274 (2005).

6. A. Kashiwagi and Y. Sakamoto, “A fast calculation method of cylindrical computer-generated holograms which perform image-reconstruction of volume data,” in Adaptive Optics: Analysis and Methods/Computational Optical Sensing and Imaging/Information Photonics/Signal Recovery and Synthesis Topical Meetings on CD-ROM OSA Technical Digest (CD) (Optical Society of America, 2007), paper DWB7, http://www.opticsinfobase.org/abstract.cfm?URI=DH-2007-DWB7.

7. T. Yamaguchi, T. Fujii, and H. Yoshikawa, “Fast calculation method for computer-generated cylindrical holograms,” Appl. Opt. 47, D63–D70 (2008), http://www.opticsinfobase.org/ao/abstract.cfm?URI=ao-47-19-D63. [CrossRef] [PubMed]

8. Y. Sando, M. Itoh, and T. Yatagai, “Fast calculation method for cylindrical computer-generated holograms,” Opt.Express 13, 1418–1423 (2005), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-13-5-1418. [CrossRef] [PubMed]

9. B. J. Jackin and T. Yatagai, “Fast calculation method for computer-generated cylindrical hologram based on wave propagation in spectral domain,” Opt.Express 18, 25546–25555 (2010), http://www.opticsinfobase.org/oe/abstract.cfm?uri=oe-18-25-25546. [CrossRef] [PubMed]

10. B. J. Jackin and T. Yatagai, “360° reconstruction of a 3D object using cylindrical computer generated holography,” Appl. Opt. 50, H147–H152 (2011), http://www.opticsinfobase.org/ao/abstract.cfm?URI=ao-50-34-H147. [CrossRef] [PubMed]

11. M. L. Tachiki, Y. Sando, M. Itoh, and T. Yatagai, “Fast calculation method for spherical computer-generated holograms,” Appl. Opt. 45, 3527–3533 (2006), http://www.opticsinfobase.org/ao/abstract.cfm?URI=ao-45-15-3527. [CrossRef] [PubMed]

12. R. J. Chien and B. K. Alpert, “A fast spherical filter with uniform resolution,” J. Comput. Phys. 136, 580–584 (1997). [CrossRef]

13. J. R. Driscoll and D. M Healy, “Computing Fourier transfroms and convolutions on the sphere,” Adv. Appl. Math. 15, 201–250 (1994). [CrossRef]

14. N. N. Lebedev, Special Functions and their Applications (Prentice Hall, 1965).

15. G. B. Arfken and H. J. Weber, Mathematical Method for Physicist (Academic Press, 2001).

16. D. M. Healy Jr., D. Rockmore, P. J. Kostelec, and S. Moore, “FFTs for the 2-sphere - improvements and variations,” J. Fourier. Anal. Appl. 9, No.4, 341–385(1998). [CrossRef]

17. P. N. Swarztrauber, “On computing the points and weights for Gauss-Legendre quadrature,” SIAM. J. Sci. Computing 24. Issue. 3, 945–954(2002). [CrossRef]

18. M. Wieczorek, SHTools, http://shtools.ipgp.fr/.

19. Z. G. Horvath, “Beyond the beam: A history of multidimensional lasers,” Opt. Photonics News 23. No. 7/8, 36–41(2012). [CrossRef]

Computation Method	Time (seconds)

Direct Integration	3730.63
Convolution	0.0573
Spherical spectrum	0.0296

Fast calculation of spherical computer generated hologram using spherical wave spectrum method

Abstract

1. Introduction

2. Theoretical background

3. Numerical computation

4. Simulation results

5. Concluding remarks

References and links

Cited By

Figures (9)

Tables (1)

Equations (38)

Optics Express