Fast calculation method of a CGH for a patch model using a point-based method Y. Ogihara and Y. Sakamoto* Graduate School of Information Science and Technology, Hokkaido University, North 14, West 9, Kita-ku, Sapporo-shi 060-0814, Japan *Corresponding author: [email protected] Received 2 September 2014; revised 27 October 2014; accepted 3 November 2014; posted 20 November 2014 (Doc. ID 222338); published 18 December 2014

Holography is three-dimensional display technology. Computer-generated holograms (CGHs) are created by simulating light propagation on a computer, and they are able to display a virtual object. There are mainly two types of calculation methods of CGHs, a point-based method and the fast Fourier-transform (FFT)-based method. The FFT-based method is based on a patch model, and it is suited to accelerating the calculations as it calculates the light propagation across a patch as a whole. The calculations with the point-based method are characterized by a high degree of parallelism, and it is suited to accelerating graphics processing units (GPUs). The point-based method is not suitable for calculation with the patch model. This paper proposes a fast calculation algorithm for a patch model with the point-based method. The proposed method calculates the line on a patch as a whole regardless of the number of points on the line. When the proposed method is implemented on a GPU, the calculation time of the proposed method is shorter than with the point-based method. © 2014 Optical Society of America OCIS codes: (090.0090) Holography; (090.1760) Computer holography. http://dx.doi.org/10.1364/AO.54.000A76

1. Introduction

Computer-generated holograms (CGHs) simulate the recording process of holography on a computer [1] and are able to display three-dimensional images by using liquid-crystal displays (LCDs). The CGHs generate the virtual scene from the data in a computer and can be used for scientific visualization. The CGHs are calculated by simulating the reflection on the surface of an object, the propagation of the object light, and the interference on the hologram. Recently, graphics processing units (GPUs) have opened up a new phase for the light propagation calculation, which requires a large amount of calculations. Numerous parallel computation algorithms for the calculation have been studied that utilize the massively parallel architecture of GPUs [2–6]. As a point-based method and a fast Fourier1559-128X/15/010A76-08$15.00/0 © 2015 Optical Society of America A76

APPLIED OPTICS / Vol. 54, No. 1 / 1 January 2015

transform (FFT)-based method are suitable for the architecture, they have been actively studied for CGHs. In the point-based method, an object is considered to be covered with point light sources. The calculation of the point-based method is simple, and this method is suited to expression by complex objects. As a result, the point-based method is commonly used to generate CGHs [7–10]. The FFT-based method is suited to expression by patch models assuming an object to be covered by small patches. Patch modeling referred to as polygonal modeling is widely used in 3D computer graphics (CGs), and there are many effective tools and software to design models. Polygonal modeling has been applied to the FFT-based method [11–15], and this makes it possible to use many CG models and modeling tools. The FFT-based method calculates light propagation very rapidly when the number of patches is small because the calculation time is proportional to the number of patches.

In this paper, we propose a fast calculation algorithm for patch modeling with a GPU, which is based on the point-based method. The method here calculates the light propagation of a line on a patch as one whole, with the calculations performed regardless of the number of points on the line, and this accelerates the calculations in patch modeling. The proposed method performs the calculations without recording the positions of the points in the memory. This is very important for GPU calculations, as accessing the memory causes latency delays because the memory on GPUs is limited. The following discusses the algorithm and the acceleration efficiency, and we also report results with optical experiments.

z3 ≫

π x − x2  yh − y2 2max : 4λ h

(3)

In Eq. (3), if the distance between a hologram plane and an object of about 1 cm is much longer than 30 cm, then the Fresnel approximation can be used. When Eq. (3) is assigned to Eq. (1), the Fresnel diffraction can be expressed as ZZ j −jkz ∞ e gx; ypxh − x; yh − xdxdy λz −∞ j  e−jkz fgxh ; yh   pxh ; yh g; (4) λz

uxh ; yh  

where pxh ; yh  is 2. Theory A.

k

2

2

pxh ; yh   e−j2zxh −yh  :

Conventional Method

1. FFT-Based Method Light propagation from an object to the hologram plane is calculated by using the Fresnel–Kirchhoff diffraction integral (Fig. 1). When the light distribution on a hologram is uxh ; yh , the Fresnel–Kirchhoff diffraction integral is expressed by ZZ uxh ; yh  

