1. Introduction

Advances in digital sensors, communications, computation, and storage have created huge collections of data. Trillions of bytes of data every day such as satellite images, driving directions, and image retrieval continually add new services. Digital data are being generated by many different sources, including digital imagers (telescopes, video cameras, MRI machines). Data from the many different sources can be collected into massive data sets via localized sensor networks, as well as the Internet. Understanding the environment requires collecting and analyzing data from thousands of sensors monitoring, this is big data environment. For example, modern medicine collects huge amounts of data about patients through imaging technology (CAT scans, MRI), genetic analysis, and other forms of diagnostic equipment [1-3]. By applying data mining to data sets for large numbers of patients, medical researchers are gaining fundamental insights into the genetic and environmental causes of diseases, and creating more effective means of diagnosis. Big-data medical image is one of important research focuses.

In medical field, embedded numerous sensing nodes (as Fig. 1) can be used to monitor and interact with physical world, such as Pulse Oximeter, Glucose Meter, Electrocardiogram (ECG), CT, MR, Social Alarm Devices, and so on. Wireless Sensor Network (WSN) can be used in all kinds of medical applications [4-7]. WSN can do telemonitoring of human physiological data, tracking and monitoring patients and doctors inside a hospital, drug administration in hospitals, and so on.

Fig. 1.Embedded numerous sensing nodes in WSN

In WSN, data coming from multiple wireless sensor nodes are aggregated and they reach the same routing node to the sink [8-9]. The application of medical image fusion (as Fig. 2) for big-data computing has strong development momentum, especially in medical therapeutics, so the fusion method study is a hot topic. Generally, the fusion for big-data medical image can be divided into 4 levels [10-11]: signal, pixel, feature and symbol. According to fusion different levels, there are relative specific fusion methods. For the pixel and the feature, the fusion methods have principal component analysis, pyramid-based, wavelet transform, etc.

Fig. 2.One of medical image sensor fusion paradigms

Discrete Wavelet Transform (DWT) can decompose the image into different scales. The wavelet coefficients have the following characteristics [12-13]: low-frequency subbands contain most of the energy of the original image, while the high-frequency coefficients take less energy and mostly in the image edge.

Said [8-9] presented encoding implementation based on set partitioning in hierarchical trees (SPIHT), which is highly efficient in coding non-textured images. This encoding algorithm can be stopped at any compressed file size or let it run until the compressed file is a representation of a nearly loss-less image. However, SPIHT coding has a serious drawback: repeated scan. This disadvantage makes encoding efficiency less. In recent years, many researchers [10-16] proposed different improved SPIHT algorithm. However, those improvement methods for SPIHT algorithm only optimized the code scheme, but neglected the transform part.

We improve the SPIHT algorithm and optimize the transform part in this paper. In transform part, after analyzing the characteristics of the wavelet coefficients, we introduce a strategy to pre-process the coefficients. Firstly, the three high-frequency wavelet coefficients are transformed to spherical coordinate domain. In this way, the correlation between wavelet coefficients in the same decomposition scale can be further reduced. However, this strategy does not affect the correlation between wavelet coefficients in the different decomposition scales. Multi-Scale Model Product (MSMP) is used to control the wavelet shrinkage functions, which can remove the coefficients that have little effect on the visual and some noise information. This pre-processing strategy can reduce the redundant ratio of big-data medical image. In coding part of medical image, the low-frequency part contains the most energy of the image, and this part does not participate into the spherical transformation. Therefore, it can be encoded separately by DPCM. The high-frequency parts in spherical coordinate domain adopt the improved SPIHT coding algorithm. We introduce a matrix of maximum pixel (MMP) to overcome the re-scan disadvantage of traditional SPIHT algorithm. When the threshold is updated, we just compare the new threshold to the coefficients in MMP. This mechanism avoids the comparison between the new thresholds to every pixel. The improved SPIHT algorithm greatly improves the fusion efficiency.

The medical image fusion with coding based on spherical coordinate domain contains the following steps: the image wavelet decomposition in different scale, the three high-frequency wavelet coefficients transforming to spherical coordinate domain, removing the coefficients by MSMP and encoding the image, transmitting to fusion center by WSN, making fusion in center center, reconstructing the image by the processed coefficient. The fusion rules of medical image are basically expressed by formula (1)

where, I(x, y) is the fused image, I1 (x, y) and I2 (x, y) are two input images, w(•) is the discrete wavelet transform method, w−1(•) is the inverse discrete wavelet transform method, ϕ is the fusion rule.

 

2. Related Works

As we know, sample-based image fusion can be divided into pixel-based image fusion and patch-based image fusion [1-3]. Generally, pixel-based image fusion can easily reflect image diversities but can’t keep image structures and it is slow, as patch-based image fusion fast and can better maintain image characteristic information within blocks. Patch-based image fusion method is based on sample-based image fusion method, which divides the input image into blocks with fixed size. The constraint rule of overlapping regions between blocks is to select matching blocks and generate images that without repeat and continuous on arbitrarily large size [4-7].

As we know, SPIHT algorithm is based on bit plane. It organizes the wavelet coefficients into spatial orientation trees. SPIHT with lists algorithm uses three different lists to store significant information of wavelet coefficients for image coding purpose [8-16]. Three lists are list of insignificant sets (LIS), list of insignificant pixels (LIP), and list of significant pixels (LSP). At first, SPIHT combines nodes into one set which is denoted as insignificant. With traveling each tree node, sets in the LIS are partitioned into four different subsets which are tested for significant state.

 

3. Spherical Coordinate Transform

3.1 Wavelet domain to spherical domain

The human visual system is characterized in observing an image. We are more sensitive to low frequency of an image, while less sensitive to the high frequency parts. Based on this characteristic of the human eyes, discrete wavelet transform decomposes the image into different scales and processes the image in different resolutions. The coefficients in the same scales are great relevance. If can reduce the correlations and reserve the main information, we will achieve high compression ratio when compressing an image. And we will reduce the image storage space and improve the efficiency of the transmission channel. From this viewpoint, the three high frequency coefficients in wavelet domain are transformed into spherical domain.

The signal f(x) denotes an image. After DWT transformation, the image decomposes into a low frequency and three high frequencies. is the horizontal high frequency, is the vertical high frequency and is the diagonal high frequency. This paper transforms the three high frequencies to spherical coordinate domain

where j, k denote the coordinate of the coefficients.

After the spherical transformation, the coefficients in the wavelet domain are mapped to spherical coordinate domain. The corresponding relationship between the various components is shown as Fig. 3.

Fig. 3.Transform from wavelet domain to spherical domain

Spherical transform can further remove the correlations of the wavelet coefficients in the same scales. The energy of the three high frequencies in the wavelet domain is totally transformed to the radial components in the spherical domain. While the phase angles just retain the pixel location information. Therefore, we just process the radial components and decrease the complexity of the algorithm. Also it avoids directly process the wavelet coefficients.

3.2 Multi-scale module products

Multi-scale norm products (MSNP) in spherical coordinate domain are defined as follows.

Definition 1 The MSNP of pixel point (x, y) in spherical coordinates is as

where k is the decomposition scales.

The MSNP has the following property: MSNP can enlarge the distance between effective signal and unimportant signal. It avoids the shortcoming that smaller wavelet coefficients are wrongly removed by traditional methods. Therefore it will effectively remove the noise information and unimportant coefficients without loss of the important feature of the original image.

Generally speaking, the image edge has the following features:

1) Has a strong gray-level mutations, a stark contrast with background.

2) The edge points may form a meaningful linear relationship, namely, there is an ordering characteristics between the adjacent edge points.

3) Processes a directional feature.

4) The relative position in the image.

5) Edge types include four kinds: pulse, step, slope and roof.

Definition 2 (the edge point for the coordinates (i, j)), the characteristic elements can be defined as the following quintuple.

where B is the contrast set; G is the gradient strength of edge point; it reflects the contrast degree between image edge and image background; L is the linear chain length of edge point, it reflects the correlation of image structure; θ is the direction angle of edge point; D is the distance between edge point and the other long chain edge (large profile), it is also the relative position in the space; S marks which type the edge is, step or roof, 0 for step type, 1 for roof type.

Theorem 1 The parameter threshold T in spherical coordinates domain for transforming from wavelet domain is

where N denotes image pixel, R is radial component in spherical coordinates.

Proof Chamboll proposed a better shrinkage threshold in Besov space [8-11] as follows

where N is the number of image pixel, R is radial component in spherical coordinates system. The equivalent relationship between Besov space norm and radial component norm is established in three steps

① Firstly, to establish equivalent relationship between Besov space norm and wavelet coefficients norm, i = 1,2,3

Suppose f is a noise image in the interval Ω = [0,1)2. It is proper to use Besov space norm to describe the regularity of image. Firstly, the image f is decomposed by compactly supported orthogonal wavelet Daubechies (db2). Suppose φj, k = 1 denotes a cycle scaling function. Ψ = {ψ(1), ψ(2), ψ(3)} is a set of orthogonal wavelet bases in the space , where

The DWT of image f is

where for φj, k = 1 , the formula (8) can rewrite as

Thus, the equivalent relationship between Besov space [8-11] norm and wavelet coefficients norm has been established.

② Secondly, to establish equivalent relationship between the threshold T in Besov space and wavelet coefficients norm , i = 1, 2,3 .

Set q = 1 , shrinkage threshold formula (6) can rewrite as follows

Therefore

Generally, dyadic wavelet is used to deal with discrete image signal, for α ≤ 1 , set , the formula (11) can rewrite as

③ Finally, to establish equivalent relationship between wavelet coefficients norm and radial component norm |R| in spherical coordinates domain.

The noise image f is decomposed on multi-scale wavelet. The high frequency parts , i = 1, 2,3 are mapped to θ,φ and radial component R in spherical coordinates system. So the equivalence relationship between radial component R and high frequency needs to be established.

Based on the spherical transform in wavelet domain, the following inequality will be established

Therefore,

where, i = 1,2,3 k ≥ o, j ∈ {0,1,…2k - 1}2.

So the following relationship will be established between wavelet coefficients norm and radial component norm.

The interval of shrinkage threshold T will be

The least upper bound of T is chose as the final shrinkage threshold

where N denotes the number of image pixel.

3.3 The processing of wavelet coefficients

In the spherical domain, R(j, k) takes the energy of the pixel point (j, k). The greater the radial components |R(j, k)| is, the more important of the image information which contains. The radial components whose absolute value is small contain the secondary information and some noise information. In order to achieve large compression ratio and do not affect the visual effect, this paper processes the radial components by soft-threshold shrinkage

Where T denotes the wavelet shrinkage threshold. We use the Donoho’s [8-11] classical threshold , σ is the noise intensity of the image, which is estimated by σ* = Median(|f|)/0.6745.

Different with the traditional shrinkage strategies, we take the threshold acting on the MSMP. This algorithm avoids killing the small amplitude coefficients and can effectively preserve the useful information of the original image.

In the spherical coordinates, angle and the radial components are independent to each other. Small angle does not always correspond to small energy. If the R(j, k) is zero, this indicate the point (j, k) signal energy is zero. Therefore, the corresponding two components θ(j, k) and φ(j, k) can be set to zero, which does not affect our recognition.

 

4. Medical Image Coding Algorithm

4.1 Low-frequency part coding

After wavelet transform, the energy carried by the low frequency almost closes to 90% of the original image. In order to achieve a higher quality of the reconstructed image, we use the DPCM coding to process the low frequency coefficients.

The size of the original image f(x, y) is M ∗ N , LL (j, k) denotes the value of the pixel point (j, k) in the low frequency wavelet coefficients. We use the row and the column to estimate other wavelet coefficients (as Fig. 4).

Fig. 4.DPCM coding in low-frequency part

Firstly, keep LL(j, N/2) and LL(M/2, k) unchanged. Other coefficients are estimated by the following formula.

where j = 1,⋯, M k = 1,⋯,N.

Secondly, estimate the row and the column :

Using the following formula to recover the low frequency coefficients

4.2 High-frequency part coding

After being pre-processed, the coefficients of the high frequencies are transformed to spherical coordinate domain. The correlation between coefficients in the same scale is reduced. The radial components R(j, k) are processed by soft threshold function. A lot of zero points are produced. This promises a high compression ratio.

4.2.1 Parameter definition

MMP: Matrix of Maximum Pixel.

H : set of coordinates of all spatial orientation tree roots (nodes in the highest pyramid level);

O(j, k): set of coordinates of all offspring of node (j, k).

D(j, k): set of coordinates of all descendants of the node (j, k).

L(j, k): set of coordinates of all non-child descendants of the node (j, k), where L(j, k)=D(j, k)−O(j, k)

T(j, k): Spatial orientation tree whose root is the node (j, k).

T : maximum scan threshold.T = 2n , n0 = ⎣log2 ( max j,k{⎢ci,j⎥})⎦ is the initial scan threshold.

4.2.2 Improved SPIHT algorithm

In traditional SPIHT algorithm, H denotes the wavelet coefficient of the low frequency. While in this paper, H represents the coefficients in the maximum scale of the high frequency. When discussing the importance of the set D(j, k) and L(j, k), once the threshold updated, all descendant coefficients of node (j, k) will compare to the new threshold. There are a lot of repeat operations that result in lower coding efficiency. We introduce the maximum matrix of pixel MMP to solve the problem of repeated comparison in scanning.

In the improved algorithm, when scanning the list LIS, the matrix MMP will be established to preserve the largest pixel nodes of descendant of the set D(j, k) and L(j, k). Once the threshold updated, we just need to compare the new threshold with the corresponding coefficients of MMP. This method can effectively reduce the number of comparisons and improve scanning efficiency. The function used for testing set significant state is defined by the following formula:

Where M is a set of coordinates.

Improved algorithm is as follows

(1) Initialization: output n0 = ⎣log2 ( max j,k{⎢ci,j⎥})⎦; set the LSP and MMP as empty lists, add the coordinates (j, k) ∈ H to the LIP, and only those with descendants also to the LIS, as type A entries.

(2) Sorting Pass:

I. For each ∀(j, k) ∈ LIP , output SM (j, k), if SM (j, k) = 1 then move (j, k) to the LSP and output the sign n of cj,k.

II. For each ∀(j, k) ∈ LIS , different from the traditional SPIHT algorithm, this part introduces MMP to preserve the largest pixel coordinates of descendants. Detailed scanning process references the algorithm flowchart.

If the type of set is D(j, k), judge whether n = n0 .

If this equation is valid, compare the node (j, k) ∈ (j, k) to the threshold and update the MMP. If this equation is not valid, just need to compare the node (j, k) ∈ MMP to the threshold.

If the type of set is L(j, k), compare the node (j, k) ∈ MMP to the threshold and update the MMP.

Refinement Pass: for each entry (j, k) in the list LSP, except those included in the last sorting pass (i.e., with the same n), output the nth most significant bit of |cj,k|.

Quantization-step Update: decremented n by 1 and go to step 2.

4.3 Fusion approach for medical image

We adopt two-layer directed graph method with multiscale edge. The multi-layer situation can be divided into two layers to handle. The method is that two images merge into an image, actually becomes two graphics fuse into a graphic. The fusion method is as follows.

1) Pre-process the original image to eliminate all interference by de-noising and smoothing process.

2) Medical image in each sensor does wavelet transformation in wavelet domain. Then the high-frequency coefficients of each image are transformed to the spherical coordinate domain to reduce the correlation in the same scale.

3) Multi-scale model product (MSMP) is used to control the shrinkage function to make the small wavelet coefficients and some noise removed. The high-frequency parts in spherical coordinate domain are coded by improved SPIHT algorithm. The low-frequency part contains the most energy of the image, which can be encoded separately by DPCM.

4) The coded image data in each sensor should be transmitted to fusion center by WSN. In fusion center, the coded image data should be decoded.

5) Do edge processing for the decoded image by the above multi-scale edge detection method. Construct the directed graph G1, G2 by the above fusion method. Based on the resulting directed graph G1, G2, mark its nodes or edges with the number (Pi , Ej), calculate the close distance of the nodes (edges) in two graphs.

6) According to the calculated close distance, combine the points within a certain threshold T1, delete the points without a certain threshold T2 (T2 >T1).

7) By the principle that the important and continuous boundary generates the longer edge curve, we check the resulting graphical boundary by step 5). Is it OK? If it is, then the resulting image by this fusion method is valid, otherwise, repeat step 5) to continue.

8) Based on the multi-scale edge of medical image, the image can be reconstructed.

Fig. 5 shows the chain of wavelet decomposition and fusion by two layers directed graph method. It expresses the fused results of the head CT image (gray edge) and MR image (black edge) under the transform scale j=3, the right graph is the enlarged chain one in this region.

Fig. 5.Fusion by two layers directed graph method

 

5. Experiment and Test

Based on the above correlation fusion algorithm, we have done tests and analysis for the fusion application of the CT & MR image, in order to test the feasibility and applicability of this method.

Fig. 6 is the original colorful and colorless medical images CT and MR for experiments. It are full of noise images. Based on this images, they can be processed according to the aforementioned fusion steps. Fig. 7 is the fusion result of the CT & MR medical images. In Fig. 7, (a) is original CT medical image, (b) is original MR medical image. (h) is fused result by our method of this paper. (c), (d), (e), (f) and (g) are fused results by the methods from reference 8, 9, 10, 11 and reference 12, respectively. From the result of the examples of medical image fusion, we can see that our method has good quality.

Fig. 6.The original medical images for experiments

Fig.7.Fusion result of the CT & MR medical image

Based on the fusion results, we can know the above presented multi-scale correlation algorithm is feasible. In order to confirm the availability of the proposed method, we have done a large number of experiments. In comparison with the belief degree of other methods, the belief degree of the proposed method in this paper is much higher. In addition, by our project, we has done an on-the-spot practical application, it has been proved that its belief degree is as high as 95% or more, it is now clear that for some application fields of demanding the high fusion quality like medical diagnosis, the method is very valuable.

5.1. Comparison of relative method

We do the comparison of relative method and give the experimental results of big-data medical images in PSNR, MSE and run-time.

Fig. 8 is the comparison of MSE for kinds of de-noising methods to Gaussian white noising CT medical image among our method and several relative methods from aforementioned references. Fig. 9 is the comparison of PSNR for kinds of de-noising methods to Gaussian white noising of the same MR medical image among our method and several relative methods from aforementioned references. Fig. 10 is fusion run-time comparison of CT & MR image among our method and several relative methods from aforementioned references. Fig. 10 shows that our method retains the advantages in fusion quality and run-time in WSN. The results indicate the novel method is effective and very useful for transmission of big-data medical image (especially, in the wireless environment). From the comparison of PSNR and MSE, we can see that our method in this paper has better results.

Fig. 8.Comparison of MSE for de-noising methods to Gaussian white noising CT medical image

Fig. 9.PSNR comparison of MR medical image

Fig. 10.Fusion run-time comparison of CT & MR image

 

6. Conclusion

We has proposed a kind of new fusion approach with coding based on spherical coordinate domain in Wireless Sensor Network for big-data medical image in this paper. In this approach, the three high-frequency coefficients in wavelet domain of medical image are pre-processed. This pre-processing strategy can reduce the redundant ratio of big-data medical image. A multi-scale model product (MSMP) is used to control the shrinkage function so as to make the small wavelet coefficients and some noise removed. The high-frequency parts in spherical coordinate domain are coded by improved SPIHT algorithm. Based on the multi-scale edge of medical image, it can be fused and reconstructed. Experimental results based on CT & MR medical image indicate the novel approach is effective and very useful for transmission of big-data medical image, which can be especially used in the wireless environment.