Recurrence formulas for fast creation of synthetic
three-dimensional holograms
Kyoji Matsushima and Masahiro Takai

A method for accelerating the synthesis of computer-generated three-dimensional 共3-D兲 holograms, basedon conventional ray tracing, is proposed. In ray tracing, computers expend almost all of their resourcesin calculating the 3-D distances between each one of the point sources composing an object and a sampling
point on the hologram. We present recurrence formulas that precisely compute the distances and reducethe computation time for synthesizing holograms to one half to one quarter, depending on the processor
type. We demonstrate that a full-parallax hologram with an area of 4800 ⫻ 4800 pixels, synthesized for
a 3-D object containing 966 point sources of light, is computed within 17 min and is optically reconstructed. © 2000 Optical Society of America
OCIS codes: 000.4430, 090.1760, 090.2870.

1. Introduction

Computer-generated holograms or digital holograms
are useful for creating three-dimensional 共3-D兲 images of virtual objects.1,2 Recent electron-beam lithography techniques enable us to fabricate synthetic
holograms with submicrometer structures.3 In suchlarge-scale holograms, however, a computational explosion in the numerical calculation of complex amplitude distributions for 3-D objects causes abottleneck. Most of the computer-generated holograms described in the literature use the discrete
Fourier transform of two-dimensional 共2-D兲 image
planes.4 –7 In this type of computer-generated hologram the required computation capacity is not solarge, since the fast-Fourier-transform algorithm isapplicable to hologram synthesis. However, it is difficult to synthesize true 3-D images with this method.
One significant method for generating the complex
amplitude distributions for 3-D objects is ray tracing.8,9 In this method objects are considered to becomposed of a large number of discrete point sourcesof light, and the complex amplitude of sphericalwaves from each point source is superimposed on thehologram plane. The ray-tracing method is simple

The authors are with the Department of Electrical Engineering,Kansai University, 3-3-35 Yamate-cho, Suita, Osaka 564-8680,
Japan. K. Matsushima’s e-mail address is matsu@kansai-u.ac.jp.
Received 6 December 1999; revised manuscript received 5 September 2000.
0003-6935兾00兾356587-08$15.00兾0© 2000 Optical Society of America

in principle and is potentially the most flexible insynthesizing holograms for true 3-D objects, whichshould be processed by standard rendering techniques such as hidden-surface elimination or surface
shading. However, there is a large computational
problem: Calculating an interference fringe pattern
requires an enormous computation time. Recent
computer technology makes it possible to synthesize
holograms on a desktop computer, but even on astate-of-the-art computer it takes several hours or as
long as a few days to create a full-parallax hologramby ray tracing. Only a horizontal-parallax-only hologram allows for the quasi-real-time synthesis of 3-D
images based on ray tracing.10,11
A typical method for reducing computation time in
ray tracing is to calculate the real-valued fringe intensity instead of the complex-valued amplitude ofthe light wave from objects.9,12 Another method
uses geometric symmetry to avoid redundant calculation of the distance between a point source of theobject and a sampling point of the hologram.13 In
addition, the use of a lookup table remarkably improves the computation time in ray tracing,12 but themain disadvantage of a lookup table is the huge sizeof the table, which increases with the number of sampling points in the object space. In another approach, computer-graphics hardware is used to assist
in calculating the diffraction pattern of point sources.14 In this vein, an attempt was made to fabricate
a special-purpose computer for the fast synthesis of
3-D holograms.15,16 These approaches may be attractive for resolving problems of the real-time synthesis of digital holograms in the future.
10 December 2000 兾 Vol. 39, No. 35 兾 APPLIED OPTICS

6587

Fig. 2. 2-D coordinates.
Fig. 1. Diagram of reference coordinate system.

