Multifocus Image Fusion and Restoration---李帅.pdf

884IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 59, NO. 4, APRIL 2010 Multifocus Image Fusion and Restoration With Sparse Representation Bin Yang and Shutao Li, Member, IEEE AbstractTo obtain an image with every object in focus, we al- ways need to fuse images taken from the same view point with dif- ferent focal settings. Multiresolution transforms, such as pyramid decomposition and wavelet, are usually used to solve this problem. In this paper, a sparse representation-based multifocus image fusion method is proposed. In the method, fi rst, the source image is represented with sparse coeffi cients using an overcomplete dictio- nary. Second, the coeffi cients are combined with the choose-max fusion rule. Finally, the fused image is reconstructed from the combined sparse coeffi cients and the dictionary. Furthermore, the proposed fusion scheme can simultaneously resolve the image restoration and fusion problem by changing the approximate criterion in the sparse representation algorithm. The proposed method is compared with spatial gradient (SG)-, morphological wavelet transform (MWT)-, discrete wavelet transform (DWT)-, stationary wavelet transform (SWT)-, curvelet transform (CVT)-, andnonsubsamplingcontourlettransform(NSCT)-basedmethods on several pairs of multifocus images. The experimental results demonstrate that the proposed approach performs better in both subjective and objective qualities. IndexTermsImagefusion,imagerestoration,sparse representation. I. INTRODUCTION N OWADAYS, image fusion has become an important sub- area of image processing. For one object or scene, mul- tiple images can be taken from one or multiple sensors. These images usually contain complementary information. Image fu- sion is the process of detecting salient features in the source images and fusing these details to a synthetic image. Through image fusion, extended or enhanced information content can be obtained in the composite image, which has many application fi elds, such as digital imaging, medical imaging, remote sens- ing, and machine vision 15. As an example of fusion that is relevant to this paper, optica imaging cameras suffer from the problem of fi nite depth of fi eld, which cannot make objects at various distances (from the sensor) all in focus. Therefore, if one object in the scene is in focus, then the other objects at different distances from the camera will be out of focus and, thus, blurred. The solution ManuscriptreceivedJanuary17,2009;revisedMay22,2009.Firstpublished October 30, 2009; current version published March 20, 2010. This work was supported in part by the National Natural Science Foundation of China under Grants 60871096 and 60835004, by the Ph.D. Programs Foundation of the Ministry of Education of China under Grant 200805320006, and by the Key Project of the Chinese Ministry of Education under Grant 2009- 120. The Associate Editor coordinating the review process for this paper was Dr. Cesare Alippi. The authors are with the College of Electrical and Information Engineering, Hunan University, Changsha 410082, China (e-mail: yangbin01420; shutao_li). Digital Object Identifi er 10.1109/TIM.2009.2026612 to get all the objects focused in one image is multifocus image fusion technique. In this technique, several images of a scene are captured with focus on different parts. Then, these images are fused with the hope that all the objects will be in focus in the resulting image 59. There are various methods available to implement image fusion. Basically, these methods can be categorized into two categories. The fi rst category is the spatial domain-based meth- ods, which directly fuse the source images into the intensity values 1013. The other category is the transformed domain-based methods, which fuse images with certain fre- quency or timefrequency transforms 1, 4, 5. Assuming that F() represents the “fusion operator,” the fusion methods in the spatial domain can be summarized as IF= F(I1,I2,.,IK).(1) The simplest fusion method in spatial domain just takes the pixel-by-pixel average of the source images. However, this method often leads to undesirable side effects, such as reduced contrast 1. If the source images are not completely registered, then a single pixel-based method, such as spatial gradient (SG)- based method 10, always results in artifacts in the fused image. Therefore, some more reasonable methods were pro- posed to fuse source images with divided blocks or segmented regions instead of single pixels 1113. However, the block- based fusion methods usually suffer from blockness in the fused image 11. For the region-based method, the source images are fi rst segmented, and the obtained regions are then fused using their properties, such as spatial frequency or SG. The segmentation algorithms, usually complicated and time consuming, are of vital importance to the fusion quality 13. A more popular method that has been explored in recent years is by using multiscale transforms. The usually used multiscale transforms include various pyramids 1417, dis- crete wavelet transform (DWT) 1, 5, 18, complex wavelet 19, 20, ridgelet 21, curvelet transform (CVT) 22, and contourlet 23. The transformed domain-based methods can be summarized as IF= T1(F (T(I1),T(I2),.,T(IK)(2) where T() represents a multiscale transform, and F() means the applied fusion operator. Pyramid decomposition is the earliest multiscale transform used for image fusion 1417. In this method, each source image is fi rst decomposed into a sequence of images (pyra- mid) in different resolutions. Then, at each position in the transformed image, the value in the pyramid with the highest 0018-9456/$26.00 2009 IEEE YANG AND LI: MULTIFOCUS IMAGE FUSION AND RESTORATION WITH SPARSE REPRESENTATION885 saliency is selected. Finally, the fused image is constructed us- ing the inverse transform of the composite images. The wavelet transform-based fusion methods employ a similarscheme to the pyramid transform-based methods. However, the performance of multiresolution transform-based methods is limited owing to that most of the multiresolution decompositions are not shift invariant, which is brought by the underlying down-sampling process 1, 5, 19. The shift-invariant extension of the DWT can yield an overcomplete signal representation, which is suitable for image fusion 2426. Further, many advanced geometric multiscale transforms, such as CVT, ridgelet, and contourlet, have been explored in recent years and shown improved results 2123, 2729. However, because the fused image obtained by transform domain-based algorithms is globally created, a little change in a single coeffi cient of the fused image in the transformed domain may cause all the pixel values to change in spatial domain. As a result, undesirable artifacts may be produced in the fusion process using the multiresolution transform-based methods in some cases 3. Obviously, effectively and completely extracting the under- lying information of the original images would make the fused image more accurate. Different from multiscale transforma- tions, the sparse representation using an overcomplete dictio- nary that contains prototype signal atoms describes signals by sparse linear combinations of these atoms 3034. Two main characteristicsofsparserepresentationareitsovercompleteness and sparsity 32. Overcompleteness means that the number of basis atoms in the dictionary exceeds the number of im- age pixels or signal dimensions. The overcomplete dictionary that contains rich transform bases allows for more stable and meaningful representation of signals. Sparsity means that the coeffi cients corresponding to a signal are sparse, that is to say, only “a few descriptions” can describe or capture the signifi cant structure information about the object of interest. Benefi ting from its sparsity and overcompleteness, sparse representation theory has successfully been applied in many practical appli- cations, including compression, denoising, feature extraction, classifi cation, and so on 3234. Recent studies have shown that common image features can also be accurately described by only a few coeffi cients or “a few descriptions” 32. Using the few coeffi cients as the salient features of images, we design a sparse representation (SR)-based image fusion scheme. In general, sparse representation is a global operation, in the sense that it is based on the gray-level content of an entire image. However, the image fusion quality depends on the accurate representation of the local salient features of source images. Therefore, a “sliding window” technique is adopted to achieve better performance in capturing local salient features and keep- ing shift invariance. In the proposed method, the source images are fi rst divided into patches, which lead to “a small size” of overcomplete dictionary to every patch. Second, the patches are decomposed by prototype signal atoms into their corresponding coeffi cients. The larger the coeffi cient is, the more salient features it contains. Third, the “choose-max” fusion rule is used to combine the corresponding coeffi cients of the source images. Finally, the result image is reconstructed using the combined coeffi cients previously obtained. In sparse representation, the dictionary (atoms) is often created by a prespecifi ed set of functions, such as discrete cosine transforms (DCT), short-time Fourier transforms, wavelet, CVT, and contourlet. More complex dictionary can be obtained by learning from a given set of signal examples 32. Note that most of the image fusion methods are based on the assumption that the source images are noise free. Therefore, these fusion algorithms can produce high-quality fused images if the assumption is satisfi ed. However, practically, the images are often corrupted by noise during acquisition or transmission processes. For multiresolution-based methods, they usually de- noise the source images, fi rst, by setting all the coeffi cients below a certain threshold to zero and keeping the remaining coeffi cientsunchanged.Then,thefi lteredimagesarefused.One advantage of our proposed method is that it can simultaneously carry out denoising and fusion of noisy source images. The rest of this paper is organized into fi ve sections. In Section II, the basic theory of sparse representation is pre- sented.InSectionIII,wepropose thefusionschemewithsparse representation theory and discuss how to simultaneously carry out image restoration and fusion. Numerical experiments and discussions are detailed in Section IV. Both advantages and disadvantages of the proposed schemes, together with some suggestions about the future work, are given in Section V. II. BASICTHEORY OFSIGNALSPARSEREPRESENTATION Sparse representation is based on the assumption that natural signals can be represented or approximately represented as a linear combination of a “few” atoms from dictionary 30. For signals Rn, sparse representation theory suggests the existence of a dictionary D RnT, which contains T pro- totype signals that are referred to as atoms. For any signal x , there exists a linear combination of atoms from D that approximates it well. More formally, it is x , s RT such that x Ds. The vector s contains “coeffi cients” of x in D. It usually assumes that T n, implying that the dictionary D is redundant. The solution is generally not unique. Finding the smallest possible number of nonzero components of s involves solving the following optimization problem: min s ?s?0subject to ?Ds x? (3) where ?s?0denotes the number of nonzero components in s. The above optimization is an NP-hard problem and can only be solved by systematically testing all the potential combina- tions of columns 31. Thus, approximate solutions are con- sidered instead. In the past decade, several pursuit algorithms have been proposed to solve this problem. The simplest algo- rithms are the matching pursuit (MP) 30 and the orthogonal MP (OMP) algorithms 32. They are greedy algorithms that sequentially select the dictionary atoms. These methods involve the computation of the inner product between the signal and dictionary columns. The MP algorithm aims to learn x Ds as follows: Let ?,? denote the inner product. Initialize the residual function r0as r0= x.(4) 886IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 59, NO. 4, APRIL 2010 Then, loop over all prototype signal atoms and select the index t l of the best atom function in D as t l = argmax t ?rl,d t ? (5) with a coeffi cient along that dimension as s(t l) = ? rl,dt l ? .(6) Then, update the residual function rl+1= rl s(t l)dt l (7) and repeat until ?rl?2= ? ? ? ? ?x l ? i=1 s(t i)dt i ? ? ? ? ? 2 (8) or the iteration reaches the priori set number. Because the computation of the MP algorithm is intensive, in this paper, the OMP is used. The OMP method updates the residual as rl+1= x Pspan ? dt 1,dt 2,.,dt l ? x(9) where P denotes the orthogonal projection onto a subspace. The OMP algorithm, in parallel, applies the GramSchmidt orthogonalization upon the chosen atoms for the effi cient com- putationofprojections.Ideally,thenumberofiterationsisequal to the number of nonzeros in s. Other well-known pursuit approaches include the basis pursuit (BP) 35 and the focal un- derdetermined system solver (FOCUSS) 33. The BP replaces the ?0-norm with an ?1-norm in (3), whereas the FOCUSS uses the ?p-norm with as a replacement for the ?0-norm. III. IMAGEFUSIONMETHOD A. Sparse Representation for Image Fusion Since the sparse representation globally handles an image, it cannot directly be used with image fusion, which depends on the local information of source images. In our method, we divide the source images into small patches and use the fi xed dictionary D with small size to solve this problem. In addition, a sliding window technique is adopted to make the sparse representation shift invariant, which is of great importance to image fusion. We assume that source image I is divided into many image patches. As shown in Fig. 1, to facilitate the analysis, the jth patch with size n n is lexicographically ordered as a vector vj. Then, vjcan be expressed as vj= T ? t=1 sj(t)dt(10) where dtis an atom from a given overcomplete dictionary, and D = d1 dt dT, which contains T atoms. sj= sj(1) sj(t) sj(T)Tis the sparse representation obtained in (3). Fig. 1.Selected image patch and its lexicographic ordering vector. Fig. 2.Schematic diagram of the proposed SR-based fusion method. Assume that the vectors responding to all the patches in image I are constituted into one matrix V. Then, V can be expressed as V = d1d2dT s1(1)s2(1)sJ(1) s1(2)s2(2)sJ(2) . . . . . . . . . . s1(T)s2(T)sJ(T) (11) whereJisthenumberofimagepatches.LetS = s1,s2,sJ. Then, (11) can be expressed as V = DS(12) where S is a sparse matrix. B. Proposed Fusion Scheme Assume that there are K registered source images I1,.,IK with size of M N. Then, the proposed fusion scheme based on image sparse representation, shown in Fig. 2, takes the following steps. Use the sliding window technique to divide each source image Ik, from left-top to right-bottom, into patches of size n n, i.e., the size of the atom in the dictionary. Then, all the patches are transformed into vectors via lexicographic ordering, and all the vectors constitute one matrix Vk, in which each column corresponds to one patch in the source image Ik. The size of Vkis (n n) (M n + 1) (N n + 1). ForthejthcolumnvectorvkjinVk,itssparserepresentation is calculated using the OMP method. The OMP iterations will stop when the representation error drops below the specifi ed tolerance. Then, we get a very sparse representation vector skj for vkj. YANG AND LI: MULTIFOCUS IMAGE FUSION AND RESTORATION WITH SPARSE REPRESENTATION887 Then, the activity level of skjresponding to the kth source image Ikis obtained as Akj= ?skj?1.(13) Fuse the corresponding columns of all sparse representation matrix S1,.,Sk,.,SKof the source images to generate SFaccording to their activity levels. The jth column of SF is obtained as sFj= sk j, k j = argmax kj (Akj).(14) The vector representation of the fused image VFcan be calculated by VF= DSF.(15) Finally, the fused image IFis reconstructed using VF. Reshape each vector vFjin VFinto a block with size n n and then add the block to IFat its responding position. This can be seen as the inverse process of Fig. 1. Thus, for each pixel position, the pixel value is the sum of several block values. Then, the pixel value is divided by the adding times at its position to obtain the reconstructed result. C. Restoration and Fusion In practice, the source images for fusion are often corrupted by noise during the acquisition or transmission process. In traditional methods, image restoration and image fusion are separately treated. Little effort has been made to combine them 36. The other benefi t of the proposed scheme is to simultaneously conduct image restoration and fusion by using the advantage of sparse representation in image restoration. Assume that image I is contaminated by additive zero-mean white Gaussian noise with standard deviation . For the jth patchs vector of image I vj, its denoising solution using the maximum a posteriori estimator is sj= min sj ? ?sj? 0 subject to ? ?vj Dsj? 2 C(16) whereCisaconstant34.Then,allthepatches sparse representation vectors constitute the sparse coeffi cient matrixS. Then, for source images I1,.,IK, the corresponding over- lapped patches and vectoring matrix are V1,V2,.,VK. Their sparse representations S1,S2,.,SKare obtained ac- cording to (16). The restored image matrices are calculated by V1= DS1, V2= DS2,.,VK= DSK. The