- Email: [email protected]
Resolving permutation ambiguity in correlation-based blind image separation
Pattern Recognition Letters 33 (2012) 559–567
Contents lists available at SciVerse ScienceDirect
Pattern Recognition Letters journal homepage: www.elsevier.com/locate/patrec
Resolving permutation ambiguity in correlation-based blind image separation Kenji Hara ⇑, Kohei Inoue, Kiichi Urahama Department of Visual Communication Design, Kyushu University, 4-9-1 Shiobaru, Minami-ku, Fukuoka-shi 815–8540, Japan
a r t i c l e
i n f o
Article history: Received 27 December 2010 Available online 8 December 2011 Communicated by G. Moser Keywords: Separation of reflection Blind image separation Permutation ambiguity Generalized multiple correlation Pruning scheme
a b s t r a c t We address the problem of permutation ambiguity in blind separation of multiple mixtures of multiple images (resulting, for instance, from multiple reflections through a thick grass plate or through two overlapping glass plates) with unknown mixing coefficients. In this paper, first we devise a generalized multiple correlation measure between one gray image and a set of multiple gray images and derive a decorrelation-based blind image separation algorithm. However, many blind image separation methods, including this algorithm, suffer from a permutation ambiguity problem that the success of the separation depends upon the selection of permutations corresponding to the orders of the update operations. To solve the problem, we improve the first algorithm above by decorrelating the mixtures while searching for the appropriate update permutation using a pruning technique. We show its effectiveness through experiments with artificially mixed images and real images. Ó 2011 Elsevier B.V. All rights reserved.
1. Introduction When a scene is photographed through a glass plate, a reflection is commonly observed superimposed over the scene. Especially in the case of through a thick grass plate or through two overlapping parallel glass plates, multiple reflections are projected upon the scene. The separation of each reflection and transparency images from the photo(s) has recently received considerable attention in the image processing and computer vision communities. The existing image separation techniques can be broadly classified into two categories: general blind source separation (Cardoso and Souloumiac, 1993; Bell and Sejnowski, 1995; Hyvaerinen and Oja, 1997; Farid and Adelson, 1999; Cardoso et al., 2002; Bronstein et al., 2005; Bedini et al., 2005; Kopriva, 2007; Gribonval and Zibulevsky, 2010) and image-specific separation (Nayar et al., 1996; Lin and Shum, 2001; Levin et al., 2004; Sarel and Irani, 2004; Diamantaras and Papadimitriou, 2005; Tonazzini et al., 2006; Levin and Weiss, 2007; Kayabol et al., 2009; Hara et al., 2009). The former is mainly based on independent component analysis (ICA) (Hyvaerinen and Oja, 1997; Cichocki and Amari, 2002) and is applied more broadly to other types of signals such as acoustics and biomedical signals. The latter does not necessarily use the ICA framework and is often restricted to separation of image mixtures. As examples of studies in the latter category, Levin et al. proposed automatic (Levin et al., 2004) and semi-automatic (Levin and Weiss, 2007) methods for separating layers from only a single mixed image based on a
⇑ Corresponding author. E-mail address: [email protected] (K. Hara). 0167-8655/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.patrec.2011.11.011
sparsity prior on natural images. Diamantaras et al. (Diamantaras and Papadimitriou, 2005) presented a robust and non-iterative method for separating image mixtures through the properties of the mixtures ratio. Sarel et al., in (Sarel and Irani, 2004), successfully separated linear mixtures of images by a decorrelation-based information exchange between mixtures. Tonazzini et al. (2006) accomplished the blind separation of images from noisy linear mixtures by formulating as a Bayesian estimation problem. Especially for the case of mixtures of more than three images, most of existing blind image separation methods suffer from a permutation ambiguity problem that the success of the separation depends upon the selection of permutations corresponding to the orders of the update operations. In this paper, we present a method for resolving this permutation ambiguity in blind image separation. First, Sarel et al.’s method (Sarel and Irani, 2004) will be extended to the blind separation of mixtures of an arbitrary number of images. We devise a generalized multiple correlation measure between one gray image and a set of multiple gray images. This multiple correlation leads us to provide a set of simultaneous linear equations for updating each mixture of images. Then source images are estimated by iterating between solving the sets of equations and cyclically permuting the mixtures of images. However, there does also exist the permutation ambiguity, just as that in the other separation methods. Thus, we further improve the first algorithm above by decorrelating the mixtures while searching for the appropriate update permutation using a pruning technique based on the discrepancy among the separation parameter values for each of the three RGB channels. We show its effectiveness through experiments with artificially mixed images and real images.
560
K. Hara et al. / Pattern Recognition Letters 33 (2012) 559–567
2. Decorrelation-based blind separation of image mixtures In this section, we describe an extension of the two-layer separation method in (Sarel and Irani, 2004), adapted to the multi-layer separation. This algorithm is included in the category of imagespecific separation, and thus we describe it in the specific framework of color image analysis. 2.1. Problem formulation Let ak = [ak,1 . . . , ak,K]T (k = 1, . . . , K) be the kth K-dimensional coefficient vector. Let IRk ; IGk and IBk (respectively LRk ; LGk and LBk ) be the three RGB component images (the image size is M = h w pixels) of the kth color mixed image Ik (respectively the kth color source h i h i image Lk). Let IVk ¼ IRk IGk IBk (respectively LVk ¼ LRk LGk LBk ) (k = 1, . . . , K) be the grayscale image (the image size N = 3M = h 3w pixels) formed by horizontally concatenating IRk ; IGk and IBk (respectively LRk ; LGk and LBk ). Then any kth grayscale mixed image IVk is given by linearly combining K independent grayscale source images LV1 ; . . . ; LVK as,
where Ci(fk, gl) is the covariance between the ith windows in fk and gl, and Vi(fk) (respectively Vi(gl)) is the variance of the ith window in fk (respectively gl). 2.3. Decorrelating multiple mixtures using cyclic iteration In this section, we derive the update equation for estimating the image sources from a given set of image mixtures. In each iterative update we select a reference image cyclically from the image mixtures and then subtruct the other image mixtures from the reference image by weighting them with the related correlation measure MGNGC introduced in the previous section. Note that the weighting parameters need to be determined optimally such that the updated reference image and the image sources are decorrelated. The exact details of how this determination of the weighting parameters is done are described below. T We will find a minimizer, rV ¼ rV1 rVK1 , of the MGNGC P V V e measure between the reference image I ¼ IVK K1 k¼1 rk Ik and V V a set of K 1 the image mixtures I1 ; . . . ; IK1 , as follows (see APPENDIX for derivation).
¼a
V 2;1 L1 ðiÞ
V 2;2 L2 ðiÞ
þa
V 2;K LK ðiÞ;
þ þ a
.. .
¼ argminr1 ;...;rK1 PN i T i 1 i QK1 V U K1 uK1 k¼1 V i Ik i¼1 uK1 ; QK1 V PN e V i¼1 V i ð I Þ k¼1 V i I k
IV1 ðiÞ ¼ a1;1 LV1 ðiÞ þ a1;2 LV2 ðiÞ þ þ a1;K LVK ðiÞ; IV2 ðiÞ
rV1 rVK1 ¼ argminr1 ;...;rK1 MGNGC K IV1 ; . . . ; IVK1 ; eI V
ð1Þ
ð5Þ
IVK ðiÞ ¼ aK;1 LV1 ðiÞ þ aK;2 LV2 ðiÞ þ þ aK;K LVK ðiÞ; ði ¼ 1; . . . ; NÞ; where IVk ðiÞ (respectively LVk ðiÞ) are the values of the ith pixels in IVk (respectively LVk ), for k = 1, . . . , K. For now, for given I1, . . . , IK with unknown a1, . . . , aK, we will estimate (constant times each of) the most likely L1, . . . , LK. 2.2. Generalized multiple correlation By extension of the generalized normalized grayscale correlation (GNGC) measure (Sarel and Irani, 2004), which is a correlation measure between two grayscale images of the same window size, we introduce a multiple correlation measure, which we will refer to as the multiple generalized normalized grayscale correlation (MGNGC) measure, between one grayscale image fK and a set of K 1 grayscale images f1, . . . , fK1 of the same window size N, as
i2 Q PN h K i i¼1 MNGC K ðf1 ; . . . ; fK1 ; fK Þ k¼1 V i ðfk Þ ; MGNGC K ðf1 ; . . . ; fK1 ; fK Þ ¼ P N QK i¼1 k¼1 V i ðfk Þ ð2Þ where MNGC iK ðf1 ; . . . ; fK1 ; f K Þ is the multiple correlation coefficient between the 5 5 sliding window centered at a pixel i in fK and the set of 5 5 sliding windows at the same pixel i in f1, . . . , fK1. It is well known that MNGCK,i(f1, . . . , fK1; fK) is expressed as
MNGC iK ðf1 ; . . . ; fK1 ; fK Þ
1=2 1 ¼ 1 i ; er KK
ð3Þ
where ~riKK is the KKth entry of the inverse matrix of the K dimensional square correlation matrix RiK whose klth entry is given by
RiK ¼ kl
C i ðfk ; g l Þ ½V i ðfk ÞV i ðg l Þ1=2
ðk; l ¼ 1; . . . ; KÞ;
ð4Þ
where U iK1 is the K 1 dimensional square correlation matrix C i ðIVk ;IVl Þ i ¼ whose klth entry is U iK1 1=2 and uK1 is the K 1 kl ½V i ðIVk ÞV i ðIVl Þ C IV ;e IV dimensional vector whose kth entry is uiK1 k ¼ i kV 1=2 ½ V i ðI k Þ Hence, setting the partial derivatives of Eq. (5) with respect to rVk , for k = 1, . . . , K 1, to zero, we get (see APPENDIX for derivation)
" K1 X N X m¼1
¼
K1 Y V C i IVk ; IVm V i Il
i¼1 N X
#
rVm
l¼1
K1 Y V ðk ¼ 1; . . . ; K 1Þ: C i IVk ; IVK V i Il
i¼1
ð6Þ
l¼1
Eq. (6) is equivalent to a set of K 1 linear simultaneous equations T with unknowns rV ¼ rV1 rVK1 . The update equations for rV is obtained analytically as,
rV ¼ ðAV Þ1 bV ;
ð7Þ
V
V
where A is the K 1 dimensional square matrix and b is the K 1 dimensional vector, as
ðAV Þkm ¼
N X i¼1
V
ðb Þk ¼
N X i¼1
1 KY C i IVk ; IVm V i IVl ;
ðk; m ¼ 1; . . . ; K 1Þ;
ð8Þ
l¼1
K1 Y V C i IVk ; IVK V i Il ;
ðk ¼ 1; . . . ; K 1Þ:
ð9Þ
l¼1
Given as input a set of K mixed images I1, . . . , IK, our presented algorithm uses a two-step iterative procedure where the K 1 dimensional vector rVk is updated, and subsequently the K image sequence is cyclicly updated (the group of K successive iterations is hereinafter referred to as a cycle), where the convergence criterion is such that all values of rV1 rVK1 are less than 103, as outlined in Algorithm 1 and Fig. 3.
561
K. Hara et al. / Pattern Recognition Letters 33 (2012) 559–567 Table 1 Interference-to-signal ratio (ISR) of image separation.
Algorithm 1. Initialize IV1 ; . . . ; IVK repeat rV Ü (AV)1bV where
ðAV Þkm ¼
N X
K1 Y V C i IVk ; IVm V i Il
i¼1 V
ðb Þk ¼
N X
l¼1
K1 Y V C i IVk ; IVK V i Il
i¼1
l¼1
ðk; m ¼ 1; . . . ; K 1Þ 3 K1 P V V V I ð r Þ I k k 7 6 K k¼1 7 6 6 V7 7 6 6 I2 7 V 7 6 6 7 I 1 7 6 . 7(6 7 6 6 . 7 .. 7 6 4 . 5 5 4 . V IK IVK1 2
IV1
3
2
until convergence return I1, . . . , IK
Condition number of mixing matrix
JADE
Relative Newton
Algorithm 1 (proposed)
5.37 (well-conditioned) 46.77 124.80 (ill-conditioned)
0.77 0.77 0.77
0.19 0.19 0.19
0.13 0.23 0.59
first step (Eq. (7)) in Algorithm 1 for each color channel image se T quence IC1 ; . . . ; ICK be rC ¼ rC1 rCK1 ðC ¼ R; G; BÞ. With the help of experiments, we confirmed the fact that the equality relation among rR, rG and rB depends upon the update permutation (note that rR rG rB is not always satisfied for all the update permutations) and that Algorithm 1 tends to fail when the permutation is chosen such that the variation among rR, rG and rB is large. This observation leads us to a pruning-based blind separation algorithm, stated in the following. At first, for each of all the possible K! update permutations, the revised algorithm, as presented in Algorithm 2, individually repeats the two-step process in Algorithm 1 for J iterations (or Q cycles with Q = J/K) and then, for the pth update permutation, records the average value
TðpÞ ¼
J 1X T j ðpÞ ðp ¼ 1; . . . ; K!Þ; J j¼1
ð10Þ
where Tj(p) is, for the jth repetition of the two-step process under the pth permutation, For comparison, Table 1 shows the numerical error measures, interference-to-signal ratio (ISR), achieved by Algorithm 1 and two of the most efficient blind source separation methods – JADE (Cardoso and Souloumiac, 1993) and Relative Newton method (Gribonval and Zibulevsky, 2010), under different condition numbers (the ratio between the maximum and minimum eigenvalues) of the mixing matrices. As we can see in Table 1, while JADE and Relative Newton method do not depend on the condition number of mixing matrix,1 Algorithm 1 deteriorates in the ill-conditioned case (condition number = 46.77, 124.80). However, in the well-conditioned case (condition number = 5.37), Algorithm 1 seems to outperform the two other methods. Fig. 1 shows the separation results corresponding to the top row (well-conditioned case) of Table 1. In each of JADE, Relative Newton, and our Algorithm 1, K! permutations are permitted for the K mixtures. For example, in Algorithm 1, it is possible for the resulting images to be negative under some permutations. Algorithm 1 often leads to improper solutions including these negative images. In our preliminary experiment we discovered and confirmed that JADE and Relative Newton also succeed or fail, depending upon the selection of permutations corresponding to the orders of the update operations. In each of Fig. 1(i)–(l), (m)–(p) and (q)–(t), shown previously, the visually determined optimum permutation is used. The failed examples are shown in Fig. 2((a), (c), (d), (g), (h) and (i) seem to be failures). It is shown in the following section that the problem in our Algorithm 1 is easily solved with a pruning scheme.
T j ðpÞ ¼krRj ðpÞ rGj ðpÞk þ krGj ðpÞ rBj ðpÞk þ krBj ðpÞ rRj ðpÞk ðj ¼ 1; . . . ; JÞ;
Algorithm 2
P ( fð1; 2; . . . ; KÞ; ð2; 1; . . . ; KÞ; . . . ; ðK; . . . ; 2; 1Þg |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl} 1
1 It is known that at low noise levels the performance of invariant algorithms such as JADE and Relative Newton method does not depend upon the mixing matrix (Cardoso, 1994).
2
K!
for all p such that p 2 P do Initialize IC1 ðpÞ; . . . ; ICK ðpÞ; C ¼ R; G; B, for p end for repeat for all p such that p 2 P do for j = 1 to J do for all C such that C 2 {R, G, B} do
rCj ðpÞ ( ðAC ðpÞÞ1 bC ðpÞ
3. Simultaneous determination of update permutation In this section we modify Algorithm 1 to decorrelate the mixtures while searching for the appropriate update permutation. Now, with the subscripts V and N(=h 3w) in Algorithm 1 replaced by C and M(=h w), where the subscript C indicates each of the three color channels (R, G, B), let rV obtained by executing the
ð11Þ
where k k denotes K 1 dimensional Euclidean distance (L2 norm) T and rCj ðpÞðC ¼ R; G; BÞ is the value of rC ¼ rC1 rCK1 obtained for the jth repetition of the two-step process under the pth permutation. If TðpÞ exceeds a certain threshold value, the corresponding update permutation p is rejected from the set of the candidate update permutations. The above process is iterated until only one permutation candidate p⁄ remains. On termination, the h i algorithm gives the estimated source images Ip ;1 ¼ IRp ;1 IGp ;1 IBp ;1 ; . . . ; Ip ;K ¼ h i IRp ;K IGp ;K IBp ;K .
3 2 K1 3 P C C C I ðpÞ ð r ðpÞÞ I ðpÞ IC1 ðpÞ K k k 7 6 j 7 6 7 6 C k¼1 7 6 6 I2 ðpÞ 7 C 7 7(6 6 I ðpÞ 1 . 7 6 6 . 7 7 6 4 . 5 .. 5 4 . C IK ðpÞ C IK1 ðpÞ 2
end for
562
K. Hara et al. / Pattern Recognition Letters 33 (2012) 559–567
if jPj > 1 then
T j ðpÞ ( krRj ðpÞ rGj ðpÞk þ krGj ðpÞ rBj ðpÞk þ krBj ðpÞ rRj ðpÞk end if end for if jPj > 1 then P TðpÞ ( 1J Jj¼1 T j ðpÞ end if end for if jPj > 1 then for all p such that p 2 P do if TðpÞ > Threshold then P Ü P {p}
end if end for end if until jPj = 1 and convergence h i I1 ( IR1 ðpÞIG1 ðpÞIB1 ðpÞ ; . . . h i ; IK ( IRK ðpÞIGK ðpÞIBK ðpÞ ; where p 2 P return I1, . . . , IK
4. Results 4.1. Artificial mixtures We present experiments with synthetic mixtures of four known images. Fig. 4(a)–(d) show the source images. We mixed them
Fig. 1. Comparison of Algorithm 1 with the existing methods: (a)–(d) source images, (e)–(h) mixed images (condition number of mixing matrix = 5.37), (i)–(l) layers separated by JADE (ISR = 0.77), (m)–(p) layers separated by Relative Newton (ISR = 0.19), (q)–(t) layers separated by our Algorithm 1 (ISR = 0.13).
K. Hara et al. / Pattern Recognition Letters 33 (2012) 559–567
563
Fig. 2. Examples of failed image separation under a permutation: (a)–(d) JADE, (e)–(h) Relative Newton, (i)–(l) our Algorithm 1.
R
Mixed image sequence
G
B
Three concatenated Color image sequence gray image sequence
Estimated source image sequence Fig. 3. Image separation through iterative image updates.
using different mixing ratios, as shown in Fig. 4(e)–(h). We show the results of our method outlined in Algorithm 2 in Fig. 4(i)–(l). One can see that our approach gives fairly good results. Fig. 5 shows the process of search and selection of the update permutation in this image separation. The vertical axis represents TðpÞ value (Eq. (10)). The horizontal axis represents the number of cycles executed for the mixture separation and selection of the update permutation. According to the criterion for accepting or rejecting update premutation candidates in Algorithm 2, six sequences {(3,1,4,2),(3,4,1,2),(4,3,2,1),(4,1,3,2),(1,3,4,2),(1,4,3,2)}, two sequences {(1,3,4,2), (1,4,3,2)}, and one sequence {(1,4,3,2)}
are selected as candidates after the first, second and third cycles, respectively, as shown by real lines in Fig. 5. Here the dotted lines represent the mixture separation processes for the rejected permutations. The labels 1, 2, 3 and 4 represent the images in Fig. 4(e), (f), (g) and (h), respectively. The circles and black crosses depict the update permutations which give perceptionally correct and incorrect separation results, respectively. We confirm that the finally selected permutation (1, 4, 3, 2) (drawn as a double circle) gives a perceptionally correct separation result, as shown at the right bottom of Fig. 5 (or Fig. 4(m)–(p)). The picture at the right top of Fig. 5 shows an unsuccessful result obtained for one
564
K. Hara et al. / Pattern Recognition Letters 33 (2012) 559–567
Fig. 4. Image separation with our Algorithm 2: (a)–(d) source images, (e)–(h) mixed images, (i)–(l) layers separated by our method.
failure
reject
reject reject success
Fig. 5. Process of mixture separation and selection of the update permutation (artificial mixtures).
of the rejected permutations. Table 2 compares the success ratio between Algorithm 1 and 2. In this experiment, we use a dataset, where each data element (each of No.1–No.10 in the top row in Table 2) contains four image mixtures synthesized using different mixing ratios. First we visually test Algorithm 1 for every permutation for each data element and then give the (number of successful permutations)/(total number of permutations) ratio in the second row in Table 2. Next we visually test Algorithm 2 for the same data element and then show success or failure of the image separation in the bottom row in Table 2. The success rates of Algorithm 1 and 2 are calculated from the above results and are shown in the right column in Table 2. Among the two algorithms, Algorithm 2 seems to achieve higher success ratio.
4.2. Real mixtures We apply our method to layer extraction from images containing multiple reflections and transparency. We photographed a picture postcard in a glass-fronted bookcase (Fig. 6(a)). One can see transparency and double reflections due to light reflected from both the surfaces of the front side and rear side glasses (region surrounded by a rectangle of Fig. 6(a)). We took three photographs under three different illumination conditions by inserting a polarization sheet at the front or back of two glasses on the bookcase, or between them (Fig. 6(c)–(e)). For acquisition of a ground truth transparency image we shot the same scene while shielding out some of the ambient light using a blackout curtain (Fig. 6(b)).
565
K. Hara et al. / Pattern Recognition Letters 33 (2012) 559–567 Table 2 Comparison between Algorithm 1 and 2. Methods
Algorithm 1
Input data number
number of successful permutations total number of permutations
Algorithm 2 (success: s; failure: )
Success ratio (%)
No. 1
No. 2
No. 3
No. 4
No. 5
No. 6
No. 7
No. 8
No. 9
No. 10
12/24
24/24
19/24
7/24
16/24
6/24
23/24
21/24
23/24
2/24
63.5
s
s
s
s
s
s
s
s
s
90.0
The separation results obtained by the JADE (Cardoso and Souloumiac, 1993), Relative Newton method (Gribonval and Zibulevsky, 2010) and our method are shown in Fig. 6(f)–(h), (i)–(k) and (l)–(n), respectively. One can see that, for both Relative Newton method and our method, the estimated transparency images (Fig. 6(i) and (l)) are similar to the ground truth transparency image (Fig. 6(b)), and that Fig. 6(j), (k), (m) and (n) are also relatively clear compared to JADE (Fig. 6(f)–(h)). In the same way as in the previous case, the process of search and selection of the permutation in this experiment is shown in Fig. 7. Here the labels 1, 2 and 3 represent the images in Fig. 6(c), (d) and (e), respectively. The update permutation (2, 3, 1) is selected and it gives a correct separation result, as shown at the center of Fig. 7 (or Fig. 6(l)–(n)). The picture at the right top of Fig. 7 shows an unsuccessful result obtained for one of the rejected permutations. 5. Conclusions In this paper, we have proposed a method of stably recovering a set of the original source images from the input mixtures. In many existing blind image separation techniques the success of the separation depends upon the selection of permutations corresponding to the orders of the update operations. By extending Sarel et al.’s bivariate generalized correlation measure (Sarel and Irani, 2004) to more than two variables and then introducing the simple pruning technique based on the discrepancy among the separation parameter values for each of the three RGB channels, we have presented a image separation method without permutation ambiguity. We have presented the separation algorithm, and shown its effectiveness through experiments with mixtures of real images. Acknowledgements This research was partially supported by the Suzuki Foundation. Appendix A. Derivation of Eq. (5) In this appendix, we derive Eq. (5) in detail. First we divide the correlation matrix RiK (Eq. (4)) into four blocks as,
ðA:1Þ where RiK1 is a K 1 dimensional square matrix and r iK1 is a K 1 dimensional vector. Hence, from the definition of the inverse matrix, ~riK;K included in Eq. (3) is expressed as
~riK;K Fig. 6. Separation of layers from images containing multiple reflections and transparency: (a) image region containing multiple reflections and transparency, (b) transparency image obtained by shielding out some of the ambient light using a blackout curtain, (c)–(e) three images taken under three different illumination conditions, (f)–(h) layers separated by JADE, (i)–(k) layers separated by Relative Newton, (l)–(n) layers separated by our method.
det RiK1 ; ¼ det RiK
ðA:2Þ
where det () denotes the matrix determinant. By first assuming detðRiK1 Þ – 0 and then applying the formula (Chen and Gu, 2000) for the determinant of a block matrix to Eq. (A.1), we have
566
K. Hara et al. / Pattern Recognition Letters 33 (2012) 559–567
failure
success
Fig. 7. Process of mixture separation and selection of the update permutation (real mixtures).
1 det RiK ¼ det RiK1 1 ðriK1 ÞT RiK1 r iK1 ;
ðA:3Þ
@ Wi @ rk
where ()1 denotes the inverse matrix. Using Eq. (A.3), we can rewrite Eq. (A.2) as,
~r iK;K ¼ 1
ðriK1 ÞT
1 : 1 RiK1 riK1
ðA:4Þ
Hence, using Eq. (A.4), we can rewrite Eq. (3) as,
T
T 1 C IV ;eI V C i IVK1 ;e IV i 1 U iK1 1=2 1=2 ½V i ðIV1 Þ ½V i ðIVK1 Þ h i K1 1=2 Q ¼2 V i IVl 0 0 V i IVk 00
ðA:5Þ n o IV1 ; . . . ; IVK1 ; eI V
MGNGC K IV1 ; . . . ; IVK1 ; eI V i2 PN h Q i V V V eV V i ðeI V Þ K1 i¼1 MNGC K I1 ; . . . ; IK1 ; I k¼1 V i Ik ¼ PN V eV QK1 i¼1 V i ð I Þ k¼1 V i Ik PN i T i 1 i QK1 V U K1 uK1 k¼1 V i Ik i¼1 uK1 ¼ ; Q PN V eV K1 i¼1 V i ð I Þ k¼1 V i Ik where
@uiK1 @ rk
l¼1
1=2 1 MNGC iK ðf1 ; . . . ; fK1 ; fK Þ ¼ ðriK1 ÞT RiK1 riK1 : Replacing {f1, . . . , fK1, fK} in Eqs. (2) and (A.5) with and then combining these equations as,
i @u@ W i K1 K1 Q C ðIV ;IV Þ C ðIV ;IV Þ ¼ 2 V i ðIVl Þ i V1 k1=2 i VK1 k1=2 ½V i ðI1 Þ ½V i ðIK1 Þ l¼1 T V V e C I ;I C IV ;e IV ðU iK1 Þ1 i 1V 1=2 i K1 1=2 ½V i ðI1 Þ ½V i ðIVK1 Þ h i h i K1 1=2 Q ¼2 V i IVl V i IVk ui1;k uiK1;k
¼
ðB:2Þ
l¼1
T V C I ;e IV C IV ;e IV i 1V 1=2 i K1 1=2 ½V i ðI1 Þ ½V i ðIVK1 Þ K1 Q V ¼ 2C i eI V ; IVk V i Il l¼1
K1 K1 P Q V ¼ 2 C i IVK ; IVk rm C i IVm ; IVk V i Il : m¼1
ðA:6Þ
U iK1
is the K 1 dimensional square correlation matrix C i ðIVk ;IVl Þ i ¼ whose klth entry is U iK1 1=2 and uK1 is the K 1 kl ½V i ðIVk ÞV i ðIVl Þ C IV ;e IV dimensional vector whose kth entry is uiK1 k ¼ i kV 1=2 , as ½V i ðIk Þ described in Section 2.3.
l¼1
P So, the partial derivatives of the numerator, Ni¼1 Wi , of Eq. (5) with respect to rk, for k = 1, . . . , K 1, are given by N @ X
@ rk
i¼1
Wi ¼ 2
N X
"
# K1 X K1 Y V C i IVk ; IVK rm C i IVk ; IVm V i Il
i¼1
m¼1
l¼1
ðk ¼ 1; . . . ; K 1Þ:
ðB:3Þ
Since one can easily see that for k = 1, . . . , K 1, the partial deriva PN V eV QK1 i¼1 V i ð I Þ k¼1 V i I k , of Eq. (5) with re-
tives of the denominator,
spect to rk is also equal to Eq. (B.3), we have
" # N K 1 1 X KY X @F @G ¼ ¼2 C i IVk ; IVK rm C i IVk ; IVm V i IVl @ rk @ rk m¼1 i¼1 l¼1
Appendix B. Derivation of Eq. (6) Denoting the numerator of Eq. (5) by
Wi ¼ uiK1
T
U iK1
1
uiK1
K 1 Y
V i IVk
PN
i¼1
Wi as,
ðk ¼ 1; . . . ; K 1Þ; ðB:1Þ
k¼1
and then taking partial derivatives of Wi with respect to rk, for k = 1, . . . , K 1, we get
ðB:4Þ
where F and G are the denominator and the numerator of Eq. (5), respectively. Thus, taking partial derivatives of (5) with respect to rk and setting them to zero, we get
K. Hara et al. / Pattern Recognition Letters 33 (2012) 559–567
@ @ F 1 @F @G 0¼ MGNGC K IV1 ; . . . ; IVK1 ; eI V ¼ F ¼ 2 G @ rk @ rk G @ rk @ rk G 1 F @F ðk ¼ 1; . . . ; K 1Þ: ðB:5Þ 1 ¼ G G @ rk Then, from the definition of MGNGC and the independency of fIVk gKk¼1 , we have F/G < 1 Hence, we get @ F/@ rk = 0 in (B.5), which leads to the following equation.
" # N K1 X K1 X Y V @F ¼2 C i IVk ; IVK rm C i IVk ; IVm V i Il @ rk m¼1 i¼1 l¼1 ðk ¼ 1; . . . ; K 1Þ:
ðB:6Þ
This equation can be rewritten as,
" K 1 X N X m¼1
¼
# 1 KY V V V C i Ik ; Im V i Il rVm
i¼1 N X
l¼1
C i IVk ; IVK
i¼1
K1 Y
V i IVl ðk ¼ 1; . . . ; K 1Þ:
ðB:7Þ
l¼1
References Bedini, L., Herranz, D., Salerno, E., Baccigalupi, C., Kuruoglu, E., Tonazzini, A., 2005. Separation of correlated astrophysical sources using multiple-lag data covariance matrices. EURASIP J. Appl. Sig. Proc. 15, 2400–2412. Bell, A.J., Sejnowski, T.J., 1995. An information maximization approach to blind separation and blind deconvolution. Neural Comput. 7 (6), 1129–1159. Bronstein, A.M., Bronstein, M.M., Zibulevsky, M., Zeevi, Y.Y., 2005. Sparse ICA for blind separation of transmitted and reflected images. Internat. J. Imag. Sci. Technol. 15 (1), 84–91. Cardoso, J.F., Souloumiac, A., 1993. Blind beamforming for non gaussian signals. In: IEE Proceedings-F, 140(6), 362–370.
567
Cardoso, J.F., 1994. On the performance of orthogonal source separation algorithms. In: Proc. EUSIPCO, pp. 776–779. Cardoso, J.F., Snoussi, H.H., Delabrouille, J., Patanchon, G., 2002. Blind separation of noisy gaussian stationary sources: Application to cosmic microwave imaging. In: Proc. EUSIPCO, vol. 1, pp. 561–564. Chen, J., Gu, G., 2000. Control-Oriented System Identification: An H-Infinity Approach. John Wiley. Cichocki, A., Amari, S., 2002. Adaptive Blind Signal and Image Processing. John Wiley. Diamantaras, K.I., Papadimitriou, T., 2005. Blind separation of reflections using the image mixtures ratio. In: Proc. IEEE ICIP, pp. 1034–1037. Farid, H., Adelson, E.H., 1999. Separating reflections from images by use of independent components analysis. J. Opt. Soc. America A 16 (9), 2136–2145. Gribonval, R., Zibulevsky, M., 2010. Sparse component analysis. In: Comon, Pierre, Jutten, Christian (Eds.), Handbook of Blind Source Separation: Independent Component Analysis and Applications. ELSEVIER. Hara, K., Inoue, K., Urahama, K., 2009. Separation of layers from images containing multiple reflections and transparency using cyclic permutation. In: Proc. IEEE ICASSP, pp. 1157–1160. Hyvaerinen, A., Oja, E., 1997. A fast fixed-point algorithm for independent component analysis. Neural Comput. 9 (7), 1483–1492. Kayabol, K., Kuruoglu, E.E., Sankur, B., 2009. Image source separation using color channel dependencies. In: Proc. ICA, pp. 499–506. Kopriva, L., 2007. Approach to blind image deconvolution by multiscale subband decomposition and independent component analysis. J. Opt. Soc. America A 24 (4), 973–983. Levin, A., Weiss, Y., 2007. User assisted separation of reflections from a single image using a sparsity prior. IEEE Trans. Pattern. Anal. Machine. Intell. 29 (9), 1647– 1654. Levin, A., Zomet, A., Weiss, Y., 2004. Separating reflections from a single image using local features. In: Proc. IEEE CVPR, pp. 306–313. Lin, S., Shum, H.Y., 2001.Separation of diffuse and specular reflection in color images. In: Proc. IEEE CVPR, pp. 341–346. Nayar, S.K., Fang, X.S., Boult, T., 1996. Separation of reflection components using color and polarization. Internat J. Comput. Vision. 21 (3), 163–186. Sarel, B., Irani, M., 2004. Separating transparent layers through layer information exchange. In: Proc. ECCV, pp. 328–341. Tonazzini, A., Bedini, L., Salerno, E., 2006. A markov model for blind image separation by a mean-field EM algorithm. IEEE Trans. IP. 15 (2), 473–482.
Contents lists available at SciVerse ScienceDirect
Pattern Recognition Letters journal homepage: www.elsevier.com/locate/patrec
Resolving permutation ambiguity in correlation-based blind image separation Kenji Hara ⇑, Kohei Inoue, Kiichi Urahama Department of Visual Communication Design, Kyushu University, 4-9-1 Shiobaru, Minami-ku, Fukuoka-shi 815–8540, Japan
a r t i c l e
i n f o
Article history: Received 27 December 2010 Available online 8 December 2011 Communicated by G. Moser Keywords: Separation of reflection Blind image separation Permutation ambiguity Generalized multiple correlation Pruning scheme
a b s t r a c t We address the problem of permutation ambiguity in blind separation of multiple mixtures of multiple images (resulting, for instance, from multiple reflections through a thick grass plate or through two overlapping glass plates) with unknown mixing coefficients. In this paper, first we devise a generalized multiple correlation measure between one gray image and a set of multiple gray images and derive a decorrelation-based blind image separation algorithm. However, many blind image separation methods, including this algorithm, suffer from a permutation ambiguity problem that the success of the separation depends upon the selection of permutations corresponding to the orders of the update operations. To solve the problem, we improve the first algorithm above by decorrelating the mixtures while searching for the appropriate update permutation using a pruning technique. We show its effectiveness through experiments with artificially mixed images and real images. Ó 2011 Elsevier B.V. All rights reserved.
1. Introduction When a scene is photographed through a glass plate, a reflection is commonly observed superimposed over the scene. Especially in the case of through a thick grass plate or through two overlapping parallel glass plates, multiple reflections are projected upon the scene. The separation of each reflection and transparency images from the photo(s) has recently received considerable attention in the image processing and computer vision communities. The existing image separation techniques can be broadly classified into two categories: general blind source separation (Cardoso and Souloumiac, 1993; Bell and Sejnowski, 1995; Hyvaerinen and Oja, 1997; Farid and Adelson, 1999; Cardoso et al., 2002; Bronstein et al., 2005; Bedini et al., 2005; Kopriva, 2007; Gribonval and Zibulevsky, 2010) and image-specific separation (Nayar et al., 1996; Lin and Shum, 2001; Levin et al., 2004; Sarel and Irani, 2004; Diamantaras and Papadimitriou, 2005; Tonazzini et al., 2006; Levin and Weiss, 2007; Kayabol et al., 2009; Hara et al., 2009). The former is mainly based on independent component analysis (ICA) (Hyvaerinen and Oja, 1997; Cichocki and Amari, 2002) and is applied more broadly to other types of signals such as acoustics and biomedical signals. The latter does not necessarily use the ICA framework and is often restricted to separation of image mixtures. As examples of studies in the latter category, Levin et al. proposed automatic (Levin et al., 2004) and semi-automatic (Levin and Weiss, 2007) methods for separating layers from only a single mixed image based on a
⇑ Corresponding author. E-mail address: [email protected] (K. Hara). 0167-8655/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.patrec.2011.11.011
sparsity prior on natural images. Diamantaras et al. (Diamantaras and Papadimitriou, 2005) presented a robust and non-iterative method for separating image mixtures through the properties of the mixtures ratio. Sarel et al., in (Sarel and Irani, 2004), successfully separated linear mixtures of images by a decorrelation-based information exchange between mixtures. Tonazzini et al. (2006) accomplished the blind separation of images from noisy linear mixtures by formulating as a Bayesian estimation problem. Especially for the case of mixtures of more than three images, most of existing blind image separation methods suffer from a permutation ambiguity problem that the success of the separation depends upon the selection of permutations corresponding to the orders of the update operations. In this paper, we present a method for resolving this permutation ambiguity in blind image separation. First, Sarel et al.’s method (Sarel and Irani, 2004) will be extended to the blind separation of mixtures of an arbitrary number of images. We devise a generalized multiple correlation measure between one gray image and a set of multiple gray images. This multiple correlation leads us to provide a set of simultaneous linear equations for updating each mixture of images. Then source images are estimated by iterating between solving the sets of equations and cyclically permuting the mixtures of images. However, there does also exist the permutation ambiguity, just as that in the other separation methods. Thus, we further improve the first algorithm above by decorrelating the mixtures while searching for the appropriate update permutation using a pruning technique based on the discrepancy among the separation parameter values for each of the three RGB channels. We show its effectiveness through experiments with artificially mixed images and real images.
560
K. Hara et al. / Pattern Recognition Letters 33 (2012) 559–567
2. Decorrelation-based blind separation of image mixtures In this section, we describe an extension of the two-layer separation method in (Sarel and Irani, 2004), adapted to the multi-layer separation. This algorithm is included in the category of imagespecific separation, and thus we describe it in the specific framework of color image analysis. 2.1. Problem formulation Let ak = [ak,1 . . . , ak,K]T (k = 1, . . . , K) be the kth K-dimensional coefficient vector. Let IRk ; IGk and IBk (respectively LRk ; LGk and LBk ) be the three RGB component images (the image size is M = h w pixels) of the kth color mixed image Ik (respectively the kth color source h i h i image Lk). Let IVk ¼ IRk IGk IBk (respectively LVk ¼ LRk LGk LBk ) (k = 1, . . . , K) be the grayscale image (the image size N = 3M = h 3w pixels) formed by horizontally concatenating IRk ; IGk and IBk (respectively LRk ; LGk and LBk ). Then any kth grayscale mixed image IVk is given by linearly combining K independent grayscale source images LV1 ; . . . ; LVK as,
where Ci(fk, gl) is the covariance between the ith windows in fk and gl, and Vi(fk) (respectively Vi(gl)) is the variance of the ith window in fk (respectively gl). 2.3. Decorrelating multiple mixtures using cyclic iteration In this section, we derive the update equation for estimating the image sources from a given set of image mixtures. In each iterative update we select a reference image cyclically from the image mixtures and then subtruct the other image mixtures from the reference image by weighting them with the related correlation measure MGNGC introduced in the previous section. Note that the weighting parameters need to be determined optimally such that the updated reference image and the image sources are decorrelated. The exact details of how this determination of the weighting parameters is done are described below. T We will find a minimizer, rV ¼ rV1 rVK1 , of the MGNGC P V V e measure between the reference image I ¼ IVK K1 k¼1 rk Ik and V V a set of K 1 the image mixtures I1 ; . . . ; IK1 , as follows (see APPENDIX for derivation).
¼a
V 2;1 L1 ðiÞ
V 2;2 L2 ðiÞ
þa
V 2;K LK ðiÞ;
þ þ a
.. .
¼ argminr1 ;...;rK1 PN i T i 1 i QK1 V U K1 uK1 k¼1 V i Ik i¼1 uK1 ; QK1 V PN e V i¼1 V i ð I Þ k¼1 V i I k
IV1 ðiÞ ¼ a1;1 LV1 ðiÞ þ a1;2 LV2 ðiÞ þ þ a1;K LVK ðiÞ; IV2 ðiÞ
rV1 rVK1 ¼ argminr1 ;...;rK1 MGNGC K IV1 ; . . . ; IVK1 ; eI V
ð1Þ
ð5Þ
IVK ðiÞ ¼ aK;1 LV1 ðiÞ þ aK;2 LV2 ðiÞ þ þ aK;K LVK ðiÞ; ði ¼ 1; . . . ; NÞ; where IVk ðiÞ (respectively LVk ðiÞ) are the values of the ith pixels in IVk (respectively LVk ), for k = 1, . . . , K. For now, for given I1, . . . , IK with unknown a1, . . . , aK, we will estimate (constant times each of) the most likely L1, . . . , LK. 2.2. Generalized multiple correlation By extension of the generalized normalized grayscale correlation (GNGC) measure (Sarel and Irani, 2004), which is a correlation measure between two grayscale images of the same window size, we introduce a multiple correlation measure, which we will refer to as the multiple generalized normalized grayscale correlation (MGNGC) measure, between one grayscale image fK and a set of K 1 grayscale images f1, . . . , fK1 of the same window size N, as
i2 Q PN h K i i¼1 MNGC K ðf1 ; . . . ; fK1 ; fK Þ k¼1 V i ðfk Þ ; MGNGC K ðf1 ; . . . ; fK1 ; fK Þ ¼ P N QK i¼1 k¼1 V i ðfk Þ ð2Þ where MNGC iK ðf1 ; . . . ; fK1 ; f K Þ is the multiple correlation coefficient between the 5 5 sliding window centered at a pixel i in fK and the set of 5 5 sliding windows at the same pixel i in f1, . . . , fK1. It is well known that MNGCK,i(f1, . . . , fK1; fK) is expressed as
MNGC iK ðf1 ; . . . ; fK1 ; fK Þ
1=2 1 ¼ 1 i ; er KK
ð3Þ
where ~riKK is the KKth entry of the inverse matrix of the K dimensional square correlation matrix RiK whose klth entry is given by
RiK ¼ kl
C i ðfk ; g l Þ ½V i ðfk ÞV i ðg l Þ1=2
ðk; l ¼ 1; . . . ; KÞ;
ð4Þ
where U iK1 is the K 1 dimensional square correlation matrix C i ðIVk ;IVl Þ i ¼ whose klth entry is U iK1 1=2 and uK1 is the K 1 kl ½V i ðIVk ÞV i ðIVl Þ C IV ;e IV dimensional vector whose kth entry is uiK1 k ¼ i kV 1=2 ½ V i ðI k Þ Hence, setting the partial derivatives of Eq. (5) with respect to rVk , for k = 1, . . . , K 1, to zero, we get (see APPENDIX for derivation)
" K1 X N X m¼1
¼
K1 Y V C i IVk ; IVm V i Il
i¼1 N X
#
rVm
l¼1
K1 Y V ðk ¼ 1; . . . ; K 1Þ: C i IVk ; IVK V i Il
i¼1
ð6Þ
l¼1
Eq. (6) is equivalent to a set of K 1 linear simultaneous equations T with unknowns rV ¼ rV1 rVK1 . The update equations for rV is obtained analytically as,
rV ¼ ðAV Þ1 bV ;
ð7Þ
V
V
where A is the K 1 dimensional square matrix and b is the K 1 dimensional vector, as
ðAV Þkm ¼
N X i¼1
V
ðb Þk ¼
N X i¼1
1 KY C i IVk ; IVm V i IVl ;
ðk; m ¼ 1; . . . ; K 1Þ;
ð8Þ
l¼1
K1 Y V C i IVk ; IVK V i Il ;
ðk ¼ 1; . . . ; K 1Þ:
ð9Þ
l¼1
Given as input a set of K mixed images I1, . . . , IK, our presented algorithm uses a two-step iterative procedure where the K 1 dimensional vector rVk is updated, and subsequently the K image sequence is cyclicly updated (the group of K successive iterations is hereinafter referred to as a cycle), where the convergence criterion is such that all values of rV1 rVK1 are less than 103, as outlined in Algorithm 1 and Fig. 3.
561
K. Hara et al. / Pattern Recognition Letters 33 (2012) 559–567 Table 1 Interference-to-signal ratio (ISR) of image separation.
Algorithm 1. Initialize IV1 ; . . . ; IVK repeat rV Ü (AV)1bV where
ðAV Þkm ¼
N X
K1 Y V C i IVk ; IVm V i Il
i¼1 V
ðb Þk ¼
N X
l¼1
K1 Y V C i IVk ; IVK V i Il
i¼1
l¼1
ðk; m ¼ 1; . . . ; K 1Þ 3 K1 P V V V I ð r Þ I k k 7 6 K k¼1 7 6 6 V7 7 6 6 I2 7 V 7 6 6 7 I 1 7 6 . 7(6 7 6 6 . 7 .. 7 6 4 . 5 5 4 . V IK IVK1 2
IV1
3
2
until convergence return I1, . . . , IK
Condition number of mixing matrix
JADE
Relative Newton
Algorithm 1 (proposed)
5.37 (well-conditioned) 46.77 124.80 (ill-conditioned)
0.77 0.77 0.77
0.19 0.19 0.19
0.13 0.23 0.59
first step (Eq. (7)) in Algorithm 1 for each color channel image se T quence IC1 ; . . . ; ICK be rC ¼ rC1 rCK1 ðC ¼ R; G; BÞ. With the help of experiments, we confirmed the fact that the equality relation among rR, rG and rB depends upon the update permutation (note that rR rG rB is not always satisfied for all the update permutations) and that Algorithm 1 tends to fail when the permutation is chosen such that the variation among rR, rG and rB is large. This observation leads us to a pruning-based blind separation algorithm, stated in the following. At first, for each of all the possible K! update permutations, the revised algorithm, as presented in Algorithm 2, individually repeats the two-step process in Algorithm 1 for J iterations (or Q cycles with Q = J/K) and then, for the pth update permutation, records the average value
TðpÞ ¼
J 1X T j ðpÞ ðp ¼ 1; . . . ; K!Þ; J j¼1
ð10Þ
where Tj(p) is, for the jth repetition of the two-step process under the pth permutation, For comparison, Table 1 shows the numerical error measures, interference-to-signal ratio (ISR), achieved by Algorithm 1 and two of the most efficient blind source separation methods – JADE (Cardoso and Souloumiac, 1993) and Relative Newton method (Gribonval and Zibulevsky, 2010), under different condition numbers (the ratio between the maximum and minimum eigenvalues) of the mixing matrices. As we can see in Table 1, while JADE and Relative Newton method do not depend on the condition number of mixing matrix,1 Algorithm 1 deteriorates in the ill-conditioned case (condition number = 46.77, 124.80). However, in the well-conditioned case (condition number = 5.37), Algorithm 1 seems to outperform the two other methods. Fig. 1 shows the separation results corresponding to the top row (well-conditioned case) of Table 1. In each of JADE, Relative Newton, and our Algorithm 1, K! permutations are permitted for the K mixtures. For example, in Algorithm 1, it is possible for the resulting images to be negative under some permutations. Algorithm 1 often leads to improper solutions including these negative images. In our preliminary experiment we discovered and confirmed that JADE and Relative Newton also succeed or fail, depending upon the selection of permutations corresponding to the orders of the update operations. In each of Fig. 1(i)–(l), (m)–(p) and (q)–(t), shown previously, the visually determined optimum permutation is used. The failed examples are shown in Fig. 2((a), (c), (d), (g), (h) and (i) seem to be failures). It is shown in the following section that the problem in our Algorithm 1 is easily solved with a pruning scheme.
T j ðpÞ ¼krRj ðpÞ rGj ðpÞk þ krGj ðpÞ rBj ðpÞk þ krBj ðpÞ rRj ðpÞk ðj ¼ 1; . . . ; JÞ;
Algorithm 2
P ( fð1; 2; . . . ; KÞ; ð2; 1; . . . ; KÞ; . . . ; ðK; . . . ; 2; 1Þg |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl} 1
1 It is known that at low noise levels the performance of invariant algorithms such as JADE and Relative Newton method does not depend upon the mixing matrix (Cardoso, 1994).
2
K!
for all p such that p 2 P do Initialize IC1 ðpÞ; . . . ; ICK ðpÞ; C ¼ R; G; B, for p end for repeat for all p such that p 2 P do for j = 1 to J do for all C such that C 2 {R, G, B} do
rCj ðpÞ ( ðAC ðpÞÞ1 bC ðpÞ
3. Simultaneous determination of update permutation In this section we modify Algorithm 1 to decorrelate the mixtures while searching for the appropriate update permutation. Now, with the subscripts V and N(=h 3w) in Algorithm 1 replaced by C and M(=h w), where the subscript C indicates each of the three color channels (R, G, B), let rV obtained by executing the
ð11Þ
where k k denotes K 1 dimensional Euclidean distance (L2 norm) T and rCj ðpÞðC ¼ R; G; BÞ is the value of rC ¼ rC1 rCK1 obtained for the jth repetition of the two-step process under the pth permutation. If TðpÞ exceeds a certain threshold value, the corresponding update permutation p is rejected from the set of the candidate update permutations. The above process is iterated until only one permutation candidate p⁄ remains. On termination, the h i algorithm gives the estimated source images Ip ;1 ¼ IRp ;1 IGp ;1 IBp ;1 ; . . . ; Ip ;K ¼ h i IRp ;K IGp ;K IBp ;K .
3 2 K1 3 P C C C I ðpÞ ð r ðpÞÞ I ðpÞ IC1 ðpÞ K k k 7 6 j 7 6 7 6 C k¼1 7 6 6 I2 ðpÞ 7 C 7 7(6 6 I ðpÞ 1 . 7 6 6 . 7 7 6 4 . 5 .. 5 4 . C IK ðpÞ C IK1 ðpÞ 2
end for
562
K. Hara et al. / Pattern Recognition Letters 33 (2012) 559–567
if jPj > 1 then
T j ðpÞ ( krRj ðpÞ rGj ðpÞk þ krGj ðpÞ rBj ðpÞk þ krBj ðpÞ rRj ðpÞk end if end for if jPj > 1 then P TðpÞ ( 1J Jj¼1 T j ðpÞ end if end for if jPj > 1 then for all p such that p 2 P do if TðpÞ > Threshold then P Ü P {p}
end if end for end if until jPj = 1 and convergence h i I1 ( IR1 ðpÞIG1 ðpÞIB1 ðpÞ ; . . . h i ; IK ( IRK ðpÞIGK ðpÞIBK ðpÞ ; where p 2 P return I1, . . . , IK
4. Results 4.1. Artificial mixtures We present experiments with synthetic mixtures of four known images. Fig. 4(a)–(d) show the source images. We mixed them
Fig. 1. Comparison of Algorithm 1 with the existing methods: (a)–(d) source images, (e)–(h) mixed images (condition number of mixing matrix = 5.37), (i)–(l) layers separated by JADE (ISR = 0.77), (m)–(p) layers separated by Relative Newton (ISR = 0.19), (q)–(t) layers separated by our Algorithm 1 (ISR = 0.13).
K. Hara et al. / Pattern Recognition Letters 33 (2012) 559–567
563
Fig. 2. Examples of failed image separation under a permutation: (a)–(d) JADE, (e)–(h) Relative Newton, (i)–(l) our Algorithm 1.
R
Mixed image sequence
G
B
Three concatenated Color image sequence gray image sequence
Estimated source image sequence Fig. 3. Image separation through iterative image updates.
using different mixing ratios, as shown in Fig. 4(e)–(h). We show the results of our method outlined in Algorithm 2 in Fig. 4(i)–(l). One can see that our approach gives fairly good results. Fig. 5 shows the process of search and selection of the update permutation in this image separation. The vertical axis represents TðpÞ value (Eq. (10)). The horizontal axis represents the number of cycles executed for the mixture separation and selection of the update permutation. According to the criterion for accepting or rejecting update premutation candidates in Algorithm 2, six sequences {(3,1,4,2),(3,4,1,2),(4,3,2,1),(4,1,3,2),(1,3,4,2),(1,4,3,2)}, two sequences {(1,3,4,2), (1,4,3,2)}, and one sequence {(1,4,3,2)}
are selected as candidates after the first, second and third cycles, respectively, as shown by real lines in Fig. 5. Here the dotted lines represent the mixture separation processes for the rejected permutations. The labels 1, 2, 3 and 4 represent the images in Fig. 4(e), (f), (g) and (h), respectively. The circles and black crosses depict the update permutations which give perceptionally correct and incorrect separation results, respectively. We confirm that the finally selected permutation (1, 4, 3, 2) (drawn as a double circle) gives a perceptionally correct separation result, as shown at the right bottom of Fig. 5 (or Fig. 4(m)–(p)). The picture at the right top of Fig. 5 shows an unsuccessful result obtained for one
564
K. Hara et al. / Pattern Recognition Letters 33 (2012) 559–567
Fig. 4. Image separation with our Algorithm 2: (a)–(d) source images, (e)–(h) mixed images, (i)–(l) layers separated by our method.
failure
reject
reject reject success
Fig. 5. Process of mixture separation and selection of the update permutation (artificial mixtures).
of the rejected permutations. Table 2 compares the success ratio between Algorithm 1 and 2. In this experiment, we use a dataset, where each data element (each of No.1–No.10 in the top row in Table 2) contains four image mixtures synthesized using different mixing ratios. First we visually test Algorithm 1 for every permutation for each data element and then give the (number of successful permutations)/(total number of permutations) ratio in the second row in Table 2. Next we visually test Algorithm 2 for the same data element and then show success or failure of the image separation in the bottom row in Table 2. The success rates of Algorithm 1 and 2 are calculated from the above results and are shown in the right column in Table 2. Among the two algorithms, Algorithm 2 seems to achieve higher success ratio.
4.2. Real mixtures We apply our method to layer extraction from images containing multiple reflections and transparency. We photographed a picture postcard in a glass-fronted bookcase (Fig. 6(a)). One can see transparency and double reflections due to light reflected from both the surfaces of the front side and rear side glasses (region surrounded by a rectangle of Fig. 6(a)). We took three photographs under three different illumination conditions by inserting a polarization sheet at the front or back of two glasses on the bookcase, or between them (Fig. 6(c)–(e)). For acquisition of a ground truth transparency image we shot the same scene while shielding out some of the ambient light using a blackout curtain (Fig. 6(b)).
565
K. Hara et al. / Pattern Recognition Letters 33 (2012) 559–567 Table 2 Comparison between Algorithm 1 and 2. Methods
Algorithm 1
Input data number
number of successful permutations total number of permutations
Algorithm 2 (success: s; failure: )
Success ratio (%)
No. 1
No. 2
No. 3
No. 4
No. 5
No. 6
No. 7
No. 8
No. 9
No. 10
12/24
24/24
19/24
7/24
16/24
6/24
23/24
21/24
23/24
2/24
63.5
s
s
s
s
s
s
s
s
s
90.0
The separation results obtained by the JADE (Cardoso and Souloumiac, 1993), Relative Newton method (Gribonval and Zibulevsky, 2010) and our method are shown in Fig. 6(f)–(h), (i)–(k) and (l)–(n), respectively. One can see that, for both Relative Newton method and our method, the estimated transparency images (Fig. 6(i) and (l)) are similar to the ground truth transparency image (Fig. 6(b)), and that Fig. 6(j), (k), (m) and (n) are also relatively clear compared to JADE (Fig. 6(f)–(h)). In the same way as in the previous case, the process of search and selection of the permutation in this experiment is shown in Fig. 7. Here the labels 1, 2 and 3 represent the images in Fig. 6(c), (d) and (e), respectively. The update permutation (2, 3, 1) is selected and it gives a correct separation result, as shown at the center of Fig. 7 (or Fig. 6(l)–(n)). The picture at the right top of Fig. 7 shows an unsuccessful result obtained for one of the rejected permutations. 5. Conclusions In this paper, we have proposed a method of stably recovering a set of the original source images from the input mixtures. In many existing blind image separation techniques the success of the separation depends upon the selection of permutations corresponding to the orders of the update operations. By extending Sarel et al.’s bivariate generalized correlation measure (Sarel and Irani, 2004) to more than two variables and then introducing the simple pruning technique based on the discrepancy among the separation parameter values for each of the three RGB channels, we have presented a image separation method without permutation ambiguity. We have presented the separation algorithm, and shown its effectiveness through experiments with mixtures of real images. Acknowledgements This research was partially supported by the Suzuki Foundation. Appendix A. Derivation of Eq. (5) In this appendix, we derive Eq. (5) in detail. First we divide the correlation matrix RiK (Eq. (4)) into four blocks as,
ðA:1Þ where RiK1 is a K 1 dimensional square matrix and r iK1 is a K 1 dimensional vector. Hence, from the definition of the inverse matrix, ~riK;K included in Eq. (3) is expressed as
~riK;K Fig. 6. Separation of layers from images containing multiple reflections and transparency: (a) image region containing multiple reflections and transparency, (b) transparency image obtained by shielding out some of the ambient light using a blackout curtain, (c)–(e) three images taken under three different illumination conditions, (f)–(h) layers separated by JADE, (i)–(k) layers separated by Relative Newton, (l)–(n) layers separated by our method.
det RiK1 ; ¼ det RiK
ðA:2Þ
where det () denotes the matrix determinant. By first assuming detðRiK1 Þ – 0 and then applying the formula (Chen and Gu, 2000) for the determinant of a block matrix to Eq. (A.1), we have
566
K. Hara et al. / Pattern Recognition Letters 33 (2012) 559–567
failure
success
Fig. 7. Process of mixture separation and selection of the update permutation (real mixtures).
1 det RiK ¼ det RiK1 1 ðriK1 ÞT RiK1 r iK1 ;
ðA:3Þ
@ Wi @ rk
where ()1 denotes the inverse matrix. Using Eq. (A.3), we can rewrite Eq. (A.2) as,
~r iK;K ¼ 1
ðriK1 ÞT
1 : 1 RiK1 riK1
ðA:4Þ
Hence, using Eq. (A.4), we can rewrite Eq. (3) as,
T
T 1 C IV ;eI V C i IVK1 ;e IV i 1 U iK1 1=2 1=2 ½V i ðIV1 Þ ½V i ðIVK1 Þ h i K1 1=2 Q ¼2 V i IVl 0 0 V i IVk 00
ðA:5Þ n o IV1 ; . . . ; IVK1 ; eI V
MGNGC K IV1 ; . . . ; IVK1 ; eI V i2 PN h Q i V V V eV V i ðeI V Þ K1 i¼1 MNGC K I1 ; . . . ; IK1 ; I k¼1 V i Ik ¼ PN V eV QK1 i¼1 V i ð I Þ k¼1 V i Ik PN i T i 1 i QK1 V U K1 uK1 k¼1 V i Ik i¼1 uK1 ¼ ; Q PN V eV K1 i¼1 V i ð I Þ k¼1 V i Ik where
@uiK1 @ rk
l¼1
1=2 1 MNGC iK ðf1 ; . . . ; fK1 ; fK Þ ¼ ðriK1 ÞT RiK1 riK1 : Replacing {f1, . . . , fK1, fK} in Eqs. (2) and (A.5) with and then combining these equations as,
i @u@ W i K1 K1 Q C ðIV ;IV Þ C ðIV ;IV Þ ¼ 2 V i ðIVl Þ i V1 k1=2 i VK1 k1=2 ½V i ðI1 Þ ½V i ðIK1 Þ l¼1 T V V e C I ;I C IV ;e IV ðU iK1 Þ1 i 1V 1=2 i K1 1=2 ½V i ðI1 Þ ½V i ðIVK1 Þ h i h i K1 1=2 Q ¼2 V i IVl V i IVk ui1;k uiK1;k
¼
ðB:2Þ
l¼1
T V C I ;e IV C IV ;e IV i 1V 1=2 i K1 1=2 ½V i ðI1 Þ ½V i ðIVK1 Þ K1 Q V ¼ 2C i eI V ; IVk V i Il l¼1
K1 K1 P Q V ¼ 2 C i IVK ; IVk rm C i IVm ; IVk V i Il : m¼1
ðA:6Þ
U iK1
is the K 1 dimensional square correlation matrix C i ðIVk ;IVl Þ i ¼ whose klth entry is U iK1 1=2 and uK1 is the K 1 kl ½V i ðIVk ÞV i ðIVl Þ C IV ;e IV dimensional vector whose kth entry is uiK1 k ¼ i kV 1=2 , as ½V i ðIk Þ described in Section 2.3.
l¼1
P So, the partial derivatives of the numerator, Ni¼1 Wi , of Eq. (5) with respect to rk, for k = 1, . . . , K 1, are given by N @ X
@ rk
i¼1
Wi ¼ 2
N X
"
# K1 X K1 Y V C i IVk ; IVK rm C i IVk ; IVm V i Il
i¼1
m¼1
l¼1
ðk ¼ 1; . . . ; K 1Þ:
ðB:3Þ
Since one can easily see that for k = 1, . . . , K 1, the partial deriva PN V eV QK1 i¼1 V i ð I Þ k¼1 V i I k , of Eq. (5) with re-
tives of the denominator,
spect to rk is also equal to Eq. (B.3), we have
" # N K 1 1 X KY X @F @G ¼ ¼2 C i IVk ; IVK rm C i IVk ; IVm V i IVl @ rk @ rk m¼1 i¼1 l¼1
Appendix B. Derivation of Eq. (6) Denoting the numerator of Eq. (5) by
Wi ¼ uiK1
T
U iK1
1
uiK1
K 1 Y
V i IVk
PN
i¼1
Wi as,
ðk ¼ 1; . . . ; K 1Þ; ðB:1Þ
k¼1
and then taking partial derivatives of Wi with respect to rk, for k = 1, . . . , K 1, we get
ðB:4Þ
where F and G are the denominator and the numerator of Eq. (5), respectively. Thus, taking partial derivatives of (5) with respect to rk and setting them to zero, we get
K. Hara et al. / Pattern Recognition Letters 33 (2012) 559–567
@ @ F 1 @F @G 0¼ MGNGC K IV1 ; . . . ; IVK1 ; eI V ¼ F ¼ 2 G @ rk @ rk G @ rk @ rk G 1 F @F ðk ¼ 1; . . . ; K 1Þ: ðB:5Þ 1 ¼ G G @ rk Then, from the definition of MGNGC and the independency of fIVk gKk¼1 , we have F/G < 1 Hence, we get @ F/@ rk = 0 in (B.5), which leads to the following equation.
" # N K1 X K1 X Y V @F ¼2 C i IVk ; IVK rm C i IVk ; IVm V i Il @ rk m¼1 i¼1 l¼1 ðk ¼ 1; . . . ; K 1Þ:
ðB:6Þ
This equation can be rewritten as,
" K 1 X N X m¼1
¼
# 1 KY V V V C i Ik ; Im V i Il rVm
i¼1 N X
l¼1
C i IVk ; IVK
i¼1
K1 Y
V i IVl ðk ¼ 1; . . . ; K 1Þ:
ðB:7Þ
l¼1
References Bedini, L., Herranz, D., Salerno, E., Baccigalupi, C., Kuruoglu, E., Tonazzini, A., 2005. Separation of correlated astrophysical sources using multiple-lag data covariance matrices. EURASIP J. Appl. Sig. Proc. 15, 2400–2412. Bell, A.J., Sejnowski, T.J., 1995. An information maximization approach to blind separation and blind deconvolution. Neural Comput. 7 (6), 1129–1159. Bronstein, A.M., Bronstein, M.M., Zibulevsky, M., Zeevi, Y.Y., 2005. Sparse ICA for blind separation of transmitted and reflected images. Internat. J. Imag. Sci. Technol. 15 (1), 84–91. Cardoso, J.F., Souloumiac, A., 1993. Blind beamforming for non gaussian signals. In: IEE Proceedings-F, 140(6), 362–370.
567
Cardoso, J.F., 1994. On the performance of orthogonal source separation algorithms. In: Proc. EUSIPCO, pp. 776–779. Cardoso, J.F., Snoussi, H.H., Delabrouille, J., Patanchon, G., 2002. Blind separation of noisy gaussian stationary sources: Application to cosmic microwave imaging. In: Proc. EUSIPCO, vol. 1, pp. 561–564. Chen, J., Gu, G., 2000. Control-Oriented System Identification: An H-Infinity Approach. John Wiley. Cichocki, A., Amari, S., 2002. Adaptive Blind Signal and Image Processing. John Wiley. Diamantaras, K.I., Papadimitriou, T., 2005. Blind separation of reflections using the image mixtures ratio. In: Proc. IEEE ICIP, pp. 1034–1037. Farid, H., Adelson, E.H., 1999. Separating reflections from images by use of independent components analysis. J. Opt. Soc. America A 16 (9), 2136–2145. Gribonval, R., Zibulevsky, M., 2010. Sparse component analysis. In: Comon, Pierre, Jutten, Christian (Eds.), Handbook of Blind Source Separation: Independent Component Analysis and Applications. ELSEVIER. Hara, K., Inoue, K., Urahama, K., 2009. Separation of layers from images containing multiple reflections and transparency using cyclic permutation. In: Proc. IEEE ICASSP, pp. 1157–1160. Hyvaerinen, A., Oja, E., 1997. A fast fixed-point algorithm for independent component analysis. Neural Comput. 9 (7), 1483–1492. Kayabol, K., Kuruoglu, E.E., Sankur, B., 2009. Image source separation using color channel dependencies. In: Proc. ICA, pp. 499–506. Kopriva, L., 2007. Approach to blind image deconvolution by multiscale subband decomposition and independent component analysis. J. Opt. Soc. America A 24 (4), 973–983. Levin, A., Weiss, Y., 2007. User assisted separation of reflections from a single image using a sparsity prior. IEEE Trans. Pattern. Anal. Machine. Intell. 29 (9), 1647– 1654. Levin, A., Zomet, A., Weiss, Y., 2004. Separating reflections from a single image using local features. In: Proc. IEEE CVPR, pp. 306–313. Lin, S., Shum, H.Y., 2001.Separation of diffuse and specular reflection in color images. In: Proc. IEEE CVPR, pp. 341–346. Nayar, S.K., Fang, X.S., Boult, T., 1996. Separation of reflection components using color and polarization. Internat J. Comput. Vision. 21 (3), 163–186. Sarel, B., Irani, M., 2004. Separating transparent layers through layer information exchange. In: Proc. ECCV, pp. 328–341. Tonazzini, A., Bedini, L., Salerno, E., 2006. A markov model for blind image separation by a mean-field EM algorithm. IEEE Trans. IP. 15 (2), 473–482.