In ray tracing, the phase of the light wave from apoint source is obtained from the distance between
the point source and the sampling point on the hologram. Most of the computation time is spent calculating this 3-D distance, because the calculation
requires an arithmetic operation on the square root.
Thus an approximation based on a polynomial expansion of the square-root function is sometimes used togenerate a fringe pattern rapidly. However, the precision of the numerical calculation in this approximation significantly decreases with an increase in theangle of incidence of the object light with respect tothe hologram.
In this paper we propose three types of recurrence
formula for fast and accurate computation of the 3-D
distance. This algorithm takes advantage of the factthat the difference between the distances from apoint source to a sampling point and to the next sampling point is much smaller than the 3-D distanceitself.
2. Ray Tracing and Its Computation Sequence

In ray tracing, we treat 3-D virtual objects as if they
were composed of a large number of point sources oflight. As is shown in Fig. 1, suppose that the position of the jth point source of light is given as 共Xjo, Yjo,
Zjo兲. A spherical wave emitted from a point source isincident upon the hologram plane. The complex amplitude of the spherical wave is sampled on the hologram plane at intervals of ␦x in the X direction and atintervals of ␦y in the Y direction. The sampling position on the hologram plane is represented by Xp ⫽p␦x and Yq ⫽ q␦y in the x and the y directions, respectively. Accordingly, the total complex amplitude sampled on the hologram plane is given as13,14
N⫺1

O共Xp, Yq兲 ⫽aj兺 r 共X , Y 兲 exp关ikr 共X , Y 兲兴,
j

j⫽0

j

p

p

q

(1)

r共xn兲 ⫽ 共 xn2 ⫹ d02兲1兾2,
d0 ⫽ 关共Yp ⫺ Y0兲2 ⫹ Z02兴1兾2,(3)

where d0, defined above, is kept constant during ascan in the X direction.

qwhere N and aj denote the total number of point
sources and the amplitude of the jth point source,
respectively. The distance between the jth point
source and the sampling point located at 共Xp, Yq兲 on
the hologram plane is written as
rj 共Xp, Yq兲 ⫽ 关共Xp ⫺ Xjo兲2 ⫹ 共Yq ⫺ Yjo兲2 ⫹ Zjo兴1兾2.

(2)

The object light is simultaneously incident upon all
the sampling points. In actual sequences in compu6588

tation, however, the location at which the amplitudeis calculated moves the hologram plane in the X or
the Y direction. This is similar to the raster of
CRT’s. The point at the coordinates 共Xp, Yq兲 shownin Fig. 1 scans the hologram in the X direction as
suffix p increases, and this scan line then moves in
the Y direction of the hologram plane. Although theexact procedures are dependent on the implementation and may not be identical, this procedure is commonly used in the calculation of the light wave ofobjects by ray tracing.
In ray tracing, since computers spend almost all of
the time computing the 3-D distance represented by
Eq. 共2兲, it is important to reduce this computation.
Therefore we take advantage of the fact that the 3-D
distance at a sampling point hardly changes at thenext sampling point.
To simplify the problem, we introduce 2-D Cartesian coordinates on the plane that contains both apoint source located at 共X0, Y0, Z0兲 and a scan line at
Yq. As is shown in Fig. 1, we define the x axis parallel to the X axis and the d axis perpendicular to thex axis on this plain. The origin of the 2-D coordinates is at the point 共X0, Yq, 0兲 in the 3-D coordinates.
Sampling points on the x axis are defined as xn ⬅ x0 ⫹n␦x 共n ⫽ 0, 1, 2, . . . 兲, and these points are associatedwith the X coordinates through xn ⫽ Xp⫹n ⫺ X0.
Hence the 3-D distance between the point source at
共X0, Y0, Z0兲 and a sampling point at 共Xp, Yq, 0兲 isrewritten as

APPLIED OPTICS 兾 Vol. 39, No. 35 兾 10 December 2000

3. Recurrence Formula

We assume that the incidence angle ␪ formed between the d axis and a straight line from a point
source to a sampling point on the hologram is smallenough to satisfy xn ⬍⬍ d0 共Fig. 2兲, and we apply anapproximation based on a polynomial expansion of
the square-root function in a calculation of r共xn兲.
The approximation, referred to as a binomial approximation in this paper, is obtained by use of the firsttwo terms of the polynomial expansion in the following:r共xn兲 ⯝ d0 ⫹