gx; y expf−jkrgdxdy; r −∞ ∞

(1)

where r is the distance from a point on an object x; y; z to a point xh ; yh  on the hologram, gx; y is the complex amplitude of an object, j is the imaginary unit, and k 2π∕λ is the wavenumber. After this, we call the light wave from the objects “object light.” The distance r is approximately described by r

(5)

The asterisk represents convolution. In Eq. (4), convolution can be calculated rapidly by using the FFT algorithm. Therefore, the calculation of the FFT-based method can be expressed as uxh ; yh  

j −jkz −1 e F fGu; vPu; vg; λz

(6)

where u; v are the coordinates of xh ; yh  in the Fourier domain, F −1 denotes an inverse Fouriertransform operator, and Gu; v and Pu; v are Fourier spectra of gxh ; yh  and pxh ; yh , respectively. In the FFT-based method, an object is created from a set of small patches, as in Fig. 2(a). The FFT algorithm performs the calculations of the light propagation from each patch separately. Calculating the light propagation of patch models is done by the FFT algorithm twice for each patch. However, the FFT-based

q xh − x2  yh − y2  z2

≈z

xh − x2 yh − y2  : 2z 2z

(2)

This approximation is the well-known Fresnel approximation which is widely used in CGHs. The Fresnel approximation can be used to satisfy the conditions [16] such as

Patch model

Hologram plane (a)

Point light sources Hologram plane (b) Fig. 1. Fresnel–Kirchhoff diffraction.

Fig. 2. Models for CGH calculation: (a) patch model; (b) point light model. 1 January 2015 / Vol. 54, No. 1 / APPLIED OPTICS

A77

method requires that the sampling pitch of an object be identical with the pixel pitch of the hologram. The order of the number of calculations with the FFT-based method is OPN 2 log N, where P is the number of patches, and N is the number of pixels on one side of the hologram. 2. Point-Based Method In the point-based method, an object is considered to be covered with point light sources as shown in Fig. 2(b). According to Eq. (1), the light distribution ul xh ; yh  on a hologram from the lth source is described as ul xh ; yh  

gl expf−jkrl g: rl

(7)

The light distribution ul xh ; yh  on a hologram is a summation of the spherical waves from all source. Supposing that the number of point light sources of an object is S, the light distribution Uxh ; yh  from the whole object is calculated by Uxh ; yh  

S X

ul xh ; yh :

(8)

l1

According to Eq. (8), computation volume of the point-based method is OSN 2 . In the conventional point-based method, the calculation time is proportional to the total number of the point light sources and the number of pixels in a hologram. When the point-based method calculates the light propagation of a large object, S becomes large and the calculation times become longer. The size of S also depends on the density of the sources, which influences the appearance of the calculated objects. A higher density of sources shows an observer more continuous surfaces, while a low source density results in spotted surfaces. B.

Proposed Method

Patch modeling has been used in the FFT-based method as described previously. Here, we propose a new faster calculation method for patch modeling using a point-based method, which is suitable for using GPUs. We assume that point lights are arranged at equal intervals Δx; Δz on a patch (Fig. 3). The position of the lth point light xl ; yl  is defined by xl  xo  lΔxzl  zo  lΔz;

(9)

where the position of the 0th point light is x0 ; y0 . By substituting Eq. (9) into Eq. (3), we obtain the following rl ≈ zo  lΔz 

xo − ξ2 lΔxxo − ξ l2 Δx2   : (10) 2zl 2zl zl

Using the approximations described subsequently and sorting in the order of l0 ; l1 ; l2, Eq. (10) is rewritten as A78

APPLIED OPTICS / Vol. 54, No. 1 / 1 January 2015

Fig. 3. Coordination of a hologram and a line on a patch.

  xo − ξ2 Δzxo − ξ2 Δxxo − ξ rl ≈ zo   l Δz −  2zo zo 2z2o 

l2 Δx2 2zl xo − ξ2 l2 Δx2  lΦ  ; 2zo 2zl

