Jalal M Fadili

Jean-Luc Starck

B Zhang

Seth W Digel

M J Fadili

J.-L Starck

Poisson noise removal in multivariate count data

Keywords: Poisson noise, multivariate data, wavelets, multi-scale variance stabilization

come

INTRODUCTION

A host of estimation methods have been proposed in the literature for non-parametric Poisson noise removal. Major contributions consist of variance stabilization: this is a classical approach which consists in preprocessing the data by applying a variance stabilizing transform (VST) such as the Anscombe transform. [START_REF] Anscombe | The transformation of Poisson, binomial and negative-binomial data[END_REF] It can be shown that the transformed data are approximately stationary, independent and Gaussian. Once we are brought to the Gaussian denoising problem, standard approaches are used before the VST is inverted to get the final estimate. Haar-Fisz transform is another widely used VST, [START_REF] Fryźlewicz | A Haar-Fisz algorithm for Poisson intensity estimation[END_REF][START_REF] Fryzlewicz | GOES-8 X-ray sensor variance stabilization using the multiscale data-driven Haar-Fisz transform[END_REF] which combines the Fisz transform 5 within the Haar transform. Jansen 6 introduced a conditional variance stabilization (CVS) approach which can be applied in any wavelet domain resulting in stabilized coefficients. However, these transformations are only valid for a sufficiently large number of counts per pixel (and of course, for a larger still number of counts, the Poisson distribution becomes Gaussian with equal mean and variance).

Hypothesis testing has also been proposed for Poisson noise removal. Kolaczyk 7 first introduced a Haar domain threshold, which implements a hypothesis testing procedure based on a user-specified false positive rate (FPR). Zhang et al. [START_REF] Zhang | Fast Poisson Noise Removal by Biorthogonal Haar Domain Hypothesis Testing[END_REF] extended the hypothesis testing to the biorthogonal Haar domain, leading to a more regular reconstruction for smooth intensities. Bijaoui and Jammal 9 derived the probability density function (pdf ) of any wavelet coefficient allowing hypothesis testing in a general wavelet domain, but is computationally complex than Haar-based methods due to the lack of closed-form pdf expressions. Kolaczyk 10 proposed "corrected" versions of the usual Gaussian-based thresholds for Poisson data. However, the asymptotic approximation adopted by [START_REF] Kolaczyk | Wavelet shrinkage estimation of certain Poisson intensity signals using corrected thresholds[END_REF] may not allow reasonable solutions in low-count situations.

Empirical Bayesian estimators for Poisson noise removal were studied in. [START_REF] Jansen | Multiscale Poisson data smoothing[END_REF][START_REF] Kolaczyk | Bayesian multiscale models for Poisson processes[END_REF][START_REF] Timmermann | Multiscale Modeling and Estimation of Poisson Processes with Application to Photon-Limited Imaging[END_REF][START_REF] Nowak | A statistical multiscale framework for Poisson inverse problems[END_REF] The low-intensity case apart, Bayesian approaches generally outperform the direct wavelet filtering. [START_REF] Nowak | Wavelet-Domain Filtering for Photon Imaging Systems[END_REF][START_REF] Antoniadis | Wavelet shrinkage for natural exponential families with quadratic variance functions[END_REF] Poisson denoising has also been formulated as a penalized maximum likelihood estimation problem [START_REF] Sardy | Automatic smoothing with wavelets for a wide class of distributions[END_REF][START_REF] Willett | Fast, Near-Optimal, Multiresolution Estimation of Poisson Signals and Images[END_REF][START_REF] Willett | Multiscale Analysis of Photon-Limited Astronomical Images[END_REF] within wavelet, wedgelet and platelet dictionaries.

In a recent paper, Zhang et al. [START_REF] Zhang | Wavelets, ridgelets and curvelets for poisson noise removal[END_REF] have proposed to combine variance stabilization and multiscale transforms, leading to the Multi-Scale Variance Stabilization Transform (MSVST). MSVST can be seen as a generalization of the Anscombe transform; see section 2 for details. This new approach is fast and easy to implement, and more importantly, works very well at very low count situations, down to 0.1 photons per pixel.

This paper

In this paper, we propose a multivalued extension of MSVST, which allows us to remove the Poisson noise in multivariate data sets, when the third dimension is not a spatial dimension, but the wavelength, the energy or the time. Such 3D data are called 2D-1D data sets in the sequel. We show that this multivalued extension can be used advantageously for approximately Gaussianizing and stabilizing the variance of a sequence of independent Poisson random vectors. This approach is shown to be fast and very well adapted to low-count situations. We use a hypothesis testing framework in the wavelet domain to denoise the Gaussianized and stabilized coefficients, and then apply an iterative reconstruction algorithm to recover the estimated vector field of intensities underlying the Poisson data. Our approach is illustrated for the detection and characterization of astrophysical sources of high-energy gamma rays, using realistic simulated observations. We show that the multivariate MSVST permits efficient estimation across the time/wavelength/energy dimension and recovery of spectral properties. Section 2 reviews the MSVST method relative to the isotropic undecimated wavelet transform and section 3 shows how it can be extended to the 2D-1D case. Section 4 presents some experiments on realistic simulated astronomical data. Conclusions are drawn in section 5.

THE 2D MULTISCALE VARIANCE STABILIZATION TRANSFORM (MSVST)

In this section, we review the MSVST method, 1 restricted here to the wavelet transform. Indeed, the MSVST can use other transforms such as the ridgelet or the curvelet transforms; see. [START_REF] Zhang | Wavelets, ridgelets and curvelets for poisson noise removal[END_REF]

VST of a filtered Poisson process

Let's start with a simple 1D case. Given X a sequence of n independent Poisson random variables X i , i = 1, • • • , n, each of mean λ i , let Y i = n j=1 h[j]X i-j be the filtered process obtained by convolving the sequence X with a discrete filter h. Y denotes any one of the Y i 's, and

τ k = i (h[i]) k for k = 1, 2, • • • .
If h = δ, then we recover the Anscombe VST 2 of Y i (hence X i ) which acts as if the stabilized data arose from a Gaussian white noise with unit variance, under the assumption that the intensity λ i is large. This is why the Anscombe VST performs poorly in low-count settings. But, if the filter h acts as an "averaging" kernel (more generally a low-pass filter), one can reasonably expect that stabilizing Y i would be more beneficial, since the signal-to-noise ratio measured at the output of h is expected to be higher.

Using a local homogeneity assumption, i.e. λ i-j = λ for all j within the support of h, it has been shown 1 that for a filter h, the transform Z = b sign(Y + c) |Y + c| with b and c defined as

c = 7τ 2 8τ 1 - τ 3 2τ 2 , b = 2 |τ 1 | τ 2 (1)
is a second order accurate variance stabilization transform, with asymptotic unit variance. By second-order accurate, we mean that the error term in the variance of the stabilized variable Z decreases rapidly as O(λ -2 ).

From (1), it is obvious that when h = δ, we obtain the classical Anscombe VST parameters b = 2 and c = 3/8. The authors in 1 have also proved that Z is asymptotically normally distributed with mean b |τ 1 |λ and unit variance. Poisson noise realizations of X, plotted as a function of the intensity λ for both Anscombe 2 (dasheddotted), Haar-Fisz (dashed) [START_REF] Fryźlewicz | A Haar-Fisz algorithm for Poisson intensity estimation[END_REF] and our VST with the 2D B 3 -Spline filter as a low-pass filter h (solid). The asymptotic bounds (dots) (i.e. 1 for the variance and √ λ for the expectation) are also shown. It can be seen that for increasing intensity, E[Z] and Var [Z] approach the theoretical bounds at different rates depending on the VST used. Quantitatively, Poisson variables transformed using the Anscombe VST can be reasonably considered to be unbiased and stabilized for λ 10, using Haar-Fisz for λ 1, and using out VST (after low-pass filtering with the chosen h) for λ 0.1.

MSVST with the Wavelet Transform

The Undecimated Wavelet Transform

The undecimated wavelet transform (UWT) uses an analysis filter bank (h, g) to decompose a signal a 0 into a coefficient set W = {w 1 , . . . , w J , a J }, where w j is the wavelet (detail) coefficients at scale j and a J is the approximation coefficients at the coarsest resolution J. The passage from one resolution to the next one is obtained using the "à trous" algorithm 20

a j+1 [l] = ( h↑(j) ⋆ a j )[l] = k h[k]a j [l + 2 j k], w j+1 [l] = (ḡ ↑(j) ⋆ a j )[l] = k g[k]a j [l + 2 j k] .
(2)

where h ↑(j) [l] = h[l] if l/2 j ∈ Z and 0 otherwise, h[l] = h[-l]
is the time-reversed version of h, and "⋆" denotes discrete circular convolution. The reconstruction is given by

a j [l] = 1 2 ( h↑(j) ⋆ a j+1 )[l] + (g ↑(j) ⋆ w j+1 )[l] .
The filter bank (h, g, h, g) needs to satisfy the so-called exact reconstruction condition [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF][START_REF] Starck | Astronomical Image and Data Analysis[END_REF] .

The so-called Isotropic UWT (IUWT) [START_REF] Starck | The Undecimated Wavelet Decomposition and its Reconstruction[END_REF] uses the filter bank (h, g = δh, h = δ, g = δ) where h is typically a symmetric low-pass filter such as the B 3 -Spline filter. The reconstruction is trivial, i.e., a 0 = a J + J j=1 w j . This algorithm is widely used in astronomical applications [START_REF] Starck | Image Processing and Data Analysis: The Multiscale Approach[END_REF] and biomedical imaging [START_REF] Olivo-Marin | Extraction of spots in biological images using multiscale products[END_REF] to detect isotropic objects.

MSVST with the UWT

Now the VST can be combined with the UWT to stabilize the wavelet coefficients of a standard separable UWT. The VST is applied to the approximation coefficients (a j ) j , leading to the following scheme:

UWT a j = h↑(j-1) ⋆ a j-1 w j = ḡ↑(j-1) ⋆ a j-1 =⇒ MS-VST + UWT a j = h↑(j-1) ⋆ a j-1 w j = ḡ↑(j-1) ⋆ A j-1 (a j-1 ) (3)
where

A j (a j ) = b (j) sign(a j + c (j) ) |a j + c (j) |. (4)
Let us define τ

(j) k = i h (j) [i]
k . Then according to (1), the constants b (j) and c (j) associated to h (j) must be set to

c (j) = 7τ (j) 2 8τ (j) 1 - τ (j) 3 2τ (j) 2 , b (j) = 2 |τ (j) 1 | τ (j) 2
.

(

) 5
The constants b (j) and c (j) only depend on the filter h and the scale level j. They can all be pre-computed once for any given h. Since these constants scale-dependent, so is the VST, hence the name MSVST.

In the case of the isotropic UWT (IUWT), given the structure of the filter bank, the stabilization procedure is given by

IUWT a j = h↑(j-1) ⋆ a j-1 w j = a j-1 -a j =⇒ MSVST + IUWT a j = h↑(j-1) ⋆ a j-1 w j = A j-1 (a j-1 ) -A j (a j ) , (6)
where A j-1 is defined similarly to A j in (4).

Asymptotic Distribution of the Detail Coefficients

The following result shows the asymptotic normality of the MSVST detail coefficients under the null hypothesis H 0 that the underlying true intensity λ is locally homogeneous. The asymptotic normal distribution has an intensity-independent variance which relies solely on the filter h and the current scale.

Theorem 2.1.

(i) MSVST+IUWT: let b (j) := sign(τ

(j) 1 )/ |τ (j) 1 |. If λ is constant within the support of the filter h (j) [k -•], then w j [k] D -→ λ→+∞ N 0, τ (j-1) 2 4τ (j-1) 1 2 + τ (j) 2 4τ (j) 1 2 - h (j-1) , h (j) 2τ (j-1) 1 τ (j) 1 , (7)
where ., . represents the scalar product.

(ii) MSVST+UWT: let b (j) := 2 |τ

(j) 1 |/τ (j)
2 , and h (j) = h↑(j-1) ⋆ • • • ⋆ h↑(1) ⋆ h for j ≥ 1 and h (0) = δ. If λ is constant within the support of the filter (ḡ ↑(j-1) ⋆ h

(j-1) )[k -•], then w j [k] D -→ λ→+∞ N (0, σ 2 j )
, where

σ 2 j = 1 τ (j-1) 2 m,n ḡ↑(j-1) [m]ḡ ↑(j-1) [n] k h (j-1) [k]h (j-1) [k + m -n] (8)
Again, the above values of b (j) , c (j) , τ

and σ j can all be pre-computed once (h, g) has been chosen.

2D-1D MSVST POISSON DENOISING

2D-1D Wavelet Transform

In the previous section, we have seen how a Poisson noise can be removed from 2D image using the wavelet transform and the MSVST. Extension to a qD data sets is straightforward by separable tensor product, and the denoising will perform very well as long as each object belonging to this q-dimensional space is roughly isotropic.

In the case of 3D data where the third dimension is either the time, wavelength or energy, we are clearly not in this situation, and the naive analysis of a 3D wavelet does not make sense. Therefore, we want to analyze the data with a wavelet transform, where the scale along the non-spatial dimension is not connected to the spatial scale. Hence, an ideal wavelet function would be defined by:

ψ(x, y, z) = ψ (xy) (x, y) ψ (z) (z) , (9)
where ψ (xy) is the spatial wavelet and ψ (z) is the wavelet along the non-spatial dimension.

In the following, for the sake of clarity and without loss of generality, we will consider only the IUWT filter bank in the spatial domain. We note j 1 the spatial resolution index (i.e. scale = 2 j1 ), j 2 resolution index along the non-spatial dimension. In order to have a fast algorithm for discrete data, we use wavelet functions associated to filter banks. Hence, our wavelet decomposition consists in applying first a 2D IUWT for each frame k z , and then, for each spatial location (k x , k y ) and for each 2D wavelet scale scale j 1 , to apply a 1D wavelet transform along z on the spatial wavelet coefficients.

Therefore, given the properties of the IUWT filter bank, the 2D-1D IUWT of the input data S is such that

S[k x , k y , k z ] = a J1,J2 [k x , k y , k z ] + J1 j1=1 w j1,J2 [k x , k y , k z ] + J2 j2=1 w J1,j2 [k x , k y , k z ] + J1 j1=1 J2 j2=1 w j1,j2 [k x , k y , k z ] . (10)
In this expression, we distinguish four kinds of coefficients:

• Detail-Detail coefficients (j 1 ≤ J 1 and j 2 ≤ J 2 ):

w j1,j2 [k x , k y , k z ] = (δ -h) ⋆ h (j2-1) ⋆ a j1-1 [k x , k y , .] -h (j2-1) ⋆ a j1 [k x , k y , .] , ( 11
)
where h is the low-pass defining the IUWT filter bank.

• Approximation-Detail coefficients (j 1 = J 1 and j 2 ≤ J 2 ):

w J1,j2 [k x , k y , k z ] = h (j2-1) ⋆ a J1 [k x , k y , .] -h (j2) ⋆ a J1 [k x , k y , .] . (12)
• Detail-Approximation coefficients (j 1 ≤ J 1 and j 2 = J 2 ):

w j1,J2 [k x , k y , k z ] = h (J2) ⋆ a j1-1 [k x , k y , .] -h (J2) ⋆ a j1 [k x , k y , .] . (13)
• Approximation-Approximation coefficients (j 1 = J 1 and j 2 = J 2 ):

a J1,J2 [k x , k y , k z ] = h (J2) ⋆ a J1 [k x , k y , .] . (14)

2D-1D MSVST

Putting all pieces together, we are now ready to plug the MSVST into the 2D-1D IUWT. The idea again is to first stabilize the (low-pass) approximation subbands, and then compute the stabilized detail coefficients. As in section 3.1, we distinguish the previous four types of coefficients for which the 2D-1D MSVST takes the following forms:

• Detail-Detail coefficients (j 1 ≤ J 1 and j 2 ≤ J 2 ):

w j1,j2 [k x , k y , k z ] = (δ -h) ⋆ A j1-1,j2-1 h (j2-1) ⋆ a j1-1 [k x , k y , .] -A j1,j2-1 h (j2-1) ⋆ a j1 [k x , k y , .] . (15)
• Approximation-Detail coefficients (j 1 = J 1 and j 2 ≤ J 2 ):

w J1,j2 [k x , k y , k z ] = A J1,j2-1 h (j2-1) ⋆ a J1 [k x , k y , .] -A J1,j2 h (j2) ⋆ a J1 [k x , k y , .] . (16)
• Detail-Approximation coefficients (j 1 ≤ J 1 and j 2 = J 2 ):

w j1,J2 [k x , k y , k z ] = A j1-1,J2 h (J2) ⋆ a j1-1 [k x , k y , .] -A j1,J2 h (J2) ⋆ a j1 [k x , k y , .] . (17)
• Approximation-Approximation coefficients (j 1 = J 1 and j 2 = J 2 ):

c J1,J2 [k x , k y , k z ] = h (J2) ⋆ a J1 [k x , k y , .] . (18)
Now, all 2D-1D wavelet coefficients w j1,j2 are stabilized. In the same vein as in Theorem 2.1, it can be shown that under the null hypothesis that the 2D-1D intensity Λ is locally homogeneous within the support of the 2D-1D wavelet at scales (j 1 , j 2 ), the noise on the stabilized coefficients is Gaussian with known scale-dependent variance that depends solely on h. Thus one can detect significant detail coefficients by binary hypothesis testing or hard thresholding at a given critical threshold. The appeal of a binary hypothesis testing approach is that it allows quantitative control of statistical significance.

Iterative Reconstruction

Following the detection step, we have to invert the MSVST scheme to reconstruct an estimate of the underlying intensity Λ. This is however not straightforward because there is no explicit reconstruction formula available. The formal reason is that the stabilizing operators A j1,j2 and the convolution operators along (x, y) and z do not commute, even though the filter bank satisfies the exact reconstruction formula. To circumvent this difficulty, we propose to reformulate the reconstruction as a sparsity-promoting convex optimization problem. In the sequel, as the noise on the stabilized coefficients is Gaussian, we assume without loss of generality that it is standardized to a unit variance.

We define the multiresolution support M j1,j2 which is determined by the set of detected significant coefficients at each scale j 1 ≤ J and j 2 ≤ J 2 and location (k x , k y , k z ), i.e.,

M j1,j2 [k x , k y , k z ] = 1 if w j1,j2 [k x , k y , k z ] is significant, 0 otherwise. (19)
We denote W the 2D-1D undecimated wavelet transform described above, R the inverse transform and S the observed Poisson noisy data set. We are seeking an estimate that preserves the significant structures in the original data by reproducing exactly the same coefficients as the wavelet coefficients of the input data S, but only at scales and positions where significant signal has been detected. At other scales and positions, we want the smoothest solution with the lowest budget in terms of wavelet coefficients. Furthermore, as Poisson intensity functions are positive by nature, a positivity constraint is imposed on the estimated intensity data Λ. Therefore our reconstruction is formulated as a constrained sparsity-promoting minimization problem over the transform coefficients w

min w w 1 subject to Mw = MWS and Rw ≥ 0 , (20)
and the intensity estimate Λ is reconstructed as Λ = R ŵ, where ŵ is a minimizer of (20). Note that the set of minimizers is not empty by coercivity. Moreover, this problem is a convex optimization problem which can be cast as a Linear Program (LP) and solved using interior-point methods. However, the computational complexity of the LP solver increases dramatically with the size of the problem. Here we propose an alternative based on the hybrid steepest descent (HSD). [START_REF] Yamada | The hybrid steepest descent method for the variational inequality problem over the intersection of fixed point sets of nonexpansive mappings[END_REF] The HSD approach allows minimizing convex functionals over the intersection of fixed point sets of nonexpansive mappings. It is much faster than LP, and in our problem, the nonexpansive mappings do have closed forms. Transposed into our context, its main steps can be summarized as follows:

Inputs: noisy data S; a low-pass filter h; multiresolution support M from the detection step; number of iterations N max . Initialization:

w (0) = MWS = Mw S , for t = 1 to N max do 1. wj1,j2 = M j1,j2 w Sj 1 ,j 2 + (1 -M j1,j2 )w (t-1)
j1,j2 , ∀j 1 , j 2 .

2.

w (t) = WP + (R ST βt [ w]),
3. Update the step β t = (N maxt)/(N max -1).

Output: Λ = Rw (Nmax) .

P + is the projector onto the positive orthant, i.e. P + (x) = max(x, 0). ST βt is the soft-thresholding operator with threshold β t . A careful convergence analysis of this algorithm is given in. [START_REF] Zhang | Wavelets, ridgelets and curvelets for poisson noise removal[END_REF] In summary, 2D-1D wavelet Poisson denoising with the MSVST involves the following three main steps:

1. Transformation : Compute the 2D-1D-MSVST using ( 15)-( 18).

2. Detection : Detect significant detail coefficients by hypothesis testing. Here, we take benefit from the asymptotic Gaussianity of the noise in the stabilized data. We can then detect the significant coefficients by classical binary hypothesis testing (corrected or not for multiple testing and dependencies).

3. Reconstruction : Reconstruct the denoised data using the knowledge of the detected coefficients in the algorithm above. We have simulated a 2D-1D data set containing seven columns and five rows of 2D Gaussian-shaped sources. Each source has a different power-law spectrum (i.e. the z-axis is energy). The data set size is 161 × 161 × 31, with an average of 0.032 photons per pixel. As the number of photon counts is very low, for visual convenience, we show in Fig. 2 the 2D image obtained after summing the noisy data set along the z-axis. Fig. 3 shows a comparison between 2D-MSVST denoising of this image, and the image obtained by first applying a 2D-1D-MSVST denoising to the input data, and summing afterward the denoised frames along the z-axis. Fig. 3(a)-(b) show denoising results for the 2D-MSVST with respectively critical threshold values τ = 3 and τ = 5, and Fig. 3(c)-(d) depict the results for the 2D-1D-MSVST using respectively τ = 4 and τ = 6 critical thresholds * . The reason for using a higher threshold level for the 2D-1D data is to correct for multiple hypothesis testings, and to get the same control over global statistical error rates. Roughly speaking, the number of false detections increases with the number of coefficients being tested simultaneously. Therefore, one must correct for multiple comparisons using e.g. the conservative Bonferroni correction or the false discovery rate (FDR) procedure. [START_REF] Benjamini | Controlling the false discovery rate -a practical and powerful approach to multiple testing[END_REF] As the number of coefficients is much higher with the whole 2D-1D data set, the critical detection threshold τ of 2D-1D denoising must be higher to have a false detection rate comparable to the 2D denoising. As we can clearly see from Fig. 3, the results are very close. This means that applying a 2D-1D denoising on the whole data set instead of a 2D denoising on the integrated image preserves the sensitivity of the MSVST. The chief advantage of the 2D-1D-MSVST is the fact that we recover the energy information for each spatial position.

EXPERIMENTAL RESULTS

MSVST-2D-1D versus MSVST-2D

Application to astronomical data

We applied our method to realistic simulated observations with the Large Area Telescope (LAT). The LAT was launched in June 2008 by NASA on the Fermi Gamma-ray Space Telescope mission. In this experiment, we generated a simulated 720 × 360 × 128 data set using the Galprop code that has a model of the diffuse gammaray emission of the Milky Way. 128 energy planes are logarithmically spaced from 30 MeV to 50 GeV. A six months LAT data set was created by multiplying the simulated data set with the exposure (6 months), and by convolving each energy band with the point spread function of the LAT instrument. Finally we have created the noisy observations assuming a Poisson noise distribution. Fig. 6 shows from top to bottom the original simulated data, the noisy data and the denoised data for energy band 171-181 Mev.

CONCLUSION

In this paper, we propose a multivalued extension of MSVST associated with an appropriate 2D-1D wavelet transform, which proved very efficient to denoise Poisson count data. The proposed algorithm is performs as well as the 2D-MSVST applied to summed frames. But unlike 2D denoising, 2D-1D extension fully exploits the information in the whole data set and recovers the information along the z-axis which is of paramount importance in many science fields such as hyperspectral imaging.

Figure 1 .

 1 Figure 1. Behavior of the expectation E[Z] (left) and variance Var [Z] (right) as a function of the underlying intensity, for the Anscombe VST, 2D Haar-Fisz VST, and out VST with the 2D B3-Spline filter as a low-pass filter h.

Fig. 1

 1 Fig.1shows the Monte-Carlo estimates of the expectation E[Z] (left) and the variance Var [Z] (right) obtained from 2 • 10 5 Poisson noise realizations of X, plotted as a function of the intensity λ for both Anscombe 2 (dasheddotted), Haar-Fisz (dashed)[START_REF] Fryźlewicz | A Haar-Fisz algorithm for Poisson intensity estimation[END_REF] and our VST with the 2D B 3 -Spline filter as a low-pass filter h (solid). The asymptotic bounds (dots) (i.e. 1 for the variance and √ λ for the expectation) are also shown. It can be seen that for increasing intensity, E[Z] and Var [Z] approach the theoretical bounds at different rates depending on the VST used. Quantitatively, Poisson variables transformed using the Anscombe VST can be reasonably considered to be unbiased and stabilized for λ 10, using Haar-Fisz for λ 1, and using out VST (after low-pass filtering with the chosen h) for λ 0.1.

Figure 2 .

 2 Figure 2. Image obtained by summing along the z-axis of the noisy data set.

Fig. 4 (Figure 3 .

 43 Fig. 4(a)-(b) display respectively a frame of the input noisy data set and the same frame after the 2D-1D-MSVST denoising. Fig. 5 displays the obtained spectra at the centers of two distinct sources. This result is rather impressive given that most of the counts in the 2D-1D data set are at most 1-2 photons (shown in white in 4(a)).

Figure 4 .

 4 Figure 4. (a): a frame from the noisy data set and (b): the same frame after the 2D-1D-MSVST denoising at threshold τ = 6.

Figure 5 .

 5 Figure 5. Pixel spectra at two different spatial locations after the 2D-1D-MSVST filtering.

Figure 6 .

 6 Figure 6. Top to bottom: simulated data of the diffuse gamma-ray emission of the Milky Way in energy band 171-181 Mev, noisy simulated data and denoised data using the 2D-1D-MSVST.

Further author information: (Send correspondence to M.J. Fadili) E-mail: Jalal.Fadili@greyc.ensicaen.fr, Telephone:


M Anwar Hasan

Christophe Negre

Low Space Complexity Multiplication over Binary Fields with Dickson Polynomial Representation

Keywords: Binary field, Dickson basis, Toeplitz matrix, multiplier, parallel, sequential

We study Dickson bases for binary field representation. Such a representation seems interesting when no optimal normal basis exists for the field. We express the product of two field elements as Toeplitz or Hankel matrix-vector products. This provides a parallel multiplier which is subquadratic in space and logarithmic in time. Using the matrix-vector formulation of the field multiplication, we also present sequential multiplier structures with linear space complexity.

INTRODUCTION

Finite field arithmetic is extensively used in cryptography. For public key cryptosystems, the size (i.e. the number of elements) of the field may be quite large, say 2 2048 . Finite field multiplication over such a large field requires a considerable amount of resources (time or space). For binary extension fields, used in many practical public key cryptosystems, field elements can be represented with respect to a normal basis, where squaring operations are almost free of cost. In order to reduce the cost of multiplication over the extension field, instead of using an arbitrary normal basis, it is desirable to use an optimal normal basis. The latter however does not exist for all extension fields, in which case one may use Dickson bases [START_REF] Dickson | The analytic representation of substitutions on a power of a prime number of letters with a discussion of the linear group[END_REF], [START_REF] Mullin | Dickson Bases and Finite Fields[END_REF] and develop an efficient field multiplier.

In this paper first we consider subquadratic space complexity bit parallel multipliers using the Dickson basis. To this end, using low weight Dickson polynomials, we formulate the problem of field multiplication as a product of a Toeplitz or Hankel matrix and a vector, and apply subquadratic space complexity algorithm for the product [START_REF] Fan | A New Approach to Sub-quadratic Space Complexity Parallel Multipliers for Extended Binary Fields[END_REF], which gives us a subquadratic space complexity field multiplier. Using the matrix-vector product formulation, we then develop sequential multipliers. For such multipliers, we consider both bit-serial and bit-parallel output formats.

The article is organized as follows. In Section 2 we present some general results on Dickson polynomials. Then in Section 3 we give a matrix-vector product approach for field multiplication using the Dickson basis representation. We use low weight Dickson polynomials and present parallel multipliers of subquadratic space complexity. In Section 4, we develop sequential multipliers that have linear space complexity. We wind up this article with a brief conclusion in Section 5.

DICKSON POLYNOMIALS

Dickson polynomials over finite fields were introduced by L.E. Dickson in [START_REF] Dickson | The analytic representation of substitutions on a power of a prime number of letters with a discussion of the linear group[END_REF]. These polynomials have several applications and interesting properties, the main one being a permutation property over finite fields. For a complete explanation on this the reader may refer to [START_REF] Lild | Dickson Polynomials[END_REF]. Our interest here concerns the use of Dickson polynomial for finite field representation for efficient binary field multiplication. There are two kinds of Dickson polynomials, and there are several ways to define and construct both of them. We give here the definition of [START_REF] Lild | Dickson Polynomials[END_REF] of the first kind Dickson polynomials. [START_REF] Lild | Dickson Polynomials[END_REF] page 9] Let R be a ring and a ∈ R. The Dickson polynomial of the first kind D n (X, a) is defined by

Definition 1: [Dickson Polynomial

D n (X, a) = n/2 i=0 n n -i   n -i n   (-a) i X n-2i . ( 1
)
For n = 0, we set D 0 (X, a) = 2 and for n = 1 we have D 1 (X, a) = X.

In this paper, we will consider only

β i = D i (X, 1) the Dickson polynomials in F 2 [X].
Theorem 1: Let P be an irreducible polynomial of degree n in

F 2 [X]. The system B = {β 1 , . . . , β n } forms a basis of F 2 n = F 2 [X]/(P ) over F 2 .
Proof: For a detailed proof we refer to [START_REF] Hasan | Subquadratic Space Complexity Multiplication over Binary Fields with Dickson Polynomial Representation[END_REF], we just give a brief explanation here. Using (1), we can see that a β i = X i + terms of lower degree . This implies that the conversion matrix from {β 1 , . . . , β n } to {X, . . . , X n }, is lower triangular with 1 on the diagonal. The conversion matrix is thus invertible and since {X, . . . , X n } form a basis of

F 2 n then B = {β 1 , . . . , β n } is a basis of F 2 n = F 2 [X]/(P ).
The following theorem will be extensively used for the construction of subquadratic multipliers in the Dickson basis.

Theorem 2: We denote

β i = D i (X, 1) the n-th Dickson polynomial in F 2 [X]
. Then for all i, j ≥ 0 the following equation holds

β i β j = β i+j + β |i-j| . (2)
April 16, 2013 DRAFT

Proof: We will show it by induction on i and j. We can easily check that equation (2) holds for i, j ≤ 1.

We suppose that the equation is true for all i, j ≤ n and we prove that the equation is true for i, j ≤ n + 1.

We first prove it for i = n + 1 and j ≤ n. We have

β n+1 β j = (Xβ n + β n-1 )β j = X(β n+j + β |n-j| ) + (β n-1+j + β |n-1-j| ),
by induction hypothesis. Now we have

β n+1 β j = (Xβ n+j + β n+j-1 ) + (Xβ |n-j| + β |n-1-j| ) = β n+1+j + β |n+1-j| .
For the other case i = n + 1 and j = n + 1, the product β n+1 β n+1 is obtained using similar tricks.

FIELD MULTIPLICATION USING LOW WEIGHT DICKSON POLYNOMIALS

In this section we consider multiplication of two elements of the binary field In Table 1

F 2 n = F 2 [X]/(P )
P = β n + β k + 1 with k ≤ n/2.
For the purpose of comparison, we mention also whether for each degree an ONB of type I or II exists (marked as NI and NII).

Our main goal here is to express the product of two elements, represented in the Dickson basis, as a Toeplitz or Hankel matrix-vector products. Recall that an n × n Toeplitz matrix T = [t i,j ] satisfies t i,j = t i+1,j+1 and a Hankel matrix H = [h i,j ] satisfies h i,j = h i+1,j-1 . We will then use the subquadratic Toeplitz matrix-vector product of [START_REF] Fan | A New Approach to Sub-quadratic Space Complexity Parallel Multipliers for Extended Binary Fields[END_REF] to design a subsquadratic multiplier.

Irreducible Dickson binomials

In this subsection we focus on finite fields

F 2 n = F 2 [X]
AB = ( n i=1 a i β i ) × n i=j b j β j =   n i,j=1 a i b j β i+j   S1 +   n i,j=1 a i b j β |i-j|   S2 (3)
Now we express each sum S 1 and S 2 as matrix-vector products. Let us begin with S 1 . We remark that S 1 has a similar expression as product of two polynomials of the same degree. In other words, S 1 can April 16, 2013 DRAFT be computed as Z A • B where

Z A =                  0 0 • • • 0 0 a1 0 • • • 0 0 . . . . . . an-1 • • • • • • a1 0 an • • • • • • a2 a1 0 an • • • a3 a2 . . . . . . 0 0 • • • 0 an                  ← β1 ← β2 . . . ← βn ← βn+1 ← βn+2 . . . ← β2n (4)
We reduce the matrix Z A modulo P = β n + 1 to get non-zero coefficients only on rows corresponding to β 1 , . . . , β n . We use the fact that β n+i for i ≥ 0 satisfies

β n+i = β i β n + β n-i = β i + β n-i .
This equation is a simple consequence of ( 2) and that β n = 1 mod P . This implies that the rows corresponding to β n+i are reduced into two rows one corresponding to β i and the other to β n-i . After performing this reduction and removing zero rows we get

S 1 = Z A • B =         a n a n-1 • • • a 2 a 1 a 1 a n • • • a 3 a 2 . . . . . . a n-1 • • • • • • a 1 a n              b 1 . . . b n      S1,1 +            0 0 • • • a n a n-1 . . . . . . 0 a n • • • a 3 a 2 a n a n-1 • • • a 2 a 1 0 • • • • • • 0 0                 b 1 . . . b n      S1,2
Finally, we get an expression of S 1 as matrix-vector product where the matrix is a sum of a Toeplitz and an essentially Hankel matrix. Now we do the same for S 2 . We split S 2 into two sums

S 2 =   n k=1 n-k j=1 a j+k b j β k   S2,1 +   n k=1 n j=k a j-k b j β k   S2,2 . (5)
We express S 2,1 and S 2,2 as matrix-vector products

S 2,1 =            a 2 a 3 • • • a n-1 a n 0 a 3 a 4 • • • a n 0 0 . . . . . . a n 0 • • • • • • 0 0 0 • • • • • • • • • 0 0            •      b 1 . . . b n      , (6)
S 2,2 =            0 a 1 a 2 • • • a n-1 0 0 a 1 • • • a n-2 . . . . . . 0 0 • • • • • • a 1 0 • • • • • • • • • 0            •      b 1 . . . b n      . ( 7
)
So now we have each of S 1 and S 2 in the required form. We finally write S 1,1 + S 2,2 = T A • B where T A is a Toeplitz matrix and S 1,2 + S 2,1 = H A • B where H A is an Hankel matrix. We obtain

A × B = T A • B + H A • B (8)
as stated at the beginning of the current subsection.

Dickson trinomials

Now we assume that the field F 2 n is defined by a three-term irreducible Dickson trinomial

P P = 1 + β k + β n , with k ≤ n/2.
The elements in F 2 n = F 2 [X]/(P ) are expressed in the Dickson basis B = {β 1 , . . . , β n }. Our aim is to express the product of two elements A and B of F 2 n as Toeplitz or Hankel matrix-vector product. We use again the expression of the product C = A × B = S 1 + S 2 given in equation [START_REF] Fan | A New Approach to Sub-quadratic Space Complexity Parallel Multipliers for Extended Binary Fields[END_REF]. Similar to the previous subsection, here we express S 1 and S 2 as matrix-vector product separately. Specifically

1) The sum S 1 is expressed as

Z A • B where Z A is given in (4)
2) For S 2 we use the expression of ( 5) and we put this expression in a matrix-vector product form.

S2 =                              a2 a3 • • • an 0 a3 a4 • • • 0 0 . . . . . . . . . an 0 • • • 0 0 0 0 • • • 0 0            +                   0 a1 • • • an-1 0 0 • • • an-2 . . . . . . 0 0 • • • a1 0 • • • • • • 0 . . . . . . 0 • • • • • • 0                                     • B. (9)
April 16, 2013 DRAFT Now we replace S 1 and S 2 by their corresponding expressions given above in A × B = S 1 + S 2 . We get

                                          0 0 • • • 0 a1 0 • • • 0 . . . . . . an-1 • • • • • • 0 an • • • • • • a1 0 an • • • a2 . . . . . . . . . 0 0 • • • an                      +                      0 a1 a2 • • • an-1 0 0 a1 • • • an-2 . . . . . . . . . 0 • • • • • • • • • a1 0 • • • • • • • • • 0 0 • • • • • • • • • 0 . . . . . . 0 • • • • • • • • • 0                                           • B +                   a2 a3 • • • an-1 an 0 a3 a4 • • • an 0 0 . . . . . . . . . an 0 • • • • • • 0 0 0 • • • • • • • • • 0 0 . . . . . . . . . 0 • • • • • • • • • 0 0                   •       b1 . . . bn       = A × B . (10)
In [START_REF] Reyhani-Masoleh | Low Complexity Sequential Normal Basis Multipliers over GF(2 m )[END_REF] the addition of two 2n × n Toeplitz matrices results in one single 2n × n Toeplitz matrix. The latter can be horizontally split in the middle to obtain two n × n Toeplitz matrices, say T up and T down , which can be then multiplied separately with vector (b 1 , . . . , b n ) with a total cost of two n × n Toeplitz matrix-vector products.

The other 2n × n Hankel matrix in [START_REF] Reyhani-Masoleh | Low Complexity Sequential Normal Basis Multipliers over GF(2 m )[END_REF] has all zero in the lower n rows, contributing nothing to the cost of the matrix-vector multiplication. Thus, the total computational cost of ( 10) is no more than three n × n Toeplitz or Hankel matrix-vector products.

The reduction:

The resulting expression of C in ( 10) is an unreduced form of A×B, since it has non zero coefficients c i on rows i = n + 1, . . . , 2n. These coefficients are obtained by multiplying T down with vector

(b 1 , b 2 , . . . , b n ),
and must be reduced modulo P = β n + β k + 1, to get an expression of C in B. We have

β i = β n β i-n + β 2n-i = (β k + 1)β i-n + β 2n-i = β i-n+k + β |i-n-k| + β i-n + β 2n-i .
We reduce the expressing of C = 2n i=1 c i β i by replacing each β i for i > n by the expression given above. Since we assume k < n 2 this process must be done two times to get a reduced expression of C. A circuit can be designed to perform this process which requires 6n -k XOR gates and is performed in time 3T X (see [START_REF] Hasan | Subquadratic Space Complexity Multiplication over Binary Fields with Dickson Polynomial Representation[END_REF] for details).

Parallel multiplier

We can design multiplier using the expression of the multiplication in F 2 n as a Toeplitz or Hankel matrixvector product (TMVP). Specifically we use the Toeplitz or Hankel matrix-vector multiplier presented April 16, 2013 DRAFT in [START_REF] Fan | A New Approach to Sub-quadratic Space Complexity Parallel Multipliers for Extended Binary Fields[END_REF] to perform these products. In Table 2, we recall the complexity of the TMVP multiplier established by Fan and Hasan [START_REF] Fan | A New Approach to Sub-quadratic Space Complexity Parallel Multipliers for Extended Binary Fields[END_REF].

Table 2 Asymptotic complexity of TMVP multiplier 2-way split method 3-way split method # AND n log 2 (3) n log 3 (6) # XOR 5.5n log 2 (3) -6n + 0.5 24 5 n log 3 (6) -5n + 1 5

Delay T A + 2 log 2 (n)T X T A + 3 log 3 (n)T X
In the case of Dickson binomials, to compute the matrix-vector products of ( 8) we need two TMVP multipliers in parallel. Each of them can use 2-way or 3-way split approach of [START_REF] Fan | A New Approach to Sub-quadratic Space Complexity Parallel Multipliers for Extended Binary Fields[END_REF]. We also need additionnal 2n XOR gates to compute the coefficient of T A and add the result of the two matrix-vector products.

In the case Dickson trinomials, as specified in Subsection 3.2, three TMVPs are done in parallel using 2-way or 3-way split approach of [START_REF] Fan | A New Approach to Sub-quadratic Space Complexity Parallel Multipliers for Extended Binary Fields[END_REF]. We also need to perform a reduction using the circuit depicted in [START_REF] Hasan | Subquadratic Space Complexity Multiplication over Binary Fields with Dickson Polynomial Representation[END_REF]. We obtain the complexities of Table 3 below where the second left most column indicates b-way splits with the value of b being either 2 or 3. (6) 48/5n log 3 (6) (3 log

11n log 2 (3) -11n (2 log 2 (n) + 1)T X +T A 3 2n log 3
3 (n) + 1)T X -11n + 3/5 +T A DT 2 7 3 n log 2 (3) 38,5 3 n log 2 (3) (2 log 2 (n) + 6)T X -6, 5n -k + 2, 5 +T A 3 2n log 3 (6) 48/5n log 3 (6) (3 log 3 (n) + 5)T X -3n -k + 7/5 +T A ONBI 2 n log 2 (3) + n 5.5n log 2 (3) (2 log 2 (n) + 1)T X -4n -0.5 +T A [4] 3 n log 3 (6) + n 24/5n log 3 (6) (3 log 3 (n) + 1)T X -3n -4/5 +T A ONBII 2 n log 2 (3) 6n log 2 (3) -3n+ (3 log 2 (n) + 1)T X [5], [9] † (8n + 2) log 2 (2n + 1) +T A ONBII 3 n log 3 (6) 72 15 n log 3 (6) -7 3 n -1+ (4 log 3 (n) + 1)T X [5], [11] ‡ (8n + 2) log 2 (2n + 1) +T A
The row of Table 3 labelled by † (resp. ‡) refers to the method of [START_REF] Zur Gathen | Efficient multiplication using type 2 optimal normal bases[END_REF] combined to the polynomial multiplication of [START_REF] Paar | A new architecture for a parallel finite field multiplier with low complexity based on composite fields[END_REF] (resp. [START_REF] Sunar | A Generalized Method for Constructing Subquadratic Complexity GF(2 k ) Multipliers[END_REF]).

April 16, 2013 DRAFT

In a recent paper Mullin et al. [START_REF] Mullin | Dickson Bases and Finite Fields[END_REF] pointed out that there were some links between the Dickson basis and the normal basis. In practice, a Dickson basis is interesting when no optimal normal basis exists for the considered field. This is the case for NIST recommended binary fields F 2 163 and F 2 283 . Using Table 1, we can remark that NIST fields can be constructed with Dickson trinomials, and thus we obtain a subquadratic multiplier in each of these cases.

SEQUENTIAL MULTIPLIERS

In this section, we present sequential multipliers. Each of these multipliers takes O(n) clock cycles but has a space complexity of O(n).

Using irreducible Dickson binomials

Multiplier with bit serial output

In the sequel, we denote the entry at (i, j) of the Toeplitz and the Hankel matrices of (8) as T i,j and H i,j , respectively. We also denote the rows of the Toeplitz matrix as T 1, , T 2, , . . . , T n, and those of the essentially Hankel matrix as H 1, , H 2, , . . . , H n, . Thus we can write

A × B =       c1 c2 . . . cn       =       T1, + H1, T2, + H2, . . . Tn, + Hn,       • B. (11)
We remark that a) T n, consists of the coordinates of input A that are rotated left one position, i.e., T n, = a n-1 a n-2 • • • a 1 a n .

On the other hand, H n, is the all zero row vector and

H n-1, = 0 a n-1 a n-2 • • • a 2 a 1 .
b) Given T i, and H i,n-1 , we can express

T i-1, for 1 < i ≤ n as T i-1, = T i,2 T i,3 • • • T i,n T i,1 + H i,n-1 .
Furthermore, given the row H i, and the entry T i+1,1 we can express

H i-1, for 1 < i ≤ n -1 as follows H i-1, = T i+1,1 H i,1 • • • H i,n-2 H i,n-1 .
The following diagram (Figure 1) corresponds to a sequential structure to realize the multiplication [START_REF] Sunar | A Generalized Method for Constructing Subquadratic Complexity GF(2 k ) Multipliers[END_REF]. In the initial clock cycle, the left side register (LR) in the diagram is loaded with T n, and the right side register (RR) with H n, . In this cycle, rows T n, and H n, are added and an inner product is performed to yield c n = (T n, + H n, ) • B. Also, in this cycle the output of MUX is a 1 (and in other cycles the MUX output is the second right most bit of RR). In the next cycle, RR is loaded with H n-1, and LR is shifted left to generate T n-1, eventually yielding c n-1 .

C = A × B in accordance with
For each of the following n -2 clock cycles, LR is shifted left, RR is shifted right, their contents are bit-wise added and an inner product is performed to produce one coordinate of C. The space and time complexity of the architecture of Fig. 1 is given in Table 4.

April 16, 2013 DRAFT Referring to [START_REF] Mullin | Dickson Bases and Finite Fields[END_REF] we denote the columns of a Toeplitz matrix as T ,i for 1 ≤ i ≤ n and those of an essentially Hankel matrix as

H ,i for 1 ≤ i ≤ n. Thus we can write A×B = T ,1 + H ,1 T ,2 + H ,2 • • • T ,n + H ,n • B, i.e., C = n i=1 b i (T ,i + H ,i ). (12)
We remark that

a) T ,1 = a n a 1 • • • a n-1 t and H ,1 = a 2 a 3 • • • a n-1 0 0 t b)
Given the column T ,i and the entry H 1,i-1 , we can express T ,i+1 as

T ,i+1 = T n,i + H 1,i-1 , T 1,i , • • • , T n-1,i t , 1 ≤ i ≤ n
, where H 1,0 is assumed to be a 1 .

Additionally, given column H ,i and entry T n,i , we can express H ,i+1 as

H ,i+1 = H 2,i , H 3,i , • • • , H n-1,i , T n,i , 0 t , 1 i ≤ n.
In the following diagram (Fig. 2), the column vectors T ,1 and H ,1 are initially loaded into the top register (TR) and the bottom register (BR) respectively. The one-bit feedback cell F is initialized with

H 1,0 = a 1 .
If TR is shifted downward and BR upward with the feedback connections as shown in the diagram, the new contents of TR and BR will be T ,2 and H ,2 respectively. Note that BR is an n -1 bits long shift register, since H n,i = 0 for 1 ≤ i ≤ n. With additional shifts on TR and BR, the remaining columns of the Toeplitz and the essentially Hankel matrices are generated.

Each corresponding pair of columns (i.e., T ,i and H ,i ) are added and the resulting columns are multiplied with b i (in the diagram these are shown using an array of XOR and AND gates).

The weighted columns are accumulated in accordance with [START_REF] Yang | Modified Sequential Normal Basis Multipliers for Type II Optimal Normal Bases[END_REF] to produce the desired output C in a total of n clock cycles. The delay and the space complexity of this architecture are given in Table 4.

April 16, 2013 DRAFT

                      0 a1 • • • an-1 a1 0 • • • an-2 . . . . . . . . . an-2 an-3 • • • a1 an-1 an-2 • • • 0            +            a2 a3 • • • an 0 a3 a4 • • • 0 0 . . . . . . . . . an 0 • • • 0 0 0 0 • • • 0 0                       • B, (13) and 
       c n+1 c n+2 . . . c 2n         =         a n a n-1 • • • a 1 0 a n • • • a 2 . . . 0 0 • • • a n         •         b 1 b 2 . . . b n         . ( 14
)
Note that c n+1 , c n+2 , • • • , c 2n can be reduced as explained at the end of Subsection 3.2.

Below we will first present a hardware structure to generate c 1 , c 2 , . . . , c n in accordance with (13). Then

Sequential multiplier with bit serial output

We denote the rows of the Toeplitz and the Hankel matrices of (13) as T i, and H i, respectively. For 1 ≤ i ≤ n, we can then write the following

T i+1, = H i-1,1 T i,1 • • • T i,n-2 T i,n-1 , H i+1, = H i,2 H i,3 • • • H i,n 0 , (15)
April 16, 2013 DRAFT where H 0,1 is assumed to be a 1 and

T 1, = 0 a 1 a 2 • • • a n-2 a n-1 , H 1, = a 2 a 3 a 4 • • • a n 0 . (16)
In Fig. 3 below, registers RR and LR are initialized with T 1, and H 1, . The feedback cell F is initialized with a 1 = H 0,1 . Then with the application of a shift to these registers together, the second rows of the Toeplitz and the Hankel matrices of ( 12) are formed in RR and LR respectively. This happens due to the fact that the shift and the feedback connection as shown in Fig. 3 essentially realize (15). The remaining rows of the two matrices are formed pair by pair with successive shifts.

Note that LR is n -1 bits long, since the right most bit of each row of the Hankel matrix is zero. The upper part of Fig. 3 is similar to that of Fig. 1 and is to add the corresponding rows of the Toeplitz and the Hankel matrices, followed by inner product operations to yield c 1 , c 2 , • • • , c n .

To generate c n+1 , c n+2 , . . . , c 2n using the structure in Fig. 3, we initialize RR with [a n , a n-1 , . . . , a 1 ],

which is the first row of the upper triangular Toeplitz matrix of (14). Register LR and cell F are initialized with all zeros. Then with successive shifts, RR will contain the remaining rows of the Toeplitz matrix and LR will have all zeros. This will result in c n+1 , c n+2 , . . . , c 2 at the output of Fig. 3.

The time and the space complexities of the structure of Fig. 3 are given in Table 4. These exclude the cost associated with the reduction of c n+1 , c n+2 , • • • , c 2n .

Bit parallel output

For (13), let T ,i and H ,i be the i-th columns of the Toeplitz and the Hankel matrices, respectively. Then (13) can be re-written as

c 1 c 2 • • • c n t = n i=1 b i (T ,i + H ,i ).
April 16, 2013 DRAFT

In order to generate T ,i and H ,i we note that in (13) T ,i = T i, and H ,i = H i, . In other words, the i-th column is the same as the i-th row for each of the matrices. Thus columns can be generated using the same system of feedback registers as shown in Fig. 3 earlier.

To obtain c 1 , c 2 , • • • , c n in bit parallel fashion the inner product unit of Fig. 3 can be replaced with an array unit of weighting AND gates and accumulators (as used in Fig. 2). The complete diagram is shown in Fig. 4, and its space and time complexities are given in Table 4.

Complexity and comparison

In Table 4 we put the resulting complexities of the different sequential multipliers based on the Dickson basis representation. For the purpose of comparison, we also give the complexity of the method of [START_REF] Fan | Subquadratic Computational Complexity Schemes for Extended Binary Field Multiplication Using Optimal Normal Bases[END_REF] using ONB of type I and II. We remark that when no ONB is available, a Dickson binomial seems to be the best choice since Dickson trinomial require an increased number of clock cycles.

CONCLUSION

In this paper we have presented new parallel multipliers based on Dickson basis representation of binary fields. The multiplier for an irreducible Dickson binomial has a complexity similar to the subquadratic multiplier for ONB II of [START_REF] Fan | Subquadratic Computational Complexity Schemes for Extended Binary Field Multiplication Using Optimal Normal Bases[END_REF]. For an irreducible Dickson trinomial, the multiplier has a slightly more space complexity, but can still be used for fields with degree of several hundreds (for example those used in today's elliptic curve cryptographic systems).

April 16, 2013 DRAFT In this paper, we have also presented sequential multipliers using the above mentioned Dickson representation. The sequential multipliers have a space complexity of O(n). We have considered both bit-serial and bit-parallel output formats for the sequential multipliers. Compared to the sequential multipliers with bit-parallel output format presented in [START_REF] Ansari | Revisiting Finite Field Multiplication Using Dickson Bases[END_REF] and [START_REF] Mullin | Dickson Bases and Finite Fields[END_REF], the sequential multipliers presented here with the same output format reduce the number of XOR and AND gates by a factor of two or more, while keeping the number of flip-flops and clock cycles about the same.

  where the polynomial P is a low weight Dickson polynomial. In particular we consider two and three-term Dickson polynomials P , i.e., Dickson binomials and trinomials. Like low weight conventional polynomials the use of low weight Dickson polynomials is expected to yield lower space complexity multipliers.

Figure 1 . 2

 12 Figure 1. Sequential multiplier with bit serial out using Dickson binomials

Figure 2 .From ( 10 )

 210 Figure 2. Sequential multiplier with bit parallel output using Dickson binomials

  we will discuss how to use part of the above hardware to generate c n+1 , c n+2 , . . . , c 2n . In practice, one can first generate c n+1 , c n+2 , . . . , c 2n . While these n bits are reduced, one can generate c 1 , c 2 , . . . , c n . This overlap of operations will effectively eliminate/hide the extra time for reduction of c n+1 , . . . , c 2n .

Figure 3 . 1 , c 2 ,

 312 Figure 3. Bit serial output sequential multiplier for Dickson trinomials

Figure 4 . 1

 41 Figure 4. Bit parallel output sequential multiplier for Dickson trinomials

  we give the degree n ∈ [160, 285] of field F 2 n = F 2 [X]/(P ) where P is a low weight Dickson polynomial. Specifically, since no irreducible Dickson binomials were available, we have looked for Almost

	Dickson binomials (ADB) with irreducibible P satisfying P × (X + 1) = β n+1 + 1. We also give Dickson
	trinomials (DT) of the form

  /(P ) where P is a two-term polynomial of the form P = β n + 1 where a β n is the n-th Dickson polynomial.The elements of F 2 n are expressed in the Dickson basis B = {β 1 , . . . , β n }. Now, our main goal is to show that the product of two elements A and B in F 2 n can be computed as a matrix-vector product M A • B where M A is a sum of a Toeplitz matrix and an essentially Hankel matrix.

				Table 1
		Irreducible Dickson binomials and trinomials
	n		n		n		n
	163	DT	195		226	DT, NI	257
	164		196	NI	227	DT	258
	165	DT	197	DT	228		259	DT
	166		198	ADB	229	DT	260	DT
	167	ADB,DT	199	DT	230	NII	261	DT,NII
	168		200	DT	231	DT, NII	262	ADB
	169		201	DT	232		263	DT
	170	DT	202	DT	233	NII	264
	171	DT	203	DT	234		265	DT
	172	DT, NI	204		235	DT	266
	173 ADB,DT,NII	205	DT	236		267
	174	NII	206		237	DT	268 ADB,DT,NI
	175	DT	207	DT	238	ADB	269	DT
	176	DT	208	DT	239	DT, NII	270	ADB,NII
	177		209	NII	240		271	DT
	178	DT,NI	210	NI,NII	241	DT	272	DT
	179	DT,NII	211	DT	242	DT	273	DT, NII
	180	NI	212	DT	243	NII	274	DT
	181	DT	213		244	DT	275	DT
	183	DT,NII	214		245	DT,NII	276
	184	DT	215	DT	246		277	DT
	185		216		247	DT	278	NII
	186	NII	217		248	DT	279	DT
	187	DT	218	DT	249		280	DT
	188	DT	219	DT	250	DT	281	NII
	189	DT, NII	220	DT	251	DT,NII	282
	190		221	DT,NII	252	ADB	283	DT
	191	DT,NII	222		253	DT	284
	192		223	DT	254	NII	285	DT
	193	DT	224	DT	255	DT
	194	DT,NII	225	DT	256	DT
	April 16, 2013							DRAFT

ADB=Dickson Binomial,DT=Dickson Trinomial,NI=ONBI and NII=ONBII If we multiply two elements A and B expressed in B and if we use Theorem 2 we get the following

Table 3

 3 Comparison of Subquadratic Space Complexity Parallel Multipliers

		b	Space	Time
			# AND	# XOR
	DB	2	2n log 2 (3)

Table 4

 4 Complexity of sequential multipliersDB=Dickson Binomial, DT=Dickson Trin., DG=General Dickson , CC=Clock Cycle.

	Archi.	#AND	#XOR	#FF	#MUX #CC	Delay
	DB	n	2n	2n + 1	1	n	T A +
	Fig 1						(1 + log 2 (n) )T X
	DB Fig 2	n	2n	2n + 1	1	n	T A + 2T X
	DT	n	2n -2	2n	0	2n	T A +
	Fig 3						(1 + log 2 (n) )T X
	DT Fig 4	n	2n -1	3n -1	0	2n	T A + T X
	DG[1]	2n	4n -3	3n	0	n	2T X + T A
	ONBI[10]	n	3n 2	2n	0	n	T A + 2T X
	ONBII[12]	n	3n-1 2	3n	0	n	T A + 2T X

ACKNOWLEDGMENT

A preliminary version of this work was presented at the WAIFI 2008 conference [6]. This work was supported in part by an NSERC research grant awarded to Dr. Hasan.


J Bobin

Y Moudden

J.-L Starck

J Fadili

SPARSITY AND MORPHOLOGICAL DIVERSITY FOR HYPERSPECTRAL DATA ANALYSIS

Keywords: source separation, hyperspectral data, sparsity, remote sensing, wavelets, astronomy

Recently morphological diversity and sparsity have emerged as new and effective sources of diversity for Blind Source Separation. Based on these new concepts, novel methods such as Generalized Morphological Component Analysis have been put forward. The latter takes advantage of the very sparse representation of structured data in large overcomplete dictionaries, to separate sources based on their morphology. Building on GMCA, the purpose of this contribution is to describe a new algorithm for hyperspectral data processing. Large-scale hyperspectral data refers to collected data that exhibit sparse spectral signatures in addition to sparse spatial morphologies, in specified dictionaries of spectral and spatial waveforms. Numerical experiments are reported which demonstrate the validity of the proposed extension for solving source separation problems involving hyperspectral data.

INTRODUCTION

Generalized Morphological Component Analysis (GMCA) is a recent algorithm for multichannel data analysis and blind source separation introduced in [START_REF] Bobin | Sparsity and morphological diversity in blind source separation[END_REF], which relies on sparsity and morphological diversity to disentangle the source signals from observed mixtures. Indeed, sparsity has been recognized recently as an effective source of diversity for BSS [START_REF] Zibulevsky | Blind source separation by sparse decomposition in a signal dictionary[END_REF]. The successes of GMCA in many generic multichannel data processing applications, including BSS, color image restoration and inpainting [START_REF] Bobin | Blind Source Separation : the Sparsity Revolution[END_REF][START_REF] Bobin | Sparsity and morphological diversity in blind source separation[END_REF], strongly motivated research to extend its applicability. In particular, there are instances where one is urged by prior knowledge to set additional constraints on the estimated parameters (e.g. equality constraints, positivity). Building on GMCA, the purpose of this contribution is to describe an algorithm for so-called hyperspectral data processing. The modified GMCA algorithm this entails is then given in section 3. Finally, numerical experiments in section 4 demonstrate the efficiency of the proposed method to deal with Hyperspectral data in component separation problems.

MODEL AND OBJECTIVE FUNCTION

Hyperspectral imaging systems collect data in a large number (up to several hundreds) of contiguous spectral intervals so that it makes sense to consider e.g. the regularity of some measurement from one channel to its neighbors. Typically, the spectral signatures of the structures in the data may be known a priori to have a sparse representation in some specified possibly redundant dictionary of template spectral waveforms. In what follows, regardless of other definitions or models living in other scientific communities, the term hyperspectral will be used generically to identify multichannel data with the above two specific properties i.e. that the number of channels is large and that these achieve a regular if not uniform sampling of some additional and meaningful physical index (e.g. wavelength, space, time) which we refer to as the spectral dimension. We further assume that hyperspectral data is structured a priori according to the simple yet common instantaneous linear mixture model given as follows :

X = AS + N = n k=1 X k + N (1)
where the measurements X ∈ R m,t are a mixture of contributions from various objects X k with different statistical and spatio-spectral properties. These induce inter-and intra-sensor structures or signal coherence in space and across sensors which we classically assume to be well represented by rank one matrices, product of a spectral signature a k ∈ R m,1 and a spatial density profile s k ∈ R 1,t . Here, s k and a k are respectively the k th row and column of S ∈ R n,t and A ∈ R m,n . The m × t random matrix N is included to account for Gaussian instrumental noise assumed, for simplicity, to be uncorrelated inter-and intra-channel with variance σ 2 . Initially, the GMCA algorithm aims at solving the following non-convex problem :

min A,S k λ k ν k 1 + 1 2σ 2 X - k a k s k 2 (2)
where s k = ν k Φ. We describe next the proposed modified GMCA algorithm for hyperspectral data processing when it is known a priori that the underlying objects of interest X k = a k s k exhibit sparse spectral signatures and sparse spatial morphologies in known dictionaries of spectral and spatial waveforms. Accounting for this prior is made by assuming that each component X k = a k s k has a sparse representation in a given multichannel signal waveform dictionnary {Φ, Ψ}. The derived optimization problem for hyperspectral data is then the following :

min A,S k λ k γ k ν k 1 + 1 2σ 2 X - n k=1 a k s k 2 (3)
where

X k = a k s k = Φγ k ν k Ψ.

GMCA ALGORITHM FOR HYPERSPECTRAL DATA

For the sake of simplicity, consider now that the Ψ and Φ reduce to respectively spectral and spatial waveforms orthogonal bases. In this case, the minimization problem ( 3) is best formulated in coefficient space as follows :

min {γ k ,ν k } 1 2σ 2 α -γν 2 + n k=1 λ k γ k ν k 1 (4)
where the columns of γ are γ k , the rows of ν are ν k and α = Ψ T XΦ T is the coefficient matrix of data X in Ω. Thus, we are seeking a decomposition of matrix α into a sum of sparse rank one matrices α k = γ k ν k . Unfortunately, there is no obvious closed form solutions to problem [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]. We propose instead a numerical approach by means of a block-coordinate relaxation iterative algorithm, alternately minimizing with respect to γ and ν. Indeed, thanks to the chosen prior, for fixed γ (resp. ν), the marginal minimization problem over ν (resp. γ) is convex. Inspired by [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF], we obtain the following system of update rules, akin to a Projected Landweber algorithm [START_REF] Blumensath | Iterative thresholding for sparse approximations[END_REF] and after some algebra (see [START_REF] Bobin | Blind Source Separation : the Sparsity Revolution[END_REF]) :

ν (+) = ∆ η γ T γ -1 γ T α (5)
γ (+) = ∆ ζ αν T νν T -1 (6)
where vector η has length n and entries η

[k] = λ k γ k 1 / γ k 2 2 , while ζ has length m and entries ζ[k] = λ k ν k 1 / ν k 2 2 .
The multichannel soft-thresholding operator ∆ η acts on each row k of ν with threshold η[k] and ∆ ζ acts on each column k of γ with threshold ζ[k]. Equations ( 5) and ( 6) rules are easily interpreted as thresholded alternate least squares solutions. Finally, in the spirit of the GMCA algorithm [START_REF] Bobin | Blind Source Separation : the Sparsity Revolution[END_REF][START_REF] Bobin | Sparsity and morphological diversity in blind source separation[END_REF], it is proposed that a solution to problem (4) can be approached efficiently using the following symmetric iterative thresholding scheme with a progressively decreasing threshold, which we refer to as hypGMCA :

1. Set the number of iterations I max and initial thresholds λ

(0) k 2.
Transform the data X into α 3. While λ (h) k are higher than a given lower bound λ min , -Update ν assuming γ is fixed using equation [START_REF] Starck | The curvelet transform for image denoising[END_REF].

-Update γ assuming ν is fixed using equation [START_REF] Zibulevsky | Blind source separation by sparse decomposition in a signal dictionary[END_REF].

-Decrease the thresholds λ (h) k . 5. Transform back γ and ν to estimate A and S.

The salient to fine process is also the core of hypGMCA. With the threshold successively decreasing towards zero along iterations, the current sparse approximations for γ and ν are progressively refined by including finer structures spatially and spectrally, alternatingly. The final threshold should vanish in the noiseless case or it may be set to a multiple of the noise standard deviation as in common detection or denoising methods. When nonunitary or redundant transforms are used, the above is no longer strictly valid. Nevertheless, simple shrinkage still gives satisfactory results in practice.

NUMERICAL EXPERIMENTS

In this section, we report on two toy BSS experiments in 1D and 2D to compare GMCA to its extension hypGMCA. First we consider synthetic 2D data consisting of m = 128 mixtures of n = 5 image sources. The sources are a set of 128 × 128 structured images. The spectra were generated as sparse process in some orthogonal wavelet domain given a priori. The wavelet coefficients of the spectra were sampled from a Laplacian probability density with scale parameter µ = 1. Finally, white Gaussian noise with variance σ 2 was added to the pixels of the synthetic mixture data in the different channels. Figure 1 pection of figure 2 allows a first qualitative assessment of the improved source recovery provided by correctly accounting for a priori spatial as well as spectral sparsity. The top images were obtained with GMCA while the bottom images, were obtained with hypGMCA. In all cases, both methods were run in the curvelet domain [START_REF] Starck | The curvelet transform for image denoising[END_REF] with the same number of iterations. The graph on figure 3 gives more quantitative results. It traces the evolution of the mixing matrix criterion C A = I n -P à † A 1 as a function of the SNR which was varied from 0 to 40dB, where P serves to reduce the scale and permutation indeterminacy of the mixing model and à † is the pseudoinverse of the estimated mixing matrix. In simulation, the true source and spectral matrices are known and so that P can be computed easily. The mixing matrix criterion is then strictly positive unless the mixing matrix is correctly estimated up to scale and permutation. Finally, as we expected since it benefits from the added a priori spectral sparsity constraint it enforces, the proposed hypGMCA is clearly more robust to noise.

In a second experiment, GMCA and hypGMCA are compared as the number n of sources is increased while the numbers of samples t and channels m are kept constant. Increasing the number of sources makes the separation task more difficult. We consider now 1D syn-Fig. 2. Left column : Estimated sources using the original GMCA algorithm. Right column : Estimated sources using the new hypGMCA. thetic source processes S generated from i.i.d. Laplacian probability density distributions with scale parameter µ = 1. The Dirac basis was taken as the dictionary of spatial waveforms Φ. The entries of the mixing matrix are also drawn from i.i.d. Laplacian distributions with scale parameter µ = 1 and the Dirac basis was also taken as dictionary of spectral waveforms Ψ. The data are not contaminated by noise. The number of samples is t = 2048 and the number of channels is m = 128. Figure 4 depicts the comparisons between GMCA and its extension to the hyperspectral setting. Each point of this figure has been computed as the mean over 100 trials. The panel on the left of Figure 4 features the evolution of the recovery SNR when the number of sources varies from 2 to 64. At lower n, the spatiospectral sparsity constraint only slightly enhances the source separation. However, as n becomes larger than 15 the spectral sparsity constraint clearly enhances the recovery results. For instance, when n = 64, the GMCA algorithm with the spectral sparsity constraint outperforms the original GMCA by up to 12dB. The right of Figure 4 shows the behavior of both algorithms in terms of C ℓ1 = n i=1 a i s i -ãi si 1 / n i=1 a i s i 1 . As expected, accounting for spectral sparsity yields sparser results. Furthermore, as the number of sources increases, the deviation between the aforementioned methods becomes wider.

CONCLUSION

We described a new algorithm, hypGMCA, for blind source separation in the case where it is known a priori that the spatial and spectral features in the data have sparse representations in known dictionaries of template waveforms. The proposed method relies on an iterative thresholding procedure with a progressively decreasing threshold. This alone gives the method true robustness to noise. As expected, taking into account the additional prior knowledge of spectral sparsity leads to enhanced performance. It was illustrated by numerical experiments that spatiospectral sparsity yields robustness to noise contamination, as well as statbility when the dimensionality of the problem increases. Current work is on enforcing other prior constraints such as positivity and on applications of the proposed method.

Fig. 1 .

 1 Fig. 1. Four 128 × 128 mixtures out of the 128 channels. The SNR is equal to 20dB.

Fig. 3 .

 3 Fig. 3. Evolution of the mixing matrix criterion C A as a function of the SNR in dB. Solid line : recovery results with GMCA. • : recovery results with hypGMCA.

Fig. 4 .

 4 Fig. 4. Abscissa : Number of sources. Ordinate -left : Recovery SNR. Right : sparsity-based criterion C ℓ1 . Solid line : recovery results with GMCA. • : recovery results with hypGMCA.


B Bleuzé

O Clatz

P Fillard

M Sermesant

N Toussaint

J Wintz

Horizontal Data Fusion for Integrative Modelling: the MedINRIA Fusion Toolbox

Introduction

The Virtual Physiological Human research aims at integrating in computational models an extensive representation of the human body. Given the numerous phenomena interacting in such a complex living system, several different modalities of observation are needed in order to gather information on these phenomena. Moreover following up on the time evolution implies comparing data at different points in time. Finally, superposition of simulated volumes and acquired data is also a needed tool. All these needs require a so called "Horizontal Data Fusion" which registers in the same spatial coordinates different datasets. The most active area of research on horizontal data fusion is medical image processing. Medical images come from a larger and larger variety of modalities, and the number of image types for a single modality increases just as fast. However, the problems linked to that explosion of data and the need to put them together is that each registration algorithm has its own input/output formats, graphical user interface and visualisation. Thus, its use by a non-expert is challenging and the combination or comparison of such algorithms is difficult. Therefore, there is a need for a standard in horizontal data fusion.

The aim of our project at INRIA within the Virtual Physiological Human Network of Excellence 1 is to provide a toolbox enabling the visualisation and the aggregation of data fusion tools and algorithms: Medinria-fusion 2 . Through the use of a core/plug-ins architecture, it explores solutions to make it easy for researchers to import their own tools, wrap them around a single application or ex- * N.

Demonstrator

The demonstrator, named MedINRIA-fusion is a subset of a larger platform called MedINRIA, developed at INRIA. MedINRIA-fusion builds on the architecture and GUI from MedINRIA and puts the framework through its paces by implementing registration tools within it. The programme imports and stores DICOM files. The intuitive and slick looking interface makes it easy to use, and also particularly appealing to clinicians thanks to its simplicity (see figure 2).

The algorithm used for this demonstrator is a simple rigid transformation estimation, based on the ITK registration framework. Different visualisation modes to evaluate the results are proposed, including the checker-board, alternating parts of each image. As an example result, figure 3 shows the multimodal registration of two MRIs (T1 and T2 modes). Three representations have been implemented so far: side by side comparison, and alpha-blending or checker-board mixing of the two images.

Perspectives

Within the VPH Toolkit, the objective of such toolbox is to aggregate the different tools existing in the community, thus it should suit the needs of other already existing algorithms to be wrapped as MedINRIA-fusion plug-ins. Moreover, by release a simple API and the software, the objective is also to attract more developers to the framework, so that any researcher can create his own plug-in.

Regarding the functionalities, the horizontal data fusion is also including other data supports, like meshes, so we plan to integrate fusion between meshes or between meshes and images in the next release.

Figure 2: Medinria-fusion: viewer interface, while using the registration tools attempt to federate researchers or developers under the same common language. More demonstrators are however needed to harden the concept.

Figure 3 :

 3 Figure 3: Checker-board fused MRIs (T1 and T2)

  Toussaint is now with the King's College London, Division of Imaging Sciences, St Thomas' Hospital The proposed framework can also help quantitative and qualitative tests of functionally similar methods. Medinria-fusion is a demonstrator that implements a registration algorithm, using the API we propose as an open framework. This software allows to visualise several formats of images (among them is DICOM).

1

http://www.vph-noe.eu/ 2 http://www-sop.inria.fr/asclepios/projects/noe/ change them between applications supporting the same framework. That attempt at a universal solution is a step towards standardisation of the exchange of algorithms. Standardisation calls for sustainability, maintainability.

http://www.itk.org

http://www.vtk.org

ConclusionStill in its infancy, this versatile framework showed an ease to import an existing algorithm, as well as an easy-to-use interface. The open API is our first


Natasha Leporé

Anand A Joshi 1⋆

Richard M Leahy

Caroline Brun

Yi-Yu Chou

Xavier Pennec

Agatha D Lee

Marina Barysheva

Greig I De Zubicaray

Margaret J Wright

Katie L Mcmahon

Arthur W Toga

Paul M Thompson

A New Combined Surface and Volume Registration

Keywords: brain, image analysis, Magnetic Resonance Imaging, image registration

3D registration of brain MRI data is vital for many medical imaging applications. However, purely intensitybased approaches for inter-subject matching of brain structure are generally inaccurate in cortical regions, due to the highly complex network of sulci and gyri, which vary widely across subjects. Here we combine a surfacebased cortical registration with a 3D fluid one for the first time, enabling precise matching of cortical folds, but allowing large deformations in the enclosed brain volume, which guarantee diffeomorphisms. This greatly improves the matching of anatomy in cortical areas. The cortices are segmented and registered with the software Freesurfer. The deformation field is initially extended to the full 3D brain volume using a 3D harmonic mapping that preserves the matching between cortical surfaces. Finally, these deformation fields are used to initialize a 3D Riemannian fluid registration algorithm, that improves the alignment of subcortical brain regions. We validate this method on an MRI dataset from 92 healthy adult twins. Results are compared to those based on volumetric registration without surface constraints; the resulting mean templates resolve consistent anatomical features both subcortically and at the cortex, suggesting that the approach is well-suited for cross-subject integration of functional and anatomic data.

INTRODUCTION

Registration of brain MRI scans is a key step in many medical imaging studies, for multimodality integration, computational anatomy studies, and longitudinal tracking of disease. In non-linear image registration, a template image T is transformed into a study image S, by aligning the geometrical features of the images. The anatomical correspondences between images can be measured and optimized in two ways, either based on intensity-derived information at each voxel such as the squared difference in intensities or mutual information between the images, or by using anatomical landmarks, which may include 3D parametric curves and surfaces extracted from the images. These feature correspondences can be used to drive the registration. The image transformation is also typically constrained via a regularizer that is chosen to ensure that the transformation is spatially smooth, invertible, and so that it satisfies other desirable properties such as inverse-consistency or transitivity. [START_REF] Leow | Statistical properties of Jacobian maps and inverse-consistent deformations in non-linear image registration[END_REF] Registration algorithms for brain MR images may also be divided into two major categories based on whether they attempt to register the entire 3D brain volume, or are restricted to the cortical surfaces. Volume-based algorithms perform quite well for subcortical structures, but have difficulty in matching cortical surfaces, due to the high complexity and variability in cortical patterns. Cortical registration is therefore used when the subsequent analysis focuses more specifically on the cortex, e.g. for population-based studies of cortical thickness. In this case, the registration is restricted to the cortical surface, and does not take into account the rest of the brain volume.

⋆: equal contribution

Cortical registrations are usually performed by first transforming the cortices to more regular surfaces such as the plane or the sphere, though some recent algorithms have been developed that directly register the cortices using 3D harmonic maps. [START_REF] Shi | Direct cortical mapping via solving partial differential equations on implicit surfaces[END_REF] Typically the surface-based registrations either use sulci or gyri as anchors to drive the transformation, [START_REF] Thompson | Genetic influences on brain structure[END_REF] or alternatively, make use of quantities defined over the whole brain such as the curvature [START_REF] Dale | Cortical surface-based analysis. I. Segmentation and surface reconstruction[END_REF] or the local conformal factor. [START_REF] Wang | Automated Surface Matching using Mutual Information Applied to Riemann Surface Structures[END_REF] To compute the flows in surface coordinates, specialized approaches using covariant PDEs, currents, or holomorphic 1-forms have been developed. [START_REF] Durrleman | Inferring Brain Variability from Diffeomorphic Deformations of Currents: An Integrative Approach[END_REF] There are many ways to perform the volume-based registrations. In this work, we focus on regularizers that assign continuum mechanical laws to the deforming image medium. The most popular model for these regularizers is the elastic one. Here a similarity function is selected, and a force is derived from it and applied to the system, while the elastic forces attempt to restore the original shape of the medium. Elastic registration is a good choice for studies where small deformations are expected, such as longitudinal studies where the same subject is imaged at two time points, and small changes are expected. However, in cases where large differences exist between the two images, the deformation can induce tearing or shearing of the elastic medium. In this case, other methods are needed. Among these, the fluid registration method was first introduced in. [START_REF] Christensen | Deformable templates using large deformation kinematics[END_REF] In this case, the image is treated as a viscous fluid that obeys a linearized version of the Navier-Stokes equation. The velocity v of the fluid is the time derivative of the deformation field v( r, t) = d u( r,t) dt . The driving force for the flow is again generated by a similarity function. The fluid paradigm allows large deformations without shearing or tearing of the image.

A new elastic regularizer was recently introduced in. [START_REF] Pennec | Riemannian elasticity: A statistical regularization framework for non-linear registration[END_REF] In this case, the regularization is based on the deformation tensors Σ = (∇ (u) + I) T (∇ u + I), with I the identity matrix. The Σ's characterize the voxelwise changes in shape and volume caused by the registration. As the Σ's are positive-definite symmetric matrices, the computations are performed in the Log-Euclidean framework, 1 which allows analytical computations on the manifold of Σ's. A fluid version of this method was developped in, [START_REF] Brun | A Tensor-Based Morphometry Study of Genetic Influences on Brain Structure using a New Fluid Registration Method[END_REF] where the rate of change of Σ was regularized instead of Σ.

In this work, we combined surface and volume based algorithms, to overcome serious weaknesses of each component code that severely limit their usefulness (one is restricted to the surface, the other cannot handle the cortex accurately). Although combination of the two is likely to be much more useful, as it allows good matching everywhere, no hybrid registration approach is yet widely used, [START_REF] Thompson | A Surface-Based Technique for Warping 3-Dimensional Images of the Brain[END_REF] . [START_REF] Thompson | Mathematical/Computational Challenges in Creating Population-Based Brain Atlases[END_REF] There are several ways to achieve a good registration over the whole brain volume, including the cortex. The simplest one, which we investigate here, is to use a cortical surface registration to obtain an initialization for the volume registration. This allows a slightly simpler implementation than one that allows the surface deformation to evolve while the volumes are matched. A similar method was designed in, [START_REF] Joshi | Surface-Constrained Volumetric Brain Registration using harmonic mappings[END_REF] where a new algorithm was introduced that performs both surface and volume registrations (see also [START_REF] Postelnicu | Combined volumetric and surface registration[END_REF] ). A landmark-based registration is first performed on the cortical surfaces, by constraining the matching between sulci. The correspondence was then extended to the whole brain using a constrained 3D harmonic mapping. The third step consists of refining the volume alignment via elastic registration.

Here we improve on the method from 13 in two ways. The algorithm in [START_REF] Joshi | Surface-Constrained Volumetric Brain Registration using harmonic mappings[END_REF] gives excellent matching in cortical areas. However, an elastic registration is used to align the subcortical structures. This method gives poor results when large deformations are required, and is not provably diffeomorphic, i.e., there is no mathematical guarantee that the interior mapping is invertible with a smooth inverse. This latter property is vital for tensorbased morphometry, a powerful approach for large-scale analysis of regional brain volumes and shape using statistics of the local deformation, 17 . 11 Here we combine a cortical registration with a fluid one for the first time, enabling precise matching of cortical folds, but allowing large deformations in the enclosed brain volume, which guarantee diffeomorphisms. For this, we will use the Riemannian fluid algorithm developed in. [START_REF] Brun | A Tensor-Based Morphometry Study of Genetic Influences on Brain Structure using a New Fluid Registration Method[END_REF] The Riemannian fluid is regularized to minimize the distortion at each voxel from the registration. Furthermore, the surface matching algorithm designed in [START_REF] Joshi | Surface-Constrained Volumetric Brain Registration using harmonic mappings[END_REF] uses sulci as landmarks to anchor the cortical registration. This method is prohibitive for large datasets, as the sulci need to be traced by hand. Here instead we use the Freesurfer software 7 (which makes use of a curvature-based registration algorithm) to align the cortical surfaces. As a result, the entire registration process is automated.

We validate our new combined algorithm using a large dataset of brain MRI images from 92 healthy adult subjects, divided into 23 pairs of monozygotic twins (MZ) and 23 same-sex pairs of dizygotic twins (DZ) from which we also obtain whole brain correlation statistics assessing the anatomical similarity between twins.

METHOD

Surface-Volume Registration

Our method aims to produce a one-to-one mapping between two brain volumes such that the subcortical structures and the cortical surfaces are aligned. Formally, let M and N represent the subject and target brain white matter volumes and let S and T denote the volumetric images such that S(x), x ∈ M , is the image intensity value at point x in the subjects volume M and T (x), x ∈ N , is the image intensity value at point x in the target volume N . The cortical surface ∂M of the subject M is the boundary of the subject white matter volume M and the cortical surface ∂N of the subject N is the boundary of the target brain's white matter volume N . There are two mains steps in our algorithm. The first step extends the surface-based registration to the enclosed 3D volume by using a 3D harmonic mapping based method [START_REF] Joshi | Surface-Constrained Volumetric Brain Registration using harmonic mappings[END_REF] and the second step performs the alignment of subcortical volumes by Riemannian fluid 4, 5 intensity registration.

Our method consists of:

1. Surface matching: A map is computed between ∂M and ∂N , the cortical surfaces of the subject and target brains, respectively. The registrations are done from each of the subjects cortices to that of the target brain.

2. Extrapolation of the surface map to the entire enclosed cerebral volume, such that the cortical surfaces remain aligned: This is done by computing a 3D harmonic map u from M to N while prescribing the deformation field defined on the cortices as a boundary condition, as in. [START_REF] Joshi | Surface-Constrained Volumetric Brain Registration using harmonic mappings[END_REF] The map is computed as follows (see [START_REF] Joshi | Surface-Constrained Volumetric Brain Registration using harmonic mappings[END_REF] for more details). First, M and N are both mapped to the surface of the unit sphere using Freesurfer. [START_REF] Dale | Cortical surface-based analysis. I. Segmentation and surface reconstruction[END_REF] Let g αβ and h αβ denote the metrics associated with the mapping of M and N to the unit sphere, respectively. The harmonic map u : M → N can be computed by minimizing the mapping energy:

E(u) = 1 2 M 3 i,j=1 3
α,β=1 h αβ (u(x)) ∂u α (x) ∂x i ∂u β (x) ∂x j dµ g . (1)
As a result of these two steps, we get a diffeomorphism between the two brain volumetric manifolds M and N . The map aligns the folding patterns of the sulci but does not align the subcortical structures because the intensity values are not used in this step. This map is refined further in the next step using Riemannian fluid-based intensity registration. [START_REF] Brun | A Tensor-Based Morphometry Study of Genetic Influences on Brain Structure using a New Fluid Registration Method[END_REF] A brief description of the method can be found in the following section.

Volume Registration

We register the outputs of the previous stage further for subcortical alignment using the Riemannian fluid registration algorithm. [START_REF] Brun | A Tensor-Based Morphometry Study of Genetic Influences on Brain Structure using a New Fluid Registration Method[END_REF] Fluid registration proceeds in a series of time steps dt. At each time t and voxel x, the velocity v is found from:

d v( x, t) dt = ∇Cost + ∇Reg( v, t) -v. (2)
Cost is the similarity cost function and Reg is the regularizer. Once the registration is completed, v is integrated over time to obtain the final displacement field u.

In particular, for the Riemannian fluid registration, 5 the regularization is done on the matrix logarithm of the rate of strain

Σ v = (∇ v + I) T (∇ v + I): Reg Riem ( v, t) = µ 4 T r(log(Σ 2 v )) + λ 8 T r(log(Σ v )) 2 . ( 3
)
Here µ and λ are fluid parameters to be selected by the user (analogous to the Lamé coefficients). Reg Riem is the fluid equivalent of the Riemannian elastic regularizer from; 21 the fluid case regularizes the rate of strain Σ v rather than the deformation tensors Σ u = (∇ u + I) T (∇ u + I) in Eq. 2. The tensor field Σ characterizes the local shape and volume deformation. These tensors live on the curved manifold of positive-definite, symmetric matrices. Taking the matrix logarithm of Σ transports it to the tangent plane at the origin of the manifold, which is a flat space (vector space) in which simple Euclidean computations can be performed. [START_REF] Arsigny | Log-Euclidean metrics for fast and simple calculus on diffusion tensors[END_REF] We chose an intensity-based cost function that penalizes the squared difference in intensity between the deforming template T and the study S:

Cost( x, u( x, t)) = (T ( x -u( x, t)) -S( x)) 2 (4)
We note that intensity matching terms based on mutual information, correlation, or other metrics could be used instead, but we do not cover these here. See 4 for details of the implementation.

As a result of this step, the diffeomorphic map generated in the previous step is refined to align the subcortical structures.

Non-linear Registration

We performed all 3 of the steps described in section 2.1 between each image in the dataset and a template, which was chosen randomly as one of the subjects images. Alternatively, as in our past work, a geometrically-centered mean template could be generated and used instead. [START_REF] Leporé | Mean template for Tensor-Based Morphometry using deformation tensors[END_REF] For each subject, we obtained a displacement field u 1 and a registered image from the surface matching and the subsequent harmonic mapping. The displacement field u 1 was used as input for the Riemannian fluid algorithm described in sec 2.2. The volume registration yielded a vector field u 2 . We used the vector fields from the surface registration as input instead of the registered images. A diffeomorphic map between the surfaces (such as those from Freesurfer and the harmonic extension to the whole brain volume) is provided as an initial input for the fluid registration; this also ensures that the final results will also be diffeomorphic, as the composition of two diffeomorphic maps. [START_REF] Christensen | Deformable templates using large deformation kinematics[END_REF] In our implementation, the final displacements were computed from the common template to each image in the data set, while the registered image was the result of transforming each subject into the template, thus allowing all comparisons to be made in the common space.

RESULTS

Data

We analyzed 3D structural brain MRI scans of 23 monozygotic (MZ) and 23 dizygotic (DZ) same-sex twin pairs (age range: 22 -25 years), as well as one of an identically scanned healthy subject used as the common target for the fluid registration. Our images were collected using a 4 Tesla Bruker Medspec whole body scanner (Bruker Medical, Ettingen, Germany) at the Center for Magnetic Resonance (University of Queensland, Australia). Three-dimensional T 1-weighted images were acquired with a magnetization prepared rapid gradient echo (MP-RAGE) sequence to resolve anatomy at high resolution. Acquisition parameters were: inversion time (T I) /repetition time (T R) /echo time (T E) = 1500 / 2500 / 3.83 msec; flip angle = 15 o ; slice thickness = 0.9 mm with a 256x256x256 acquisition matrix.

Extracerebral (non-brain) tissues were manually deleted from the MRI images using the 3D interactive program Display (Montreal Neurological Institute, McGill University, Canada). All scans were then aligned to the ICBM53 template using 9-parameter registration (i.e., translational and rotational alignment, allowing scaling in 3 independent directions) found in the FMRIB's Linear Image Registration Toolbox, FLIRT. [START_REF] Jenkinson | Improvised optimization for the robust and accurate linear registration and motion correction of brain images[END_REF]

Twin Statistics

For each subject, at each voxel, the displacement u was computed, as well as the deformation matrix Σ and the determinant of the Jacobian matrix J = (Id + ∇u). det(J) measures local volume expansions, where det(J) > 1, and shrinkage when det(J) < 1. In the full set of MZ and DZ twins, the resemblance between twin pairs at each voxel can be computed from det(J) to assess the Intraclass Correlation:

23 ICC = σ 2 b (σ 2 b +σ 2 w )
, where σ 2 b and σ 2 w are the pooled variance between pairs and within pairs, respectively. We did not want to assume that the data on volumetric differences (here det(J)) were normally distributed, so we computed p values at each voxel with a voxelwise permutation test (5000 permutations), to establish a null distribution for the ICC statistics at each voxel (the null hypothesis was ICC = 0 (no correlation)).

Details of the Numerical Implementation

We used Freesurfer for segmentation of the brains and extraction of the cortical surface meshes. The surface registration step in Section 2.1 was performed by Freesurfer by using a spherical parameterization of the cortex and a metric representing cortical 'depth' 7 that is a function of the surface curvature. The extension of the surface registration to the full volume, described in Section 2.1, and the Riemannian fluid algorithm, in Section 2.2, were implemented in MATLAB. To test our method, we selected one brain from the database as a target brain volume and registered the rest of the brains to it. On an Intel Pentium IV processor the registration of each subject to a target image takes approximately 1 hour.

Registration Results

Figure 1 shows a comparison of the registration from the fluid code alone, without cortical surface constraints, compared to the combined code, which also matches the cortex. The fluid registration alone gives poor results in some areas, e.g., most of the cortex. This is a serious problem for all fluid algorithms, hence the initialization with the cortical registration used here. Initializing with the displacements from the cortical registration noticeably improves the cortical results. Fig. 1 also shows that the subcortical structures are now also better aligned; this may be due to deep cortical surfaces passing near the them, especially at the corpus callosum.

The top row in Figure 2 shows the results of averaging the intensities of all 92 registered images. For the results using the proposed hybrid algorithm, features at the cortical boundary and the ventricles are still clearly defined after averaging. This is in marked contrast with those from the Riemannian registration alone, where many key anatomical features are blurred away by the averaging. We can see that alignment of the surface gives sharper images near the surface boundary. We also observed less blurring for the subcortical structures such as the ventricles. The second row show the variance in intensity at each voxel for the images in our dataset. The variance is much lower for the combined algorithm. Again this is true even in subcortical areas, except within the ventricles.

DISCUSSION

We investigated a new combined algorithm for the registration of brain MR images. Our algorithm was tested on a large data set of 92 healthy adult twins, and showed significant improvement over using fluid registration alone, both for subcortical areas and in the vicinity of the cortex. Intraclass correlation statistics were significantly improved over using fluid registration alone. The premise of using twin data to help validate registration is that any errors in registration will tend to be uncorrelated between twins, and will tend to deplete the natural correlation in brain anatomy that is due to genetic factors.

There are several ways to achieve a good registration over the whole brain volume, including the cortical surface. The simplest one, which we investigated here, is to use a cortical surface registration to initialize a 3D registration. In that case, as we showed here, the matching converges to better local minimum compared to the one found from the fluid registration alone. However, as the cortex is allowed to move during the fluid registration, it can potentially reduce the accuracy of the matching found from the cortical matching algorithm. Thus, as an alternative, one could choose instead to re-introduce the cost function from the first cortical matching algorithm into the volume registration, either as a hard constraint or with a Lagrange multiplier. As the matching was already quite good with the first method, this second method was not pursued here, though it would be an interesting extension of this work. Finally, in cases where both the above methods fail, a third option would be to force the displacements from the Riemannian fluid algorithm to match those found from the cortical mapping.

As a next step, a large-scale group validation is needed to compare the performance of our algorithm to other state-of-art registration algorithms. In the near future, we plan to use a validation such as the one in [START_REF] Klein | Evaluation of 14 nonlinear Deformation Algorithms Applied to Human Brain MRI Registration[END_REF] to further validate our method. Because our approach produces diffeomorphic mappings while enforcing cortical constraints, it is likely to improve results in tensor-based morphometry studies, which are now being conducted on a large scale to assess the effects of disease, medication, and specific genetic polymorphisms on brain structure.

.

Figure 1 .

 1 Figure 1. Comparison of the results of the combined algorithm to those of the 3D Riemannian algorithm with no surface constraints. Registration results are shown mapping one subject to a template in coronal (left column), sagittal (middle column) and horizontal (right column) slices. The top row shows the moving image and the second row shows the fixed target image. The last two rows show the results of the registration for the volumetric registrations without surface constraints (third row) and with constraints (fourth row); the fourth row shows accurate cortical and subcortical registration. Results are more accurate for the combined code both on the cortex (see boxed area in figure) and for subcortical structures.

Figure 2 .

 2 Figure 2. The top row shows the result of averaging the intensities for all 92 registered images, with the colorbar showing the mean intensity values. Anatomical features near the cortex are much sharper when the combined code is used. The bottom row shows the voxelwise variance in intensity in these same registered images, with the variance values shown according to the colorbar on the right. The top and bottom left images are for our combined algorithm, while the ones on the right are for the Riemannian algorithm alone. The intensity variance in the registered images is much lower with the combined algorithm, which suggests that misalignment of cortical anatomy is reduced.

Figure 3 .

 3 Figure 3. Maps of p-values for the Intraclass correlation for the monozygotic twins. Sulcal patterns are known vary naturally even among monozygotic twins, but subcortical anatomy is under strong genetic control; areas of highly significant correlation include the anterior cingulate and paralimbic regions.

Acknowledgments

This work was supported by NIH grant R01 HD050735 and the National Health and Medical Research Council, Australia grant 496682.