xn2
.
2d0

(4)

This approximation is suitable for display devices
with low spatial resolution, such as liquid-crystal
panels, because the incidence angle in these devices
is limited to a low value by the sampling theorem toavoid aliasing. In devices not less than approximately 10兾␭ in spatial resolution, the incidence angle
can become so large that the binomial approximation
is no longer appropriate. In that case the numericalerror of approximation 共4兲 would significantly increase with increasing xn. This problem can be
avoided by use of a recurrence formula.
Suppose that the sampling pitch ␦x is sufficiently
smaller than the distance r共xn兲. The distance r共xn⫹1兲
can be expanded into a Taylor’s series around xn:
r共 xn⫹1兲 ⫽ r共 xn兲 ⫹ ␦xr⬘共xn兲 ⫹冉 冊冋 册
冋 册⫽ r共xn兲 ⫹ xn ⫹⫻␦x2
r⬙共xn兲 ⫹ . . .
2

␦x
r共xn兲

␦x
xn2
␦x
⫺
2 r共xn兲2r共xn兲

2

⫹....

(5)

Assuming that ␦x兾r共xn兲 is much smaller than unity,
the first two terms on the right-hand side are a good
approximation for r共xn⫹1兲. Therefore the distanceat xn⫹1 is given approximately by the distance at xn:
r共xn⫹1兲 ⯝ r共xn兲 ⫹共xn ⫹ ␦x兾2兲␦x
.
r共xn兲

(6)

A simultaneous recurrence formula for calculatingr共 xn兲 is obtained by substitution of xn ⫽ x0 ⫹ n␦x into
approximation 共6兲 as follows:rn⫹1 ⫽ rn ⫹ sn兾rn,
sn⫹1 ⫽ sn ⫹ c1.

(7)

Fig. 3. Numerical error of recurrence formula and binomial approximation; d0 ⫽ 105␭.

cal error is defined as Ecal ⫽ 兩rn ⫺ r共xn兲兩, and thiserror is indicated in units of wavelength ␭, wherer共xn兲 denotes an exact distance.
In Fig. 3 the numerical errors of the recurrence
formula and the binomial approximation of Eq. 共4兲
are plotted on a logarithmic scale as functions of xn
with sampling pitches of the hologram from ␭ to 10␭.
The recurrence formula was evaluated step by stepfrom x ⫽ 0 to xmax ⫽ d0 tan ␪max, where ␪max is the
maximum incidence angle, to avoid aliasing. For
in-line-type holograms, this angle is given by ␪max ⫽tan⫺1共␭兾2␦x兲.
As is shown in Fig. 3, the error of the recurrence
formula is dependent on ␦x, unlike in the binomial
approximation. This error is less than ␭兾10 at thexmax with a sampling pitch of ␭, whereas at the same
location of xmax the polynomial approximation gives
rise to an error in excess of 103␭.
Figure 3 shows that the precision of the recurrence
formula is superior to that of the binomial approximation, and its error is dependent on the sampling
pitch as well as on the incidence angle. Therefore in

Here sn is an arithmetical series equivalent to n␦x2,
and it is introduced to eliminate the multiplicationfrom Eqs. 共7兲. The starter value and the constant c1
are given asr0 ⬅ 共d02 ⫹ x02兲1兾2,
s0 ⬅ 共x0 ⫹ ␦x兾2兲␦x,
c1 ⬅ ␦x2.

(8)

In this recurrence formula only two additions anda division are required for each step. Moreover, thenumerical error of Eqs. 共7兲 has significantly less dependence on the incidence angle ␪, as is described in
the next section.
4. Numerical Errors of the Recurrence Formula

Figures 3 and 4 show numerical errors of distances
calculated by the recurrence formula. The numeri-

Fig. 4. Error contours given by numerical calculations of recurrence formula for sampling pitches of 共a兲 ␦x ⱕ ␭兾2 and 共b兲 ␦x ⬎ ␭兾2.
10 December 2000 兾 Vol. 39, No. 35 兾 APPLIED OPTICS

