Argheesh Bhanot
email: bhanot@unistra.fr

Céline Meillier
email: meillier@unistra.fr

Fabrice Heitz
email: fabrice.heitz@unistra.fr

Laura Harsan
email: harsan@unistra.fr

C Meiller

Spatially Constrained Online Dictionary Learning for Source Separation

Keywords: Source separation, dictionary learning, fMRI, scintigraphy, hyperspectral data

Whether in medical imaging, astronomy or remote sensing, the data are increasingly complex. In addition to the spatial dimension, the data may contain temporal or spectral information that characterises the different sources present in the image. The compromise between spatial resolution and temporal/spectral resolution is often at the expense of spatial resolution, resulting in a potentially large mixing of sources in the same pixel/voxel. Source separation methods must incorporate spatial information to estimate the contribution and signature of each source in the image. We consider the particular case where the position of the sources is approximately known thanks to external information that may come from another imaging modality or from a priori knowledge. We propose a spatially constrained dictionary learning source separation algorithm that uses e.g. high resolution segmentation map or regions of interest defined by an expert to regularise the source contribution estimation. The originality of the proposed model is the replacement of the sparsity constraint classically expressed in the form of an 1 penalty on the localisation of sources by an indicator function exploiting the external source localisation information.

The model is easily adaptable to different applications by adding or modifying the constraints on the sources properties in the optimisation problem. The performance of this algorithm has been validated on synthetic and quasi-real data, before being applied on real data previously analysed by other methods of the literature in order to compare the results. To illustrate the potential of the approach, different applications have been considered, from scintigraphic data to astronomy or fMRI data.

I. INTRODUCTION

The issue of source separation, or unmixing, is well known to the signal and image processing community. It concerns a very large number of applications and can occur under different conditions of source mixing. A large part of the literature is devoted to blind source separation (BSS) [START_REF] Wold | Principal component analysis[END_REF]- [START_REF] Comon | Independent component analysis, a new concept?[END_REF]. BSS methods allow to solve cocktail party problems for which P signals (or images) composed of a mixture of R sources are observed, without any a priori on the properties of the sources. If the number of sources involved in the mixtures is not known a priori, then they must also be estimated [START_REF] Bioucas-Dias | A variable splitting augmented Lagrangian approach to linear spectral unmixing[END_REF], [START_REF] Šmídl | Automatic regions of interest in factor analysis for dynamic medical imaging[END_REF]. Early BSS methods mainly comprised Independent Component Analysis (ICA) and sparse decomposition analysis.

Many variants of the ICA approach have been proposed in the literature to solve BSS problems. All of them are based on the general principle of spatial independence of the sources, which makes it possible to estimate their temporal (or spectral) signatures.

For instance, in brain functional networks detection in functional Magnetic Resonance Imaging (fMRI) data, ICA is widely used to separate spatial sources by assuming the independence of the temporal signals associated with each spatial source, i.e. functional network. Spatial ICA has proven effective in [START_REF] Mckeown | Spatially independent activity patterns in functional MRI data during the Stroop color-naming task[END_REF], [START_REF] Xu | Spatial ICA reveals functional activity hidden from traditional fMRI GLM-based analyses[END_REF] for fMRI data, but the main drawback of ICA approaches is the unknown number of sources which is set arbitrarily and may lead to a large number of nuisance sources, that must be screened manually or by a semi-automatic method [START_REF] Sourty | Towards an automated selection of spontaneous coactivity maps in functional magnetic resonance imaging[END_REF].

In contrast to the BSS problem, many unmixing problems involve a dictionary of pre-defined sources.

For example, in hyperspectral imaging for remote sensing, libraries of light spectra corresponding to the different materials that may be observed in the scene are available, so only the proportion of the different materials in each pixel is estimated. Between these two extreme cases, there are a large number of unmixing problems where some information on the form or location of the sources or the type of mixture is known. Sum-to-one and positivity constraints on the coefficients of the mixing matrix are classic in signal and image processing. In remote sensing applications, hyperspectral data linear unmixing is for instance carried out by methods based on nonnegative matrix factorization [START_REF] Yang | Blind spectral unmixing based on sparse nonnegative matrix factorization[END_REF], [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF]. In recent years, sparse decomposition methods have been widely used to solve source separation problems. The sparsity constraint is another way to reduce the set of solutions. It can be combined with the two latter constraints.

The sparsity may concern the mixing itself, i.e. for a given observed signal, only a few number of sources is involved, or the decomposition of the sources on a dictionary (wavelet, discrete cosine transform, or custom atoms) [START_REF] Abolghasemi | Blind separation of image sources via adaptive dictionary learning[END_REF]. Recently for fMRI application, sparse analysis based on dictionary learning methods have proven to be promising [START_REF] Abraham | Extracting brain regions from rest fMRI with total-variation constrained dictionary learning[END_REF]- [START_REF] Mensch | Compressed online dictionary learning for fast resting-state fMRI decomposition[END_REF]. Dictionary learning methods take into account the spatial sparsity of the functional networks in the form of 1 constraints on the mixing matrix in the minimisation problem.

