Space-Variant Gabor Decomposition for Filtering 3D Medical Images

Space-Variant Gabor Decomposition for Filtering 3D Medical Images Darian Onchis1,2(B), Codruta Istin3, and Pedro Real4 1University of Vienna, Vienna, Austria [email protected] 2West University of Timisoara, Timisoara, Romania 3Politehnica University of Timisoara, Timisoara, Romania 4University of Seville, Seville, Spain Abstract. This is an experimental paper in which we introduce the possibility to analyze and to synthesize 3D medical images by using multivariate Gabor frames with Gaussian windows. Our purpose is to apply a space-variant filter-like operation in the space-frequency domain to correct medical images corrupted by different types of acquisitions errors. The Gabor frames are constructed with Gaussian windows sampled on non-separable lattices for a better packing of the space-frequency plane. An implementable solution for 3D-Gabor frames with non-separable lattice is given and numerical tests on simulated data are presented. 1 Introduction The study of noise removal in medical images was approached in different ways by numerous authors. From the first studies involving only convolutional filters [12], to the rank algorithms [13,14] and most recently to methods to average pixels depending on their neighborhood statistics [10,11], the search was both for maximizing the filter capabilities and to propose fast algorithms. While very successful for various signal analysis applications in medicine or telecommunications, time-frequency analysis and especially the Gabor frames expansions for image processing in 2D or 3D was of limited interest due to the large indexing problems i.e. a 3D image is analysed with a 6D lattice, and to the problems related to time consuming implementation. In a fairly recent paper [15], the first author introduced a procedure to efficiently perform an nD Gabor frames decomposition. Based on that result, we propose in this paper a method to decompose, to filter and to reconstruct 3D medical images that overcomes both problems. We are using a tensor product decompositions [4] to reduce the decomposition to only 2D lattice case and also fast algorithms to implement this case. The result is an experimental framework for fast Gabor frames construction used to analyze 3D-data images sampled on quincunx-type lattices. The research is in the simulation stage and we present in here only tests done on the MATLAB MRI data set. Therefore the main contribution of this conference paper is to introduce a completely invertible 3D Gabor transform applied to filter 3D medical images, c Springer International Publishing AG 2017 M. Felsberg et al. (Eds.): CAIP 2017, Part II, LNCS 10425, pp. 455–461, 2017. DOI: 10.1007/978-3-319-64698-5 38 456 D. Onchis et al. different from existing examples in literature that present only 3D Gabor filters [7,9]. This paper is structured as follows: In the second section we recall the necessary results from Gabor analysis, while in the third section we present the 3D Gabor frames construction. The application to 3D medical imaging is given in the fourth section. In the last section, we present timings for the 3D Gabor construction and the conclusions are drawn. 2 Theoretical Preliminaries Frames (gi)i∈Igeneralize the idea of a basis in a Hilbert space Hand they are formally defined as: Definition 1. Afamily(gi)i∈Iin a Hilbert space His called a frame if there exist constants A, B > 0such that for all g∈H Ag2≤ i∈I |g,gi|2≤Bg2(1) Every element f∈Hhas an expansion of the form: f=SS−1f= i∈I S−1f,gigi= i∈I f,S−1gigi where Sdenotes the invertible frame operator [3]: Sf =i∈If,gigi.The family (γi)i∈I=(S−1gi)i∈Iis again a frame with frame bounds B−1and A−1 and is called a canonical dual frame. The main tool for time-frequency analysis is the Short-Time Fourier Transform in short STFT, defined for functions in L2(R)as Vgf(λ)=Vgf(a, b)=f,MbTag=f,π(λ)g(2) where Taf(t)=f(t−a) is the translation (time shift) and Mbf(t)=e2πib·tf(t) is the modulation (frequency shift), for λ=(a, b)inR2. The operators MbTaare called time-frequency shifts. Their composition is denoted by π(λ). In order to obtain Gabor frames, the STFT is sampled over a time-frequency lattice. In the standard 1D case the regular lattice is of the form Λ=aZ×bZwith the condition ab < 1. Therefore a Gabor frame is defined as: G(g,a,b):={gk, l =MblTakg, k,l ∈Z},(3) i.e. the elements of the Gabor frame are translated and modulated versions of one atom g. Space-Variant Gabor Decomposition for Filtering 3D Medical Images 457 Theorem 1 (Dual Gabor Frames [5]). If G(g,a,b)is a frame for L2(R), then the canonical dual frame takes the form G(γ,a,b)for γ=S−1g. where the Gabor frame operator S is defined as: S:=  k, l∈Z f,MblTakgMblTakg(4) We introduce also the result of Bourouihiya [2], which extends the classical results of Lyubarski and Seip to higher dimensions [6]forg0(x)=2 −1/4e−πx2 with x∈R: Lemma 1. Let g=g0⊗···⊗g0(nfactors) and Λab =(a1Z×···×anZ)×(b1Z× ···×bnZ).ThenG(g,Λab)is a frame if and only if ajbj<1for 1≤j≤n. Based on this result, another extension to non-separable lattices of the form NΛa,b is possible, which besides a multi-variate sampling in each dimension, give us a better packing of the time-frequency plane by using a (non-separable) hexagonal lattice that match with the circular contour lines of the Gaussian. The representation of the non-separable lattice is based on the rectangular lattice via a shear operation. Therefore, we can give the following lemma: Lemma 2. Given a window g=g0⊗···⊗g0(nfactors) and a lattice of the form NΛa,b,whereΛis a rectangular lattice and Nis a shear lattice, the system G(g,Λab)is a frame if and only if ajbj<1for 1≤j≤n. Proof. The shear matrix N is a symplectic matrix hence the determinant det(N) = 1. Therefore the volumes vol(NΛa,b)=vol(Λa,b) are equal and the Lemma 1extends to non-separable lattices. In applications sampled data of finite lengths are analysed; the process of sampling and periodization are employed [16]. In this way also the number of shifts in time and frequency becomes finite. The redundancy of a discrete system, not necessarily a Gabor system, is defined as the fraction of the number of used discrete function over the length of the domain, #shift L. 3 3D Gabor Frame Construction In this section we present how to construct numerically a 3D Gabor transform for 3D image filtering. We will exploit the possibilities given by the Lemma 2for choosing a non-separable lattice of quincunx-type. This situation can be easily expanded to more dimensions or can be reduced to less dimensions following the same procedure. We write n for data cube length in one direction and the modulation and translation operator are defined on Zn. For further discretization details [8] is a comprehensive source. 458 D. Onchis et al. For the 3D case, we will use a multi-variate generalized Gaussian window, obtained as a tensor product of 1D windows in the form: G3=g1⊗g2⊗g3 We consider our 3D Gaussian window of size n1×n2×n3as a complexvalued function on the additive Abelian group G=n1×n2×n3. The joint position-frequency space is G× G=n1×n2×n3× n1×n2×n3. Now, let’s pay attention to the lattice. We would like to use a non-separable lattice obtained as in the hypothesis of the Lemma 2by applying a shear matrix to the matrix generating the regular standard lattice. There are other symplectic matrices like the rotation matrix using the fractional Fourier transform [1] that can be used to transform a rectangular lattice into a quincunx-like lattice. For our 3D case using a non-separable lattice, we will consider the following quincunx-like lattice: Λ3=I3×3Q3 O3I3×3A6·Z6 where Q3=⎛ ⎝ a1/2b100 0a2/2b20 00a3/2b3 ⎞ ⎠ and A6=diag(a1,a 2,a 3,b 1,b 2,b 3) This 6D non-separable lattice can be reduced to the case of 2D nonseparability in time and frequency by the expansion in the canonical basis. Moreover, we can write the following tensor product: 1a1/2b1 01 ⊗⎛ ⎝ 100 000 000 ⎞ ⎠+1a2/2b2 01 ⊗⎛ ⎝ 000 010 000 ⎞ ⎠+1a3/2b3 01 ⊗⎛ ⎝ 000 000 001 ⎞ ⎠ In this conditions it is feasible to perform a Gabor analysis followed by a Gabor synthesis at level of the constituting vectors, and therefore reducing the case of applying 3D matrices to the faster case of applying the transform over each dimension. Therefore, we can compute the coefficients of the expansion (e.g. the 3D Gabor filter) by using any fast algorithms developed for the 1D Gabor transform [8]. For a complete Gabor transform, we use the Theorem 1and we invert the frame operator in order to obtain the dual frame. Using the computed dual frame, we can recover the 3D data cube. Space-Variant Gabor Decomposition for Filtering 3D Medical Images 459 4 Applications to 3D Medical Image Processing We have tested our implementation for the MRI data set that comes with MATLAB2017. Loading mri.mat adds two variables to the workspace: D (128-by-128- by-1-by-27, class uint 8) and a grayscale colormap, map (89-by-3, class double). D comprises 27 128-by-128 horizontal slices from an MRI data scan of a human cranium. Values in D range from 0 through 88, so the colormap is needed to generate a figure with a useful visual range. The procedure takes as inputs D = 3D medical image and the G = the Gaussian analysing atoms. The parameters a = time shift and b = frequency shift are defined depending of the application and under the liniar independence contraints. The output is DG = data cube recovered after Gabor analysis. The complete algorithm is summarized below: Algorithm 1. 3D Gabor transform filtering with quincunx lattice Input: G- Gaussian atom, D- 3D medical image Output: DG - 3D Gabor filtered medical image 1: Generate Gabor matrix Gfor a quincux sampling lattice; 2: Compute the Gabor coefficients GC over each dimension 3: Remove the low energy Gabor coefficients 4: Recover DG over each dimension with the use of the dual atom In the Figs. 1and 2, we present the results of the applying the Gabor filtering. Fig. 1. Unfiltered 3D image We observed that due to the reconstruction after the removal of the low energy coefficients the image becomes much smoother. The next step is to corroborate this results with the needs of medical practitioners or radiologists. 460 D. Onchis et al. Fig. 2. Gabor filtered 3D image 5 Conclusions In this paper, we described the construction of Gabor frames beyond the standard Gabor processing on regular lattices for signal and image processing. In this way, one has the freedom to choose a non-separable lattice and to sample the nD-Gaussian window in a multi-variate way, while still being assured that the result will be a Gabor frame. This added flexibility will be useful for application like 3D space-variant filtering or decomposition in the 3D Gabor domain for different features identifications (e.g. 3D plane waves). The matrix we used for generating the lattice in our example, verifies the conditions under which the hypothesis of the Lemma 2are true (i.e. shear matrix). The lattice parameters aj,b jshould be chosen according to the time-frequency concentration of the corresponding windows gjto obtain well concentrated systems. Therefore, under the conditions ajbj<1forj=1,2,3, we obtain Gabor frames for decompositions and analysis. The timings for the implementation in MatlabTM on a notebook with Intel core i5 at 2.3 GHz and with 8.00 GB of RAM memory are given in Table 1. Table 1. Timings comparison. The case 1. With redundancy 2, ab =1 2thecase2., with redundancy 4, ab =1 4. Timings for 3D MRI data 1. n = 128 2. n = 128 Quincunx lattice Coefficients 12.438 s 72.103 s Synthesis 13.134 s 78.363 s Acknowledgments. The first author gratefully acknowledge the support of the Austrian Science Fund (FWF): project number P27516. Space-Variant Gabor Decomposition for Filtering 3D Medical Images 461 References 1. Bastiaans, M.J., van Leest, A.J.: From the rectangular to the quincunx Gabor lattice via fractional Fourier transformation. IEEE Signal Process. Lett. 5(8), 203– 205 (1998) 2. Bourouihiya, A.: The tensor product of frames. Sampl. Theory Signal Image Process. 7(1), 65–76 (2008) 3. Christensen, O.: An Introduction to Frames and Riesz Bases. Applied and Numerical Harmonic Analysis. Birkh¨auser, Boston (2003) 4. Christensen, O., Feichtinger, H., Paukner, S.: Gabor Analysis for Imaging. Handbook of Mathematical Methods in Imaging. Springer, Berlin (2010) 5. Gr¨ochenig, K.: Foundations of Time-Frequency Analysis. Birkh¨auser, Boston (2001) 6. Lyubarskii, Y.I.: Frames in the Bargmann space of entire functions. In: Entire and Subharmonic Functions, pp. 167–180. American Mathematical Society, Providence (1992) 7. Mikula, K., Sgallari, F.: Semi-implicit finite volume scheme for image processing in 3D cylindrical geometry. J. Comput. Appl. Math. 161(1), 119–132 (2003) 8. Qiu, S., Feichtinger, H.G.: Discrete Gabor structures and optimal representation. IEEE Trans. Signal Process. 43(10), 2258–2268 (1995) 9. Wang, Y., Chua, C.-S.: Face recognition from 2D and 3D images using 3D Gabor filters. Image Vis. Comput. 23(11), 1018–1028 (2005) 10. Buades, A., Morel, J.M.: A Non-local algorithm for image denoising. In: IEEE Computer Society Conference on Computer Vision and Pattern Recognition, vol. 2, pp. 60–65, 20–26 June 2005 11. Dabov, K., Foi, A., Katkovnik, V., Egiazarian, K.: Image denoising by sparse 3D transform-domain collaborative filtering. IEEE Trans. Image Process. 16(8), 2080– 2095 (2007) 12. Pratt, W.K.: Digital Image Processing: PIKS Scientific inside, 4th edn. Wiley, Los Altos (2007) 13. Yaroslavsky, L.P., Kim, V.: Rank algorithms for picture processing. Comput. Vis. Graph. Image Process. 35, 234–258 (1986) 14. Storozhilova, M., Lukin, A., Yurin, D., Sinitsyn, V.: 2.5D extension of neighborhood filters for noise reduction in 3D medical CT images. In: Gavrilova, M.L., Tan, C.J.K., Konushin, A. (eds.) Transactions on Computational Science XIX. LNCS, vol. 7870, pp. 1–16. Springer, Heidelberg (2013). doi:10.1007/978-3-642-39759-2 1 15. Onchis, D.M.: Optimized frames and multi-dimensional challenges in timefrequency analysis. Adv. Comput. Math. 40(3), 703–709 (2014) 16. Søndergaard, P.L.: Gabor frames by sampling and periodization. Adv. Comput. Math. 27(4), 355–373 (2007)