OSA's Digital Library

Optics Express

Optics Express

  • Editor: Andrew M. Weiner
  • Vol. 21, Iss. 1 — Jan. 14, 2013
  • pp: 935–948
« Show journal navigation

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

Boaz Jessie Jackin and Toyohiko Yatagai  »View Author Affiliations


Optics Express, Vol. 21, Issue 1, pp. 935-948 (2013)
http://dx.doi.org/10.1364/OE.21.000935


View Full Text Article

Acrobat PDF (3325 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

A fast calculation method for computer generation of spherical holograms in proposed. This method is based on wave propagation defined in spectral domain and in spherical coordinates. The spherical wave spectrum and transfer function were derived from boundary value solutions to the scalar wave equation. It is a spectral propagation formula analogous to angular spectrum formula in cartesian coordinates. A numerical method to evaluate the derived formula is suggested, which uses only N(logN)2 operations for calculations on N sampling points. Simulation results are presented to verify the correctness of the proposed method. A spherical hologram for a spherical object was generated and reconstructed successfully using the proposed method.

© 2013 OSA

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

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

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

] 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

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

]. 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(N4) or O(N3) 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

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]

] developed a fast computation method where the hologram is a tilted plane with respect to the object surface. Sakamoto et al.[5

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).

] 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

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.

]. Yamaguchi et al.[7

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]

] 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

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]

] 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

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]

], 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

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]

].

Fast computation solutions for spherical computer generated hologram emplyoing PSF(convolution method) was proposed by Tachiki et al.[11

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]

]. 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

