- Email: [email protected]
Optical multiresolution analysis with spatial localization
Optics Communications 283 (2010) 2056–2060
Contents lists available at ScienceDirect
Optics Communications journal homepage: www.elsevier.com/locate/optcom
Optical multiresolution analysis with spatial localization Javier Mazzaferri 1, Silvia Ledesma * Image Processing Laboratory (LPI), Physics Department, ‘J.J. Giambiagi’, FCEyN, University of Buenos Aires, Argentina
a r t i c l e
i n f o
Article history: Received 29 September 2009 Received in revised form 7 January 2010 Accepted 10 January 2010
Keywords: Image processing Fourier optics Multiresolution analysis Feature extraction
a b s t r a c t Multiresolution analysis is very useful for characterization of textures, segmentation tasks, and feature enhancement. The development of optical methods to perform such procedures is highly promissory for real-time applications. Usually, the optical implementations of multiresolution analysis consist in the decomposition of the input scene in different frequency bands, obtaining various filtered versions of the scene. However, under certain circumstances it could be useful to provide just one version of the scene where the different filters are applied in different regions. This procedure could be specially interesting for biological and medical applications in situations when the approximate localization of the scale information is known a priori. In this paper we present a fully optical method to perform multiresolution analysis with spatial localization. By means of the proposed technique, the multi-scale analysis is performed at once in a unique image. The experimental set-up consist of a double-pass convergent optical processor. The first stage of the device allows the multiple band decomposition, while the second stage confines the information of each band to different regions of the object and recombines it to achieve the desired operation. Numerical simulations and experimental results, which prove the very good performance of the method, are presented. Ó 2010 Elsevier B.V. All rights reserved.
1. Introduction The interest in studying and developing optical architectures to perform real-time image processing, has notoriously grown over recent years; relevant examples of these applications can be found in recent publications [1–3]. For instance, these kind of devices have been implemented in combination with optical microscopes to carry out specific operations, by only optical methods [4,5]. Likewise, the optical processing techniques could improve the in vivo measurements, present in medical applications. In particular, feature extraction of medical images could help for the diagnosis of some diseases [6]. Moreover, in biological images it could be of interest to follow in real-time the evolution of a specific feature. The wavelet transforms have been extensively employed for multiresolution analysis, showing excellent capabilities. Furthermore, these techniques have been widely used for different image processing applications, such as: texture analysis, feature extraction and image segmentation among others [7–10]. This method has been employed in diverse technological and scientific applications related to biological [11–13] and medical [14,6] image treatment.
* Corresponding author. E-mail address: [email protected] (S. Ledesma). 1 Present address: Instituto de Ciencia de Materiales de Aragón (ICMA), Consejo Superior de Investigaciones Científicas-Universidad de Zaragoza, Facultad de Ciencias, Plaza San Francisco s/n, 50009 Zaragoza, Spain. 0030-4018/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.optcom.2010.01.018
It is a common practice, in multiresolution analysis, to apply a bank of spatial filters and obtain a set of images at different spatial resolutions. In these cases the whole function is filtered to each scale. In fact, when wavelet transforms are employed to perform this processing, it is usual to present the results according to typical arrangements of the different filtered images [15]. The information is displayed separately in different images, each one corresponding to an individual scale parameter. However, for certain applications it could be of interest to analyze one region of the image in one frequency band, and other regions of the image in other frequency bands, displaying the result in a unique final image. The implementation of such a procedure could be simple to perform digitally, however in this case, it would not be done in real-time. Several experimental architectures have been proposed to implement the multiresolution analysis optically. Clear examples of this are two pioneering works presented by Mendlovic and Ouzieli [16,17], where the Dammann grating and computer-generated multireference matched filters were employed to achieve the direct and inverse 2D wavelet transform optically. In all these cases, the information corresponding to each scale is obtained in separated images. As it was mentioned above, until now the procedures needed to obtain a unique image with different scales in different regions of the object should be achieved with digital operations. The novelty of the method presented here is that it performs this procedure in real-time with a fully optical process. We propose a device based on the wavelet transform that can achieve optical
2057
J. Mazzaferri, S. Ledesma / Optics Communications 283 (2010) 2056–2060
spatially-localized multiresolution analysis. The design consists in a double-pass optical correlator, in which the frequency decomposition is done with a bank of tuning filters which are binarized versions of the mexican hat wavelet. The paper is organized as follows: in Section 2 the proposed method is presented as well as the theoretical fundamentals of the technique are exposed. The experimental details are explained and both the numerical simulation and experimental results are presented in Section 3. In Section 4 the conclusions are summarized. 2. Localized multiresolution analysis The optical multiresolution analysis method presented in this work is based on a convergent multichannel processor [2]. The design proposed here, consists in a double-pass correlator which is sketched in Fig. 1. The monochromatic light source S and lens L1 shape a convergent wave-front that impinges the object f. The amplitude of the optic field in the plane conjugated with S through L1 is a scaled version of the Fourier transform F of f. The phase of the field in this plane is quadratic and represents a defocus that can be compensated by properly placing the lens L2 and the set of binary masks BM to achieve a focused result. The filter H, designed to match the scale of F, is placed in the Fourier plane mentioned above. Lens L2 images the filtered version I just in front of BM, completing the first stage of the processor. This step allows to separate the information of the object in individual frequency bands. A set of replicas of the input object, each one filtered at a particular scale, is obtained in different regions of I. By means of BM, a different part of each replica is conserved while the rest is blocked. The light direction is reversed with mirror M1 and the second stage of the device, composed by lens L2, filter H, beam splitter BS and mirror M2, recombines the masked components of I to compose the final image s. It should be pointed out that the masking operation introduces some distortions in the final image. However we will demonstrate with some approximations, that in typical working conditions these effects may be neglected. The multiple filter H, which is analogous to the filter proposed in Ref. [2], is composed by N different filters Hj(u, v), where u and v are the spatial frequencies corresponding to the spatial Cartesian coordinates x and y of the object [2]. Each filter Hj is multiplied by a linear phase factor as is explicitly shown below
Hðu; v Þ ¼
N X
i/jL
Hj ðu; v Þe ;
ð1Þ
j¼1
where f/jL g is a set of linear phases defined as /jL ¼ 2pðukj þ v nj Þ with different values of the pair (kj,nj) for each filter Hj. If an object f is processed with such a filter, the image in the correlation plane (x2,y2) can be described as
Iðx2 ; y2 Þ ¼ F1 ½Fðu; v Þ Hðu; v Þ ¼
N X
f hj dðx2 kj ; y2 nj Þ;
where F1 ½ stands for the inverse Fourier transform, and denotes the convolution operation. In each term of the sum in Eq. (2) there is the convolution between the input scene and the impulse response hj of each filter Hj. The convolutions with the Dirac’s delta functions d(x, y) indicate that each filtered image is displaced off-axis in the intermediate plane. In this work, we have used a set of filters {Hj} based on the Mexican hat wavelet
Wðu; v Þ ¼ 4p2 ðu2 þ v 2 Þ exp 2pðu2 þ v 2 Þ ;
ð3Þ
expressed in the Fourier domain, which is frequently employed for feature extraction and multiresolution analysis. The filters of the bank are based on different scaled versions of Eq. (3). As we are mainly interested in the pass-band feature of these functions, we have used the binarized versions WB(u, v) described in Eq. (4). To physically represent H we choose the phase modulation mode of a SLM, and to this end the amplitude of the filter was phase-encoded [2]. The binary filters are described by
( W B ðu; v Þ ¼
1 ! 0
C D2
2
2 < u2 þ v 2 < C þ D2
! elsewhere
;
ð4Þ
where D and C are the width and centre of the annular-shaped frequency band. As we have mentioned above, we employed the bank of filters {Hj} described by Eq. (4), where the frequency bands of each filter has no intersection with the frequency bands of the other filters of the bank. As it was stated in Eq. (2), the function I(x2, y2) has the distribution sketched in Fig. 2a. A mirror M1, placed in the intermediate plane (x2, y2), conjugates the wave-front and inverts the direction of propagation. This plane acts as the input of the second stage of the device. The wavefront emerging from M1 goes through L2, and its Fourier transform is obtained before the filter. For the second stage, the filter transmission can be described by H(u, v). The optical field emerging from the filter, can be expressed as:
F½I Hðu; v Þ ¼ F ðu; v ÞH ðu; v Þ Hðu; v Þ ¼ F ðu; v ÞjHðu; v Þj2 :
ð5Þ
As H(u, v) is a phase only function, it has constant unitary amplitude, therefore the field distribution of Eq. (5) can be simplified to F*(u, v). By means of the beam splitter BS and a second mirror
a Fig. 1. Proposed optical device for spatially-localized multiresolution analysis.
ð2Þ
j¼1
b
Fig. 2. (a) Sketch of the intermediate image obtained with the first stage of the double-pass correlator, (b) arrangement of binary masks in the intermediate plane.
2058
J. Mazzaferri, S. Ledesma / Optics Communications 283 (2010) 2056–2060
M2, the inverse Fourier transform of Eq. (5) is obtained at the end of the second stage of the device. Therefore, if no binary mask is placed in M1, the intensity at the output plane results:
jsðx3 ; y3 Þj2 ¼ jF1 ½F ðu; v Þj2 ¼ jf ðx1 ; y1 Þj2 :
ð6Þ
It should be mentioned that the coordinates (x1,y1) and (x3, y3) are linearly related by a constant factor, which characterizes the magnification of the system. In what follows, and for the sake of clarity, we will employ the coordinates (x, y) to describe the spatial domain. Eq. (6) shows that the second stage counteracts the effect of the first stage. Let us now consider the effect of positioning a binary mask BM, composed by the elements {mj}, in the intermediate plane (x2, y2). The stops are placed in such a way that each convolution f hj in I is masked by a binary element mj which, in turn, does not overlap with the others. The spatial arrangement of the masks is sketched in Fig. 2b. If BM is placed, the field distribution after the mask, in the second stage, can be described by N nh X
i o f hj mj dðx2 kj ; y2 nj Þ;
ð7Þ
j¼1
N X
hj ðx; yÞ ðmj ðx; yÞf ðx; yÞÞ;
ð13Þ
j¼1
where we have used Eq. (4) to obtain that hj hj = hj. Employing the same approximation for hj ðx; yÞ, it is possible to write the field distribution in the zero order of the final plane as N X
mj hj f :
ð14Þ
j¼1
Finally, as the set of functions {mj} are binary and do not overlap, the intensity of the optical field in the final plane can be described as
jsj2 ¼
N X
2 mj f hj :
ð15Þ
j¼1
where it was used that both d and mj are real functions. The optical field in the Fourier plane, is described as N n h io X M j ðu; v Þ F ðu; v Þ Hj ðu; v Þ
ð8Þ
j¼1
e
support region of mj, and mj does not have variations at the scale of hj. Therefore the integrand of Eq. (12) is not null only if r x and s y. Within this approximation, we can replace mj(x, y) by mj(r, s) and obtain
j
i/L ðu;v Þ
:
3. Results
In Eq. (8), Mj denotes the Fourier transform of mj. The optical field after passing through the filter can be written as N X N n h io X M j ðu; v Þ F ðu; v Þ Hj ðu; v Þ j¼1
ð9Þ
l¼1
Hl ðu; v Þe
ið/jL ðu;v Þþ/lL ðu;v ÞÞ
:
By employing the definition of /L, the exponential factor becomes 1 in those terms where j = l. Therefore, the double summation in Eq. (9) can be split in two summations. One of them corresponding to j = l and the other to j – l. The terms with j = l contribute to the zero order in the final plane, meanwhile the terms with j – l yield offaxis field distributions due to the linear phase factors in the exponentials. Besides, as the filters Hj ðu; v Þ and Hl(u,v) do not intersect, the terms with j – l vanish to zero. Then, the terms that contribute to the final plane are only in the zero order region and can be described as N n h io X M j ðu; v Þ F ðu; v Þ Hj ðu; v Þ
ð10Þ
j¼1
Hj ðu; v Þ: By means of the inverse Fourier transform, the complex amplitude of the field in the final plane can be calculated as: N X
Each term in the sum of Eq. (15) is the product between a binary mask mj and a different squared convolution jf hjj2. Therefore, the whole sum represent different regions of the input object, where each one is filtered to a different scale. This is the desired result.
n h io hj mj f hj :
We have tested the performance of the proposed method with an input scene consisting in a variation of the first three steps of generation of the Sierpinsky carpet, as is shown in Fig. 3a. This scene has the most important features at spatial frequen1 1 and 36 . To analyze this function, cies (in units of 1/pix) around 12 we have employed a filter bank {H1(u, v),H2(u, v)} composed by two binary functions as described by Eq. (4). H1 filters a frequency 1 1 , and H2 filters a frequency band centred at 36 . The band centred at 12 binary filters were obtained by choosing a threshold value that produces frequency bands as wide as possible without overlap. The resulting D values, in units of 1/pix, for H1 and H2 are approx1 1 and 36 , respectively. With these functions, we have deimately 12 signed a multiple phase filter as described in Eq. (1), which was physically implemented by encoding the amplitude information in the phase [2,18]. The linear phases employed to build the filter H were selected to redirect the information filtered by H1 and H2 off-axis. The light corresponding to H1 was redirected in the horizontal direction meanwhile the light filtered with H2 was deviated in the vertical direction. The information corresponding to spatial frequencies outside the bands of H1 and H2 is obtained in the zero order of the intermediate plane. The resulting filter H is shown in Fig. 3b. By means of numerical simulations we have computed the output field distribution of the first stage of the device, which is shown in Fig. 4a. In the figure, three orders can be distinguished.
ð11Þ
j¼1
By explicitly writing the convolution f hj , and introducing mj(x, y) in the integral, it is obtained that N Z Z X
1
1
j¼1
mj ðx; yÞf ðr; sÞhj ðx r; y sÞdrds
ð12Þ
hj ðx; yÞ: As it was previously pointed out, in most of the cases of interest, Eq. (12) can be approximated by a simpler expression. Indeed, generally, the support region of hj is very much smaller than the
Fig. 3. (a) Input scene, (b) multiple phase filter (gray levels from black to white denote phase values ranging from 0 to 2p).
J. Mazzaferri, S. Ledesma / Optics Communications 283 (2010) 2056–2060
2059
m2
m0
m1
a
b
Fig. 4. (a) Field distribution in the intermediate plane of the device, (b) binary mask.
The order on the left is the input scene filtered with H1, the upper order is the same scene filtered with H2, and in the zero order there is the information of the input scene that does not belong to the frequency bands of the filters. To test our method, we have designed a binary mask that blocks the different regions of the intermediate field of Fig. 4(a). The mask is designed to filter in different frequency bands the information of two different parts of the objects. In the smaller upper part of the scene, no filter is applied. Besides, in the bigger lower part of the object, H2 is used. The designed mask is shown in Fig. 4b. As it can be observed, the mask blocks the lower region of the left and zero orders. Therefore, the only information in this part of the object that reaches the final plane, belongs to the upper order. Within the upper part of the scene, no order is blocked, therefore the object is fully reconstructed there. Fig. 6a shows the numerical simulation of the final plane where the desired results can be observed. In order to test the experimental performance of the method, the optical device of Fig. 5 was mounted. As can be noted, the system follows the sketch presented in Fig. 1. The light of an Argon laser is filtered and expanded to act as the point light source of the device. The input scene f(x, y) and the filter H(u, v) are displayed in amplitude (SLM1) and phase (SLM2) spatial modulators respectively. Lens L1 is positioned to obtain the Fourier transform of f(x,y) in the plane where H(u, v) is placed. The SLM’s consist of LCD panels and a set of wave quarter plates (WQP) and linear polarizers (P), that manages the incident and emerging elliptical polarization, as is sketched in the right of Fig. 5. By means of controlling the incoming and emerging polarization states, the desired modulation is obtained [19,20]. The mirror M1 on the right of the system inverts the direction of propagation of the light, and the binary mask BM is positioned in the image plane of the input object, in contact with M1. The second stage of the double correlator is built with the beam splitter BS and
L3 M2
CCD SLM
SLM1
BS
SLM2
WQP1 WQP2
BM
Ar Laser
M1 L1
L2 Video projector
P1
PC
Fig. 5. Double-pass experimental device.
LCD
P2
Fig. 6. (a) Numerical simulation result obtained with the binary mask of Fig. 4b. Experimental results obtained: (b) without binary mask, (c) blocking the 0 and 2 orders, (d) blocking the 0 and 1 orders, (e) with the binary mask of Fig. 4b.
the mirror M2. Lens L3 is introduced to make the output image fit the size of the CCD sensor. It should be noted that to have the filter positioned in the Fourier planes of both the first and the second stage of the device, lens L2 was placed to have its focal plane in the Fourier plane of the first stage. With the input scene and the Fourier filter shown in Fig. 3, experimental results were obtained. As a first test, the output image without any binary mask in the intermediate plane was obtained. As it was predicted in Eq. (6), the image of the input scene was obtained in the final plane. Fig. 6b shows the experimental image. By means of a binary mask, the upper and zero order in the intermediate plane were blocked. In this case, the output of the system showed the input scene filtered with H1, as can be observed in Fig. 6c. Similarly, with another mask, the left and zero orders of the intermediate plane were blocked obtaining, in this case, the input scene filtered with H2. The corresponding output of the system is shown in Fig. 6d. Finally, to perform the spatially-localized multiresolution analysis, the mask presented in Fig. 4b was placed in the intermediate plane of the system. Fig. 6e shows the experimental result of the localized multiresolution analysis, which coincides very well with the numerical simulations presented in Fig. 6a. Despite the optical implementation could be further optimized, excellent experimental results were obtained. It was demonstrated that the presented system permits the full-optical realization of spatially localized multiresolution analysis, which has promissory applications. 4. Conclusions In this work, we propose an optical arrangement to perform spatially-localized multiresolution analysis. The proposed device allows to process different subregions of the input scene at different scales, avoiding to split the scale information in separated images. The whole process is performed optically, and therefore, in real-time. The fundamental interest of this kind of devices is the fact that they can be integrated in other optical systems, as for example optical microscopes, allowing the realization of fulloptical real-time image processing. The proposed method is based in a double-pass multichannel correlator that is based on wavelet multiresolution analysis. The first stage of the device performs the sub-band decomposition and, by means of a set of binary masks, the spatial localization of scale information is obtained at the end of the second stage. Numerical simulations and experimental
2060
J. Mazzaferri, S. Ledesma / Optics Communications 283 (2010) 2056–2060
tests were provided to demonstrate the high performance of the method. Very good results were obtained, suggesting promissory applications of the proposed architecture. Acknowledgement This work was supported with grants from Consejo Nacional de Investigaciones Cientı´ficas y Técnicas (CONICET), Universidad de Buenos Aires (UBA) and Agencia Nacional de Promoción Cientı´fica y Tecnológica (ANPCyT) from Argentina. S. Ledesma is a member of CONICET. References [1] B. Javidi, I. Moon, S. Yeom, E. Carapezza, Optics Express 13 (12) (2005) 4492. [2] J.E. Mazzaferri, S.A. Ledesma, C.C. Iemmi, Journal of Optics A: Pure and Applied Optics 5 (2003) 425. [3] H. Zhang, C.M. Cartwright, M.S. Ding, W. Gillespie, Optics Communications 185 (2000) 277. [4] S. Fürhapter, A. Jesacher, S. Bernet, M. Ritsch-Marte, Optics Express 13 (2005) 689.
[5] A. Jesacher, S. Fürhapter, S. Bernet, M. Ritsch-Marte, Physical Review Letters 94 (2005) 233902. [6] W. Qian, L. Li, L.P. Clarke, Medical Physics 26 (3) (1999) 402. [7] C.S. Lu, P.C. Chung, C.E. Chen, Pattern Recognition 30 (5) (1997) 729. [8] N. Fatemi-Ghomi, P. Palmer, M. Petrou, in: Proceedings of the Workshop on Performance Characteristics of Vision Algorithms (ECCV 96), 1996. [9] M. Unser, IEEE Transactions on Image Processing 4 (11) (1995) 1549. [10] T. Chang, C. Kuo, Image Processing IEEE Transactions on 2 (4) (1993) 429. [11] G. Papari, P. Campisi, N. Petkov, A. Neri, EURASIP Journal on Applied Signal Processing 2007 (2007) 119. [12] J.-C. Olivo-Marin, Pattern Recognition 35 (9) (2002) 1989. [13] M. Unser, A. Aldroubi, Proceedings of the IEEE 84 (4) (1996) 626. [14] C. Wählby, I. Sintorn, F. Erlandsson, G. Borgefors, E. Bengtsson, Journal of Microscopy 215 (2004) 6776. [15] S.G. Mallat, S. Zhong, IEEE Transactions on Pattern Analysis and Machine Intelligence 14 (1992) 710. [16] D. Mendlovic, I. Ouzieli, I. Kiryuschev, E. Marom, Applied Optics 34 (1995) 8213. [17] I. Ouzieli, D. Mendlovic, Applied Optics 35 (1996) 5839. [18] I. Moreno, E. Ahouzi, J. Campos, M.J. Yzuel, Applied Optics 36 (1997) 7428. [19] A. Márquez, J. Campos, M.J. Yzuel, I. Moreno, J.A. Davis, C. Iemmi, A. Moreno, A. Robert, Optical Engineering 39 (2000) 3301. [20] A. Márquez, C. Iemmi, I. Moreno, J.A. Davis, J. Campos, M.J. Yzuel, Optical Engineering 40 (2001) 2558.
Contents lists available at ScienceDirect
Optics Communications journal homepage: www.elsevier.com/locate/optcom
Optical multiresolution analysis with spatial localization Javier Mazzaferri 1, Silvia Ledesma * Image Processing Laboratory (LPI), Physics Department, ‘J.J. Giambiagi’, FCEyN, University of Buenos Aires, Argentina
a r t i c l e
i n f o
Article history: Received 29 September 2009 Received in revised form 7 January 2010 Accepted 10 January 2010
Keywords: Image processing Fourier optics Multiresolution analysis Feature extraction
a b s t r a c t Multiresolution analysis is very useful for characterization of textures, segmentation tasks, and feature enhancement. The development of optical methods to perform such procedures is highly promissory for real-time applications. Usually, the optical implementations of multiresolution analysis consist in the decomposition of the input scene in different frequency bands, obtaining various filtered versions of the scene. However, under certain circumstances it could be useful to provide just one version of the scene where the different filters are applied in different regions. This procedure could be specially interesting for biological and medical applications in situations when the approximate localization of the scale information is known a priori. In this paper we present a fully optical method to perform multiresolution analysis with spatial localization. By means of the proposed technique, the multi-scale analysis is performed at once in a unique image. The experimental set-up consist of a double-pass convergent optical processor. The first stage of the device allows the multiple band decomposition, while the second stage confines the information of each band to different regions of the object and recombines it to achieve the desired operation. Numerical simulations and experimental results, which prove the very good performance of the method, are presented. Ó 2010 Elsevier B.V. All rights reserved.
1. Introduction The interest in studying and developing optical architectures to perform real-time image processing, has notoriously grown over recent years; relevant examples of these applications can be found in recent publications [1–3]. For instance, these kind of devices have been implemented in combination with optical microscopes to carry out specific operations, by only optical methods [4,5]. Likewise, the optical processing techniques could improve the in vivo measurements, present in medical applications. In particular, feature extraction of medical images could help for the diagnosis of some diseases [6]. Moreover, in biological images it could be of interest to follow in real-time the evolution of a specific feature. The wavelet transforms have been extensively employed for multiresolution analysis, showing excellent capabilities. Furthermore, these techniques have been widely used for different image processing applications, such as: texture analysis, feature extraction and image segmentation among others [7–10]. This method has been employed in diverse technological and scientific applications related to biological [11–13] and medical [14,6] image treatment.
* Corresponding author. E-mail address: [email protected] (S. Ledesma). 1 Present address: Instituto de Ciencia de Materiales de Aragón (ICMA), Consejo Superior de Investigaciones Científicas-Universidad de Zaragoza, Facultad de Ciencias, Plaza San Francisco s/n, 50009 Zaragoza, Spain. 0030-4018/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.optcom.2010.01.018
It is a common practice, in multiresolution analysis, to apply a bank of spatial filters and obtain a set of images at different spatial resolutions. In these cases the whole function is filtered to each scale. In fact, when wavelet transforms are employed to perform this processing, it is usual to present the results according to typical arrangements of the different filtered images [15]. The information is displayed separately in different images, each one corresponding to an individual scale parameter. However, for certain applications it could be of interest to analyze one region of the image in one frequency band, and other regions of the image in other frequency bands, displaying the result in a unique final image. The implementation of such a procedure could be simple to perform digitally, however in this case, it would not be done in real-time. Several experimental architectures have been proposed to implement the multiresolution analysis optically. Clear examples of this are two pioneering works presented by Mendlovic and Ouzieli [16,17], where the Dammann grating and computer-generated multireference matched filters were employed to achieve the direct and inverse 2D wavelet transform optically. In all these cases, the information corresponding to each scale is obtained in separated images. As it was mentioned above, until now the procedures needed to obtain a unique image with different scales in different regions of the object should be achieved with digital operations. The novelty of the method presented here is that it performs this procedure in real-time with a fully optical process. We propose a device based on the wavelet transform that can achieve optical
2057
J. Mazzaferri, S. Ledesma / Optics Communications 283 (2010) 2056–2060
spatially-localized multiresolution analysis. The design consists in a double-pass optical correlator, in which the frequency decomposition is done with a bank of tuning filters which are binarized versions of the mexican hat wavelet. The paper is organized as follows: in Section 2 the proposed method is presented as well as the theoretical fundamentals of the technique are exposed. The experimental details are explained and both the numerical simulation and experimental results are presented in Section 3. In Section 4 the conclusions are summarized. 2. Localized multiresolution analysis The optical multiresolution analysis method presented in this work is based on a convergent multichannel processor [2]. The design proposed here, consists in a double-pass correlator which is sketched in Fig. 1. The monochromatic light source S and lens L1 shape a convergent wave-front that impinges the object f. The amplitude of the optic field in the plane conjugated with S through L1 is a scaled version of the Fourier transform F of f. The phase of the field in this plane is quadratic and represents a defocus that can be compensated by properly placing the lens L2 and the set of binary masks BM to achieve a focused result. The filter H, designed to match the scale of F, is placed in the Fourier plane mentioned above. Lens L2 images the filtered version I just in front of BM, completing the first stage of the processor. This step allows to separate the information of the object in individual frequency bands. A set of replicas of the input object, each one filtered at a particular scale, is obtained in different regions of I. By means of BM, a different part of each replica is conserved while the rest is blocked. The light direction is reversed with mirror M1 and the second stage of the device, composed by lens L2, filter H, beam splitter BS and mirror M2, recombines the masked components of I to compose the final image s. It should be pointed out that the masking operation introduces some distortions in the final image. However we will demonstrate with some approximations, that in typical working conditions these effects may be neglected. The multiple filter H, which is analogous to the filter proposed in Ref. [2], is composed by N different filters Hj(u, v), where u and v are the spatial frequencies corresponding to the spatial Cartesian coordinates x and y of the object [2]. Each filter Hj is multiplied by a linear phase factor as is explicitly shown below
Hðu; v Þ ¼
N X
i/jL
Hj ðu; v Þe ;
ð1Þ
j¼1
where f/jL g is a set of linear phases defined as /jL ¼ 2pðukj þ v nj Þ with different values of the pair (kj,nj) for each filter Hj. If an object f is processed with such a filter, the image in the correlation plane (x2,y2) can be described as
Iðx2 ; y2 Þ ¼ F1 ½Fðu; v Þ Hðu; v Þ ¼
N X
f hj dðx2 kj ; y2 nj Þ;
where F1 ½ stands for the inverse Fourier transform, and denotes the convolution operation. In each term of the sum in Eq. (2) there is the convolution between the input scene and the impulse response hj of each filter Hj. The convolutions with the Dirac’s delta functions d(x, y) indicate that each filtered image is displaced off-axis in the intermediate plane. In this work, we have used a set of filters {Hj} based on the Mexican hat wavelet
Wðu; v Þ ¼ 4p2 ðu2 þ v 2 Þ exp 2pðu2 þ v 2 Þ ;
ð3Þ
expressed in the Fourier domain, which is frequently employed for feature extraction and multiresolution analysis. The filters of the bank are based on different scaled versions of Eq. (3). As we are mainly interested in the pass-band feature of these functions, we have used the binarized versions WB(u, v) described in Eq. (4). To physically represent H we choose the phase modulation mode of a SLM, and to this end the amplitude of the filter was phase-encoded [2]. The binary filters are described by
( W B ðu; v Þ ¼
1 ! 0
C D2
2
2 < u2 þ v 2 < C þ D2
! elsewhere
;
ð4Þ
where D and C are the width and centre of the annular-shaped frequency band. As we have mentioned above, we employed the bank of filters {Hj} described by Eq. (4), where the frequency bands of each filter has no intersection with the frequency bands of the other filters of the bank. As it was stated in Eq. (2), the function I(x2, y2) has the distribution sketched in Fig. 2a. A mirror M1, placed in the intermediate plane (x2, y2), conjugates the wave-front and inverts the direction of propagation. This plane acts as the input of the second stage of the device. The wavefront emerging from M1 goes through L2, and its Fourier transform is obtained before the filter. For the second stage, the filter transmission can be described by H(u, v). The optical field emerging from the filter, can be expressed as:
F½I Hðu; v Þ ¼ F ðu; v ÞH ðu; v Þ Hðu; v Þ ¼ F ðu; v ÞjHðu; v Þj2 :
ð5Þ
As H(u, v) is a phase only function, it has constant unitary amplitude, therefore the field distribution of Eq. (5) can be simplified to F*(u, v). By means of the beam splitter BS and a second mirror
a Fig. 1. Proposed optical device for spatially-localized multiresolution analysis.
ð2Þ
j¼1
b
Fig. 2. (a) Sketch of the intermediate image obtained with the first stage of the double-pass correlator, (b) arrangement of binary masks in the intermediate plane.
2058
J. Mazzaferri, S. Ledesma / Optics Communications 283 (2010) 2056–2060
M2, the inverse Fourier transform of Eq. (5) is obtained at the end of the second stage of the device. Therefore, if no binary mask is placed in M1, the intensity at the output plane results:
jsðx3 ; y3 Þj2 ¼ jF1 ½F ðu; v Þj2 ¼ jf ðx1 ; y1 Þj2 :
ð6Þ
It should be mentioned that the coordinates (x1,y1) and (x3, y3) are linearly related by a constant factor, which characterizes the magnification of the system. In what follows, and for the sake of clarity, we will employ the coordinates (x, y) to describe the spatial domain. Eq. (6) shows that the second stage counteracts the effect of the first stage. Let us now consider the effect of positioning a binary mask BM, composed by the elements {mj}, in the intermediate plane (x2, y2). The stops are placed in such a way that each convolution f hj in I is masked by a binary element mj which, in turn, does not overlap with the others. The spatial arrangement of the masks is sketched in Fig. 2b. If BM is placed, the field distribution after the mask, in the second stage, can be described by N nh X
i o f hj mj dðx2 kj ; y2 nj Þ;
ð7Þ
j¼1
N X
hj ðx; yÞ ðmj ðx; yÞf ðx; yÞÞ;
ð13Þ
j¼1
where we have used Eq. (4) to obtain that hj hj = hj. Employing the same approximation for hj ðx; yÞ, it is possible to write the field distribution in the zero order of the final plane as N X
mj hj f :
ð14Þ
j¼1
Finally, as the set of functions {mj} are binary and do not overlap, the intensity of the optical field in the final plane can be described as
jsj2 ¼
N X
2 mj f hj :
ð15Þ
j¼1
where it was used that both d and mj are real functions. The optical field in the Fourier plane, is described as N n h io X M j ðu; v Þ F ðu; v Þ Hj ðu; v Þ
ð8Þ
j¼1
e
support region of mj, and mj does not have variations at the scale of hj. Therefore the integrand of Eq. (12) is not null only if r x and s y. Within this approximation, we can replace mj(x, y) by mj(r, s) and obtain
j
i/L ðu;v Þ
:
3. Results
In Eq. (8), Mj denotes the Fourier transform of mj. The optical field after passing through the filter can be written as N X N n h io X M j ðu; v Þ F ðu; v Þ Hj ðu; v Þ j¼1
ð9Þ
l¼1
Hl ðu; v Þe
ið/jL ðu;v Þþ/lL ðu;v ÞÞ
:
By employing the definition of /L, the exponential factor becomes 1 in those terms where j = l. Therefore, the double summation in Eq. (9) can be split in two summations. One of them corresponding to j = l and the other to j – l. The terms with j = l contribute to the zero order in the final plane, meanwhile the terms with j – l yield offaxis field distributions due to the linear phase factors in the exponentials. Besides, as the filters Hj ðu; v Þ and Hl(u,v) do not intersect, the terms with j – l vanish to zero. Then, the terms that contribute to the final plane are only in the zero order region and can be described as N n h io X M j ðu; v Þ F ðu; v Þ Hj ðu; v Þ
ð10Þ
j¼1
Hj ðu; v Þ: By means of the inverse Fourier transform, the complex amplitude of the field in the final plane can be calculated as: N X
Each term in the sum of Eq. (15) is the product between a binary mask mj and a different squared convolution jf hjj2. Therefore, the whole sum represent different regions of the input object, where each one is filtered to a different scale. This is the desired result.
n h io hj mj f hj :
We have tested the performance of the proposed method with an input scene consisting in a variation of the first three steps of generation of the Sierpinsky carpet, as is shown in Fig. 3a. This scene has the most important features at spatial frequen1 1 and 36 . To analyze this function, cies (in units of 1/pix) around 12 we have employed a filter bank {H1(u, v),H2(u, v)} composed by two binary functions as described by Eq. (4). H1 filters a frequency 1 1 , and H2 filters a frequency band centred at 36 . The band centred at 12 binary filters were obtained by choosing a threshold value that produces frequency bands as wide as possible without overlap. The resulting D values, in units of 1/pix, for H1 and H2 are approx1 1 and 36 , respectively. With these functions, we have deimately 12 signed a multiple phase filter as described in Eq. (1), which was physically implemented by encoding the amplitude information in the phase [2,18]. The linear phases employed to build the filter H were selected to redirect the information filtered by H1 and H2 off-axis. The light corresponding to H1 was redirected in the horizontal direction meanwhile the light filtered with H2 was deviated in the vertical direction. The information corresponding to spatial frequencies outside the bands of H1 and H2 is obtained in the zero order of the intermediate plane. The resulting filter H is shown in Fig. 3b. By means of numerical simulations we have computed the output field distribution of the first stage of the device, which is shown in Fig. 4a. In the figure, three orders can be distinguished.
ð11Þ
j¼1
By explicitly writing the convolution f hj , and introducing mj(x, y) in the integral, it is obtained that N Z Z X
1
1
j¼1
mj ðx; yÞf ðr; sÞhj ðx r; y sÞdrds
ð12Þ
hj ðx; yÞ: As it was previously pointed out, in most of the cases of interest, Eq. (12) can be approximated by a simpler expression. Indeed, generally, the support region of hj is very much smaller than the
Fig. 3. (a) Input scene, (b) multiple phase filter (gray levels from black to white denote phase values ranging from 0 to 2p).
J. Mazzaferri, S. Ledesma / Optics Communications 283 (2010) 2056–2060
2059
m2
m0
m1
a
b
Fig. 4. (a) Field distribution in the intermediate plane of the device, (b) binary mask.
The order on the left is the input scene filtered with H1, the upper order is the same scene filtered with H2, and in the zero order there is the information of the input scene that does not belong to the frequency bands of the filters. To test our method, we have designed a binary mask that blocks the different regions of the intermediate field of Fig. 4(a). The mask is designed to filter in different frequency bands the information of two different parts of the objects. In the smaller upper part of the scene, no filter is applied. Besides, in the bigger lower part of the object, H2 is used. The designed mask is shown in Fig. 4b. As it can be observed, the mask blocks the lower region of the left and zero orders. Therefore, the only information in this part of the object that reaches the final plane, belongs to the upper order. Within the upper part of the scene, no order is blocked, therefore the object is fully reconstructed there. Fig. 6a shows the numerical simulation of the final plane where the desired results can be observed. In order to test the experimental performance of the method, the optical device of Fig. 5 was mounted. As can be noted, the system follows the sketch presented in Fig. 1. The light of an Argon laser is filtered and expanded to act as the point light source of the device. The input scene f(x, y) and the filter H(u, v) are displayed in amplitude (SLM1) and phase (SLM2) spatial modulators respectively. Lens L1 is positioned to obtain the Fourier transform of f(x,y) in the plane where H(u, v) is placed. The SLM’s consist of LCD panels and a set of wave quarter plates (WQP) and linear polarizers (P), that manages the incident and emerging elliptical polarization, as is sketched in the right of Fig. 5. By means of controlling the incoming and emerging polarization states, the desired modulation is obtained [19,20]. The mirror M1 on the right of the system inverts the direction of propagation of the light, and the binary mask BM is positioned in the image plane of the input object, in contact with M1. The second stage of the double correlator is built with the beam splitter BS and
L3 M2
CCD SLM
SLM1
BS
SLM2
WQP1 WQP2
BM
Ar Laser
M1 L1
L2 Video projector
P1
PC
Fig. 5. Double-pass experimental device.
LCD
P2
Fig. 6. (a) Numerical simulation result obtained with the binary mask of Fig. 4b. Experimental results obtained: (b) without binary mask, (c) blocking the 0 and 2 orders, (d) blocking the 0 and 1 orders, (e) with the binary mask of Fig. 4b.
the mirror M2. Lens L3 is introduced to make the output image fit the size of the CCD sensor. It should be noted that to have the filter positioned in the Fourier planes of both the first and the second stage of the device, lens L2 was placed to have its focal plane in the Fourier plane of the first stage. With the input scene and the Fourier filter shown in Fig. 3, experimental results were obtained. As a first test, the output image without any binary mask in the intermediate plane was obtained. As it was predicted in Eq. (6), the image of the input scene was obtained in the final plane. Fig. 6b shows the experimental image. By means of a binary mask, the upper and zero order in the intermediate plane were blocked. In this case, the output of the system showed the input scene filtered with H1, as can be observed in Fig. 6c. Similarly, with another mask, the left and zero orders of the intermediate plane were blocked obtaining, in this case, the input scene filtered with H2. The corresponding output of the system is shown in Fig. 6d. Finally, to perform the spatially-localized multiresolution analysis, the mask presented in Fig. 4b was placed in the intermediate plane of the system. Fig. 6e shows the experimental result of the localized multiresolution analysis, which coincides very well with the numerical simulations presented in Fig. 6a. Despite the optical implementation could be further optimized, excellent experimental results were obtained. It was demonstrated that the presented system permits the full-optical realization of spatially localized multiresolution analysis, which has promissory applications. 4. Conclusions In this work, we propose an optical arrangement to perform spatially-localized multiresolution analysis. The proposed device allows to process different subregions of the input scene at different scales, avoiding to split the scale information in separated images. The whole process is performed optically, and therefore, in real-time. The fundamental interest of this kind of devices is the fact that they can be integrated in other optical systems, as for example optical microscopes, allowing the realization of fulloptical real-time image processing. The proposed method is based in a double-pass multichannel correlator that is based on wavelet multiresolution analysis. The first stage of the device performs the sub-band decomposition and, by means of a set of binary masks, the spatial localization of scale information is obtained at the end of the second stage. Numerical simulations and experimental
2060
J. Mazzaferri, S. Ledesma / Optics Communications 283 (2010) 2056–2060
tests were provided to demonstrate the high performance of the method. Very good results were obtained, suggesting promissory applications of the proposed architecture. Acknowledgement This work was supported with grants from Consejo Nacional de Investigaciones Cientı´ficas y Técnicas (CONICET), Universidad de Buenos Aires (UBA) and Agencia Nacional de Promoción Cientı´fica y Tecnológica (ANPCyT) from Argentina. S. Ledesma is a member of CONICET. References [1] B. Javidi, I. Moon, S. Yeom, E. Carapezza, Optics Express 13 (12) (2005) 4492. [2] J.E. Mazzaferri, S.A. Ledesma, C.C. Iemmi, Journal of Optics A: Pure and Applied Optics 5 (2003) 425. [3] H. Zhang, C.M. Cartwright, M.S. Ding, W. Gillespie, Optics Communications 185 (2000) 277. [4] S. Fürhapter, A. Jesacher, S. Bernet, M. Ritsch-Marte, Optics Express 13 (2005) 689.
[5] A. Jesacher, S. Fürhapter, S. Bernet, M. Ritsch-Marte, Physical Review Letters 94 (2005) 233902. [6] W. Qian, L. Li, L.P. Clarke, Medical Physics 26 (3) (1999) 402. [7] C.S. Lu, P.C. Chung, C.E. Chen, Pattern Recognition 30 (5) (1997) 729. [8] N. Fatemi-Ghomi, P. Palmer, M. Petrou, in: Proceedings of the Workshop on Performance Characteristics of Vision Algorithms (ECCV 96), 1996. [9] M. Unser, IEEE Transactions on Image Processing 4 (11) (1995) 1549. [10] T. Chang, C. Kuo, Image Processing IEEE Transactions on 2 (4) (1993) 429. [11] G. Papari, P. Campisi, N. Petkov, A. Neri, EURASIP Journal on Applied Signal Processing 2007 (2007) 119. [12] J.-C. Olivo-Marin, Pattern Recognition 35 (9) (2002) 1989. [13] M. Unser, A. Aldroubi, Proceedings of the IEEE 84 (4) (1996) 626. [14] C. Wählby, I. Sintorn, F. Erlandsson, G. Borgefors, E. Bengtsson, Journal of Microscopy 215 (2004) 6776. [15] S.G. Mallat, S. Zhong, IEEE Transactions on Pattern Analysis and Machine Intelligence 14 (1992) 710. [16] D. Mendlovic, I. Ouzieli, I. Kiryuschev, E. Marom, Applied Optics 34 (1995) 8213. [17] I. Ouzieli, D. Mendlovic, Applied Optics 35 (1996) 5839. [18] I. Moreno, E. Ahouzi, J. Campos, M.J. Yzuel, Applied Optics 36 (1997) 7428. [19] A. Márquez, J. Campos, M.J. Yzuel, I. Moreno, J.A. Davis, C. Iemmi, A. Moreno, A. Robert, Optical Engineering 39 (2000) 3301. [20] A. Márquez, C. Iemmi, I. Moreno, J.A. Davis, J. Campos, M.J. Yzuel, Optical Engineering 40 (2001) 2558.