6589

Fig. 4 we plotted contours of the numerical error on
the parameter plane spanned by the sampling
pitches and the incidence angle. Contours for the
sampling pitches that are less or more than ␭兾2 areshown in the separate Figs. of 4共a兲 and 4共b兲, respectively, because the incidence angle is limited to avoid
aliasing in the case of Fig. 4共b兲 but not in Fig. 4共a兲.
As is shown in Fig. 4共a兲, the error does not exceed␭兾10 for incidence angles less than ⬃60°. In Fig.
4共b兲, since the error of the recurrence formula neverexceeds ␭兾10, there is no contour line for an error of␭兾10.
The numerical error of the recurrence formula
given by Eqs. 共7兲 does not exceed ␭兾10 in almost all of
the area indicated in Fig. 4. However, if a precisionof ␭兾10 or higher is required for a certain fringe calculation, it is possible to use the algorithm for error
control that we present in Section 5.
5. Error Estimation and Error Control

Errors for both the recurrence formula and the binomial approximation can be estimated from remainders of the series expansion.
The polynomial expansion of the square root iswritten as
r共xn兲 ⫽ d0 ⫹

xn2
xn 4
⫺
⫹....
2d0 8d03(9)

The binomial approximation error is given as thethird term on the right-hand side of Eq. 共9兲. Whenx0 ⫽ 0, the error expression is rearranged by substitution of xn ⫽ n␦x:
␦x4 4
EBA ⯝ ⫺n.
8d03

1 ␦x2xk2.
2 r共xk兲3(11)

Since the starter value of r0 has no error, the total
error after repetition n can be obtained by integrationof ⑀共k兲 from 0 to n. Therefore the total error is written as
ERF共n兲 ⫽兰n⑀共k兲dk,(12)

␦x关共x0 ⫹ n␦x兲3 ⫺ x03兴,
6r03

(13)

0

⯝⫺

where r共 xk兲 in the denominator of relation 共11兲 isapproximated by the constant distance r0. In particular, the error of the recurrence formula when x0 ⫽0 is given as
ERF ⯝ ⫺
6590

␦x4 3
n.
6d03

The ratio of the binomial approximation error to therecurrence formula error is given as EBA兾ERF ⫽ 3n兾4by relations 共10兲 and 共14兲. When the incidence angle
increases, n becomes quite large. Therefore the recurrence formula has a remarkable advantage over a
simple binomial approximation.
In the recurrence formula, numerical errors areaccumulated for every single step and increase as therepetition times increase. However, it is possible toreset the accumulated error to zero by recalculationof the exact distance before the error exceeds a permissible limit. The error estimation in relation 共14兲
can be used to predict the critical repetition limit.
The critical repetition limit is obtained asnlim ⫽(10)

In the recurrence formula a single-step error ⑀共k兲,
introduced when xk⫹1 is calculated from xk, is given
as the term 关␦x兾r共xk兲兴2 in Eq. 共5兲:
⑀共k兲 ⬵ ⫺Fig. 5. Error of recurrence formula with error control; d0 ⫽ 105␭.
The sampling pitches and permissible errors are 共a兲 ␦x ⫽ 2␭, Elim ⫽10⫺3␭; 共b兲 ␦x ⫽ 3␭, Elim ⫽ 5 ⫻ 10⫺3␭; 共c兲 ␦x ⫽ ␭, Elim ⫽ 10⫺4␭.

(14)

APPLIED OPTICS 兾 Vol. 39, No. 35 兾 10 December 2000

冋冉

x03 ⫺

6r03Elim
␦x

冊

1兾3

册

⫺ x0 ␦x⫺1(15)

when relation 共13兲 is solved for n, where Elim denotes
the permissible error of the distance.
Figure 5 shows the errors in distance calculations
based on the error-control procedure with permissible
errors of, Fig. 5共a兲, 10⫺3␭; 5共b兲, 5 ⫻ 10⫺3␭; and 5共c兲,
10⫺4␭. All error curves are controlled with less than
the given permissible errors. The overhead for computing nlim, of course, becomes large for small Elimand should not be ignored in terms of the computation time.
6. Variations of the Recurrence Formula
A.

Divisionless Recurrence Formula

In most processors division is one of the hardest
arithmetic operations. Division is usually several
times as long as multiplication. Thus it is significant to eliminate division from the recurrence formula, but it is not worthwhile if such elimination
leads to an increase in the error. To eliminate division, we have introduced a variable pn, defined as thereciprocal of distance rn. Moreover, pn⫹1 is calculated by a single step of the Newton–Raphson
method, which is an excellent algorithm for calculating a reciprocal value.
If ␰ is the reciprocal of ␹, a single step of Newton–

Raphson method for improving ␰ from an estimatedvalue ␰0 is described with a function defined as f 共␰兲 ⫽␹ ⫺ ␰⫺1:␰ ⫽ ␰0 ⫺f 共␰0兲
⫽ ␰0共2 ⫺ ␹␰0兲.
f⬘共␰0兲

(16)

Thus the following recursion relation is given by substitution of pn⫹1, rn⫹1, and pn into ␰, ␹, and ␰0, respectively:pn⫹1 ⫽ pn共2 ⫺ rn⫹1pn兲.

(17)

The recurrence formula of Eqs. 共7兲 is rewritten, withpn, asrn⫹1 ⫽ rn ⫹ sn pn,
sn⫹1 ⫽ sn ⫹ c1.

(18)

The starter value of pn is given as p0 ⫽ 1兾r0. Other
constants are the same as that of the original basic
recurrence formula. This modified divisionless procedure involves three multiplications, two additions,and a subtraction.
The error of the divisionless recurrence formula isactually equivalent to that of the basic one, because
the Newton–Raphson method has quadratic convergence, and pn as the estimated value is quite close topn⫹1 in conventional holograms. When ␦x ⫽ ␭ in
Fig. 3, for example, the difference between distances
calculated by basic and divisionless recurrence formulas does not exceed 2 ⫻ 10⫺7␭ at the maximum
incidence angle. Therefore, Eqs. 共18兲 can be treatedin practice as another expression of Eqs. 共7兲.
B.

Scaled–Divisionless Recurrence Formula

Here we employ birefined numerical operations aswell as variables encoded as double-precision
floating-point numbers to investigate the numericalerror of the recurrence formula. This is necessary toretain the accuracy of the calculated distance and to
avoid rounding errors. However, electronic circuitsthat handle double-precision floating-point operations are too large and too complex to be used in
hardware fabricated especially for the real-time synthesis of holograms.15,16 Additionally, current microprocessors tend to support some special
instruction sets for multimedia data processing;these instructions simultaneously process multiple
data encoded in single-precision floating-point operations. These special functions are expected to be
used to synthesize holograms in the future.
The number of significant digits must be reduced to
let us take advantage of single-precision floatingpoint operations. Therefore the following scaled
variable should be introduced into the recurrence formula:qn ⫽ rn兾r0 ⫺ 1.

(19)

Fig. 6. Error of scaled-divisionless recurrence formula. d0 ⫽105␭.

By substituting r0 from Eq. 共19兲 into relation 共18兲, we
simplify the recurrence formula by using qn to obtainqn⫹1 ⫽ qn ⫹ sn pn,
pn⫹1 ⫽ pn共2 ⫺ qn⫹1pn ⫹ pn兲,sn⫹1 ⫽ sn ⫹ c2.

(20)

The starter formulas and a constant are given as
q0 ⫽ 0,p0 ⫽ 1,s0 ⫽ 共x0 ⫹ ␦x兾2兲␦x兾共d02 ⫹ x02兲,

c2 ⬅ ␦x2兾共d02 ⫹ x02兲.

(21)

