Blind image separation using pyramid technique

Signal and image separation is an important processing step for accurate image reconstruction, which is increasingly applied to many medical imaging applications and communication systems. Most of the conventional separation approaches are based on frequency domain and time domain. These approaches, however, are sensitive to noise and thus often produce undesirable results.

In this paper, we propose a novel method of image separation. It incorporates the property of pyramid component extracted from the image and a finite ridgelet transform (FRT) to obtain a precise analysis of the images and thus correctly separate the images even in a highly noisy environment. We obtain the multiple components of the target images by employing a pyramid processing, which operates in the various domains and thus can decompose the image into multiple components.

In addition, the pyramid decomposition in the proposed method can eliminate information redundancy in the target image and thus can substantially enhance the quality of image separation. We have conducted extensive simulations, which demonstrate that the proposed pyramid structure with FRT outperforms the conventional methods based on time domain and trigonometric transforms.

Blind source separation (BSS) has been one of the major research areas for over a decade and is receiving growing attention due to its processing applications in image and signal processing. It aims at extracting a set of source signals from an observed mixture of signals with little or no information about either the mixing environment or the mixing process and sources. The applications of BSS range from medical engineering to neuroscience and also from telecommunications to financial time series analysis. For example, its recent applications include astronomical imaging, remote sensing, medical imaging, biological data analysis, image and speech signal processing, etc. [1,2,3].

Independent component analysis (ICA) has been often regarded as an attractive solution to the BSS problem. Its process is based on non-Gaussianity method, and so, it can utilize its statistical independence of the sources to calculate the de-mixing matrix and extract the source signals with a scaling factor and permutation [4, 5].

In biomedical applications, ICA has been applied to the functional magnetic resonance imaging (fMRI) data analysis applies ICA. For example, in the article of [6], temporal dynamics and their spatial sources have been successfully recognized by real-valued ICA. In addition, the ICA has been applied to [7] for classification in electroencephalography (EEG) which has a two-state output (fatigue state vs. alert state). ICA has been also used in gait activity analysis, which usually relies on multiple sensors such as pressure, gyroscope, and accelerometer. The multiple sensors often incur crosstalk problem sensors where each sensor interferes with another sensor. ICA has also been exploited in [8] to enhance an automated classification technique to recognize toe walking gait from normal gait in idiopathic toe walking (ITW) children.

The simplest BSS model assumes the existence of n independent sources s₁, s₂,..., s _n , and the same number of linear and instantaneous mixtures of these sources, x₁, x₂,..., x _n , that is,

In vector-matrix notation, the above mixing model in the presence of noise can be expressed as

Here, A is an N × N square mixing matrix.

Equation (2) can be expressed in matrix form as follows:

The model described above is represented in Fig. 1.

The mixing and de-mixing models

Full size image

The de-mixing process [9,10,11] can be represented by calculating the separating matrix W, which is the inverse of the mixing matrix A, and computing the independent sources, which are obtained by

In this paper, we introduce a novel (differential) image separation algorithm that separates mixed images by extracting the components of the images using a pyramid technique. In this way, the image structure can be decomposed into multiple images of different scales. The proposed method creates different levels of scaled-down images in a pyramid structure. We therefore conduct the separation process on each level of the pyramid. While the lowest scale image at the top of the pyramid has same features, it incurs lower redundancy than the original image at the bottom of the pyramid. Our separation process conducted on the scaled images of the pyramid can, therefore, lead to better separation performance with lower redundancy in the resulting separated images.

Our method has the following two advantages over the most of ICA methods. Its first advantage is the high performance under noisy condition. Most of the ICA techniques consider only noiseless data; hence, they often lead to poor results in the presence of noise [12, 13]. In contrast, our method can separate the mixed image under a noisy condition and still provide high peak signal-to-noise ratio (PSNR). The second advantage is its fast processing and yet accurate separation results. Since it removes the redundancy in the image information, it can obtain the estimated image sources faster and more accurately than the ICA methods. The key contribution of proposed algorithm extracts the scaled-down images of the pyramid, in a way that maintains the important information of the original image, while reducing the redundant information.