(11)

 Δzxo − ξ2 Δxxo − ξ :  Φ  Δz − zo 2z2o

(12)

≈ zo  where



As lΔz ≪ zo the factor 1∕zl in Line 2 of the Eq. (10) can be approximated by   1 1 1 Δz ≈ 1−  l : (13) zl zo  lΔz zo zo The denominator of the third line in Eq. (10) is approximated as z0 because Line 3 of Eq. (10) is smaller than Line 2 of Eq. (10). The use of the different approximations causes small errors which are discussed in the next section. By substituting the approximated r into Eq. (2), we obtain an object light uξ; 0 as n

2

−jk zo xo2z−ξ o

uξ; 0  e

o

L−1 X gl l0

rl

2

2

Δx −jkl 2z

e

l

e−jkΦl :

(14)

Note that the l2 Δx2 ∕2zl  is dependent on l and not on ξ; this means that the term is the phase of the point light source. Here, we assume that the 2 2 gl ∕rl e−jkl Δx ∕2zl  is a complex constant C. This approximation results in some errors and limitation for the arrangement of the point light source, as described in the next section. By using this, the object light is expressed by n o 2 L−1 X −jk zo xo2z−ξ o uξ; 0  Ce e−jkΦl : (15) l0

The summation indicates L terms of a geometric series with a common ratio e−jkΦ . The series is rewritten by the formula for a geometric series as

uξ; 0  Ce−jk



2

zo xo2z−ξ o



1 − e−jkLΦ : 1 − e−jkΦ

(16)

Equation (16) shows that the object light from point light sources arranged linearly is calculated as one independent of the number of points. As a result, Eq. (16) provides very fast calculation of the object light from a large number of points. C.

Triangular Patch

In the patch model, an object is in many cases expressed by triangular patches. This section describes the method of generation of a triangular patch with the proposed method. Figure 4 shows the coordinates of a triangular patch with the proposed method of this subsection; a triangular patch is expressed by integrating with the lines on the patch. An object light uξ; 0 in a one-dimensional hologram from each line is described by Eq. (11). By extending to two dimensions, an object light uξ; η with M lines is defined as uξ; η 

M−1 X

n

2

2

mo −η jk zo xmo −ξ2zy mo

Cm e

m0

o

1 − e−jkLm Φxy ; 1 − e−jkΦxy

(17)

where xmo ; ymo  is the starting point of the mth line, Cm is the constant at the mth line, Δx; Δy; Δz is the interval of point light sources, and Δzfxmo − ξ2  ymo − η2 g 2z2o fΔxxmo − ξ  Δyymo − ηg  : zo

Φxy  Δz −

(18)

Equation (17) shows that the calculations of the light wave of a triangular patch is performed by linecalculations of the number of lines on the patch.

The computation volume of the proposed method is OPLN 2 , where L is the number of lines on one patch, p asa line is calculated as 1. This is similar to OP QN 2 , where the average number of point light sources per patch is L2 . The proposed method will be very much faster than the conventional point-based method when the number of points is large, and the computation volume of the number is larger than that of the FFT-based method with patches. It may be assumed that computation volume of the proposed method is intermediate between the conventional point-based method and the FFT-based method. B. Interval of the Point Light Source 2

3. Theoretical Discussion A.

δx >

Volume of Calculations

The computation volume of the conventional pointbased method and FFT-based method is defined by OSN 2  and OPN 2 log N, respectively, as detailed in the previous section. The computation volume of the point-based method can be rewritten as OPQN 2  using the relation S  PQ, where Q is the average number of point light sources per patch.

λR ; A

(19)

where R is the distance between the point light source and the lens, and A is the diameter of the pupil. When R  50 cm and A  5 mm, then δx is about 0.1 mm. The intervals Δx and Δy of the point light sources must be wider than this δx. Figures 5(a) and 5(b) show the reconstructed images when the point light sources interval satisfies the condition of Eq. (19). The images are high contrast and clear, and we can see the details of the triangles. On the other hand, when the interval in Fig. 5(c) is 0.05 mm, which is too small to satisfy the condition, the reconstructed image is unclear and blurred. C.

Fig. 4. Coordination of a triangular patch.

2