Figure 6 shows the error curves of Eqs. 共20兲 with a␦x of 10␭. The original basic recurrence formula
gives no significant values for monorefined operations, whereas the modified version yields effective
distances. However, its error curve involves irregularity owing to the rounding errors.
Unfortunately, because of rounding errors, eventhe modified version is not useful when ␦x is less than⬃10␭.
7. Measurement of Computation TimeThe time required for synthesizing holograms is directly proportional to two factors: the number ofpoint sources for the object and the number of pixels
in the hologram. Therefore the computation time ofa hologram with M pixels and for an object with N
point sources is given as
Ttotal ⫽ ␶a MN,(22)

where ␶a is an elemental time coefficient that is dependent on the arithmetic algorithm; ␶a is defined as
the time necessary for calculating the fringe intensityof a single point source in a single pixel. We employed the coefficient ␶a in units of nanoseconds per
pixel per point to evaluate the performance of therecurrence formulas.
Two processors were used in the performance mea10 December 2000 兾 Vol. 39, No. 35 兾 APPLIED OPTICS

6591

surements. One was an Intel Pentium II with aclock frequency of 450 MHz; this is the most common
processor and is representative of complexinstruction-set computers 共CISC’s兲. Pentium-series
processors have an instruction set for calculating thesquare root and execute it in almost the same time as
the division operation. The other processor was an
Alpha 21164A with a clock frequency of 600 MHz.
This is the typical processor for reduced-instructionset computers 共RISC’s兲 and has no instruction set forthe square root. The basic design of the Alpha processor emphasizes the fast execution of simple arithmetic operations instead of the computation of higher
functions such as square roots or trigonometric functions.
In the benchmark program we computed the realvalued fringe intensity given as Re关OR*兴, where R
denotes the distribution of the complex amplitude of
the reference light. Since we assume in-line holograms, R is a constant. Thus the bipolar fringe intensity was calculated as
I共Xp, Yq兲 ⫽ Re关O共Xp, Yq兲兴
⬵

N⫺1

aj

j⫽0

j

兺 Z cos kr 共X , Y 兲,
o

j

p

q

(23)

where we replaced the factor aj 兾rj 共Xp, Yq兲 in Eq. 共1兲
with aj 兾Zjo for simplification.
The distance rj 共Xp, Yq兲 was computed by use of thesquare root and three types of recurrence formula,i.e., basic, divisionless, and scaled– divisionless. Birefined operations were employed in the square root,
the basic, and the divisionless recurrence formulas,and monorefined operations were used in the scaled–divisionless recurrence formula. The benchmark
program for both processors was implemented with
Microsoft Visual C⫹⫹ 6.0. This program uses a table lookup method for the trigonometric function but
does not employ any reduction method based on geometric symmetry.13
We measured the computation time TN that is necessary for calculating the interference pattern for N
point sources in a hologram with 2048 ⫻ 2048 pixels.
The time coefficient was obtained as ␶a ⫽ 共T10 ⫺T5兲兾共20482 ⫻ 5兲. The time coefficient was determined only for calculating the fringe pattern and did
not include overhead time for image scaling, file
input– output, or any initialization.
Figure 7 shows the measured time coefficients.
The longest time coefficient was 294 共ns兾pixel兲兾point
obtained for the square root on the Pentium II,
whereas the shortest was 47 共ns兾pixel兲兾point in thedivisionless recurrence formula on the Alpha 21164A.

Fig. 7. Performance of recurrence formulas.

rections and a dimension of 3 cm ⫻ 3 cm. We used
the wire frame of a cone as an object, which contains
966 point sources of light. The object has an actual
width of 2 cm in the X direction and was set at a Z
position of 20 cm behind the hologram.
Figure 8 shows the optically reconstructed image ofthe fabricated in-line hologram. The total calculation time required for synthesizing the hologram was
1005 s on the Alpha processor. This leads to a timecoefficient of 45 共ns兾pixel兲兾point, which is a little
shorter than the coefficient shown in Fig. 7. This
shorter time coefficient is attributed to our limiting
the incidence angle to avoid aliasing, which meansthat light emitted from a point source would not always cover the whole hologram plane in this setup.
9. Discussion