The remainder of this paper is organized as follows. In Section 2, we provide the related work. Section 3 illustrates various techniques including the principle of pyramid image, finite ridgelet transform, and trigonometric transforms. In Section 4, the proposed image separation approach is presented. To demonstrate the effectiveness of the proposed technique, an extensive set of simulation experiments and performance comparison is reported in Section 5. Finally, Section 6 presents the concluding remarks.

In literature, there have been several papers published, which propose various approaches to the source image separation problem. The method in [14] considers a nonlinear real-life mixture of document images that occur when a page of a document is scanned and the back page shows through. It used a separation method based on the fact that the high-frequency components of the images are sparse and are stronger on one side of the paper than on the other one. Astrophysical image separation has been considered for a blind source separation method in [15]. In the work of [16], feedback sparse component analysis of image mixture was developed to extract the image sources by utilizing a feedback mechanism and sparse component analysis (SCA). In [17], a wavelet packet transform method was proposed in combination with a geometric de-mixing algorithm. It decomposes the mixed images by a wavelet transform (WT) and then uses the most relevant component as an input to its de-mixing geometric algorithm.

In the article of [18], columns or rows of mixed images were concatenated to arrange them into a 1-D mixed image. Then, a source separation of frequency-time approach with mutual diagonal was introduced to enhance these 1-D signals to resolve their components. Then, two-dimensional (2-D) astrophysical image components were achieved by segmenting separated 1-D original signals and rearranging these segments as columns or rows. Recently, researchers implemented sparse component analysis (SCA) to improve the method of blind image separation [19, 20]. These approaches could accurately separate the image mixtures using linear clustering when the linear clustering has less run time than super-plane clustering techniques, and the image sources are sparse enough [21].

The work of [22] applied the discrete cosine transform (DCT) as an approach to get the information in the frequency domain. It uses a block-segmented DCT reorganization to get the information in the segmented blocks while selecting the sparsest block by comparing the linear strength in each block. Moreover, the authors of [22] used the geometric characteristic of sparse blocks to study the linear orientations that match with the mixing matrix columns.

3.1 Principle of pyramid image enhancement

An image can be decomposed and analyzed in a form of a pyramid with a few levels of scaled-down images. The pyramid places the original image at the first level and adds scaled-down images at the higher levels [23] as illustrated in Fig. 2.

Pyramid technique effect

Full size image

The pyramid scales down an image using the low-pass filter with a Gaussian mask expressed by Eq. 5.

The motivation behind our pyramid technique is that the surrounding pixels within a certain area often have the similar characteristics, and thus, they are highly correlated with each other. To estimate the inverse matrix from Eq. (4) directly, the mixed image is converted from 2-D signal to 1-D signal. It is, therefore, inefficient to apply ICA algorithm to the pixel values, since most of the information around neighboring pixels is redundant, and the entropy of the pixels in the same area is low. Therefore, as a more efficient method to calculate the inverse matrix, we proposed a new technique that can remove the redundancy without affecting the information, and that can increase the entropy, while maintaining the features. The pyramid technique has been proposed in this paper, which scales down the image while maintaining the main features such as salient features and removing the redundant information. In the presented work, we use three levels to construct pyramid levels as shown in Fig. 2. Level 1 is the input image. Level 2 is the output after applying the filter based on Eq. (5), followed by a downsampling step. The above steps are then repeated to produce level 3. The ratio between the image outputs of the two consecutive levels determines the scale of the pyramid levels with respect to the original image. We can use these scales for further processing using ICA separation.

3.2 Transform techniques

3.2.1 Ridgelet transform

Ridgelet transform (RT) is a highly effective approximation approach to represent an image object as described by Candes and Donoho [24, 25]. It has a discontinuity across a line, and a curvelet, which is adopted in their papers as a type of RT, and is an effective transform for objects with discontinuities across curves. The approximation quality of RT is very close to the ideal Lagrangian condition and is in general better than any other algorithms such as Fourier transform (FT) and wavelet transform (WT). Due to such advantage, RT is widely used in image analysis, such as watermarking, image enhancement, image de-noising, and texture classification [26, 27].