We assumed that gl ∕rl e−jkl Δx ∕2zl  is a complex constant in the previous section. This assumption means that each point light source has a constant amplitude, and that it has the phases −kl2 Δx2 ∕2zl . 1∕rl is almost constant in the Fresnel region, and the constant amplitude will not be a cause of problems because a set of point light sources is on the same patch. But what is the situation when the phase changes? Imagine that there is only a single point light source in a space, and the phase of this point light source changes without any amplitude change. An observer watching the point light source with constant amplitude does not detect any difference in the phase change. If another point light source is closing to the source, the light waves from these two sources will give rise to interference. As an approximation of a line gets finer, Δx approaches 0 and L increases accordingly; the resulting wave starts to resemble a wave converging to a particular point in the space due to the phase term −kl2 Δx2 ∕2zl . As a result, the assumption requires that each point light source must form the isolated image on the viewers’ optical system. The resolution δx of a human eye is defined as

Error

The approximation of the denominator of the Line 3 in Eq. (10), to z0 , causes errors in the calculated hologram data. Line 3 can be rewritten as 1 January 2015 / Vol. 54, No. 1 / APPLIED OPTICS

A79

Fig. 6. Calculation time in terms of the number of points.

A. Comparison of Calculation Times

Fig. 5. Reconstructed images affected by the phase approximation; (a) Δx  0.5 mm; (b) Δx  0.2 mm; (c) Δx  0.05 mm.

  lΔxxo − ξ lΔxxo − ξ Δz 1− ≈ l : zl zo zo

(20)

This makes the distance error Δr Δr 

ΔxΔzxo − ξ 2 l : z2o

(21)

This indicates that the error increases with the distance from the equating position xo · zo . For the error Δr  0.004 mm, when Δx  0.1 mm and a patch size of 10 mm (l  100), xo − ξ  10 mm and zo  50 cm. The error increases in proportion to l2 , but the error is vanishingly small for l < 1000. This condition is satisfied in the Fresnel region. The error causes the degradation of image resolution, but Δr is smaller than that of a human eye. Under the condition of Eq. (19), the errors are not significant for many point light sources because the focused images of the point light sources do not interfere each other. 4. Experiments

We carried out optical experiments and compared the experiments in terms of calculation time to clarity.

We compared the calculation time of the object light by the point-based method, the FFT-based method, and the proposed method. Details of the CPU and the GPU used for the calculation are listed in Table 1. Parallel computation on the GPU was implemented by CUDA version 6.0 [17]. The FFTs of the Fresnel diffraction on the GPU were implemented by the CUFFT library in CUDA. In the FFT-based method, three FFTs were operated using the CUFFT, and the multiplication of complex amplitudes was calculated to generate an object light from one patch. We measured the calculation time with three FFTs and one multiplication. 1. Comparison of the Number of Points Forming a Patch We measured the calculation time for object light to compare mainly the point-based method and the proposed method. The hologram size was 2048 × 2048. The object to be calculated in this experiment is one square, parallel to the hologram plane, and Δx and Δy are 0.2 mm. Figure 6 and Table 2 show the results of this experiment. The calculation time of the point-based method is proportional to the number of points. In comparison, the calculation time of the proposed method increases in proportion to the square root of the number of points. These results show that the theoretical analysis detailed in the previous section is valid. The differences in the calculation time between the point-based method and the proposed method become longer when the number of points expressing the patch increases. The proposed method is about 120 times faster than the point-based method when calculating 16,384 points Table 2.

Table 1.

GPU Parameters

Calculation Time on the GPU, in Terms of the Number of Points

Time (ms)

[GPU] CPU Memory capacity GPU Core clock of GPU Capacity of global memory The number of streaming processors (SP)

A80

Intel Core i7 8 GB Tesla K20 5.2 GHz 5.0 GB 2496

APPLIED OPTICS / Vol. 54, No. 1 / 1 January 2015

Number of Points 1 256 1024 4096 16,384

Point-Based

FFT-Based

Proposed

Fast calculation method of a CGH for a patch model using a point-based method.

Holography is three-dimensional display technology. Computer-generated holograms (CGHs) are created by simulating light propagation on a computer, and...
650KB Sizes 2 Downloads 7 Views