PROCEEDINGS OF THE IEEE

• Document Type: CREST
• Collection: General CIA Records
• Document Number (FOIA) /ESDN (CREST): CIA-RDP88B00553R000100280003-6
• Release Decision: RIFPUB
• Original Classification: K
• Document Page Count: 148
• Document Creation Date: December 19, 2016
• Document Release Date: February 14, 2006
• Sequence Number: 3
• Publication Date: January 1, 1981
• Content Type: OPEN

W-ryf9e~~ 3F
~~??9 4# CI 8 F ~ I A
$0003;:
THE I`NSTI7U~
32. D ELECTRONICS ENGINEERS
ISSND018 821'q) a'
JANUARY 1981
MONOMC)DE OPTICAL IBERS
ILy HADIsNG.& REFLECTA,N?CE MAPS
i~:.k -i -r ., ..? ... ;, . ,.. 6~.hnr?.... . bx+^ ,$r .~?. r1y..7 . J/ e
Hy = (e/p) I/2 Ex
while the other field components Ey, Ez, Hx, and Hz are neg-
ligible, >V specifies the spatial variation in the plane perpendic-
ular to the fiber axis, p is the permeability of the medium, e =
eon2, where n nco = nci and co is the dielectric constant of
vacuum.
Because nco = nci , the fields are only weakly influenced by
the polarization properties of the fiber structure [51, [9]. If
this is not obvious, then recall that plane wave reflection from
a dielectric interface is nearly insensitive to the polarization of
the incident wave when the two dielectrics are similar. Accord-
ingly, the spatial dependence, 0, of the fields must be insensi-
tive to polarization effects so that t 11 is a solution of the scalar
wave equation, i.e., 1l ''II,,
{aY2 +---+k2(r) r aY ( W(r)=Q2D(r) (5)
k(r) = 27rn(r)/A. (6)
The solution of (5) corresponding to the fundamental mode is
that with the largest 0 and with ii independent of the polar
angle 0. Although polarization effects are very small, they
must be considered in determining 0 for HE,, modes of a non-
circular fiber (Section IX).
In summary, the fundamental mode is approximately a trans-
verse electromagnetic wave given by (4) with the spatial de-
pendence 0(r) a solution of the scalar wave equation. If our
arguments are difficult to follow, then the reader may prefer
the derivation by formal perturbation methods [51, [91,
[10].
III. GAUSSIAN APPROXIMATION FOR FIELDS
OF THE FUNDAMENTAL MODE
Our main objective is to find a good approximation for the
field i/i(r) and the propagation constant 0 of the fundamental
mode on fibers with refractive index profiles, n(r), like those
of Fig. 1(b) where in is defined in Section VII. For such pro-
files, we know that i(i(r) must be maximum at r = 0, decreasing
to zero with increasing r. Furthermore, numerical solution of
(5) for the step and power-law core profiles show that i1(r) is
approximately Gaussian in appearance [ 11 ] - [ 15 ].Accordingly
we assume that the field of the HEII mode has the form
iji(r) = exp [- z (r/ro )2 ]
where ro, the spot size, is determined by the following ele-
mentary, variational method [ 14] . If (7) is an accurate ap-
proximate solution of (5), then it can be used as a trial function
in a stationary expression for the propagation constant, (3, with
the spot size found as the value of ro giving rise to the largest 0
[141, [161. Recall that the fundamental mode is defined as
that with the largest 0. The stationary expression for R of (5)
is (see Additonal Note on (8) at the end)
fo- {- d ) z +k2(r)02} rdr
dr
R2 =-
J'
rV~2 dr
Thus the procedure for finding the spot size ro is straight-
forward. We simply substitute the approximation of iJi(r)
given by (7) into (8) and then find ro as that value satisfying
(a(32/aro) = 0. The approximation for the propagation con-
stant 0 is found by substituting ro into (8). Knowing ro and a,
the field is fully specified by (4) and (7). It remains only to
specify the refractive index profile n(r) and perform the algebra
discussed above. To this end we find it convenient to express
n(r) as
n2(r) = n.I +s(r2/pz)i2
where A, the refractive index difference parameter, is defined
as
(n 2 2 1/2
co ncl)
while s(r2/p2) characterizes any profile shape (s = 1 at maxi-
mum) with p a scaling parameter. For the profiles of Fig. 1(b),
s = 1 at r = 0 and s = 0 at r = ??. Our definition of A differs
from some authors, e.g., [7], but it is more convenient for
this presentation.
IV. EXAMPLE 1: THE GAUSSIAN REFRACTIVE
INDEX FIBER
We first consider the refractive index profile of Fig. 1(b) for
i*pD Bo 9553]fkoob 1/?bk8b6o3~6i Gaussian func-
1 The fields of higher order modes are also nearly TEM but they are
not uniformly polarized iAppi O IFdr Release 2006/03/10
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
s(r2/P2) = exp [-(r/P)2 ]. (11)
The Gaussian profile is unique because it has the same shape as
the field intensity t/ ,2 given by (7) but with p playing the role
of ro. Furthermore, it approximates those fibers for communi-
cation that undergo significant diffusion at the core-cladding
interface during the manufacturing process. Most importantly
however, the Gaussian profile provides us with a very simple
model for understanding all aspects of propagation on mono-
mode fibers and one where our approximate solution is highly
accurate.
Substituting (7) and (11) into (8) leads to an expression for
Q in terms of ro given by
(PQ)2 = (Pkeo)2 - (P/ro)2 - V2 {(P/ro)2 + I}-' (12)
where koo = 27rnco/X and we have introduced the usual dimen-
sionless fiber parameter V defined as
The spot size ro t:o be substituted into (7) and (12) is found
by solving (a/32/aro) = 0 leading to
ro =p2/(V- 1). (14)
This expression is physically meaningful only when V> 1
(ro positive) but we find below that this does not detract from
(14) providing a comprehensive description of communica-
tion fibers. Substituting ro into (12) leads to an expression for
/3 of the form
(PQ)2 = (Pkco)2 - 2V+ 1 (15)
which is restricted to V 2, since from (3) we know that /3
kc0 .
The spot size and propagation constant fully specify the field
of the fundamental mode and thus the light transmission prop-
erties of monomod.e fibers as we discuss below. The usefulness
of the parameter V is evident from the fact that ro /p, and hence
i,li(r), depend on V only.
Range of V for Single Moded Operation: When V is increased
above a certain value other modes can propagate, e.g., this value
is V=2.41 for the step profile fiber of Fig. 1(b) [4] and
V = 2.59 for the Gaussian profile fiber.2 However, in practice
fibers tend to be effectively single moded for larger values of
V, say V < 3 for the step profile, because the higher order
modes suffer radiation losses due to fiber imperfections [ 17 ] .
From this discussion we can appreciate one reason why
monomode fibers are weakly guiding. To be single moded,
V < 3 for the profiles of Fig. 1(b), which from (13) restricts
the physical dimensions of the fiber cross section to p < 3X/2z.
Supposing that X =- I ?m, then we see that p is minute unless
A 2. Because this is undesir-
able in practice, it explains why we are unconcerned by the
restriction of (14) to V > 1.
B. Profile Dimension for Maximum Light Concentration
It is interesting to determine the refractive index profile
dimension p for which the modal light intensity will be most
densely concentrated near to the fiber axis, assuming a given
A and wavelength of excitation. In other words, we ask what
p gives the smallest spot size ro. By differentiating (14) with
respect to p, recalling from (13) that V is proportional to p,
we find that the optimum p occurs when V = 2, i.e., p = X/7r0.
At V = 2, we also observe that ro = p so that the intensity
distribution, S(r), exactly matches the shape of the refractive
index profile s(r2/p2) = exp [-(r/p)2 1.
This last result can also be appreciated in a qualitative fashion
from elementary physics. Suppose we want to confine the light
intensity to the profile shape exp [-(r/p)2 1. Now we know that a
light beam has an innate spread 66 in ray directions due to diffrac-
tion. Thus, if the Gaussian beam is to maintain its exp [- z (r/p)2 ]
field distribution as it propagates along the fiber, it is necessary
that this 6B diffraction spread just equals the angle Bc necessary
for rays to be bound to the fiber, where both 50 and 0. are
inclined to the fiber axis. It is well known from geometric optics
that light inclined at angles 0 > Alneo, where 0 is given by (10),
is lost by radiation (refraction), so we take 0c = Ilnco. The
spread in ray directions 60 due to diffraction is determined by a
two-dimensional Fourier transform of the field distribution,
where the transform variables are r and 2an0/X, with 61 inclined
to the fiber axis. This leads to exp [-2(rrnpe/X)2] for the angular
(0) dependence of the transform of exp [- 2 (r/p)2 ] . Now we
(arbitrarily) define the characteristic diffraction spread So as
that value of 0 for which the transform falls to e-2, so that 60 =
X/2irnp. Demanding that 80 = 1./nco to prevent radiation loss,
is equivalent from (13) to demanding that V = 2. Obviously it is
only fortuitous that this result is identical to that found above,
since our arguments are qualitative. Nevertheless, the derivation
does convey some insight into the physics of propagation on
monomode fiber as a balance between the spread in ray directions
due to diffraction on one hand and their containment by the re-
fractive index gradient on the other hand.
C. Pulse Spreading
Because information is normally transmitted along optical
fibers by pulses, i.e., by digital envelope modulation, we are
interested in pulse distortion. Now a pulse of finite duration
has a characteristic spread in frequencies as does any practical
s
u
th t th 1 d
o
rce so a e pu se sprea s as It travels due to material
'The "exact" results reported ing,tis pap~r or the srrppyy~~qq~~ e,~ ~y~~~yyt~ ~~{{e
profiles of Fig. 1(b) AkW 6~1d~~te itQXet7~ ~~6jZn6klbl t341~liltibe CIAsR13R88B3 53 W(Ha &28>1tI3d#pends on frequency,
of (5). and also due to waveguide (intramodal) dispersion, because 0
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
SNYDER: UNDERSTANDING MONOMODE OPTICAL FIBERS
TABLE I
PARAMETERS OF THE FUNDAMENTAL MODE
(Parameters of the fundamental mode expressed in terms of the spot size
ro of the modal field as defined by (7) or by 11. 1 of Table II. Explicit
expressions are given in Table II for both ro and /3 for the profiles
considered in this paper. The expressions for V,-' and I) are derived
using II.4 and 11.5 of Table II.)
ted
-di
ower density
1 1 .z -(r/ro)2
H
S(r) =
= 2 E/u) e
rec
z
p
y
2 X
Total power
(r
P = 2nr J rS (r) dr
0
( 1
= 2 luJ ~r o
Jt rS(r)dr
Jr
-(r/r
Fraction
of power
n(r) =
= 1- e o
within 0
to r
f rS(r)dr
0
1 = dm = -"Co. ( lI2
{ 1 +
v
n
P
+ y
S
g
l
V
l l
f l
2
(Group velocity)
o
o
y = (pf3) 2 - (Pkco
)2
z
2l
n3 dv
Wavegvide dispersion
D- 1 Co
(
,
p
1 d p
}
{V ad
parameter
aw
A
put
v
o
(neglecting terms in A)
Smoothed Out"
Gaussian Profile
Step Profile
Step Profile
S(r)
l~el~'
-(r/P)2(V-1)
?~#I'e-(r/p)2RnV2
1(El'e-(r/P)2(m+l)(Vm - 1)
for V>1
e
2~u
2 u
lJ
2 u
~J
2
P
(l i 2
v
l x z 7f E
(
TV
F2/(m+l)(Vm+2 _l)
for V>1
1
ZIUJ
2
uJ knVS
2`uJ
2
)
-(r/P)2(V-1)
1-e-(r/p)29,nV2
- (r/p)2 (m+1)(Vm+2_1)
i~(r
1_e
1-e
for V>l
1- (1/V)2,r=p
2
2 m+
-1
W}1
~2~
((I
A
co 1I
A
LJ1l col+ (
,
}
co D [ m+2 ~~Ir
LUUl- I
1)
(
Vg
I
I
1_ `]131C0
VnCO
V2
neO
for V>1
2
D
1/V3
2(1- in V)/V3
(( l
1
mV3 I- 1mm21 Vm+2 L
l J J
for V>l
11
11
(a) (b)
Fig. 2. (a) The normalized power density versus radial distance where S
and P are defined on Table I. (b) The fraction of power flowing within
the cross section 0 to r versus the dimensionless parameter V of (13)
where n is defined on Table I. Both (a) and (b) are for the HE,, mode
of a Gaussian profile fiber. p is shown in Fig. 1(b).
has a nonlinear dependence on frequency. These two effects entiating II.4 of Table II using 11.5 to simplify. Pulse spread
can combine in various ways depending on the nature of the depends on the spread in the group velocity dvg due to the
material dispersion and normally cancel each other at some spread in frequencies dw associated with the pulse. In prac-
wavelength [ 18 ] , [191. Here we examine the influence of tice, however, it is common to quantify pulse distortion by
waveguide dispersion only. the spread in the mode "transit time" or group delay [201,
The group veloci .,,,, for vg iv Y~ra~g C} }~~ C ul ? Differentiating
the pulse. The expression for ug1 m a ~`e' Ii's found` y i fer- g 1 a d1i' j3f ri`t r lse distortion that
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
TABLE II
FIELD OF THE FUNDAMENTAL MODE
(Outline of the procedure for finding the fields, i.e., spot size ro and propagation constant /3, of the
HE? mode on fibers with a profile shape s(r2/p) . Explicit results are given for the class of
profiles depicted in Fig. 1(b).)
2
iQz
~ro1
Ex-e e z (E/U)/ E
II.1
x
with E
, H
, E
and H
negligible
x
Z
,,
Y
n2(r/p) = nc1+ s(r2/P2),S2 I2.2
}'
a=f 2 -ni
l
11.3
2
(PS) 2 = (Pk
) 2 - P J +v2 e-%Ox ds (x) dx+ s (o) - 1~ II.4
co
ro l dx
where kco = 2Trnco/X and V = 2Trpo/I with x = (r/p) 2 and
co = (p/ro) 2
co = (p/r?) 2 is found by solving
1 = -V2 J xe ?Ox ds(x) dx 11.5
-
_
o
a
X
Gaussian Profile
Step Profile
"Smoothe
d Out" Step Profile
s(r2/pl)
e-(r/p)2
1 r< p
?
P-1(m + 1)
tme_t dt
0 r>p
J
a a= (m+l)r2/p2
(P/r
?
V-1
Zn V2
2
(m+ 1) (V m+2-.k) kn V2+ (knV2 -2)knV2
for V> 1
2
(P1)
-2V+1
(pkco)2
(pkco)2-P.o V2-1
2
z +2
) - (m+2)Um
(pk
f
for V> 1
+ (
co
2
m ->s (pkco) 2- In v2 - 1 - (kn
is proportional to DA3, where A is given by (10) and D is a
dimensionless distortion parameter that depends on V only
and is defined in Table I. The greater I D 103, the greater the
pulse spreading due to. waveguide dispersion. Thus by com-
paring fibers with the same V value, we again see an advantage
of the weakly guiding A p, where s is defined by (9). By following the identical
procedure discussed in Section IV for the Gaussian profile
fiber we find that
ro = p2/ln V2
(PQ)2 = (pkco )2 - In V2_ 1
(17)
r
a = p-1g \p , Vl exp I - 2 (0) 2 R (V - 1)V > I `which leads to the'other parameters shown in Table Iwhere
/ L 3 P V 2 J kco = 2rrnco /X.
Approved For Release 2006/Q3bJ0 ~d I ~t~~Q R AQ Q~ ~~Q also to the step
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
SNYDER: UNDERSTANDING MONOMODE OPTICAL FIBERS
+0.06
+0.04
Fig. 3. The waveguide dispersion parameter D, which is a measure of
the group delay or pulse spread in the absence of material dispersion,
as is defined on Table I. m is shown in Fig. 1(b).
profile. However, the radius p for maximum confinement of
light power is found from ro to occur when V = e0.5 = 1.65,
corresponding to ro = p = 1.65 X/27r0 compared to V = 2 with
ro = p = X/ir. for the Gaussian profile. Accordingly, light of a
given wavelength can be more highly confined by a step pro-
file by the fact that ro for the step is 0.83 that of the Gaussian
refractive index profile, assuming both profiles have the same
A value of (10).
The step refractive index profile and the Gaussian profile
differ significantly in their waveguide dispersion characteristics.
The waveguide distortion parameter D(V) for the step fiber
equals 2(1 - In V)/V3 as shown in Table I. This function is
compared to the Gaussian profile result in Fig. 3 where we
note that the step profile has a V at which the waveguide
dispersion is exactly zero.
VII. EXAMPLE 3: FIBERS WITH DIFFUSE
CORE-CLADDING BOUNDARIES
The procedure we describe in Section III for finding the spot
size is applicable to profiles of arbitrary shape so that, e.g., the
fields of 11.1 in Table II can represent any profile n2(r) as spec-
ified by II.2. Substituting this n2(r) into (8) and integrating
by parts leads to an expression for 0, in terms of the spot size,
as given by II.4. The spot size is found by solving (aa2/aro) = 0
or equivalently from 11.5 of Table II. Thus Table II gives
simple expressions for finding the spot size and 0 of arbitrary
profile shapes. However, only special profiles will lead from
II.5 to an analytical expression for spot size, e.g., the Gaussian
profile fiber with (ds(x)ldx) = e-x, with x = (r2/p2) and the
step profile fiber for which (ds(x)ldx) = -S (x - 1) using the
familiar Dirac delta function notation. Having analyzed both
the step and the Gaussian profiles, it is of interest to study the
characteristics of fibers with "smoothed out" step profiles
intermediate between these two extreme cases. The step and
Gaussian profile turn out to be extreme limits of the general
class of profile shapes illustrated in Fig. 1(b) and defined by
shape functions s? x (x) of the form,
sm(x) _ F(m'+ 1) f tme-t dt
(m+1)x
where x = (r/p)2 and F is the usual Gamma function and
For integer values of m, F(m + 1) = m! and (18) becomes
i /
sm(r2/p2)=exp [-(m+1)(r/P)2] M M (m+1) I r)2i
i=0 P
(19)
We then see that m = 0 is the Gaussian profile while it can also
be shown that m = 00 is the step profile. Each profile has the
same value.
The convenient property of these "smoothed out" step
profiles is that asm(x)/ax is a simple analytical expression as
seen from (18). Substituting this into II.5 of Table II leads to
an analytical expression for the spot size given by
ro = p2/(m + 1) {V2/(m+2)- 1 } (20)
which is valid for the continuum of values m > 0 and is identical
to that of the Gaussian profile when m = 0 and the step pro-
file when m = 00, since (x< -1)/a - In x as a - 0. Substituting
(20) into II.4 of Table II leads to the simple expression for (3
in Table II and the other parameters of Table I. As we men-
tioned in Section IV, the maximum V for single mode operation
is about 2.41 for the step and 2.59 for the Gaussian so that the
greater m, the closer the maximum value of V is to 2.41.
Because of our simple analytical expressions for the field, (3
and the other parameters of Tables I and II we can quantify
the role of profile shape on propagation in monomode fibers
without recourse to numerical methods. We consider wave-
guide distortion, since it is much more sensitive to profile
shape than the other parameters discussed above. From the
expression for waveguide dispersion D given in Table I, we see
that there exists a value of V at which D = 0, occuring when
V={(m+2)/m}(m+2)l2-exp (1+1/m) when m is large.
Thus, the value of V for zero waveguide dispersion decreases
from V=-as m-0 to V= e 1 as m--. In Fig. 3 we com-
pare D for the m = 0, 2 and 00 profiles. It is clear that wave-
guide dispersion is highly sensitive to profile shape [26] for
small values of m, say m < 10.
VIII. APPLICATION TO REFRACTIVE INDEX PROFILES
OF MORE GENERAL SHAPE
While the procedure of this paper can be applied, at least
formally, to profiles of arbitrary shape, we can obtain ana-
lytical expressions for a few special cases only. Accordingly,
it is often useful to have a qualitative description of propaga-
tion. Such a description is suggested from the results of the
smoothed out profiles, discussed in the last section. By com-
paring profiles with identical maximum height A and cross
sectional area 27r f o rs (r) dr, we have separated the dependence
of modal parameters on profile shape from that of profile area.
Our results show that spot size ro is comparatively insensitive
to dramatic changes in profile shape, e.g., to within 4.6 per-
cent, all profiles (m = 0 to 00) have the same ro when V = 3.
Furthermore, to within 3.5 percent, the lowest V for propaga-
tion of the second mode is V = 2.5 for all smoothed out pro-
files. These conclusions hold also for profiles of more general
shape, including power law profiles and profiles with depres-
sions about r = 0. Thus profiles with the same maximum height,
A and same area in cross section 27r fo rs(r) dr taking s = 1 at
maximum, have approximately the same spot size and cutoff
V for the second mode. While it is possible to construct ex-
in rs? 1 dr = p2/2.
rule when operated near cutoff of the second mode. Suppose
we are given a profile sha e s (r) with s = 1 at maximum
~
01 db2$0003-6 1
Approved For Release 2006/03/10: CIA-RDP88B01 553R00
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
and a maximum refractive index difference Do defined by
(2) and (10). Propagation of the HE,, mode on this fiber can
then be modeled by one of the smoothed out step profiles
characterized by A and p. Following our rule, we take A =
Do and equate profile areas so that p2 = 2 fo rso(r) dr. It
is not critical which smoothed out profile we use as far as
ro and the cutoff of the second mode are concerned. However,
the waveguide dispersion parameter D varies significantly for
the smaller m values (m < 10) of Fig. 1(b). Accordingly, we
select a smoothed out profile which most nearly resembles
the given profile. For example, if the given profile is step like
with a dip near r = 0, then we use the step profile results listed
in Tables I and II, V = 2.41 for cutoff of the second mode,
where in all parameters p and A are defined as above. In this
way we have a qualitative description of HE,, mode propaga-
tion for a more general class of profiles.
IX. PHILOSOPHY OF THE METHOD AND ITS EXTENSION
TO NONCIRCULAR CROSS SECTIONS
By recognizing that (7) is an exact solution of (5) for a para-
bolic profile, i.e., for s = I - (r/p)2, we appreciate the essence
of the approximation presented in this paper-namely that the
fields of the parabolic profile fiber can be "fit" to the fields of
arbitrary profile fibers. The optimum fit is found from the
stationary expression for 0. Thus the approach can be general-
ized to fibers of noncircular cross section by fitting the HE,,
fields of the noncircular parabolic profile fiber to the noncircu-
lar fiber of arbitrary refractive index grading. For example,
the HE,, fields of fibers with elliptically deformed cross sec-
tions are approximated by
t(x,Y)=exp{-i[(x/xo)2+(Y/Yo)2]} (21)
where the intensity distribution is now elliptical being charac-
terized by lengths x0 and yo. This is the exact solution of the
scalar wave equation for a parabolic profile s = 1 - (x/px)2 -
(y/py)2 but, with the Cartesian form of (8), can be fit to arbi-
trary profiles s(x2/pl , y2/py). From the procedure of Section
III we then derive two equations, each similar to II.5, for deter-
mining x0 and yo. The propagation of HE,, modes on non-
circular waveguides is like that of plane waves in an anisotropic
crystal in that the Q's of the x- and y-polarized HE,, modes
differ [ 9] . This behavior is not displayed by the scalar theory
presented here but is accounted for by perturbation formulas
once the solution for 0 and O(x, y) of the scalar wave equation
are known [5], [9].
In summary, we have used a powerful conceptual tool for
finding the modes of optical fibers where the solution of the
parabolic fiber is the building block from which the modes
of arbitrary profile fibers are constructed.
X. ACCURACY OF THE APPROXIMATION
To appreciate the accuracy of the results presented here, we
again note that (7) is an exact solution of (5) for the para-
bolic profile s = I - (r/p)2. Thus the more the profile deviates
from 1 - (r1p)2, the greater potential error. Now the Gaussian
profile is significantly closer to 1 - (r/p)2 than is the step pro-
file. Accordingly, the results for the Gaussian profile are sig-
ificantly more accurate than those of the step profile. While
this is true [ 14 ] , we emphasize that the results for the step
profile are themselves highly accurate but just not as accurate
as those of the Gaussian profile. For example, the parameter
i7(p) of Table I for the step profile is less than 1.2 percent in
error within the range of greatest practical interest, i.e., (2 <
V < 3). Although the error in 7l(p) increases as V decreases,
even at V = 1.25 it is only 3.7 percent. The accuracy of a is
even better, e.g., the dimensionless parameter p{kco - (32}1/2
is less than 1 percent in error within the range 2 < V < 3.
Obviously, the error depends on the parameter in question,
e.g., while our approximation for 0 is highly accurate, its
second derivative is susceptible to errors. Even so, the distor-
tion parameter D which is proportional to d 2 (p(3 )/d V 2 ex-
hibits the correct qualitative features but the zero crossing in
Fig. 3 is slightly in error, e.g., it should be V = 3.01 rather
than 2.72 for m = -. On the other hand, the error2 is less
than 9.5 percent for D of the Gaussian profile when V ? 1.5.
Field Far from Fiber Axis: The Gaussian approximation for
0(r) is deficient far from the axis, r = 0, where the exact fields
of profiles with circular symmetry have a Ko (wr) dependence
rather than Gaussian, where K0 is a modified Bessel function,
w = rp2 + (27rnC1 /X)2 11/2 and K0(x) - (7./2x)1/2 e-x for large
x. The far (evanescent) field is usually of minor interest, but it
is necessary to accurately evaluate evanescent coupling phe-
nomena such as cross talk between fibers. We can use the
results of the Gaussian approximation to generate an accurate
expression for the far field. To do this we rearrange (5) by
writing k2(r) = kc1 + {k2(r) - kc1 }, leading to
a2 1 a
{aY2 + (Y) aY - w2 0(r) = - {k2(r) - kcl } 1G(r)
(22)
with w defined as above. We then approximate Vl(r) on the
right-hand side by (7), because k2(r) - kel is significant only
for small values of r when the Gaussian approximation is rea-
sonable. The resulting equation is solved by standard Green's
function methods [16], [27] leading, for example, to
fi(r)~(V2/V+1)Ko{(V- 1)r/p}exp[2(V- 1)2/(V+1)]
for the Gaussian profile (11). When V = 2.5, this result is less
than 4 percent in error for r > 1.5 p, taking modal power to be
the same as in the Gaussian approximation (Table I).
Additional Note on (8): This result follows after multiplying
(5) by rV(r) and using the identity
a2 W a>~ a 1 a) - ~aw)2.
r>.~ 2 + rly ar r
ar ar ar ` ar Then, by integrating from 0 to 00, we obtain (8) plus another
term rlJiao/ar evaluated at r = 0 and oo. This term is zero be-
cause i is finite at r = 0 and decays to infinity exponentially
fast.
ACKNOWLEDGMENT
C. Pask conceived the "smoothed out" step profiles given
by (18) and, in addition to T. Burkitt, C. Hussey, J. Love, and
R. Sammut, provided constructive criticism.
REFERENCES.
M. Boerner and S. Maslowski, "Single-mode transmission systems
for civil telecommunication," Proc. Inst. Elec. Eng., vol. 123,
pp. 627-632, June 1976.
W. A. Gambling and H. Matsumura, "A comparison of single-mode
and multimode fibres for long distance telecommunications," in
Fiber and Integrated Optics, D. B. Ostrowsky, Ed. New York:
Plenum, 1978, pp. 333-344.
A. Hondros and P. Debye, "Electromagnetic waves in dielectric
wires," Ann. Phys., vol. 32, pp. 465-476, May 1910.
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
SNYDER: UNDERSTANDING MONOMODE OPTICAL FIBERS
141 D.Marcuse, Light Transmission Optics. New York: Van Nostrand
Reinhold, 1972, ch. 8.
[51 A. W. Snyder and J. D. Love, Optical Waveguide Theory. London,
England: Chapman and Hall, 1981, ch. 12, 13, and 30.
[61 A. W. Snyder, "Asymptotic expressions for eigenfunctions and
eigenvalues of a dielectric or optical waveguide," IEEE Trans.
Microwave Theory Tech., vol. MTT-17, pp. 1130-1138, Dec.
[71
[81
[91
[101
1969.
D. Gloge, "Weakly guiding fibers," Appl. Optics, vol. 10, pp.
2252-2258, Oct. 1971.
J. A. Arnaud, "Transverse coupling in fiber optics. Part II: Cou-
pling to mode sinks," Bell Syst. Tech. J., vol. 53, pp. 675-696,
Apr. 1974.
A. W. Snyder and W. R. Young, "Modes of optical waveguides,"
J. Opt. Soc. Amer., vol. 68, pp. 297-309, Mar. 1978.
D. L. A. Tjaden, "First order corrections to weak-guidance ap-
proximations in fibre optics theory," Philips J. Res., vol. 33, pp.
103-112, 1/2 1978.
[ill W. A. Gambling and I-I. Matsumura, "Simple characterization fac-
tor for practical single-mode fibres," Electron. Lett., vol. 13, pp.
691-693, Nov. 1977.
[121 D. Marcuse, "Loss analysis of single-mode fiber splices," Bell
Syst. Tech. J., vol. 56, pp. 703-718, May 1977.
[131 -, "Gaussian approximation of the fundamental modes of
graded fibers," J. Opt. Soc. Amer., vol. 68, pp. 103-109, Jan.
1978.
[14] A. W. Snyder and R. A. Sammut, "Fundamental (HE11) modes
of graded optical fibers," J. Opt. Soc. Amer., vol. 69, pp. 1663-
1671, Dec. 1979.
[151 A. W. Snyder, "Acuity of compound eyes: Physical limitations
and design," J. Conip. Physiol., vol. 116, pp. 161-182, Jan.
1977.
[161 J. Mathews and R. L. Walker, Mathematical Methods of Physics.
New York: Benjamin, 1965, pp. 315-323; 261-262.
[171 W. A. Gambling, H. Matsumura, and C. M. Ragdale, "Zero mode
dispersion in single-mode fibres," Electron. Lett., vol. 14, pp.
618-620, Sept. 1978.
[181 F. Kapron, `Maximum information capacity of fibre optic wave-
guides," Electron. Lett., vol. 13, pp. 96-97, Feb. 1977.
[191 W. A. Gambling, H. Matsumura and C. M. Ragdale, "Mode disper-
sion, material dispersion and profile dispersion in graded-index
single-mode fibers,"IEE. MOA,vol. 3, pp. 239-246, Nov. 1979.
[20] S. E. Miller and A. G. Chynoweth, Optical Fiber Communications.
New York, Academic Press, 1979, ch. 4.
[211 W. A. Gambling, H. Matsumura, and C. M. Ragdale, "Joint loss in
single-mode fibers," Electron. Lett., vol. 14, pp. 491-493, July
1978.
[22] K. Peterman, "Fundamental mode microbending loss in graded-
index and W-fibres," Opt. Quantum Electron., vol. 9, pp. 167-
175, Mar. 1977.
[231 C. Pask and R. A. Sammut, "Development in the theory of fibre
optics," Proc..IREE Aust., vol. 40, pp. 89-101, June 1979.
[ 241 A. W. Snyder, I. A. White, and J. D. Mitchell, "Radiation from
bent optical waveguides," Electron. Lett., vol. 11, pp. 275-277,
July 1975.
[251 I. A. White, "Radiation from bends in optical waveguides: The
volume current method," IEE J. Microwave, Opt., Acoust., vol.
3, pp. 186-188, Sept. 1979.
[26] A. W. Snyder and R. A. Sammut, "Dispersion in graded single-
mode fibers," Electron. Lett., vol. 15, pp. 269-271, May 1979.
[271 P. M. Morse and H. Feshbach, Methods of Theoretical Physics.
New York: McGraw-Hill, 1953, p. 804.
[28] R. S. Anderssen, F. de Hoog, and J. D. Love, "Modes of Gaussian
profile optical fibres," Opt. Quantum. Electron. Submitted.
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
Hill Shading and the Reflectance Map
BERTHOLD K. P. HORN
Abstract-Shaded overlays for maps give the user an immediate ap-
preciation for the surface topography since they appeal to an important
visual depth cue. A brief review of the history of manual methods is
followed by a discussion of a number of methods that have been pro-
posed for the automatic generation of shaded overlays. These
techniques are compared using the reflectance map as a common
representation for the dependence of tone or gray level on the orienta-
tion of surface elements.
1. INTRODUCTION
0 F THE SEVERAL ways of depicting surface form on
maps, hill shading has the most immediate appeal and
provides for quick comprehension of the topography. In
this sense, hill shading is complementary to the use of con-
tours, which provide accurate terrain elevations but require
careful scrutiny if one is to ascertain the surface form. Shaded
maps are most important when the interpreter's time is limited,
as in aviation, for users that are not trained cartographers, and
for small scale maps, where contours degenerate into messy
tangles of lines.
Why then do we not see more shaded maps? One reason is
the expense of the present manual methods of production,
which require skilled artists with good insight into cartography.
Working from existing contour maps, ridge and stream lines
extracted from such. maps, and at times aided also by aerial
photography, they wield airbrushes, in what is a slow, tedious,
and imprecise operation. Different individuals called upon to
create such images by manual methods will inevitably produce
different results because of the inherent subjective judgment.
The resulting differences in expression of the terrain character-
istics of the same surface at the same scale provide a particular
problem for a map series, where adjoining sheets should match
in terms of hill-shading symbology. This justifies investigation
of an objective system which makes the treatment of all terrain
forms comparable and repeatable.
Attempts at automation began with the notion that the gray
levels used in the shading should derive from a model of how
light might be reflected from a surface. Ignoring shadowing
and mutual illumination effects, it seems clear that the re-
flected intensity will be a function of the local surface inclina-
tion. The choice of a method for calculating the gray tone
based on the orientation of each surface element has however
been the subject of occasionally bitter controversy for almost
two centuries. Much of the difficulty stems from a lack of a
common representation that would allow comparison of
methods which appear at first glance to be incomparable.
The recently developed reflectance map constitutes such a
common denominator. It is a simple device developed origi-
nally for work in machine vision where one is interested in
Manuscript received December 26, 1979; revised July 22, 1980. This
work was supported in part by the Advanced Research Projects Agency
of the Department of Defense under Office of Naval Research Contract
N00014-75-C-0643.
The author is with the Artificial Intelligence Laboratory, Massachu-
setts Institute of Technology, Cambridge, MA 02139.
Digital
Terrain
Model
Gradient
Estimation
Reflectance
Calculation
R
R
Graphic
Output
Processing
Fig. 1. Block diagram of a system for the generation of relief shading.
The gray-value is calculated by applying the reflectance map to the
gradient estimate obtained by sampling neighboring points in the
digital terrain model.
calculating surface shape from the gray levels in an image.
This is clearly just the inverse of the problem of producing
shaded pictures from a surface model. The reflectance map is
a plot of apparent brightness versus two variables, namely the
slope of the surface element in the west-to-east direction and
the slope in the south-to-north direction. Producing a shaded
overlay for a map then is simply a matter of calculating these
two slopes for each surface element and looking up the
appropriate gray level in the reflectance map (see Fig. 1). This
is a very simple, local computation that can.be carried out
efficiently even on enormous databases. The resulting gray
levels can then be fed to a graphic output device that will pro-
duce a continuous tone or halftone photographic transparency
from the given stream of numbers.
What reflectance map is to be used? Careful comparison of
more than a dozen proposed shading methods shows that some
of the simplest provide a good impression of the shape of the
surface. These experiments also show that the most com-
monly used assumptions about surface reflectance do not lead
to the best results, while simple monotonic functions of the
surface slope in the direction away from the assumed light
source work admirably. What matters is the visual impression,
not theoretical rules [ 1 ] . One goal of this paper is a review of
various hill-shading methods that have been proposed in the
past. Much can be learned from these efforts when they are
evaluated in terms of the corresponding reflectance maps.
H. EARLY HISTORY OF HILL SHADING
Chiaroscuro, the technique of using light and shade in
pictorial representation of three dimensional shapes, has been
used by artists for many centuries. Leonardo da Vinci put it
to good effect in his maps of Toscana, drawn in 1502 and
1503, that contained oblique shaded views of relief forms
illuminated from the left [ 1 ] . Woodcuts of the area. around
Zurich in Switzerland drawn half a century later by Maurer use
shaded sideviews as well. Overhead views using relief shading
appear for the first time in maps of the same area drawn a
century after that by Gygers, but these then gave way to less
desirable forms [ 1 ] .
The choice of the representation for relief forms depend to a
great extent on the available reproduction technology.. Wood-
cuts and engraving methods lend themselves to linear forms,
where brightness of an area in the reproduction is controlled
by the spacing and width of darkened lines. Useful directional,
textural effects can be generated by orienting these line frag-
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
0018-9219/81/0100-0014$00.75 ? 1981 IEEE
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
HORN: HILL SHADING AND THE REFLECTANCE MAP
ments, or hachures (Schraffuren), along lines of steepest had been used previously, unguided by his careful analysis
descent. Crowding of such lines in steep areas may have given [271-[29]. While line-based methods give rise to beautiful,
rise to notions of "steeper implies darker." easy to interpret maps, they cannot show the fine detail of sur-
Lehmann proposed the first rigorous relationships [21, [31 face topography possible with halftones and must be based on
between surface slopes and quantities measurable on the smoothed, generalized information such as contours. These
printed map. In 1799, when his method (B6schungsschraffen) lines also tend to interfere with others used to portray
was published anonymously, the techniques for measuring the planimetric information.
surface accurately at a large enough number of points did not A shaded overlay can also be produced by photographing an
exist. Results of this first method of illustrating shape are in appropriately illuminated scaled model of the surface. If this
some ways analogous to those one might obtain by illuminat- model has a matte or diffusely reflecting surface, a map over-
ing a model of the surface from above, an arrangement that lay of high quality will result provided attention is paid to the
gives rise to images that are difficult to interpret. projection geometry. While this was an approach taken early
Partly as a result of this, an alternate form (Schatten- on [27], it really only became practical in the 1950's with the
schraffen) evolved [4]-[6], in which the line thickness is introduction of milling machines that allow an operator to
varied according to the orientation of the local surface patch carve a model by tracing contours on an existing map
with respect to a light source, usually assumed to be near the [30]-[37]. This is still an expensive, slow process however, in
top left of the map when it is oriented properly for viewing. part because of the manual work required to smooth out the
For maps with North at the top this corresponds to northwest. resulting "terraced" model.
Surface patches sloping downward in that direction are por- The Swiss school of cartography improved on earlier forms
trayed with a light tone, while those sloping upward in that [281-[30], [381, [391 and developed shading to a fine art,
direction get a dark tone. Since flat areas have no lines of producing numerous outstanding maps in this time [401-[48].
descent, they remain white. Aside from this defect, this Imhof argues that automated methods, such as relief model
method produces an image similar to one obtained by photography, cannot produce results nearly as impressive,
obliquely illuminating a diffusely reflecting model of the sur- since the cartographer cannot easily influence the process [ 1 ] .
face. Having flat areas appear white makes maps produced by The manual shading method is however slow and expensive,
this method a little difficult to interpret. They are neverthe- and consequently has not been used except for small areas and
less superior to those made by the earlier method, as evidenced those of particular interest or military importance. One can-
for example by the "Dufourkarte" of Switzerland made be- not expect, with significant areas of the world still not mapped
tween 1842 and 1864 using this approach [ 11. These methods at large scales, and the rising cost of labor, that shaded overlays
for portraying surface shape preceded the widespread use of produced this way will be used in many maps.
contours [7], in part because the latter require detailed surface Yoeli [49]-[57] saw the potential of the digital computer in
measurements that were not available before the advent of dealing with this dilemma. It is possible to implement Wiechel's
photogrammetry. method based on oblique illumination of a diffusely reflecting
While lithography was invented by Alois Sneefelder in 1796, surface if terrain elevations can be read into a computer and
it found little application in cartography until around 1850. It suitable continuous tone output devices are available. Yoeli
permitted the production of multicolored maps, but more im- was hampered by the lack of such devices at that time.
portantly, led to the use of halftones, destined to ultimately Blachut and Marsik tried to simplify the required calculations
replace lines as a means of modulating the average reflectance to the point where a computer might not even be required
in the printed map. W. H. Fox Talbot invented a photome- [58], [591. Peucker helped popularize the whole idea of
chanical halftone process in 1852, but commercial success computer-based cartography [191, [60]-[62], and found a
came only years after the patenting of the halftone screen by piecewise linear approximation to the equation for the bright-
Frederick von Egloffstein in 1865, and the crossline screen of ness of a diffuse reflector that works well [ 61 ] . Many other
William A. Leggo in 1869. interesting reports appeared during this time on the subject of
Having access to these new reproduction schemes, Wiechel hill-shading, too numerous to mention individually [63]-[701.
[8] developed shading methods (Schraglichtschummerung) to Brassel [71]-[74] took Imhof's admonitions seriously and
replace the use of hachures as described above. His funda- tried to implement as much as had been formalized of the
mental paper, based in part on work by Burmester [91 on "Swiss manner." With the output devices available to him at
shaded pictures of regular surfaces, placed the field of hill that time it was not easy to judge whether the added com-
shading on a sound foundation. Wiechel discovered the error plexity was worth the effort. All of these computer based
regarding flat surfaces, for example, and developed a graphic methods require detailed digital terrain models. The storage
method for determining the gray value from contour interval capacity and techniques for handling this kind of information
and direction. Unfortunately, the means for controlled genera- now exist [32], [75]-[81 1 as do the photographic output de-
tion of halftones as a function of surface orientation did not vices needed. There has been significant progress, too, in the
then exist and his work was ignored for a long time. automatic generation of digital terrain models directly from
aerial photographs [811-[871, partly as a byproduct of work
III. HILL SHADING IN THIS CENTURY on orthophoto generation [881-[91 ]. More compact and ap-
Two methods based on lines, this time contours instead of propriate representations for these terrain models are under
lines of steepest descent, were explored by Kitiro Tanaka in investigation [921-[95], as are alternate methods for relief
the 1930's. His first method used the lines of intersection of portrayal such as block diagrams [961-[102].
the terrain surface with uniformly spaced, parallel, inclined Considerable progress has been made recently in the com-
planes [ 101, [ I 1 ] . Tanaka's initiative gave rise to considerable puter graphics area in the portrayal of regular objects with
discussion [121-[19], partly in the form of an acrimonious simple surfaces [103]-[116]. Early models for the reflection
debate [201-[231 . Q,His other, Who Qw~s base al Q li t e j 7,.1Q~120 are being elaborated,
of a terraced model 2Sf~i~~epf1'aihQFl`['g~dnw~lt Ci'fi I~owWail MX1~~1t1'Q, Anar surface [ 1211
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
[ 130] . In this context, work on models of the microstructure
of surfaces is relevant [1311-[1361. In a recent effort in the
machine vision area, a method was developed for portraying
the dependence of brightness on surface orientation using the
so-called reflectance map [ 1371-[ 1411. The reflectance map
can be determined if the detailed geometric dependence of re-
flection from the surface [ 1421, [ 143 ] and the distribution of
light sources are known. Alternatively, it can be found
empirically, or derived directly by analyzing the interaction of
light rays with the surface microstructure.
As a result of the development of the reflectance map, the
availability of detailed digital terrain data, small computers
able to perform the simple calculations required, and geo-
metrically accurate gray-level output devices, we may say that
automatic hill shading has come of age-
IV. DIGITAL TERRAIN MODELS
For many applications of cartographic data it is useful to
have machine-readable surface representations. Such terrain
models are used for example in the design of roads and in
order to determine the region irradiated by a radio frequency
antenna. Initially, digital terrain models were generated
manually by interpolation from existing contour maps. This is
a tedious error-prone process producing a digitized version of
the surface represented by the contours, which in turn is a
smoothed, generalized version of the real surface.
The contour information on topographic maps is produced
by manual scanning of stereo pairs of aerial photographs.
Today, fortunately, stereocomparators often come equipped
with coordinate readouts that allow the extraction of informa-
tion needed for the generation of digital terrain models [ 144].
Conveniently taken during orthophoto generation [ 881-[ 911,
the data tends to be accurate and detailed. Even more exciting
is the prospect for machines that achieve stereo fusion without
human help [8l]-[871, since they will lead to the automatic
production of digital terrain models. In the past such ma-
chines had difficulties dealing with uniform surfaces such as
lakes, featureless surfaces, large slopes, and depth disconti-
nuities, as well as broken surfaces, such as forest canopies.
This is apparently still true when aerial photographs are used
with disparities large enough to ensure high accuracy.
Various representations can be chosen for the surface eleva-
tion information. Series expansion, a weighted sum of mathe-
matical functions such as polynomials, Gaussian hills or
periodic functions may be used. These tend to be expensive
to evaluate however and not accurate in approximating sur-
faces that have slope discontinuities. This is important for
many types of terrain, at all but the largest scales. Perhaps the
simplest surface representation is an array of elevations {ztj }
based on a fixed grid, usually square. Determining the height
at a particular point is simple and the interchange of terrain
models between users is easy since the format is so trivial.
One disadvantage of this kind of surface representation is the
high redundancy in areas where the surface is relatively
smooth. The illustrations in this paper are based on digital
terrain models consisting of arrays of elevation values.
Methods that achieve considerable data compression by
covering the surface with panels stretched between specially
chosen points have been developed [ 921-[ 95 1. These exploit
the fact that real geographical surfaces are not arbitrary sets
of elevations but have definite structure and regularity. Such
representations mayApo Xeld FrolrlRelo sei2 /03M'0
voluminous ones, if users can be persuaded to accept the
greater programming complexities involved.
Digital terrain models may also be referred to as digital
elevation models if they contain no information other than
the elevation values.
V. THE REFLECTANCE MAP
The human visual system has a remarkable ability to deter-
mine the distance to objects viewed, as well as their shape,
using a variety of depth cues. One such cue is shading, the de-
pendence of apparent brightness of a surface element on its
orientation with respect to the light source(s) and the viewer.
Without this particular depth cue we would be hard pressed to
interpret pictures of smooth, opaque objects such as people,
since other cues like stereo disparity and motion parallax are
absent in a flat, still photograph. It can be shown that
shading contains enough information to allow the observer to
recover the shape. In fact, a computer program has been de-
veloped that can do this using a single digitized image [ 137 1.
Such work in the area of machine vision has led to a. need to
model the image-forming process more carefully 11381. The
input to the visual sensing system is image irradiance, which is
proportional to scene radiance (here loosely called apparent
brightness) [ 140]. Scene radiance in turn can be related to the
underlying geometric dependence of reflectance of the surface
material and the distribution of light sources [1421, [1431.
Here we concentrate on the dependence of scene radiance on
the orientation of the surface element. Shaded overlays for
maps are interpreted by the viewer using the same mechanism
normally employed to determine the shape of three-dimensional
surfaces from the shading found in their images. Thus shaded
overlays should be produced in a way that emulates the
image-forming process, one in which brightness depends on
surface orientation. This is why the reflectance map, which
captures this dependence, is useful in this endeavor.
Consider a surface z(x, y) viewed from a great distance above
(see Fig. 2). Let the x-axis point to the east, the y-axis north,
and the z-axis straight up. The orientation of a surface element
can be specified simply by giving its slope p in the x (west-to-
east) direction and its slope q in the y (south-to-north) direc-
tion. The slopes p and q are the components of the gradient
vector, (p, q). The apparent brightness of a surface element
R(p, q) depends on its orientation, or equivalently, the local
gradient. It is convenient to illustrate this dependence by
plotting contours of constant apparent brightness on a graph
with axes p and q. This reflectance map [ 138] provides a
graphic illustration of the dependence of apparent brightness
on surface orientation. The pq-plane, in which the reflectance
map is drawn, is called the gradient space, because each point
in it corresponds to a particular gradient.
Surface orientation has two degrees of freedom. We have
chosen here to specify the orientation of a surface element by
the two components of the gradient. Another useful way of
specifying surface orientation is to find the intersection of the
surface normal with the unit sphere. Each point on the sur-
face of this Gaussian sphere again corresponds uniquely to a
particular surface orientation. If the terrain is single-valued,
with no overhangs, all surface normals will point more or less
upwards and pierce the Gaussian sphere in a hemisphere lying
above an equator corresponding to the horizontal plane.
Gradient space happens to be the projection of this hemisphere
from the center of the sphere onto a plane tangent at the upper
C1A-RDP88B00553R000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
HORN: HILL SHADING AND THE REFLECTANCE MAP
x (East)
Fig. 2. Coordinate system and viewing geometry. The viewer is actu-
ally at a great distance above the terrain so that the projection is
orthographic.
While we will not use this representation in the calculation of
relief shading, it is helpful in understanding previous attempts
at. graphical portrayal of the dependence of apparent brightness
on surface orientation. The first such method was developed
by Wiechel more than a century ago [8]. His brilliant analysis
appears to have been largely ignored partly because it de-
pended on mathematical. manipulations that may have been
inaccessible to many of the intended users. Later, Kitiro
Tanaka invented another method showing the variation of
apparent brightness with. surface gradient [101, [ 11 ] , [241,
[251. This second precursor of the reflectance map also ap-
pears to have found little following.
Since the reflectance map gives apparent brightness as a
function of local surface gradient only, it does not take into
account effects dependent on the position of the surface
element. One such effect is illumination of one surface ele-
ment by another. Fortunately this mutual illumination effect
is small unless surface reflectance is quite high [ 138] . It is not
known whether mutual illumination effects aid or hinder the
perception of surface shape. They are difficult to calculate
and so have not been emulated in work on hill-shading.
Another position dependent effect on apparent brightness is
the blocking of light by one portion of the surface before it
reaches another. Cast shadows can be calculated by determin-
ing which surface elements are not visible from the point of
view of the light source [ 1391. Shadows cast by one compli-
cated shape on another are hard to interpret however and ap-
parently detract from the visual quality of shaded overlays
[ 1 ] , [351, [361. They are thus rarely included.
Scattering of light by air molecules and aerosol particles
changes the apparent brightness of a surface element viewed
through the atmosphere. The brightness is shifted towards a
background value equal to the brightness of an infinitely thick
layer of air. The difference between the brightness and the
background value decreases with the thickness of the gaseous
layer through which the surface is viewed [ 145 ]. The resulting
reduction in contrastAPP;8z it'P.c' bWsi6r2VN ft/lb
aerial perspective and can be a useful depth cue, although there
is no general agreement that it aids the perception of surface
shape. It has been used at times by map-makers and can be
modeled easily [1], [71], [73], [74]. The effect has not been
added to any of the hill-shading schemes presented here in
order to simplify comparisons.
VII. WHERE DO REFLECTANCE MAPS COME FROM?
A reflectance map may be based on experimental data. One
can mount a sample of the surface in question on a goniometer
stage and measure its apparent brightness from a fixed view-
point under fixed lighting conditions while varying its orienta-
tion. Instead, one can take a picture of a test object of known
shape and calculate the orientation of the corresponding sur-
face element for each point in the image. The reflectance map
is then obtained by reading off the measured brightness there.
Alternatively, one may use even more detailed information
about light reflection from the surface. The bidirectional re-
flectance distribution function (BRDF) describes how bright a
surface will appear viewed from one specified direction when
illuminated from another specified direction [142], [143].
By integrating over the given light source distribution one can
calculate the reflectance map from this information [1401.
Crudely speaking, the reflectance map is like a "convolution"
of the BRDF and the source-radiance distribution.
Most commonly, reflectance maps are based on phenomeno-
logical models, rather than physical reality. The so-called
Lambertian surface, or perfect diffuser, for example, has the
property that it appears equally bright from all viewing direc-
tions. It also reflects all light, absorbing none. It turns out
that these two constraints are sufficient to determine uniquely
the BRDF of such a surface, and from it the reflectance map,
provided the positions of the light sources are also given.
Some reflectance maps are based on mathematical models of
the interaction of light with the surface. Such models tend to
be either too complex to allow analytic solution or too simple
to represent real surfaces effectively. Nevertheless some have
come quite close to predicting the observed behavior of
particular surfaces [l34]-[136].
Here, new reflectance maps will be determined, based on
proposed methods for producing shaded overlays for maps.
Their derivation will not depend on an understanding of the
image-formation process or the physics of light reflection.
Instead, they will require an analysis of how the brightness of
a point in the overlay depends on the gradient of the underly-
ing geographical surface.
Which reflectance map should be used? The answer to this
question must depend on the quality of the impression a
viewer gets of the shape of the surface portrayed. Various
methods for producing shaded overlays can be compared by
evaluating sample products and classified according to the
corresponding reflectance maps. It will become apparent that
in this way general conclusions can be drawn about a new
method just by inspecting its reflectance map.
VIII. NORMALIZATION OF GRAY TONE
A picture made by applying varying amounts of light absorb-
ing substances, such as ink, to an opaque, diffusely reflecting
material like paper, has a limited dynamic range. Reflectance
is limited at the low end by the properties of the ink and at
the high end by the paper, which will at most reflect all the
light incident upon it, unless it fluoresces. The diffuse reflec-
td1K f '$`8~Y~65~ 618bjg6d jL9 bone. Similarly, if 'UU U12-
 Approved For Release 2006/03/10 : CIA-RDP88BOO553ROO0100280003-6
absorbing substances are used on a transparent substrate, a
limit applies, since transparency cannot be larger than one.
The problem of fitting a given image into the available dy-
namic range is fundamental to all methods of reproduction. A
normalization is applied so that the maximum apparent bright-
ness to be reproduced is represented by a reflectance of one
(or whatever the maximum is for the paper being used). This
scaling will have to be applied whenever relief shading is based
on models of image-formation by light reflected from the
terrain surface.
IX. GRADIENT ESTIMATION
The apparent brightness of a surface element depends on its
orientation with respect to the viewer and the light source.
The orientation of the surface element is described fully by a
surface normal, or equivalently by the gradient. The compo-
nents of the gradient are the slopes p (in the west-to-east
direction) and q (in the south-to-north direction). These
slopes have to be estimated from the array of terrain eleva-
tions. It is convenient to use a short hand here for elevations
in the neighborhood of a particular point (see Fig. 3)..In the
context of a single point at discrete coordinate (i, j), we will
denote the elevation at that point by zoo, while elevations of
the adjacent grid points to the west and east will be called z_0
and z+o , respectively. Similarly, elevations at the points to the
south and north will be denoted zo_ and zo+?
The simplest estimates for the slope p might be
P. = (z+o - zoo )/Ax
p- = (zoo - z_o )/Ox
where Ax is the grid interval in the west-to-east direction, ex-
pressed in the same units as the terrain elevations. These esti-
mates are biased, actually estimating the slope half a grid
interval to the right and left of the central point, respectively.
Their average however, the central difference, is unbiased,
pc = (z+o - z-o )/2 Ax.
Numerical analysis [146]-[149] teaches us that for certain
classes of surfaces an even better estimate is obtained using a
weighted average of three such central differences,
pw = [(z+++2z+o +z+_)- (z_++2z_o +z__)]/80x.
Symmetrically, one can estimate the south-to-north slope,
qw = [(z++ + 2zo+ + z-+) - (z+_+2zo-+z--)]/80y.
These expressions produce excellent estimates for the com-
ponents of the gradient of the central point. The results de-
pend on elevations in a 3 X 3 neighborhood, with individual
elevation values weighted less than they are in the simpler ex-
pression for the central difference. This has the advantage that
local errors in terrain elevation tend not to contribute as
heavily to error in slope. At the same time, more calculations
are required and three rows of the digital terrain model have to
be available at one time.
Care has to be taken to avoid corruption of the slope
estimates by quantization noise in the elevation values. Nu-
merical problems due to the division of small integers may
result when a terrain model is too finely interpolated, with
limited vertical resolution. If it is necessary to generate many
pixels in the output, it is better to interpolate the gray values
produced
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
Z_+
Zo+
Z++
Z_o
Zoo
Z+o
Z-_
Z0_
Z+_
X. GRADIENT SMOOTHING EFFECTS
More complicated slope estimators than the ones described
tend to introduce a smoothing effect, as can be seen by apply-
ing them near points of discontinuity in slope. To illustrate
this more clearly, consider two horizontal smoothing opera-
tions H+ and H- that modify the terrain model as follows:
I
H+: zoo = (zoo +z+o)/2
H-: zoo = (z-o +zoo)/2.
It can now be seen that the central difference slope estimate
p, on the original terrain model, equals the biased estimate p+
calculated from the terrain model smoothed using H-, or,
equivalently, the biased estimate p_ calculated from the terrain
model smoothed using H+. Next consider two vertical smooth-
ing operation V+ and V- in which the terrain model is modi-
fied as follows:
V+: zoo = (zoo + zo+)/2
V-: zoo = (zo - + zoo )/2.
The complicated slope estimate pw can be shown to produce
the same result as the first difference p+ calculated from a
terrain model smoothed by applying H-, V+, and V-. Simi-
larly the slope estimate qw equals q+ calculated from a terrain
model smoothed by applying V-, H+, and H- (actually, since
all of these operations are linear, their order can be arbitrarily
rearranged). Perhaps any "smoothing" desired should be done
as a separate editing operation, combined with the removal of
"glitches" from the digital elevation model, rather than as part
of the slope estimation. Also for terrain models of relatively
limited size this smoothing may be undesirable. Some other
slope estimators are simpler and introduce less smoothing. For
example one can combine two biased estimates of the slope to
get,
PI/2 = [(z+++z+o)- (zo++zoo)]/20x
and symmetrically,
g112 = [(z+++zo+)- (z+o +zoo)]/20y.
Here the average gradient in the top-right quadrant (zoo, z+o,
z++, zo+) rather than at the central point is being estimated,
using elevations in a 2 X 2 neighborhood only. For the graphic
illustrations presented here, the expressions for p1/2 and q1/2
were used to estimate the gradient.
At this time some terrain models are still produced by hand
and have rather limited size. Rather than smoothing the ter-
rain, one may wish to increase apparent resolution by some
means. This can be done quite effectively by combining biased
slope estimates (see Fig. 4). For every point in the terrain
by the shA '3 1F1c4r Release 2006/03/10 : C41R 188B 053 ROO1 1 QII2 A responding to the
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
HORN: HILL SHADING AND THE REFLECTANCE MAP
(p-, q+) (p+, q+)
(p_,q_) (p+,q-)
O 6
Fig. 4. Combinations of biased slope estimates can be used to plot
four times as many gray-tones as there are elevations values in the ter-
rain model. The limited amount of data in a small terrain model
may be stretched this way to produce reasonably detailed hill-shading
output.
four quadrants around it. Each is based on a different com-
bination of the slope estimates (p- or p+) and (q_ or q+) as
appropriate for that quadrant. No miracles should be antici-
pated; this method cannot create information where there is
none, but it can stretch what is available to its limits.
More complicated slope estimators than those discussed here
do not seem called for, since the simple ones shown produce
excellent results. Furthermore, estimators having wider sup-
port, while known to be more accurate for certain classes of
functions such as polynomials, may perform worse on typical
terrain with its discontinuities in slope along ridge and stream
lines.
It has been cartographic practice to assume a light source in
the northwest at a 45? elevation above the horizon. It is help-
ful in this case to introduce a rotated coordinate system as
described in Appendix A.
XI. EXAGGERATION OF TERRAIN ELEVATIONS
Compared to objects of a size that allow for easy manipula-
tion by a human observer, the surface of the earth is in many
places, though not everywhere, rather flat. The range of slopes
is often so small as to cause disappointment with correctly
proportioned models, so that height is often exaggerated in
physical models. Similarly, shading based on models of light
reflection from a surface tends to have undesirably low con-
trast. Here too terrain elevations may be exaggerated for all
but the most mountainous regions. This is equivalent to
multiplication of the components of the gradient by a constant
factor, and corresponds to a simple transformation of the re-
flectance map. For reflectance maps based on reflection of
light originating from an assumed source, a similar effect can
often be achieved by a decrease in the elevation of the source.
For flat surfaces the source may be lowered to a mere 10? or
20? above the horizon, where normally it might be at 45?.
XII. PRODUCING SHADED OVERLAYS
The generation of shaded images from a digital terrain model
using the reflectance map is straightforward (see Fig. 1). For
each point in the terrain model the local gradient (p, q) is
found. The reflectance :map then provides the appropriate
brightness R(p, q), to be plotted on a suitable gray-level out-
put device. All computations are local and can be accom-
plished in a single pass through the image.
To illustrate these ideas a simple program
procedure SHADING(N, M, DX, DY); integer N, M; real DX, DY;
begin array Z[0:N,0:M1, B[0:N-1,0:M-1];
real procedure PE(I, J); integer I, J;
PE :_ (Z[I,.J) + Z[I-I,J] - Z[I,J-1] - Z[I-1,J-1]) / (2.0 X DX);
real procedure QE(I, J); integer 1, J;
QE := (Z[I,J] + Z[I,J-1] - Z[1-1,J] - Z[I-1,J-11) / (2.0 X DY);
real procedure R(P, Q); real P, Q;
RM:=MAX(0.0,MIN(1.0,(1.0-FP-Q)/2.0));

for J := 1 step 1 until M-1 do
for I := 1 step 1 until N-1 do
B[I-1,J-1] := RM(PE(I, J), QE(I, J));

end
Fig. 5. Simple program to generate shaded output from aterrain model.
described later on. Two arrays are used, Z to store the terrain
elevations and B to store the calculated brightness values. The
latter has one row and one column fewer, since its entries
correspond to points lying between those in the elevation array
(the formulas for p1/2 and q1/2 are used). The spacing of the
underlying grid is DX in the west-to-east direction and DY in
the south-to-north direction. The procedures PE(I, J) and
QE(I, J) estimate the slopes, while the procedure RM(P, Q)
calculates the corresponding brightness using a particularly
simple reflectance map. The resulting values range from
0.0 (black) to 1.0 (white) and have to be scaled appropriately.
before they can be fed to a particular gray-level output device.
Typical terrain models are quite large and may exceed allow-
able array storage limits or even the address space of a com-
puter. Fortunately only two (or three) rows of the terrain
model are needed for the estimation of the local slopes. The
program given can be easily modified to read the terrain
model, and to write the calculated gray values, one line at a
time. This makes it possible to deal with terrain models of
essentially arbitrary size.
Next one should note that terrain models typically are stored
using integer (fixed point) representation for elevations to
achieve compactness and because elevations are only known
with limited precision (an elevation may be given in meters as
a 16-bit quantity for example). Similarly, gray values to be
sent to a graphic output device are typically quantized to
relatively few levels because of the limited ability of the
human eye to discern small brightness differences and the
limited ability of the device to accurately reproduce these
(a typical output device may take values between 0 and 255).
The calculations can thus be carried out largely in integer
(fixed point) arithmetic and even a simple computer is
adequate.
XIII. USE OF LOOKUP TABLES
Some of the formulas for reflectance maps discussed later on
are quite elaborate and it would seem that a lot of computa-
tion is required to produce shaded output using them. Fortu-
nately it is possible to make the amount of computation
equally small in all cases by implementing the reflectance map
_tlgcp~bUg beginning.
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
Since elevations are quantized, so are the estimates of slope.
It is therefore not necessary that one be able to determine the
apparent brightness for all possible values of the gradient
(p, q). Further, it is reasonable to place an upper limit on
slope, so that only a finite number of possible values can occur
(For example, if slopes between -1.55 and +1.60 are con-
sidered, in increments of 0.05, then there are only 64 possibili-
ties for p and 64 for q, and a lookup table with 4096 entries
can be used). A second justification for the use of a lookup
table is the quantization of the gray values produced. It makes
little sense to calculate the apparent brightness with very high
precision only to coarsely quantize the result. A convenient
rule of thumb is that the number of possible discrete values for
each gradient component need not be more than the number
of gray levels available from output device. The final choice of
quantization must take into account both of the above
considerations.
One can separate the estimation of slope from the calcula-
tion of gray value, and produce an intermediate file of coded
surface gradient values. This file need not be larger than the
original terrain model if the gradient is quantized properly
(if p and q can each take on 64 values, each gradient can be
encoded as a 12 bit value). The code in the lookup table can
be based on ways of expression surface orientation other than
in terms of components of the surface gradient. In any case, a
file of surface orientation codes can be fed through a lookup
table procedure to produce the final output. In this fashion
different reflectance maps, encoded as different lookup tables,
can be applied to a terrain model easily, with little more effort
than reading and writing a file. The illustrations here were
produced this way.
Many gray-level raster displays have a translation table be-
tween the image memory and the digital-to-analog converter
driving the cathode ray tube intensity control. The quantized,
packed reflectance map can be loaded into this lookup table,
while the image memory is loaded with the coded slope
matrix. This allows one to view the same terrain with a
variety of assumed reflectance properties simply by reloading
the translation table, which is small compared to the image
memory.
XIV. TAXONOMY OF REFLECTANCE MAPS
Here we have discussed some of the issues one is likely to
encounter when developing a program that produces shaded
output. In the remainder of this paper we will analyze a num-
ber of proposed hill-shading methods in terms of their equiva-
lent reflectance maps. Notational tools will be introduced as
they are needed. Rather than proceed in strict historic order,
we will discuss relief shading methods in the following groups:
rotationally symmetric reflectance maps-gray tone de-
pends on slope only;
2) methods based on varying line spacing or thickness to
modulate average reflectance;
3) ideal diffuse reflectance and various approximations
thereto;
gray tone depends only on the slope of the surface in the
direction away from the assumed light source;
methods depending on more elaborate models of diffuse
reflectance' from porous material, such as that covering
the lunar surface;
models for gloss and lustrous reflection-smooth surface,
XV. AVERAGE REFLECTANCE OF EVENLY
SPACED DARK LINES
Some early methods for hill shading achieve the desired con-
trol of gray tone by varying the spacing between printed lines.
One advantage of this approach is the ease with which such
information can be printed, since it is not necessary to first
screen a continuous tone image. One disadvantage is the con-
fusion created when the lines used for this purpose are layed
on top of others portraying planimetric information. While
the directional textural effects of the lines are important in
conveying information about shape, we concentrate here on
the average reflectance.
Consider inked lines with reflectance rb covering an area of
paper with reflectance r,.,, (see Fig. 6). The ratio of the area
covered by ink to the area not covered is the same as the ratio
of the width of the lines to the width of the uninked spaces.
This in turn equals b/w, where b is the width of the inked line
and w the width of the uninked space measured along any
direction not parallel to the lines. If we ignore diffusion of
light in the paper, then the average reflectance of the surface is
R = (wr,,, + brb)/(w + b)
R = rw - b(r,,, - rb)/(w + b).
If, for example, the paper reflects all the incident light, and
the ink none, then rW = 1 and rb = 0, so that R = 1 - b/(w + b).
XVI. SLOPE OF THE SURFACE IN AN
ARBITRARY DIRECTION
In the calculation of gray value produced by some methods
of hill shading it is necessary to know the slope of the surface
in an arbitrary direction, given the slope p in the west-to-east
direction and the slope q in the south-to-north direction. Note
that p and q are the first partial derivatives of the elevation z
with respect to x and y, respectively. Consider taking an
infinitesimal step dx in the x direction and an infinitesimal
step dy in the y direction. The change in elevation dz is given
by
dz=pdx+qdy.
Along a contour line for example, the elevation is constant, so
that for a small step dx = a ds and dy = b ds, we can. write:
(p, q) ? (a, b) ds = 0
where "?" denotes the dot-product. The local direction of the
contours, (a, b) is of course perpendicular to the local gradient
(p, q)?
Now consider taking a small step in an arbitrary direction
(po, qo) say. That is let dx = po ds and dy = qo ds. The length
of the step, measured in the xy-plane is,
,/p0 2-2
qo ds.
While the change in elevation is,
dz = (pop + qoq) ds.
Consequently the slope, change in elevation divided by length
of the step, is
s=(pop +goq)I po +qo?
extended sorpproveoe -`or e se t2s~t ~c/03/10 : CIA-RDP88B00553RO00100280003 6 for
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
HORN: HILL SHADING AND THE REFLECTANCE MAP
.- b --a--o w --
Fig. 6. Magnified portion of surface covered with lines. The average
tone depends on the fractional area covered by the lines, as well as
the reflectance of the paper and the ink.
x-axis, then, the above can also be written,
s=pcosct +q sin a.
The direction in the xy-plane in which the slope is maximal
can be found by differentiating with respect to a. The direc-
tion of steepest ascent is (p, q) and the maximum slope
equals p2 + q2 .
XVII. LEHMANN'S BOSCHUNGSSCHRAFFEN (A)
One of the earliest methods for depicting surface shape using
a form of shading is that of Lehmann [2] , [3 1. Illustrations
based on ad hoc scales of increasing darkness as a function of
slope ("Schwarzegradscalen") had been published before, but
there was no systematic analysis of this approach until the
appearance of an anonymous publication attributed to Leh-
mann. In his method, short lines in the direction of steepest
descent, called hachures, are drawn with spacing and thickness
specified by rules that.ensure that the fractional area darkened
is proportional to the angle of inclination of the surface 0.
That is, steeper implies darker. The lines merge, producing a
continuous black area, when 0 exceeds some maximum value
00, typically 450 or 60?. The slope of the surface equals the
tangent of the angle of inclination or "dip." Using the expres-
sion for the slope in the direction of steepest ascent, we get,
tan0= p2+ q2.
Consequently, the average reflectance is,
R (p, q) = rw - (rw - rb) tan-1 p2 + q2 /00.
When the angle of inclination exceeds the maximum, the lines
coalesce and R (p, q) = rb. We can also write the above in
another form,
R'(0, 0) = rw - (rw - rb) (Oleo).
Here, 0, the azimuth of the direction of steepest descent, does
not appear in the formula. on the right, since apparent bright-
ness in this case depends only on the magnitude of the slope.
The direction and magnitude of the surface gradient can be
found from a map prepared according to Lehmann's rules.
The direction of steepest descent lies along the hachures, while
the slope is directly related to the average tone that results
from the width and spacing of these lines. In analyzing his
method we have concentrated on calculating the average re-
flectance produced in the printed product. It should be pointed
out that this method also gives rise to textural effects that will
not be discussed.
Another interesting aspect of Lehmann's method is that the
lines or hachures were drawn starting on one contour and end-
ing on the next. This greatly contributed to the later develop-
ment of the contour representation (Isohypsen) for terrain
Fig. 7. Spacing between successive contour lines along a given direction
on the topographic map.
surfaces, that was to ultimately replace most of these early
attempts at portraying surface shape [7] .
XVIII. CONTOUR DENSITY (B)
Another method is based on the observation that lines on a
contour map are more crowded in steep areas and that this
crowding leads to darkening of tone or average gray value.
This side effect may be helpful in visually conveying informa-
tion about the nature of the surface. In order to calculate the
dependence of the average local reflectance on the gradient
(p, q), we have to determine the spacing of contour lines on
the map. We assume that the surface is locally smooth and can
be approximated by a plane, at least on the scale of the spacing
between contour lines (if this is not the case, aliasing, or under-
sampling problems occur in any case).
Consider a portion of the surface with slope s in some direc-
tion not parallel to the contour lines (see Fig. 7). Assume that
the map scale is k and the vertical contour interval 8. Then it
is clear that the spacing between contours on the map d can be
obtained from the formula for slope,
s = S /(d/k).
If we take the cross section of the surface in the direction of
steepest ascent, then s = p2 + q2 . As a result we can write,
d=kS/ p2+q2.
On the map, d = b + w. That is, the spacing between contours
is the sum of the width of the contour lines and the width of
the blank spaces between them. The average reflectance then
is,
R (p, q) = rw - (b/k5) (rw - rb) p2 + q2.
The result can also be expressed as,
R'(0, 0) = rw - (b/kb) (rw - rb) tan 0
where 0 is the inclination of the surface. The above expres-
sions only hold if w is not negative. When the slope is too
steep, contour lines overlap, and the average reflectance is
simply equal to rb. In the special case that rw = 1 and rb = 0,
the above simplifies to,
R (p, q) = 1 - (b/k5) p2 + q2.
Typically (b/kS) may equal 1 or 1 / / .
XIX. DIFFUSE SURFACE UNDER VERTICAL .
ILLUMINATION (C)
The methods discussed so far produce tones that depend on
the magnitude of the gradient only, not its direction. This is
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
similar to the effect one would obtain if a physical model of
the terrain was illuminated vertically, with the light source
placed near the viewer. An ideal diffusing surface has an ap-
parent brightness that is proportional to the cosine of the inci-
dent angle i as discussed later. This is the angle between the
direction of the incident rays and the local normal, which, in
the case of vertical illumination, is just 0. Therefore,
R'(0,0)=cos0
R(p,q)=1/ l+p2+q2.
Instead of illumination from a point source, one may consider
the effect of a distributed source. A uniform hemispherical
source illuminating a diffusely reflecting surface leads to a
result of the following form [ 1401,
Rl(0, ~) = cost (0/2) = 2 + 1 cos 0
R(p,q)=(I+1/ l+p2+q2)/2.
This reflectance map leads to flatter, even less interpretable
pictures, since the range of reflectances has been halved and all
reflectances have been shifted upwards by a half. In the deri-
vation of the formula above, reflection from the surrounding
terrain surface is ignored. If the terrain surface diffusely re-
flects a fraction p of the incident light, the constant term in
the above expression is increased from 2 to (1 + p)/2, while
the coefficient of cos 0 decreased from to (1 - p)/2. It is at
times suggested that a component of surface brightness due to
distributed illumination from the sky be added to that result-
ing from oblique illumination. This however typically detracts
from the shaded result, rather than improving it.
The methods discussed so far give rise to rotationally sym-
metric reflectance maps, that can be described adequately by a
single cross-section, showing tone versus slope [ 11 , [351 , [36].
This representation has sometimes been misused for asymmetric
reflectance maps, where it does not apply. Rotationally sym-
metric reflectance maps produce shaded images that are diffi-
cult to interpret. Moving the assumed light source away from
the overhead position gives rise to better shaded map overlays,
but forces us to introduce some new concepts.
XX. THE SURFACE NORMAL
The surface normal is a vector perpendicular to the local
tangent plane. The direction of the surface normal n can be
found by taking the cross-product of any two vectors parallel
to lines locally tangent to the surface (as long as they are not
parallel to each other). We can find two such vectors by re-
membering that the change in elevation when one takes a small
step dx in the x-direction is just dz = p dx, while the change in
elevation corresponding to a step dy in the y-direction is dz =
q dy. The two vectors, (1, 0, p) dx and (0, 1, q) dy, are, there-
fore, parallel to lines tangent to the surface and so their cross-
product is a surface normal
n = (1, 0, p) X (0, 1, q) = (-p, -q, 1)?
Note that the gradient (p, q) is just the (negative) projection
of this vector on the xy-plane. A unit surface normal N can
be obtained by dividing the vector
%F1 +p +q2.
n by its magnitude n =
time ilhelnfulCt~vUCt~1R ft; Co~Yd1Ti'f25 Si Yl~c 1
Fig. 8. Definition of the azimuth angle 0 and the zenith angle 0. Here,
azimuth is measured counter-clockwise from the x-axis in the xy-
plane, while the zenith angle is measured from the z-axis.
can then be given as an azimuth angle 0 measured. anticlock-
wise from the x-axis, and a polar or zenith angle 0 (see Fig. 8).
(In navigation, the azimuth angle is usually measured clock-
wise from north, and the elevation angle is given instead of the
zenith angle. These are just the complements of the angles
used here.) The unit vector in the direction so defined equals,
N = (cos 0 sin 0, sin 0 sin 0, cos 0).
To find the azimuth and zenith angle of the surface normal we
identify components of corresponding unit vectors. Then,
sin-q/+q2
cos 0 _ -p/ p2 + q2
sin0= p2+q2J 1+p2+q2
cos0=1/ 1+p2+q2.
Conversely,
p=-COS 0tan0
and
q=-sin 0tan0.
We will find it convenient to use both vector and spherical
coordinate notation to specify direction.
XXI. POSITION OF THE LIGHT SOURCE
The reflectance maps discussed so far are rotationally sym-
metric about the origin, only the magnitude of the gradient,
not its direction affecting the resulting gray value. This cor-
responds to a situation where the light source is at the viewing
position. Most hill-shading methods have the assumed light
source in some other position, typically in the northwest, with
a zenith angle of around 45?(00 = 45?, 00 = 135?). The unit
vector,
S = (cos 00 sin 00, sin 00 sin 00, cos 00)
points directly at the light source. A surface element will be
illuminated maximally when the rays from the light source
strike it perpendicularly, that is, when the surface normal
t source. By identifying components in the
gh
l
i
points at the
urce
t
l
~
~
iorat~ exRPPPAO AW99 1049 @R l' (nti ng a at tio, the 1) source
~,`fi I
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
HORN: HILL SHADING AND THE REFLECTANCE MAP
one finds that the components of the gradient of such a sur-
face element are,
Po = -cos 00 tan 00
qo = -sin ?o tan 00.
When the source is in the standard cartographic position, this
means,
Po = i/VT
qo 1/VT?
This standard position for the assumed light source was proba-
bly chosen because we are used to viewing objects lighted from
that direction [ 1 ] . When we look at nearby objects in front of
us, our body blocks the light arriving from behind us. Further,
when writing on a horizontal surface, many of use find our
right hand blocking light coming from that direction. We thus
often arrange for light sources to be to the left, in front of us.
While we can certainly interpret shading in pictures where the
light source is not in this standard position, there seems to be
a larger possibility of depth reversal in that case, particularly if
the object has a complex, unfamiliar shape.
Returning now to the specification of the position of the
light source, we find two identities that will be helpful later.
cos(0- 00)=(pop+qoq)/[ P2 +q2 Po +qo]
PoP+qoq=tan0 tanOo cos(0- 00).
It also follows that the slope of the surface in the direction
(Po, qo) away from the light source is,
s = tan 0 cos (0 - 00).
XXII. TANAKA'S ORTHOGRAPHICAL RELIEF METHOD (D)
A method proposed by Kitiro Tanaka in 1930 [101, [ 11 ] ,
involves drawing the lines of intersection of the surface with
evenly spaced inclined planes. These planes are oriented so
that their common normal points towards an equivalent light
source (see Fig. 9). Thus slopes tilted away from this direction
have contours spaced closely, giving rise to heavier shading
than that on horizontal surfaces, while surfaces lying parallel
to the inclined planes are lightest. As in Lehmann's method,
some information may be conveyed by the directional texture
of the contours. Here we concentrate on the average reflec-
tance only.
A contour is the intersection of the terrain's surface z =
z(x, y) with a plane. The equation z = z0 applies to a hori-
zontal plane appropriate for ordinary contours. For "inclined
contours" an inclined plane is used with an equation of the
form
1) ? (x,Y, z)/l1 + Po + qo = zo.
(-Po, -qo,
The vector (-po, -q0, 1.) is perpendicular to the inclined
planes. Ordinary contours represent the locus of the solution
of z (x, y) = z0 , while inclined contours are the loci of solu-
tions of the equation,
.
[-pox - q0y+z(x,Y)]/ 1 +Po +qo =z0r
Fig. 9. Side-view of a hill cut by inclined planes. Viewed from above,
the lines of intersection crowd together where the surface slopes away
from the equivalent source. Conversely, there are no lines where the
terrain surface is parallel to the inclined planes.
of this equation! All we need are the slopes of this new sur-
face. Differentiating the above expression with respect to x
and y, we get,
q'=(q-qo)l 1+Po+qo?
(P - Po)2 + (q - qo)2
R(P,4)=rw-(b/kS)(rW rb) l+P2+q2 0
We obtain the expression for contour density, derived earlier,
when po = qo = 0. Also, in the special case that rb = 0, rx, = 1,
Po = 1/NF2 , and Q0 -- -1/-\r2 ,
R(P, q) = 1 - (b/k&) (P - 1/VT)2 + (q + 1/VT)2/VT?
It is sometimes useful to express the apparent brightness as a
function of the azimuth 0 and zenith angle 0 of the surface,
normal. If we let 00 be the azimuth and 00 the zenith angle
of the normal to the inclined planes, then the formula can be
rewritten as follows:
R'(0, 0) = rW - (b/ks) (r,,, - rb) cos 00
tang 0 - 2 tan 0 tan 00 cos (0 700) + tan 200.
When 00 = 45?, rb = 0 and rti? = 1, then, as Tanaka showed
[101, [111,
R'(0, 0) = 1 - (b/k6),\11 - sin 20 cos (0 - 0o)/(N12 cos 0).
How does one choose the parameter (blkS)? Tanaka felt
that the shading produced by his method should match that
seen on a surface covered with an ideal material called a per-
fect diffuser. The apparent brightness of such a surface varies
with the cosine of the incident angle, between the surface
normal and a vector pointing at the light surface. He intro-
duced a parameter called the line factor. It is the ratio of the
width of the inked line b to the interval between inclined con-
tours for a horizontal surface k5 /sin 00. The line factor is just,
(b/kS) po+qo/ l+po+qo.
Ta ka ro osed varying the line width b in order to produce
We can now apply oar 88f~k ~ ~g~r deW ~}R 90
to the modified surface z (x, y) de med y e e n e sh ~~t8$@O~i4~0f#: t diffuser, but
 Approved For Release 2006/03/10 : CIA-RDP88BOO553ROO0100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
Fig. 10. "Block-diagram" representation of terrain surface. This is an isometric projection of a series
of uniformly spaced vertical profiles of the surface viewed from the southeast. Note the shading
effect due to the variation in line spacing.
realized the impracticality of this approach for all but poly- nate system of the surface rather than one oriented with
hedral surfaces [ 10] , [ II ] . Resigned to using a fixed line respect to the viewer. Details may be found in Appendices
width, he chose to optimize the line factor by considering the B and C, where contour density shading and Tanaka's inclined
brightness distribution on a spherical cap extending to 45? contour method are shown to be special cases of this more
slope. With the source at 45? elevation, the least deviation general situation.
from the brightness distribution one would see if the surface
was a perfect diffuser is obtained when the line factor equals
0.3608. Consequently, (b/ks) = 0.3608. Finally then,
R (p, q) = 1 -- 0.3608 (p - 1/\)2 + (q +l /V)2 .
It is unfortunate that this method later gave rise to some mis-
understanding as well as a less rigorous hybridized form [ 15 1.
A common representation for relief form is the block dia-
gram, an oblique view of a series of equally spaced vertical
profiles [971-[1021. The projection typically is orthographic,
although at times a perspective projection is utilized. Surfaces
not visible to the viewer are eliminated (see Fig. 10). Shading
can of course be applied to oblique views as may be done in
sophisticated flight simulators of the future. We concentrate
here on map forms that provide for superposition of plani-
metric information however, and digress only to point out that
part of the appeal of block diagrams lies in their implicit shad-
ing, due to the variation in the spacing of lines.
Following the discussion in the last section, it is clear that
the equivalent light-source position is in the horizontal plane
at right angles to the vertical cutting surfaces. The analysis
just presented then applies directly. Things are a little more
difficult if the result is to be expressed in terms of the coordi- terrace l1~a
ff d h~
Approved For Release 2006/03/10 : CIA-RDP8 >3 Od~~~O~~b1~Ob~3- the light source has
XXIII. WIECHEL'S CONTOUR-TERRACE MODEL (E)
Imagine a three-dimensional model of the terrain built by
stacking pieces of some material cut according to the shape of
the contours on a topographic map [8]. If the thickness of
the material is chosen correctly the model will be a scaled
approximation of the terrain, looking a little like a tiered cake.
Illuminating this construction with a distant point source will
give rise to a form of shading since each contour "terrace"
casts a shadow on the one beneath it (see Fig. 11). Wiechel
[8] was the first to analyze the reflectance properties of such
a surface. In order to calculate the average brightness of a
portion of the model, when viewed from above, we must
determine the width of the shadow relative to the width of the
terrace.
The width of the shadow, measured perpendicular to the
contours, varies, depending on the orientation of contours
relative to the direction of the rays from the source. For ex-
ample, when measured this way, the width is zero where the
contour is locally parallel to the projection of the rays on the
xy-plane. Measured in a vertical plane containing the light
source however, the width of the shadow is constant, since the
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
5)
Fig. 11. Shadows cast in the contour terrace model. The width of the
shadows, measured perpendicular to the contours, varies with the di-
rection of the contours relative to the direction of the incident rays.
/ \
Fig. 12. Section of the contour terrace model in a vertical plane con-
taining the light-source. The width of the shadow b measured in this
plane is constant, while the width of the terrace b + w depends on the
slope of the surface in a direction parallel to the projection of the in-
cident rays on the ground plane.
a zenith angle 00, the contour interval is S, and the map scale
k, then,
tan 00 = (b/kS )
tan 00 = po + q o
To calculate the average brightness we must know the width
d of the terrace in the model, measured in the same vertical
plane. The slope in this plane evidently is just
s = -kS/d.
We know that the slope of a surface in the direction (po, qo)
s=(pop +goq)l po+qo?
Solving for d from the last two equations and for b from the
two before them, we get
bld =-(pop + qo q)?
For example, when the local surface normal (-p, -q, 1) is per-
pendicular to the direction to the source (-po, -q0, 1), their
dot-product is zero and b/d = 1. The terrace is then covered
tour interval and the map scale have cancelled, as one might
have predicted.
When (pop + q o q) < -1, shadows coalesce and no further
increase in bld is possible. When, on the other hand, (pop +
qo q) > 0, the slope is facing towards the light source. This
means that no shadow is cast. In this model, shading only
occurs on slopes facing away from the source, while those
facing towards it are all uniformly bright. This is certainly
not what one would expect of a real surface and suggests that
the contour-terrace model has some shortcomings. This is not
surprising since apparent brightness depends on surface orien-
tation, not height, and while the model represents height with
reasonable accuracy it does a poor job of modeling surface
orientation. Indeed the surface of the model is mostly hori-
zontal, with some narrow strips of a vertical orientation. The
latter are not even visible from above.
Wiechel noted that light would be reflected from these verti-
cal surfaces onto the terraces [8]. The surface thus appears
brighter, viewed from above, near vertical surfaces facing
towards the light source. He made the simplifying assumption
that reflection produces uniformly bright patches with the
same shape as shadows that would be cast were a source to be
placed opposite the actual light source. This is not a reason-
able assumption unless the vertical surfaces are made of narrow
mirror facets, each oriented perpendicular to the direction of
the incident light! In this case, surfaces illuminated by reflec-
tion as well as by direct light have a brightness twice that of
those illuminated only by direct light. This version of the
model is fortunately simple enough to be amenable to analysis.
First note that, if we assume the surface to be an ideal diffuser,
then the brightness of horizontal surfaces that are neither
shadowed nor illuminated by reflection equals the cosine of
the zenith angle of the source. Therefore, let rb = 0 and rw =
cos 00, where
cos00=l/ 1+p0+qo
and so
R(p,q)=(l+pop+qoq)/ 1+po+qo
R'(0, 0) _ [ 1 + tan 0 tan 00 cos (0 - 00)] cos 00.
When the source is in the standard position (northwest at 450)
this becomes,
R(p,q)_[I+(p-q)/ ]/ / .
Note that here apparent brightness already becomes equal to
one when the angle of inclination is about 30.36? towards the
light source. This may be contrasted with the case of the ideal
diffuser, to be discussed later, where it reaches one only for an
inclination of 45?. Wiechel used this model as the second
approximation to the ideal diffuser (the first will be discussed
later) and expressed his result as [8]
(cos i/cos e)
where i is the incident angle, and e is the emittance angle, here
equal to 0. These angles will play an important role in the dis-
cussion of more recent methods later on.
According to Raisz and Imhof [ 1 ] , [271-[29] terraced con-
tour models were used in the late 1800's. An early example is
an alpine excursion map published in 1865 that employed
exactly by the shadow. In the above expression both the con- "contour shadows" [ 1 ] . The first attempts at photography of
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
obliquely illuminated surfaces also used terraced terrain models
[27]. Wiechel probably was influenced by these early efforts
when he chose to develop this method for hill shading.
XXIV. WIECHEL'S HELLIGKEITSMAASSSTAB
Wiechel based his method for irregular surfaces on that de-
veloped earlier by Burmester for regular surfaces [9]. In order
to make his approach practical he needed a graphical device
for translating measurements of contour interval and direction
of steepest descent into gray tones. The "Helligkeitsmaass-
stab" (his spelling) is arranged so that these measurements can
be transferred directly, and the correct tone determined from
a series of isophotes, contours of constant brightness. Steep
slopes, with small contour intervals correspond to points near
the origin of this diagram, while those of gentle slope map into
points further away.
His diagram therefore is a sort of inside-out reflectance map.
The main difference is that radial distance from the origin in
gradient space is proportional to tan 0, while it is proportional
to cot 0 in this early precursor. This corresponds to a con-
formal mapping operation referred to as inversion with respect
to the unit circle. Wiechel showed that his diagram corre-
sponded to the image of an appropriate illuminated logarith-
moid made of the desired material. The equation of this sur-
face is z = -log x2 +y2. The reflectance map, by the way,
can be thought of as the image of a paraboloid [ 1381.
It is indeed unfortunate that Wiechel's construction was
ignored. Wiechel developed two shading methods that did not
require this two-dimensional diagram. In each case apparent
brightness depended only on the slope of the surface in the
direction away from the light source. This property manifests
itself in the reflectance map in the form of parallel straight-line
contours. The effect is less apparent in Wiechel's diagram,
where isophotes become nested circles through the origin, with
centers along the line in the direction of the light source.
XXV. TANAKKA'S RELIEF CONTOUR METHOD (F)
Kitiro Tanaka, in 1939, developed an ingenious method
[24]-[26] for drawing the shadows one would see if one
looked at a contour-terrace model. His method is based on the
observation that the length of the shadow, measured in the
direction of the incident rays, is constant. Using a pen with a
wide nib one can trace the contours, while maintaining the
orientation of the nib parallel to the direction of the incident
rays (as in roundhand writing). Only those portions of the
contours are traced that correspond to slopes facing away
from the assumed light source. Tanaka used black ink on gray
paper for reasons that will become apparent. If the reflectance
of this paper is rg then,
R(p,q)=rg+(rg-rb)(pep+qoq)
provided (pop + q0 q) < 0, otherwise R (p, q) = rg.
Tanaka also came up with a way of modulating the average
reflectance of the paper in areas that corresponded to slopes
facing towards the source. His approach is somewhat analogous
to taking the negative of a picture of the contour-terrace model
obtained by illuminating it from the other side. Thus white
"shadows" are cast in the opposite direction to the black shad-
ows. These can be drawn with white ink on gray paper using
reflectance will be,
R(p,q)=rg- (rg- rw)(pop+qoq)
where r,,, is the reflectance of the white ink. When (pop +
q0 q) < 0, no "shadows" appear and R (p, q) = rg. Tanaka
combined the two methods, tracing contours using both white
and black ink. The corresponding reflectance map R (p, q)
equals one of the expressions above depending on whether the
slope locally faces away from or towards the assumed source.
He apparently also experimented with nibs of different width
for white and black ink. This corresponds to changing the ele-
vation of the assumed sources. If the width of the nib is b,
then the relationship is,
(b/kb) = tan 0o = pa + qo.
The results of this tedious manual method are most impressive
[241-[261. One can write the above expressions in the alter-
nate notation,
R'(0, 0) = rg + (rg - rb) tan 0 tan 00 cos (0 - (fio ),
when cos (0 -
00) < 0
R'(0,
0) = rg - (rg - r,,,) tan 0 tan 00 cos (0 - 00),
when cos (0 -
00) > 0.
Tanaka preferred a reflectance for the gray background half-
way between that of the black ink and the white ink. Placing
the light source in the standard position we get,
R(p,q)= [I +(p- q)/J]/2
R'(0, 0) _ [ 1 + tan 0 cos (0 - 00) ] /2.
This result can also be expressed as, (cos i cos g)/cos e? where g
is the phase angle, here equal to 0o . Note that except for scaling
by cos g, this is the same result as that obtained by Wiechel for
his contour-terrace model. One effect of this scaling is that
apparent brightness rises to one only when the angle of inclina-
tion is 45?, on the other hand, horizontal surface now have a
gray value of only 0.5.
XXVI. TANAKA'S HEMISPHERICAL BRIGHTNESS
DISTRIBUTION
Tanaka needed a way to display the dependence of tone on
surface orientation to permit comparison of the results pro-
duced by his two methods and what would be seen if the sur-
face modeled were an ideal diffuser. He chose an oblique view
of the brightness distribution on a spherical cap extending to
45? inclination [ 101 , [ 1 1 ] , [24] -[26 J. If the cap is increased
until it is a hemisphere, one obtains something like the reflec-
tance map. One difference is that radial distance from the
origin in gradient space is proportional to tan 0, while here it is
proportional to sin 0. Thus, while the reflectance map is a
central projection of the Gaussian sphere onto a horizontal
plane, this is a parallel projection. Put another way: we are
dealing here with an image of a hemisphere, while the reflec-
tance map is the image of a paraboloid.
Tanaka's oblique views of the distribution of brightness
versus surface orientation do not provide the quantitative
information available in a contour representation such as
the same method as before except that now the section of the Wiechel's. His method is nevertheless very helpful and it is
contours that correspond to slopes facing towards they light unfortunate that few seem to have paid any attention to it,
source are traced. Xis pprO ed I-or e1ease~Goro r j c]gips}~y$8B08onti gdO8s~OQf~ ppr p,Fiate forms. It is
 Approved For Release 2006/03/10: CIA-RDP88B00553R000100280003-6
HORN: HILL SHADING AND THE REFLECTANCE MAP
not uncommon for example to see the dependence of tone on
surface orientation shown as a curve depending on one vari-
able, slope, when it clearly depends on two, slope and the
direction of steepest descent, or equivalently, the two compo-
nents of the gradient.
XXVII. LAMBERTIAN SURFACES (G)
We now turn from graphical methods using variation in line
spacing and line thickness to those utilizing continuous tone or
halftone techniques. These are often based on a model of
what the terrain would look like were it made of some ideal
material, illuminated from a predetermined direction. The
result differs from an aerial photograph, since no account is
taken here of varying terrain cover, the light source is often
placed in a position that is astronomically impossible, and the
terrain model has been smoothed and generalized. Not being
like an aerial photograph is an advantage, since aerial photo-
graphs, taken with the sun fairly high in the sky, often do not
provide for easy (monocular) comprehension of surface
topography.
The amount of light captured by a surface patch will depend
on its inclination relative to the incident beam. As seen from
the source the surface is foreshortened, its apparent (or pro-
jected) area equal to its true area multiplied by the cosine of
the incident angle. Thus the irradiance is proportional to cos i.
Strangely, it is commonly assumed that the radiance (apparent
brightness) of the surface patch is also proportional to cos i.
This is generally not the case since light may be reflected dif-
ferently in different directions, as can be seen by considering a
specularly reflecting material.
One can however postulate an ideal surface that reflects all
light incident on it and appears equally bright from all viewing
directions. Such a surface is called an ideal diffuser or Lam-
bertian reflector and has the property that its radiance equals
the irradiance divided by 7r [ 1421, [ 143 1. In this special case
the radiance is proportional to the cosine of the incident angle.
No real surface behaves exactly like this, although pressed
powders of highly transparent materials like barium sulfate
and magnesium carbonate come close. Matte white paint, opal
glass, and rough white paper are somewhat worse approxima-
tions, as is snow [ 1311. Most proposed schemes for automatic
hill-shading are based on models of brightness distribution on
ideally diffusing surfaces [8], [10], [11], [49]-[57], [71],
[731, [741, even though there is no evidence that perception
of surface shape is optimized by this choice of reflectance
model. As we will see, reflectance calculations based on this
model are not particularly simple either.
The cosine of the incident angle can be found by considering
the appropriate spherical triangle (see Fig. 13) formed by the
local normal N, the direction towards the source S, and the
vertical V. One then finds the following, as Wiechel already
showed [81,
R'(6, 4)) = cos 00 cos 0 + sin 6o sin 0 cos (0 - 00).
Alternatively one can simply take the dot-product of the unit
vector N normal to the surface and the unit vector S pointing
towards the source [ 1381, [ 1401
(--p, -q, 1)(-po, -qo, 1)
( 1 +-p2+ q2 I +po +qo)
Fig. 13. Spherical triangle used in calculating the incident angle i from
the azimuth and elevation of the light-source and the azimuth and
elevation of the surface normal. The direction towards the viewer is
V, the direction to the source is S, while the surface normal is N.
then is
R(p,q)=(l+pop+g0q)I( 1+p2 +q2 l+po+qo).
When (1 + pop + q0 q) < 0 the surface element is turned away
from the source and is self-shadowed. In this case, R (p, q) = 0.
In the case of a point source of light at 450 zenith angle in
the northwest, the reflectance map becomes
R(p,q)1 +(p-q)/J I 1+p2+q2]?
XXVIII, PEUCKER'S PIECEWISE LINEAR
APPROXIMATION
(H)
The computation of gray value using the equation for the
cosine of the incident angle is complicated and slow because of
the appearance of the square root. Peucker [611 experimented
with a number of approximations that are easier to compute.
He found that an adequate, piecewise linear approximation for
slopes less than one, is
0.3441p
- 0.5219q + 0.6599,
for p + q > 0
0.5129p
-- 0.3441 q + 0.6599,
forp+q 0, R (p, q) > 1 and so all surfaces facing
towards the light source are white. No information is avail-
able to the viewer regarding surface shape in these areas. If
the assumed light source is in the standard position we get the
simple formula,
R (p, q) = 10(p-q)/Y `
Marsik also limited the density to a maximum of 0..7 to avoid
interference with planimetric information on the map.
XXXV. LOMMEL-SEELIGER LAW (N)
Many surfaces have reflectance properties that differ greatly
from those of an ideal diffuser. The photometry of rocky
planets and satellites has intrigued astronomers for many years
[ 121 ] -[ 130 1. Several models have been proposed to explain
the observed behavior. One of the earliest, developed by
Lommel [ 119 ] and modified by Seeliger [ 120 1, is based on an
analysis of primary scattering in a porous surface [ 1261, [ 128 ] .
Their model consists of a random distribution of similar par-
ticles suspended in a transparent medium and results in a
reflectance function that is given here in its simplest form,
I/[ 1 + (cos a/cos i )]
unless,cos i < 0, when the surface is self shadowed. Here i is
the incident angle, and e is the emittance angle, the angle
between the local surface normal and the direction to the
viewer, here equal to 0. The expression equals 1 /(1 + cos g)
when i = 0, where g is the phase angle, here equal to 00. Using
this value for normalization and remembering the expression
for cos i one finds
) 2(1 +pop + qoq)
R (
=-
R'(0, 0) = (1 + cos 00)/[ l + cos 0/(cos d cos 00
p, q
(1+ 1+p2
0 +q0)(1+p2+q2)
+ sin 0 sin 00 cos (0 - 00))]
Incidentally, this function does not satisfy Helmholtz's reci-
procity law [124], and therefore cannot correspond to the
[1+1/ l+po+qo]
reflectance of any real surface illuminated by a point source.
[1+
l+po+qo/(l+pop+qoq)]
XXXIV. MARSIK'S AUTOMATIC RELIEF SHADING (M)
Blachut and Marsik further modified Wiechel's approxima-
tion, partly as a result of their dissatisfaction with the fact that
a horizontal surface does not appear white when a perfectly
diffusing material is assumed [58], [59]. This may have
stemmed in part from early conventions in map-making where
horizontal surfaces were portrayed without hachures [41-[6].
Marsik also aimed for simpler calculations and considered the
slope in the direction towards the source. For some reason,
he proposed making the density of the printed result equal to
the tangent of the projected slope angle 0' (see Fig. 15). Den-
sity is the logarithm (base 10) of the reciprocal of the reflec-
tance. Applying the analogue rule to the upper spherical tri-
angle (see Fig. 16) one can show that,
0 = cos 0 sin 0' - sin 0 cos 0' cos (0 - 00).
tan. 0' = tan 0 cos (0 - Oo),
unless (1 +pop +qoq) < 0, when R (p, q) = 0. When the
source is in the standard position,
(1+1/J)[1+(p-q)/5]
(1+/2-)+(p-q)l-,r2-
The Lommel-Seeliger law has been used in automated relief
shading by Batson, Edwards, and Eliason [70].
Based on detailed measurements and modeling, Fesenkov
[ 1231, [127] and later Hapke [ 1281-[1301 further improved
the equations for the reflectance of the material in the maria
of the moon. Hapke imagined the surface as an open porous
network into which light can penetrate freely from any direc-
tion. His result has three components: the Lommel-Seeliger
formula for reflection from a surface layer containing many
scattering points of low reflectance, Schonberg's formula
[122] for reflection from a Lambertian sphere and a compli-
cated factor resulting from mutual obscuration of the particles.
The results of such investigations are often expressed in terms
of angles other than the ones introduced so far The I omm 1
RI WpVrbv~t Vor Re1?dase 2006/03/10: 44 P88BOG 5* 00400e8? in a way which
 Approved For Release 2006/03/10 : CIA-RDP68B00553R000100280003-6
HORN: HILL SHADING AND THE REFLECTANCE MAP
simplifies the problem of calculating the shape of the lunar
surface from shading in a single image [ 137] , [138 1, [ 151 ] .
The angles needed, luminance longitude and luminance lati-
tude, are defined in Appendix D.
XXXVI. MINNAERT'S REFLECTANCE FUNCTION (0)
Minnaert discusses a large variety of models for the reflec-
tion of light from rough surfaces [126]. He also proposed a
class of simple functions of the form,
cosK i cosK -1 e
intended to fit observations of the radiance of lunar material
while obeying the reciprocity law [ 124] . Here K is a parameter
to be chosen so that the best fit with experimental data is
obtained. This parameter is meant to he between zero and
one, with the above expression becoming equal to that for the
perfect diffuser when K := 1. We can normalize this expression
so it equals one when i = 0,
R(p, q) = cos' i cos'' -a/cosK-1 g
R(P, q)= [(1 +pop+qoq)/(l +p2 +q2)]K
(~ + p2 + q2 / l + Po + qo ).
XXXVII. PARTICULARLY SIMPLE REFLECTANCE
MAPS
Several methods discussed here have reflectance depending
only on the slope in the direction away from the assumed light
source, leading to parallel straight line contours in the reflec-
tance map. These include Wiechel's first and second "approxi-
mation," Tanaka's relief contour method, the "law" of Lom-
mel and Seeliger, Minnaert's formula when K = 2, as well as
Marsik's automatic relief shading. These methods are quite
effective in producing overlays that are easy to interpret. One
can construct more such reflectance maps, including some that
are even easier to calculate. One possibility, for example, is,
R(p,q)=2+2(p'+a)/b
P'=(PoP+goq)l Po+qo
is the slope in the direction away from the source. Values less
than or equal to zero correspond to black, while values greater
than or equal to one correspond to white. The parameters a
and b allow one to chose the gray value for horizontal surfaces
and the rapidity with which the gray values changes with sur-
face inclination. The simple program shown earlier (see Fig.
5) uses this form with a = 0, b = 1/VT and po = 1 /VT, qo =
-1/VT.
A simple alternative, where all possible slopes are mapped
into the range from zero to one is,
R(p,q)= 2+2 (p'+a)/ b2 +(p'+a)2.
This has the advantage that the reflectance does not saturate
for any finite slope and all changes of inclination in the verti-
cal plane including the source translate into changes in gray
level.
Another way to achieve this effect is the following, some-
what reminiscent of Lehmann's approach,
R (P, q) = 2 + (I /7T) tan-1 [(7T/2) (P'+ a)/b].
These three formulas are given in a form where the rate at
which the gray value changes with surface inclination is the
Fig. 17. Spherical triangles used to calculate the first off-specular angle
s. It is the angle between S, the center of the source, and S, the direc-
tion from which light is specularly reflected towards the viewer. Equiv-
alently, it is the angle between V, the direction of the viewer, and V',
the direction in which light from the center of the source is specularly
reflected.
S
XXVIII. GLOSSINESS-THE FIRST OFF-SPECULAR ANGLE
Not all surfaces are matte. Some are perfectly specular or
mirror-like. Since smooth, specularly reflecting surfaces form
virtual images of the objects around them, patches of high
brightness will appear when such a surface is illuminated by an
extended source, like a fluorescent light fixture, or by light
streaming in through a window. The size of the patches de-
pends on the solid angle subtended by the source as well as the
surface curvature, while the brightness distribution is that of
the source.
To study reflection of an extended source in a specular sur-
face, it is useful to introduce the "off-specular" angle s between
the direction S to the center of the source and the direction
S', of the point that is specularly reflected to the viewer (see
Fig. 17). This, incidentally, is also the angle between the
direction to the viewer V and the direction V' in which the
rays from the center of the source are specularly reflected.
We assume a circularly symmetric source, with brightness
L(s) at eccentricity s. This is the brightness the viewer observes
in. the specularly reflecting surface. Calculating the first off-
specular angle s is simple using the appropriate spherical
triangles:
(P)
cos s = cos 2i cos g - sin 2i sin g cos X
cos e = cos i cos g - sin i sin g cos X.
Here, i is the incident angle, between the local normal and the
direction to the source, e = 0, is the emittance angle, between
local normal and the direction to the viewer, while g = 0o is
the phase angle, between source and viewer. Eliminating
X from the two equations and expanding the sine and cosine
of 2i, one gets,
Substituting expressions in p and q for cos i, cos e and cos g
one can rewrite this as,
coss=[2(1+pop+qoq)/(1+p2+q2)- 1]J l+po+qo.
This result can also be obtained simply by finding the direc-
tion S' from which a ray must come to be specularly reflected
to the viewer V, by a surface element with normal N,
same at (p' + a) = 0. Approved For Release 2006/03/10: CIA-RDP88B0055 R0~('1 6 0(r03-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
where V= (0, 0, 1). The off-specular angle is the angle be-
tween S' and the center of the source S so
coss=S?S'=2(S?N)(V?N)-(S?V).
Note that the cosine of the first off-specular angle can be
calculated easily, without using trigonometric functions. The
contours of constant cos s turn out to be nested circles in
gradient space, with centers lying on the line from the origin
to the point (po, q,)). This can be seen by noting that the
locus of the point S', for constant s, is a circle about the point
S and that circles on the Gaussian sphere give rise to circles in
gradient space when projected stereographically [ 138] .
The cosine of the off-specular angle s equals one when con-
ditions are right for specular reflection, that is, when e = i and
g = i + e. This can be seen by setting e = i = g12 in the trigono-
metric expression for cos s.
XXXIX. Bul-TUONG'S FORMULA-SPECULAR SURFACE,
EXTENDED SOURCE (Q)
Having seen how to calculate the off-specular angle s, we can
now make a reflectance map, by assigning the distribution of
source brightness L(s). This function should be nonnegative
monotonically decreasing with s, and equal to one when s = 0.
For ease of calculation one choice might be
L(s) = cos" (s/2) = [!-(I (1 + cos s)]" l2
where n is a number that defines how compact the bright patch
is (a useful value might be around 20). So far, we have devel-
oped the reflectance map for a specular surface and a circu-
larly symmetric source. Many surfaces, such as glazed pottery
or smooth plastic, have both glossy and diffuse components
reflection. Specular reflection takes place at the smooth inter-
face between two materials of different refractive index, while
the matte component results from scattering of light that
penetrates some distance into the surface layer.
We can combine these two components as follows
R(p, q) _ [(1 - a) + aL(s)] cos i/cos (g/2)
where a determines how much of the incident light is reflected
specularly. The expression is scaled so that its maximum is
(approximately) equal to one. Here we have assumed the
source, while distributed, is compact enough so that the diffuse
reflection component can be approximated as cos i. The above
expression obeys the reciprocity law of Helmholtz [ 124]
which applies to real surfaces illuminated by a point source.
Bui-Tuong used a reflectance function similar to the one de-
rived above in his computer graphics work [ 113 1. He appar-
ently tried to model reflection from a surface that is not
perfectly smooth. This requires a different off-specular angle
however, as will be seen in the next section.
XL. LUSTER-THE SECOND OFF-SPECULAR ANGLE
Refulgency, gloss or shine can also appear when a point
source is reflected in a surface that is not perfectly smooth.
When a slightly uneven surface, of a material that gives rise
to metallic or dielectric reflection, is illuminated by a point
source, bright patches will be seen surrounding points where
the local tangent plane is oriented correctly for specular reflec-
tion. The size of these patches will depend on the roughness
of the surface and the surface curvature, while the distribution
of brightness will depend to some extent on the texture of
the microstructure of the surface.
Approved For Release 2006/03/10 :
Fig. 18. Spherical triangles used to calculate the second off-specular
angle s'. It is the angle between the actual surface normal N and a
surface normal N' oriented to specularly reflect rays from the source
towards the viewer.
In this case we will need to calculate the second off-specular
angle s' between the local normal N and the normal N' oriented
for specular reflection of rays from the source S towards the
viewer V (see Fig. 18). By considering the appropriate spheri-
cal triangles one finds,
cos s' = cos i cos (g/2) - sin i sin (g12) cos X
cos e = cos i cos g - sin i sing cos X.
Eliminating X from the two equations and expanding the sine
and cosine of the phase angle g, one finds,
cos s' = (cos i + cos e)/(2 cos (g/2))
cos s' = (cos i + cos e)/(NF2 1 + cos g).
This result can also be obtained by finding the vector N',
normal to a surface element oriented to specularly reflect a ray
from the source in the direction of the viewer V. That is,
N'=(S+V)/IS+VI.
The off-specular angle is the angle between the actual surface
normal N, and this vector N'
coss'=N?N'= [(S ?N)+(V ?N)]/-,r2- 1 +(S ? V).
The surface microstructure of an uneven surface can be
modeled by many randomly disposed mirror-like facets, too
small to be optically resolved, each turned a little from the
average local surface orientation. One can define a distribu-
tion P(s') describing what fraction of these microscopic facets
are turned away from the average local normal by an angle s'.
For ease of calculation one choice might be,
P(s') = cos" s'.
XLI. BLINN'S FORMULA-ROUGH SURFACE, POINT
SOURCE (R)
One can use the fact that a normal N' oriented for specular
reflection of the point source towards the viewer, lies in the
direction (-pl, -ql, 1), where
p 1 = - cos 00 tan (0o /2)
ql = -sin 00 tan (90/2).
We can also find N' by normalizing the vector (S + V), so that
CIA-RDP88B00553R000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
HORN: HILL SHADING AND THE REFLECTANCE MAP
its third component equals 1:
P1 =Po/[1+ l+po+qo]
q1 =q0/[1+ l+pa+qo].
A surface with gradient (pl, q1) is oriented just right to specu-
larly reflect a ray from the source to the viewer. This can be
seen by noting that when. p = pi and q = q1,
cosi=cose=l/ _p22
cos g = 2/(l + p2 + q2) - 1.
In any case,
toss'=(1 +plp+qlq)/( 1 +-p2+ 42 1 +pl +qi)?
Note that s' will tend to be (roughly) half of s when both
angles are small. Combining matte components of surface
reflection with those from the rough outer surface we get,
R(p, q) = [(1 - a) + aP(s')] cos i/cos (g/2).
The above reflectance map also obeys Helmholtz's reciprocity
law and is normalized so that its maximum is (approximately)
equal to one. Blinn and Newell give a similar reflectance func-
tion, claiming it was what Bui-Tuong had proposed [114 1.
The two are not the same however since the two off-specular
angles are different; in fact, the contours of constant s' are
nested ellipses in gradient space, while, as mentioned earlier,
the contours of constant s are nested circles. Indeed, Bui-
Tuong's model corresponds to reflection of an extended, rota-
tionally symmetric source in a specular surface, while the
model presented in this section applies to reflection of a point
source in a rough surface.
XLII. BLINN AND NEWELL'S MODEL FOR SPECULAR
SURFACES
One of the methods described by Blinn and Newell [ 114]
assumes a perfectly specular surface in which the world sur-
rounding the object is reflected. To make computations
feasible, they imagine the surrounding objects at a distance
great enough so that each part of the surround appears to lie in
essentially the same direction from every point of the surface
of the object. In this case one can imagine the brightness
distribution of the surrounding objects projected onto the
inside of a large sphere. The gray value used for a particular
surface patch then is found by computing the direction S' from
which a ray must come to be specularly reflected to the viewer
V, by a patch with surface normal N. We have already seen
that,
S' = 2(V ? N)N - V.
The appropriate gray value is then determined from the spher-
ical distribution of brightness. In practice the sphere is mapped
onto a plane by calculating the zenith angle, Bo and azimuth,
00 of S' [ 114] . The brightness distribution can be equally
well specified in gradient space [ 138], since it is also a projec-
tion of the Gaussian sphere.
Surface models incorporating randomly dispersed mirror-like
facets were first studied in the 1700's by Bouguer [ 1181. This
type of microstructure has been investigated extensively since
then, despite the difficulties of reasoning about the three-
dimensional nature of reflection from such surfaces. Recently,
Torrance and Sparrow further elaborated on these models
[1341, [135] in order to match more closely experimental
data showing maximum brightness for angles of reflection
larger than the incident angle. They included in their consid-
erations the effects of obstruction of the incident and emer-
gent rays by facets near the one reflecting the ray. Blinn sim-
plified and explained their calculations [115] and used them
in producing shaded images of computer models of various
objects. The overall result can be broken into a product of
three terms, one dependent on the distribution of facet orien-
tations, the second being the formula for Fresnel reflection
from a flat dielectric surface, while the third is the geometric
attenuation factor accounting for partial occlusion of one
facet by another. We will not discuss these models in any
more detail here.
Models for glossy or lustrous reflection have been used with
great success in computer graphics to increase the impression
of realism the viewer has when confronted with a synthetic
picture of objects represented in the computer. Unfortunately,
these methods do not seem to improve the presentation of
surface shape for cartographic purposes.
XLIII. COLORED SHADING
It is often said that quantitative information about the sur-
face cannot be obtained from relief shading [ 1 ] . Contour
lines on the other hand do allow measurements of elevation
and estimation of the gradient. Shading does provide some
information about the gradient too, but cannot be used to
determine both of its components locally, since only one mea-
surement is available at each point. Since we can perceive the
shape of objects portrayed by shaded pictures, it seems that
these local constraints do lead to a global appreciation of
shape, apparently based on our assumption that the surface is
continuous and smooth.
If two shaded images, produced with the assumed light
source in different positions, were available however, two
measurements could be made at each point allowing one to
determine the gradient locally [ 1411. It is inconvenient to
work with two shaded overlays; fortunately though, they can
be combined by printing them in different colors. In fact, yet
another overlay can be added in a third color, but it adds no
new information, since the two components of the gradient
are already fully determined by the first two.
Colored shading corresponds to illumination by multiple
sources, each of a different color. The exact color at each
point in the printed result is uniquely related to the gradient
at that point. Thus quantitative information is available in
this new kind of map overlay. Further, ambiguities present
in black and white presentations disappear. By positioning the
light sources properly, one can avoid problems occasioned by
the accidental alignment of ridge or stream lines with the
direction of incident light. Thus the need for ad hoc adjust-
ments of the azimuth of the assumed light source is removed.
Colored shading is easy to interpret in terms of surface shape
and effective in portraying surface form. It is unlikely how-
ever that it will be widely used because of the added expense
of printing and conflict with existing uses of color in cartog-
raphy to distinguish various kinds of planimetric information.
Amongst other things, color is now used to code height and
surface cover. Further, yellow is used in ordinary shading for
sun-facing slopes, while violet is used for shaded regions [ 1521.
This is thought to simulate the increased sky illumination
component in areas turned away from the sun.
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
 7
Fig. 19. Reflectance maps in the order in which they were introduced. The letter codes on subfigures match the letters found after the head-
ings of corresponding sections of this paper. Small crosses mark points in gradient space where slope-components are integer multiples of 1.
Where appropriate, a small square marks the gradient of a surface element that is perpendicular to the rays from the assumed light-source.
XLIV. SUMMARY AND CONCLUSIONS Methods that have been proposed for mechanizing the genera-
After a brief review of the history of hill-shading an efficient tion of relief shading were also treated. The automated method
method for providing shaded overlays was described. It de- described here is very flexible, since it can use any reflectance
pends on a lookup table containing sampled values of the map.
reflectance map. Traditional, manual methods were explored Eighteen of the reflectance maps described were plotted as
in terms of their equivalent reflectance maps, as were phenom- contour diagrams (see Fig. 19). The letter under a subfigure
enological model t ra ty c s onds to the letter apppearing after the heading of the 70- ro'veu ro ease'' bb!OJ! VT1: CIA DP88B00553R000100280003-6
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
 Approved For Release 2006/03/10 : CIA-RDP88BOO553ROO0100280003-6
HORN: HILL SHADING AND THE REFLECTANCE MAP
section in which the corresponding reflectance map is dis-
cussed. The first three (A, B, and C) are independent of the
direction of the gradient, depending only on the slope. These
give rise to rotationally symmetric diagrams. Six other dia-
grams (E, F, K, M, N, and P) show parallel straight lines. These
correspond to reflectance maps which depend only on the
slope in the direction away from the source. Reflectance maps
for perfectly diffusing (G) and glossy (Q and R) surfaces are
also included.
Shaded images of a mathematically defined surface were
then created using these reflectance maps (see Fig. 20). Several
of the subfigures give one a good appreciation for the shape of
the object. Assumption of a perfectly diffusing surface (G)
and Wiechel's modified brightness function (L) lead to good
results, while the images corresponding to glossy surfaces
(Q and R) are perhaps the most vivid.
The same set of reflectance maps was then used to make
shaded images of a re 'o S 't
- hick a digital 80
Approved For Release 2006/03/10: CIA-RDP88BOO553ROOO1O(~Z dOv IM g
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
Fig. 20. ShSrl:'wi yadfFaQ*Re sebWQMM& P? BO~~ AO~1G0g8 YO tl previous figure.
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
HORN: HILL SHADING AND THE REFLECTANCE MAP
Approved For Release 213OM3Mt(n4U -RDP88B00553R000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
ge Bppr jBoQowQ2 used for the previous
Fig. 21. Shaded images of a digital terrain model. 2Ru W jjq s W s
figure.
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
HORN: HILL SHADING AND THE REFLECTANCE MAP
Approved For Release 2004/O /10,a&J.}RDP88B00553R000100
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
Fig. 22. Shaded images of several digital terrain models. A particularly simple reflectance map was used. The lines in the subfigures correspond
to a length of 10 km on the surface. (a) Lake Louise, Alberta, Canada. (a) Gulf Islands, British Columbia. (c) Les Diablerets. Switzerland.
(d) Dent de Morcles, Switzerland. (3) Mexico City, Mexico. (f) Jewell Ridge, Virginia. (g) White Tail Butte, Wyoming. (h) Tehacnapi Moun-
tains, California. (i) Mount Index, Washington. (j) Mauritius, Indian Ocean.
terrain model was available (see Fig. 21). Some reflectance
maps appear much better than others in conveying an imme-
diate impression of surface shape. Rotationally symmetric
reflectance maps (A, B, and C), corresponding to overhead
illumination of the terrain, are not very good for example.
Perfectly diffuse reflectance (G) is not optimal either. In fact,
various approximations to the formula for a Lambertian re-
flector (H and K) seem to produce better results. Glossy
reflectance components (Q and R), while very useful in the
portrayal of regular objects, result in tones that are too dark
to be useful in a map overlay. We may also not be used to
Marsik's method (M), in which half of the surface is a fea-
tureless white, is clearly not very effective. Several of the
other methods require careful scrutiny before conclusions
about their adequacy can be made. Amongst the best are
Wiechel's modified brightness method (L) and the modifica-
tion of Brassel's method presented here (J). Several of the
methods depending on the slope in the direction away from
the light source appear to be quite adequate (F, K, N, and P).
These are to be recommended unless there are good reasons to
prefer one of the other methods.
Shaded images were created from several other digital terrain
seeing a geographiAp eorFsmRadetl sef2006/03/10 : Cl 86BD0@59M i 1928b6mimple methods (P).
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
HORN: HILL SHADING AND THE REFLECTANCE MAP
Fig. 22. (Continued.)
The terrain models differ widely in their quality, resolution metric detail [ 15:31. This is an area that has not received much
and origin. They do show the utility of the methods described attention so far. Another important issue relates to the appro-
here in presenting the information in a digital terrain model to priate scale for shaded overlays. Shaded overlays are useful
a human observer. for large scale neaps. For small scale maps it is necessary to
Shading is an important depth cue. The choice of reflectance generalize the surface to avoid the appearance of complex tex-
map should not be based on some ad hoc model of surface tures that may be difficult to interpret [ 11, [481, [731, [741 ,
behavior, experimental measurement of reflectance of some [ 1541, [ 1551. This nonlinear process of removing small hills,
material, or formulas that happen to be easy to calculate. In- ridges and valleys has not yet been satisfactorily automated.
stead, one should use a reflectance map that gives rise to an An as yet unexplored possibility depends on finely sampled
immediate, accurate perception of surface shape. terrain elevations. This is the ability of shading to show fine
It is important to arrange for the range of gray tones in the detail. Contour maps have to be carefully generalized or
A-le o th~t to avoid showin confusing detail on a scale smaller
shaded overlays to lJP8VA'Fur k9eib e? JIlO : C -~tDP88B00553R000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
than the contour interval. This is not the case with shading,
although historically the manually produced maps have always
shown only quite coarse features. We do not yet know whether
the textures produced by the shading method when working
from really fine terrain models will be confusing, or of great
value in identifying different types of terrain.
APPENDIX A
ROTATED GRADIENTS
It has been cartographic practice to assume a light source in
the Northwest at a 450 elevation above the horizon. It is help-
ful in this case to introduce a rotated coordinate system (see
Fig. 23) with
p'= (p - q)/P
Fig. 23. Rotated coordinate system that may be convenient when the
assumed light-source is in the northwest. The reflectance map is
symmetrical about the p'-axis.
4~=(p+q)lV~_2_.
and q,,
If Ox = Ay = A say, then the slopes in the Northwest to South-
east -,r2 and in the Southwest to Northeast direction, can be esti- pw [(z+o + z+- + za-) - (z-0 + z-+ + zo+)1 /4 A
mated particularly easily by combining the formulas for pw qw = [(zo+ + z++ + z+o) - ((z - + z-- + z - )] /4 \ A.
Approved For Release 2006/03/10 : CIA-RDP88B00553R000180280003-b?
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
HORN: HILL SHADING AND THE REFLECTANCE MAP
If one wishes to estimate the slopes for the center of the top-
right quadrant (in the unrotated coordinate system) rather
than the central point one may combine the expressions for
p1/2 and q1/2 to get the simple formulas,
p'12 = (z+o - zo+)/. A
q'1/2 = (z++ - zoo)/N/-2-
oo )/ / A.
One advantage of the rotated coordinate system stems from
the fact that models of surface reflectance considered here are
symmetric with respect to a line pointing towards the source.
That is, a surface element with slopes p' = po and q= qo say,
has the same apparent brightness as one with slopes p' = po
and q' = -qo. Thus a lookup table based on the rotated coor-
dinate system can be smaller, since only that half of the table
corresponding to q' > 0 need be stored.
So far we have assumed that the grid of the terrain model is
aligned with the geographical coordinates. If instead the whole
model is rotated anticlockwise by an angle 0, then slopes
p" and q" can first be estimated from the model as described
and then transformed as follows:
p=p" Cos 0-q"sin 0
q=p"sinb+q" cos 0.
Alternatively, the model can be resampled to produce a new
version on a grid aligned with the axes.
APPENDIX B
SHADING APPARENT IN BLOCK DIAGRAMS
We can analyze the shading apparent in block diagrams by
calculating the spacing between lines as a function of the sur-
face orientation. Let a local surface normal ben = (-p, -q, 1).
A series of parallel planes, with common normal s, cuts the
terrain surface. The intersections of these planes with the sur-
face are viewed from a direction specified by the vector v. It is
assumed that the viewer is at a great distance so that the pro-
files are projected orthographically along lines parallel to v (see
Fig. 24).
The line of intersection of one of the cutting planes with the
local tangent plane will be parallel to the vector n X s, since
the line lies in both planes and is therefore perpendicular to
the normals, n and s. Now construct a plane through the line
of intersection and the viewer. This plane, called the viewing
plane, contains both n X s and v. The normal e of the viewing
plane must therefore be perpendicular to both and can be
defined as,
e = (n?v)s-(s?v)n.
If we let p = (x, y, z), then the equation for the local tangent
plane can be written,
n ?p=Cn
for some value of the constant cn. Similarly, the equation of a
particular cutting plane is,
Fig. 24. The viewing plane contains the viewer and the line of intersec-
tion of the slicing plane with the terrain surface. Line spacing in the
block diagram equals the spacing between successive viewing planes.
The dotted line is, parallel to the vector v.
Different values of cs correspond to different cutting planes.
The plane corresponding to the value cs + dcs is separated
from the plane corresponding to the value cs by a distance
dcs/s, where s is the magnitude of the vector s. The equation
for the viewing plane is just
e?p=Ce.
Successive cutting planes will intersect the tangent plane in
parallel lines. These give rise to parallel viewing planes corre-
sponding to different values of the constant ce. The spacing of
these viewing planes is of interest, since it equals the spacing
of the lines in the orthographic projection. The plane corre-
sponding to the value ce + dce is separated from the plane
corresponding to the value ce by a distance of dce/e, where
e is the magnitude of the vector e. In order to relate the
spacing of lines in the block diagram to the spacing of the
cutting planes we need to find the relationship between dce
and dcs.
A point p on the line of intersection lies in all three planes
and therefore simultaneously satisfies the three equations
given above for these planes. Expanding the last one of these,
e ? p = cei we obtain,
(n-v)(s-p)-(s-v)(n' p)=Ce
Here, cn is fixed and
Ce and es is simply
(n - v) cg - (s . v) cn = Cc .
so the relationship between changes in
dce = (n ? v)dcs.
If the interval between cutting planes is S and the map scale
is k, then dcs/s = (kS). Consequently the spacing between
lines in the block diagram dcele is
d = kS (n - v)(s/e)
where e is the magnitude of the vector e = (n ? v)s - (s ? v)n.
Finally, we remember that
R(p, q) = rw - (rw - rb)(bld)
where b is the thickness of the lines. Thus
R(p, q) = r,,, - (b/kS)(r, - rb)(e/s)(1/n ? v).
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
The view vector is tangent to the surface when n ? v = 0. When
this dot product becomes negative, the surface is turned away
from the viewer and should not be visible. Also note that
d = kS, when s ? v = 0. One should therefore chooses and u so
that they are not orthogonal, to avoid getting only evenly
spaced parallel lines.
In the case of perspective projection, line density will increase
with distance, and the resulting reflectance will be lowered
because of a change in the effective scale factor k. If the
projected profiles are plotted on a raster device, one has to
also take into account the fact that the number of dots per
unit line length is not constant. The dot density varies as
max [ I cos 0 I , I sin 01 ] , where 0 is the angle between the line
and the direction of the raster. This variation should be in-
cluded if an accurate reflectance map is to be derived for out-
put of this form.
APPENDIX C
ISOMETRIC VIEWS OF VERTICAL PROFILES
The transformation between the terrain coordinate system
and that of an observer viewing the terrain obliquely can be
found by multiplying a rotation matrix corresponding to rota-
tion by 0, about the x-axis with a matrix corresponding to
rotation by (7r/2 + O,) about the z-axis, where 0, is the azimuth
and 0? is the zenith angle of the direction specified by the
vector v. If the coordinates in the observer's system are x', y',
and z', one finds,
x' = -sin O0x + cos O ,y
y' = -cos 0v cos 0?x - sin 0? cos 00y + sin 0?z
z' = +cos 00 sin 0?x + sin 0? sin 0,y + cos 0?z.
In the case of orthographic projection, the values of x' and
y' are simply multiplied by the map scale k, to determine
coordinates in the block diagram.
The general formula derived in Appendix B applies to all
combinations of viewpoint and cutting plane orientation. It
is interesting to look at a few special cases however. We can,
for example, check our result for the contour interval in an
ordinary contour map. Here n = (-p, -q, 1), as always, and
s = (0, 0, 1), since we are considering the intersection of
the surface with horizontal planes. Further, v = (0, 0, 1)
since the viewer is vertically above the surface. Here then
s = 1, n ? v = 1, and e (p, q, 0). The line interval is, therefore,
d = (kS)l P2 + q2
The same reflectance map is obtained as before. Slightly more
complicated is the case of Tanaka's inclined contours, where
s = (-po, -qo, 1). Here, again, n ? v = 1, while,
. e = (p -Po,q-qo,0)
S= l+po+qa.
The line interval is, therefore,
d =(k6)[ 1 +Po +qol (P - Po)2 +(q - qo)2 ].
A result leading to the same reflectance map as the one derived
before.
Finally, consider profiles running West to East, that is,
s = (0, 1, 0). The resulting traces may be viewed isometrically
Fig. 25. Luminance longitude a and luminance latitude 0 of a surface
element are defined as the longitude and latitude of a patch on a
sphere with the same orientation. Longitude and latitude are measured
relative to the luminance equator through the light source S and the
viewer V.
diagram. Then v = (1, -1, 1). Consequently, n ? v = (1 - p + q)
ands ? v = -1. Further, e = (-p, I - p, 1) and hence,
d=(0)(1 - p+q)/(N[2 1 -p+p2
So, if rb = 0 and rw = 1,
R(p, q) = 1 - \ (b/kS) 1 - p + p2 p + q).
Similarly, for profiles running South to North, s = (1, 0, 0),
and,
R(p, q)= 1 - (b/k6) 1 +q+q2/(1 - p + q).
At times two orthogonal sets of slicing planes will be used,
producing a mesh on the surface. The reflectance map corre-
sponding to this case can be found by adding the last two
formulas and subtracting one from the result.
APPENDIX D
LUMINANCE LONGITUDE AND LUMINANCE LATITUDE
A convention for specifying the orientation of the surface
element relative to the direction of a light source and the
viewer has become established in the work on planetary and
lunar photometry. Imagine a sphere illuminated by a light
source above the point S, viewed by an observer above the
point V (see Fig. 25). These two points define a great circle
which is called the luminance equator. Points on the sphere
can be referenced using the longitude a measured from the
point V along the equator, and the latitude 0.
All possible surface orientations can be found on the sphere,
and each surface orientation can be identified with some point
N say. The luminance longitude and luminance latitude corre-
sponding to a particular surface orientation are the longitude
and latitude of N. It is not difficult to show that,
from the Southeast, A bV@t 1F H RI 2 W6pW'0 : CFA1? W8$BI 3FW0041 O2 00[i26- 2IEG + E2 ]
 HORN: HILL SHAD1 PANDeTdHE REFLREleaseE 2006/03/10 : CIA-RDP88B00553R000100280003-6 MAP
where we have used the shorthand notation, I = cos i, E = cos e,
and G = cos g. These results can also be expressed in terms of
the components of the gradient:
tans=(pop +404)/ Po+qo.
So tan a is simply the slope in the direction away from the
source. Now,
1+2IEG-(I2+E2+G2)= (gop-Poq)2
[(I +p2 +q2) (l +po +q2
0)]
I2 2IEG+E2 = [(PoP+gog)2 +(P0 +qo)]
[(1 +P2 +q2) (1 +P0 +qo)]
tan(3=(qop- Poq)l (Pop+goq)2 +(po +qo).
The Lommel-Seeliger law can be expressed in terms of lumi-
nance longitude and luminance latitude as,
cos (a +g)/[cos a + cos (a + g)]
and it is clear from this form that scene radiance is indepen-
dent of luminance latitude.
ACKNOWLEDGMENT
The author would like to thank Kurt Brassel, Thomas
Peucker, George Lukes, Robert McEwen, Marsha Jo Hannah,
and James Mahoney for generously supplying digital terrain
model. Blenda Horn helped in the preparation of the text and
Karen Prendergast created the figures. Helpful comments were
provided by Robert Sjoberg, Katsushi Ikeuchi, and William
Silver. Encouragement by Thomas Peucker, Kitiro Tanaka, and
Kurt Brassel was instrumental in the generation of this paper.
[151 A. H. Robinson and N. J. W. Thrower, "A new method for ter-
rain representation," Geographical Rev., vol. 47, no. 4, pp. 507-
520, Oct. 1957.
1 161 A. H. Robinson, "The cartographic representation of the statis-
tical surface," Int. Yearbook of Cartography, vol. 1, pp. 53-63,
1961.
[171 N. J. W. Thrower, "Extended uses of the method of orthogonal
mapping of traces of parallel, inclined planes with a surface,
especially terrain," Int. Yearbook of Cartography, vol. 3, pp.
26-28, 1963.
[181 C. M. King, Techniques in Geomorphology. New York: St.
Martin's Press, 1966, pp. 255-256.
[191 T. K. Peucker, M. Tichenor, and W.-D. Rase, "Die Automa-
tisierung der Methode der schragen Schnittflachen," Kar-
tographische Nachrichten, vol. 22, no. 4, pp. 143-148, 1972.
[201 T. M. Oberlander, "A critical appraisal of the inclined contour
technique of surface representation," Ann. Assoc. Amer. Geog-
raphers, vol. 58, no. 4, pp. 802-813, 1968.
[211 A. H. Robinson, and R. D. Sale, Elements of Cartography, 3rd
ed. New York: Wiley, 1969, pp. 189-196.
[22 ] A. H. Robinson, and N. J. W. Thrower, "On surface representa-
tion using traces of parallel inclined planes," Ann. Assoc. Amer.
Geographers, vol. 59, no. 3. pp. 600-603, 1969.
[23] T. M. Oberlander, "Reply to Robinson-Thrower commentary,"
Ann. Assoc. Amer. Geographers, vol. 59, no. 3, pp. 603-605,
1969.
[24] Kitiro Tanaka, "The relief contour method of representing to-
pography on maps," Geographical Rev. Jap., vol. 15, nos. 9 &
10, pp. 655-671, 784-796 (Japanese) p. 797 (English abstract),
1939.
[251 -, "The relief contour method of representing topography on
maps," Geographical Rev., vol. 40, no. 3, pp. 444-456, 1950.
[26] -, "The relief contour method of representing topography on
maps," Surveying and Mapping, vol. 11, p. 27, 1951.
[27] C. Kopcke, "Ueber Reliefs und Relief-Photogramme," Civilin-
genleur, vol. 31, 1885.
[28] Y. Pauliny, "Memoire fiber eine neue Situations Plane- und
Landkartendarstellungsmethode," Streffen's Oesterreichischer
Militdrische Zeitschrift, vol. 36, pg. 177, 1895.
[29] Erwin Raisz, General Cartography. New York and London,
1938.
[301 W. Pillewizer, "Gelandedarstellung durch Reliefphotographie,"
Kartographische Nachrichten, vol. 7, p. 141, 1957.
[311 O. C. Stoessel, "Photomechanische Reliefschummerung," Nach-
richten aus dem Karten- und Vermessungswesen, vol. 1, no. 10,
pp. 53-55, 1959.
[321 A. A. Noma
and M
G
Misulia "Pro
rammi
t
h
,
.
.
g
ng
opograp
ic
, REFERENCES maps for automatic terrain model construction," Surveying and
Mapping, vol. 19, no. 3, p. 335, Sept. 1959.
[1) E. Imhof, Kartographische Gelandedarstellung. Berlin, Ger- [331 H. R. Wilkerson, "Reliefschummerung durch Photographic, von
many: W. de Gruyter & Co., 1965. Gelandemodellen," Nachrichten aus dem Karten- und Vermes-
[2] J. G. Lehmann, Darstellung einer neuen Theorie der Bezeichnung sungswesen, vol. 1, no. 10, pp. 60-62, 1959.
der schiefen Flachen im Grundriss oder der Situationzeichnung [341 H. Friedemann, "Von neuen Erfindungen: Anordnung und Ver-
der Berge. Leipzig, Germany, 1799. fahren zur Herstellung von Schummerungen fur kartographische
[3] -, Die Lehre der Situations-Zeichnung oder Anweisung zum Zwecke," Kartographische Nachrichten, vol. 12, p. 150, 1962.
richtigen Erkennen und genauen Abbilden der Erd-Oberfldche [35] P. Richarme, "L'estompage photographique," Bulletin du
in topographischen Charten und Situation-Planen. Dresden, comite francais de cartographie, vol. 17, p. 188, 1963.
Germany: Arnoldische Buch- und Kunsthandlung, 1816. [36] P. Richarme, "The photographic hill shading of maps," Survey-
[4] F. Chauvin, Die Darstellung der Berge in Karten und Planen, ing and Mapping, vol. 23, no. 1, pp. 47-59, 1963.
mit besonderer Riicksicht auf ihre Anwendbarkeit im Felde. [37] C. R. Gilman, "Photomechanical experiments in automated car-
Berlin, Germany: Nauck'sche Buchhandlung, 1852. tography," in Proc. ASM Fall Meeting (San Francisco, CA),
[5] -, Das Bergzeichnen rationell entwickelt. Berlin, Germany: 1971.
Nauck'sche Buchhandlung, 1854. [381 H. G. Lyons, The Representation of Reliefs on Maps. Ministry
[61 C. Vogel, "Die Terraindarstellung auf Landkarten mittels Schraf- of Finance, National Printing Department, Egypt, 1909.
fierung,"Petermanns Geogr. Mitt., vol. 39, p. 148, 1893. [39] L. J. Harris, Hill-shading for Relief-depiction in Topographical
[7] H. Bach, Die Theorie der Bergzeichnung in Verbindung mit Maps. London, 1959.
Geognosie. Stuttgart, Germany, 1853. [401 E. Imhof, "Gelandedarstellung in Karten grosser und mittlerer
[8] H. Wiechel, "Theorie und Darstellung der Beleuchtung von nicht Massstabe," Vortrag Natf. Ges. Zurich, January 1947.
gesetzmassig gebildeten Fl" h
a
[9]
[10]
c en mit Rucksicht auf die [41] B. Carlberg, "Schweizer Manier und wirklichkeitsnahe Karte,
Bergzeichnung," Civilingenieur, vol. 24, pp. 335-364, 1878. Probleme der Farbgebung," Kartographische Nachrichten, vol.
L. Burmester, Theorie und Darstellung gesetzmdssig gestalteter 4, pp. 8-14, 1954.
Fldchen. Leipzig, Germany, 1875. [42] G. Pohlmann, "Heutige. Methoden und Verfahren der Gelan-
Kitiro Tanaka, "A new method of topographical hill delinea- dedarstellung," Kartographische Nachrichten, vol. 8, no. 3, pp.
tion," Memo. of the College of Engineering, Kyushu Imperial 71-78, 1958.
Univ., Fukuoka, Japan., vol. 5
no. 3
pp
121-143 1930 [43] H Mi
t
"D
,
,
.
e
zner
ie Schummerung unt Ahi
,..,ernname ener naturge-
[ 11 ] -, "The orthographical relief method of representing hill mi ssen Beleuchtung," Kartographische Nachrichten, vol. 9,
features on a topographical map," Geographical J., vol. 79, no. no. 3, pp. 73-79, 1959.
3, pp. 213-219, Mar. 1932. [44] E. Imhof, "Probleme der Kartographischen Gelandedarstel-
[ 12 ] H. StJ. L. Winterbotham, "Note on Professor Kitiro's method of lung," Nachrichten aus dem Karten- und Vermessungswesen, vol.
orthographical relief," Geographical /., vol. 80, pp. 518-520, 1, no. 10, p. 9--31, 1959.
1932. [451 J. S. Keates, "Techniques of relief representation," Surveying
[131 P. Wilski, "Eine neuee Japanische Darstellung der Hohen auf and Mapping, vol. 21, no. 4, pv. 459-463, Dec. 1961.
Landkarten,"Petermanns Geogr. Mitt., vol. 80, p. 359, 1934. [46] F. Holzel, "Die Gelandeschummerung in der Krise?" Karto-
[ 14 ] Jadwiga Remiszewska, " Metoda ciec pochylych w kartografic- graphische Nachrichten, vol. 12, no. 1, pp. 17-21, Feb. 1962.
znym obrazie urzezbienia" (The method of inclined cut in the [47] E. Imhof, Ed., International Yearbook of Cartography. Lon-
cartographic presentation of the form of the land), Polish Geo- don, England: George Philip and Son, 1963.
graphical Rev., vol. 27, pp. 125-134, 1955. [48] F. Holzel G
Approved For Release 2006/03/10 : CIA-RDP8gB068% ooTu }' 0'u0- shading," Nach-
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
richten aus dem Karten- and Vermessungswesen, vol. 5, no. 5,
pp. 23-33, 1963.
[49] p Yneli, "Relief shading," Surveying and Mapping, vol. 19, p.
229, 1959.
[50] -, "Analytische Schattierung," Kartographische Nachrichten,
vol. 14, no. 4, pp. 142-148, 1965.
[511 - "Analytical hill shading," Surveying and Mapping, vol. 25,
no. 4, pp. 573-579, Dec. 1965.
[52] -, "Analytische Schattierung and Dichte," Kartographische
Nachrichten, vol. 16, no. 1, pp. 17-23, 1966.
[53] -, "Analytical hill shading and density," Surveying and Map-
ping, vol. 26, no. 2, pp. 253-259, June 1966.
[54] -, "Die Mechanisierung der Analytischen Schattierung,"
Kartographische Nachrichten, vol. 16, no. 3, pp. 103-107, 1966.
[551 -, "The mechanisation of analytical hill shading," Cartographic
J., vol. 4, no. 2, Dec. 1967.
[561 - "Die Richtung des Lichtes bei analytischer Schattierung,"
Kartographische Nachrichten, vol. 17, no. 2, pp. 37-44, 1967.
[571 - "An experimental electronic system for converting con-
tours into hill-shaded relief," International Yearbook of Car-
tography, vol. 11, pp-'I I 1- 114, 1971.
[581 T. J. Blachut, Z. Marsik, and D. Makow, "Relief shading pro-
cess," Canadian Patent No. 051-R14, 1969.
[591 Z. Marsik, "Automatic Relief Shading," Photogrammetrica, vol.
27, no. 2, pp. 57-'70, 1971.
[601 T. K. Peucker, "Computer cartography," Association of Amer-
ican Geographers, Washington, DC, Commission on College
Geography, Resource Paper No. 17, pp. 41-54, 1972.
[611 T. K. Peucker and D. Cochrane, "Die Automation der Relief-
darstellung-Theorie and Praxis," International Yearbook of
Cartography, vol. 14, pp. 128-139, 1974.
[62] T. K. Peucker, M. Tichenor, and W: D. Rase, "The computer
version of three relief representations," in Display and Analysis
of Spatial Data, J. C. Davis, and M. McCullagh, Eds. New York:
Wiley, 1974.
[63] M. Eckert, Die Kartenwissenschaft: Forschung and Grundlagen
zu einer Kartographie als Wissenschaft - Vol. 1. Berlin, Ger-
many: Walter de Gruyt, 1962.
[641 L. D. Carmichael, "Experiments in relief portrayal," Carto-
graphic J., vol. 1, pp. 11-17, 1964.
[651 M. Eckert, Die Kartenwissenschaft: Forschung and Grundlagen
zu eener Kartographie als Wissenschaft - Vol. 2. Berlin, Ger-
many: Walter de Gruyt, 1965.
[661 M. S. Monmonier, "The production of shaded maps on the digi-
tal computer," Professional Cartographer, vol. 17, no. 5, pp. 13-
14, Sept. 1965.
[671 P. K. Koldayev, "Plastic colour and shadow relief representa-
tion," Academy of Sciences of the USSR, Council of Soviet
Cartographers, Moscow, 1967.
[68] B. F. Sprunt, "Computer-generated halftone images from digital
terrain models," M.Sc. dissertation, Dep. Mathematics, Univ.
Southampton, Southampton, England. 1969.
[69] B. Anda, "Automatic hill-shading using an automatic flatbed
drafting machine with a standard photohead," ITC Journal
(Enschede, The Netherlands). no. 2, pp. 212-216, 1974.
[70] R. M. Batson, E. Edwards, and E. M. Eliason, "Computer gener-
ated relief images," J. Res., U. S. Geological Survey, vol. 3, no.
4, pp. 401-408, July-Aug. 1975.
[711 K. Brassel, "Modelle and Versuche zur automatischen Schrdg-
lichtschattierung," Ph.D. dissertation, Geography Dep., Univ.
Zurich, Klosters, Switzerland, 1973.
[72] - "Ein- and mehrfarbige Printerdarstellungen," Karto-
graphischeNachrichten, no. 5, pp. 177-183, 1973.
[731 -, "Ein Modell zur Automatischen Schraglichtschattierung,"
International Yearbook of Cartography, pp. 66-77, 1974.
[74] - "A Model for Automated Hill Shading," Amer. Cartog-
rapher, vol. 1, no. 1, pp. 15-27, Apr. 1974.
[75] W. Blascke, "Le Modele Digital M.I.T.," Societe Francaise de
Photogrammetrie, Bull. 27, pp. 37-40, July 1967.
[761 B. W. Boehm, "Tabular representations of multi-variate func-
tions-with applications to topographic modelling," in ACM,
22nd National Conf. Proc., pp. 403-415, 1967.
[77] W. Aumen, "A new map form: Numbers," International Year-
book of Cartography, vol. 10, pp. 80-84, 1970.
[78] M. W. Grist, "Digital ground models: An account of recent re-
search," Photogrammetric Rec., vol. 70, no. 4, pp. 424-441,
Oct. 1972.
[79] F. Silar, "Das digitale Gelandemodell-Theorie and Praxis,"
Vermessungstechnik, vol. 20, no. 9, pp. 327-329, 1972.
[80] K. Torlegard, "Digital terrain models-General survey and
Swedish experiences," Bildmessung and Luftbildwesen, vol. 40,
no. 1, pp. 21-30, 1972.
[81] American Society of Photogrammetry, Proc. Digital Terrain
Models (DTM) Symp. (St. Louis, MO, May 9-11, 1978).
[82] Wild Heerburg and Raytheon, "B8 Stereomat Automated Plot-
ter," company sales literature, 1968.
Approved For Release 2006/03/10 :
[ 831 S. Bertram, "The UNIMACE and the automatic photomapper,"
Photogrammetric Eng., vol. 35, pp. 569-576, 1969.
[841 R. H. Seymour and A. E. Whiteside, "A new computer-assisted
stereocomparator," Bendix Tech. J., pp. 1-5, Spring 1972.
[851 B. G. Crawley, "Gestalt contours," Canadian Surveyor, vol. 28,
no. 3, pp. 237-246, Sept. 1974.
[86] Bendix Research Laboratories, "AS-11B-X automatic stereo
mapper," RADC-TR-76-100, Rome Air Development Center,
Griffiss A F B, N Y, Apr. 1976.
[87] D. J. Panton, "Digital stereo mapping," Countermeasures, p. 12,
May 1976.
[88] W. Loscher, "Some aspects of orthophoto technology," Photo-
grammetric Rec., vol. 6, no. 30, pp. 419-432, 1967.
[891 T. J. Blachut, "Further extension of the orthophoto technique,"
Canadian Surveyor, vol. 22, no. 1, pp. 206-220, 1968.
[90] T. J. Blachut and M. C. Van Wijk, "3-D information from ortho-
photos," Photogrammetric Eng., pp. 365-376, April 1970.
[91] T. A. Hughes, A. R. Shope, and F. S. Baxter, "USGS automatic
orthophoto system," Photogrammetric Eng., vol. 37, pp. 1055-
1062, 1971.
[921 A. Beyer, "Zur Erfassung flachen Gelandes durch willkiirlich
verteilte Hohenpunkte," Vermessungstechnik, vol. 20, no. 6,
pp. 204-207, 1972.
[93] T. K. Peucker and N. Chrisman, "Cartographic data structures,"
American Cartographer, vol. 2, no. 1, pp. 55-69, Apr. 1.975.
[94] E. Keppel, "Approximating complex surfaces by triangulation
of contour lines," IBM J. Res. Develop., Jan. 1975.
[95] T. K. Peucker, R. J. Fowler, J. J. Little, and D. M. Mark, "Digital
representation of three-dimensional surfaces by triangulated
irregular networks (TIN)," Tech. Rep. 10 (Revised), Dep.
Geography, Simon Fraser University, Barnaby, BC, Canada,
1976.
[96] A. K. Lobeck, Block Diagrams. New York: Wiley, 1924.
Reprinted by Emerson-Trussell, Amherst, MA, 1958.
[97] M. Schuster, Das Geographische and Geologische Blockbild.
Berlin, Germany: Akademie Verlag, 1954.
[981 D. A. Goosen, "Blockdiagrams," ITC Information (Delft, The
Netherlands), no. 3, Spring 1962.
[991 G. F. Jenks and D. A. Brown, "Three-dimensional map construc-
tion," Science, vol. 154, no. 3750, pp. 857-864, 18 Nov. 1966.
[100] B. Kubert, J. Szabo, and S. Giulieri, "The perspective represen-
tation of functions of two variables," J. A.C.M., vol. 15, no. 2,
pp. 193-204, Apr. 1968.
[1011 D. Douglas, "Viewblok: A computer program for constructing
perspective view block diagrams," Rev. Geographie de Montreal,
vol. 26, pp. 102-104, 1971.
[1021 T. J. Wright, "A two-space solution to the hidden line problem
for plotting functions of two variables," IEEE Trans. Comput.,
vol. C-22, pp. 28-33, Jan. 1973.
[1031 C. Wylie, G. Romney, and D. Evans, "Half-tone perspective
drawings by computer, " in Proc. Fall Joint Computer Conf.,
pp. 49-58, 1967.
[1041 A. Appel, "Some Techniques for shaded machine rendering of
solids," in Proc. Spring Joint Computer Conf., pp. 37--45, 1968.
[105] J. E. Warnock, "A hidden-surface algorithm for computer
generated half-tone pictures," TR 4-15, Dep. Computer Science,
Univ. Utah, Salt Lake City, 1969.
[106] G. S. Watkins, "A real-time visible surface algorithm," Report
UTEC-CSC-70-101, Dep. Computer Science, Univ. Utah, Salt
Lake City, June 1970.
[107] W. J. Bouknight, "A procedure for generation of three-dimen-
sional half toned computer graphics presentations," Commun.
A.C.M., vol. 13, no. 9, pp. 527-536, Sept. 1970.
[108] R. A. Goldstein and R. Nagel, "3-D visual simulation," Simula-
tion, vol. 16, pp. 25-31, 1971.
[1091 H. Gouraud, "Computer display of curved surfaces," IEEE
Trans. Comput., vol. C-20, pp. 623-629, June 1971.
[110] M. E. Newell, R. G. Newell, and T. L. Sancha, "A new approach
to the shaded picture problem," in Proc. ACM National Conf.
(Boston, MA), vol. 1, pp. 443-450, 1973.
[1111 J. Staudhammer and D. J. Odgen "Computer graphics for half-
tone three-dimensional object images," Computers and Graphics,
vol. 1, no. 1, pp. 109-114, 1975.
[1121 E. A. Catmull, "Computer display of curved surfaces," in Proc.
IEEE Conf. Computer Graphics, Pattern Recognition and Data
Structures (Los Angeles, CA, May 1975), IEEE Cat. No.
75CH0981-1C, pp. 11-17, 1975.
[1131 Phong Bui-Tuong, "Illumination for computer-generated images,"
Commun. A.C.M., vol. 18, no. 6, pp. 311-317, June 1.975.
[114] J. F. Blinn and M. E. Newell, "Texture and reflection in com-
puter generated images," Commun. A.C.M., vol. 19, no. 10, pp.
542-547, Oct. 1976.
[115 ] J. F. Blinn, "Models of light reflection for computer synthesized
pictures," in SIGGRAPH '77,Proc. 4th Conf. Computer Graphics
and Interactive Techniques, A.C.M., pp. 192-198, 1977.
[116] - "A scan line algorithm for displaying parametrically de-
CIA-RDP88B00553R000100280003-6
 HORN: HILL SHADING AND THd For Release 2006/03/10 : CIA-RDP88BOO553R000100280003-6
fined surfaces," in SIGGRAPH '78, Proc. 5th Conf. on Com-
puter Graphics and Interactive Techniques, A.C.M., 1978.
[117] J. H. Lambert, Pholometria sive de mensura de gratibus luminis,
colorum et umbrae, Eberhard Klett, Augsburg, 1760. Trans-
lated by W. Engelman "Lambert's Photometrie," Ostwald's
Klassiker der exacten Wissenschaften, no. 31-33, Leipzig, Ger-
many, 1892.
[118] l'Abbe de Lacaille, Traite d'optique sur la gradation de la
lumiere, (Ouvrage posthume de M. Bouguer), L. F. Delatour, A
Paris, 1760. Translated by W. E. K. Middleton, Optical Treatise
on the Gradation of Light. Toronto, Ont., Canada: University
of Toronto Press, 1961.
[119] E. Lommel, "Ueber Fluorescenz," Ann. Phys. (Leipzig, Ger-
many), vol. 10, pp. 449-472, 1880.
[1201 H. Seeliger, "Die Photometrie von diffus reflektierenden
Flachen," S. B. Bayer. Akad. Wiss., vol. 18, p. 20, 1888.
[1211 A. Markov, "Les particularites dans le reflexion de la lumiere
par la surface de la lune," Astronomische Nachrichten, vol. 221,
pp. 65-78, 1924.
[1221 E. Schonberg, "Untersuchungen zur Theorie der Beleuchtung
des Mondes auf Grund photometrischer Messungen," Acta Soc.
Sci. Fennicae, vol. 50, pp. 1-70, 1925.
[123] V. G. Fesenkov, "Photometric investigations of the lunar sur-
face," Astronomicheskii Zhurnal, vol. 5, pp. 219-234, 1929.
Translated by Redstone Scientific Information Center. Apr.
1968.
[124] M. Minnaert, "The reciprocity principle in lunar photometry,"
Astrophys. J., vol. 93, pp. 403-410, 1941.
[125] V. A. Fedoretz, "Photographic photometry of the lunar surface,"
Publ. Kharkov Obs., vol. 2, pp. 49-172, 1952.
[1261 M. Minnaert, "Photometry of the moon," in Planets and Satel-
lites, G. P. Kuiper, and B. M. Middlehurst, Eds. Chicago, IL:
Univ. Chicago Press, 1961, vol. 3, ch. 6, pp. 213-248.
[127] V. Fesenkov, "Photometry of the Moon," in Physics and As-
tronomy of the Moon, Z. Kopal, Ed. New York: Academic
Press, 1962, pp. 99-:130.
[1281 B. W. Hapke, "A Theoretical Photometric Function for the
Lunar Surface," J. Geophys. Res., vol. 68, no. 15, pp. 4571-
4586, Aug. 1963.
[129] B. Hapke and H. Van Horn, "Photometric studies of complex
surfaces, with applications to the moon," J. Geogphys. Res.,
vol. 68, no. 15, pp. 4545-4570, Aug. 1963.
[130] B. Hapke, "An improved theoretical lunar photometric func-
tion," Astronomical J., vol. 71, no. 5, pp. 333-339, June 1966.
[131] W. E. Middleton and A. G. Mungall, "The luminous directional
reflectance of snow," J. Opt. Soc. Amer., vol. 42, no. 3, pp.
572-579, 1952.
[1321 M. Planck, The Theory of Heat Radiation. New York: Dover,
1959.
[133] P. Beckmann and A. Spizzichnio, The Scattering of Electromag-
netic Waves from Rough Surfaces. New York: Pergamon Press,
1963.
[1341 K. E. Torrance, E. M. Sparrow, and R. C. Birkebak, "Polariza-
tion, directional distribution, and off-specular peak phenomena
in light reflected from roughened surfaces," J. Opt. Soc. Amer.,
vol. 56, no. 7, pp. 916-925, July 1966.
[135] K. E. Torrance and E. M. Sparrow, "Theory for off-specular
reflection from roughened surfaces," J. Opt. Soc. Amer., vol.
57, no. 9, pp. 1105-1114, Sept. 1967.
[1361 T. S. Trowbridge and K. P. Reitz, "Average irregularity represen-
tation of a rough surface for ray reflection," J. Opt. Soc. Amer.,
vol. 65, no. 5, pp. 531-536, May 1975.
[137] B. K. P. Horn, "Determining shape from shading," in The Psy-
chology of Computer Vision, P. H. Winston, Ed. New York:
McGraw-Hill, 1975, ch. 4.
[138] - "Understanding image intensities," Artificial Intelligence,
vol. 8, no. 11, pp. 201-231, 1977.
[1391 -, and B. L. Bachman, "Using synthetic images to register
real images with surface models," Commun. A.C.M., vol. 21, no.
11, pp. 914-924, Nov. 1978.
[140] and R. W. Sjoberg, "Calculating the reflectance map," Appl.
Opt., vol. 18, no. 11, pp. 1770-1779, June 1979.
[141] R. J. Woodham, "Photometric stereo: A reflectance map tech-
nique for determining surface orientation from image intensity,"
Image Understanding Systems and Industrial Applications, in
Proc. S.P.I.E., Aug. 1978, vol. 155.
[142] F. E. Nicodemus, J. C. Richmond, and J. J. Hsia, I. W. Ginsberg,
T. Limperis, "Geometrical considerations and nomenclature for
reflectance," NBS Monograph 160, National Bureau of Stan-
dards, U.S. Dep. Commerce, Washington, DC, Oct. 1977.
[143] F. E. Nicodemus, Ed., "Self-Study manual on optical radiation
measurements," NBS Technical Notes 910-1, 910-2, & 910-3,
National Bureau of Standards, U.S. Dep. Commerce, Washing-
ton, DC, 1976, 1977, & 1978.
[1441 H. C. Babcock, "Evaluation of a stereocompilation digitizer, "
in Congress on Surveying and Mapping, 30th Annual Meeting,
1970, pp. 338-347.
[1451 E. J. McCartney, Optics of the Atmosphere: Scattering by
Molecules and Particles. New York: Wiley, 1976.
[1461 S. D. Conte and C. de Boor, Elementary Numerical Analysis,
New York: McGraw-Hill, 1972.
[147] R. W. Hamming, Numerical Methods for Scientists and Engineers.
New York: McGraw-Hill, 1962.
[148] R. D. Richtrneyer and K. W. Morton, Difference Methods for
Initial-Value Problems. New York: Wiley, 1967, pp. 136-143.
[1491 F. B. Hildebrand, Introduction to Numerical Analysis. New
York: McGraw-Hill, 1956, 1974.
[150] P. M. Bridge and J. L. Inge, "Shaded relief of Mars," Atlas of
Mars, MH 25 M IR, JPL Contract WO-8122, USGS, U.S. Dep.
Interior, 1972.
[1511 T. Rindfleisch, "Photometric method for lunar topography,"
PhotogrammetricEng., vol. 32, pp. 262-276, 1966.
[1521 W. Bantel, "Der Reproduktionsweg vom einfarbigen Reliefori-
ginal zur mehrfarbigen Reliefkarte," International Yearbook of
Cartography, vol. 13, pp. 134-136, 1973.
[153] A. DeLucis, "The effect of shaded relief terrain representation
on map information accessibility," in American Congress on
Surveying and Mapping, 31st Annual Meeting (Washington, DC),
pp. 641-657, 1971.
[1541 J. Neumann, "Begriffsgeschichte and Definition des Begriffes
'Kartographische Generalisierung'," International Yearbook of
Cartography, vol. 13, pp. 59-67, 1973.
[155] F. Topfer, Ed. Kartographische Generalisierung. Gotha,
Leipzig, Germany: Geographisch-Kartographische Anstalt, 1974.
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
Acousto-Optics-A Review of Fundamentals
Abstract-The paper first reviews the historical development of
acousto-optics and its applications. Following this, a heuristic explana-
tion of acousto-optic effects is presented, with the emphasis on the
plane wave model of interaction. Finally, there is a discussion of some
basic configurations of relevance to signal processing.
I. INTRODUCTION
A OUSTO-OPTICS deals with the interaction of sound
and light. The existence of such an interaction was
predicted by Brillouin [ 1 ] in 1922. Experimental
verification followed in 1932, by Lucas and Biquard [2] in
France, and Debye and Sears [31 in the U.S. Brillouin's
original theory predicted a phenomenon closely analagous
to X-ray diffraction in crystals. In the latter, the atomic
planes cause multiple reflections of an incident electromag-
netic plane wave. These reflections interfere constructively
for certain critical angles of incidence, to cause enhanced
overall reflection (also called diffraction or scattering). In
acoustic diffraction, the role of the atomic planes is assumed
by planes of compression and rarefaction, induced by ultra-
sonic waves with frequencies between 1 MHz and 2 GHz.
As a result, similar diffraction effects occur as in crystals.
The similarity between atomic planes and sound wave
fronts should, however, not be carried too far. For one thing,
atomic planes are sharply defined in location and regular in
structure. Sound waves, on the other hand, are essentially
sinusoidal and, when limited in the transverse direction, spread
as they propagate, thereby giving rise to very complex density
distributions and wavefronts. To a first approximation how-
ever, the analogy is very useful.
The fact that the sound wavefronts move causes the dif-
fracted light to be Doppler shifted. Brillouin predicted a basic
Doppler shift equal to the sound frequency. This phenomenon
of frequency shifting became of importance only in recent
times and forms, in fact, the basis of heterodyning techniques
in modern signal processing applications.
Long before such sophisticated techniques became practical
through the invention of the laser, the application of acousto-
optics as such to (parallel) signal processing had become evi-
dent, and experiments were already in progress using arc lamps
and other conventional light,sources. The pioneer in this field
is probably Okolicsanyi who, in the early 1930's used acousto-
optic techniques in the Scophony television projector [4].
In the course of this work Okolicsanyi made a thorough investi-
gation of possible sound cell applications (for instance time
compression) which he published in 1937 [ 5 ] .
Manuscript received August 14, 1980.
The author is with the Department of Electrical and Computer Engi-
neering, University of Iowa, Iowa City, IA 52242.
It was not until the early 1960's that the subject was taken
up again when signal processors were investigated by Rosenthal
[6] , Liben [7] , Slobodin [8] and Arm et al. [9] . By the end
of the 1960's lasers were beginning to be used routinely which
in turn led to the development of coherent heterodyne detec-
tion techniques by King et al. [ 10] and Whitman et al. [ 11 ] .
Thus, almost 50 years after Brillouin's analysis, all the phe-
nomena he predicted had been put to use.
During the past decade even more refined techniques of
signal processing were developed, in particular ones involving
crossed sound cells (two-dimensional processing) and time
integration. A detailed discussion of this development is
outside the scope of this introduction. The interested reader
will find abundant information in a number of review articles
[ 121 -[171 , as well as in subsequent papers in this issue.
In what follows, the author will first give a brief, heuristic
explanation of acousto-optic effects, emphasizing the plane
wave interaction model but tracing, wherever practical, the
connection with other, perhaps more familiar models. Next,
some configurations that form the building blocks of most
signal processors will be reviewed. This paper will be limited
to the basic elements and avoid system details such as lens
configurations, dynamic range considerations, intermodulation
effects, etc. A discussion of such details is not appropriate
for this brief review; they can best be learned from the papers
that follow.
II. HEURISTIC THEORY
As pointed out in the Introduction, Brillouin's analysis of
acousto-optic interactions predicted an effect, now called
acoustic Bragg diffraction, similar to X-ray diffraction in
crystals. The principle of acoustic Bragg diffraction is illus-
trated in Fig. 1 which shows a beam of light directed at a
high frequency sound wave. At certain critical angles of
incidence ?aB, the incident beam generates a new one whose
direction differs by 2aB. The angle aB is called the Bragg
angle and is given by
A
sin aB = 2A
where A is the wavelength of light in the acoustic medium and
A is the acoustic wavelength.
The value given by (1) refers to the angle observed inside the
medium of sound propagation. Outside the medium, Snell's
law changes the value of aB to aB
nX X?
I sin aB = =
A A
Q018 9 z866/3?PUq0CIAORDPc~8B00553R000100280003-6
Approved For Release
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
KORPEL: ACOUSTO-OPTIC-A REVIEW OF FUNDAMENTALS
Fig. 3. Diffraction in the Debye-Sears region showing many generated
orders.
Fig. 1. Diffraction in the Bragg region showing downshifted (top) and
upshifted (bottom) interaction.
kin
kin
Fig. 2. Wave vector diagrams illustrating diffraction in the Bragg region
for downshifted (top) and upshifted (bottom) interaction.
where A? denotes the vacuum wavelength. In the drawings,
this refraction effect is not shown; for simplicity beams are
shown propagating in the direction appropriate to the medium.
If the wavefronts of the sound move away from the incident
light (Fig. 1, top), they also move away from the diffracted
beam which is called the -1 order in this case. It is easily
shown that this results in a negative Doppler shift, i.e., the
diffracted light is downshifted in frequency by the sound
frequency co. The case of upshifted Bragg diffraction is also
shown in Fig. 1 (bottom); the diffracted light is now called
the +1 order.
The essential properties of Bragg diffraction can be explained
in many different ways [181-[221. This author's preference
lies with the model in which the interaction is thought of as
a collision between photons and phonons. The momenta of
the interacting particles are given by hk and hK for the pho-
tons and phonons respectively, where h = h/27r, h is Planck's
constant, k and K are the propagation vectors of light and
sound. Conservation of momentum is then illustrated by the
diagrams of Fig. 2, from which the Bragg angle conditions can
be derived easily [ 181
sin aB = 2 K/k = 2 X/A.
Fig. 4. Multiple scattering in the Debye-Sears region (bottom) made
possible by the presence of many plane waves of sound in the radia-
tion pattern of the transducer (top).
If the width of the transducer L is decreased, as in Fig. 3, the
column of sound in the medium will look less and less like a
plane wave. More accurately, the angular plane wave spectrum
of the soundfield (i.e., the radiation pattern) broadens so as
to make available additional plane waves of sound for a wider
range of interaction. The basic idea is shown in Fig. 4. The
first effect of shortening L is that the up-and-down-shifted
orders kl and k_1 are now both generated due to the simul-
taneous presence of soundwaves K1 and K_1 with amplitudes
S1 and S-1. If L is decreased further, additional soundwaves
K2 and K_2 appear, with amplitudes S2 and S-2, which inter-
act with already scattered lightwaves k1 and k_1 to generate
the second orders k2 and k_2. If we continue to decrease L,
more and more orders are generated, as shown in Fig. 3. (Note
that the nth order is shifted in frequency by an amount nw.)
The extent to which this multiple scattering process can oper-
Note that, classically, the diagrams of Fig. 2 represent
phase synchronous interaction, a concept which was already
implicit in Brillouin's early paper.
Frequency up- or downshift follow from the same diagrams ate is determined by the ratio of the angular width A/L, of
by considering energy conservation in terms of photon-and the sound radiation pattern to the separation X/A between
phonon energies htwi and 1iw. orders. Customarily a parameter Q, inversely proportional
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
K2L AL
Q= k =27r 1 .
Thus the condition for Bragg diffraction (the Bragg region)
is stated as
Q >> 1 (5)
while the opposite condition denotes the multiple scattering
regime where many orders may be generated. This latter
regime is called the Raman-Nath [241, Lucas-Biquard, or
Debye-Sears region after early investigators.
Conventionally, Debye-Sears diffraction is analyzed by
using a moving phase grating model as was first done by
Raman and Nath [241. This approach can be shown to lead
to the same results as are obtained by using multiple scatter-
ing calculations [25], [26]. The amplitude of the diffracted
orders is given by
An = (-j)11EinJn(kCISIL/2). (6)
where Ei?c is the amplitude of the incident light, ISM represents
the peak condensation or strain in the soundfield S(x, t) _
ISM cos (wt - kx), and C is a material constant in terms of
which the refractive index variation is given by
An(x, t) = z nCS(x, t).
With (7) it is easily shown that (6) may also be written as
An = (-j)nEincJn(&) (8)
where a = k? &L denotes the peak phase shift of the light
due to the peak refractive index variation On, and k? = 21r/Au.
The constant C is related to the so called strain-optic coeffi-
cient p by
C= -n2p. (9)
More generally [ 18] , n is a tensor and strictly speaking, equa-
tion (9) applies only to liquids.
Fig. 5 shows the dependence of the diffracted orders on a,
leaving out the (-j)n form. Quite often a modest sound pres-
sure is used (weak interaction, a > 1 /v0. The third term, for AT sufficiently
large, is the desired correlation R12 (?) as a function of x/ V.
D. One-Cell System with Electronically Inserted Reference
We note briefly a method suggested by Kellman for increas-
ing the versatility and convenience of time-integration acousto-
optic processors. The system employed is a simple modifica-
tion of that of Fig. 11, in which input signal v1(t) is replaced
by signal
v1(t) = bl(t) cos [27r(vo + ve)t +(31(t)] +A cos 27rvot (57)
the phase of the zeroth order by 90? and Tblocks the frequency- obtained from v1(t) by shifting the carrier frequency by an
downshifted diffra tproved For Release t 00t6/03/10 e : CIA RDP?8B00553 0001002g80003-6 same frequency.
+u2(t)[1 +(4)bi(t - x/V)} at
 Approved For Release 2006/03/10 : CIA-RDP88BOO553ROO0100280003-6
RHODES: AO SIGNAL, PROCESSING
With v (t) as the input, the acousto-optic cell transmittance,
assuming Bragg regime operation, is given by
t(x,t)=I+j[A+(2)u"1(t-x/V)] exp[j21rv,(t-x/V)]. (58)
The cell is illuminated as before by the modulated diode and
imaged, this time with a zeroth-order stop. The resultant
detector plane intensity is given by
Id(x, t) = [B+u2(t)JJA+(Z)v1(t- x/V)12
= [B + u2(t)] [A2 + (4)bi(t - x/V) +Av1 (t - x/V)].
(59)
Assuming AT>> I/vo, the associated time-integrated intensity
is approximated by
EAT(x) = A2BAT + A J u1 (t - x/V) u2 (t) dt (60)
DT
which is similar to the results for the previous system (which
corresponds to setting A equal to unity), but now the ratio of
the correlation term to the bias term is proportional to A-1
and can therefore be adjusted for optimum system perfor-
mance [40].
E. Linear Intensity Modulation Method
As our final example, we describe a time-integration cor-
relator/convolver that represents a significant departure from
the other systems-both time-integrating and space-integrating-
that we have considered to this point. The scheme, described
by Sprague and Koliopoulos [33], relies on characteristics of
acousto-optic cells operating in the Bragg regime at high dif-
fraction efficiencies. For an acousto-optic modulator with
sinusoidal driving signal v(t) = b cos 27rvot, the intensity of the
diffracted wave Idiff for Bragg regime operation is given by
Idiff = Iillum sin2 (Kb), where Iillum is the intensity of the
incident light and K is a constant. For low diffraction effi-
ciency, Kb -l ms for a large array) and
the low dynamic range (typically 20-30 dB). Recent work
aimed at improving the device yield and quantum efficiency of
large frame transfer imagers has resulted in a device which uses
a virtual phase clocking [ 13 1. The device has a simpler gate
construction than other types of area CCD arrays and a larger
effective collection area per pel.
The two basic categories of linear array photodetectors
useful for AO signal processing are the photodiode and the
metal-oxide-semiconductor (MOS) depletion mode sensor.
The mechanisms for accessing photodiode line arrays are
either self-scanned MOS shift registers, random address MOS
switches, or CCD parallel-in/serial-out shift registers. MOS
depletion mode sensors are usually coupled directly to CCD
shift registers by direct charge injection.
An example of a commercially available randomly address-
able photodiode linear array is the Reticon CP1006 device
[ 141. The device consists of two interdigitated rows of 256
Approved
PRI = 240U3
Fig. 3. Signal output from the integrated optical real-time spectrum
analyzer when two simultaneous RF pulses are applied to its input.
photodiodes with a pitch of 50 pm (so called 2P configura-
tion). Its dynamic range is 20 dB for noncoherent detection
processors. The array can be sequentially scanned or any pel
can be read at the strobe rate of 1 MHz [151. The device
has been used in bulk AO RF spectrum surveillance and identifi-
cation systems [ 16) . The random access feature enables pre-
determined frequency slots to be addressed without the time
delay of sequentially scanning all pels. An example of a com-
mercially available photodiode array with CCD PI/SO shift
registers is the Reticon CP 1023 device [171. The device
has 256 pels on an 18-?m pitch. The device has four indepen-
dent outputs of 64 pels each, capable of a 6-MHz data rate
(10.6-?s access time). Its specified dynamic range is 30 dB
for noncoherent detection systems.
A self-scanned photodiode line array specifically designed
for use with integrated and bulk optical real time spectrum
analyzers has been built and tested both separately [181,
[ 191 and as an integral part of a working integrated optical
spectrum analyzer [201. The sensor has 140 pels on a 1.2-pm
pitch, OP architecture (pels side by side) and a single photo-
diode for zeroth-order beam dumping. The array is composed
of seven groups of 20 pels each. Each group is addressed in
parallel by a ten-stage dynamic PMOS shift register. Each
group has two electrometer output circuits. One circuit out-
puts odd pels while the other circuit outputs even pels during
each clock cycle. Using this scheme, fourteen parallel samples
are available every clock cycle. The device operates at a clock
For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
Fig. 4. Photomicrograph of the 5050 silicon photodetector designed
specifically for AO spectrum analyzers,
rate of 5 MHz resulting in an access (or integration) time of 2
ps. Large charge handling capacity is achieved by having large
parasitic and photodiode junction capacitance. The measured
dynamic range of the device is 41 dB for 2-ps access time. The
pel to pel electrical crosstalk is limited by isolating the N+ photo-
diode elements in a boron tub made by ion implantation (a
tub is a region of opposite carrier concentration to that of the
substrate material in which device structures are fabricated.)
Holdover crosstalk is reduced by using a common source reset
switch. A cross section of the device is shown in Fig. 2. The
tub junction is 6 pm deep. This type of tub structure limits
crosstalk by cutting off the diffusion of both deeply absorbed
photons and thermally generated carriers, at the expense of
reduced quantum efficiency. An oscillograph of the output of
the detector when two simultaneous RF pulses of 1.0-ps pulse
duration, 4.17-kHz repetition rate, and 680- and 700-MHz
frequencies are applied to the integrated optical spectrum
analyzer is shown in Fig. 3 [ 211. A photomicrograph of the
sensor is shown in Fig. 4.
A MOS depletion mode linear array photodetector which is
accessed by PI/SO CCD shift registers has been described (221.
The device is intended for use with an integrated optical spec-
trum analyzer. It consists of four groups of 25 pels on a pitch
of 8?m with a 2P architecture. It is designed to be readout in
four parallel channels, each channel addressing 25 pels or in
two parallel channels each channel addressing 50 pels.
IV. DEVICE PHYSICS
The fundamental process in a photodetector is the genera-
tion of free charge carriers by a direct or indirect band-gap
transition due to absorption of photons in a semiconductor
material. The three classes of photodetectors are photocon-
ductors, depletion-layer devices, and avalanche photodiodes.
For the purpose of line or area photodetector arrays for AO
signal processing, only depletion-layer devices fabricated using
silicon integrated circuit technology have been developed ex-
tensively (although a GaAs photodiode line array for AO has
been made [23] and is discussed in Section VI-B.) These de-
vices form their photocharge collecting depletion regions
either by reverse biased p-n junctions or by electric fields
(i.e., CCD's). Excellent review papers and books exist in the
literature describing the device physics of these photodetectors
[ 101, [ 241. Therefore, only those physical processes that are
of critical importance for AO signal processing photodetectors
are presented here in detail. They include quantum efficiency,
responsivity, crosstalk, charge handling capacity, data rate,
optical to electrical transfer function, noise sources, dynamic
range, and sensitivity. As an example, a generalized expression
which shows how these various sensor parameters interrelate
to optical parameters has been derived for an AO space-in-
tegrating spectrum analyzer [4] and is given here in detail
i
Charge
This expression describes the parameters associated with a given
pel of the sensor. The term S is the efficiency of the opto-
electric interaction expressed as the fraction of light deflected
per applied watt per instantaneous system bandwidth, A f Ps
exposure
responsivity time
[RD X]
deflected transfer
optical function
power density
8P8PoI
(7rHP/4)Af
dynamic
range
QT QT(NES).
QT(NES)~
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
 I
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
is the RF signal power of a single tone applied to the AO de-
vice which is required to cause saturation of a given pet. P0
is the power of the optical source. The term i is the transmis-
sion loss in the optics including that due to optical apodiza-
tion. The term (7rHP/4) is the area of a resolvable spot where
H is the height of the pet element and P is its pitch. This term
implicitly includes the optical modulation transfer function
(MTF) of the system. The criteria for choosing H and P for a
specified optical resolution are considered in Section V-A for
bulk systems and V-B for integrated optical systems. The term
MD is the product of the modulus of the electrical transfer
function and the sensor MTF. The electrical transfer function
is a pleasure of the fidelity of the sensor's signal"processorread-
out with respect to frequency (data rate), amplitude, and phase
of the sampled data. The combined MTF is a measure of the sen-
sor electrical crosstalk performance and its optical-to-electrical
transfer characteristic.
A. Quantum Efficiency
Quantum efficiency is defined as the number of signal carriers
collected by the detector divided by the number of photons
that have impinged upon the detector. The factors which in-
fluence quantum efficiency are the optical transmission of
the insulating layers that cover the detectors, junction depths
for photodiode arrays, depletion layer thicknesses, bulk life-
time, surface recombination velocity, optical absorption
coefficient, and pel-to-pel isolation techniques. The absolute
quantum efficiency of a detector is determined by the collec-
tion efficiency of the individual regions of the semiconductor
that make up the device multiplied by the transmission ef-
ficiency of the external optical layers and isolation techniques.
A representative cross section of a silicon photodiode sensor
array for AO signal processing is shown in Fig. 2 for example.
The first semiconductor region is an N+ layer. In this region,
several mechanisms contribute to quantum efficiency. Surface
recombination of minority carriers, characterized by the surface
recombination velocity, causes a loss in efficiency since carriers
created by absorption of photons at or very near the surface
can be lost by hole-electron recombination. This loss is im-
portant for short wavelength photons which have a high ab-
sorption coefficient. Similarly, recombination in the bulk
volume of the N+ layer causes loss of carriers. Bulk recombi-
nationloss is characterized by the minority carrier lifetime [ 25 1.
However, in the N+ region a gradient of ionized donar atoms
gives rise to an electric field that aids in separating electron-
hole pairs once created. The second region is the depletion
layer which spreads from the N+ junction into the P-type layer.
Carriers created within this region experience a high electric
field and are almost immediately collected. The carriers in
the undepleted P tub experience no built-in electric field if
this region is considered to be of uniform doping. Only dif-
fusion and recombination affect quantum efficiency in this
region. The thickness of this layer is important because it
limits pel-to-pel electrical crosstalk. Photons which are ab-
sorbed below the tub-substrate junction are not collected by
any photodiode and therefore cannot contribute to crosstalk.
However, the price paid for this crosstalk improvement is a
loss in quantum efficiency. The equations which describe the
quantum efficiency in these three regions are quite complex.
Gary [26] and more recently Van de Wiele [27] describe both
numerical and simplified analytical expressions for the quantum
efficiency in these regions.
The quantum efficiency of the photodetector is also in-
fluenced by the technique used to isolate or define pels. In
the example shown in Fig. 2, adjacent pet isolation is provided
by the potential maximum formed between N+ diffusions. In
this example, all photons absorbed between pels are collected.
However, two other techniques are possible. The first is a
heavily doped channel stop diffusion or ion implant. This type
of diffusion is necessary in MOS-CCD type structures to define
pet elements but results in substantial recombination of ab-
sorbed photogenerated charge in this region. Another tech-
nique is to isolate pels by the preferential etching of narrow
grooves in (110) silicon [ 28 1. This technique also results in an
effective decrease in quantum efficiency since photons which
impinge in the area of the grooves are not collected. The de-
gradation in quantum efficiency for both cases sited above can
be described by a factor nlsj, which is simply the ratio of
photoactive device area to total pet area.
The loss of efficiency due to optical reflection and absorp-
tion in the insulating layers covering the active region of
photodiodes can be minimized by choosing their thickness
to be (I) X and requiring their refractive index to be the
square-root value of the bulk semiconductor. The refractive
index of silicon at optical wavelengths is 3.4 while that of
Si3 N4 , a commonly deposited dielectric layer, is 1.9. There-
fore, silicon nitride is a reasonable material to use for anti-
reflection coating of silicon sensors. In practice, however,
silicon nitride is never placed directly on silicon because of
surface passivation problems. Usually a 100 A layer of thermal
oxide is grown before a silicon nitride deposition. To avoid
spurious optical reflections, the surface of the photodetector
active region should also be made planar. This is particularly
important if the sensor is to be used as a coherent detector.
A total quantum efficiency which combines all the effects
above can be expressed as
n = (n150 nREF)(nN+ + nDep + n,.). (6)
Consider, for example, the following key silicon photodiode
device parameters
Xi 18 pm--tub-substrate junction depth;
S 5 cm/s--surface recombination velocity;
a 600 cm-1 @ 8500 A absorption coefficient or
3800 cm-1 @ 6328 A absorption coefficient.
Sensors with these parameters can have a total quantum effi-
ciency of 0.48 for X = 8500 A light and 0.82 for X = 6328 A.
The lower quantum efficiency for longer wavelength light is
due to a large fraction of photons being absorbed below the
tub substrate junction.
The quantum efficiency of a MOS-CCD can be determined
using the terms in (6) but noting that there exists no junction so
nN+ is equal to zero. The quantum efficiency of CCD's is in
general lower than for photodiodes due to the lower trans-
missivity of the gate material covering the depletion collec-
tion area. Polysilicon is usually used as the transparent gate
although tin oxide is sometimes used and is of higher trans-
missivity [293.
B. Responsivity
Silicon photodiodes and CCD's have a narrow-band response
in the near infrared and visible region of the optical spectrum
given by the relationship
_ In
R~ by
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
where h is Planck's constant, v is the optical frequency, q is
the electronic charge and n is the total quantum efficiency.
The units of spectral response are amperes per watt or current
density per irradiance input. The detector responsivity at a
specific wavelength is given by the relationship
Aqn
RD= ARC = by
where A is the active area of a pet. The units of responsivity
are charge per exposure input or coulombs per microjoule per
square meter. As an example, consider the pet elements of a
detector to be 12 j. m X 20 um and n to be 0.5 at X = 8500 A
and 0.80 at 6328 A. The responsivity is
RD>, = 8.2 X 10-17 coulombs/?J/m2 - 513 e-/?J/m2
at 8500 A and
RDX = 9.8 X 10-17 coulombs/pJ/m2 - 610 e-/pJ/m2
at 6328 A.
C. Crosstalk Between Pels
An important consideration in the use of photodetectors for
AO signal processing is the relationship between the sensor
resolution and the resolution of the optical system. The key
physical property of the semiconductor which determines to
what extent photons incident upon a given pet contribute to
signal charge within that pel or adjacent pels is the absorption
coefficient. Fig. 5 shows the measured intrinsic absorption
coefficient for silicon and germanium as a function of wave-
length [30] . Also shown are the emission wavelength of some
laser sources used for AO signal processing. The photon flux
density 4)(x) as a function of the depth into the semiconductor
is given by an exponential function
'F(x) = 4)(0) a-> Vbi, the maximum charge handling capability of a
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
BORSUK: PHOTODETECTORS FOR AO SIGNAL PROCESSORS
photodiode in units of charge reduces to the expression
QT = 2(VbiVR)112Cio +CpVR. (18)
Typical realizable values for the terms in (18) for a silicon
photodiode array are
Cj0 = 0.01 pF
Vbi 0.6 V
Cp=0.1 pF
VR=10V
which yields a charge capacity of 5 to 10 X 106 e-.
The charge capacity of a shallow n-type buried channel CCD
is given by the relationship [ 10]
QT = 0.75 NDtjA (19)
where ND is the donor concentration of the channel, tj is the
n-p junction thickness, and.A is the area of the photon collect-
ing CCD gate. Using values of ND = 3 X 1016 cm, ti = 0.3 ?m,
and A = 10-6 cm2 , the number of photocarriers which the
typical shallow buried channel CCD can hold is ' 7 X 105
electrons.
E. Intrinsic Detector Transfer Characteristics
The use of photodetectors for AO signal processing place a
premium on the fidelity of the sensor's amplitude and time
transfer characteristics. The intrinsic time response of both
CCD's and photodiode arrays to pulsed optical signals is de-
pendent upon the time necessary for photogenerated carriers
to be collected. In the depletion region of both CCD's and
photodiodes, the free carriers are collected quickly. The aver-
age transit time to for an electron in the depletion region is
given by the relationship
Xd _ 2Kse0(V+ Vbi)IRNA
t0 2VS 2V5
where Xd is the depletion width, Vs is the saturation velocity
of elections, Ks is the relative dielectric constant of silicon, eo
is the permittivity of free space, Vbi is the built-in junction
voltage, V is the reverse bias on the diode, and NA is the ac-
ceptor impurity concentration. As an example, for NA = 1016
cm-3 and V = 10 V, the average electron transit time in the
depletion region is 5.6 ps. However, in the volume outside the
depletion region the free carriers must diffuse to the depletion
region. In imaging arrays which do not have a tub-substrate
junction to cutoff deeply absorbed photons, the diffusion time
of charge to the pel can be many microseconds. Photocarriers
In the case of CCD's, this capacitance is at the output sensing
node. In the case of photodiodes, every pel has both a junc-
tion and a parasitic capacitance.
The linearity of a photodiode's response is dependent upon
the ratio of the nonlinear junction capacitance to the parasitic
capacitance. The response of photodiodes that are intensely
illuminated within one exposure interval is quite complex.
The diodes can be driven into forward bias. In the ideal case,
their voltage response to light intensity then has a logarithmic
transfer function similar to solar cells. However, in realizable
photodiodes they act like leaky integrators in this mode.
The MOS CCD converts charge to voltage at the electrometer
output by discharging junction depletion capacitance. This
junction capacitance is made small by design to achieve high
signal charge sensitivity. The junction is strongly reversed
biased so that signal charge discharges it only slightly and it
thus approximates a linear response. If the incident light in-
tensity is very high, the MOS CCD register will spill charge
into adjacent pels causing blooming. This phenomenon can be
limited by adding to the sensor structure a drain such that
excess charge is removed before blooming can occur.
IV. DETECTOR DYNAMIC RANGE AND NOISE
The linear dynamic range of photodetectors for AO signal
processing can be expressed as the ratio of the maximum
charge capacity of a pel to the rms number of noise electrons,
both. expressed here in units of charge
QT
DDET __ [QT(NES)] 112
The relationship between system dynamic range and sensor
dynamic range is defined by
Dsystem KDDET (22)
where K is a function related to the system transfer function.
The inequality indicates that the system dynamic range may be
limited by the detector. The noise equivalent signal (NES) is
defined as the input exposure density, E(NES), which will
make the SNR equal. to unity at the sensor output
2 1/2
E(NES) _ RDX (p1/m2) (23)
where [QT(NES)]1I2 is the rms number of noise photoelec-
trons in coulombs and RD, is the responsivity previously de-
fined. The NES is separable into two components: temporal
noise and fixed pattern noise. This separation can be expressed
in terms of rms noise electrons by the expression
from a short optical pulse will then be collected over this time
interval. If the exposure time is comparable to the optical
pulse width, a temporal signal degradation will result. This
time delay can cause degradation in AO spectrum analyzer sys-
tems whose processing period is less than several microseconds.
This form of signal degradation is one contributor to the so
called holdover crosstalk. The use of the tub-substrate junction
to enhance pel-to-pel crosstalk performance is thus also impor-
tant to limit temporal crosstalk in both CCD's and photodiode
arrays.
The physical mechanism of electron-hole creation by the
absorption of photons is identical for both MOS-CCD's and
photodiodes. However, the linearity of the charge-to-voltage
conversion can be substantially different depending upon the
design of the sensor. The major reason for a nonlinear charge-
to-voltage ca APPa~?evPgtph?oxo }ode ~an;NpNgwt10
junction on depletion on capa u 1 ~ ~
LVTl1-4nol-1 __ - [QT
emporal + Q patial11/2 (24)
q q
where q is the electronic charge in coulombs. A summary of
these noises sources is presented in Table I for CCD's and
photodiode arrays and discussed in detail below. The choice
between selecting a photodiode, a CCD, or a photodiode/CCD
for a given AO system application is dependent principally
upon tradeoffs between speed of operation, sensitivity, and
dynamic range. In general, for low temporal noise perfor-
mance (e.g., QT(NES)/q < 200 a-) and moderate dynamic
range (102), CCD's are superior to photodiode arrays princi-
pally because of the low output sensing capacitance at the
electrometer as opposed to the photodiode capacitance. On
ilQl sQB o5153iRQ0&1((1) 3118z), sensitivity is
 0-IT, OHNSON.JYSUIST
MOS ELECTROMETER
SIGNAL PROC. AMP
FIXED PATTERN
AID QUANTIZING NOISE
principally limited by amplifier noise common to both CCD's
and photodiode arrays to about 500-1000 a-, making the
higher dynamic range (105) obtainable with photodiodes
attractive.
Several optoelectrical techniques exist to extend the dy-
namic range of the photodetector. One technique utilizes two
linear photodetector arrays and a beam splitter [33]. Light
intensity is divided unequally between the two detectors.
After the photodetector which receives the majority of light
intensity saturates, the second detector array output is uti-
lized. In this way the dynamic ranges of the two detector ar-
rays can be combined to yield a total detection range which
is the sum of the dynamic range of the individual photodetec-
tor arrays. Another elegant optical scheme for extending the
detector dynamic range of a one-dimensional processor uses a
two-dimensional detector array illuminated by means of a
cylindrical lens through a stepped density filter [34]. The
rows of the detector parallel to the direction of changing opti-
cal transmissivity become a measure of the optical intensity
while the columns of the array represent the one dimensional
spatial transform of the processor.
A. Detector Temporal Noise
There are several fundamental sources of temporal uncorre-
lated noise in photodetector arrays. They include the Johnson-
Nyquist thermal noise associated with the charging of capaci-
tors through a resistance, the shot noise due to leakage current,
the shot noise of surface and bulk traps in MOS-CCD's, and the
signal shot noise arising from quantum fluctuations of im-
pinging photons. Other mechanisms that contribute to the
overall temporal noise properties of detectors are peculiar to
the specific circuit implementation of the detector. They in-
clude electrometer amplifier noise, noise due to transfer ineffi-
ciency in MOS-CCI)'s, and analog and digital processor noise.
The noise associated with the charging of capacitors is of
importance since this mechanism is used repeatedly in photo-
detector arrays for the conversion of signals from charge to
voltage or vice versa. Referring to Fig. 8, a capacitor whose
initial charge is zero is to be charged to a voltage Vo through
a resistor R by closing a noiseless switch. The voltage across
the capacitor (assuming the capacitance is independent of
voltage) can be expressed as a function of time by the
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
BULK TRAPS
LEAKAGE
PHOTON SHOT
TRANSFER INEFFICIENCY
k.` .I DHNSON-NYQUIST IN
OUTPUT RESET OPERATION
V. kTRAI T
1
Fig. 8. Schematic representation of the charging of a capacitor through
a noiseless switch showing the Johnson-Nyquist thermal noise source.
V(t) = Vo (1 - e-'/RC) (25)
or in terms of the mean charge on the capacitor
Qm(t) = CVo(1 - e-t/RC). (26)
Barbe has shown that the variance from the mean value of Q(t)
is given by the relationship [ 10]
Q(t)2 = kTC(l - e-2t/RC)
or in terms of voltage as:
V(t)2 = C (1 - e-2t/RC)
Signal detection of photodetector arrays usually involve two
operations. The first is the fast charging of a capacitor through
the action of a low resistance switch. Ideally the second
operation is the discharging of the capacitance by signal charge
generated by light when the switch is off. In fact, thermally
generated leakage current and the emptying of surface and
bulk traps also contribute to the process and thus add shot
noise. The RC/2 time constant associated with the variance of
the voltage (or charge) is different for the two conditions. For
example, in the charging mode the time constant associated
with the variance of the voltage for a typical MOS switch
operating in the subthreshold (or triode) region charging a
total capacitance of 0.2 pF (representing depletion and para-
sitic capacitance) is (104 E2 X 0.2 pF)/2 = 1 ns. The time nec-
essary for the capacitor to reach better than 99 percent of Vo
is 5 RC = 10 ns. Consequently, the mean deviation at the end
of the charging time, t = 5 RC, is
-rkT11/z
Vn C /J
Qn = (kTC)112 (30)
in charge.
Assuming that the low resistance switch is turned. off at the
end of the charging period (Roff > 1012 2), the RC/2 time
constant of the variance is now 0.1 s. This phenomenon can
be used to obviate the noise source by the sequential sampling
and differencing of signals commonly called correlated double
sampling (CDS) [351. Physically, an ideal capacitor has no
sources of noise and hence the voltage set across it is as exact
as is the ability to measure it.
The junction and depletion capacitances associated with
semiconductors are not ideal because of the presence of leak-
age current and intraband traps. Leakage (or dark) current can
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
BORSUK: PHOTODETECTORS FOR AO SIGNAL PROCESSORS
be described by the relationship [361.
iL = gni rA Bw/ +AsS + iD (31)
`TG.R /
where n1 is the intrinsic carrier density (ni ^' exp (-EG/kT)).
AB and AS are the areas of the bulk junction and surface, w/ is
the width of the depletion region, TG.R is the bulk generation-
recombination lifetime, and S is the surface recombination
velocity. Typical values for silicon integrated circuits are
TG.R = 100-300 ?s and S:= 1-5 cm/s. The diffusion current,
iD, is due to the diffusion. of carriers to the depletion region
which are thermally generated in the neutral bulk within a
diffusion length L = DT(,. A typical value for L in silicon
integrated circuits is ' 250 ?m. Therefore, without a junction
structure for limiting this diffusion of carriers to the depletion
region, the entire bulk semiconductor can contribute leakage
current to the collecting region. The effect of the leakage cur-
rent is two fold. First, the semiconductor depletion capaci-
tance is discharged in time. Second, the generation process is
random, and this contributes a temporal noise associated with
the sensing of the charge or voltage on the capacitance. The
process obeys Poisson statistics and can therefore be described
as shot noise
Q3z (Leakage) = q iLT (32)
where 7 is the sample time interval. The sample time interval
can correspond to the exposure time of the detector as photo-
generated charge is integrated or it can correspond to the time
interval associated with the processing of photogenerated
signals.
Another source of shot noise associated with photodetectors
which shift charge signals by using buried channel CCD shift
registers is bulk trapping noise. This source of circuit noise can
be very small because the trapping occurs with deep lying
energy levels with time constants greater than 10 is. An ex-
pression for high signal level bulk trapping noise is given
by [101
Q,, (trap) = MNS Nt (1 - e z/re) a-r/Te (33)
ND
where M is the number of analog stages, NS is the number of
signal charge, and Nt is the density of bulk traps in the charge
collecting volume. The term ND is the doping density of the
buried channel, r is the time the charge packet is under a gate,
and Te is the emission constant for a bulk trap. If the bulk
trapping noise is computed for the MOS sensor element then
M = 1 and r is the exposure time. If it is being computed for
the CCD shift register, then M is equal to the number of trans-
fers and r is the transport clock period. An MOS-CCD detector
for use with a high-speed AO processor (i.e., short exposure
time and high-speed readout) suffers little noise contributed
by bulk traps. However, bulk traps can be the predominant
detector noise source in a well designed, cooled CCD for use
in long integration time AO systems such as time-integrating
correlators and radiometers.
The shot noise due to the random fluctuations in signal
photons incident upon the detector is given by the relationship
Q2 (hv) = gRDX(DT (34)
One Set Per Pixel I Common
To All Pixels
Address
Reset
Switch
Switch
OR VR
(b) CCD Read-Out Circuit
Fig. 9. Schematic representation of photodiodes and CCD arrays
showing the key difference between circuit elements.
Temporal circuit noise sources common to both CCD and
photodiode arrays arise principally at the output electrometer
where charge is converted to voltage, and in the subsequent
analog signal processing amplifiers. The first temporal noise
source of the detector output circuitry is the (kTC)112 noise
created by the discrete resetting of the combined gate and
parasitic capacitances on the MOS electrometer. The values
and sources of these capacitances are not the same for CCD's
and photodiode arrays. Fig. 9 shows the critical difference
between a CCD and. a photodiode array.
There are three contributors to the total capacitance in the
photodiode case. They are the junction capacitance C1 the
photodiode parasitic oxide capacitance CP and the gate and
parasitic capacitance of the electrometer gate CI. The sum
(Cp + C/) is directly proportional to the charge handling
capacity of the photodiode and is sometimes made large by
design to achieve large signal handling capability. In this
case, signal detection is achieved by the discharging the ca-
pacitor by photon-generated carriers. The signal charge is then
the difference between the original number of coulombs
placed on the capacitor by the reset action and that number
less those removed by recombination with photogenerated
charge. A net increase in dynamic range of Ci + Cp can be
achieved at the expense of low light level sensitivity by adding
parasitic capacitance at each diode. This method of increasing
diode capacitance is preferred to increasing junction capaci-
tance because the latter also adds leakage current. The detec-
tor sensitivity suffers however since the (kTC) 112 noise also
goes up with increased capacitance. The MOSFET electrometer
capacitance CI is directly proportional to the white noise den
sity for both CCD's and photodiodes. Therefore, if 1/f noise is
neglected, it is always of advantage to make CI as small as
where RD is the responsivity previously defined and 4) is the possible.
irradiance~or photon flux). In the case of CCD sensors, the signal charge is directly pro-
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1,.JANUARY 1981
portional to the number of photons incident upon the detec-
tor. The sensing of this signal charge is then performed by a
diode whose capacitance is independent of the charge handling
capacity of the array. The smaller is this capacitor, the greater
is the sensitivity of the device to low light (charge) levels and
the smaller is the (kTC)1/2 noise associated with the resetting
of the combined. gate and capacitance of the electrometer
(CG + Cp). Therefore, if (kTC)1/2 is the dominant noise
source, CCD arrays will have greater sensitivity than photo-
diode arrays unless CDS is used.
The second temporal noise source is due to the MOSFETS
which are typically used as both the on-chip electrometer and
the signal processing amplifiers. Their noise variance can be
expressed in terms of an equivalent input noise charge by the
relationship
Qn = En'fnCf (35)
where 1s fn is the noise bandwidth and CI is the input capaci-
tance of the transistor. The term E,2, is the input noise density
in (V2/Hz) which is composed of white noise and 1/f noise [37]
a
En = 4kT 2 1 + KoIDs
3 gm gemfb
The 1 If noise term (last term on RHS) is an empirical relation-
ship. The term ILLS is the drain current of the MOSFET in the
high current pentode regime of operation. Ko is a device-
related constant inversely proportional to gate area. The terms
a and b are experimentally determined such that the data fits
the equation and have a range typically of 0.5 to 2 for a and
about unity for b. The term gm is the transconductance of the
transistor defined here for the pentode operating regime as
W 11/2
gm = 2/1C0. L Its 1 (37)
where ji is the carrier mobility, W/L is the transistor width-to-
length ration, Co, is the gate capacitance per unit area, and
IDS is the source-to-drain current. The frequency fa at which
the 1/? noise density equals the white noise density level is
given by
3 KoIDs
fa 8 kT gm. (38)
This frequency usually ranges from 10 kHz to 100 kHz for
MOSFET's. In particular, measurements made at Westinghouse
on p-channel enhancement mode MOSFET'S used for photo-
sensor amplifiers are typically no more than 60 kHz [38].
The amplifier input capacitance CI consists of the parallel
combination of parasitic and gate capacitance. It can be
shown that for minimum noise, the parasitic capacitance
should be made equal to the gate capacitance. Substituting
(36) and (37) into (35), and using the above condition, the
noise variance in terms of charge at the input of the amplifier
can be written as
2 ~8kT(2.CooW)3/2L5/2 + KoIa (2C0xW)L31
Qn(Elect)_ DS I
The substitution of only (36) into (35) yields the result
Fig. 10. Sampled data output response showing the effect of inadequate
amplifier bandwidth.
This formulation indicates the importance of having a high
g,,,/Cr at the electrometer. The higher the gr?ICI, the more
sensitive is the detector output to photo generated charge as
well as the inherent charge noise sources (e.g., leakage, traps,
etc.) and the less sensitive to off-chip amplifier noise.
The sensor's output amplifier bandwidth B0 and its relation-
ship to noise bandwidth Afn is an important subject because
of the discrete sample nature of the analog voltage or current
signals reconstructed at the output of the sensor. Although
a bandwidth of only half the sample rate is required to prevent
aliasing, an effective loss in gain and accuracy is introduced if
the signal processing amplifiers lack sufficient bandwidth to
reconstruct the discrete samples. Referring to Fig. 10, the
relationship between amplifier gain, bandwidth, and sample
rate is given by
AV(eff)=AV[l - exp (-irB0/ff)] (41)
when B0 is the -3-dB bandwidth of the sensor amplifiers and
ff is the clock frequency. If the amplifier has a bandwidth
much greater than the sample rate, the discrete output samples
will be an excellent representation of photogenerated charge
packets collected by the sensor during an exposure time, but
the noise power is increased due to the larger bandwidth.
However, if the amplifier bandwidth is insufficient, the rela-
tionship between discrete charge packets and output voltage
samples will not be perfect (for a general description. of this ef-
fect for periodic waveforms see [39] ). Thus a sample accuracy,
SA can be defined as the ratio of the output sample value ob-
tained for a specific amplifier bandwidth to that obtained for
the case B0 >>fe,. An exact expression for sample accuracy
depends upon considering all previous samples in time since
the output response is composed of a transient (continuously
changing optical signals) and steady state response. However, a
reasonable simplification is to consider only the last previous
sample in time. This leads to the following expression for the
accuracy SA,
V 1
SA = 1 - Vl exp (-IrB0/f,)] (42)
2
where V1 is the amplitude of the Nth signal sample at the time
t1 and V2 is the next sample at time t3. An error, E, can be
defined for achieving a specified accuracy by the relationship
E= (1 - SA) = Vl exp (-irB0/f,). (43)
2
Therefore, the bandwidth as a function of clock frequency and
sample error can be written as
Qn2 al~ 1 (Elect.) = 4kT + CI 0 DS K I lAfn (40) B f` to V 1 0
g
(44)
+ B005~3RO01 U(1280g
A pro ~ le s~PA06/03/10: CIA-RD
P88
~U 1
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
As an example, consider an error of 0.01 percent is required
and that V1 / V2 is equal to 0.5. The required bandwidth is
then 2.7 times the clock rate, or 5.4 times the maximum signal
frequency.
The relationship between signal processor amplifier noise
and bandwidth can be obtained as an equivalent input noise by
evaluating the integral equation
Qre(P) = CI fo I T(if) I2 En df (45)
here T(lf) is the transfer function of the electronics. The
w
equivalent input noise density, En, has been shown earlier to
be composed of white and 1/f noise components while the
transfer function is dependent upon the amplifier gain-
frequency characteristics and characteristic response of the
sampled data. The transfer function can be written as
jrsf
T(if) TO (1 n, ) 1 (1 - lIT')
if ] 1+7(f/Bo e
(46)
The first term on the right is due to periodic sampling while
the second term is due to the relationship between the data
sample hold time rs and the clock period To = 1 /fc. The third
term is the characteristic rolloff of a one pole filter response
(typically achieved in operational amplifiers by the addition of
a feedback capacitance) where B0 is the -3 dB bandwidth.
The last term is the characteristic frequency response of a cor-
related double sampled filter circuit [341. The term r' is the
time between the clamp and sample operation of a given pel in
the sample interval. Some reasonable assumptions can be
made to evaluate the combined transfer function of (46). The
frequency characteristic of a one-pole filter can be written as
{I/[ 1 + (f/B0 )2 11. The frequency characteristic of a CDS filter
with r' _ (2)To is given by the relationship (1 - cos7rf/f f). The
frequency characteristic of a sample and hold circuit inherent
in any analog-to-digital (A/D) converter is (sinc2 rrf/fc) for
rg = To. The CDS filter has double zeros at f = 2nff,
n = 0, 1, 2, 3, ? ? ? which suppress 1/f noise, while the sample-
and-hold circuit has zeros at f = kff, k = 1, 2, 3, ? ? . , which
suppress clock noise. The relationship between these functions
is shown in Fig. 11.
The effectiveness of a CDS circuit in removing the (kTC)1/2
noise of the output electrometer is dependent upon the off
resistance of the resetting switch and the time interval [40]
Qeff(CDS) = kTCI[2(l - eT uRoffC7] (47)
This relationship shows that effectively all of the (kTC)112 can
be removed if R0 ffC1 >> T'. However, if the opposite is the
case, (kTC)112 noise will not be removed but will increase by
the factor Therefore, for AO radiometry applications
which require long integration times, the off-resistance of the
reset switch must be very high to effectively use CDS. Values
of 1012 12 for this type switch are realizable. On the other
hand, for AO spectrum analysis systems in which fast access
time is required, the integration time is usually in the micro-
second range. Therefore, CDS can be used effectively to re-
move (kTC)112 noise in these applications. However, for high
speed detection, the noise bandwidth must be large. The
amplifier noise can then be the dominant noise component
over the (kTC)112 noise associated with sensor reset operations.
In such cases, the implementation a CDS circuit is justified
T2
0
[_+(f/B0)2]
2T
02
T
0 2
sinc2( of/f(;)
T
02
Fig. 11. The first cycle components of the transfer function of sampled
data sensors.
only if the 1 If noise density is large. It can be shown that the
effectiveness of CDS in removing 1/f noise goes as fa/ff which
means that at high clock rates CDS filters 1/f noise very ef-
fectively. A recent treatment of the response of a CDS circuit
to 1 If noise is presented in [ 41 J.
An important source of temporal circuit noise which cannot
be ignored is that which arises from the A/D quantization of
signals at the sampled data output of the detector. The vari-
ance of this noise is given by the relationship [42]
2
Q(A/D) 2
where QI is the quantizing interval. For photodetectors which
have a linear transfer function between exposure density and
output signal, QI is given by the relationship
QI = AE (/1J
QB /m2 /bit) (49)
where AE is the mean change in exposure density and AB is
the mean change in output bits of the A/D.
B. Fixed Pattern Noise
Fixed pattern (or spatial) noise is due to unequal gain and/or
different dc offsets between pels of photodetectors and corre-
lated noise due to digital feedthrough from switching circuits.
The effect of fixed pattern noise is to reduce the dynamic
range of the system since any serial block signal processing
performed on the detector's output must include the variance
between pels due to nonuniform response.
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO, 1, JANUARY 1981
The characterization of fixed pattern noise is best performed
by experimental. techniques. In the dark, the standard devia-
tion about the mean of the sampled signal voltage of the array,
Vdc, can be measured. Measuring the gain nonuniformity is
somewhat more tedious. The array is uniformly illuminated
with light at several levels of irradiance. At each level, the out-
put voltage of each pel is measured and its dc offset is sub-
tracted to obtain a normalized response. The standard devia-
tion of the normalized response can then be obtained for each
irradiance level and used to set an upper boundary on input
optical irradiance. Using this largest value, a worst case noise
voltage due to nonuniform gain is determined. The two noise
voltages can be converted into an effective noise charge at the
sensor input by the relationship
[U,2(Spatial)11/2 = CT [ Vac + VG2 11/2 (50)
AV (eff)
where CT is the effective sensor input capacitance and A V (eff)
is the voltage gain from sensor element to output. For detec-
tor systems with linear transfer functions, spatial nonuniformi-
ties are dominated by dc offsets which are typically in the
range of 10-20 mV.
The most common sources of unequal offsets in photodetec-
tors are
1) nonuniform leakage current in the sensor region;
2) nonuniform leakage current in the analog shift register;
3) threshold voltage offsets;
4) nonuniform reset switch feed through.
The most common sources of unequal gain in photodetectors
1) nonequal volumes for photogenerated carriers;
2) differences in :integrating capacitance of the photocollect-
ing node;
3) nonuniform amplifier gains in the analog signal processing
chain.
The effect of fixed pattern noise can be removed at some
expense. A signal processor capable of adding or subtracting a
predetermined offset level and multiplying by a predetermined
weighting coefficient to each pel can be used to improve
photodetector fixed pattern performance. The processor must
have sufficient capacity to store all offset levels and multiplica-
tion coefficients. It must also be able to operate at the maxi-
mum data rate of the detector or suffer the loss of information.
Signal processors of this type have been built for use with bulk
AO spectrum analyzer systems but are substantially larger in
volume than the optical processing part of the system.
C. Signal-to-Noise Ratio
The SNR of the detector for noncoherent detection can be
written as
S RDx d>T
~N)DET IQT (NES) + gR
DX r1112'
The contributions to the equivalent noise charge may be sum-
marized in terms of rms noise electrons by the relationship
[QT ( gBS)]'/2 = 8 [kTC
T+Qn (Leakage) + Qn (traps)
+ Q,2~ (Elect,) + Qn (P) + Q;~ (AID)
+ Qn (spatial)].
is only present in CCD's and can be limited by careful device
processing. The noise contribution of the post electrometer
amplifier processor may be reduced by having a high gm/Cj at
the electrometer. Section IV-B described an experimental
method for measuring the spatial contribution of noise to
the NES. The combined contribution of all temporal noise
sources can also be related to a measured sensor output voltage
by the relationship
[Qn (temporal)1'12 _ CTVO (temporal)
q qAv (eff)
where Vo is the rms noise voltage over the noise bandwidth
measured for a specific pel of the sensor. This measurement
can be made by commanding a sample-and-hold device to
sample the same pel of the sensor every exposure time, low-
pass filtering of the sample and hold signal to the noise band-
width, and then measuring the result with a true rms voltmeter
whose bandwidth is at least the noise bandwidth. Since there
is no viable mechanism for achieving photoelectric gain in
present solid-state photodiode and CCD sensor arrays, single
photon shot-noise limited performance cannot be achieved.
(Present silicon avalanche photodiodes do not have uniform
gain characteristics as a function of reset voltage and tempera-
ture and are thus not presently considered for large arrays.)
However, examination of (51) shows that an improvement in
SNR by a factor of (T)112 can be achieved when QRDX4'T >> Q7'
(NES). The detector is then photon shot noise limited. This
mode of operation is possible by adjustment of either the irra-
diance 4 or the integration time r. The practical limit to this
improvement is set by the finite charge handling capacity of
the array (i.e., the sensor output saturates) and by the photon
shot noise due to optical scattering (which has been neglected
in the expressions presented here).
V. PITCH AND DETECTOR PEL REQUIREMENTS
A. Bulk Acousto-Optic Processors
The number of resolvable spots N can be expressed by:
N= Bf = (f_)kh(x) (54)
a
where ka is a factor which accounts for optical aperture weight-
ing and beam shape while h(x) is a factor which specifies reso-
lution. The term d is the width of the limiting optical aperture
in the system while Va is the velocity of sound in the Bragg
device. The term Af is the instantaneous AO bandwidth of the
Bragg device while B is the bandwidth of a resolvable element.
The term ko is equal to WtWa/l.22 for a circularly symmetri-
cal beam. Wt and Wa are the weighting functions for truncation
of the Gaussian beam and aperture apodization respectively.
The term h(x) is unity if resolution is specified by the Rayleigh
criterion. This specification is arbitrary and defines as resolv-
able two adjacent equal intensity spots separated by a dip of
-0.97 dB. If greater resolution between adjacent spots is re-
quired, h will be less than unity. For example, if a -3 dB dip
in intensity is required between adjacent spots of a uniformly
illuminated circular aperture, a factor of 1.15 (h(x) = .87)
greater resolution than is obtained using the Rayleigh criterion
is necessary. The requirement to measure the dip between
diffracted spots sets the minimum number of pels Mp at twice
the number of resolvable spots for AO processor sensors. The
pitch P (center-to-center pel spacing) of the array can be deter-
mined by equating the following geometric expression:
Techniques to imp ove detec~to SON a e the f die M P
Johnson-Nyquist 4WYJ 00 V i A$ ~fnbYse CIA-RDP88B00553R0 0128603-6 (55)
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
BORSUK: PHOTODETECTORS FOR AO SIGNAL PROCESSORS
Zero Beam
Detector
140 Element
Focal Plane
Array
Lo &
Amplifiers
Fig. 12. Schematic representation of a real-time IOC spectrum analyzer.
to the half angle of optical deflection obtained over a band- for use with an IOC system is in the detector pitch and the
width Of, given by 0/2 = (X/2u)Af where F is the focal length mounting of the detector to the IOC. The pitch of the detec-
of the transform lens in the system. Using the small angle tor, as in the bulk case, is proportional to the focal length
approximation, the expression for P is then given by of the transform lens and to the frequency resolution given
by [451
1
P=(AF/2dk) h(x).
In this analysis, the limiting optical aperture is set by the
beam expander rather than by the rectangular aperture of the
AO device. The beam expander output is circularly symmetri-
cal. Optical truncation and apodization of the beam cross
section can be performed at this point before cylindrical optics
channel the light through the AO device. The beam is then
passed through additional cylindrical and spherical optics
which reconstruct circularly symmetric spots and far-field
patterns. An advantage of this optical arrangement over the
use of an anamorphic beam expander is that apodization of
the circularly symmetric beam expander output by modifica-
tion of the Gaussian cross-section intensity profile results in
circularly symmetric spots at the photosensor. A circular spot
shape is advantageous because it minimizes the height H of the
pel elements for a given spot diameter. The minimum height
H of the pel elements is set equal to the diameter of the spot
specified as twice the radius of the first zero of the far field
diffraction pattern and can be expressed by
H= 2FX/dko. - (57)
The area of a pel is simply the product of H and P while the
aspect ratio is the ratio of H to P equal to 4h(x). The pitch
of the array can be matched to the required spot dimension
by the choice of the focal length of the final lens. However,
the optimum sensor design for AO systems need match the
aspect ratio as well. For example, if resolution is specified by
the Rayleigh criterion, h(x) equals one. The resulting aspect
ratio is four. Typical values for commercially available photo-
sensor arrays are of the order of 1-1.4. Note however, that
the more stringent the requirements on resolution, the lower
the value of the aspect ratio.
Neff Vs
where Neff is the effective index of refraction for a given opti-
cal mode within the waveguide and Vs is the acoustic velocity
of the surface acoustic wave (SAW) across the optical aperture.
However, unlike bulk optical systems, the focal length F is
limited principally by waveguide attenuation and the ineffi-
ciency of waveguide reflectors. The smaller the pitch the more
advantageous is the detector for IOC. Typical commercially
available photodetectors described in Section III have a pitch
18 pm, while the two detectors designed specifically for use
with IOC, have a pitch of 12 ?m and 8 ?m. Reducing pel
pitch below these dimensions must await new technology
developments associated with integrated circuit lithography
and processing.
The mounting of the detector to the IOC requires precision
butt coupling. The sensor must be brought to within several
micrometers of the edge of the integrated optical planar guide
to avoid light loss and optical crosstalk due to the rapid
divergence of the beam. The divergence of the beam is deter-
mined by the numerical aperture of the waveguide. To prevent
light from being reflected back into the waveguide where it
will undergo internal reflections and possibly return to the
detector array in adjacent pels, the detector array used for
the IOC spectrum, analyzer can be placed at a 450 angle to the
edge. Light reflected off the surface of the detector cannot
then re-enter the waveguide. Fig. 13 shows a schematic repre-
sentation of this arrangement, while Fig. 14 shows the actual
photodetector as an integral part of a working IOC spectrum
analyzer.
VI. ADVANCED MICROELECTRONICS FOR ACOUSTO-OPTIC
SIGNAL PROCESSING
B. Integrated-Optic Processors The technology challenge of applying advanced microelec-
Integrated optical circuits (IOC) consist of both active and tronics to AO signal processing is just beginning. Substantial
passive optical devices formed on planar dielectric waveguides. innovation in the design and fabrication of microelectronic
Optical signal processing functions can be performed by the circuits to perform high-speed processing over a wide dynamic
combination of these elements. Because of their inherent range of optical signals in both linear and area arrays are re-
planer construction they are well suited to performing one- quired at this time. Optical processing systems have demon-
dimensional real-time spectrum analysis and spatial or time strated the capability of performing 1015 complex multiplica-
integrating correlation/convolution functions [ 201, [431 , [441. tions per second.. Much of this computational power is useless
Fig. 12 shows a schematic diagram of a real-time IOC spec- without the development of advancedy}photodetector systems
trum analyzer. ThAPpGQle*iFotfi tea$etl~ 0 k?-: Qt RDf 88 (m553Rbfb1 $6% y6amic range optical
detector intended for use with a bulk system and one intended detectors but have as an integral part of their design high speed
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
Fig. 13. Schematic representation of a butt coupled sensor to an IOC.
Fig. 16. Photomicrograph of an experimental groove isolated photo-
diode array fabricated in (110) silicon.
INTEGRATED
LIGHT
INTENSITY
Fig. 14. Photograph showing a sensor butt coupled to a working IOC
spectrum analyzer.
NEW VHSIC
PHASE I
CANDIDATE
Fig. 15. Block diagram of advanced AO processors utilizing VHSIC
technology.
data processing. Two potential technology candidates for
these advanced "smart" sensors are silicon very high speed
integrated circuits (VHSIC) and GaAs microelectronic inte-
grated circuits.
A. Very High Speed Integrated Circuits Technology
The VHSIC program is a De t f
t
D
TIME
Fig. 17. Optical response of one pel of the array showing the linear to
logarithmic response during an integration interval.
Cco
CIRCULATING
SIR
Fig. 18. Schematic representation of a high-speed detector preprocessor
for use with AO instantaneous spectrum analyzers.
par men
o
efense initiative to circuits. Specifically, the capability to accurately define lines
develop high speed, low power digital silicon integrated circuits and spaces of I Am will allow fabrication of small pitch, high
composed of greater than 104 digital logic gates per die. This density photosensor arrays. This capability can also support
transistor density requires the technological development of the fabrication of high-speed peristaltic analog CCD's and
new lithographic and processing techniques. A figure of merit amplifiers to perform the high-speed scanning and subsequent
for VHSIC is 5 X 1011 (Gate-Hz)/cm2. The effective utiliza- processing of large arrays. Finally, the ability to accurately
tion and implementation of VHSIC circuits requires accurate control the concentration, depth, and lateral diffusion of im-
device and circuit models and computer aided design (CAD). purity dopants in the circuit will invite the development of
VHSIC technology will impact optical signal processing in novel detector and signal processing structures. The second
two ways as illustrated in Fig. 15. The first benefit is the us , e ' nt of a high-speed
of VHSIC techniquApprGyQd p9CtfiiR4ll ?gf2 (eIs'6( : li iE o s o w ft i pperro~r~~m~s~~ many specialized functions
Si"111 CCU
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
unique to the AO signal processing requirement. Such an
AO VHSIC processor could be used to perform detection
algorithms and pattern (signature) recognition algorithms of
optically processed data as well as remove fixed pattern noise
and other dilatorious properties of the sensor.
An advanced high speed detector preprocessor silicon micro-
electronic circuit for use with integrated and bulk real time
AO spectrum analyzers is being developed [461. The circuits
will be capable of analyzing and sorting by amplitude fre-
Approved For Release 2006/03/10
quency and time of arrival wide dynamic range optical signals
presented at its focal plane. The photosensor part of the cir-
cuit utilizes narrow grooves to isolate diodes dielectrically.
This unique geometry serves to limit adjacent pel electrical
crosstalk and to allow operation of the diode sensors into their
nonlinear forward biased region so that high optical dynamic
range can be achieved. A photomicrograph of an experimental
groove isolated photodiode array fabricated on (110) silicon
is shown in Fig. 16 while the linear to log response of a partic-
: CIA-RDP88B00553R000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
ular element of this array to HeNe 6328 A laser light illumina-
tion is shown in Fig. 17 [471. The signal processor utilizes
high-speed peristaltic CCD's with on-chip bipolar transistor
amplifiers for analog implementations of comparators and
analog buffers, low noise short channel MOS charge to voltage
amplifiers and CCD internal nondestructive readout charge
to current amplifiers [48], [49]. A line to line differencing
scheme is utilized to cancel fixed pattern noise in the sensor
as well as remove CW signals. A schematic diagram of this
circuit is shown in Fig. 18.
Another example of an advanced custom silicon microelec-
tronic circuit which enables a unique AO phenomenon to be
utilized is the AO exciser processor [50). The processor uti-
lizes the coherent property of laser light to allow reconstruc-
tion of time domain signals that have been transformed by the
Bragg AO interaction. A schematic representation is shown in
Fig. 19. Laser light incident on the Bragg device is diffracted
in spatially resolved spots as in the ubiquitous AO spectrum
analyzer. However, since the optical frequency of each de-
flected spot is doppler shifted in frequency proportional to
the input electrical signal the superposition of reference laser
light upon these spots results in optical mixing [ 51 ] , [52).
A photodetector array of sufficient instantaneous bandwidth
placed to intercept these spots will have a time domain output
signal of identical frequency spectrum content but possibly
amplitude shaded by a weighting function.
If an array of photodiodes is placed at the image plane and
a means for either adding or deleting the response from any
one diode to a common output bus is implemented, excision
of particular frequencies in the time domain can readily be
achieved. Fig. 20 is a representation of such a switching photo-
array along with a schematic. Fig. 21 shows the actual silicon
device built to perform this function [531, [541. The posi-
tions of the switches, determining to which output bus each
photodiode is connected, are set by loading serially a binary bit
stream into the device's shift registers. The switch positions
are then held in this state until new data is entered. Using this
type of AO system, made possible by a unique photodetector,
narrow-band signals are separated from broad-band signals
over a wide instantaneous bandwidth. This function is use-
ful as a prewhitening filter in spread spectrum communica-
tion systems; as well as for ECM signal processors to increase
sensitivity in a jamming and a dense emitter radio frequency
environment.
B. GaAs Technology
The use of GaAs integrated circuits for advanced high-speed
sensor/processors is possible. A line array of GaAs detectors
for use with an integrated optical spectrum analyzer has been
fabricated and tested [231. The detectors were fabricated as
electroabsorption avalanche photodiodes. Measurements on
individual diodes yielded noncoherent dynamic ranges of
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
BORSUK: PHOTODETECTORS FOR AO SIGNAL PROCESSORS
I
57 dB. These detectors consisted of a line array of individually
readout photodiodes. The ability to marry GaAs detectors
with high-speed GaAs CCD's and digital logic functions may
eventually be feasible as this technology matures [55]. A
schematic representation of a 16-pel photodiode/CCD linear
GaAs array is shown in Fig. 22 [56]. The device utilizes a
Schottky diode as the sensor connected to a Schottky gate
mesa-type GaAs CCD. Mesa-type GaAs CCD's have been
operated at high speeds (500 MHz) in a continuous clocking
mode [57].
VII. CONCLUSION
A description of the properties and use of solid state photo-
detectors for AO signal processing has been presented. Special
attention has been given to the device physics which describe
the noise properties of these sensors because of their impor-
tance in determining the ultimate performance of the entire
system. The different categories of device architectures have
performance of AO systems. New silicon VHSIC and GaAs
sensor/processors will be developed in the future that allow
the full potential. of these unique optical signal processors to
be achieved.
ACKNOWLEDGMENT
The author acknowledges the many useful discussions and
contributions made by the staff of the Westinghouse Advanced
Technology Laboratory. Specifically, the critical review and
discussion of this paper by M. H. White, D. R. Lampe, N.
Bluzer, and E. C. Malarkey is greatly appreciated. The support
of the Navy in developing this technology has been significant.
Specifically, the technical contribution of T. Giallorenzi, G.
Anderson and A. Spezio from U.S. NRL and program support
from J. Koenig and N. Butler of Navelex in pursuing photo-
sensor microelectronics for acousto optics is acknowledged.
also been explained because of their differences to those 111 W. K. Middleton, Vision Through The Atmosphere.
h there have Canada: University of Toronto Press, 1952.
Althou
s
i
i
d i
g
ng sensor
.
n scene
mag
typically use
121 been some new sensors specifically made for AO applications, D. Hecht, "Spectrum analysis using acousto-optic
Opt. Eng. vol. 16, , pp pp. . 461-466, Oct. 1977.
the photodetector i( pi b?pottRgibAi e2 6y6y1V : _~~ . B %6llf66t1 6t6 f -10 nA/cm2.
liE
~fl
niliinliinitfflniinRi[ni[flniinMnTPMMNM 19
31
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
[5]
[61
[71
optic radiometry," Proc. SPIE Optics in Radar Systems, vol. 128,
pp. 344-352, Sept. 1977.
R. A. Spague, "A review of acousto-optic correlations," Opt.
Eng., vol. 16, pp. 467-474, Oct. 1977.
T. M. Turpin, "Spectrum analysis using optical processing," this
issue, pp. 79-92.
"Time integrating optical processing," in Proc. SPIE Symp.
Real Time Signal Processing (San Diego, CA), vol_ 154, pp. 196-
203, Aug. 1978.
[8] - "Real time,input transducer for coherent optical process-
ing," Int. Optical Computing Conf. (Apr. 9-11, 1974, Zurich,
Switzerland).
J. B. Barton, J..1. Cuny, D. R. Collins, "Performance analysis of
EBS-CCD imaging tubes/status of ICCD development," in Proc.
1975 Int. Conf Application of Charge Coupled Devices, pp.
133-146, Oct. 1975.
single crystal germanium and silicon at 77? K and 300? K," Phys.
Rev., vol. 99, p. 1151, 1955.
[311 A. Turley, unpublished, 1978.
D21 R. H. Dyck, W. Steffe, "Effects of optical crosstalk in CCD image
sensors," in Proc. Int. Conf. Application of CCD's, pp. 1-55,
Oct. 1978.
1331 U.S. Patent 3 962 5.77.
[341 U.S. Patent 3 934 153.
[351 M. H. White, D. R. Lampe, F. Blaha, and I. Mack, "Characteriza-
tion of surface channel CCD image arrays at low light levels,"
IEEESC, vol. 9, p. 3, Feb. 1974.
[361 M. H. White, "Photodiode sensor array," in Solid State Imaging,
NATO Advanced Study Institute Series, Series E, NO1 6.
[371 P. R. Gray and R. G. Meyer, Analysis and design of integrated
circuits. New York: Wiley, 1977.
[38] A. Turley, private communication, 1979.
[391 C. D. McGillen and G. R. Cooper, Continuous and Discrete
Signal and System Analysis. New York: Holt, Rinehart, and
Winston, 1974, pp. 238-241.
[401 R. W. Brodersen and S. P. Emmons, "The measurement of noise
in buried channel charge-coupled devices," in 1975 Int. Conf.
Application of CCD's (Oct. 29, 1975), p. 331.
[411 R. J. Kansy, "Response of a correlated double sampling circuit
to 1/f noise," IEEE J. Solid-State Circuits, vol. SC-15, p. 373,
June 1980.
[421 B. Gold and C. M. Rader, Digital Processing of Signals. New
York: McGraw-Hill, 1969, p. 104.
[43] S. H. Chang and J. T. Boyd, "An integrated optical channel wave-
guide-CCD transversal filter," IEEE Electron Device Lett., vol.
EDL-1, Mar. 1980.
[441 C. S. Tsai, C. C. Lee and B. Kim, "A review of recent progress on
guided wave acousto-optics with applications to wideband com-
munication and signal processing," Tech. Dig. 1980 Topical Meet.
Integrated and Guided-Wave Optics (Incline Village, NV), vol.
ME1-1, 1980.
[45] M. C. Hamilton, D. A. Wille, and W. J. Miceli, "An integrated
optical RF spectrum analyzer," Opt. Eng., vol. 16, no. 5, p. 475,
Sept.-Oct. 1977.
[46] Navy Contract N00173-79-C-0485.
1471 J. Kim, unpublished, Apr. 1980.
[481 G. M. Borsuk, M. H. White, N. Bluzer and D. R. Lampe, "Design
and analysis of new high speed peristaltic CCD's," in Proc. 1978
Conf. Application of CCD's (Oct. 25-27), 1978.
[491 R. J. Brewer, "The low light level potential of a CCD imaging
array," IEEE Trans. Electron Devices, vol. ED-27, p. 401, Feb.
1980.
150 1 Navy Contract N00039-78-R-0417.
[51] W. J. Thaler, "Frequency modulation of a HeNe laser beam by
ultrasonic waves in quartz," Appl. Phys. Lett., vol. 5, no. 2, p. 29,
July 16, 1964.
[521 G. M. Borsuk and W. J. Thaler, "Frequency-moduated laser com-
munication system," IEEE Trans. Sonics Ultrason., vol. SU-17,
Oct. 1970.
[531 F. Kub, unpublished, Mar. 1980.
[54] N. Bluyer, F. Kub, and G. M. Borsuk, "Steering of high frequency
low level photocurrents between summing busses beyond the
gm/21rc limit," Proc. IEDM, Dec. 1980.
[55] R. C. Eden and I. Deyhimy, "Application of GaAs integrated cir-
cuits and charge-coupled devices (CCD's) for high speed signal
processing," in Proc. SPIE Acousto Optic Bulk Wave Devices,
vol. 214, p. 39, 1979.
[561 F. Kub, Unpublished, May 1980.
[571 I. Deyhimy at al., "An ultra high speed GaAs CCD," in Tech.
Dig. 19 79 IEDM, p. 619, 19 7 9.
[10] D. F. Barbe, "Imaging devices using the charge-coupled concept,"
Proc. IEEE, vol. 63, pp. 38-67, Jan. 1975.
[11] "Advances in CCD and imaging, Session II," in Proc. ISSCC,
1980.
[12 ] M. M. Blouke, J. F. Breitzmann and J. E. Hall, "Three Phase,
backside illuminated 500 X 500 CCD imager," in Proc. 1978
ISSCC, p. 36, 1978.
[131 J. Hynecek, "Virtual Phase CCD Technology," in Proc. 1979
IEDM, p. 611.
[141 Reticon data sheet "512 random access line scanner-CP 1006,"
Feb. 2, 1980.
[15] G. Weckler, private communication, Apr. 1980.
[16] R. Ekstein, "ELINT sensor design goals," in Proc. SPIE Conf.
Workshop Acousto-Optic Bulk Wave Devices (Monterey, CA),
Nov. 27-29, 1979.
[171 Reticon Data Sheet "CP 1023 tapped 256 CCPD," Mar. 18, 1980.
[181 G. M. Borsuk, A. Turley, G. E. Marx, E. C. Malarkey, "Photo-
sensor array for integrated optical spectrum analyzer systems,"
Proc. SPIE, vol. 176, p. 109, 1979.
[19] G. E. Marx and G. M. Borsuk, "Evaluation of a photosensor for
integrated optics spectrum analyzers," in Tech. Dig. 1980 Topical
Meet. Integrated and Guided Wave Optics (Incline Village, NV),
vol. ME 6, 1980.
[201 D. Mergerian et al., "Working integrated optical RF spectrum
analyzer," Appl. Opt., vol. 19, no. 18, Sept. 15, 1980.
[21 ] D. Mergerian, unpublished, Mar. 1980.
[221 J. Y. M. Lee, B. Chen, "Detector array for an integrated optic
spectrum analyzer," Extended abstracts, vol. 80-1, Spring Meet.
Eleetro-Chemical Soc., May 11-16, 1980.
[23] "GaAs waveguide detector array for A IOC spectrum analyzer,"
Semiconductor Res. Lab, Washington University, St. Louis, MO,
Final Rep. N064351-1 July 31, 1979.
[241 Solid State Imaging, Jespers, Van De Wiele, and White, Eds.,
NATO Advanced Study Institute Series, Series E, Applied Sci-
ences, no. 16 (Leyden, The Netherlands: Noordhoff, 1976).
[251 A. S. Grove, Physics and Technology of Semiconductor Devices.
New York: Wiley, 1967, p. 121.
[26] P. A. Gary, "Modeling and optimization of a silicon photosensor
for a reading aid," Stanford Electronics Lab., Tech. Rep. No.
4822-1, May 1967.
[27] F. Van DeWiele, "Photodiode quantum efficiency," in Solid
State Imaging, NATO Advanced Study Institute Series, Series E,
no. 16 (Leyden, The Netherlands: Noordhoff, 1976).
[281 N. Bluzer, private communication, 1979.
[291 D. McCann et al., "Buried-channel CCD imaging arrays with tin-
oxide transparent gates," in Proc. 1978 ISSCC, p. 30, 1978.
[301 W. C. Dash and R. Newman, "Intrinsic optical absorption in
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88BOO553ROO0100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
Contributors
Gerald M. Borsuk (SM'78) was born in Newark,
NJ, on December 15, 1944. He received the
B.S., M.S., and Ph.D. degrees in physics from
Georgetown University, Washington, DC, in
1966, 1970, and 1972, respectively.
Since 1977, he has been employed by the
Westinghouse Electric Corporation, Advanced
Technology Laboratories, Baltimore, MD, and is
Manager of the Solid State Device Physics De-
partment primarily involved in research and
development concerning charge-coupled devices
and optical techniques for signal processing. Prior to joining West-
inghouse, he worked for the ITT Electro Physics Laboratory, Colum-
bia, MD.
Berthold K. P. Horn received the B.Sc. (elec.
eng.) degree from the University of the Witwa-
tersrand, Johannesburg, South Africa, and the
S.M. and Ph.D. degrees from the Massachusetts
Institute of Technology (M.I.T.), Cambridge,
MA, in 1965, 1968, and 1970 respectively.
He has taught at both the University of the
Witwatersrand and M.I.T. He has been Associ-
ate Professor in the Department of Electrical
Engineering and Computer Science of M.I.T.
since June 1976. He directs the machine vision
and machine manipulation work of the M.I.T. Artificial Intelligence
Laboratory. His present interests also include locomotion, tomography,
image understanding, exact color reproduction, and the protection of
wildlife.
Peter Kellman was born in New York City, NY,
on November 29, 1952. He received the B.S.
degree in electrical engineering from Carnegie-
Mellon University, Pittsburgh, PA, in 1975, and
the M.S..and Ph.D. degrees in electrical engi-
neering from Stanford University, Stanford,
CA, in 1977 and 1979, respectively.
He joined ESL, Incorporated, a subsidiary of
TRW, Sunnyvale, CA, in 1975 and has worked
on the development of automatic signal search
and collection systems using acousto-optic sig-
nal processing technology. He is currently manager of Acousto-Optic
System Development.
Dr. Kellman is a member of the IEEE Acoustics Speech and Signal
Processing and Communications Societies, Society of Photo-Optical
Instrumentation Engineers, and the Optical Society of America.
Adrianus Korpel (SM'67-F'75) received the
B.S., M.S., and Ph.D. degrees in 1953, 1955,
and 1969, respectively, from the University of
Technology, Delft, The Netherlands, all in
electrical engineering.
From 1955 to 1960 he was with the Postmas-
ter's General Research Laboratories, Melbourne,
Australia. In 1960 he joined Zenith Radio
Corporation, Chicago, IL, where in 1975, he
became Director of research in engineering
physics. Since 1977 he has been with the Uni-
Dr. Korpel is a member of both the Optical Society of America and
the Acoustical Society of America.
J. William Murray (S'59-M'63) received the
B.S.E.E. degree in 1960 from the University of
Washington, Seattle, and the M.E.E. degree in
1962 from New York University, New York.
As a Member of Technical Staff at Bell Tele-
phone Laboratories, Whippany, NJ, from 1960
to 1966, he participated in development of
signal and data processing equipment, and
associated computer software, for operation,
testing, and simulation of antiballistic missile
radar systems, and for measurement of antenna
patterns of large multifunction phased-array radars. From 1966 to
1977 he was employed by the Applied Technology Division of Itek
Corporation, Sunnyvale, CA. During this period he was primarily con-
cerned with development of system designs for radar warning systems
for use on military tactical aircraft. This included development of
computer-based interferometer systems for providing highly accurate
azimuth and elevation direction finding information. Since 1977 he
has been a Member of Technical Staff at ESL, Incorporated, a sub-
sidiary of TRW, Sunnyvale, CA. His current responsibilities are primar-
ily in connection with system design of specialized computer-based
signal processing systems. This includes development of systems for
performing real-time analysis, display, and collection of the multichan-
nel output signals available from acousto-optic channelized receivers.
Mr. Murray is a member of Tau Beta Pi.
William T. Rhodes (S'70-M'71) was born in
Palo Alto, CA, on April 14, 1943. He received
the B.S. degree in physics in 1966 and the M.S.
and Ph.D. degrees in electrical engineering in
1968 and 1971, respectively, from Stanford
University, Stanford, CA.
In 1971 he joined the faculty of the Georgia
Institute of Technology where he is currently
an Associate Professor teaching in the areas of
modern optics, communications, and signal
processing. He is coauthor of a book on lasers
and their applications and is editing a book on optical information pro-
cessing. In 1976 he was a Humboldt Research Fellow at the University
of Erlangen in West Germany. His research is in the areas of image
formation and processing and in hybrid optical/electronic signal pro-
cessing. He is an associate editor for Optics Letters and serves on the
editorial boards of Optica Acta and Optical Engineering.
Dr. Rhodes is currently the Chairman of the IEEE Computer Society
Technical Committee on Optical Processing and the Vice Chairman of
the Optical Society of America Technical Group on Information Pro-
cessing and Holography.
Harry N. Shaver (S'53-M'57) was born in Bis-
bee, AZ, on June 25, 1935. He received the
B.S.E.E. and M.S. degrees in electrical engineer-
ing from the University of Arizona, Tucson, in
1957 and 1959, respectively, and the degree of
Engineer from Stanford University in 1965.
He joined ESL, Incorporated, a subsidiary of
TRW, Sunnyvale, CA, in 1972 and has been
active in the areas of signal search, automatic
signal recognition, and optical signal processing.
He is currently Manager of the Signal System
versity of Iowa, Iowa City, as a faculty member of the Department of
Electrical and Computer Engineering. He has worked in the areas of
coherent optics, acousto-optics, and acoustic imaging. His present
research interests include nonlinear phenomena.
Approved For Release 2006/03/10 : CIA-RDP88BOO553ROO0100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
Mr. Shaver is a member of Tau Beta Pi, Pi Mu Epsilon, Sigma Pi
Sigma, and Phi Kappa Phi.
Allan W. Snyder received the B.Sc. degree from
the Pennsylvania State University, University
Park, the S. M. degree from the Massachusetts
Institute of Technology, Cambridge, MA, the
M.S. degree from Harvard University, Cam-
bridge, MA, the Ph.D. degree from University
College London, England, and the D. Sc. degree
from the University of London, London, Eng-
land, in the fields of biology, electrical engineer-
ing, and applied physics.
At present he holds the chair of optical phys-
ics and vision research in the Institute of Advanced Studies at the Aus-
tralian National University. He is Head of the Department of Applied
Mathematics, a department devoted to interdisciplinary fields including
membrane biophysics, optics and vision. He is also a professor in the
Department of Neurobiology. During the academic years of 1974
and 1977 he was a Visiting Professor at the Technische Hochschule,
Darmstadt, Germany, and at the Yale University School of Medicine,
New Haven, CT, respectively. His industrial appointments include nine
months with Peter Kiewitt as a Communication Engineer on the ice cap
in Greenland (1960), a Research Scientist with the Sylvania Applied
Research Laboratory (1963-1967), Consultant to the British Post Office
(1967-1969), the Standard Telecommunications Laboratory, England
(1969-1970), Telecom Australia (1972) and the French Centre National
d'Etudes des Telecommunications (1976) all concerned with optical
fibers. He conceived and jointly organized the first international work-
shop on photoreceptor optics in Darmstadt, Germany, 1974, and the
first international workshop on optical waveguide theory in Lannion,
France, 1976. He was one of the four invited members on the 1977
Marconi International Fellowship Committee. In addition to more than
150 papers, he is the co-author of Photoreceptor Optics (Springer Verlag,
1975) and Optical Waveguide Theory (Chapman and Hall, 1981) and
was the Guest Editor of the 1977 special issue of Optical and Quantum
Electronics.
Professor Snyder has been the recipient of a National Science Founda-
tion fellowship in biophysics, the John Simon Guggenheim Fellowship
in visual neurobiology, the research medal of the Royal Society of
Victoria and the Edgeworth David Medal.
Terry M. Turpin received his B.S.E.E. degree
from the University of Akron, Akron, OH, in
1966, and a M.S. degree in electronic engineer-
ing from Catholic University, Washington, DC,
in 1970.
From 1967 to the present, he has been
actively involved in optical processing research
and development with the Department of De-
fense. He has had extensive experience with
one- and two-dimensional optical processing
systems and hands-on experience with a wide
Shi-Kay Yao (S'68-M'74-SM'80) was born in
China on June 30, 1945. He received the B.S.
degree in electrical engineering from National
Taiwan University, Taipei, Republic of China,
in 1967, and the M.S. and Ph.D. degrees in elec-
trical engineering from Carnegie-Mellon Univer-
sity, Pittsburgh, PA, in 1969 and 1974,
respectively.
He is currently Senior Scientist at the TRW
Technology Research Center, Torrance, CA,
contributing to the development of optoelec-
tronic devices and fiber optical communication. Prior to this, he was
manager of Optical Signal Processing at the Electronics Research Center
of Rockwell International, Anaheim, CA, from 1976 to 1980, where his
major contributions were in Integrated Optics, waveguide optical lenses,
waveguide photodetectors, and surface acousto-optic devices using ZnO
piezoelectric thin film. He was with Harris Corporation from 1974 to
1976 working on acousto-optic devices and surface acoustic wave de-
vices. He has been working on optical image processing (1968-1970),
laser communication (1970-1973) as part of his Ph.D. dissertation re-
quirement, and integrated optical switching devices on LiNbO3 (1973-
1974) as post doctorate research, at Carnegie-Mellon University. He is
the author or co-author of over 40 journal articles, technical presenta-
tions and holds several U.S. patents.
Dr. Yao is currently the Treasurer of the Los Angeles chapter of IEEE
Quantum Electronics and Applications Society. He is a member of the
Optical Society of America, Eta Kappa Nu, and Sigma Xi.
Eddie H. Young, Jr. received the B.S. degree
(with honors) in electrical engineering from the
University of Illinois, Urbana, in 1961 and the
M.S. and Ph.D. degrees in electrical engineering
from the University of Texas, Austin, in 1963
and 1968, respectively.
He is currently Director of the Optical Sys-
tems at Harris GSG leading engineering activi-
ties in optical systems and optical devices
development. Prior to this, he was involved in
the study of acoustic and optical interaction in
-bulk and surface waves for applications in signal processing. He first
became interested in surface wave acoustics while on a Summer Na-
tional Science Foundation postdoctorate residence at Stanford Micro-
wave Laboratory. Subsequently, he was funded a NSF Grant to con-
tinue research at Old Dominion University. He also received a grant
from NASA to work on a computer-aided network and analysis pro-
gram. He spent 2 years as an Assistant Professor of Electrical Engineer-
ing at Old Dominion University, Norfolk, VA. At the University of
Texas, his doctoral dissertation and graduate studies primarily involved
plasmas and electromagnetic theory. In particular, he is involved in the
fabrication of acoustic cells, and in the analysis and experimentation of
the cells for correlation and spectrum analysis. In addition, some basic
research is carried out in surface wave devices. He has written numerous
technical papers.
Dr. Young is a member of the Eta Kappa Nu, Tau Beta Pi, and Sigma
Xi.
variety of electro-optical components.
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
Proceedings Letters
This section is intended primarily for rapid dissemination of brief
reports on new research results in technical areas of interest to IEEE
members. Contributions are reviewed immediately, and acceptance is
determined by timeliness and importance of the subject, and brevity
and clarity of the presentation. Research letters must contain a clear
concise statement of the problem studied, identify new results, and
make evident their utility, importance, or relevance to electrical engi-
neering. Key references to related literature must be given.
Contributions should be submitted in triplicate to the Editor,
PROCEEDINGS OF THE IEEE, 345 East 47 Street, New York, N. Y. 10017.
The length should be limited to five double-spaced typewritten pages,
counting each illustration as a half page. An abstract of 50 words or
less and the original figures should be included. Instructions covering
abbreviations, the form for references, general style, and preparation
of figures are found in "Information for IEEE Authors," available on
request from the IEEE Editorial Department. Authors are invited to
suggest the categories in the table of contents under which their letters
best fit.
After a letter has been accepted, the author's company or institution
will be requested to pay a voluntary charge of $70 per printed page,
calculated to the nearest whole page and with a $70 minimum, to cover
part of the cost of publication.
Quasi- Optics of the Waves Guided by a Slab of
Uniaxially Anisotropic Dielectric
Abstract-A ray optical method is given for deducing the dispersion
relation and the group velocity of the TM wave guided by a slab of
uniaxially anisotropic dielectric with particular emphasis on confirming
the correctness of the reverse lateral displacement of the ray on total
reflection at an interface separating the dielectric from the free space.
A planar slab of uniaxially anisotropic magnetoionic medium occupy-
ing the region -o < x, y < o, and -d < z < d is situated in free space as
shown in Fig. 1. The magnetization is in the z direction and is normal
to the surfaces of the slab. All the field quantities are assumed to be
independent of y and have the harmonic time dependence of the form
exp (-iwt). With suitable normalizations, inside the slab, the TM waves
(Hy, Ex, Ez) are governed by the following relations:
1 a2 a2
C 2 + az 2 +w2 Hy(x,z)=0
8x
Ex iw az ' Ez =
1 ally 1 8Hy
and e = 1 - 1/w2. The fields in the free space outside the slab are also
governed by (1) and (2) with e = 1. We restrict our attention to fre-
quencies below the plasma frequency, that is w < 1 and e < 0.
Assuming a plane wave solution of the form Hy (x, z) - exp [i (Oxx +
0,z)], from (1) we obtain the dispersion relation
2O2 -- w2) - 0z + w2 = 0. (3)
Also, we assume that Ox is positive real. For w < 1, Oz is also real. For
the sake of convenience, we take oz to be also positive. The group
velocity vg = z vgx + z vgz in the unbounded uniaxially anisotropic di-
electric can be deduced from (3) as
-wOx
Vgx (1 w2) [1 +02x/(1 - w2)2] w and az is positive
real leading to a total reflection in the slab. Applying the boundary
conditions that Hy and. Ex are continuous at the interface z = d yields
A2/A1 = exp (-i24) (8)
tan 4=-= (9)
kz [w2 + w20/(1 - w2)] 112.
Since v$z > 0, A2/A1 given in (8) is the reflection coefficient both for
the phase and the ray.
By considering a space-like wave packet and a time-like wave packet
separately, we can establish [1] , [2] that the ray on total reflection at
a plane interface between the uniaxially anisotropic dielectric and the
free space undergoes a lateral displacement in the x direction as given
by
2x5 = 2a$/39x = 2w2/0x0zaz > 0 (10)
and an associated time delay as given by
2w az
2ts = -2ao/8w = 1 + 2 > 0. (11)
0zaz (1 w )
In (8)-(11), 0 is a function of the independent variables Ox and w. The
incident ray is directed in the -x direction and the lateral displacement
determined in (10) is in the +x direction. Thus the ray on total reflec-
tion has a reverse lateral displacement. Felsen [3] has previously ob-
tained this result by a systematic evaluation of the integral representation
of the field excited by a simple source in a half space of the uniaxially
anisotropic dielectric. The lateral displacement can be either in the
forward or in the reverse direction but the associated time delay should
always be positive. The positive time delay obtained in (11) confirms
the correctness of the reverse lateral displacement deduced in (10).
The dispersion relation of the TM wave guided along the slab can be
deduced using the zig-zag model for the propagation of the phase fronts
of a homogeneous plane wave inside the slab [4]. The zig-zagging
wave-normal directions of the successively and totally reflected homo-
geneous plane wave inside the slab are shown by the line PIP2P3P4 in
Fig. 1. The dashed line PIP4 represents a phase front. The closure
condition is that the resultant phase change between any two points
on a phase front along the zig-zagging wave-normal directions including
the phase changes occurring on total internal reflection at the two inter-
faces is equal to an integral multiple of 27r. Symmetry considerations
show that the phase change introduced by the total internal reflec-
tion is the same for the two interfaces. Since the distance P1P2P3P4 =
4d cos 6 the total Base chap a from P to P alon the zi z,, in
(OX - w2)1/2
Manuscript received August 1, 1980. P' p g r 4 g g gg g
The author is with the Department of Electrical and Computer En- wave-normal directions is obtained as 40d cos op - 40 where 0 is the
gineering, University of Wisconsin, Madison, WI 53706. magnitude of the wavenumber in the wave-normal directions. Since
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
0018-9219/81/0100-0121$00.75 ? 1981 IEEE
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
pp. 174-185, Feb. 1974.
[3] L. B. Felsen, "Lateral waves on an anisotropic plasma interface,"
IRE Trans. Antennas Propagat., vol. AP-10, pp. 347-349, May
1962.
[41 S. R. Seshadri, Fundamentals of Transmission Lines and Electro-
magnetic Fields. Reading, MA: Addison-Wesley, 1971, pp. 428-
436.
Comments on "Resolution of Coherent Sources Incident
on a Circular Antenna Array"
Fig. 1. Geometry of a planar slab of the uniaxially anisotropic dielec-
tric showing the zig-zagging wave-normal directions and the associated
zig-zagging ray directions.
a cos Op = (iZ, the closure condition is given by
4d[w2 + w2i3 /(1 - w2)] 1/2 - 4(p = 27rn
where n is an integer. This closure condition is the dispersion relation
of the TM wave guided by the slab. The allowable values of ax and the
associated discrete values of the angle of incidence op are given by the
dispersion relation (12), together with (9) and (5a).
For the guided wave, in view of (12), w becomes an implicit function
of px. The group velocity Vg of the guided wave is in the x direction
and the x component Vgx of Vg is given by dw/d(3x as obtained from
(12). Taking the derivative of (12) with respect to ox and using (4a)
and (5b), we determine the group velocity of the guided wave to be
given by
d tan 0g - a(p/a/3x
-Vgx =
d tan og/(-vgx) - ao/aw
Associated with the zig-zagging wave-normal directions are the zig-
zagging ray directions. The direction of the ray is different from the
corresponding wave-normal direction, as already indicated. On total
reflection at an interface between the uniaxially anisotropic dielectric
and the free space, a ray undergoes a reverse lateral shift. In Fig. 1
P1P2P3P4 is a part of the zig-zagging wave-normal directions incident
on the interface at an angle Bp. RIR2R3R4 represents the correspond-
ing ray directions incident on the interfaces at an angle 0g. R 1R4 is on
the midplane z = 0 of the slab. Since the reflections at the two inter-
faces are identical, RIR2R3R4 forms a unit cell in the periodic pattern
of the zig-zagging ray directions. The times taken for the ray to advance
from R I to R 2, from R2 to R 3, and from R 3 to R4 are given by d tan og/
(-vgx), -2at/aw, and d tan 611(-v gx), respectively. Hence, the time
taken for the ray associated with the guided wave to advance from R 1
to R4 along the actually ray direction is given by 2 [d tan og/(-vgx) -
ao/aw]. In this time the ray associated with the guided wave advances
in the -x direction from R1 to R4. The distance from R1 to R4 is
given by 2 [d tan o8. - ar/a px] . The minus sign in front of a0/a(3x
specifically takes account of the fact that the lateral displacement is in
the reverse direction. Hence the ray method [2] gives the group veloc-
lity of the guided wave as
d tan 0g - ao/agx
-Vgx d tan og/(-ugx) - ao/aw
which is identical to that obtained in (13) from the dispersion relation.
The correctness of the group velocity deduced via the ray method con-
firms the reverse direction of the lateral displacement of the ray result-
ing from the total reflection at the interface between the anisotropic
dielectric and the free space.
The reverse lateral displacement was first obtained by Felsen [3].
Since the geometry is uniform in the longitudinal direction, the reversal
of the ray in the longitudinal direction is somewhat intriguing. But, the
positive time delay associated with the reverse lateral displacement and
the correctness of the group velocity of the wave guided by a slab de-
duced by the ray method which explicitly includes the reverse lateral
displacement confirm the correctness of the reverse lateral displacement.
[ 1 ] H. K. V. Lotsch, "Reflection and refraction of a beam of light at a
plane interface," J. Opt. Soc. Amer., vol. 58, no. 4, pp. 551-561,
Apr. 1968.
[2] H. Kogelnik and H. P. Weber, "Rays, stored energy, and power
flow in dielectric waveguides," J. Opt. Soc. Amer., vol. 64, no. 2,
EDMUND K. MILLER
In the above letter,I a technique is presented for using a circular array
of N equally spaced antennas to obtain the azimuths and strengths of
M (M < N/2) far-field sources. The procedure involves two basic steps:
1) use of a FFT (assuming N is a power of 2) to transform the N in-
dividual antenna-element strengths; 2) solution of an Mth-order linear
system composed of the transformed data to find the coefficients of an
Mth-order polynomial whose roots provide the azimuth angles.
It is the purpose of this brief letter to mention that step 2 above is
identical to Prony's Method [ 1 ] , a procedure which underlies much of
modern signal processing. Although the approach in the above letter
differs from the more direct approach presented in [1], both lead to
the same result as can be seen by comparing their respective starting
equations, which are sums of exponentials, and the subsequent linear
system and polynomial solutions. For example, the (12), in the above
letter, which has the form (redefining the Q1's to eliminate the al-
ternating minus signs of the author)
M
Wm-i= L Q1 W1-n+i,
1=1
can be seen to be the same as the linear predictor which arises in
Prony's method. In (1), the Wk come from step 1 above, and the Q1
are the linear predictor coefficients, solved for from (1), which form
the polynomial
PM+QIPM-1 +...QM= 0
whose roots provide the source directions (via their complex logarithms).
As Moody points out, basically the same technique works for finding
the azimuth angles of far-field sources from linear arrays of equally
spaced elements [2]. Furthermore, the azimuth and elevation angles
of far-field sources can be found from a similar approach using planar
arrays and receiving antennas [ 3 ]. Other applications for which Prony's
Method has been used include the synthesis of nonuniformly spaced
linear arrays [4] ; the imaging of linear source distributions from their
far fields [5] ; and the extraction of the complex resonance frequencies
(poles) of objects from the impulse or spectral responses [6), [7].
Reply2 by M. P. Moody3
I wish to thank E. K. Miller for his very useful comments and for the
list of references supplied relating to Prony's Method, of which I have
not been hitherto aware.
Although he has shown that the method derived was not, after all,
novel, it is hoped that the application of the method with deconvolu-
tion in solving the resolution problem (epoch detection) of overlapping
pulses may be found to be of use.
Since submitting the original correspondence, I have become aware of
another closely related method which, from the literature available,
does not appear to have been linked with Prony's Method. I refer to
Burg's so called maximum entropy method [ 8], [91 (also known as the
autoregressive power spectral density estimate or Linear Prediction
Spectral Estimate) and have included some recent references for the
benefit of readers in this field [10]-[12]. The method was derived by
Burg based on maximizing entropy, but results in a set of equations
(also known as modified Yule-Walker equations) which bear a re-
Manuscript received June 5, 1980.
E. K. Miller is with the Lawrence Livermore Laboratory, P.O. Box
808, Livermore, CA.
I M. P. Moody, Proc. IEEE, vol. 68, pp. 276-277, Feb. 1980.
'Manuscript received July 18, 1980.
3 M. P. Moody is with the Electrical Engineering Department, Queens-
land Institute of Technology, George Street, Brisbane 4000, Australia.
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
0018-9219/81/0100-0122$00.75 ? 1981 IEEE
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
semblance to Prony's matrix equation. Since Burg's method is equiv-
alent to modeling the data as the impulse response of an Nth-order
recursive digital filter (IIR filter) whose Z coefficients are the un-
knowns in the matrix equation [131, the positions of poles can be
found on the Z plane by solution of the resulting polynominal in Z.
Some work has also been carried out combining autoregressive with
zero-producing moving average (ARMA) processes with good results
in noisy conditions [ 141.
[1]
[21
REFERENCES
F. Hildebrand, Introduction to Numerical Analysis. New York:
McGraw-Hill, 1956, pp. 378-382.
J. M. Kelso, "Measuring the vertical angles of arrival of HF sky-
wave signals with multiple modes," Radio Sci., vol. 7, pp. 245-
250, 1972.
C. J. Benning, "Adaptive array beamforming and steering using
Prony's method," in Prot. Annu. Symp. USAF Antenna R&D,
AFAL (Univ. of Illinois, Allerton, IL), 1969.
syntheses using Prony's method," Electron. Lett., vol. 14, pp.
180-181, 1978.
E. K. Miller, D. L. Lager, and J. T. Okada, "Imaging of linear
source distributions," Lawrence Livermore Laboratory, Rep.
UCRL-52822, Nov. 1979.
-Evaluation or a processing technique for transient data," IEEE
Trans. Antennas Propagat., vol. AP-26, pp. 165-173, Jan. 1978.
J. N. Brittingham, E. K. Miller, and J. L. Willows, "Pole extrac-
tion from real-frequency information," Proc. IEEE, vol. 68, pp.
263-273, Feb. 1980.
[91
[10]
sented at the 37th Anna. Meet. Society of Exploration Geolo-
gists 4Oklahoma City, OK), Oct. 31, 1967.
Maximum entropy spectral analysis," Ph.D. dissertation,
Stanford University, Stanford, CA, 1975.
D. N. Swingler, "A comparison between Burg's maximum entropy
method and a nonrecursive technique for the spectral analysis of
deterministic signals," J. Geophys. Res., vol. 84, no. B2, pp. 679-
685, Feb. 1979.
[11] S. Holm and J. M. Hovem, "Estimation of scalar ocean wave
spectra by the maximum entropy method," IEEE J. Oceanic
Eng., vol. Oe-4, pp. '76-83, July 1979.
[12] S. M. Kay, "The effects of noise on the autoregressive spectral
estimator," IEEE Trans. Acoust., Speech and Signal Processing,
vol. ASSP-27, pp. 478-485, Oct. 1979.
[13] A. van den Bosr "Alternative interpretation of maximum entropy
spectral analysis," IEEE Trans. Inform. Theory (Corres.), vol.
IT-17, pp. 493-494, July 1971.
[14] M. Kaveh, "High resolution spectral estimation for noisy signals,"
IEEE Trans. Acoustics, Speech and Signal Processing (Corres),
vol. ASSP-27, pp. 286-288, June 1979.
Electrostatic Fields Inside Two Planar Distributions
of Potential
Abstract-The potentials and fields inside two parallel planar distribu-
tions of potential can be written as Green's function integrals. Each
Green's function can be evaluated numerically by decomposing it into
two series of terms (one for large and one for small distances) which
converge rapidly within their intervals of use. The accuracy of trunca-
tion after a few terms is evaluated.
1. INTRODUCTION
The calculation of electron trajectories in electron devices requires
knowledge of interior electric fields determined by boundary voltages.
In many such devices, the shape of the boundary voltage plates is too
complicated for analytical determination of the interior fields. With
the advent of digital electronic computers, it has become feasible to
generate accurate numerical solutions to potentials and fields in the
interior of electron devices, but the cost of generating enough field
points, or of regenerating trajectory points when needed, for calculating
electron trajectories can increase rapidly with increasing volume. For
many applications, especially for devices which approximate a sandwich
of planar distributions of boundary voltages over large areas, the cost
can become prohibitive. In this work, we propose a semianalytical,
seminumerical technique to overcome this difficulty. We shall assume
that two parallel planar distributions of voltage form the boundary (if
necessary, two planar grids of voltages can be generated numerically to
Manuscript received March 28, 1980; revised August 18, 1980.
The author is with RCA Laboratories, Princeton, NJ 08540.
bound a convenient part of the interior). Unlike a previous communi-
cation [ 1 ] , the two voltage distributions may be different.
A Green's function connecting an arbitrary interior point to all grid
points of the two planes can be calculated for interior potential and
fields. It can be expressed as either of two series expansions (one for
large, the other for small distances). Only a few terms of each are
needed for reasonable accuracy. The problem is essentially reduced to
a two-dimensional integral over the two planar boundary distributions.
II. ANALYSIS
Let the two bounding plates by located at yi = ? Y (plus sign for
i = 1, minus sign for i = 2). Let Vi(xo, zo) be the voltage distribution
on plate i. Consider an interior point r = (x, y, z) and define Ri =
(x - x0, y - yi, z - zo) as the distance vector between r and a point on
plate i. The methods of [ 11 then yield for potential 0 (r),
2
'D (r) Jfdxo dzo Vi (xo, zo) G (Ri)
i=1
1 sinh [K (Y ? y)]
G(R1)= 27r 0 dKKJ0(KRp) (1)
sink (2KY)
where Rp = [(x -.x0)2 + (z - z0)2] 1/2 is the parallel component of Ri.
Further development of the Green's function G (Ri) can be pursued, as
in [I I, to obtain,
G(R')
2 YZ [G1(, n) - G2(E, n)]
G1(l:,n)
{(4m+1
+-n)/[(4m+1
+n)2
+12]3/2}
m=0
G2(,)
- :
{(4m + 3
? n)/[(4m + 3
? n)2
+ 1.2 ] 3/2 }
(2)
1n=0
where t = [(x - 1:0)2 + (z - z0)2] 1/2/Y and n = y/Y are dimension-
less versions of variables Rpp and y, and where we retain the convention
that the upper sign holds for i = 1 and the lower for i = 2. The con-
vergence of (2) requires at least as many terms as are given by the
integer closest to the ratio Rp/4Y. Hence distant source points require
relatively many terms in the Green's function. It is desirable to supple-
ment (2) by an expansion that handles distant source points more
economically. Fortunately an alternative exists. Rewrite G(Ri) as,
[cosh (nt) sinh (nt)1
G(Ri)=-- dttJ0Qt) I ? J. (3)
4"Y2 J o L cosh (t) sinh (t)
Utilizing the fact that the Bessel function of the first kind and zero
order J0(t) can be written as the real part of a Hankel function that
goes to zero sufficiently rapidly when t -> o, we may deform the
integration in the complex t plane to one along the positive imaginary
axis that contains infinitesimal half circles around the singularities of
1/cosh (t) and 1/sinh (t).1 The resulting series is
_ 1
G(Rt) -2 [Gc(f, n) - Gs(, n)]
27r y
Gc(,n)=ir 2- (-1)'n(m+i)cos[(m+1)irn]K0[(+n+2)7rf]
m==0
Gs(>;, n) = ?7r E (-1)'m sin (mrrn)Ko(mat) (4)
m=o
with the upper sign. for i = 1 and the lower for i = 2. The smallest argu-
1 It may be helpful to bear in mind that the real part of an integral of
a complex function along the positive imaginary axis results. The infin-
itesimal half circles yield real contributions, whereas the connecting
pieces along the imaginary axis only yield purely imaginary contribu-
tions to the integral. Hence only the former contribute to the real part
of the integral.
Approved For PMp#91i? 169 6161812~.W~~WPWR000100280003-6
111
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
TABLE I
POTENTIAL
1
0.5
1
1.5
2
A
.04709
.03817
.02227
.01076
.00478
-0.5
E
.04735
.03843
.02232
.01076
.00478
A
.08358
.06420
.03402
.01522
.00648
-0.25
B
08397
06459
03399
01522
00648
.
.
.
.
_
A
.1452
.1006
.04519
.01820
.00735 1
0
E
.1458
.1011
.04520
.01820
'00735
A
.2727
.1532
.05305
.01873
.00713
0.25
E
.2734
.1539
.05309
.01873
.00713
A
.6297
.2184
.05178
.01582
.00569
0.5
E
.6305
.2192
.05173
.01582
.00569
TABLE II
HORIZONTAL FIELD
0
0.5
1
1.5
2
A
0
.0308
.0290
.0168
.07952
0.5
E
0
.0308
.0294
.0168
.07452
A
0
.0638
.0513
.0255
.01117
0.25
E
0
.0638
.0511_
.0255
.01117
A
0
.1342
.0605
.0333
.01320
0
E
0
.1342
.0804
.0333
.01320
A
0
.2993
.1142
.0375
.01335
0.25
E
0
.2993
.1147
.0375
.01335
A
0
.6743
.1356
.0344
.01105
0.5
E
0
.6743
.1350
.0344
.01105
ment of a Bessel function in (4) is 7rt/2 - 1.6 t. Since K0(x) goes as
x-1/2 exp (-x) for large x, we find that the sum of all these terms con-
verges rapidly with increasing m when t > 1. Hence (4) is the desired
alternative to (2) for large C.
Similar expansions for the electric field are easily worked out [E(r) is
the negative gradient of the potential 4) (r)] . Each component of the
electric field defines a new Green's function,
1
Gy (R1) 2rrY3 g,7(t, r7)
(x - xO)
Gx(Ri) = 2rnY4 g(, rl)
(Z - ZO)
Gz (RI) = 27rY4 g (, TI)
and we shall simply give the expansions of gt and gn for t < 1 and i; > 1.
Small ~(t < 1):
112-2(4m+1-++71)2 12-2(4m+3?n)2
gQ.n) +
m=0 [(4m + 1 + 1)2 +2]5/2 [(4m + 3 ? n)2 + X21 5/2
gt(t,'n)=3y- j[(4m+1+77)2+ 215/2
m=0 l
4m+3?s
[ (4m + 3 ? +7)2 + l2 ] 5/2 (6)
TABLE III
VERTICAL FIELD
~
0
0.5
1
1.5
2
A
-.1184
-.0902
-.0461
-.01974
-.00812
-0.5
-.1189
-.0907
-0462
-.01974
-.00812
A
-.1822
-.1212
-.0471
-.01548
-.00535
-0.25
E
1827
-
1217
-
-.0468
-.01548
=.00535
.
.
A
-.3343
-.1744
-.0406
-.00775
-.00143
0
E
-.3348
-.1749
-.0410
_00775
-.00143
A
-7683
-.2471
.0184
.00423
.00331
0.25
E
-.7688
-.2476
.0182
.00423
.00331
A
L
-.2559
-.2369
.0346
.01934
.00810
0.5
E
-.2559
.2375
.0347
.01934
.00810
gt(t,7n)=7r2 (-1)m {(m+2)2 cos [(m+2)7rr7]KI[m+2)7rl]
m=O
+ m2 sin (mrrrl)KI(mrrt)}.
The sign convention is retained in (6) and (7): upper sign for i = 1 and
lower sign for i = 2.
Calculations of the three Green's functions of (4) and (5) for y = 1
are shown in Tables I-III in the rows marked "E" for values of t
and for 1771 5 ? The rows marked "A" give approximations for the
m < 3 partial sums in (4), (6), and (7). The description of this half
space 1711 < 2 by just a few terms appears to be accurate to well within
a percent.
REFERENCE
[1 ] D. A. de Wolf, "A useful potential solution V(x, y, z) for sym-
metric rectangular stepwise constant-potential boundaries at
y = ?Y," Proc. IEEE, vol. 66, pp. 85-86, Jan. 1978.
Correction to "Criteria for the Separation of Real and
Complex Natural Modes of Dynamical Systems"
In the above letter,1 (15) and (16) contain misprints. The correct
version is
II[C-R rCII2j411\CLfII.
Manuscript received June 6, 1980.
The author is a Visiting Professor with the Faculty of Engineering,
Concordia University, and McGill University, Montreal, P.Q., Canada.
1 F. M. Reza, Proc. IEEE, vol. 68, pp. 174-175, Jan. 1980.
Frequency Domain Least-Mean-Square Algorithm
S. SHANKAR NARAYAN AND A. M. PETERSON.
Abstract-Frequency domain adaptive filtering can be performed by
Fourier transforming the input-signal vector and weighting the contents
of each frequency bin. By reducing the eigenvalue spread of the data
autocorrelation matrix, frequency domain filtering promises great
improvements in convergence rate over the conventional time-domain
adaptive filtering.
Adaptive filters are being widely used in a variety of applications,
including noise cancellation [1], linear prediction filtering [2], and
spectral estimation [3] . Filters of this type are generally implemented
in the time domain in tapped-delay-line (TDL) form, and the Widrow-
Hoff adaptive least-mean-square (LMS) algorithm [4] is used to obtain
the filter parameters. The purpose of this letter is to describe a fre-
quency domain equivalent of the time domain adaptive filter. This
g77(5,+i)=rr2 (1)'n {(m+2)2 sin [(m+2)7rnlKo[(m+2)17 l
m=0
? m2 cos (m7r17)K0(mrrt)}
Manuscript received August 21, 1980.
The authors are with the Center for Radar Astronomy, Department
of Electrical Engineering, Stanford University, Stanford, CA 94305.
Approved For F gS1-12 ~o/b 1b12 1 P8WAbMR000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
INPUT
SIGNAL
(Knl
BAND PASS DIGITAL COMB
FILTERS. IMPLEMENTED
BY OFT
COMPLEX
ADAPTIVE LMS
ALGORITHM
Fig. 1. A conventional TDL adaptive filter.
method promises great improvements in convergence rate over the time
domain approach.
The general form of the TDL adaptive filter [4] is shown in Fig. 1.
Its function is to weight and sum the delayed samples of the input
signal xn to form an adaptive output. The input signal vector Xn and
the weight vector An at the time of nth iteration are defined as
Xn = [xn xn-J ... xn-(N-1)] T
An = [ano an1 ... an(N-1)] T (2)
respectively. The corresponding filter output is
Yn = X nAn
The error en is the difference between the desired response do and
the output yn
en=dn -yn. (4)
In the LMS adaptive algorithm, the mean-square error (en) is mini-
mized by continuously updating the weight vector An as each new data
sample is received, according to the expression
An-1 = An + 2,uXn?n (5)
where p is a positive constant that governs the rate of convergence, and
proper choice of this constant ensures stability of the adaptive process.
It can be shown [5] that for stationary input and sufficiently small ?,
the speed of convergence of the algorithm is dependent on the ratio
of maximum to minimum eigenvalues, \max/Nmin, of the data auto-
correlation matrix RXX,
RXX = E(XnXn). (6)
Slow convergence rate can, be expected when this ratio is large. This
eigenvalue spread is bounded by the relation [2],
max I X (e] `'')12
Amax 06wou; rev1sea neptember 10, 1980. ^ ^^
The author is with H. R. .-Singer, Science Park Road, State College, h = A P (8)
PA 16801. where (hi) is the impulse response given by the Maclaurin series expan-
0018-9219/81/0100-0126$00.75 ?1981 IEEE
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
N-1
H(Z) = 57 hiZ-i. (9)
1=0
Equation (8) can be solved by the least squares normal equations given
by
P=(ATA)-- ATh. (10)
As is well known, the ordinary least squares estimators of P of the gen-
eralized linear regression model are unbiased and consistent. In addition,
the maximum likelihood estimators of the regression coefficients, P, are
equivalent to the least squares estimators.
Equation (10) can be readily solved for r leading to the following
equations for b 1 and b 2.
N-2 N-3
br k3h1 +ks Y hihi+I +k6 Y hihi+2
i=0 i=0
N-2
b2 -- k4hi +k6 > hihi+l +k
i=0
k3 =_ [h0 L hi lID
{13 ~+
k4 = Y hihi+r D
FN-3 1
ks = ! hIJ D
k6 = hihi+l] (L)
rN-2 l
k7= I hl -h'J D
Dt htj [ haj - [L hihi+I] - [hi2] ha. (12)
III. COMPUTATIONAL CONSIDERATIONS
Equations (11) and (12) have been expressed in a form that is suitable
for exposition but not for computation. It can be observed that many
variables are simple functions of other, already computed, variables.
For instance, if EN- 3 hI has already been computed, then
N-2 N-3
h? = hi + hN2-2
i=0 i=0
which involves just a single multiplication and addition.
If a measure of an algorithm's complexity is the total number of
operations involved, then the complexity of this algorithm is linear in
N. Depending on the method used to compute operations, the total
number of operations is given by K * N where K varies between 3 and 5
and N is the length of the input data record. For the DFT using the
FFT, the number of computations varies as N log2 N. For large data
records (N > 32) the autoregressive-moving average (ARMA) method
has an advantage.
IV. EXPERIMENTAL RESULTS
A computer program was written to test the accuracy of three methods
of estimating a pure tone. The methods are, DFT processing (11, the
method of moments [4] and the method presented in this paper. A
pure sinusoidal was generated located approximately at the center of
the normalized frequency spectrum (frequency normalized with respect
to the sampling rate.) The sinusoid was added to a Gaussian random
S/S IN DECIBELS
Fig. 1. Error performance of three frequency estimators. (a) ARMA
estimation (b) DFT estimation (c) method of moments.
process of variable intensity. The error was computed as the magni-
tude of the difference between the true and estimated frequency
divided by the true frequency. The results of this experiment are shown
in Fig. 1. It can be noticed that on a logarithmic scale the error is linear
for the ARMA method. In addition, the errors for large and small SIN
ratios is smaller for the ARMA than for the DFT method.
V. CONCLUSION
A method has been presented which can efficiently and accurately
estimate the peak frequency in an arbitrary spectrum. Simulated test-
ing has shown the viability of the method.
Many alterations to the ARMA method of frequency estimation can
be suggested. One modification is to use a weighted least squares
formulation where the weight are determined by the errors of the
estimates. The estimation is then a two step process.
Alternately, following the development of [7], the problem can be
solved recursively. The frequency can then be estimated adaptively
permitting a continuous tracking of the peak frequency.
REFERENCES
D. C. Rife and G. A. Vincent, "Use of the discrete Fourier trans-
form in the measurement of frequencies and lends of tones,"
Bell Syst. Tech. J., vol. 49, pp. 197-228, Feb. 1970.
[2] R. B. Blackman and T. W. Tukey, The Measurement of Power
Spectra from the Point of View of Communications Engineering.
New York: Dover, 1959.
L. C. Palmer, "Course frequency estimation using the discrete
Fourier transform," IEEE Trans. Inform. Theory, vol. IT-20, pp.
104-109, Jan. 1974.
to tone detection," Bell Syst. Tech. J., vol. 55, pp. 143-155, Feb.
1976.
J. P. Burg, "Maximum entropy spectral analysis," presented at the
37th Meet. Soc. Exploration Geophysicists, Oklahoma City, OK,
Oct. 31, 1967.
"Predictive filtering and smoothing of short records by using
maximum entropy," J. Geophys. Res., vol. 78, pp. 4959-4964,
1973.
C. R. Guarino, "Adaptive signal processing using FIR and IIR fil-
ters," Proc. IEEE, pp. 957-9 59, June 1979.
Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
 Approved For Release 2006/03/10 : CIA-RDP88B00553R000100280003-6
PROCEEDINGS OF THE IEEE, VOL. 69, NO. 1, JANUARY 1981
A Method for Improving the Classification Speed
of Clustering Algorithms Which Use a
Euclidean Distance Metric
Abstract-Many pattern recognition computer programs use one of
the clustering algorithm techniques. Often these algorithms use a
Euclidean distance metric as a similarity measure.
A scheme is proposed where both the Euclidean metric and a more
simple city-block metric are utilized together to reduce overall classifi-
cation time. The relation between the Euclidean and city-block dis-
tances is introduced as a scalar function. The bounds of the function
are given and used to decide whether classification of each pattern
vector is to be achieved by the computationally slow Euclidean distance
or the faster city-block distance. The criteria is that the classification
should be identical to the original Euclidean only scheme.
1. INTRODUCTION
Given a specific n-dimensional pattern vector, one of the important
operations required of a clustering algorithm is to select from the K
class references, the class reference which is closest to the pattern. The
closeness can be expressed by one of the many forms of distance
metric [1) .
Many clustering algorithms produce optimum results if the Euclidean
metric is used [2]. This metric can be expressed as
E(X, Z1) = [(X - Zj)T(X _Z j)] 112
where X represents the n-dimensional vector to be classified and Zj
represents the /th class reference point. The task is to compute and
select
MIN E(X, Zj), / = 1, 2, ? - ? , K.
For each pattern vector there are K computations.
The unweighted city-block distance metric can be computed more
quickly than the Euclidean metric because no multiplications are in-
volved (this is especially true for computers using software multiplica-
tion). The unweighted city-block metric is given by
n
C(X,Zj)=T IXk - ZjkI
k=1
where k is the kth component of the n-dimensional vector.
However, for a scheme which is optimum using a Euclidean metric,
inferior classification accuracy would result if a city-block metric was
used, albeit with an improvement in classification speed.
Under certain conditions on the data, the minimum of the K Euclid-
ean distances corresponds to the minimum of the K city-block dis-
tances. Under these conditions identical classification of the pattern
vector is produced for both measures. The fact that the magnitudes of
the two distances may be different does not affect classification.
If it could be deduced for a particular pattern vector, that the small-
est of the set of Euclidean distances corresponded to the smallest of the
set of city-block distances, then the city-block distances could be com-
puted instead of the Euclidean and a saving in classification time would
result.
A method of deciding which one of the distances to employ is given
below.
11. A MINIMUM EUCLIDEAN DISTANCE FROM A CITY-BLOCK RATIO
Consider two city-block distances, designated Cl and C2. A logic
decision is to be made to either use the smaller of Cl and C2 to obtain
a classification, or to reject them, and compute and utilise the smallest
Euclidean distance. Lemmas 1 and 2 will be given without proof since
this is beyond the scope of this letter.
Lemma 1: The relationship between the Euclidean and city-block
metrics is given by
C1= Ej ? gj(9n) (1)
Manuscript received June 11, 1980.
J. D. Curie is with Plessey Electronic Systems Research, Eastleigh
Road, Havant, Hampshire P09 2PE, England.
J. J. Hill is with the Department of Electronic Engineering, University
of Hull, Hull, HU6 7RX England.
where gj (en) is a scalar function of the orientation of Ej with respect to
the axes of the n-dimensional space containing X and Z,.
Lemma 2: gj (Bn) is a bounded function
I