Suppose that there is a function ψ:R → R satisfying the admissibility condition

where \( \widehat{\psi} \) stands for the Fourier transform of the function ψ. For each a > 0, b ∈ R, and θ ∈ [0, 2π], a bivariate RT ψ_a,b,θ:R² → R² is defined as

For a fixed θ, ψ_{a, b, θ}(x₁, x₂) is constant along the line x₁ cos θ + x₂ sin θ = constant.

Given an integrable signal f(x₁, x₂), the RT is defined as

It follows that

where ψ_{a, b}(t) = a^−1/2. ψ((t − b)/a).

Then the RT can be expressed as

As a result, the formula is represented by

From signal f(x₁,x₂), we can calculate Radon transform (RAT). Thus, the RT space can be expressed as an implementation of a 1-D wavelet transform to the slices of the Radon space.

It is known that approximate RAT for an image can be effectively computed with the fast Fourier transform (FFT). This approach is summarized by [28,29,30]:

1. (a)

   2-D FFT step: calculate the 2-D FFT of the image.

2. (b)

   Cartesian to polar conversion step: obtain samples on the recto polar as shown in Fig. 3.

3. (c)

   1-D inverse FFT step: calculate the 1-D inverse FFT on each angular line.

Ridgelet transform calculation using FFT in rectopolar

Full size image

For the implementation of the Cartesian to polar conversion, we use a rectopolar coordinate plane. The geometry of this coordinate plane is presented in Fig. 3, where the data points are marked with circles. Here, for an image of size n × n, there are 2n radial lines in the frequency plane selected by connecting the origin to the vertices lying on the boundary of the array. The grid lines of the rectopolar coordinate plane are the intersections between the set of radial lines and that of Cartesian lines parallel to the axes. Thus, there are 2n × n points (marked with circles) on the rectopolar grid lines, and the corresponding data structure is a rectangular format with n × 2n elements.

To complete the RT, we perform a 1-D wavelet transform along the radial variable in Radon space. Figure 4 displays the flow graph of the RT.

Ridgelet transform flowchart

Full size image

Do and Vetterli supposed a different execution of the ridgelet transform called finite ridgelet transform.

It has numerical exactness like the RT with little computational complications. As supposed above, a separate RT can be realized via a Radon transform and a 1-D discrete wavelet transform (DWT) as presented in Fig. 4.

The finite Radon transform (FRAT) is simply an addition of image pixels over a certain set of lines. Those lines are known in a limited geometry in a similar scheme to the lines for the constant Radon transform (RAT) in the Euclidean geometry [12,13,14,15,16].

We denote Z _p  = {0, 1, 2 ... p − 1}, where p is a prime number. Note that Z _p is a limited field with modulo p processes.

Then, the FRAT of a real function f on the limited grid Z _p ² is given by

Here, L_k,l denotes the group of points that make up a line on the lattice Z _p ².

The lines of the FRAT show a wrap-around effect in the transform. This means that the FRAT deals with the input image as one period of a periodic picture. In the FRAT domain, the energy is best compressed if the mean is removed from the image f(i, j) prior to the transform.

In Eq. (12), the factor \( 1/\sqrt{p} \) is supplied in order to normalize l₂ standard between the result and input of the FRAT. With an invertible FRAT and by using Eq. (13), we can have an invertible separate finite ridgelet transform (FRT) by taking the separate wavelet transform on each repetition of FRAT projection repetition, (r _k [0], r _k [1], … r _k [p−1]), where the trend k is constant. The total recording is known as the FRT as shown in Fig. 4.

3.2.2 Discrete wavelet transform