We used double-precision floating-point variables formost of numerical calculations reported here, be-

8. Construction and Reconstruction of a Hologram

A hologram was synthesized by use of the divisionlessrecurrence formula. The hologram was an in-line
type with binary transmittance. It was fabricatedby use of an image setter with a resolution of 4064
dpi. The hologram was composed of 4800 ⫻ 4800
pixels with sampling pitches of 6.25 ␮m in both di6592

APPLIED OPTICS 兾 Vol. 39, No. 35 兾 10 December 2000

Fig. 8. Photograph of optical reconstruction of an in-line hologram calculated with divisionless recurrence formula. For an object with 966 point sources, calculation takes 16 min 45 s on theAlpha processor.

cause it is expected that the rounding errors may beignored in double precision. In a rough estimation,
the variables for the 3-D distance must be capable ofencoding the minimum distance variation of rmin ⫽r共␦x兲 ⫺ r共0兲 as well as the distance d0. Therefore thenecessary significant figures are estimated aslog10共d0兾rmin兲. This is rewritten as log10共2d02兾␦x2兲
by use of the binomial approximation. For example,
the significant figures necessary for a case in whichd0 ⫽ 105␭ and ␦x ⫽ ␭ are approximately 10.3. In
this case, single-precision floating-point computation,
which has approximately seven significant figures,
does not have enough precision to encode the 3-D
distance. For a distance of d0 ⫽ 105␭, the rounding
errors in double-precision floating-point operation
with 15 significant figures are not significant for ␦x ⬎10⫺2␭.
As is seen from Fig. 7, the performance of the arithmetic algorithm was strongly affected by the architecture of the processors. On the Pentium IIprocessor the shortest ␶a was obtained with thescaled– divisionless recurrence formula, which wasapproximately half the time of the worst case, thesquare root. However, on the Alpha processor the
shortest time was obtained with the divisionless recurrence formula, which was approximately a quarter of the worst time. The arithmetic operation of
the square root in RISC processors is usually implemented as a library provided with the compiler supplier. Although such a library is sophisticated andhighly optimized, many steps of the arithmetic operation are needed to guarantee that the calculated
value is precision enough for most of the numericalanalysis. Therefore the recurrence formula, composed of only a few steps of the arithmetic operation,is more effective in RISC processors than in CISC
processors.
We expected the scaled– divisionless recurrence
formula with monorefined arithmetic operation to
be faster than that of birefined operation only on theAlpha processor, because the Alpha processor has a
special instruction set for single-precision floatingpoint operations. However, the results show thatthe computation time for the scaled– divisionless
formula was approximately 4% longer than that forthe divisionless formula on the Alpha processor.
This probably corresponds to optimization of thecompiler, because the computation time recorded
for the scaled– divisionless recurrence formula issensitive to the setting of the compiler’s optimization switches.
10. Conclusion

We have proposed an arithmetic acceleration algorithm that includes recurrence formulas for synthesizing computer-generated 3-D holograms. In theconventional ray-tracing method it is important to
reduce the computation time for 3-D distances between a point source of the object and a sampling
point of the hologram. The 3-D distance can be calculated with the recurrence formulas in a few simple
arithmetic steps.

Three types of recurrence formula have been presented. The basic recurrence formula requires only
two additions and a division to obtain the distance.
The divisionless recurrence formula, in which the division operation is eliminated by the Newton–
Raphson method, is a modification of the basic
formula. The precision of both recurrence formulas
is sufficient for synthesizing holograms. In addition, the error estimation makes it possible to control
the error. In the scaled– divisionless recurrence formula single-precision floating-point operation is allowed for a sampling pitch greater than 10␭. This
type of recurrence formula is expected to facilitate thefabrication of special hardware for the real-time synthesis of 3-D holograms.
The performance of recurrence formulas, which
was defined as the computation time of the fringe
intensity of a single point source for a single pixel,
was measured on Pentium II 共450 MHz兲 and Alpha21164 共600 MHz兲 processors. The divisionless recurrence formula reduced the time required for synthesizing holograms to approximately a half and aquarter on the Pentium II and the Alpha processors,
respectively.
Finally, we constructed an in-line binary hologram with 4800 ⫻ 4800 pixels for an object containing 966 point sources as a demonstration of theperformance of the divisionless recurrence formula.
The hologram was synthesized on the Alpha processor in 16 min 45 s and was optically reconstructed.
The authors thank Y. Maeda and Y. Yasuda forvaluable discussions of mathematics.
References
1. T. S. Huang, “Digital holography,” Proc. IEEE 59, 1335–1347
共1971兲.
2. G. Tricoles, “Computer generated holograms: an historicalreview,” Appl. Opt. 26, 4351– 4360 共1987兲.
3. N. Yoshikawa, M. Itoh, and T. Yatagai, “Binary computergenerated holograms for security applications from a synthetic
double-exposure method by electron-beam lithography,” Opt.
Lett. 23, 1483–1485 共1998兲.
4. F. Wyrowski and O. Bryndahl, “Iterative Fourier-transform
algorithm applied to computer holography,” J. Opt. Soc. Am. A
5, 1058 –1065 共1988兲.
5. F. Wyrowski, “Iterative quantization of digital amplitude holograms,” Appl. Opt. 28, 3864 –3870 共1989兲.
6. M. A. Seldwitz, J. P. Allebach, and D. W. Sweedney, “Synthesisof digital holograms by direct binary search,” Appl. Opt. 26,
2788 –2798 共1987兲.
7. L. Legeard, P. Réfrégier, and P. Ambs, “Multicriteria optimallyfor iterative encoding of computer-generated holograms,” Appl.
Opt. 36, 7444 –7449 共1997兲.
8. J. P. Waters, “Holographic image synthesis utilizing theoretical methods,” Appl. Phys. Lett. 9, 405– 407 共1966兲.
9. A. D. Stein, Z. Wang, and J. J. S. Leigh, “Computer-generated
holograms: a simplified ray-tracing approach,” Comput.
Phys. 6, 389 –392 共1992兲.
10. P. St.-Hilaire, S. A. Benton, M. Lucente, J. Underkoffler, andH. Yoshikawa, “Electronic display system for computational
holography,” in Practical Holography IV, S. A. Benton, ed.,SPIE Proc. 1212, 174 –182 共1990兲.
10 December 2000 兾 Vol. 39, No. 35 兾 APPLIED OPTICS

6593

11. P. St.-Hilaire, S. A. Benton, M. Lucente, J. Underkoffler, andH. Yoshikawa, “Real-time holographic display: improvements using a multichannel acoustooptic modulator and holographic optical elements,” Practical Holography V, S. A.
Benton, ed., SPIE Proc. 1461, 254 –261 共1991兲.
12. M. Lucente, “Interactive computation of holograms using alook-up table,” J. Electron. Imag. 2, 28 –34 共1993兲.
13. J. L. Juárez-Pérez, A. Olivares-Pérez, and R. Berriel-Valdos,“Nonredundant calculation for creating digital Fresnel holograms,” Appl. Opt. 36, 7437–7443 共1997兲.

6594

APPLIED OPTICS 兾 Vol. 39, No. 35 兾 10 December 2000

14. A. Ritter, J. Böttger, O. Deussen, M. König, and T. Strothotte,“Hardware-based rendering of full-parallax synthetic holograms,” Appl. Opt. 38, 1364 –1369 共1999兲.
15. T. Ito, T. Yabe, M. Okazaki, and M. Yanagi, “Special purpose
computer holography HORN-1 for reconstruction of virtual
image in three dimensions,” Comput. Phys. Commun. 82, 104 –110 共1994兲.
16. T. Ito, H. Eldeib, K. Yoshida, S. Takahashi, T. Yabe, and T.
Kunugi, “Special purpose computer for holography HORN-2,”
Comput. Phys. Commun. 93, 13–20 共1996兲.