In this paper, we are interested in the problem of source separation in spatially structured data: 2D or 3D images that contain temporal information (fMRI, scintigraphy) or light spectrum information (hyperspectral imaging). We consider two categories of this kind of unmixing problems. In the first one, for a given pixel/voxel, different sources contribute to the mix in the sense that the spatial resolution is not fine enough to allow spatial separation of the sources. This is the case, for example, with unmixing problems in remote sensing [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF] or fMRI applications [START_REF] Mensch | Compressed online dictionary learning for fast resting-state fMRI decomposition[END_REF], [START_REF] Bhanot | Online dictionary learning for single-subject fMRI data unmixing[END_REF]. In this case, the mixing matrix is a matrix of proportions where for a given pixel, the sum of the contributions of each source is equal to 1.

In the second category of unmixing problems, the mixing is additive, the signals of the different sources are superimposed and their sum forms the observed mixing signal. Decomposition of scintigraphic image sequences into tissue images and their time-activity curves, or unmixing of light sources in hyperspectral data in astronomy are examples of this second category of problems. In this case, the sum-to-one constraint is not relevant, the coefficients of the mixing matrix are the intensity of the contribution of each source in the mixture. Since the observed signals are observed in the form of images, the constraints that can be defined in the optimisation problem should be related to the location of the sources and not to their shape. We consider the case where no information on the temporal or spectral signature of the different sources, or on their dependence is available.

We propose a dictionary learning method that introduces sparsity constraints on the spatial localisation of sources from external knowledge. Additional constraints on the mixing matrix (positivity and sum-toone constraints) can be added or removed depending on the application. To illustrate the potential of the approach, different applications have been considered, from scintigraphic data to astronomy or fMRI data.

All these data are of very different natures, as well as the a priori information available on the location of the sources. We thus show that our algorithm is adaptable to different types of data and different types of a priori knowledge on the location of sources. In the case of multimodal observations, information regarding the possible spatial location of sources is usually derived from a high spatial resolution image that does not provide the second dimension, namely temporal or spectral information. Unlike multimodal image fusion problems, such as pansharpening [START_REF] Loncan | Hyperspectral pansharpening: A review[END_REF], our goal is not to produce a spatially and spectrally or temporally well-resolved image. We rather aim at exploiting a segmentation information from a high spatial resolution image in order to improve the unmixing of spectral or temporal sources at a lower resolution image level. In some applications, such a segmentation map is not available, but approximate spatial location information can be provided by an expert who can define regions of interest (ROIs) (see for example the unmixing method for highly realistic simulated renography dataset in [START_REF] Filippi | Robust unmixing of dynamic sequences using regions of interest[END_REF], [START_REF] Benali | Foundations of factor analysis of medical image sequences: a unified approach and some practical implications[END_REF]). This paper is organised as follow: section II presents the observation model and the dictionary learning approach for unmixing spatially structured data based on the spatial sparsity constraint. The performance and adaptability of the algorithm to different applications are illustrated on synthetic data in section III and on real data in section IV. Comparison of our approach with application-specific state-of-the-art methods are proposed on three different applications: unmixing in scintigraphic image sequences, fMRI and hyperspectral astronomic datacubes.

II. MODEL AND METHOD

The classical linear model used in source separation may be written as:

Y UA, (1)
where Y ∈ R N ×P is the observed data. In spatio-spectral separation, N can be interpreted as the spectral length and in spatio-temporal separation N is the length of the temporal signals. P is the number of voxels or pixels, depending upon the dataset. Matrix U ∈ R N ×R contains the temporal/spectral signatures where R is the number of sources. Matrix A ∈ R R×P , usually called mixing or abundance matrix, codes the fraction of the R components contributions at each voxel or pixel.

A. Constrained optimisation formulation

Given the observation model ( 1), the following minimisation problem:

min A,U 1 2 Y -UA 2 F ( 2
)
does not have a unique solution because of the joint estimation of A and U, and the ill-posedness of the problem. In order to restrain the number of solutions, we introduce some standard constraints on matrix

A such as the positivity constraint A ∈ R + and the sum-to-one constraint A[., i] = 1, with i = 1 : P , as these are the proportions voxel by voxel. The form of matrix A may also be constrained by some extra information from high resolution (HR) segmentations of sources or source locations knowledge:

we know which ROIs may contribute to a given voxel, i.e. present a non-zero proportion at this voxel.

When the total number of sources R is high, this a priori knowledge allows to constrain the possible solutions of the minimisation problem. In the most general setting, the unmixing problem is recast as:

min A,U 1 2 Y -UA 2 F + µ σ 2 U 2 F + I R + (A)+ I S (A) + I M ( Ã) (A), (3)
where the first term is the data fidelity term, the second term is a Tikhonov regularisation controlled by parameter µ σ set to 10 -4 to prevent bad conditioning (see section II-B). The third term is a positivity constraint where I R + (A) = ∞ if at least one of the elements of A is negative, and 0 otherwise. The fourth term in eq. ( 3) codes an optional sum-to-one constraint on each column of matrix A, I S (A) = ∞ if at least one column of A does not sum to one, and 0 otherwise. The last term I M ( Ã) (A) is the indicator function on the set M ( Ã) of matrices having a structure similar to a given binary "structure matrix" Ã, i.e. A ∈ M ( Ã) if and only if A ∈ R R×P and coefficient A i,j = 0 if Ãi,j = 0. Ã is a binary matrix, where element ( Ã) r,i = 1 if, according to a priori knowledge about spatial localisation of the sources, the r th region of interest could exist in the i th voxel, and 0 otherwise. This results in I M ( Ã) (A) = ∞ if at least one element of A is non-zero while it is zero in Ã, and 0 otherwise.

