Ferréol Soulez

A "LEARN 2D, APPLY 3D" METHOD FOR 3D DECONVOLUTION MICROSCOPY

This paper presents a 3D deconvolution method for fluorescence microscopy that reached the first place at the "the 3D Deconvolution Microscopy Challenge" held during ISBI 2013. It uses sparse coding algorithm to learn 2D "high resolution" features that will be used as a prior to enhance the resolution along depth axis. This is a three steps method: (i) deconvolution step with total variation regularization, (ii) denoising of the deconvolved image using learned sparse coding, (iii) deconvolution using denoised image as quadratic prior. Its effectiveness is illustrated on both synthetic and real data.

INTRODUCTION 1.Motivation

Thanks to the widespread of fluorescent labeling techniques, fluorescent microscopy is one of the most used imaging modality in biology and epifluorescence (a.k.a. wide field) microscopes are the simplest set up for such observations. It consists on imaging at its emission wavelength a cellular structure marked by fluorescent dye excited by uniform illumination. On the resulting 2D image, structures are more or less defocalized according to their distance to the focal plane. Moving this focal plane through the sample produces a 3D representation of the object. However, its characteristics both in terms of resolution and noise are often too coarse to fit in constraints of current biological research. Furthermore the resolution is not isotropic and the actual resolution along the depth axis is far more coarser than the resolution along lateral dimensions. Beside more advanced microscope designs such as the confocal or multiphotons microscopes, one can use only computational methods numerically invert the blurring process and enhance epifluorescence micrographs as demonstrated by Agard and Sedat [START_REF] Agard | Optical sectioning microscopy: cellular architecture in three dimensions[END_REF]. Since, many deconvolution algorithms were proposed (see [START_REF] Sibarita | Deconvolution microscopy[END_REF][START_REF] Sarder | Deconvolution methods for 3-d fluorescence microscopy images[END_REF] for reviews) and these "deconvolution microscopy" methods lead to various commercial and open-source software [START_REF] Griffa | Comparison of deconvolution software in 3d microscopy. a user point of view part 1[END_REF][START_REF] Griffa | Comparison of deconvolution software in 3d microscopy. a user point of view part 2[END_REF] but still suffer of several drawbacks: computational and memory costs, restoration artifacts, nonquantitative results. . . In the 2013, the first edition of the "3D Deconvolution Microscopy Challenge" was organized to stimulate the community to look for novel, global and practical approaches to this problem. I present in this paper the 3D deconvolution method that was awarded in this challenge.

Notations

In this document I use the following notations:

• y ∈ R N is the observed data;

• N = Nx × Ny × Nz is the size of the observed data; work supported by projects MiTiV and POLCA funded by the French National Research Agency (ANR DEFI 09-EMER-008-01 and ANR-10-BLAN-0511).

• g ∈ R K is the processed data;

• K = Kx × Ky × Kz is the size of the processed data; • k = 1, . . . , K is the voxels index; • H ∈ R K×K is the forward operator; • w ∈ R K is the vector of statistical weights; • u + ∈ R K is the solution of the first step; • u ++ ∈ R K is the solution of the second step; • u +++ ∈ R K is the solution of the last step; • s ∈ {1, 2, 3}
is the 2D subspace index corresponding to the sections {(x, y), (x, z), (y, z)} respectively;

• m = 1, . . . , M [s] is the pixel index in patches of size M [s] ;

• p ∈ R (M [1] +M [2] +M [3]
)×K is the patch vector with p m,k,s the value of the m th pixel of the patch of the orientation s of the slice centered on pixel k; • P [s] is the set of all overlapping patches laying in the subspaces s; • = 1, . . . , L is the atom index in dictionary containing L atoms;

• D ∈ R L×(M [1] +M [2] +M [3]
) is the dictionary; • D is the set of dictionary with L atoms of size M with unit

2 norm; • α ∈ R L×K×3 is the vector of coefficients; • R ∈ R K×K×M ×3 is a reconstruction operator; 2. ALGORITHM 2.1. Principle
As mentioned in the introduction, the resolution along the depth axis is coarser (at best 2 times) than the resolution along lateral dimensions. In most deconvolution method, the resolution of the deconvolved data is still coarser in depth than laterally. The key idea of this work is to extract "high resolution" features from lateral (x, y) sections and use them as a prior to enhance the lower resolution along depth axis. It means that I implicitly rely on that object structures are statistically isotropic (no preferred orientation).

This "copy -paste" method is achieved using sparse coding on a learned dictionary [START_REF] Mairal | Online dictionary learning for sparse coding[END_REF][START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF]. In sparse coding, every patch of the denoised image is a linear combination of few dictionary elements. In learned sparse coding, this dictionary is learned on a training set or directly using noisy image patches. The main contribution of this work is to learn the dictionary only using 2D patches extracted from lateral (x, y) sections of the 3D image of the previously deconvolved image. This learned dictionary is then used to denoise every 2D images of every sections, namely (x, y), (x, z) and (y, z), which are then averaged.

The algorithm we propose here has three steps: data

Step 1

Step 2

Step 3 (x, y) section (x, z) section Fig. 1. Results of the three steps of our method on parts of the green channel of the first "3D Deconvolution Microscopy Challenge" data set.

1. deconvolution using Total Variation regularization, 2. denoising of the deconvolved image using learned sparse coding, 3. deconvolution using denoised image as quadratic prior.

The results of these three steps are illustrated on figure (1) on a zoom of the green channel (with filaments) of the deconvolution challenge data set.

Step 1: Total variation deconvolution

This first step is exactly identical to the deconvolution step described in [START_REF] Soulez | Blind deconvolution of 3d data in wide field fluorescence microscopy[END_REF]. The maximum a posteriori estimates u + obtained by minimizing a cost function under a non-negativity constraint:

u + = arg min u≥0 (J data (u) + µ1 J TV (u)) , (1)
with J data (u) a data fitting term, J TV (u) a regularization term and µ1 an hyper-parameter.

Data-fitting term

Under uncorrelated Gaussian noise assumption, the data fitting term writes:

J data (u) = k w k [H • u] k -g k 2 , (2)
where w k is the inverse of the noise variance at voxel k. This model can cope with non-stationary noise and can be used to express confidence on measurements on each voxel of the data. Thus it can deal with unmeasured voxels (due to saturation, borders. . . ) by setting w k = 0 for such voxels. With such a formulation, it is possible to take rigorously the borders into account. In order to apply H using FFT-based computations, most works suppose that H is circulant. Doing so, to keep a rigorous direct model without the unrealistic periodic assumption, both data and estimated object must be extended. This field extension is mandatory for two reasons: (i) to prevent spatial aliasing and (ii) to estimate parts of the scene outside of the field of view that contributes to the measurements inside the instrument field of view [START_REF] Soulez | Heterogeneous Multidimensional Data Deblurring[END_REF]. For that reason, as the PSF has the same size as the data, we extend the reconstructed image u and the data g to at least K = (2 Nx -1) × (2 Ny -1) × (2 Nz -1), padding g with zeros and setting w k = 0 in the extended area. Furthermore, except for very low detector noise, this formulation can account for mixed Poisson + Gaussian noise by approximating it as a non-stationary uncorrelated Gaussian noise [START_REF] Mugnier | Mistral: a myopic edge-preserving image restoration method, with application to astronomical adaptive-optics-corrected long-exposure images[END_REF]:

w k def = γ max(g k , 0) + σ 2 k -1 if g k is measured, 0 otherwise, ( 3
)
where γ is the quantization factor of the detector and σ 2 k is the variance of other approximately Gaussian noises (e.g., read-out noise) at voxel k.

TV regularization term

As most observed objects are smooth with few sharp structures (e.g., edges and spikes), I use as regularizing prior an hyperbolic approximation of the classical 3D total variation [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]:

J TV (u) = k ∇ k u 2 2 + 2 , (4)
with ∇ k the finite difference operator to approximate spatial gradient at voxel k. Parameter > 0 ensures differentiability of J TV at 0.

Step 2: Denoising

The solution u + of the TV deconvolution step is then "denoised" using learned sparse coding, cleaning artifacts and enhancing resolution especially along the depth axis. Both dictionary learning and sparse coding are done using the SPAM toolbox 1 [START_REF] Mairal | Online dictionary learning for sparse coding[END_REF][START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF].

Dictionary learning

The dictionary D [s=1] is learned on the set P [1] of all overlapping patches extracted from lateral sections of the 3D deconvolved image u + :

D [s=1] = arg min D∈D min α λ0 |α ,k | + k,m q k,m - D [s=1] ,m α k, 2    , (5)
1 an open source toolbox: http://spams-devel.gforge.inria.fr/index.html where q is a vector of 2 normalized patches:

q k,m = p [s=1] k,m m p [s=1] k,m 2 .

Sparse coding denoising

The learned dictionary is then used to denoise every 2D sections using sparse coding. If the sampling rate differs from one dimension to another, the dictionary is resampled in order to preserve atom physical size along each dimension using interpolation functions f1→2 and f1→3:

D [s=2] = f1→2(D [s=1] ) and D [s=3] = f1→3(D [s=1] ).
The denoised image is

u ++ k = k ,m,s R [s] k,k ,m l D [s] l,m α [s] k ,l (6)
α [s] = arg min α m p [s] k,m - D [s] ,m α [s] ,k 2 +λs α [s] ,k . (7)
R is a reconstruction operator which simply copies a weighted version of every synthetic patch k at its proper position in the 3D image. For instance, if only the central pixels m0 of the reconstructed patches are averaged:

R k,k ,m,s = 1/3 if k = k and m = m0, 0 otherwise, , (8)
as each pixels are at the center of only 1 patch in each of the 3 subspaces.

Step 3: Deconvolution using quadratic prior

As expected, the solution of the step 2 u ++ may be negative and yields to a model (H u ++ ) which is too far from the actual data as attested by the value of J data (u ++ ). As a consequence, a final deconvolution step using u ++ as a quadratic prior was added:

u +++ = arg min u≥0 J data (u) + µ2 u -u ++ 2 2 , (9)

CHOICE OF THE PARAMETERS

For both deconvolution steps (step 1 and 3), few parameters were needed, namely the Gaussian noise variance in each channel, the hyper-parameters µ1 and µ2 and the parameter of the TV hyperbolic approximation.

For the learned sparse coding step, the main parameters are set according to Mairal prescriptions [START_REF] Mairal | Online dictionary learning for sparse coding[END_REF][START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF]:

• size of the patch M [1] = 9 × 9,

• number of elements of the dictionary L = 500,

• hyper-parameter λs = 1/ √ M [s] .

In addition to these parameters, there is however a lot of possible variations of the presented method. For instance any other deconvolution method can be used as step 1, orthogonal matching pursuit can be preferred to the LASSO to solve Eq. ( 6) and/or Eq. ( 5). Such questions should be assessed in further more extensive studies.

RESULTS

Simulations

This method was first tested on synthetic micrographs simulated for "the 3D Deconvolution Microscopy Challenge" held during ISBI 2013. This dataset is composed of four channels with 1024 × 512 × 256 voxels. The four channels of this dataset were deconvolved independently using an aberration free diffraction PSF. The deconvolution results, presented in Fig. 2, show the strong improvement both in term of resolution and denoising. The increase of resolution is particularly spectacular along the depth axis for sharp structures like the filaments in green. With these results, the presented method reached the first rank at the "the 3D Deconvolution Microscopy Challenge". Unfortunately, ground truth and the quantitative figures of merit of each participant of this challenge are still undisclosed.

I have demonstrated the robustness of the presented method by processing experimental data shared by Cazares et al. on the cell image library website2 . This dataset consists on an observation of microtubules in a Drosophila S2 cell labeled Alexa Fluor 488 and observed with a Zeiss Elyra structured illumination microscope (SIM) with a NA = 1.4 objective. It is composed of a 3D widefield image and a super-resolution SIM image that can be used as a ground truth. Both micrographs have (1024 × 1024 × 44) voxels of size 40 × 40 × 110 nm. The PSF, parametrized by its aberrations, was estimated by deconvolving the widefield data by the super-resolution SIM data using the PSF estimation stage of the blind deconvolution algorithm presented in [START_REF] Soulez | Blind deconvolution of 3d data in wide field fluorescence microscopy[END_REF].

The data, the deconvolution result and the structured illumination image are shown on Fig. 3. Compared to widefield, the deconvolution result confirms the effectiveness of my method to increase the resolution especially along the depth dimension. Furthermore, even if much more details can be seen on SIM observation that still have a better resolution, the deconvolution result seems less noisy.

Fig. 2 .

 2 Fig. 2. (x, y) and (x, z) sections of the first "3D Deconvolution Microscopy Challenge" data set compared with its deconvolution using the proposed method.

  y) section (x, z) section

Fig. 3 .

 3 Fig. 3. Widefield observation and its deconvolution compared to the super-resolution image given by structured illumination (SIM) on Drosophila S2 cell.

http://www.cellimagelibrary.org/images/36797


Louis Thibon

Ferréol Soulez

Éric Thiébaut

Fast automatic myopic deconvolution of angiogram sequences *

We present a fast unsupervised myopic deconvolution method dedicated to quasi-real time processing of video sequences such as angiograms. Our method is based on a Bayesian approach of which the tuning parameters are automatically set thanks to the marginalized likelihood of the observed image. We demonstrate the effectiveness of our approach on simulated and empirical images.

CONTEXT AND PREVIOUS WORK

Coronary angiography is a medical imaging technique used to visualize heart vessels (the coronaries) in order to diagnose and prevent potential heart failures. Its principle consists in injecting a radiocontrast agent by catheter into an artery and capturing an image sequence of the cardiovascular system thanks to X-ray irradiation.

In a previous work [START_REF] Soulez | Blind deconvolution of video sequences[END_REF], we shown that the blur of angiogram sequences could be approximated by a convolution by a shift-invariant point spread function (PSF). We found that the actual PSF is constant for a given video sequence but depends on the operating conditions and on the patient. We therefore proposed to use multi-frame blind deconvolution to improve the quality of the videos both in terms of signal to noise ratio and resolution. We observed that the PSF are approximately isotropic and bell shaped with a profile similar to a Lorentzian distribution:

h(r) ≈ η 1 + 2 r/γ 2 , (1)
with r the 2-D position on the detector, γ the full width at halfmaximum (FWHM) of the PSF and η a normalization factor. This multi-frame blind deconvolution method is however lengthy and requires the tuning of a number of control parameters which make it impracticable for non-specialists and incompatible with the quasi real-time requirement. In this paper, we propose several improvements to achieve a fast and unsupervised method which works well in practice.

UNSUPERVISED MYOPIC DECONVOLUTION

Data model

For a given frame of the sequence, a raw image writes [START_REF] Soulez | Blind deconvolution of video sequences[END_REF]:

y = Hγ • x + n (2)
with y ∈ R M the data (yj is the value measured by j-th pixel and M is the number of pixels), Hγ ∈ R M ×N the blur operator, x ∈ R N the object of interest in the form of a pixelized image with N pixels and n ∈ R M a term which accounts for noise and approximation errors. The notation Hγ indicates that the blur depends on some parameters such as the width γ of the PSF.

Bayesian approach

Our goal is to recover the object x given the data y. This is a myopic deconvolution problem as it requires to solve a joint problem of deblurring the data and estimating the PSF parameters. Deconvolution is well known to be an ill-conditioned inverse problem which has to be regularized to obtain a robust and stable solution [START_REF] Titterington | General structure of regularization procedures in image reconstruction[END_REF]. The maximum a posteriori (MAP):

x MAP = arg max x Pr(x | y, θ) (3)
provides such a solution. In words, MAP aims at maximizing the posterior probability of the unknowns x given the data y and socalled hyper-parameters θ (among others the ones which characterize the PSF). Using Bayes rule:

Pr(x | y, θ) = Pr(y | x, θ) Pr(x | θ) Pr(y | θ) ,
and realizing that the evidence Pr(y | θ) does not depend on x yields:

x MAP = arg max x Pr(y | x, θ) Pr(x | θ) , (4)
where the likelihood of the data Pr(y | x, θ) depends on the noise statistics and where the prior distribution Pr(x | θ) implements prior knowledge about the sought image x. This however assumes that the hyper-parameters θ are known. In order to have an unsupervised method, we want to be able to automatically determine θ. To that end, we propose to use maximum likelihood estimator for the hyperparameters:

θ GML = arg max θ Pr(y | θ) . (5)
Using again Bayes rule, Pr(y | θ) is obtained by marginalizing the joint distribution of y and x given θ:

Pr(y | θ) = Pr(y, x | θ) dx = Pr(y | x, θ) Pr(x | θ) dx .
The problem now amounts to integrate the joint distribution. Owing to the number of variables (N the number of pixels in the restored image), the only realistic cases are when the integral can be carried out analytically. This implies to restrict our-self to Gaussian distributions.

Generalized maximum likelihood

We introduce the, possibly improper, Gaussian distribution of the centered variables u with symmetric positive semi-definite weighting matrix W:

G(u; W) def = W 2 π 1/2 + exp - 1 2 u 2 W , (6)
with u 2 W = u •W• u and |W| + the product of the strictly positive eigenvalues of W. The weighting matrix W is also called the precision matrix by some authors [3]. If the covariance matrix of u is non-singular, taking W = Cov(u) -1 and noting that |W| + is then the determinant of W, the usual Gaussian distribution is retrieved. However, the definition in Eq. ( 6) let us cope with improper statistics.

In practice, the statistics of the error term n in Eq. ( 2) is well approximated by a centered Gaussian distribution, the likelihood of the data y knowing the object x and the tuning parameters θ is then given by:

Pr(y | x, θ) = G(y -Hγ • x; B/σ 2 ) , (7)
with Hγ • x = E(y | x, θ) and B/σ 2 = W y|x,θ the expectation and precision matrix of the data given x and θ respectively. The scaling term σ 2 characterizes the noise level. It is, for instance, an estimate of the noise variance after some whitening implemented by the matrix B which is, at least, positive semi-definite and may be different from the identity to account for non-i.i.d. noise. Since we restrict our-self to Gaussian distributions, the most general expression of the a priori distribution of x knowing θ is:

Pr(x | θ) = G(x -x; µ C) , (8)
with x = E(x | θ) and µ C = W x|θ the expected value and the precision matrix of x knowing θ respectively. The matrix C is symmetric positive semi-definite and µ > 0 is a tuning parameter. In what follows, in order to factorize the noise level characterized by σ 2 , we introduce α = σ 2 µ as the tuning parameter for the priors.

In the considered case (affine direct model and Gaussian distributions) and providing H γ • B • Hγ + α C is invertible, the MAP solution is unique and has a closed form:

x MAP = x + R θ • (y -Hγ • x) , (9)
with R θ the reconstruction matrix given by:

R θ = (H γ • B • Hγ + α C) -1 • H γ • B , (10)
which does not depend on the noise level σ 2 but solely on:

θ = {α, γ} (11)
which are thus the tuning parameters of our myopic problem: γ let us tune the PSF, while α let us set the relative weight of the priors. Providing that B/σ 2 = W y|x,θ does not depend on x, it can be easily shown that the evidence is also a Gaussian distribution:

Pr(y | θ) = G(y -Hγ • x; Q θ /σ 2 ) , (12)
with Hγ • x = E(y | θ) the expectation of y knowing θ and Q θ /σ 2 = W y|θ the corresponding precision matrix with:

Q θ = B • (I -A θ ) , ( 13
)
where I is the identity and A θ def = Hγ • R θ is the so-called influence matrix [START_REF] Craven | Smoothing noisy data with spline functions[END_REF]. Now the evidence Pr(y | θ) can be maximized with respect to σ 2 to yield the generalized maximum likelihood (GML) estimate of the noise variance:

σ 2 GML = 1 M+ y -Hγ • x 2 Q θ (14)
with M+ the number of strictly positive eigenvalues of Q θ . Replacing the noise variance by its maximum likelihood value σ 2 GML , maximizing the Pr(y | θ) with respect to the other tuning parameters amounts to minimize:

GML(θ) = y -Hγ • x 2 Q θ M+ |Q θ | 1/M + + (15)
which is a generalization of the criterion considered by Wahba [5] and Trouvé et al. [3] as it takes into account non-i.i.d. noise and the a priori expectation of the parameters x = E(x | θ).

Practical formulation

In our case, the noise is approximately uniform and uncorrelated but we have to cope with unmeasured pixels (see leftmost image in Fig. 1) so B is the diagonal matrix B = diag(b) with b the mask of measured pixels:

bj = 1 if j-th pixel is valid 0 else. (16)
In order to regularize the inverse problem, we impose the smoothness of the object via its prior distribution as it is typically done in image reconstruction with x = 0 and:

C = D • D , ( 17
)
where D is a finite difference operator. That is to say that maximizing Pr(x | θ) amounts to minimizing

x 2 C = D • x 2
2 which is a quadratic measure of the roughness of x. Using 2-D indices for the sake of clarity, our regularization typically writes:

x 2 C = j 1 ,j 2 (xj 1 +1,j 2 -xj 1 ,j 2 ) 2 + j 1 ,j 2 (xj 1 ,j 2 +1 -xj 1 ,j 2 ) 2 ,
(18) with N1 and N2 the restored image dimensions, thus N = N1 × N2 is the number of restored pixels.

Myopic estimation

To summarize, we propose to solve the myopic part of our problem by minimizing GML(θ) with respect to θ. This is very like the strategy proposed by Blanc et al. [START_REF] Blanc | Marginal estimation of aberrations and image restoration by use of phase diversity[END_REF]. Taking into account that x = 0, the GML criterion is finally given by:

GML(θ) = y • B • (I -A θ ) • y M+ |B • (I -A θ )| 1/M + + , (19)
with M+ the number of strictly positive eigenvalues of B • (I-A θ ).

For comparison, we also consider the generalized cross validation criterion [START_REF] Craven | Smoothing noisy data with spline functions[END_REF] which has been extensively used as the method of choice for choosing the best hyper-parameters:

GCV(θ) = y • (I -A θ ) • B • (I -A θ ) • y M+ [tr(I -A θ )/M+] 2 . ( 20
)
Fig. 1: Blind and myopic deconvolution of a single coronarographic image. Left: original image. Center: blind deconvolution by Soulez et al. [START_REF] Soulez | Blind deconvolution of video sequences[END_REF]. Right: proposed myopic deconvolution (here with a circular mask).

Circulant approximation

The criteria (GML or GCV) which have to be minimized to determine the tuning parameters θ are not in general applicable in large dimension because they involve to compute the determinant or the trace of a very large matrix. There exist means to approximate the quantities for large matrix [7], but our experience is that the achieved precision is insufficient for our purposes. Since H is a convolution, it can be approximated by a circulant block circulant matrix which is diagonalized by the discrete Fourier transform (DFT):

Hγ ≈ F -1 • diag( ĥ) • F , (21)
with F the DFT operator and ĥ = F • h the DFT of the shift invariant PSF h. Note that the corresponding adjoint of the PSF operator is:

H γ ≈ F -1 • diag( ĥ ) • F.
The same trick can be applied to the regularization operator C = D • D which is also shift-invariant:

C ≈ F -1 • diag(ĉ) • F . ( 22
)
For instance, the regularization term in Eq. ( 18) is approximated by:

x 2 C ≈ k1,k 2 ĉk 1 ,k 2 |x k 1 ,k 2 | 2 ,
with x = F • x the DFT of x and:

ĉk 1 ,k 2 = 4 sin 2 (π k1/N1) + 4 sin 2 (π k2/N2) , (23)
which is real and non-negative. Finally, by only considering a subimage where all pixels are measured:

B ≈ I . ( 24
)
Putting all together, we obtain the following approximate expressions in Fourier domain for the generalized maximum likelihood and generalized cross-validation criteria:

GML(θ) ≈ k qk |ŷ k | 2 k∈K + qk 1/M + , (25)
GCV(θ) ≈ N k q2 k |ŷ k | 2 k qk 2 , (26)
with ĥ = F • y the DFT of the raw image y and:

qk = α ĉk | ĥk | 2 + α ĉk , (27)
and where K+ is the set of frequels for which qk = 0 and M+ is the size of this set.

Algorithm summary

To process a sequence of images, we first estimate the hyperparameters on a sub-region of a selected image of the sequence by minimizing GML(θ) under circulant approximation. To carry out this optimization, we use Nelder-Mead simplex algorithm [8] as described in [9] with simple modifications to implement bound constraints on the parameters θ. Given the GML estimate of the hyper-parameters, we perform a regularized deconvolution for each image of the sequence. This amount to solve the linear equations:

(H γ • B • Hγ + α C) • x = H γ • B • y . (28)
At least due to the non-stationary operator B, these equations cannot be easily diagonalized by FFT, so we use an iterative method such as the linear conjugate gradients [10].

RESULTS

We tested our method on simulated images and on a set of empirical angiogram sequences provided by the cardiac center of Hôpital de la Croix Rousse (Lyon, France). Figures 1 and2 show typical results for empirical and simulated images respectively.

Validation of our approximations

In practice, the tuning parameters obtained from the minimization with respect to α and γ of the two versions of the GML criterion, the exact one in Eq. ( 19) and the approximated one in Eq. ( 25), are very similar. Following Trouvé et al.

[3], we used the generalized singular value decomposition (GSVD) of the operators H and D 29), but with GML parameter estimation.

to speed-up the computations of the exact criterion. The DFT approximation is however much faster1 , so we favor the use of our DFT-based method.

Comparison with other approaches

When the problem was solely to estimate the regularization level α, Wahba [5] found that GML and GCV give equally good results. However, in our myopic context where both α and the width γ of the PSF are unknown, we observed that GML provides much better estimates of the tuning parameters than GCV. This confirms the conclusions of [3] who also favored GML rather than GCV to estimate the depth from defocus of real images.

Compared to the original method by Soulez et al. [START_REF] Soulez | Blind deconvolution of video sequences[END_REF], our method gives similar results in terms of quality (see Fig. 1) in spite of the simplifications made (myopic instead of blind deconvolution and quadratic regularization instead of total variation). The myopic approach is however much faster: on an Intel Core i7 CPU at 2.94 GHz, it takes about 0.4 second per 512 × 512 image instead of 80 seconds for the blind method. Including the myopic estimation step, a typical angiogram sequence of 80 images is processed in only 30 seconds. Our myopic algorithm is therefore a quasi real-time unsupervised method.

Finally, we also compared our strategy to the deterministic approach of [11]. They considered a quadratic regularization such that:

ĉk 1 ,k 2 = 1 1 + k 2 1 + k 2 2 /k0 p , (29)
where k0 is the cut-off frequency (in frequels) of the object and p is the exponent of the decreasing rate. Figure 2 shows that their regularization which requires to fit 2 more parameters (k0 and p) and have therefore more flexibility to adapt to the actual object structure yields an improved restoration compared to our simple regularization, given in Eq. (18). However, this is only true if parameter estimation is done via GML, not by fitting the power spectrum of the data as suggested by [11].

CONCLUSION

We have presented a fast unsupervised algorithm for myopic deconvolution of images. Our algorithm was tested on simulations and empirical sequences of medical images (angiograms). The first stage of the method consists in estimating the tuning parameters by means of the generalized maximum likelihood (GML) while the second stage is a regularized deconvolution. To estimate quickly the GML criterion, we rely on circulant approximations of the operators which are thus diagonalizable by discrete Fourier transform (DFT). We have validated our assumptions on simulations and by comparing the approximated criterion to the exact one.

Our algorithm is currently under testing at the Hôpital de la Croix Rousse (Lyon, France). Processing a typical coronarographic sequence of 80 images with 512 × 512 pixels takes only 30 seconds. This is about 4 times the observing duration. Thus our algorithm can be used right after the acquisition to quickly enhance the quality of the image and help the physician to make an improved diagnostic. A graphical interface have been developed and the software will be made freely available at http://mitiv.univ-lyon1.fr.

Fig. 2 :

 2 Fig. 2: Myopic deconvolution of an image blurred by a Gaussian with FWHM γ = 2 pixels and with SNR = 8 dB. A: result with finite difference regularization, cf. Eq. (23), and GML parameter estimation. B: result following method in [11]. C: result with same regularization as in B, cf. Eq. (29), but with GML parameter estimation.

Solving the exact myopic problem by GSVD can take up to

hours of computation on a 16 × 16 region of interest; compared to a few milliseconds when using the DFT, even on a larger, say 64 × 64, region.

This work is part of the MiTiV project funded by the French National Research Agency (ANR DEFI 09-EMER-008-01).


Xiaojun Bao

Hongfei Zhang
email: zhanghongfei@lzu.edu.cn

Haifei Zhang

G Royer

Junqing Li

October

Systematical calculation of α decay half-lives with a generalized liquid drop model

Keywords: Preformation factor, α decay, Half-lives, Shell effects

A systematic calculation of α decay half-lives is presented for even-even nuclei between Te and Z=118 isotopes. The potential energy governing α decay has been determined within a liquid drop model including proximity effects between the α particle and the daughter nucleus and adjusted to reproduce the experimental Q value. The α decay half-lives have been deduced from the WKB barrier penetration probability. The α decay half-lives obtained agree reasonably well with the experimental data.

Introduction

The α decay is one of the most important decay channels of the heavy and superheavy nuclei.

Measurements on the α decay can provide reliable information on the nuclear structure such as the ground state energy, the ground state half-life, the nuclear spin and parity, the nuclear deformation, the nuclear clustering, the shell effects, and the nuclear interaction [1][2][3][4].

Experimentally α decay of nuclei is used to identify the new nuclides and new elements through α decay chain from unknown parent nucleus to a known nuclide [5]. Recently, the interest in the α decay has been renewed because of the development of radioactive beams and new detector technology under low temperature. Some newly synthesized superheavy elements have recently been identified using α decay [6][7][8][9].

The process of α decay is fundamentally a quantum tunneling effect, which was first explained by Gamow and by Condon and Guerney in the 1920s [10,11]. Later on, theoretical calculations were performed to predict the absolute α decay width, to extract nuclear structure information, and to pursue a microscopic understanding of the α decay phenomenon.

These studies are based on various theoretical models such as the shell-model, fission-like model, and cluster-model [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27]. The simple empirical relations between α decay half-lives and decay energies are also discussed [28][29][30][31]. Generally tunneling penetration is used to describe the α decay, in which the penetration probability was calculated using WKB approximation assuming α particle tunneling through the potential barrier between α cluster and the daughter nucleus in the parent nucleus. In the unified fission approach [32,33] the decay constant λ is simply the product of the barrier penetrability P and of a constant assault frequency ν 0 . Then, the height, position and width of the potential barriers are the main ingredients determining the half-lives. In the cluster-model [12,34], the cluster is assumed to form before it penetrates the barrier and a preformation factor is included in the calculation. The decay constant λ is defined as the product of the preformation factor, the assault frequency and the penetration probability. Usually, computing the α formation amplitude is a difficult task because the actual wave functions involved cannot be well defined.

The α preformation factor is very important from the viewpoint of the nuclear structure.

Numerous studies of the α decay have been concentrated on this problem [35,36]. In recent studies, the α preformation factor is extracted by dividing the experimental α decay width by the barrier penetration probability, which can easily be obtained from the WKB approximation [25,37].

The generalized liquid drop model (GLDM) was employed to calculate the nuclear potential. It has been successfully used to obtain the potential barriers for fusion reactions [38] and nuclear decays [39][40][41]. In our previous calculations [16], within a superasymmetric fission picture allows us to reproduce the experimental α decay half-lives when the experimental Q α values are used. The α decay half-lives have been calculated using the GLDM [16], in which the assault frequency ν 0 is fixed as 1.0 × 10 20 s -1 . In this work, the α decay halflives are calculated in the preformed cluster-model and the assault frequencies ν 0 have been calculated using a classical method. This paper is organized in the following way. In Section 2 we present the framework of the Generalized Liquid Drop Model. The numerical results are presented and discussed in Section 3. Finally, a brief summary of the present work is given in Section 4.

Theoretical framework

The α decay constant is defined as, λ = P 0 ν 0 P.

(

The assault frequency is calculated using the classical method,

ν 0 = 1 2R 2E α M α , (2.2)
where R is the radius of the parent nucleus and E α is the energy of the α particle, corrected for recoil; M α being its mass.

The penetration probability P is calculated within the WKB approximation. The GLDM energy of a deformed nucleus is defined as [38]:

E = E V + E S + E C + E prox , ( 2.3)
where the different terms are respectively the volume, surface, Coulomb, nuclear proximity energies.

For one-body shapes, the volume E V , surface E S and Coulomb E C energies are given by

E V = -15.494(1 -1.8I 2 )A M eV, ( 2.4
)

E S = 17.9439(1 -2.6I 2 )A 2/3 (S/4πR 2 0 ) M eV, ( 2.5
)

E C = 0.6e 2 (Z 2 /R 0 )B C M eV.
(2.6)

B C is the Coulomb shape dependent function, S is the surface of the one body deformed nucleus and I is the relative neutron excess.

B C = 0.5 (V (θ)/V 0 )(R(θ)/R 0 ) 3 sin θdθ, ( 2.7)
where V (θ) is the electrostatic potential at the surface and V 0 the surface potential of the sphere. The effective sharp radius R 0 has been chosen as

R 0 = 1.28A 1/3 -0.76 + 0.8A -1/3 fm. (2.8)
This formula proposed in Ref. [42] is derived from the droplet model and the proximity energy and simulates rather a central radius for which R 0 /A 1/3 increases slightly with the mass. It has been shown that this selected more elaborated expression can also be used to reproduce accurately the fusion, fission and cluster and alpha decay data. To ensure volume conservation, the radii R 1 and R 2 of the daughter and α nuclei are given by

R 1 = R 0 (1 + β 3 ) -1/3 , ( 2.9
)

R 2 = R 0 β(1 + β 3 ) -1/3 , ( 2.10)
where,

β = 1.28A 1/3 1 -0.76 + 0.8A -1/3 1 1.28A 1/3 2 -0.76 + 0.8A -1/3 2 . (2.11)
When the fragments are separated [43],

E V = -15.494[(1 -1.8I 2 1 )A 1 + (1 -1.8I 2 2 )A 2 ] MeV, (2.12
)

E S = 17.9439[(1 -2.6I 2 1 )A 2/3 1 + (1 -2.6I 2 2 )A 2/3 2 ] MeV, (2.13
)

E C = 0.6e 2 Z 2 1 /R 1 + 0.6e 2 (Z 2 2 /R 2 ) + e 2 Z 1 Z 2 /r MeV. (2.14)
The surface energy comes from the effects of the surface tension forces in a half space. When a neck or a gap appears between separated fragments an additional term called proximity energy must be added to take into account the effects of the nuclear forces between the close surface. It moves the barrier top to an external position and strongly decreases the pure Coulomb barrier:

E prox (r) = 2γ hmax h min Φ[D(r, h)/b]2πhdh, ( 2.15)
where

γ = 0.9517 (1 -2.6I 2 1 )(1 -2.6I 2 2 ) M eV f m -2 , (2.16)
r is the distance between the mass centres, h is the transverse distance varying from the neck radius or zero to the height of the neck border, D is the distance between the opposite surfaces in consideration and b is the surface width fixed at 0.99 fm. Φ is the proximity function. The surface parameter γ is the geometric mean between the surface parameters of the two fragments.

The barrier penetrability P is calculated within the action integral

P = exp[- 2 Rout R in 2B(r)(E(r) -E(sphere))dr], (2.17)
with

E(R in ) = E(R out ) = Q exp . B(r) = µ, in which µ is the reduced mass.
To calculate the absolute α decay width, the α preformation factor P α is indispensable based on the Gamow picture, which measures the probability that an α cluster is present in the decaying nucleus. Within a superasymmetric fission picture the preformation factor P 0 has been taken as 1 in previous studies and that allows us to reproduce the experimental α decay half-lives when the experimental Q α values are used. However, there are still small differences between the calculated and experimental values and these discrepancies may be used to determine the α preformation probability. The microscopic calculation gives a value of 0.3 for the α cluster preformation factor of even-even nucleus 212 Po [12]. However, development on the microscopical description of the α cluster preformation factor is still slow due to the complexity of the nuclear many-body problem. Experiments have shown that the preformation factor varies smoothly in the open shell region and has a value smaller than 1.0 [1]. As a result, it is reasonable and appropriate to take the preformation factor as a constant for all even-even nuclei. This means both the medium mass α emitters and heavy ones can be well described in a consistent way by the GLDM. We fix the value of the preformation factors P α =0.38 for even-even nuclei, which is consistent with both the microscopic calculations and the experimental data of open shell nuclei [1,12]. Finally, the partial half-life is related to the decay constant λ by T 1/2 = ln 2/λ. The half-lives of α decay from the ground state to ground state for even-even nuclei with proton number Z=52-118 have been calculated within the GLDM. The ground state spin and parity of all even-even nuclei is 0 + . This means that the α decay of even-even nuclei mainly proceeds to the ground state of the daughter nucleus. Actually, the parent nucleus can also decay to the excited states of the daughter nucleus, this probability is very small in normal cases, and it can be neglected for a systematic calculation of half-lives.

Results and discussions

The numerical results for the even-even nuclei from Te to 118 are listed in Table 1. The first and second columns denote the parent nucleus and neutron number N, respectively.

The third and forth columns are respectively the experimental decay energies and partial half-lives of α decay. The calculated half-lives are listed in the last column. The meanings of columns 6-10 are similar to those of columns 1-5. It is known experimentally that the magnitude of α decay half-lives of the even-even nuclei varies in a very wide range from 10 -7 to 10 22 s. Although the amplitude of the variation of half-lives is as high as 10 29 times, we can see from Table 1 that the experimental α decay half-lives of many even-even nuclei are reproduced within a factor of 3 by GLDM. But for the four regions the vicinity of N=50, N=126, N=152 and superheavy region the deviations between experimental and calculated values are relatively large. The big deviation occurring for Te isotope chains are mainly due to the shell effect 100 Sn nucleus. In the superheavy nuclei region, one can see from the Table 1 that the calculated values deviate obviously from the experimental ones. This may be understood as a consequence of the following two sides. On the theoretical side, the superheavy nuclei α preformation factor is smaller than the one of medium and heavy nuclei [25] in our previous study. Theoretical calculation should give better agreement with the experimental data if the preformation factor is considered as a variable with different parent nuclei. On the experimental side, it would be interesting to improve the precision of experimental α decay half-lives. In order to show the systematic behavior of the agreement between model and data clearly, we also plot the comparison of experimental half-lives and theoretical ones in Fig. 2.

The theoretical points (stars) almost coincide with the experimental ones (black circles).

First, Po, Rn, Ra and Th isotope chains half-lives are show in Fig. 2 as a function of neutron number N. On the one hand, the half-lives increase with increasing neutron number up to the spherical shell closure N=126, and then decrease rapidly with the neutron number.

The minimum half-lives appears at N=128 for all the four selected isotopes, indicating α emission is much easier for the nuclei with two neutrons outside the shell closure. Moreover, the values of the minimum half-lives are nearly identical at N=128, after which the half-lives increase again with increasing N. On the other hand, in Fig. 2 The penetration probabilities and the Q α values for Cm, Cf, Fm, No and Rf isotopes are also shown in Fig. 3. We can see from this figure that the change in Q α is similar to that in penetration probability. They decrease with increasing neutron number up to N=152, then increase with N up to N=154. After that, they decrease again. It is well known that the shell effects for α radioactivity are related to the Q α value, which is maximum when the daughter nuclei has a magic number of neutrons and protons. The changes in Q α and P are similar. The lower is Q α the more difficult it will be for an α cluster to penetrate the potential barrier in the parent nucleus in accordance with quantum tunneling penetration theory.

The calculated half-lives are compared with the experimental ones for isotopes of Th, U, Pu, Cm, Cf, Fm, No and Rf in Fig. 4. We can see from this figure that the experimental data are well reproduced by GLDM in these isotopes considered here. The theoretical points (stars) almost coincide with the experimental ones (black circles). Fig. 4 shows that the αdecay half-lives decrease as the value of Z increases, which reflects the stability gained by these nuclei when two protons are removed. The half-lives do not increase monotonously with N and a maximum exists at N=152 due to the shell effects. We may conclude that N=152 is a deformed magic number. Around the shell N=152, the variation of half-lives is approximately 10 times for Cm, Cf, Fm and No isotopic chains. For instance, the αdecay half-life of 250 Cf (N=152) is only ten times larger than that of 252 Cf (N=154). This demonstrates that the influence of the deformed shell on half-lives is less than that of the spherical one.

Summary

A global calculation of the α decay half-lives are presented for the even-even nuclei from Z=52 to Z=118, in which the penetration probability was calculated using WKB approximation assuming α particle tunneling through the potential barrier between α cluster and the daughter nucleus in the parent nucleus. The barriers are constructed with the GLDM and the decay energies used in our calculations are extracted from experimental data. The assault frequency ν 0 has been calculated using a classical method. The value of the preformation factors P α has been fixed to 0.38 for even-even nuclei. Finally, the decay constant is the product of the preformation, the penetration, and the assault frequency. The obtained decay half-lives agree reasonably well with the available experimental data. The change in α decay half-lives with neutron number shows that shell effects play an important role in the behavior of α decay half-lives around magic numbers. It is worth noting that the aim of this work is not only to reproduce the experimental data, but also to extend our understanding of α decay half-lives around shell closures.

Figure 1 :

 1 Figure 1: The penetration probabilities and experimental decay energies Q α of α decay even-even Pb, Po, Rn, Ra and Th isotopes.

Figure 2 :

 2 Figure 2: Comparison of the calculated α decay half-lives with the experimental data for even-even nuclei ranging from Z=82 to Z=90, showing the shell effect at the neutron magic numbers N=126.

Fig. 1

 1 Fig.1shows the penetration probability in the left hand side as a function of neutron number N. The penetration probabilities decrease with increasing neutron number up to the spherical shell closure N = 126, and then increase rapidly with the neutron number. The maximum probability appears at N = 128 for all the four selected isotopes. For the nuclei

  with two neutrons outside the closed shell, the α particle emission is easier than that of the other nuclei of the same isotopes. The closed shell structures play also the key role for the penetrability mechanism. Moreover, the values of the maximum penetration probabilities are nearly identical at N = 128 for Po, Rn, Ra, and Th isotopes, after which the penetration probability decreases again with increasing N. At the same neutron number the penetration probability increases with increasing proton number because proton number of Po isotopes are closer to a proton magic number Z = 82 than those of Rn, Ra, and Th. The closer the nucleon number is to a magic number, the more difficult the α particle penetrates from the parent nuclei. The experimental decay energy Q α is shown in the right hand side also as a function of N in Fig.1. One can see from the figure that Q α decreases with the increasing neutron number N for the Pb, Po, Rn, Ra, and Th isotopes before the neutron number N=126. The Q α values increase sharply after N=126 and reach the maximum at N=128, then they decrease again quickly for Po, Rn, Ra and Th isotopes. The changes in Q α and penetration probability P are similar.

  there is a clear decrease in the half-life before the neutron number N=128, shown in the Po, Rn, Ra and Th isotopes. This is attributed to the strong N=126 shell effect: the main effect of the N=126 shell is included in the decay Q value, which is closely related to the nuclear structure, and the remaining effect is largely absorbed into the α preformation factor. Because the constant preformation factor cannot completely describe the detailed features of nuclear structure, the strong shell effects are clearly shown from the increased deviations in the neighborhood of N=126. At the same neutron number, half-lives decrease with increasing proton number because proton number of Po isotopes are closer to a proton magic number Z=82 than those of Rn, Ra and Th. Both N=126 and Z=82 are well known as magic numbers. The closer the nucleon number is to shell closure, the more difficultly the α cluster forms in the parent nuclei. The dramatic change of the half-lives around the magic number indicates that the shell effects play an important role on α formation mechanism in the parent nuclei. There is no doubt that the theoretical calculation should give better agreement with the experimental data if the preformation factor is considered as a variable with different parent nuclei instead of a constant, especially for the closed shell region nuclei.

Figure 3 :

 3 Figure 3: The penetration probabilities and experimental released energies Q α of α decay even-even Cm, Cf, Fm, No and Rf isotopes.

Figure 4 :

 4 Figure 4: Comparison of the calculated α decay half-lives with the experimental data for even-even nuclei ranging from Z=90 to Z=104, showing the shell effect at the neutron magic number N=152.

Table 1 :

 1 Comparison between experimental and theoretical α decay half-lives and of the even-even nuclei with proton number Z=52-118. The units of the α decay energies and half-lives are MeV and second respectively, and the experimental results are taken from the recently data[44,45].

	Nucleus N	Q	T 1/2 (exp.)	T 1/2 (cal.) Nucleus N	Q	T 1/2 (exp.)	T 1/2 (cal.)
	106 52 Te	54 4.290 8.0 × 10 -5 65.0 × 10 -5 162 76 Os	86 6.767 2.1 × 10 -3 3.6 × 10 -3
	108 52 Te	56 3.420 4.3 × 10 0	28.5 × 10 0 164 76 Os	88 6.479 4.1 × 10 -3 31.4 × 10 -3
	110 52 Te	58 2.699 6.2 × 10 5	99.0 × 10 5 166 76 Os	90 6.139 3.0 × 10 -1 4.9 × 10 -1
	110 54 Xe	56 3.875 3.9 × 10 0	0.87 × 10 0 168 76 Os	92 5.816 4.9 × 10 0	9.1 × 10 0
	112 54 Xe	58 3.330 3.0 × 10 2	20.3 × 10 2 170 76 Os	94 5.537 7.8 × 10 1	13.9 × 10 1
	146 62 Sm	84 2.528 2.2 × 10 15 12.7 × 10 15 172 76 Os	96 5.224 1.7 × 10 3	3.9 × 10 3
	148 64 Gd	84 3.271 2.2 × 10 9	6.4 × 10 9	174 76 Os	98 4.870 1.8 × 10 5	2.8 × 10 5
	150 64 Gd	86 2.808 5.6 × 10 13 22.0 × 10 13 186 76 Os	110 2.820 6.3 × 10 22 10.4 × 10 22
	150 66 Dy	84 4.351 1.2 × 10 3	2.4 × 10 3	166 78 Pt	88 7.286 3.0 × 10 -4 4.5 × 10 -4
	152 66 Dy	86 3.726 8.6 × 10 6	25.3 × 10 6 168 78 Pt	90 6.990 2.0 × 10 -3 3.3 × 10 -3
	154 66 Dy	88 2.945 9.5 × 10 13 19.6 × 10 13 170 78 Pt	92 6.707 6.0 × 10 -3 26.4 × 10 -3
	152 68 Er	84 4.934 1.1 × 10 1	2.1 × 10 1	172 78 Pt	94 6.464 1.0 × 10 -1 1.6 × 10 -1
	154 68 Er	86 4.279 4.8 × 10 4	7.4 × 10 4	174 78 Pt	96 6.183 1.2 × 10 0	1.7 × 10 0
	156 68 Er	88 3.483 6.7 × 10 9	38.5 × 10 9 176 78 Pt	98 5.885 1.6 × 10 1	2.6 × 10 1
	154 70 Yb	84 5.474 4.4 × 10 -1 6.9 × 10 -1 178 78 Pt	100 5.573 2.7 × 10 2	5.8 × 10 2
	156 70 Yb	86 4.811 2.6 × 10 2	8.7 × 10 2	180 78 Pt	102 5.240 1.9 × 10 4	2.3 × 10 4
	158 70 Yb	88 4.170 4.3 × 10 6	4.5 × 10 6	182 78 Pt	104 4.951 4.2 × 10 5	7.9 × 10 5
	156 72 Hf	84 6.028 2.4 × 10 -2 3.5 × 10 -2 184 78 Pt	106 4.629 5.9 × 10 7	10.0 × 10 7
	158 72 Hf	86 5.404 6.4 × 10 0	10.9 × 10 0 186 78 Pt	108 4.598 5.3 × 10 9	7.1 × 10 9
	160 72 Hf	88 4.902 1.9 × 10 3	2.7 × 10 3	188 78 Pt	110 4.320 3.3 × 10 12 1.6 × 10 12
	162 72 Hf	90 4.416 4.9 × 10 5	14.8 × 10 5 190 78 Pt	112 3.252 2.0 × 10 19 1.8 × 10 19
	158 74 W	84 6.613 1.3 × 10 -3 2.0 × 10 -3 172 80 Hg	92 7.524 2.3 × 10 -4 3.9 × 10 -4
	160 74 W	86 6.065 1.1 × 10 -1 1.6 × 10 -1 174 80 Hg	94 7.233 2.0 × 10 -3 2.6 × 10 -3
	162 74 W	88 5.677 3.0 × 10 0	5.9 × 10 0	176 80 Hg	96 6.899 2.3 × 10 -2 2.9 × 10 -2
	164 74 W	90 5.279 1.7 × 10 2	3.1 × 10 2	178 80 Hg	98 6.577 5.0 × 10 -1 3.5 × 10 -1
	166 74 W	92 4.856 5.5 × 10 4	4.1 × 10 4	180 80 Hg	100 6.258 5.4 × 10 0	5.2 × 10 0
	168 74 W	94 4.500 1.6 × 10 6	4.6 × 10 6	182 80 Hg	102 5.996 7.8 × 10 1	5.6 × 10 1

4 × 10 -3 0.14 × 10 -3

Acknowledgements

The work is supported by the Natural Science Foundation of China (Grants 11175074, 11105035 and 10975064), the Fundamental Research Funds for the Central Universities (grants lzujbky-2012-5), by the CAS Knowledge Innovation Project NO.KJCX-SYW-N02.