Fig. 1 Geometry of the problem.

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
2uεμ2ut2=0
(1)
where r is the radius of the spherical surface of interest and θ, ϕ represent the azimuthal and polar angles in the surface. The Laplacian operator ∇2 in spherical co-ordinates is defined as
2=1r2r(r2r)+1r2sinθθ(sinθθ)+1r2sin2θ2ϕ2
(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
1r2r(r2ur)+1r2sinθθ(sinθuθ)+1r2sin2θ2uϕ21c22ut2=0
(3)
The solution of Eq. (3) can be found by separation of variables [14

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

, 15

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

], which can be expressed as shown below
u(r,θ,ϕ,t)=R(r)Θ(θ)Φ(ϕ)T(t)
(4)
The separable variables obey the following four ordinary differential equations.
d2Φdϕ2+m2Φ=0
(5)
1sinθddθ(sinθdΘdθ)+[n(n+1)m2sin2θ]Θ=0
(6)
1r2ddr(r2dRdr)+k2Rn(n+1)r2R=0
(7)
1c2d2Tdt2+k2T=0
(8)
The solution to all the above equations are derived in by Arfken [15

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

]. Only the final results are used in this research work. The solution to azimuthal Eq. (5) is
Φ(ϕ)=Φ1eimϕ+Φ2eimϕ
(9)
where m must be an integer so that there is continuity and periodicity in Φ(ϕ). Φ1 and Φ2 are constants.

The solution to polar Eq. (6) is
Θ(θ)=Θ1Pnm(cosθ)+Θ2Qnm(cosθ)
(10)
where, Pnm and Qnm are the associated Legendre polynomials of first and second kind respectively and Θ1 and Θ2 are constants. Qnm is not finite at the poles where cos(θ) = ±1, so this solution is discarded(Θ2 = 0).

For the radial differential equation Eq. (7) the solutions are
R(r)=R1hn(1)(kr)+R2hn(2)(kr)
(11)
where, hn(1) and hn(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 (R2 = 0)

The solution to Eq.(8) is
T(ω)=T1eiωt+T2eiωt
(12)
Again here we assume T2 = 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 Ynm[14

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

, 15

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

] defined by
Ynm(θ,ϕ)(1)m(2n+1)(nm)!4π(n+m)!Pnmcos(θ)eimϕ
(13)
where the quantity
P¯nm=(2n+1)(nm)!4π(n+m)!Pnm(cosθ)
(14)
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)
Ynm(θ,ϕ)=P¯nm(cosθ)eimϕ
(15)
Combining all the above equations, the travelling wave solution to Eq. (3) can be represented as
u(r,θ,ϕ,ω)=n=0m=nnAmn(ω)hn(kr)Ynm(θ,ϕ)
(16)

The radiated field is completely defined when the coefficient Amn 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 Ynm(θ,ϕ)* and integrating over the sphere, gives
Amn=1hn(ka)π2π202πu(a,θ,ϕ)Ynm(θ,ϕ)*sin(θ)dθdϕ
(17)
where, dΩ = sin(θ)dθdϕ′, is the solid angle for integrating on a sphere. Inserting the expression for Amn back into Eq. (16), we get,
u(r,θ,ϕ)=n=0m=nnYnm(θ,ϕ)([π2π202πu(a,θ,ϕ)Ynm(θ,ϕ)*dΩ]hn(kr)hn(ka))
(18)

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

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

] which is given by
u(x,y,z)=14π2ei(kxx+kyy)dkxdky([u(x,y,0)ei(kxx+kyy)dxdy]eikzz)
(19)
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 (kx, ky). The propagation of the spectrum is defined by the quantity eikzz which is known as the transfer function. The propagated spectrum is recomposed into the wavefield at destination by the integral over kx, ky. 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
    Umn(a)=u(a,θ,ϕ)Ynm(θ,ϕ)*dΩ
    (20)
  • The decomposed wave components(spectra) are expressed by Eq. (13) which is composed of a travelling wave component in ϕ and defined by eimϕ and a standing wave component in θ given by Pnm(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 kx and ky 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 Ynm(θ,ϕ) as the basis function. In other words, the spectral space for spherical system is spanned by Spherical harmonic functions( Ymn(θ,ϕ)). 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
    Umn(r)=hn(kr)hn(ka)Umn(a)
    (21)
  • Hence the quantity hn(kr)/hn(ka) can be refered to the Transfer function(TF) for spherical system, as opposed to the quantity eikzz for planar system.
  • The inverse spherical harmonic transform(ISHT) that recomposes the wave field back from the spectrum is given by
    u(r,θ,ϕ)=n=0m=nnUnm(r)Ynm(θ,ϕ)
    (22)

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.
hn(1)(x)=(i)n+1eixix
(23)
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).

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.

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 Ynm 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
Umn(r)=π/2π/2(ππu(r,θ,ϕ)eimϕdϕ)P¯nm(cosθ)dθ
(26)
First, the quantity within the round brackets alone is to be computed which is nothing but a Fourier transform operation. The Fourier coefficients Um(θ) are evaluated for m = −N,...,N as shown below
Um(θ)=ππu(r,θ,ϕ)eimϕdθ
(27)
=1Ii=1Iu(r,θ,ϕi)eimϕi
(28)
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 Um(θ) is to be evaluated for |m| ≤ nN. This is done using the Gaussian-Legendre quadrature as shown below,
Unm=π/2π/2Um(θj)P¯nm(cosθ)sinθdθ
(29)
=j=|m|NUm(θj)P¯nm(cosθj)wj
(30)
where θj and wj are respectively the Gauss nodes and weights and are calculated using the Fourier-Newton method as described by Swarztrauber [17

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]

]. 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(θ,ϕ)=m=NN(n=|m|NUnmPnm(cosθ))eimϕ
(31)
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
Um(θ)=n=|m|NUnmP¯nm(cosθ)
(32)
u(θ,ϕ)=m=NNUm(θ)eimϕ
(33)
Thus using this numerical procedure fast computation of wave propagation in spherical computer generated holograms is achieved, which uses only O(N(logN)2 operations for the SHT computation.

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

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]

] 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)×TF])+ISHT[SHT(Reference)×TF]|2
Then a hologram for the same object was simulated using the the well known direct integration formula defined in Eq. (34).
H(r,θ,ϕ)=O(θ,ϕ)exp(ikL)Ldxdy
(34)
where,
L=r2+a22ra[sin(θ)sin(θ)+cos(θ)cos(θ)cos(ϕϕ)]
(35)
Hologram=|Hobject(r,θ,ϕ)+Hreference(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

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]

] 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.
Fig. 4 Computed hologram(intensity) using a)proposed method and b)direct integration.

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.
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).

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.
Fig. 8 Hologram(intensity).

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×Conjugate[Reference])×TF]|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))2 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.

Table 1. Comparision of calculation speed

table-icon
View This 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))2 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

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

] 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]

OCIS Codes
(050.1960) Diffraction and gratings : Diffraction theory
(090.1760) Holography : Computer holography
(090.2870) Holography : Holographic display

ToC Category:
Holography

History
Original Manuscript: October 16, 2012
Revised Manuscript: November 19, 2012
Manuscript Accepted: November 22, 2012
Published: January 9, 2013

Citation
Boaz Jessie Jackin and Toyohiko Yatagai, "Fast calculation of spherical computer generated hologram using spherical wave spectrum method," Opt. Express 21, 935-948 (2013)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-21-1-935


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  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,” Science280, 1544–1545 (1998). [CrossRef]
  2. N. M. Lawandy and R. M. Balachandran “Random laser?,” Nature373, 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. A10, 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 ApplicationsTung H. Jeong and Hans I. Bjelkhagen, eds., Proc. SPIE5742, 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.Express13, 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.Express18, 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, 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. Computing24. 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 News23. No. 7/8, 36–41(2012). [CrossRef]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited