- Email: [email protected]
Neighborhood operation binary image algebra for optical morphological image processing
I February1996 OPTICS COMMUNICATIONS ELSEWIER
Optics Communications
123 ( 1996 ) 705-7
15
Full length article
Neighborhood operation binary image algebra for optical morphological image processing Shifu Yuan, Guofan Jin, Minxian Wu, Yingbai Yan Tsinghua University, Department
of Precision Instruments, Beijing 100 084, Chinu
Received 27 March 1995
Abstract
A neighborhood operation binary image algebra (NOBIA) which has only one basic operation is presented and developed. The basic operation of NOBIA is a convolution followed by a nonlinear filtering function and then an intersection operation. The parallel architecture of an optical neighborhood operation digital image processor is designed to efficiently perform the parallel morphological image processing algorithm provided by NOBIA, and the optical implementation hardware of the processor is discussed. Using an incoherent optical convoluter as a three-dimensional free space interconnection device, a smart LCLV as an optical nonlinear device, and a simple optical circuit for intersection operation, an optical hardware is constituted and used to realize the basic operation of NOBIA experimentally. Any image transformation can be performed with this hardware by executing the basic operation repeatedly according to proper iterative programming.
1. Introduction With optical parallel architecture a massive number of optical computing systems have been proposed. The most important and efficient application of optical computing is in the field of image processing, and many structures of parallel image algebra suited for optical implementation have been presented. Mathematical morphology [ 1,2] is a widely-used image algebra in which a broad class of image processing operations are represented by morphological transformations (MT) with the combination of dilations and erosions. On the basis of mathematical morphology, Huang et al. have developed a binary image algebra (BIA) , in which any image transformation can be implemented by the use of appropriate reference images and the three fundamental operations of complement, union and dilation [ 3,4]. Along this direction, an image logic algebra (ILA) has been proposed which has more operations [ 5,6]. Recently, Liu et al. have proposed another parallel binary image algebra, called one-operation image algebra (OIA), with a logic operation followed by a dilation [ 71. This paper presents our attempts at searching a simpler algebraic representation for optical implementation of parallel image processing. According to BIA, three hardwares are needed to respectively perform the three fundamental operations. In OIA, an operation is constructed with a logic operation followed by a dilation. In fact, if we can construct a basic operation that can be programmed to perform at least the three fundamental operations of BIA, another binary image algebra can be developed with this basic operation and thus all of the image transformation can be represented by such an image algebra. The dilation of a binary image X with a reference image R can be performed by the convolution of images X and R followed by a thresholding [ 81. We have defined a neighborhood 0030-4018/96/$12.00
0 1996 Elsevier Science B.V. All rights reserved
SSDlOO30-40 I8( 95 )003 17-7
706
S. Yuan et al. /Optics
Communications
123 (1996) 705-715
operation as a convolution followed by a nonlinear filtering with a point nonlinear function f9]. The combination of a neighborhood operation and an intersection operation can be regarded as a basic operation. The intersection, complement and dilation can be efficiently performed with such a basic operation by selecting different formats of the point nonlinear function. Based on the basic operation, we develop a new binary image algebra that is called neighborhood operation binary image algebra (NOBIA) . All morphological image transformations such as erosion, opening, closing and hit-miss transform can be represented easily with the NOBIA. Thus NOBIA provides a simple parallel algorithm for optical morphological image processing. Also, we construct a simple processor, i.e. the optical neighborhood operation binary image processor (ONODIP) , to perform the algorithm of NOBIA. The suggested ONODIP is in fact a single instruction multidata stream (SIMD) type of parallel processor. The ONODIP has the advantage of the inherent parallelism of optics. Because optical free space interconnection is adopted, the processor has a very flexible interconnection style, thus the structuring elements for morphological image processing can not only be very large, but also be altered easily.
2. Neighborhood
operation binary image algebra (NOBIA)
The binary image algebra (BIA) , extending from mathematical morphology, is the basis of NOBIA. For convenience, we here give a short review of BIA. In BIA, a binary image is regarded as a set of image points (pixels with value 1). The domain of images is defined as the universal image, which is the set W= ((x, y) IxEZ,,~EZ,), where Z,= {0, f 1, . . . . kn), and n is a positive integer. The image space S is the power set (the set of all subsets) of the universal image, i.e., S = P( W). Any image is a subset of the universal image. The three fundamental operations of BIA are defined as follows: (1) Complement of an image X:
rz=1(x, Y> I(.& y) E WA (4 Y) e-v .
(1)
(2) Union of two images X and R: Xu R= ((x, y> I (x, y) EXV (x, y) ERJ .
(2)
(3) Dilation of two images X and R: XOR=
{xi +x2,
y1 fy,)
E
WI (x,, Y,) =f,
0,
(xzt ~2) ER
(X+0) A (R+0) I (X=0)
V (R=0)
7
(3)
.
Remarks: E denotes belong to, A denotes AND, V denotes OR, and (21is the null image having no image point. The erosion of an image X by a structuring element R can be written as XeR=%ri,
(4)
where the reflection image R is defined as i={(-.x,
-y)I(x,y)ER].
(5)
The morphological dilation and erosion are in fact a convolution followed by a point nonlinear filtering in the spatial domain [ 81. We define a neighborhood processing operation NP(X, R) as the convolution of a binary image X and a binary image R followed by a point nonlinear function [ 91, i.e., NP(X, R) =NL[ U(X, R)] ,
(6)
where NL(LI) denotes a point nonlinear filtering function, and u is the gray level of a pixel in the convolution resultant image U(X, R) _The image V(X, R) can be expressed as U(X, R) =X*R,
(7)
S. Yunn et al. /Optics
I
Communications
707
I23 (1996) 705-715
_____ L._--.--~ 1
2
N
u
(a) Fig, I. Nonlinear functions function NL’( I( ).
for the neighborhood
operation.
(a) The first-type nonlinear
function
NL’( u). (b) The second-type
nonlinear
where * denotes convolution. The value of pixel (p, q) in the image X can be expressed as x(p, q) which has two possible states 0 and 1, where p=o, +1 1 . . .. 5 II and q = 0, f 1, . . , 5 n. The value of pixel (j, k) in the structuring element R is expressed as r(j, k) which has two possible states 0 and 1, where j= 0, + 1, . .., +n and k=O, * 1, . ., f R. According to Eq. (7), the value of a pixel (p, q) in the image U(X, R) is expressed as U(P. 4) =
C C-Q-j, I
s-k)r(.L
k)
(8)
T
h
where.r(p, q) =0 when ]p] >n or 141 >n. Dilation and the complement of dilation can be performed directly by NP(X, R) with different formats of the point nonlinear function. Figs. la and b respectively give the point nonlinear functions for the dilation and its complement. The function shown in Fig. la is called the first-type nonlinear function and is expressed as NL’( u), while the function shown in Fig. 1b is called the second-type nonlinear function and is expressed as NL*( U) . Then the neighboring processing function with the nth nonlinear function is represented as NP”(X. R) =NL”[ U(X, R)] =NL”[X*R] where II = written as
1,
,
2. The relations between the neighboring
(9) processing
function and the morphological
dilation can be
NP’(X. R) =NL’[
U(X, R)] =XOR,
( 10)
NP’(X. R) =NL’[
U(X, R)] =XOR.
(11)
When we select the reference image R to be an element image Z= { (0, 0) } that has only one image point at the origin, Eqs. ( 10) and ( 11) respectively give the original image X and its complement image X, i.e., NP’ (X. I) = NL’ [ U(X, I) ] =X ,
(12)
NP’(X, I) =NL’[
(13)
U(X, Z)] =X.
According to Eqs. ( 10) and ( 13), the neighborhood processing function NP(X, R) can perform dilation and complement operations. In BIA, the three operations (complement, union and dilation) are chosen as the three elements of the complete set of fundamental operations. Because union can be performed with complement and intersection, the three fundamental operations of the complete set can also be complement, intersection and dilation. Thus the question is how to construct a basic operation that is capable of performing at least complement, intersection and dilation. The basic operation F,,( . ) of the neighborhood operation binary image algebra is then defined as the neighborhood processing function NP(X, R) followed by an intersection, i.e., F,,{X,R,NF”‘,
Y)=NP”(X,R)nY=NL”[X*R]nY,
(14)
S. Yuan et al. /Optics
Communications
:.
, : Binary
- Gray Image Input -
Binary Image
: ;
Nonlinear
Convolution
i Image
Filtering
x
I
j i
123 (1996) 705-715
Image R
Neighborhood
:
Function NL(u)
j
Processing
Binary Image
I
I Nonlinear
Reference
. Output Intersection
Binary Image Y
Operation
Fig. 2. Structure frame of the basic operation of NOBIA.
where n = 1, 2. X and Y are two input binary images, R is a reference image. The structure of the basic operation is shown in Fig. 2. Selecting different forms of the nonlinear function, the basic operation F,,( . } has different forms, i.e., F, (X, R, NP’, Y) =NP’(X,
R) n Y= (XOR)
17 Y,
(15)
F2 (X, R, NP*, Y) =NP’(X,
R) I’-’Y= (XOR)
n Y.
(16)
Complement,
intersection
and dilation can all be performed by the operation F,( . ) :
Fz(X,Z,NP2,
W)=NP’(X,Z)nW=X,
(17)
F,(X,Z,NP’,
Y)=NP(x,z)nY=xnr,
(18)
F, (X, R, NP’, W) =NP’(X,
R) f7 W=XOR,
(19)
where W is the universal image. Because any image transform can be implemented by the iteration of complement, intersection and dilation with appropriate reference images, it can also be represented by the iteration of the basic operation F,( . ) . On this basis, we define NOBIA as an algebra with an image space S= P( W) and a basic operation F,,( . }, which is represented as NOBIA=
(P(W),
3. NOBIA representation All morphological
(20)
F,,) .
of morphological
image transformations
image transformations can be derived from F,,( . ). Here we give some examples:
(1) Union of two images X and Y: XUY=XnY=F,(F,(X,Z,NP2,
F,(Y,Z,
NP*, W}),Z,NP*,
W) .
(21)
(2) Difference of X by K XIY=XnY=F,(Y,Z, (3)
Erosion of image X with a reference image R: XeR=F,(F,(X,Z,
(4)
(22)
NP2,X).
Symmetric
difference
NP2, W),d,
NP2, W) .
of two images X and Y:
(23)
S. Yuan et al. /Optics Communications 123 (1996) 705-715
709
XAY= (X/Y) u (Y/X) = F, ( Y, I, NP’, X) U F2 {X, I, NP’, Y) = F2 ( F2 ( F2 { Y, Z, NP*, X), I, NP’, Fz ( F2 {X, I, NP’, Y), I, NP*, W} }, I, NP*, W) . (24)
(5)
Opening of image X by a reference image R: XoR==(XeR)OR=F,{F,(F,{X,
NP’, I, W},k,NP2,
W), R,NP’,
W).
(25)
(6) Closing of X by a reference image R: X.R=(XOR) (7)
R, NP’, W), li, NP*, W) .
(26)
transform of image X by an image pair R( Rn Rb) :
Hit-miss XOR=
eR=F2(F2{X,
(XeR,-)
fl (&R,)
=J?@Z?,~IX@Z&,
= F, (F, (X, I, NP’, W}, dr, NP*, F2 (X, d,, NP2, W) } . (8)
Edge detection. The out-object
edge detection can be denoted as
ED,,(X) = (XOR)/X=F,{X,
Z, NP2, F, {X, R, NP1, W}} .
(27)
(28)
The in-object edge detection can be denoted as ED,(X) =X/(XBR)
= (%Ori) nX=F,
{F,(X,
Z, NP2, W}, d, NP’, X) .
(29)
The central edge detection can be denoted as ED,(X) (9)
= (XOR)
nXBR=F,
(F2{X, Z, NP2, W), d, NP’, F, (X, R, NP’, W)} .
Shifting an image X from the coordinate X@((i,j)l=F,(X,
( 10) Noise remoual.
(30)
(x, y) to (x + i, y -kj) :
l(i,j)),NP’,WJ.
(31)
The removal of noise points in an image X can be denoted as
NR(X) =XoR@R =F2(Fz(F,(F2{F2{X,NP2,Z,
W},ri,NP2,
W), R,NP’,
W}, R, NP’, W},i,NP’,
W}
(32)
= F, [ F2 { F2 ( F2 { F2 (X, R, NP2, W}, i, NP*, W), NP’, Z, WI, I?, NP’, W), R, NP’, W}
(33)
or NR(X)X.RoR
where R is a small structuring element. ( 1 I ) Skeltonization. In mathematical morphology, SK(X)
=
the skeleton of an object X is expressed as
U S,(X), oCn
S,(X) = (XenR)
n (XenR)oR,
wheren=1,2,...,N,andN=MAX{kJX8kR#0}.InNOBIA,itcanbedenotedas
710
S. Yuan et al. /Optics
SK(X) = f-l
f-l
J,(X) = Fz (
O<“=GN
J,(x)
J,(X),
Z, NP’, W) ,
1, IQ”
, Jz(x)
I23 (1996) 705-7I5
(35)
O
= FI
J,(X) =S,(X)
f-l
Communications
{F,
1..
J’,
= ((X8nR)
{.f,
(x>,
1,
. . ., 1,
NP’, JN- ,(x) }, 1, NP’, JN(x)
n (X8nR) eR0R=X0nzh--lX0(n+
wheren=l,2
4. Comparison
,...,
Z, NP’, W), (n+l)ii,
(364
l)I;i0R
= F2 { F2 ( F2 {X, Z, NP*, W}, ni, NP*, F2(F2[F2{X,
} ,
(36b) NP’, W), R, NP*, W)}, Z, NP’, W),
N,andN=MAX{k~F2(F2(X,Z,NP2,W),kZ?,NP2,W}#0}.
of NOBIA with BIA and OIA
In comparison with the existing image algebra, NOBIA has an obvious advantage in the parallel algorithm of image processing. In NOBIA, there is only one basic operation which has dual processing modes by selecting two nonlinear functions. Then only one unit of optical hardware is required by NOBIA, and the control and programming of the processing system can be very simple. But in BIA there must be three hardware units to respectively perform three fundamental operations. An efficient morphological image processing structure should have less basic operations and less iteration steps for a morphological transform. NOBIA has the smallest number of basic operations (only one) ; furthermore, the suggested NOBIA may also simplify greatly the iteration steps of morphological image processing. For example, a closing operation requires four iteration steps in BIA, but it needs only two steps in NOBIA; the hit-miss transform needs six iteration steps in BIA, while in NOBIA it needs only three steps. We also compare the iteration steps of other morphological transforms in BIA with that in NOBIA, and this gives the comparing results in Table 1. NOBIA also has obvious advantages comparing with OIA, though there is also only one operation in OIA. The basic operation of OIA has at least 16 processing modes (for 16 binary logic functions) which must be programmed to select. Moreover, the complete system of OIA has two hardware units because the basic operation is divided into two parts (one for logic operations, the other for dilation) by the image output port. Then the control and programming will be very complicated. But in NOBIA there is truly one hardware unit, and it only requires to programmably select the formats of a nonlinear function. Table 1 The comparison
results of iteration steps in NOBIA and BIA Image transform
Iteration steps in BIA
Iteration steps in NOBIA
complement union intersection dilation erosion difference symmetric difference opening closing hit-miss transform central edge detection
1
1
I
3
4 I 3 3 7 4 4 6 7
I 1 2
I 5 3 2 3 3
S. Yuan et al. /Oprics Communications 123 (1996) 705-715
711
Control Unit Image Data Flow CxSignal
Clock
Flow
Program Counter
F,(‘} Basic Operation
c
Test and Feedback Controller
:
Image Output
DeVice
Parallel Processing Unit Fig. 3. Architecture of the ONODIP 5. Optical neighborhood
operation digital image processor (ONODIP): Architecture
NOBIA provides a parallel algorithm for morphological image processing. The architecture for implementing the algorithm should be parallel. Fig. 3 shows the architecture of the optical neighborhood digital image processor (ONODIP) for executing the parallel algorithms of NOBIA. ONODIP is in fact a cellular SIMD machine and consists of a parallel processing unit under the supervision of a control unit. The control unit includes a clock, a program counter, a test and feedback controller, and an instruction memory. The parallel processing unit includes an image memory device for storing multiple images, three frame buffers, and a central processing unit (CPU) for executing the basic operation F,{ . ). The three buffers X, R and Y are respectively used to store the three input images for the central processing unit. The images stored in the image memory device can be respectively sent to the three frame buffers. The CPU can be controlled by the controlling unit to perform proper functions. The ONODIP operates as follows: ( 1) the corresponding images are sent to the frame buffers X, R and Y from the image memory device; (2) the basic operation F,{ ) is executed by selecting the proper formats of the nonlinear function; (3) the result is output, tested for program control, or sent back to the image memory device. Then the system can be fed back to step ( 1) .
6. Optical implementation
hardware of the ONODIP
The optical implementation hardware of the ONODIP is the CPU for optically performing the basic operation F,,( ) . According to Eq. ( 14), we need an optical convoluter, an optical spatial nonlinear device with two operation modes, and an optical parallel intersection processor.
712
S. Yuan et al. /Optics
Communications 123 (1996) 705-715
ML
Pl
p2
P”
Fig. 4. An incoherent optical convoluter.
6.1. Incoherent optical convoluter The convolution can be performed with the incoherent optical convoluter [9] as shown in Fig. 4. It is a generalization of the incoherent optical correlator which has been used to implement morphological image processing
1101. The input plane P, is illuminated by a uniform incoherent light. A mask M is put close in front of an imaging lens L with a focal lengthf. The distance between plane P, and mask M is D. The distance between Mask M and lens L is very small and can be neglected. The mask M is used to display the structuring element image R. The shortest distance between two near pixels (called pixel sampling period) in plane P, is a. The pixel sampling period in mask M is d. Plane P, is the imaging plane in which a clear image of the input plane Pi can be obtained. The distance v between lens L and plane P y can be written as
Df
(37)
ll=o-f’
The distance between lens L and output plane P, isf+ z < v. The pixel sampling period in the output plane P, is b. When ‘=
f 2(d-a> d(D-f) +fa’
the intensity distribution I,(mb,
(38) Z2 in output plane P, has the form
nb) = C C I,]( m-j)a, .i k
(n-k)a]r(jd,
kd) ,
(39)
where m, n, j, k are integers, I, is the input image light intensity distribution in input plane Pi, and rf’jd, kd) the binary distribution in the mask M. That is to say, in the output plane P, we can get the convolution result. 6.2. Liquid crystal light valve - a smart optical nonlinear device ONODIP requires a smart optical nonlinear device which has electrically-controllable nonlinear characteristics. The liquid crystal light valve (LCLV) is one of such optical devices. Because of its nonlinearities of input-output
S. Yuan et al. /Optics Communications 123 (I 9%) 705-715
I output
1 output
0
0
I
713
1input (b)
Fig. 5. The LCLV input-output characteristic curves.
x
OP
I
Collimated Beam ~ Fig. 6. Optical circuit for performing the intersection operation.
P I
CB
Fig. 7. Experimental setup for the implementation of the basic operation.
characteristic, LCLV has been used to perform optical parallel logic gates [ 11,121. The LCLV used in our experiments is a 4.5”twisted nematic reflective LCLV. The device operates between two cross polarizers at room temperature, and it has the capacity of low light-intensity response. The LCLV has the spatial resolving power of 40 lp/mm, the contract ratio of 100 : 1. The working voltage is 4-16 V rms with 200 Hz-5.5 kHz audio frequency. Fig. 5a gives the LCLV input-output characteristic curve when the working voltage is 5.6 V rms with 5 10 Hz, and Fig. 5b gives another input-output characteristic curve when the working voltage is 6.5 V rms with 510 Hz. For clarity, the physical units are omitted in Fig. 5. The first- and second-type nonlinear functions can be respectively simulated with the input-output characteristic curves shown in Figs. 5a and b. The types of nonlinear function can be selected just by changing the size of the working voltage. 6.3. Optical implementation of parallel intersection The image intersection operation is a parallel logic operation which can be performed very easily by using the optical method. That is also an important reason why we select intersection to construct the basic operation of NOBIA. Using a collimated beam to illustrate an image mask X and then another image mask Y, the intersection result can be obtained in the output plane OP, as shown in Fig. 6.
714
S. Yuan et al. /Optics
Communications
123 (1996) 705-715
Fig. 8. Experimental results of the implementation of the basic operations. (a) Image X, (b) structuring Experimental result of F, (X, R, NP’, Y): (e) Experimental result of F,{X, R,NP’, I’).
6.4. Experimental
demonstration
element R, and (c) Image Y. (d)
of the basic operation F,,( . }
The experimental arrangement for the optical implementation of the basic operation F,( . ) is shown in Fig. 7. The diffuse screen G is put just before the input plane IP. The incoherent uniform light (S) is used to illuminate the diffuse screen G. The mask M, of the structuring element R is put close in front of the lens L with a focal length f= 184 mm, and the distance D between the plane IP and lens L is 400 mm. The distance between lens L and the writing plane of LCLV is f+z= 320 mm. Plane IP, mask Mi, lens L and LCLV construct an incoherent optical convoluter. When the image shown in the input plane IP passes through mask M, and lens L, the convolution resultant image V(X, R) is formed in the writing plane of LCLV. A collimated HeNe laser beam (CB) is used as the reading light to read the written-in image of LCLV, where P is a polarizer and A is an analyzer. Selecting different nonlinear characteristics of LCLV to simulate the corresponding formats of the nonlinear function, the corresponding binary image is read out from LCLV. The output binary image passes through the mask M2 coded by binary image Y and executes the intersection operation with image Y. Then in the output plane OP the operation result of F,,(X, R, NP”, Y) can be obtained. Using different nonlinear characteristic curves of LCLV to perform different type nonlinear filterings, the different formats of the basic operation are implemented experimentally. Fig. 8 gives the experimental results. Fig. 8a is an input binary image X. Fig. 8b shows the structuring element R with 3 X 3 image points, and Fig. 8c gives the image Y. Figs. 8d and e respectively give the experimental results of F, {X, R, NP’, Y) and F2{X, R, IQ’*, Y}.
7. Discussion and conclusion
We have proposed a new neighborhood operation binary image algebra to provide a unified theory of the parallel algorithm for morphological image processing. NOBIA has only one basic operation that is a convolution followed by a nonlinear filtering operation and then an intersection operation. All possible morphological image transformations can be expressed readily by NOBIA. NOBIA greatly simplifies the hardware structure and iteration steps of morphological image processing. Then we have presented a new architecture of the optical neighborhood operation digital image processor. The suggested ONODIP is of a high parallelism, simple structure and potential fast processing ability, and can be used to efficiently execute the algorithm of NOBIA. We have also discussed the implementation of the basic operation of NOBIA. Using an incoherent optical convoluter as a three-dimensional free space interconnection device, a smart LCLV as an optical nonlinear device, and a simple optical circuit for the intersection operation, an optical hardware is constituted and used to experimentally demonstrate the basic operation. The experimental results are given. Because the two kinds of nonlinear functions are both thresholding operations which can be performed directly and easily with a smart spatial light modulator, the optical hardware can be used as the CPU of ONODIP. All morphological image processing operations can be performed by executing the CPU repeatedly. NOBIA can be extended to the case of gray scale images by threshold decomposition [ 131, so that ONODIP can also be suitable for the processing of gray scale images.
S. Yuan et al. /Optics
Communications
123 (1996) 705-715
715
Acknowledgements This work was supported by the National Natural Science Foundation Lixue Chen for his helpful discussions throughout the course.
of China. The authors wish to thank Dr.
References 1 I I J. Sena, Image analysis and mathematical morphology (Academic Press, New York, 1982). 12I P. Maragos. Opt. Eng. 26 ( 1987) 623. 131KS. Huang, B.K. Jenkins and A.A. Sawchuk, Comp. Vision Graph. Image Proc. 45 (1989) 295. [ 4 I KS. Huang, A.A. Sawchuk. B.K. Jenkins, P. Chavel, J.M. Wang, A.G. Weber, C.H. Wang and I. Glaser, 151M. Fukui and K. Kitayama, Appl. Optics 31 (1992) 581. 16 1M. Fukui and K. Kitayama, Appl. Optics 31 ( 1992) 4645. 17 1L. Liu. Z. Zhang, and X. Zhang, J. Opt. Sot. Am. A 11 (1994) 1728. 18I L. Liu, Optics Lett. 14 (1989) 482. 19 I Shifu Yuan, Shijie Zhao. Xueru Zhang and Lixue Chen, Proc. SPIE 1989 ( 1993) 402. 1 101 KS. O’Neill and W.T. Rhodes, Proc. SPIE 638 (1986) 41. I I1 I B.K. Jenkins, A.A. Sawchuk, T.C. Strand and B.H. Soffer. Appl. Optics 23 (1984) 3455. ( I2 I Shifu Yuan. Shijie Zhao, Xueru Zhang, Lixue Chen and Jing Hong, Optik 97 (1994) 149. [ 13I F.Y.C. Shih, O.R. Mitchell, IEEE Trans. Pattern Anal. Mach. Intell. 11 (1989) 31.
Appl. Opt. 32 ( 1993) 166
Optics Communications
123 ( 1996 ) 705-7
15
Full length article
Neighborhood operation binary image algebra for optical morphological image processing Shifu Yuan, Guofan Jin, Minxian Wu, Yingbai Yan Tsinghua University, Department
of Precision Instruments, Beijing 100 084, Chinu
Received 27 March 1995
Abstract
A neighborhood operation binary image algebra (NOBIA) which has only one basic operation is presented and developed. The basic operation of NOBIA is a convolution followed by a nonlinear filtering function and then an intersection operation. The parallel architecture of an optical neighborhood operation digital image processor is designed to efficiently perform the parallel morphological image processing algorithm provided by NOBIA, and the optical implementation hardware of the processor is discussed. Using an incoherent optical convoluter as a three-dimensional free space interconnection device, a smart LCLV as an optical nonlinear device, and a simple optical circuit for intersection operation, an optical hardware is constituted and used to realize the basic operation of NOBIA experimentally. Any image transformation can be performed with this hardware by executing the basic operation repeatedly according to proper iterative programming.
1. Introduction With optical parallel architecture a massive number of optical computing systems have been proposed. The most important and efficient application of optical computing is in the field of image processing, and many structures of parallel image algebra suited for optical implementation have been presented. Mathematical morphology [ 1,2] is a widely-used image algebra in which a broad class of image processing operations are represented by morphological transformations (MT) with the combination of dilations and erosions. On the basis of mathematical morphology, Huang et al. have developed a binary image algebra (BIA) , in which any image transformation can be implemented by the use of appropriate reference images and the three fundamental operations of complement, union and dilation [ 3,4]. Along this direction, an image logic algebra (ILA) has been proposed which has more operations [ 5,6]. Recently, Liu et al. have proposed another parallel binary image algebra, called one-operation image algebra (OIA), with a logic operation followed by a dilation [ 71. This paper presents our attempts at searching a simpler algebraic representation for optical implementation of parallel image processing. According to BIA, three hardwares are needed to respectively perform the three fundamental operations. In OIA, an operation is constructed with a logic operation followed by a dilation. In fact, if we can construct a basic operation that can be programmed to perform at least the three fundamental operations of BIA, another binary image algebra can be developed with this basic operation and thus all of the image transformation can be represented by such an image algebra. The dilation of a binary image X with a reference image R can be performed by the convolution of images X and R followed by a thresholding [ 81. We have defined a neighborhood 0030-4018/96/$12.00
0 1996 Elsevier Science B.V. All rights reserved
SSDlOO30-40 I8( 95 )003 17-7
706
S. Yuan et al. /Optics
Communications
123 (1996) 705-715
operation as a convolution followed by a nonlinear filtering with a point nonlinear function f9]. The combination of a neighborhood operation and an intersection operation can be regarded as a basic operation. The intersection, complement and dilation can be efficiently performed with such a basic operation by selecting different formats of the point nonlinear function. Based on the basic operation, we develop a new binary image algebra that is called neighborhood operation binary image algebra (NOBIA) . All morphological image transformations such as erosion, opening, closing and hit-miss transform can be represented easily with the NOBIA. Thus NOBIA provides a simple parallel algorithm for optical morphological image processing. Also, we construct a simple processor, i.e. the optical neighborhood operation binary image processor (ONODIP) , to perform the algorithm of NOBIA. The suggested ONODIP is in fact a single instruction multidata stream (SIMD) type of parallel processor. The ONODIP has the advantage of the inherent parallelism of optics. Because optical free space interconnection is adopted, the processor has a very flexible interconnection style, thus the structuring elements for morphological image processing can not only be very large, but also be altered easily.
2. Neighborhood
operation binary image algebra (NOBIA)
The binary image algebra (BIA) , extending from mathematical morphology, is the basis of NOBIA. For convenience, we here give a short review of BIA. In BIA, a binary image is regarded as a set of image points (pixels with value 1). The domain of images is defined as the universal image, which is the set W= ((x, y) IxEZ,,~EZ,), where Z,= {0, f 1, . . . . kn), and n is a positive integer. The image space S is the power set (the set of all subsets) of the universal image, i.e., S = P( W). Any image is a subset of the universal image. The three fundamental operations of BIA are defined as follows: (1) Complement of an image X:
rz=1(x, Y> I(.& y) E WA (4 Y) e-v .
(1)
(2) Union of two images X and R: Xu R= ((x, y> I (x, y) EXV (x, y) ERJ .
(2)
(3) Dilation of two images X and R: XOR=
{xi +x2,
y1 fy,)
E
WI (x,, Y,) =f,
0,
(xzt ~2) ER
(X+0) A (R+0) I (X=0)
V (R=0)
7
(3)
.
Remarks: E denotes belong to, A denotes AND, V denotes OR, and (21is the null image having no image point. The erosion of an image X by a structuring element R can be written as XeR=%ri,
(4)
where the reflection image R is defined as i={(-.x,
-y)I(x,y)ER].
(5)
The morphological dilation and erosion are in fact a convolution followed by a point nonlinear filtering in the spatial domain [ 81. We define a neighborhood processing operation NP(X, R) as the convolution of a binary image X and a binary image R followed by a point nonlinear function [ 91, i.e., NP(X, R) =NL[ U(X, R)] ,
(6)
where NL(LI) denotes a point nonlinear filtering function, and u is the gray level of a pixel in the convolution resultant image U(X, R) _The image V(X, R) can be expressed as U(X, R) =X*R,
(7)
S. Yunn et al. /Optics
I
Communications
707
I23 (1996) 705-715
_____ L._--.--~ 1
2
N
u
(a) Fig, I. Nonlinear functions function NL’( I( ).
for the neighborhood
operation.
(a) The first-type nonlinear
function
NL’( u). (b) The second-type
nonlinear
where * denotes convolution. The value of pixel (p, q) in the image X can be expressed as x(p, q) which has two possible states 0 and 1, where p=o, +1 1 . . .. 5 II and q = 0, f 1, . . , 5 n. The value of pixel (j, k) in the structuring element R is expressed as r(j, k) which has two possible states 0 and 1, where j= 0, + 1, . .., +n and k=O, * 1, . ., f R. According to Eq. (7), the value of a pixel (p, q) in the image U(X, R) is expressed as U(P. 4) =
C C-Q-j, I
s-k)r(.L
k)
(8)
T
h
where.r(p, q) =0 when ]p] >n or 141 >n. Dilation and the complement of dilation can be performed directly by NP(X, R) with different formats of the point nonlinear function. Figs. la and b respectively give the point nonlinear functions for the dilation and its complement. The function shown in Fig. la is called the first-type nonlinear function and is expressed as NL’( u), while the function shown in Fig. 1b is called the second-type nonlinear function and is expressed as NL*( U) . Then the neighboring processing function with the nth nonlinear function is represented as NP”(X. R) =NL”[ U(X, R)] =NL”[X*R] where II = written as
1,
,
2. The relations between the neighboring
(9) processing
function and the morphological
dilation can be
NP’(X. R) =NL’[
U(X, R)] =XOR,
( 10)
NP’(X. R) =NL’[
U(X, R)] =XOR.
(11)
When we select the reference image R to be an element image Z= { (0, 0) } that has only one image point at the origin, Eqs. ( 10) and ( 11) respectively give the original image X and its complement image X, i.e., NP’ (X. I) = NL’ [ U(X, I) ] =X ,
(12)
NP’(X, I) =NL’[
(13)
U(X, Z)] =X.
According to Eqs. ( 10) and ( 13), the neighborhood processing function NP(X, R) can perform dilation and complement operations. In BIA, the three operations (complement, union and dilation) are chosen as the three elements of the complete set of fundamental operations. Because union can be performed with complement and intersection, the three fundamental operations of the complete set can also be complement, intersection and dilation. Thus the question is how to construct a basic operation that is capable of performing at least complement, intersection and dilation. The basic operation F,,( . ) of the neighborhood operation binary image algebra is then defined as the neighborhood processing function NP(X, R) followed by an intersection, i.e., F,,{X,R,NF”‘,
Y)=NP”(X,R)nY=NL”[X*R]nY,
(14)
S. Yuan et al. /Optics
Communications
:.
, : Binary
- Gray Image Input -
Binary Image
: ;
Nonlinear
Convolution
i Image
Filtering
x
I
j i
123 (1996) 705-715
Image R
Neighborhood
:
Function NL(u)
j
Processing
Binary Image
I
I Nonlinear
Reference
. Output Intersection
Binary Image Y
Operation
Fig. 2. Structure frame of the basic operation of NOBIA.
where n = 1, 2. X and Y are two input binary images, R is a reference image. The structure of the basic operation is shown in Fig. 2. Selecting different forms of the nonlinear function, the basic operation F,,( . } has different forms, i.e., F, (X, R, NP’, Y) =NP’(X,
R) n Y= (XOR)
17 Y,
(15)
F2 (X, R, NP*, Y) =NP’(X,
R) I’-’Y= (XOR)
n Y.
(16)
Complement,
intersection
and dilation can all be performed by the operation F,( . ) :
Fz(X,Z,NP2,
W)=NP’(X,Z)nW=X,
(17)
F,(X,Z,NP’,
Y)=NP(x,z)nY=xnr,
(18)
F, (X, R, NP’, W) =NP’(X,
R) f7 W=XOR,
(19)
where W is the universal image. Because any image transform can be implemented by the iteration of complement, intersection and dilation with appropriate reference images, it can also be represented by the iteration of the basic operation F,( . ) . On this basis, we define NOBIA as an algebra with an image space S= P( W) and a basic operation F,,( . }, which is represented as NOBIA=
(P(W),
3. NOBIA representation All morphological
(20)
F,,) .
of morphological
image transformations
image transformations can be derived from F,,( . ). Here we give some examples:
(1) Union of two images X and Y: XUY=XnY=F,(F,(X,Z,NP2,
F,(Y,Z,
NP*, W}),Z,NP*,
W) .
(21)
(2) Difference of X by K XIY=XnY=F,(Y,Z, (3)
Erosion of image X with a reference image R: XeR=F,(F,(X,Z,
(4)
(22)
NP2,X).
Symmetric
difference
NP2, W),d,
NP2, W) .
of two images X and Y:
(23)
S. Yuan et al. /Optics Communications 123 (1996) 705-715
709
XAY= (X/Y) u (Y/X) = F, ( Y, I, NP’, X) U F2 {X, I, NP’, Y) = F2 ( F2 ( F2 { Y, Z, NP*, X), I, NP’, Fz ( F2 {X, I, NP’, Y), I, NP*, W} }, I, NP*, W) . (24)
(5)
Opening of image X by a reference image R: XoR==(XeR)OR=F,{F,(F,{X,
NP’, I, W},k,NP2,
W), R,NP’,
W).
(25)
(6) Closing of X by a reference image R: X.R=(XOR) (7)
R, NP’, W), li, NP*, W) .
(26)
transform of image X by an image pair R( Rn Rb) :
Hit-miss XOR=
eR=F2(F2{X,
(XeR,-)
fl (&R,)
=J?@Z?,~IX@Z&,
= F, (F, (X, I, NP’, W}, dr, NP*, F2 (X, d,, NP2, W) } . (8)
Edge detection. The out-object
edge detection can be denoted as
ED,,(X) = (XOR)/X=F,{X,
Z, NP2, F, {X, R, NP1, W}} .
(27)
(28)
The in-object edge detection can be denoted as ED,(X) =X/(XBR)
= (%Ori) nX=F,
{F,(X,
Z, NP2, W}, d, NP’, X) .
(29)
The central edge detection can be denoted as ED,(X) (9)
= (XOR)
nXBR=F,
(F2{X, Z, NP2, W), d, NP’, F, (X, R, NP’, W)} .
Shifting an image X from the coordinate X@((i,j)l=F,(X,
( 10) Noise remoual.
(30)
(x, y) to (x + i, y -kj) :
l(i,j)),NP’,WJ.
(31)
The removal of noise points in an image X can be denoted as
NR(X) =XoR@R =F2(Fz(F,(F2{F2{X,NP2,Z,
W},ri,NP2,
W), R,NP’,
W}, R, NP’, W},i,NP’,
W}
(32)
= F, [ F2 { F2 ( F2 { F2 (X, R, NP2, W}, i, NP*, W), NP’, Z, WI, I?, NP’, W), R, NP’, W}
(33)
or NR(X)X.RoR
where R is a small structuring element. ( 1 I ) Skeltonization. In mathematical morphology, SK(X)
=
the skeleton of an object X is expressed as
U S,(X), oCn
S,(X) = (XenR)
n (XenR)oR,
wheren=1,2,...,N,andN=MAX{kJX8kR#0}.InNOBIA,itcanbedenotedas
710
S. Yuan et al. /Optics
SK(X) = f-l
f-l
J,(X) = Fz (
O<“=GN
J,(x)
J,(X),
Z, NP’, W) ,
1, IQ”
, Jz(x)
I23 (1996) 705-7I5
(35)
O
= FI
J,(X) =S,(X)
f-l
Communications
{F,
1..
J’,
= ((X8nR)
{.f,
(x>,
1,
. . ., 1,
NP’, JN- ,(x) }, 1, NP’, JN(x)
n (X8nR) eR0R=X0nzh--lX0(n+
wheren=l,2
4. Comparison
,...,
Z, NP’, W), (n+l)ii,
(364
l)I;i0R
= F2 { F2 ( F2 {X, Z, NP*, W}, ni, NP*, F2(F2[F2{X,
} ,
(36b) NP’, W), R, NP*, W)}, Z, NP’, W),
N,andN=MAX{k~F2(F2(X,Z,NP2,W),kZ?,NP2,W}#0}.
of NOBIA with BIA and OIA
In comparison with the existing image algebra, NOBIA has an obvious advantage in the parallel algorithm of image processing. In NOBIA, there is only one basic operation which has dual processing modes by selecting two nonlinear functions. Then only one unit of optical hardware is required by NOBIA, and the control and programming of the processing system can be very simple. But in BIA there must be three hardware units to respectively perform three fundamental operations. An efficient morphological image processing structure should have less basic operations and less iteration steps for a morphological transform. NOBIA has the smallest number of basic operations (only one) ; furthermore, the suggested NOBIA may also simplify greatly the iteration steps of morphological image processing. For example, a closing operation requires four iteration steps in BIA, but it needs only two steps in NOBIA; the hit-miss transform needs six iteration steps in BIA, while in NOBIA it needs only three steps. We also compare the iteration steps of other morphological transforms in BIA with that in NOBIA, and this gives the comparing results in Table 1. NOBIA also has obvious advantages comparing with OIA, though there is also only one operation in OIA. The basic operation of OIA has at least 16 processing modes (for 16 binary logic functions) which must be programmed to select. Moreover, the complete system of OIA has two hardware units because the basic operation is divided into two parts (one for logic operations, the other for dilation) by the image output port. Then the control and programming will be very complicated. But in NOBIA there is truly one hardware unit, and it only requires to programmably select the formats of a nonlinear function. Table 1 The comparison
results of iteration steps in NOBIA and BIA Image transform
Iteration steps in BIA
Iteration steps in NOBIA
complement union intersection dilation erosion difference symmetric difference opening closing hit-miss transform central edge detection
1
1
I
3
4 I 3 3 7 4 4 6 7
I 1 2
I 5 3 2 3 3
S. Yuan et al. /Oprics Communications 123 (1996) 705-715
711
Control Unit Image Data Flow CxSignal
Clock
Flow
Program Counter
F,(‘} Basic Operation
c
Test and Feedback Controller
:
Image Output
DeVice
Parallel Processing Unit Fig. 3. Architecture of the ONODIP 5. Optical neighborhood
operation digital image processor (ONODIP): Architecture
NOBIA provides a parallel algorithm for morphological image processing. The architecture for implementing the algorithm should be parallel. Fig. 3 shows the architecture of the optical neighborhood digital image processor (ONODIP) for executing the parallel algorithms of NOBIA. ONODIP is in fact a cellular SIMD machine and consists of a parallel processing unit under the supervision of a control unit. The control unit includes a clock, a program counter, a test and feedback controller, and an instruction memory. The parallel processing unit includes an image memory device for storing multiple images, three frame buffers, and a central processing unit (CPU) for executing the basic operation F,{ . ). The three buffers X, R and Y are respectively used to store the three input images for the central processing unit. The images stored in the image memory device can be respectively sent to the three frame buffers. The CPU can be controlled by the controlling unit to perform proper functions. The ONODIP operates as follows: ( 1) the corresponding images are sent to the frame buffers X, R and Y from the image memory device; (2) the basic operation F,{ ) is executed by selecting the proper formats of the nonlinear function; (3) the result is output, tested for program control, or sent back to the image memory device. Then the system can be fed back to step ( 1) .
6. Optical implementation
hardware of the ONODIP
The optical implementation hardware of the ONODIP is the CPU for optically performing the basic operation F,,( ) . According to Eq. ( 14), we need an optical convoluter, an optical spatial nonlinear device with two operation modes, and an optical parallel intersection processor.
712
S. Yuan et al. /Optics
Communications 123 (1996) 705-715
ML
Pl
p2
P”
Fig. 4. An incoherent optical convoluter.
6.1. Incoherent optical convoluter The convolution can be performed with the incoherent optical convoluter [9] as shown in Fig. 4. It is a generalization of the incoherent optical correlator which has been used to implement morphological image processing
1101. The input plane P, is illuminated by a uniform incoherent light. A mask M is put close in front of an imaging lens L with a focal lengthf. The distance between plane P, and mask M is D. The distance between Mask M and lens L is very small and can be neglected. The mask M is used to display the structuring element image R. The shortest distance between two near pixels (called pixel sampling period) in plane P, is a. The pixel sampling period in mask M is d. Plane P, is the imaging plane in which a clear image of the input plane Pi can be obtained. The distance v between lens L and plane P y can be written as
Df
(37)
ll=o-f’
The distance between lens L and output plane P, isf+ z < v. The pixel sampling period in the output plane P, is b. When ‘=
f 2(d-a> d(D-f) +fa’
the intensity distribution I,(mb,
(38) Z2 in output plane P, has the form
nb) = C C I,]( m-j)a, .i k
(n-k)a]r(jd,
kd) ,
(39)
where m, n, j, k are integers, I, is the input image light intensity distribution in input plane Pi, and rf’jd, kd) the binary distribution in the mask M. That is to say, in the output plane P, we can get the convolution result. 6.2. Liquid crystal light valve - a smart optical nonlinear device ONODIP requires a smart optical nonlinear device which has electrically-controllable nonlinear characteristics. The liquid crystal light valve (LCLV) is one of such optical devices. Because of its nonlinearities of input-output
S. Yuan et al. /Optics Communications 123 (I 9%) 705-715
I output
1 output
0
0
I
713
1input (b)
Fig. 5. The LCLV input-output characteristic curves.
x
OP
I
Collimated Beam ~ Fig. 6. Optical circuit for performing the intersection operation.
P I
CB
Fig. 7. Experimental setup for the implementation of the basic operation.
characteristic, LCLV has been used to perform optical parallel logic gates [ 11,121. The LCLV used in our experiments is a 4.5”twisted nematic reflective LCLV. The device operates between two cross polarizers at room temperature, and it has the capacity of low light-intensity response. The LCLV has the spatial resolving power of 40 lp/mm, the contract ratio of 100 : 1. The working voltage is 4-16 V rms with 200 Hz-5.5 kHz audio frequency. Fig. 5a gives the LCLV input-output characteristic curve when the working voltage is 5.6 V rms with 5 10 Hz, and Fig. 5b gives another input-output characteristic curve when the working voltage is 6.5 V rms with 510 Hz. For clarity, the physical units are omitted in Fig. 5. The first- and second-type nonlinear functions can be respectively simulated with the input-output characteristic curves shown in Figs. 5a and b. The types of nonlinear function can be selected just by changing the size of the working voltage. 6.3. Optical implementation of parallel intersection The image intersection operation is a parallel logic operation which can be performed very easily by using the optical method. That is also an important reason why we select intersection to construct the basic operation of NOBIA. Using a collimated beam to illustrate an image mask X and then another image mask Y, the intersection result can be obtained in the output plane OP, as shown in Fig. 6.
714
S. Yuan et al. /Optics
Communications
123 (1996) 705-715
Fig. 8. Experimental results of the implementation of the basic operations. (a) Image X, (b) structuring Experimental result of F, (X, R, NP’, Y): (e) Experimental result of F,{X, R,NP’, I’).
6.4. Experimental
demonstration
element R, and (c) Image Y. (d)
of the basic operation F,,( . }
The experimental arrangement for the optical implementation of the basic operation F,( . ) is shown in Fig. 7. The diffuse screen G is put just before the input plane IP. The incoherent uniform light (S) is used to illuminate the diffuse screen G. The mask M, of the structuring element R is put close in front of the lens L with a focal length f= 184 mm, and the distance D between the plane IP and lens L is 400 mm. The distance between lens L and the writing plane of LCLV is f+z= 320 mm. Plane IP, mask Mi, lens L and LCLV construct an incoherent optical convoluter. When the image shown in the input plane IP passes through mask M, and lens L, the convolution resultant image V(X, R) is formed in the writing plane of LCLV. A collimated HeNe laser beam (CB) is used as the reading light to read the written-in image of LCLV, where P is a polarizer and A is an analyzer. Selecting different nonlinear characteristics of LCLV to simulate the corresponding formats of the nonlinear function, the corresponding binary image is read out from LCLV. The output binary image passes through the mask M2 coded by binary image Y and executes the intersection operation with image Y. Then in the output plane OP the operation result of F,,(X, R, NP”, Y) can be obtained. Using different nonlinear characteristic curves of LCLV to perform different type nonlinear filterings, the different formats of the basic operation are implemented experimentally. Fig. 8 gives the experimental results. Fig. 8a is an input binary image X. Fig. 8b shows the structuring element R with 3 X 3 image points, and Fig. 8c gives the image Y. Figs. 8d and e respectively give the experimental results of F, {X, R, NP’, Y) and F2{X, R, IQ’*, Y}.
7. Discussion and conclusion
We have proposed a new neighborhood operation binary image algebra to provide a unified theory of the parallel algorithm for morphological image processing. NOBIA has only one basic operation that is a convolution followed by a nonlinear filtering operation and then an intersection operation. All possible morphological image transformations can be expressed readily by NOBIA. NOBIA greatly simplifies the hardware structure and iteration steps of morphological image processing. Then we have presented a new architecture of the optical neighborhood operation digital image processor. The suggested ONODIP is of a high parallelism, simple structure and potential fast processing ability, and can be used to efficiently execute the algorithm of NOBIA. We have also discussed the implementation of the basic operation of NOBIA. Using an incoherent optical convoluter as a three-dimensional free space interconnection device, a smart LCLV as an optical nonlinear device, and a simple optical circuit for the intersection operation, an optical hardware is constituted and used to experimentally demonstrate the basic operation. The experimental results are given. Because the two kinds of nonlinear functions are both thresholding operations which can be performed directly and easily with a smart spatial light modulator, the optical hardware can be used as the CPU of ONODIP. All morphological image processing operations can be performed by executing the CPU repeatedly. NOBIA can be extended to the case of gray scale images by threshold decomposition [ 131, so that ONODIP can also be suitable for the processing of gray scale images.
S. Yuan et al. /Optics
Communications
123 (1996) 705-715
715
Acknowledgements This work was supported by the National Natural Science Foundation Lixue Chen for his helpful discussions throughout the course.
of China. The authors wish to thank Dr.
References 1 I I J. Sena, Image analysis and mathematical morphology (Academic Press, New York, 1982). 12I P. Maragos. Opt. Eng. 26 ( 1987) 623. 131KS. Huang, B.K. Jenkins and A.A. Sawchuk, Comp. Vision Graph. Image Proc. 45 (1989) 295. [ 4 I KS. Huang, A.A. Sawchuk. B.K. Jenkins, P. Chavel, J.M. Wang, A.G. Weber, C.H. Wang and I. Glaser, 151M. Fukui and K. Kitayama, Appl. Optics 31 (1992) 581. 16 1M. Fukui and K. Kitayama, Appl. Optics 31 ( 1992) 4645. 17 1L. Liu. Z. Zhang, and X. Zhang, J. Opt. Sot. Am. A 11 (1994) 1728. 18I L. Liu, Optics Lett. 14 (1989) 482. 19 I Shifu Yuan, Shijie Zhao. Xueru Zhang and Lixue Chen, Proc. SPIE 1989 ( 1993) 402. 1 101 KS. O’Neill and W.T. Rhodes, Proc. SPIE 638 (1986) 41. I I1 I B.K. Jenkins, A.A. Sawchuk, T.C. Strand and B.H. Soffer. Appl. Optics 23 (1984) 3455. ( I2 I Shifu Yuan. Shijie Zhao, Xueru Zhang, Lixue Chen and Jing Hong, Optik 97 (1994) 149. [ 13I F.Y.C. Shih, O.R. Mitchell, IEEE Trans. Pattern Anal. Mach. Intell. 11 (1989) 31.
Appl. Opt. 32 ( 1993) 166