Combining sets R + ∩ S = S + , the optimisation problem can be rewritten as follows:

min A,U 1 2 Y -UA 2 F + µ σ 2 U 2 F +I S + (A) +I M ( Ã) (A). (4)
Estimating jointly U and A in eq. ( 4) is a typical problem of dictionary learning (DL). But, unlike conventional DL algorithms, there is no sparsity regularisation term in the form of an 1 penalty: it is the sources localisation information coded in the structural term I M ( Ã) (A) which enforces the sparse decomposition of each voxel. A classical way to solve the joint estimation problem is to optimise alternatively the cost function eq. ( 3) along U and A.

B. Estimation of the temporal / spectral signatures matrix U

Considering that A is fixed, problem ( 3) becomes:

min U 1 2 Y -UA 2 F + µ σ 2 U 2 F . (5)
The Tikhonov regularisation term µσ 2 U 2 F is introduced to improve the conditioning of problem ( 2), µ σ is set to 10 -4 to prevent collinearity between columns of U.

The solution of ( 5) is the ridge estimator defined by:

Û = YA T (AA T + µ σ I R ) -1 , (6)
where I R is the R × R identity matrix.

C. Estimation of the abundance / mixing matrix A

Consider that U is fixed, then problem eq. ( 3) becomes min A f (A), where:

f (A) = 1 2 Y -UA 2 F +I S + (A) +I M ( Ã) (A). (7)
Note that this function is separable according to the pixels/voxels i ∈ {1,. . . , P}, which leads to:

min ai f (a i ), with: f (a i ) = 1 2 y i -Ua i 2 F +I S + (a i ) +I M (ãi) (a i ), (8)
where a i is a column vector from the matrix A (and with an abuse of notation for f (.)). The set of all the vectors with a structure similar to a i is given by ãi , where ãi is a column of Ã. The regularisation terms in eq. ( 8) can be summarised as:

g(a i ) = I M (ã)∩S + (a i ). (9)
Note that the objective function eq. ( 8) is convex since the first term is convex and differentiable and g(a i ) is convex but non differentiable. The proof of convexity for g(a i ) is given in Appendix A. Minimisation of the objective function given by eq. ( 8) belongs to the class of problems on which the proximal gradient methods can be applied. Different algorithms are available, for example, alternating direction method of multipliers (ADMM) [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF], projected gradient, also known as iterative shrinkage-thresholding algorithm (ISTA) or FISTA (Fast ISTA) [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]. Algorithm FISTA was preferred for its rapid convergence: its implementation is given in Algorithm 1.

1 Initialisation of A (0) , ω (1)
2 while STOPPING CRITERIA = TRUE do 3 Minimisation problem w.r.t U 4 U (l+1) = YA (l) T (A (l) A (l) T + µ σ I R ) -1 5
Parallel minimisation w.r.t. the columns a i of A by employing FISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] 6

for k ← 1 to proxsteps do 7 a (k) i = prox g (ω (k) -λ∇f (ω (k) )) 8 t (k+1) = 1+ √ 1+4(t (k) ) 2 2 9 ω (k+1) = a (k) i + t (k) -1 t (k+1) (a (k) i -a (k-1) i ) 10 end 11 a (l+1) i = a (proxsteps) i 12 end
13 return A, U Algorithm 1: Alternate optimisation algorithm to estimate A and U

In Algorithm 1, ∇f (a i ) is the gradient of f (a i ), given by U T (Ua iy i ). The step size λ is set equal to the inverse of the Lipschitz constant of ∇f (a i ) i.e. 1/L, where L = U T U F . t (k+1) is an auxiliary variable which helps in the fast convergence of FISTA, ω calculates intermediate values based on a special linear combination of the last two points and prox refers to the proximal operator [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]. In our case, the proximal operator is just the projection of a i in the positive orthant, with the vector normalised to-sum-to-one. This projection also forces the elements of abundance matrix (A) r,i to be non-zero only at positions where the region of interest r projects on pixel/voxel i (I M (ã) constraint). The proximal operator of the function g is:

prox g (y) = argmin x∈M (ã)∩S + x -y 2 = P M (ã)∩S + (y), ( 10
)
where P is the projection operator on set M (ã) ∩ S + . The orthogonal projection of a vector y ∈ R R on M (ã) ∩ S + is obtained using the projection onto convex sets (POCS) method [START_REF] Boyd | Alternating projections[END_REF]. POCS algorithm alternates projection onto the simplex S + = R + ∩ S and projection onto the set M (ã) of vectors having the same structure as ãi . Only a few iterations are required for convergence of the POCS algorithm.

Convergence towards a global minimum of DL algorithms cannot be proven. In practice, a good initialisation of A and the presence of pure pixels (as in remote sensing applications) in each region guarantee a good joint estimation of U and A. Previous work [START_REF] Bhanot | Online dictionary learning for single-subject fMRI data unmixing[END_REF] has demonstrated the importance of well-defining the spatial constraint on abundance I M ( Ã) (A) to ensure an acceptable estimate of abundances and spectral or temporal signatures.

III. EVALUATION ON SYNTHETIC DATASETS