Wavelets have become an efficient tool in several signal processing areas such as signal de-noising, image fusion, and signal restoration and compression. The conventional discrete wavelet transform (DWT) may be regarded as the result of filtering the input signal with a bank of band-pass filters whose impulse responses are all approximately given by scaled versions of a mother wavelet. The scaling factor between adjacent filters is usually 2:1, which leads to octave bandwidths and center frequencies that are one octave apart from each other as illustrated by Fig. 5 [31]. The outputs of the filters are usually maximally decimated so that the number of DWT output samples equals the number of input samples and the transform is invertible as shown in Fig. 5.

Reconstruction and decomposition based on wavelet. a Full decomposition-reconstruction of two band filter bank. b Decomposition tree of wavelet packet. c Reconstruction tree of wavelet packet

Full size image

The art of calculating a good wavelet lies in the design of appropriate filters, H₁, H₀, G₁, and G₀, to realize various trade-offs between frequency and spatial space characteristics while satisfying the condition of perfect reconstruction (PR) introduced by [31]. In Fig. 5, the procedure of interpolation and decimation by 2:1 as the result of H₁ and H₀ defines all odd components of these signals to zero.

For the low pass branch, this is equivalent to multiplying x₀(n) by \( \frac{1}{2}\left(1+{\left(-1\right)}^n\right) \).

Hence, X₀(z) is converted to\( \frac{1}{2}\left\{{X}_0(z)+{X}_0\left(-z\right)\right\} \). Similarly, X₁(z) is converted to \( \frac{1}{2}\left\{{X}_1(z)+{X}_1\left(-z\right)\right\} \).

Thus, the expression for Y(z) is given by the equation below [31]:

The first PR condition requires aliasing cancelation and forces the above term in X(−z) to be zero [31].

Hence, {H₀(−z)G₀(z) + H₁(−z)G₁(z)} = 0, which can be achieved if:

Here, k is limited to odd numbers (usually k = ± 1).

From X(z) to Y(z), the transfer function need to be unity in the second condition of PR:

3.2.3 Trigonometric transform

The two primary trigonometric transforms are the discrete cosine transform (DCT) and the discrete sine transform (DST). Trigonometric transform has an energy compaction feature. The properties of these transforms are described below.

DCT

The DCT is a 1-D transform with the capability of energy compaction. For a 1-D signal x(n), an application example of DCT is given by [32].

where

DST

The DST is another transform and can be calculated by Eq. (17). Application examples of the DST can be found in [31]:

As described in the prior sections, we merge the benefits of the pyramid technique and FRT. First, the mixed images are decomposed into frequency bands using different transforms. Then, each frequency band is handled, separately, using the pyramid technique to extract its details.

A flow diagram of the proposed method is depicted in Fig. 6, which is also described by the following steps:

• Step 1: Decompose the mixed image into different transforms using the finite ridgelet transform (FRT), wavelet transform (WT), discrete sine transform (DST), and discrete cosine transform (DCT).

• Step 2: Apply a pyramid construction on each transform to obtain the different scale components of each transform in each pyramid level. We chose three levels of pyramid construction in the present work, while it can be extended to a larger number of levels.

• Step 3: Conduct a separation operation on all the pyramid levels in each transform and calculate the inverse matrix (un-mixing matrix). The operation proceeds from level 3 (the smallest scale) towards level 1 (the largest scale). In level 3 of pyramid component, we start with a random matrix to calculate the values of inverse matrix. The output values of estimating inverse matrix from level 3 are used as the input matrix to update the inverse matrix values for level 2. The updated output of the inverse matrix from level 2 is in turn used as the input matrix for level 1 to calculate the final values of the inverse matrix. The final estimated values of the inverse matrix are applied to the original mixed image to extract accurate separated images in step 4.

• Step 4: Calculate an estimate of the separated image using the mixed image with the calculated inverse matrix.

Block diagram of the proposed image separation algorithm

Full size image

ICA has been regarded as one of the most efficient approaches reported in different fields. It provides advantages of fast convergence and straightforward implementation. Figure 7 illustrates a flowchart that realizes the fast independent component analysis (FastICA) approach [32]:

Flowchart of the FastICA algorithm

Full size image

Next, we introduce several performance metrics to evaluate the image separation results for each given method, which are also described in [33–35].

1. (i)

   Signal-to-noise ratio (SNR):

Here, the size of the image is N × M, while f(x, y) represents the original image and \( \tilde{f}\left(x,y\right) \) an estimated image.

1. (ii)

   Root mean square error (RMSE):

RMSE is to measure the square error between two images. It considers image degradation as perceived variation in information.

1. (iii)

   Peak signal-to-noise ratio (PSNR):

PSNR is one of the most widely used metrics for evaluating the quality of estimated image. The higher the PSNR values are, the higher quality the estimation output provides.

1. (iv)

   Normalized cross-correlation (NCC):

NCC is another common performance metric that is useful to compare the estimation results from different source images.

A computer simulation is presented in this section to evaluate the performance of the proposed approach after the image was mixed. In all experiments, test images were used which were extracted from a standard image database. We assumed that the mixed images are corrupted by an additive white Gaussian noise (AWGN) with zero mean and unit variance to illustrate the visual aspect of the various mixed images; we reported in Fig. 9 one of each noisy mixture in several noise levels.

We have conducted experiments on the images of Fig. 8 and obtained better image separation results from the proposed separation method compared with other separation methods. Due to the space restriction of the paper, we summarize the detailed experiment results with Cameraman and Baboon images. Figure 10 shows the experimental results with Cameraman and Baboon images using the proposed separation method and various other methods.

Original images, images from left are Cameraman and Baboon

Full size image

These test images are created by a convolutional mixing process using a set of mixing matrices generated randomly by MATLAB, and the criteria of this matrix are normally distributed random numbers. Also, Fig. 9 shows the result of mixing process at different noise level.

Mixing results at several noise level. a 4 dB. b − 5 dB. c − 10 dB. d − 15 dB

Full size image

As a result of the experiments, the separated images are shown in Fig. 10. The numerical results of these experiments are included in Tables 1, 2, 3, and 4 at noise level 4 dB. These tables give an image quality comparison between separation algorithms, revealing the superiority of the proposed algorithm at noise level 4 dB. We used SNR, PSNR, RMSE, and NCC to evaluate the estimation quality of separated images.

Estimated results of separated image at noise level 4 dB. a Proposed method (FRT with pyramid). b FRT without pyramid. c DWT with pyramid. d Time with pyramid. e DWT without pyramid. f Time without pyramid. g DCT with pyramid. h DST with pyramid. i DCT without pyramid. j DST without pyramid

Full size image

As illustrated by Fig. 10 and from Tables 1, 2, 3, and 4, the proposed approach performs preliminary separation better than the conventional separation methods based on the time domain, wavelet transform, and trigonometric transform. We can observe that the resulting images of Fig. 10a separated by the proposed RT with homomorphic method has better quality than other images of Fig. 10b–j separated by various different methods.

From Tables 1, 2, 3, and 4, it can be observed that the proposed RT with pyramid operation produces higher quality and efficiency compared with all other approaches tested in our experiments. Tables 1 and 2 prove that the SNR and PSNR of the images separated by the proposed method are higher than those of all the other methods. Table 3 indicates that the RMSE of the resulting images separated by the pyramid operation with FRT is the best compared with all other methods considered in our experiments. In addition, as illustrated in Table 4, the NCC of the proposed technique shows a value closer to 1 than any other methods do. Here, an NCC result of 1 is the best possible result. From the experimental results of Tables 1, 2, 3, and 4, it is observed that the separation quality of the time domain-based method is relatively low. This result follows from the fact that the sources must satisfy statistical independence to allow FastICA methods to achieve high-quality separation results.

Figures 11, 12, 13, 14, 15, 16, 17, and 18 plot an extensive set of simulation results measured with a wide range of noise levels. These results compare the separation quality of the tested images using various evaluation metrics including SNR, PSNR, RMSE, and NCC. Figures 11 and 12 present the SNR output of the separated images for Cameraman and Baboon, respectively. These figures demonstrate the performance comparison of the proposed BSS method with respect to the other methods at different input noise levels. It shows that the proposed method provides the highest performance. For example, the proposed FRT with pyramid method achieves an SNR of higher than 7 dB for the Baboon image, whereas the DCT without pyramid method gives an SNR as low as 0.5 dB. Figures 13 and 14 illustrate the PSNR measurement of the separated images, where the proposed method can obtain PSNR values of 15 dB or higher, whereas other methods provide much poorer PSNR values in the range of 4~14 dB.

Output SNR vs. input SNR for Cameraman image overall separation performances

Full size image

Output SNR vs. input SNR for Baboon image overall separation performances

Full size image

Output PSNR vs. input SNR for Cameraman image overall separation performances

Full size image

Output PSNR vs. input SNR for Baboon image overall separation performances

Full size image

Output RMSE vs. input SNR for Cameraman image overall separation performances

Full size image

Output RMSE vs. input SNR for Baboon image overall separation performances

Full size image

Output NCC vs. input SNR for Cameraman image overall separation performances

Full size image

Output NCC vs. input SNR for Baboon image overall separation performances

Full size image

On the other hand, Figs. 15 and 16 demonstrate the output of RMSE. In the RMSE results, the lowest curve indicates the best result. We can observe that the proposed method produces an RMSE value of 0.25 or lower, which is 0.2~0.6 lower than other methods. Figures 17 and 18 illustrate the simulation results of NCC. It is observed that the proposed method produces an NCC value very close to 1, while other methods provide much NCC values in the range of − 0.1~ 4.

This paper addressed the blind image separation problem by introducing a new image separation technique based on a novel concept of pyramid processing and ridgelet transform. The proposed approach first uses FRT domain coefficients to obtain the frequency components. It then applies a pyramid processing to estimate the mixing matrix by constructing the different level scales to extract more details in information and remove redundant information. We conducted an extended set of simulation experiments using various image separation methods that employ the proposed FRT as well as other methods including DWT, time domain, DCT, and DST along with pyramid and non-pyramid operations, respectively. The experimental results demonstrate that the proposed method outperforms all other methods that we tested. In summary, it presents PSNR values of 12~16 dB under a wide range of noise condition, while all other methods provide much poor PSNR values of 4~10 dB under the same noise condition. The proposed method, therefore, appears to be an efficient approach to separate mixed images even under noisy conditions.

BSS:

Blind source separation

DCT:

Discrete cosine transform

DST:

Discrete sine transform

DWT:

Discrete wavelet transform

EEG:

Electroencephalography

FastICA:

Fast independent component analysis

fMRI:

Functional magnetic resonance imaging

FRAT:

Finite Radon Transform

FRT:

Finite ridgelet transform

ICA:

Independent component analysis

ITW:

Idiopathic toe walking

NCC:

Normalized cross-correlation

PR:

Perfect reconstruction

PSNR:

Peak signal-to-noise ratio

RAT:

Radon Transform

RMSE:

Root mean square error

RT:

Ridgelet transform

SNR:

Signal-to-noise ratio

WT:

Wavelet transform

This work was supported by IITP grant through the Korean Government, development of wide area driving environment awareness and cooperative driving technology which are based on V2X wireless communication under grant R7117-19-0164, and it was also supported by the Center for Integrated Smart Sensors funded by the Ministry of Science of Korean Government, ICT and Future Planning as Global Frontier Project (CISS-2016).

Authors and Affiliations

1. Department of Electronic Engineering, College of Electrical and Computer Engineering, Chungbuk National University, Cheongju City, South Korea

   M. Y. Abbass & HyungWon Kim

2. Engineering Department, Nuclear Research Center, Atomic Energy Authority, Cairo City, Egypt

   M. Y. Abbass

Contributions

MYA and HWK designed the proposed algorithm together. MYA implemented it with MATLAB. Both authors wrote and approved the final manuscript. The corresponding author is HWK (hwkim@cbnu.ac.kr).

Corresponding author

Correspondence to HyungWon Kim.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Reprints and permissions