In this section we evaluate the unmixing performance of our algorithm on two different synthetic datasets. Dataset I was created to show unmixing of signals/spectra taking into account different situations that could occur in real applications such as fMRI or astronomical data unmixing. Dataset II is an example of realistic synthetic data in scintigraphy used in [START_REF] Filippi | Robust unmixing of dynamic sequences using regions of interest[END_REF], for which the authors have proposed an unmixing method based on prior knowledge of the location of the regions of interest. This method, called Robust Unmixing of Dynamic Sequences Using Regions of Interest (RUDUR), has been compared in [START_REF] Filippi | Robust unmixing of dynamic sequences using regions of interest[END_REF] to 2) Algorithm details: The ground truth is given by the localisation map in Fig. 1a. To initialise A (0) , each region was dilated with a 7 pixels square structuring element and then the proportion for each region over each pixel was calculated, respecting the sum-to-one condition. The dilatation was done to introduce the uncertainty in the localisation of regions; as the localisation is seldom precise when dealing with real data. The algorithm used for unmixing is given in Algorithm 1. 400 steps were adopted for FISTA, in combination with 50 steps of alternate optimisation. The weighting parameter in the Tikhonov regularisation was set to 10 -4 as no more smoothing was required. A standard normalisation was applied to the data before processing:

different
Y (i) = Y (i) -µ (i) σ (i) Y
, where µ is the mean of the temporal signal of the i th pixel and σ is the standard deviation of the timecourse of the i th pixel.

3) Results and discussion: We observe that the timecourses and the abundances for the seven regions are well estimated even if the abundances are not perfectly initialised. In Fig. 2 the -20dB and -10dB curves converge to values less than 1 since some noise remains in the estimated timecourses. The 0 th iteration shows the ratio Y-UA F η F calculated with the initial enlarged A and U, estimated with least squares.

B. Dataset II

In an effort to objectively evaluate the performances of our approach, we propose to test and compare our method on a more realistic synthetic data set of scintigraphy images created for the evaluation of the performances of the state-of-the-art RUDUR method [START_REF] Filippi | Robust unmixing of dynamic sequences using regions of interest[END_REF]. We have reused the RUDUR code, as it is distributed by the authors [START_REF] Filippi | Robust unmixing of dynamic sequences using regions of interest[END_REF].

1) Data description:

This second dataset provides a first test case in which the method is confronted with a real application, on physical model-based simulations of scintigraphic images, with ground truth. It enables a comparison with a recent reference method in the field [START_REF] Filippi | Robust unmixing of dynamic sequences using regions of interest[END_REF]. Region of interest-based unmixing methods are common in scintigraphy [START_REF] Šmídl | Automatic regions of interest in factor analysis for dynamic medical imaging[END_REF], [START_REF] Benali | Foundations of factor analysis of medical image sequences: a unified approach and some practical implications[END_REF], [START_REF] Nijran | The importance of constraints in factor analysis of dynamic studies[END_REF]. All standard methods are source separation methods.

They estimate the time activity curves (TACs) and emissions of a tracer (a radioactive element) in the different body organs. A dataset of scintigraphic data has been made available at [START_REF] Šámal | Database of dynamic renal scintigraphy[END_REF]. This dataset is based on a Monte Carlo simulation of scintillation camera imaging [START_REF] Brolin | Dynamic (99m)Tc-MAG3 renography: images for quality control obtained by combining pharmacokinetic modelling, an anthropomorphic computer phantom and Monte Carlo simulated scintillation camera imaging[END_REF]. The datacube comprises images of size 21 × 26, with N = 60. The dataset, containing R = 3 regions, is shown in Fig. 3. In Fig. 3 the first row (a) shows the ground truths for the spatial maps. In the second row (b) the ROI initialisation is presented. This initialisation is not the same as in [START_REF] Filippi | Robust unmixing of dynamic sequences using regions of interest[END_REF] as our method needs strictly greater initial ROIs than the regions which need to be unmixed.

2) Algorithm details: RUDUR algorithm has a soft constraint on the source locations which allows the regions to unmix data even if the ROIs selected lie in the interior of the actual locations of the regions. In our model, we use a hard constraint on the locations of the regions, through the regularization term I M ( Ã) , so initial ROIs should be stricly enclosing the regions for which we want to estimate the timecourses. To achieve this, the binary mask of ROIs used in [START_REF] Filippi | Robust unmixing of dynamic sequences using regions of interest[END_REF] have been dilated with a 5 pixels square structuring element.

We should note that this application corresponds to a an additive case of unmixing, so the sum-toone constraint was dropped off in our algorithm. Further, as scintigraphy timecourses should be strictly positive (representing the emission of the tracer), eq. ( 3) was changed to:

min A,U 1 2 Y -UA 2 F + µ σ 2 U 2 F + I R + (A)+ I M ( Ã) (A)+ I R + (U).
Due to the addition of I R + (U), the constraint of positivity on the TACs, the least squares given in step 4 of Algorithm 1 to solve for U had to be replaced by FISTA steps to estimate U in each alternate step. The initialisation of the algorithm was done with the help of least squares using the initial dilated ROIs.

As the sum-to-one condition is dropped, the proximal operator of the function g in Algorithm 1 changes to:

prox g (y) = argmin x∈M (ã)∩R + x -y 2 = P R + ∩M (ã) (y), (11)
where P is now the projection operator on the set R + ∩ M (ã).

3) Results and discussion: We ran our algorithm with 500 steps of alternate optimisation. At each iteration l, convergence is monitored by the optimisation gain κ l defined as:

κ l = Y -U (l) A (l) F -Y -U (l-1) A (l-1) F Y -U (l-1) A (l-1) F . ( 12
)
which decreases to 10 -16 at the 500th alternate step. Estimated temporal signals and spatial maps were normalised by the criteria given in [START_REF] Filippi | Robust unmixing of dynamic sequences using regions of interest[END_REF] for comparison with the provided ground truth. The normalised mean squared error (NMSE) and normalised mean absolute error (NMAE) for the estimated spatial maps and time activity curves were calculated. In addition, we also provide spectral angle distances (SAD) for the estimated time activity curves. The spatial maps unmixed by our method are presented in Fig. 3 (c) and those estimated by RUDUR are displayed in Fig. 3 (d). Visually the maps are very similar although the errors are lower in the case of RUDUR. The last row(e) shows the TACs estimated by the proposed method (in blue). We observe in Fig. 3 that our solution is near the target solution (in dashed black) for all the sources and is close to the solution provided by RUDUR.

IV. APPLICATION ON SEMI-REAL AND REAL DATA

In this section we discuss the application of our algorithm to sources unmixing of (semi-real) resting state (rs) 3D fMRI data and on a real dataset of hyperspectral astronomic images. In these two applications, in addition to the 3D + time or hyperspectral data, a high resolution segmentation map is available as a standard, with a registration on the low resolution images containing the sources to be separated.

A. Semi real rs-fMRI dataset

Resting state functional Magnetic Resonance Imaging (rs-fMRI) has been widely used for studying brain functional connectivity [START_REF] Lee | Resting-state fMRI: a review of methods and clinical applications[END_REF]. Rs-fMRI allows the observation of changes in cerebral activity by analysing the blood-oxygen-level-dependent (BOLD) signal [START_REF] Logothetis | Interpreting the BOLD signal[END_REF]. At rest, only spontaneous activity is measured and a set of anatomical regions with the same fluctuations are considered part of a common network. We apply the algorithm on a single subject rs-fMRI. The motivation to work on single subject has been detailed in [START_REF] Bhanot | Online dictionary learning for single-subject fMRI data unmixing[END_REF]. Benefit of integrating a high resolution (HR) anatomical atlas in the single-subject case has also been demonstrated in this previous work.

In order to test the proposed unmixing method, we use data acquired in a preclinical study with an

Alzheimer mouse model at ICube laboratory under the IRIS platform, the studies on this data are given in [START_REF] Degiorgis | Longitudinal alterations of resting-state functional connectivity in alzheimer's disease in a tauopathy mouse model[END_REF], [START_REF] Degiorgis | MRI analysis of brain connectivity in a mouse model of Alzheimer's disease[END_REF]. The data consist of a 3D+t rs-fMRI and a 3D anatomical image registered to the rs-fMRI image.

The anatomical image has a dimension of 256 × 256 × 34 and 0.08299 × 0.07812 × 0.4 mm resolution. Mouse Brain Atlas [START_REF]A mesoscale connectome of the mouse brain[END_REF]. Each colour represents a label associated with an anatomical region.

1) Validation dataset:

A validation dataset is created by introducing a set of synthetic temporal signatures in seven small regions of the real Azheimer mice dataset. The seven regions, labelled ACAd1, ACAd5, ACAd6a, ACAv1, ACAv5, ORBl1 and PL1, have been arbitrarily chosen in the prefrontal cortex.

A first synthetic signal is obtained by averaging the real signals of the regions ACAd1 and ACAd5 which were already highly correlated in the real data. This signal is then modified to create signals with arbitrary high correlation or anti-correlation for the regions ACAd1, ACAd5, ACAd6a, ACAv1, ACAv5, ORBl1 and PL1 (see blue lines in plots of figure Fig. 6). These correlations do not have a physical significance, they are used as a ground truth for evaluation of the proposed algorithm performances [START_REF] Bhanot | Online dictionary learning for single-subject fMRI data unmixing[END_REF]. The seven synthetic signals are shown in Fig. 6.

The classical fMRI preprocessing pipeline of slice timing and co-registration is applied on the rs-fMRI dataset. Also, the confounding signals are regressed before analysing the data. The next step consists in registering the spatially well-resolved ABA template to the artificially augmented anatomical image (which is already perfectly aligned with the rs-fMRI data). The registration of the ABA mouse template to the anatomical images provides the deformation field that is applied to the HR segmentation map to transport the different labelled regions on the augmented rs-fMRI data. The spatial resolution of rs-fMRI data is augmented by subdividing each original voxel into 3×6×2 high resolution voxels. The registration of the anatomical image to the augmented rs-fMRI leads to an increase in its own resolution. Large Telescope, which produces hyperspectral observations of the deep sky. In these hyperspectral images, we can observe hundreds or even thousands of galaxies. Depending on their age, chemical composition, type, distance, etc, these galaxies have different spectra. These spectra may contain emission lines, continuous components, etc. One of the main objectives of MUSE data analysis is the detection of very distant galaxies, which therefore emit very low light flux. Spectrum of distant galaxies consists of a single emission line, the Lyman-α emission line, which is a marker of the strong presence of hydrogen in the galaxy. They are difficult to detect due to their distance and their very faint intensity compared to closer galaxies, and a strong noise affects the data. Moreover if two galaxies are aligned in the direction of observation, their spectra will blend inside pixels of MUSE images.

Recently two articles [START_REF] Bacon | The MUSE Hubble Ultra Deep Field Survey-I. Survey description, data reduction, and source detection[END_REF], [START_REF] Inami | The MUSE Hubble Ultra Deep Field Survey-II. Spectroscopic redshifts and comparisons to color selections of high-redshift galaxies[END_REF] were published around a MUSE dataset called Ultra Deep Field (UDF) which corresponds to an area of the sky previously observed by the Hubble Space Telescope (HST). The HST observation is a spatially well-resolved image, of spatial resolution 0.1 × 0.1 arcsec, for which there is a segmentation map presented in [START_REF] Rafelski | UVUDF: ultraviolet through near-infrared catalog and photometric redshifts of galaxies in the Hubble Ultra Deep Field[END_REF]. Due to the difference in resolution of the MUSE data, which is only 0.7 × 0.7 arcsec, two distinct sources in the HST image, may overlap in the MUSE data creating a mix in the spectra. In [START_REF] Bacon | The MUSE Hubble Ultra Deep Field Survey-I. Survey description, data reduction, and source detection[END_REF], [START_REF] Inami | The MUSE Hubble Ultra Deep Field Survey-II. Spectroscopic redshifts and comparisons to color selections of high-redshift galaxies[END_REF], the information provided in the Rafelski catalogue is exploited to perform the deblending and prove that MUSE, despite its lower spatial resolution, allows, thanks to the spectral information, to unmix two spatially close or even superimposed sources. We have selected the same portion of the image that is presented in figure 21 of [START_REF] Bacon | The MUSE Hubble Ultra Deep Field Survey-I. Survey description, data reduction, and source detection[END_REF] where the objects identified by For defining the spatial constraint I M ( Ã) required by our model, the HST segmentation map provided by [START_REF] Rafelski | UVUDF: ultraviolet through near-infrared catalog and photometric redshifts of galaxies in the Hubble Ultra Deep Field[END_REF] is used. This map is perfectly registered on the MUSE data. Then by degrading the spatial resolution from 0.1 arcsec to 0.7 arcsec, we obtain binary masks for all the objects present in the field.

For the three central sources, the obtained binary masks are shown on the first line of Fig. 8. The algorithm used for unmixing galaxy spectra corresponds exactly to the version described in section II. One hundred alternated optimisation steps allows to reach a gain eq. ( 12) equal to 10 -10 . It should be noted that the background is considered here as a source, its mask is available in the segmentation map and it is processed in the same way as for the galaxies, to degrade its resolution to the resolution of the MUSE data.

It is impossible to compare quantitatively the results obtained with a ground truth, since no such information exists for the MUSE data. However, same conclusions as in [START_REF] Bacon | The MUSE Hubble Ultra Deep Field Survey-I. Survey description, data reduction, and source detection[END_REF] can be drawn about the spectrum estimated by our algorithm for source #4451: at wavelength λ = 6242.5 Ang, there is an emission line corresponding to object #4451 of Rafelski's catalogue. This emission line has the characteristics of the Lyman-α line (Lyα), namely an asymmetric profile as illustrated in Fig. 7. These result are very similar to the ones presented in Figure 21 in [START_REF] Bacon | The MUSE Hubble Ultra Deep Field Survey-I. Survey description, data reduction, and source detection[END_REF] that is reproduced in Fig. 9 by courtesy of the authors. Similarity between results presented in Fig. 7 and Fig. 9 confirms the interest of our generic approach to solve this particular type of unmixing problem.

C. Discussion

The last application shows the ability of the proposed algorithm to unmix hyperspectral data with a spatially constrained dictionary learning algorithm. It also underlines the necessity of having some external knowledge about the localisation of sources in higher spatial resolution to improve classical dictionary learning unmixing algorithm. In the case of MUSE data, such information coming from the segmentation of another telescope observation (here the Hubble Space Telescope) is rarely available.

Moreover, a good segmentation of this additional image is also required. Besides, this a priori knowledge of sources spatial localisation must be precisely registered to the spectral/temporal data to unmix. As for the fMRI application, in [START_REF] Bhanot | Online dictionary learning for single-subject fMRI data unmixing[END_REF] the authors showed that the estimated values were sensitive to the quality of the localisation maps provided. As these maps are used to define the regularisation term I M (ãi) , their precision influences the estimation of the temporal/spectral signatures and their abundances in strongly mixed voxels: the better the precision of the localisation map, the better the estimation of the spectral/temporal information associated to each source.

From an algorithmic point of view, the strong points of the proposed approach are the small number of parameters to be set and its genericity. The algorithm has a single intrinsic parameter µ σ used in the Tikhonov regularisation, it has only to be set to a very low arbitrary value, as explained in section II-B.

With the scintigraphy application, it has been shown that it is easy to adapt the unmixing problem by adding or removing constraints. In some applications, the proportions of each region have a natural smooth evolution from one pixel to another. A constraint for spatial smoothing within the regions could be introduced in a future variant of the algorithm. However, depending upon the nature of the added constraint, it would increase the complexity of the algorithm. For the current version, the constraints that are in the form of indicator functions, would lead to the introduction of nested loops and consequently to a significant increase in calculation time.

Concerning the execution time of the algorithm, the main factors are the stopping criteria of the different nested iterative algorithms and the size of the images to be unmixed. Furthermore, the calculation time depends much more on the spatial dimension of the image than on the temporal/spectral dimension as the complexity lies in the estimation of A, and to a lesser extent on the number of regions. The algorithm approached convergence for the astronomical data in a few minutes on a laptop, whereas 3 days CPU time was necessary for the fMRI data on a server with 32 cores, with a processing speed of each core around 1603 Mhz. The applications we presented do not require real-time processing but a GPU implementation would lead to a significant gain in speed as the code is highly parallelisable.

V. CONCLUSION

A method to unmix data consisting of an image and a temporal or spectral dimension has been presented. The proposed method proves effective in unmixing problems where some prior information related to the ROIs is available or where such information can be obtained from a registered high resolution image. It should be noted that the method does not take into consideration the morphology, the local structure or the texture of the sources, but only their approximate locations. The spatial constraints are classically expressed as an 1 -penalty to promote sparsity of the mixture in each voxel. Problems with such constraints are generally solved by dictionary learning algorithms. The originality of the proposed approach lies in the replacement of this penalty by a constraint on localisation of the different regions of interests. In some of the examples presented in the article, it has been shown that the algorithm can be easily adapted for specific applications by introducing or dropping constraints of sum-to-one on spatial maps or positivity on timecourses. State-of-the-art results have been obtained on three very different applications. The comprehensiveness of the method makes it possible to easily adapt it to other fields of application such as remote sensing applications for which incorporation of spatial information has been proved to enhance hyperspectral image unmixing performance [START_REF] Wang | A survey of methods incorporating spatial information in image classification and spectral unmixing[END_REF]. The code will be available soon online for the community.

APPENDIX

The convexity of I M (ãi)∩S + (a i ) can be proven if we can prove that the set M (ã i ) ∩ S ∩ R + is convex. This can be easily verified by following the assumption that each voxel contains the contribution of at

where U GT r is the ground truth temporal signal for r th region, U r is the estimated signal for the r th region, N is the length of the temporal signal. NMAE and NMSE can be written in a similar manner to evaluate spatial maps.

Fig. 1 :

 1 Fig. 1: Localisation map and temporal signatures used to build synthetic data.

Fig. 2 :

 2 Fig. 2: Convergence curves plotted for different values of SNR for Dataset I. The curves represent Y-UA F η F plotted against the number of iterations.

FunctionalFig. 3 :

 3 Fig. 3: (a) Scintigraphic data. Ground truth for spatial maps, (b) Initial ROIs , (c) Spatial maps estimated by the proposed algorithm, (d) spatial maps estimated by RUDUR, (e) TACs estimated by RUDUR and our method.

Fig. 4 :

 4 Fig. 4: Rs-fMRI data unmixing. 3D representation of the segmentation map associated with the Allen

FFig. 5 :

 5 Fig. 5: Rs-fMRI data. The figure shows that Y -UA F converges; as is seen by the flattening of the curve in the last iterations.

Fig. 6 :Fig. 7 :

 67 Fig. 6: Rs-fMRI data. The plot shows samples corresponding to the first 500 seconds of the synthetic signals (in blue) and their corresponding estimated timecourses (in dashed red)

Fig. 8 :

 8 Fig. 8: Hyperspectral astronomic data. On top, from left to right, binary mask of sources ID#4451, ID#4460 and ID#4465. Bottom, from left to right, estimated abundance map sources ID#4451, ID#4460 and ID#4465

LyαFig. 9 :

 9 Fig.9: Hyperspectral astronomic data. Reproduction of figure21from paper[START_REF] Bacon | The MUSE Hubble Ultra Deep Field Survey-I. Survey description, data reduction, and source detection[END_REF] with the authors' permission and pending approval from A&A.

  The mean squared errors (MSE) for the estimation of the timecourses are displayed in table I, The MSE in estimating the abundances are given in table III. From table I and table III we see that for most of the regions the errors decrease as SNR increases. The spectral angle distances (SAD)(formula in Appendix B) given in table II follow a similartrend, proving the effectiveness of the unmixing method. The estimated timecourses for each region were normalised by standard deviation before calculating the SAD and MSE. We also generated synthetic data where region 6 was completely included in region 5. In this case, due to noise it was impossible to correctly estimate the timecourse of the region included in the other (and therefore its abundance).

	-20dB -10dB 0dB 10dB 20dB
	reg1 9.5e-02 1.0e-02 1.1e-03 1.1e-04 1.1e-05
	reg2 1.3e-01 1.3e-02 1.2e-03 1.2e-04 1.3e-05
	reg3 2.2e-01 2.1e-02 2.4e-03 3.2e-03 4.0e-03
	reg4 2.3e-01 2.2e-02 2.5e-03 3.2e-03 4.0e-03
	reg5 1.3e-01 1.3e-02 1.3e-03 1.3e-04 1.3e-05
	reg6 1.4e+00 3.4e-01 4.2e-02 4.3e-03 6.3e-04
	reg7 9.6e-03 9.6e-04 9.6e-05 9.6e-06 1.0e-06

TABLE I :

 I Region wise mean squared errors for U for different SNRs

	-20dB -10dB 0dB	10dB 20dB
	reg1 1.8e+01 5.8e+00 1.9e+00 5.9e-01 1.9e-01
	reg2 2.0e+01 6.5e+00 2.0e+00 6.4e-01 2.0e-01
	reg3 2.7e+01 8.3e+00 2.8e+00 3.2e+00 3.6e+00
	reg4 2.7e+01 8.6e+00 2.9e+00 3.2e+00 3.6e+00
	reg5 2.0e+01 6.5e+00 2.1e+00 6.5e-01 2.1e-01
	reg6 7.1e+01 3.4e+01 1.2e+01 3.7e+00 1.4e+00
	reg7 5.6e+00 1.8e+00 5.6e-01 1.8e-01 5.7e-02

TABLE II :

 II Region wise spectral angles for U for different SNRs Fig. 2 illustrates the convergence of the algorithm for different SNRs. Fig. 2 plots Y-UA F

	-20dB -10dB 0dB 10dB 20dB
	reg1 5.8e-04 6.1e-05 4.7e-06 3.2e-07 2.8e-08
	reg2 2.2e-03 2.1e-4 3.8e-05 1.0e-04 1.3e-04
	reg3 3.2e-03 3.4e-04 4.9e-05 5.6e-05 7.0e-05
	reg4 3.4e-03 3.8e-04 5.4e-05 5.9e-05 7.1e-05
	reg5 2.1e-03 2.0e-04 4.6e-05 1.1e-04 1.3e-04
	reg6 9.2e-04 8.9e-05 7.0e-05 2.1e-04 2.6e-04
	reg7 5.4e-03 5.4e-04 6.3e-05 9.4e-05 1.2e-04

η F , as a function of the number of steps in the alternate minimisation and η denotes the white Gaussian noise

TABLE III :

 III Region wise mean squared errors for A for different SNRs present in the dataset at different SNRs. In the very first steps of the alternate optimisation we see that the curves decrease sharply and ultimately settle around a particular value when convergence is achieved.

TABLE IV :

 IV Errors (NMAE and NMSE) and spectral angle differences between the estimated timecourses and the ground truth for the scintigraphy dataset.

TABLE V :

 V Errors between the estimated spatial maps for the different regions and the ground truth. The errors given here were calculated by restricting the pixels of the estimated spatial maps to the initial ROIs.

ACKNOWLEDGMENT

The authors would like to thank IRIS platform [36] at ICube for the Alzheimer's data and equally R.

Bacon for providing MUSE dataset.

2) Algorithm details: Finally, DL is performed at the (low) resolution of the initial fMRI data Y ∈ R N ×P , where N = 490 temporal samples and P = 21024 voxels after extracting the brain. The initial abundance matrix A (0) ∈ R R×P is constructed as follows. Let's say that each voxel i ∈ {1,. . . , P} was subdivided into J high resolution voxels during the artificial augmentation step. For each voxel i of Y and all regions r ∈ {1, . . . , R}, the element (A (0) ) r,i will contain the proportion of high resolution voxels in voxel i, occupied by region r. If region r is not transported to the low-resolution voxel i then (A (0) ) r,i = 0. Matrix à which supports the spatial constraint I M ( Ã) (A) in eq. ( 4) is defined as:

( Ã) r,i = 1 when (A (0) ) r,i > 0 and 0 elsewhere.

3) Results and discussion:

The proposed DL method is applied to the validation data set. Empirically, the algorithm convergences to an acceptable solution for A and U after 500 iterations, see Fig. 5, corresponding to a gain on the optimization κ l < 10 -3 (eq. ( 12)). For the estimation of A, the FISTA algorithm requires a stopping criterion or a maximum number of iterations. In our implementation, FISTA is stopped when a

Fig. 6 shows the estimated timecourses in dashed red against the ground truth signals in blue. The mean squared error (MSE) of the estimated timecourses are given on the plots. Despite the strong mixing in the voxels of the seven regions, our algorithm provides a very good estimate of the synthetic timecourses introduced in the data. In neurosciences applications, these timecourses are then used to build functional brain networks.

B. Hyperspectral astronomic data

In order to confront our algorithm with an unmixing problem on real data, we finally consider the problem of unmixing sources in hyperspectral astronomic data. Some datasets in this domain comprise of hyperspectral datacubes and an external information on the spatial location of the sources. We studied datacubes from the MUSE instrument, the Multi Unit Spectroscopic Explorer [START_REF] Bacon | The MUSE second-generation VLT instrument[END_REF], installed at the Very least one region of interest. By definition:

= +∞ if a i has a non-zero value where ãi is 0 = +∞ if the sum of a i is not equal to 1

= +∞ if at least one element in a i is negative = 0 otherwise where (ã i ) is a binary vector where element (ã i ) r = 1 if the r th region of the segmentation map intersects the i th voxel, and 0 otherwise. A set formed by the intersection

and for all v, w ∈ M (ã i ) ∩ S ∩ R + , a line segment defined :

Proof that z ∈ R + and z ∈ S is straightforward.

Proof that the line segment z belongs to subspace M (ã i ) is detailed in the following.

Let u ∈ R + . By introducing the complementary vector (1 -(ã i )) of ãi , where the notation 1 is a vector of ones of the same size as ãi , we have the equivalence:

Implication u ∈ M (ã i ) =⇒ u T (1 -(ã i )) = 0 is straightforward (vector u lying in M ( ãi ) has its coefficients (u) r = 0 when 1 -(ã i ) r = 1). It is easy to prove that the reverse is true provided that u ∈ R + . As a consequence, vectors v and w lying in M (ã i ) verify:

which yields:

which, since z ∈ R + and according to [START_REF] Varoquaux | Multi-subject dictionary learning to segment an atlas of brain spontaneous activity[END_REF], implies that z ∈ M (ã i ).

The formulae for the different criteria used to analyse the experimental results are given below: