Xiuheng Wang
email: xiuheng.wang@mail.nwpu.edu.cn

Jie Chen
email: dr.jie.chen@nwpu.edu.cn

Cédric Richard
email: cedric.richard@unice.fr

David Brie
email: david.brie@univ-lorraine.fr

LEARNING SPECTRAL-SPATIAL PRIOR VIA 3DDNCNN FOR HYPERSPECTRAL IMAGE DECONVOLUTION

Keywords: Hyperspectral image deconvolution, ADMM, spectral-spatial prior, 3D convolution, deep learning

Hyperspectral image (HSI) deconvolution is an ill-posed problem aiming at recovering sharp images with tens or hundreds of spectral channels from blurred and noisy observations. In order to successfully conduct the deconvolution, proper priors are required to regularize the optimization problem. However, handcrafting a good regularizer may not be trivial and complex regularizers lead to difficulties in solving the optimization problem. In this paper, we use the alternating direction method of multipliers (ADMM) to decompose the optimization problem into iterative subproblems where the prior only appears in a denoising subproblem. Then a 3D denoising convolutional neural network (3DDnCNN) is designed and trained with data for solving this problem. In this way, the hyperspectral image deconvolution is then solved with a framework that integrates the optimization techniques and deep learning. Experimental results demonstrate the superiority of the proposed method with several blurring settings in both quantitative and qualitative comparisons.

INTRODUCTION

Hyperspectral imaging captures a series of images of the same scene over many continuous narrow spectral bands. The high spectral resolution provided by hyperspectral images enables researchers to conduct analyses that cannot be done with conventional imaging techniques. Hyperpectral imaging has been used in many applications such as remote sensing [START_REF] Bioucas-Dias | Hyperspectral remote sensing data analysis and future challenges[END_REF], medical science [START_REF] Lu | Medical hyperspectral imaging: a review[END_REF], atmospheric monitoring [START_REF] Keith | Monitoring effects of a controlled subsurface carbon dioxide release on vegetation using a hyperspectral imager[END_REF]. However, observed hyperspectral images are sometimes blurred during the imaging process, leading to degraded performance in subsequent analyses. Thus, it is desirable to restore images by deconvolution (inversion of the degradation process) techniques before further processing.

Multichannel images contain abundant spectral information across neighboring wavelengths, making image restoration more complex than ordinary 2D images [START_REF] Bongard | 3d deconvolution of hyper-spectral astronomical data[END_REF][START_REF] Sarder | Deconvolution methods for 3-d fluorescence microscopy images[END_REF]. Deconvolution of multichannel (multispectral) images involves Wiener filter [START_REF] Bongard | 3d deconvolution of hyper-spectral astronomical data[END_REF][START_REF] Galatsanos | Digital restoration of multichannel images[END_REF], Kalman filter, and regularized least-squares [START_REF] Galatsanos | Least squares restoration of multichannel images[END_REF]. For hyperspectral deconvolution, in [START_REF] Gaucel | Adaptive-3d-wiener for hyperspectral image restoration: Influence on detection strategy[END_REF] the authors use adaptive 3D Wiener filter, and in [START_REF] Bongard | 3d deconvolution of hyper-spectral astronomical data[END_REF] the authors use filter-based linear methods for astronomic hyperspectral images. 2D Fast Fourier Transforms (FFTs) [START_REF] Thiébaut | Introduction toimage reconstruction and inverse problems[END_REF] and Fourier-wavelet technique [START_REF] Neelamani | Forward: Fourierwavelet regularized deconvolution for ill-conditioned systems[END_REF] are also considered for hyperspectral image restoration, to compute solution efficiently in Fourier and wavelet domains respectively. Since the deconvolution problem is ill-posed, it is important to incorporate prior information on images to regularize the solution. To this end, the fast hyperspectral restoration algorithm in [START_REF] Henrot | Fast positive deconvolution of hyperspectral images[END_REF] performs deconvolution under positivity constraints while accounting for spatial and spectral correlation. In [START_REF] Song | Online deconvolution for industrial hyperspectral imaging systems[END_REF], an online deconvolution algorithm is devised by considering sequentially collected data by push-broom devices.

Similar to conventional image deconvolution, properly defining priors and designing regularizers play an important role in improving the performance and enhancing the stability of the deblurring processing. However, it is a non-trivial task to handcraft a powerful regularizer, and complex regularizers may introduce extra difficulties in solving optimization problems. The variable splitting technique is one of the efficient ways to solve optimization problems with a data fidelity term and a regularization term, especially for the case with non-differentiable regularizers such as 1-norm and TV-norm regularizations. Within this framework, various methods have been proposed based on the ADMM [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] or the half quadratic splitting method [START_REF] Geman | Nonlinear image recovery with halfquadratic regularization[END_REF] to solve a variety of 2D image inverse problems like denoising, deblurring and superresolution [START_REF] Venkatakrishnan | Plugand-play priors for model based reconstruction[END_REF][START_REF] Brifman | Turning a denoiser into a super-resolver using plug and play priors[END_REF][START_REF] Zhang | Learning deep cnn denoiser prior for image restoration[END_REF][START_REF] Chang | One network to solve them all-solving linear inverse problems using deep projection models[END_REF].

Recently, learning the prior from the data is used for 2D image processing and shows its advantageous over using the handcrafted regularizers. Benefiting from the variable splitting technique, this approach plugs the data-based image denoising methods as a module in the optimization iterations to solve various inverse problems [START_REF] Zhang | Learning deep cnn denoiser prior for image restoration[END_REF][START_REF] Chang | One network to solve them all-solving linear inverse problems using deep projection models[END_REF]. This efficient strategy has not been employed in hyperspectral image deconvolution problems, though similar difficulties of designing regularizers are encountered therein. In this work, we propose a hyperspectral deconvolution technique that uses spectralspatial priors learnt from data by a deep neural network. The optimization problem of multichannel image deconvolution is addressed by the ADMM algorithm. A new effective 3DDnCNN is trained as a denoiser to learn priors of hyperspectral images. Experimental results show the effectiveness of the proposed strategy.

PROBLEM FORMULATION

We denote a degraded hyperspectral image and its latent clean counterpart by Y ∈ R P ×Q×N and X ∈ R P ×Q×N respectively, where P , Q, and N are the numbers of rows, columns and spectral dimension of the image. The degraded image and the clean image of the ith spectral band are denoted by Yi ∈ R P ×Q and Xi ∈ R P ×Q . For ease of mathematical formulation, the columns of Yi and Xi Following the linear degradation model in [START_REF] Henrot | Fast positive deconvolution of hyperspectral images[END_REF] and supposing that the convolution is separable over spectral bands, vectors yi and xi are related by

yi = Hixi + ni, for i = 1, • • • , N. ( 1
)
where Hi is the 2D blurring matrix, ni is the additive independent and identically distributed (i.i.d.) Gaussian noise with standard deviation σ for all channels. Under the assumption that {Hi} N i=1 are the same over all channels, i.e., Hi = H for i = 1, • • • , N , we can estimate xi by seeking the minimum of the following objective function:

x = arg min x N i=1 1 2 yi -Hxi 2 + λΦ(x) (2)
where x = col{xi} N i=1 with col{• • • } stack its vector arguments to from a connected vector, the first squared-error term is the data fidelity term, and Φ(x) is the regularizer that enforces desirable properties of the solution with λ ≥ 0 being the regularization parameter. In hyperspectral image processing, spatial and spectral correlations are often encoded in Φ(x).

PROPOSED METHOD AND NETWORK DESIGN

Designing a good regularizor Φ(x) along with efficient solving method is not trivial task. Instead, we propose to learn priors from hyperspectral data and incorporate it into model-based optimization to tackle the regularized inverse problem in [START_REF] Lu | Medical hyperspectral imaging: a review[END_REF]. More specifically, using the variable splitting technique, we transform problem (2) into two sub-problems, namely, a simple quadratic problem and a 3D-image denoising problem. These sub-problems are iteratively solved, using a linear method and a deep neural network, respectively, until the convergence criterion is met.

Variable splitting based on ADMM

ADMM is adopted to decouple the data fidelity term and the regularization term in [START_REF] Lu | Medical hyperspectral imaging: a review[END_REF]. By introducing an auxiliary variable z, problem (2) can be written in the equivalent form:

x = arg min x N i=1 1 2 yi -Hxi 2 + λΦ(z) s.t. z = x with x = col{x1, • • • , xN }. (3)
The associated augmented Lagrangian function is given by where v is the dual variable, and ρ ≥ 0 is the penalty parameter. Scaling v as u = 1 ρ v, (4) can be iteratively solved by repeating the following steps:

Lρ(x, z, v) = arg min x N i=1 1 2 yi -Hxi 2 + λΦ(z) + v T (x -z) + ρ 2 x -z 2 (4)
x k+1 = arg min x N i=1 1 2 yi -Hxi 2 + ρ 2 x -xk 2 (5a)
z k+1 = arg min z λΦ(z) + ρ 2 zk -z 2 (5b)
u k+1 =u k + x k+1 -z k+1 (5c)
where xk = z k -u k and zk = x k+1 + u k . In this way, the data fidelity term and the regularization term in (2) are decoupled into sub-problems (5a) and (5b). Note that (5a) is a least square problem, which can be solved analytically with the solution:

x k+1 = col{(H T H + ρI) -1 (H T yi + ρx k,i )} N i=1 (6)
and (5b) can be reformulated as

z k+1 = arg min z 1 2( λ/ρ) 2 zk -z 2 + Φ(z) (7)
which is actually performed in the 3D image domain. From a Bayesian viewpoint, (7) can be regarded as a denoising problem, removing Gaussian noise with noise-level λ/ρ from the noisy HSI zk to obtain the clean HSI z k+1 . In other words, a denoising operator can be used for the regularization term Φ(x).

Learning spectral-spatial prior via 3DDnCNN

As discussed above, instead of using handcrafted regularizers, we propose to train a CNN-based denoiser for (7) by directly learning the spectral-spatial prior from hyperspectral image datasets. A 3D denoising convolutional neural network is specifically designed for this task. Unlike 2D convolution resulting in spectral information distortion, 3D convolution extracts spatial feature of neighboring pixels and spectral feature of adjacent bands simultaneously without reducing spectral resolution. 3D convolution also involves fewer parameters and is more appropriate for hyperspectral image processing due to the difficulty in capturing a large number of hyperspectral data (see Fig. 2).

K do xk = z k -u k for i = 1 to N do x k+1,i = (H T H + µI) -1 (H T yi + ρx k,i ) end for xk = x k+1 + u k z k+1 = F(z k , Θ) u k+1 = u k + x k+1 -z k+1 k = k + 1 end while
The structure of the proposed 3DDnCNN is illustrated in Fig. 1. Let us define a 3D-block by the composition of a 3D convolution layer (3DConv), a batch normalization (BN) layer and a Rectifier Linear Unit (ReLU) activation function layer. Besides the input and output layers, a 3D convolution layer (3DConv), a ReLU activation function layer, B 3D-blocks and a last 3D convolution layer are sequentially connected to form the proposed network. Batch normalization is also used to speed up the training process as well as to boost the denoising performance [START_REF] Zhang | Beyond a gaussian denoiser: Residual learning of deep cnn for image denoising[END_REF]. The last convolutional layer contains one 3D-filter while the others are composed of 32 3D-filters. The kernel size of each 3D-filter is 3×3×3, which implies that the depth of the kernel in spectral dimension and size of the kernel in spatial dimension are 3 and 3×3 respectively.

The input of the proposed 3DDnCNN is the noisy hyperspectral image z = z + v, where v is the Gaussian noise. Inspired by 2D image denoising algorithm [START_REF] Zhang | Beyond a gaussian denoiser: Residual learning of deep cnn for image denoising[END_REF], we adopt the residual learning to predict the residual error v ≈ v in our denoising network, then we can achieve the estimated clean image by ẑ = zv. For training the network, we use the following loss function:

(Θ) = M m=1 F(zm; Θ) -(zm -zm) 1 (8)
where {(zm, zm)} M m=1 represents M generated noisy-clean HSI (patch) pairs used to train the network function F parameterized by Θ. Note that 1-norm is used as we find it leads to better performance than 2-norm in denoising.

After the 3DDnCNN is trained with massive noisy-clean HSI pairs, it is incorporated into the ADMM-based framework as a denoiser, yielding Algorithm 1.

EXPERIMENTAL RESULTS

In this section, we validate the proposed method to show its effectiveness. We used two public hyperspectral image databases, namely, CAVE [START_REF] Yasuma | Generalized assorted pixel camera: postcapture control of resolution, dynamic range, and spectrum[END_REF] database with 32 indoor HSIs recorded under controlled illuminations, and Harvard [START_REF] Chakrabarti | Statistics of real-world hyperspectral images[END_REF] database with 50 indoor and outdoor HSIs captured under daylight illumination. In the CAVE database, the images consist of 512 × 512 pixels, with 31 spectral bands of 10 nm, covering the visible spectrum 400-700 nm. The images from the Harvard database have a spatial resolution of 1392 × 1040 pixels and 31 channels, ranging from 420 nm to 720 nm at a wavelength interval of 10 nm. The top left 1024 × 1024 pixels were For training and evaluating the proposed 3DDnCNN, the first 20 images of the CAVE database were used as the training set and the others were used as the test set. In the Harvard database, the first 30 hyperspectral images were used as the training set while the others were used as the test set. We implemented our denoising network based on Pytorch and Adam optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] with an initial learning rate 0.0002 and a mini-batch of 64 to minimize the loss function [START_REF] Gaucel | Adaptive-3d-wiener for hyperspectral image restoration: Influence on detection strategy[END_REF] in 500 epochs. The weights were initialized by the method in [START_REF] He | Delving deep into rectifiers: Surpassing human-level performance on imagenet classification[END_REF]. In every epoch of the training phase, each original hyperspectral image was randomly cropped into patches of size 64× 64, and each patch was randomly flipped and rotated for data augmentation.

In our proposed method, the regularization parameter λ was set to 9.6×10 -5 and the penalty factor ρ was set to 0.06 except for setting (c) λ = 1.28×10 -3 and ρ = 0.8 due to its different noise level. The influence of the number of 3D-blocks B is examined in Fig. 4, where the average root-mean-square error (RMSE) value based on the test set in the CAVE database was used to evaluate the denoising performance. As shown in Fig. 4, large values of B generally leads better results. In our deconvolution experiment, we set B = 8 as a larger value by considering the computational cost and memory demand.

Once the denoiser 3DDnCNN was trained, it was plugged into the ADMM framework. The number of iterations K in Algorithm 1 was set to 20 which was sufficient to ensure the convergence. To evaluate the quality of the deblurred images, four quantitative metrics including the root mean-square error (RMSE), peak-signal-tonoise-ratio (PSNR), spectral angle mapper (SAM) [START_REF] Yuhas | Discrimination among semi-arid landscape endmembers using the spectral angle mapper (sam) algorithm[END_REF] and structural similarity (SSIM) [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF]. We compared the proposed method with two deconvolution methods that can be used for hyperspectral images, namely the algorithms with well-designed regularizers in [START_REF] Krishnan | Fast image deconvolution using hyper-laplacian priors[END_REF] and [START_REF] Henrot | Fast positive deconvolution of hyperspectral images[END_REF]. The first method (denoted by HLP), considering the spatial priors, i.e., the hyper-Laplacian priors of images. The second method (denoted by SSP), considering both the spatial and spectral priors of hyperspectral images. The average quantitative results of these methods on the test sets are reported in Table 1. The proposed method achieved satisfactory results in various blur scenarios. For visual comparison, we take scenario (a) for example. Fig. 5 illustrates blurred image, deblurred images, ground truth of real and fake peppers and superballs (two images from CAVE dataset) at the 15h and 24th bands respectively. Fig. 6 shows the spectra difference between recovered data and the clean data (sponges from CAVE dataset) on several pixels. We observe that our proposed method exhibited better spatial and spectral visual results.

CONCLUSION

This paper presented hyperspectral image deconvolution method based on the ADMM framework. Instead of using handcraftted priors, we designed a denoiser based on 3DDnCNN to learn the spectral-spatial information of hyperspectral images from data, and Fig. 6. Spectral difference between restored data and clean data on several pixels. plugged it to ADMM based optimization. Experimental results on two public databases demonstrated that proposed method can effectively handle scenarios with various blurring setting. In the future, we will further improve the structure of 3DDnCNN and extend the proposed method to blind hyperspectral deconvolution.

Fig. 1 .

 1 Fig. 1. Architecture of the proposed 3DDnCNN for hyperspectral image denoising.

Fig. 2 .

 2 Fig. 2. Difference between 2D and 3D convolution.

Algorithm 1

 1 HSI deconvolution with prior learnt from 3DDnCNN. Input: Network parameters Θ, blurred observation y, number of spectral bands N , regularization parameter λ, penalty factor ρ, number of iterations K. Output: Deblurred HSI x. Initialize x = x0, auxiliary variable z0 = x0, scaled dual variable u0 = 0, k = 0. while Stopping criteria are not met and k ≤

Fig. 3 .

 3 Fig. 3. Blurring kernels used in the experiments.

Fig. 4 .

 4 Fig. 4. The denoising performance of 3DDnCNN with different B. extracted in our experiments. The following settings of blur kernels and i.i.d. white Gaussian noise with stand deviation σ were used to generate degraded images in the experiment. (a) 15×15 Gaussian kernel with bandwidth σ k = 1.6, and σ = 0.01. (b) 15×15 Gaussian kernel with bandwidth σ k = 2.4, and σ = 0.01. (c) 15×15 Gaussian kernel with bandwidth σ k = 1.6, and σ = 0.03. (d) Circle kernel with diameter of 7, and σ = 0.01. (e) Motion kernel from [22] of size 15×15, and σ = 0.01. (f) Square kernel with side length of 5, and σ = 0.01. See Fig. 3 for these kernels.For training and evaluating the proposed 3DDnCNN, the first 20 images of the CAVE database were used as the training set and the others were used as the test set. In the Harvard database, the first 30 hyperspectral images were used as the training set while the others were used as the test set. We implemented our denoising network based on Pytorch and Adam optimizer[START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] with an initial learning rate 0.0002 and a mini-batch of 64 to minimize the loss function[START_REF] Gaucel | Adaptive-3d-wiener for hyperspectral image restoration: Influence on detection strategy[END_REF] in 500 epochs. The weights were initialized by the method in[START_REF] He | Delving deep into rectifiers: Surpassing human-level performance on imagenet classification[END_REF]. In every epoch of the training phase, each original hyperspectral image was randomly cropped into patches of size 64× 64, and each patch was randomly flipped and rotated for data augmentation.In our proposed method, the regularization parameter λ was set to 9.6×10 -5 and the penalty factor ρ was set to 0.06 except for setting (c) λ = 1.28×10 -3 and ρ = 0.8 due to its different noise level. The influence of the number of 3D-blocks B is examined in Fig.4, where the average root-mean-square error (RMSE) value based on the test set in the CAVE database was used to evaluate the denoising performance. As shown in Fig.4, large values of B generally leads better results. In our deconvolution experiment, we set B = 8 as a larger value by considering the computational cost and memory demand.Once the denoiser 3DDnCNN was trained, it was plugged into the ADMM framework. The number of iterations K in Algorithm 1 was set to 20 which was sufficient to ensure the convergence. To evaluate the quality of the deblurred images, four quantitative metrics including the root mean-square error (RMSE), peak-signal-tonoise-ratio (PSNR), spectral angle mapper (SAM)[START_REF] Yuhas | Discrimination among semi-arid landscape endmembers using the spectral angle mapper (sam) algorithm[END_REF] and structural similarity (SSIM)[START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF]. We compared the proposed method

Fig. 5 .

 5 Fig. 5. Visual results for comparison of the methods.

Table 1 .

 1 Average RMSE, PSNR, SAM, SSIM of different methods on CAVE and Harvard database in 6 blur scenarios.

	Scenarios Methods	CAVE database RMSE PSNR SAM	SSIM	Harvard database RMSE PSNR SAM	SSIM
		HLP	3.781	37.427 11.19 0.9257	3.562	38.181	8.41	0.9133
	(a)	SSP	4.299	36.536	8.24	0.9474	3.687	38.657	5.06	0.9353
		proposed	2.559	41.173	6.27	0.9602	2.722	41.116	4.84	0.9486
		HLP	4.947	35.340 10.98 0.9058	4.356	36.885	8.07	0.8903
	(b)	SSP	5.285	34.697	8.49	0.9255	4.431	37.225	5.10	0.9110
		proposed	4.774	35.961	7.77	0.9199	3.766	38.852	4.77	0.9202
		HLP	7.778	30.042 25.83 0.6140	7.729	30.515 23.24 0.5947
	(c)	SSP	4.606	35.655 12.23 0.9161	4.060	37.079	8.22	0.9046
		proposed	3.898	37.424	6.46	0.9395	3.737	38.958	4.58	0.9165
		HLP	4.486	36.102 11.52 0.9116	4.004	37.354	8.58	0.8992
	(d)	SSP	4.788	35.609	8.41	0.9350	3.992	38.085	5.08	0.9232
		proposed	3.017	39.887	6.80	0.9477	2.987	40.500	4.82	0.9384
		HLP	4.312	36.183 12.97 0.9035	3.891	37.248	9.88	0.8945
	(e)	SSP	4.473	36.151	8.51	0.9427	3.807	38.372	5.19	0.9317
		proposed	2.024	43.042	6.30	0.9638	2.306	42.258	4.77	0.9568
		HLP	3.933	37.054 11.77 0.9211	3.642	37.928	8.86	0.9090
	(f)	SSP	4.305	36.539	8.33	0.9463	3.646	38.771	5.08	0.9347
		proposed	2.078	42.718	6.29	0.9669	2.458	41.807	4.87	0.9541

This work was funded in part by NSFC 61811530283, 111 project B18041, the seed Foundation of Innovation and Creation for Graduate Students in Northwestern Polytechnical University. The work of C. Richard was funded in part by the CoopIntEER program CNRS-NSFC (DIALOG project).


Aneline Dolet
email: aneline.dolet@grenoble-inp.fr

Ri Ammanouil
email: ritaammanouil@gmail.com

Virginie Petrilli
email: virginie.petrilli@lyon.unicancer.fr

Cédric Richard
email: cedric.richard@unice.fr

Piero Tortoli
email: piero.tortoli@unifi.it

Didier Vray
email: didier.vray@creatis.insa-lyon.fr

François Varray
email: francois.varray@creatis.insa-lyon.fr

In Vitro and In Vivo Multispectral Photoacoustic Imaging for the Evaluation of Chromophore Concentration

Keywords: photoacoustic imaging, spectral unmixing, blood oxygen concentration evaluation

Multispectral photoacoustic imaging is a powerful noninvasive medical imaging technique that provides access to functional information. In this study, a set of methods is proposed and validated, with experimental multispectral photoacoustic images used to estimate the concentration of chromophores. The unmixing techniques used in this paper consist of two steps: (1) automatic extraction of the reference spectrum of each pure chromophore; and (2) abundance calculation of each pure chromophore from the estimated reference spectra. The compared strategies bring positivity and sum-to-one constraints, from the hyperspectral remote sensing field to multispectral photoacoustic, to evaluate chromophore concentration. Particularly, the study extracts the endmembers and compares the algorithms from the hyperspectral remote sensing domain and a dedicated algorithm for segmentation of multispectral photoacoustic data to this end. First, these strategies are tested with dilution and mixing of chromophores on colored 4% agar phantom data. Then, some preliminary in vivo experiments are performed. These consist of estimations of the oxygen saturation rate (sO 2 ) in mouse tumors. This article proposes then a proof-of-concept of the interest to bring hyperspectral remote sensing algorithms to multispectral photoacoustic imaging for the estimation of chromophore concentration.

Introduction

Photoacoustic imaging is a hybrid medical imaging technique that provides access to functional information. Nanosecond laser pulse illumination can provide optical properties of biological tissues that are related to molecular composition. The transmitted optical energy is absorbed by the optical absorbers in the tissues, which then undergo thermal expansion. At each laser pulse, this phenomenon generates ultrasonic waves that can be detected at the tissue surface by a conventional ultrasound system. The use of several laser wavelengths provides the multispectral absorption characteristics that can be used to distinguish the imaged tissues from each other [START_REF] Beard | Biomedical photoacoustic imaging[END_REF][START_REF] Van Veen | Determination of VIS-NIR absorption coefficients of mammalian fat, with time-and spatially resolved diffuse reflectance and transmission spectroscopy[END_REF]. In addition, as the absorption by a chromophore is linearly related to its concentration, multispectral photoacoustic imaging can also be used to determine the concentrations of chromophores present in the imaged region [START_REF] Glatz | Blind source unmixing in multi-spectral optoacoustic tomography[END_REF][START_REF] Razansky | Multispectral photoacoustic imaging of fluorochromes in small animals[END_REF]. These two properties, discrimination of tissues and concentration evaluation, are of major interest, depending on the application at hand. They can be exploited to determine different concentrations of the same chromophore (e.g., estimation of the concentration of a contrast agent in the body [START_REF] Mienkina | Evaluation of Ferucarbotran (Resovist) as a Photoacoustic Contrast Agent[END_REF]), or to distinguish one particular chromophore from all of the other imaged ones without considering its dilution (e.g., determination of the level of vascularization and calculation of the concentration of oxygenated and deoxygenated blood, for oxygenation rate evaluation [START_REF] Mercep | Combined pulse-echo ultrasound and multispectral optoacoustic tomography with a multi-segment detector array[END_REF][START_REF] Singh | In vivo demonstration of real-time oxygen saturation imaging using a portable and affordable LED-based multispectral photoacoustic and ultrasound imaging system[END_REF], cancer tissue evaluation [START_REF] Jnawali | Automatic cancer tissue detection using multispectral photoacoustic imaging[END_REF], or imaging connectivity in brain [START_REF] Nasiriavanaki | High-resolution photoacoustic tomography of resting-state functional connectivity in the mouse brain[END_REF]).

In photoacoustic images, each pixel can be interpreted as a mixture of pure chromophores, where the weights represent their average concentrations. The spectra of pure chromophores are called endmembers, and their weights are the abundance coefficients. Unmixing techniques for multispectral photoacoustic data are presented in this article, with the objective being to estimate the concentrations of the chromophores present in the imaged area. In the field of photoacoustic imaging, several unmixing methods have already been used to detect chromophores from the body [START_REF] Tzoumas | Unmixing Molecular Agents From Absorbing Tissue in Multispectral Optoacoustic Tomography[END_REF] or contrast agents injected into tissues. For instance, principal component analysis (PCA) and independent component analysis were used in Reference [START_REF] Glatz | Blind source unmixing in multi-spectral optoacoustic tomography[END_REF]. However, the abundance matrices calculated by these two methods do not allow accurate estimation of concentrations of chromophores. In particular, these matrices can contain negative values that cannot be interpreted as concentrations even after scaling. This issue is one of the challenges addressed in this paper. Indeed, this study highlights experimental results as a proof of concept of the usefulness of photoacoustic imaging to evaluate chromophore concentration.

The concentration of an endmember should be equal to zero if the chromophore is absent in the considered pixel, and equal to one if it is a pure pixel. The endmember will lie in the interval (0, 1) if the pixel is partly composed of the chromophore (i.e., diluted chromophore, or a mix of chromophores). The abundance matrix must then be subject to the following constraints: (1) each abundance coefficient must be in the range [0, 1]; and

(2) for a given pixel, the sum of all K abundance coefficients, which correspond to the K endmembers, must be less than 1 (i.e., diluted chromophore) or equal to 1 (i.e., pure chromophore, or a mixture of chromophores). These constraints have never been investigated in the field of photoacoustics, while they are usually encountered in remote sensing, where imaged areas might contain, for example, roads, buildings, rivers, or forests, which correspond to the chromophores to be unmixed. The abundances are, thus, related to the local concentrations of the chromophores in the scene. Spectral unmixing is also a widely explored field for various applications, like food safety [START_REF] Gowen | Hyperspectral imaging-An emerging process analytical tool for food quality and safety control[END_REF][START_REF] Bermana | Classification of Sound and Stained Wheat Grains Using Visible and near Infrared Hyperspectral Image Analysis[END_REF], pharmaceutical process monitoring [START_REF] Gendrin | Pharmaceutical applications of vibrational chemical imaging and chemometrics: A review[END_REF], and industrial and forensic applications [START_REF] Picon | Fuzzy Spectral and Spatial Feature Integration for Classification of Nonferrous Materials in Hyperspectral Data[END_REF][START_REF] Brewer | Forensic analysis of bioagents by X-ray and TOF-SIMS hyperspectral imaging[END_REF]. However, unmixing of remote sensing data is the most closely related area to our multispectral photoacoustic data. As many efficient algorithms have been devised and analyzed over the last decade, we focus on this particular field in this study. The main challenge of this paper is then to show that photoacoustic imaging can benefit from these innovations, principally considering the constraints presented above, because of the similarities between these two fields. However, photoacoustic imaging contains less spectral bands than remote sensing and is subject to spectrally dependent light attenuation which makes the unmixing even more challenging. Moreover, the methods presented in this study are based on relative comparison between spectra. Even if PCA is also a relative method, as mentioned before, it does not allow proper evaluation of concentration, based on authors preliminary experiments. The study explores, then, other field to adapt unmixing methods that better suit the considered goal. As all spectra are impacted in the same way by the light absorption and ultrasound attenuation for each wavelength, and these methods allow us to avoid taking into account these parameters.

Unmixing of Photoacoustic Data

A multispectral photoacoustic image consists of a two-dimensional area that is imaged at L wavelengths. Each pixel x i is characterized by its coordinates s i ∈ IN 2 and is endowed with a spectrum A i ∈ IR L . The N samples, x i , of the region of interest are expressed as follows:

x i = [s i , A i ] ∈ X, (1)
with i ∈ {1, . . . , N} as the sample index; see Figure 1a. Before any other processing, the data need to be normalized over the range [0,255]. This allows the same range of values for data from different acquisition systems. Two types of pixels are illustrated in Figure 1: pixel x 1 contains only noise, called the background, with low and flat photoacoustic amplitudes at all wavelengths (Figure 1b, green curve A 1 ), and pixel x 2 provides information on the optical absorbers, with a high photoacoustic signal (Figure 1b, purple curve A 2 ).

Linear Mixing Model

The acquisition of photoacoustic data over a region of interest at several wavelengths provides a multispectral (i.e., three-dimensional) image, where each pixel is endowed with a spectrum. This spectrum is either pure, and is then considered as a single endmember, or mixed if it is composed of a mixture of endmembers. The linear mixing model (LMM) [START_REF] Keshava | Spectral unmixing[END_REF] considers that each mixed pixel is a convex combination of the spectra of the endmembers. More formally, this is defined as follows:

A i = K ∑ k=1 u ki E k + g i , ∀ i ∈ {1, . . . , N}, (2)
where A i ∈ IR L is the spectrum at the i-th pixel, L is the number of wavelengths, K is the number of endmembers, u ki is the abundance of the k-th endmember in the i-th pixel, E k is the L-dimensional spectrum of the k-th endmember, and g i is a vector of zero-mean white noise that defines the sensor noise and the error of the model. All of the vectors are column vectors. Then, U is used for the matrix of abundances with the generic (i, j)-th entry u ij . As abundances are relative contributions, they must be positive and their sum has to be equal to one:

   u ki ≥ 0, ∀ i ∈ {1, . . . , N}, ∀k ∈ {1, . . . , K} K ∑ k=1 u ki = 1, ∀ i ∈ {1, . . . , N}. (3)
The LMM is a simple but effective model that is widely used in remote sensing [START_REF] Keshava | Spectral unmixing[END_REF]. For medical applications, the LMM can be a powerful tool for estimation of concentration of chromophores. In this context, it is assumed that the observed spectra can either correspond to a fully concentrated chromophore or a diluted chromophore, or to a mixture of several chromophores. The original LMM with the constraints of Equation (3) does not consider dilution. It can, however, be extended to this scenario by relaxing the sum-to-one constraint, which can be performed by considering a zero endmember with the LMM, called a shadow endmember.

Unmixing Strategy

Data unmixing requires inversion of a mixing model that can be either linear or nonlinear, depending on the hypotheses. The study focus on unsupervised methods to be robust to the various photoacoustic systems that exist. Unsupervised algorithms have been devised in remote sensing to extract endmembers and to estimate the abundance matrix. Unsupervised methods, such as group lasso with unit sum and positivity constraints (GLUP) [START_REF] Ammanouil | Blind and fully constrained unmixing of hyperspectral images[END_REF] and vertex component analysis (VCA) [START_REF] Nascimento | Vertex Component Analysis: A Fast Algorithm to Unmix Hyperspectral Data[END_REF], can do both of these simultaneously, subject to constraints [START_REF] Glatz | Blind source unmixing in multi-spectral optoacoustic tomography[END_REF]. Other unsupervised methods, such as N-FINDR [START_REF] Winter | An Algorithm for Fast Autonomous Spectral End-Member Determination in Hyperspectral Data[END_REF], only extract the endmembers. A supervised unmixing algorithm is then required to calculate the abundances. The fully constrained least-square algorithm (FCLS) [START_REF] Heinz | Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF] usually provides satisfactory performance as long as the endmembers are accurately extracted. FCLS can be used with GLUP, VCA, and N-FINDR [START_REF] Ammanouil | Blind and fully constrained unmixing of hyperspectral images[END_REF], or with any other strategy for extraction of the endmembers. In this study, we use FCLS to estimate the abundance matrices.

In the following, we use GLUP, VCA (that was already used in photoacoustic imaging in Reference [START_REF] Deán-Ben | Fast unmixing of multispectral optoacoustic data with vertex component analysis[END_REF]), and N-FINDR to estimate the endmembers in experimental multispectral photoacoustic data to evaluate the interest of remote sensing methods to process these data. We compare their performances with the spatio-spectral mean-shift (SSM-S) algorithm [START_REF] Dolet | Spatial and Spectral Regularization to Discriminate Tissues Using Multispectral Photoacoustic Imaging[END_REF] devised by some of the present authors.

Method

Pre-Processing

The pre-processing introduced in this section aims to discriminate pixels of interest from the background. Among the possible strategies, background pixels can be determined by thresholding filtered data with a Sobel filter. The Sobel threshold λ S can be set as follows:

λ S = 2 * ∇ S , (4)
where ∇ S is the mean of the Sobel gradient magnitude [START_REF] Pratt | Digital Image Processing[END_REF]. For our application, only the edges detected by the Sobel filter were used to calculate ∇ S . Threshold λ S was applied to the sum over all of the wavelengths of the multispectral photoacoustic image, to design a mask where the background pixels were set to 0. The estimation methods introduced hereafter were only applied to pixels out of the background.

Endmember Extraction

GLUP Algorithm

Group lasso with unit sum and positivity constraints assumes that the endmembers are not known but are present in the form of some pure (unmixed) pixels in the image [START_REF] Ammanouil | Blind and fully constrained unmixing of hyperspectral images[END_REF]. Based on this assumption, the LMM (2) can be reformulated as follows:

A i = N ∑ j=1 u G ji A j + g i , ∀ i ∈ {1, . . . , N}, (5)
where u G ji is the abundance of A j in A i . On the one hand, if A j is an endmember, the j-th row U G j of matrix U G calculated with GLUP should have nonzero entries that represent the abundance map of A i . On the other hand, if A j is a mixed pixel, U G j should have all of its entries equal to zero. This means that U G should have N -K rows of zero, with the other K rows equal to the rows of U.

The premise in GLUP is that U G allows the identification of the endmembers in observation A through its nonzero rows. The resulting unmixing problem requires that U G only has a few rows different from zero, in addition to the non-negativity and sum-to-one constraints. This leads to the following convex optimization problem:

min U G 1 2 N ∑ j=1 A -AU G 2 F + µ N ∑ j=1 U G j 2 (6)
subject to:

     u G ji ≥ 0, ∀ i, j ∈ {1, . . . , N} N ∑ j=1 u G ji = 1, ∀ i ∈ {1, . . . , N},
where µ > 0 is a small regularization parameter that is set by the user, and

A = [A 1 , . . . , A N ]
is the matrix of all pixels to unmix. The first term in Eqtuation [START_REF] Mercep | Combined pulse-echo ultrasound and multispectral optoacoustic tomography with a multi-segment detector array[END_REF] ensures that the observations match the model in Equation ( 5), and the second term is the group lasso regularizer that induces sparsity by possibly driving several rows of U G to zero [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF]. The minimization problem has constraints to ensure that the abundances obey the positivity and the sum-to-one constraints. This optimization problem can be solved with a primal-dual method; see Reference [START_REF] Ammanouil | Blind and fully constrained unmixing of hyperspectral images[END_REF] for details.

To conclude, GLUP allows the identification of the endmembers in A by identification of the nonzero rows in U G . GLUP also provides the estimated abundances that correspond to the nonzero rows in the estimated matrix U G . In this paper, similarly to Reference [START_REF] Ammanouil | Blind and fully constrained unmixing of hyperspectral images[END_REF], we did not use the abundances estimated by GLUP. For a fair comparison with other algorithms, we re-estimated the abundances with the FCLS algorithm from the endmember spectra estimated using GLUP.

VCA Algorithm

Vertex component analysis also assumes the presence of pure pixels in the data [START_REF] Nascimento | Vertex Component Analysis: A Fast Algorithm to Unmix Hyperspectral Data[END_REF]. The driving principle of the VCA algorithm is to project the data onto a direction orthogonal to the subspace spanned by the endmembers already extracted. The new extracted endmember is the farthest signal in this projection domain. Considering this endmember, a new subspace is calculated, and the same procedure is iteratively performed until the preset number K of endmembers to extract is reached.

The first step of the VCA algorithm consists of determination of the initial subspace. Two methods are recommended, which depend on the signal-to-noise ratio. If the signal-to-noise ratio is greater than λ V in Equation [START_REF] Singh | In vivo demonstration of real-time oxygen saturation imaging using a portable and affordable LED-based multispectral photoacoustic and ultrasound imaging system[END_REF], this first subspace is calculated using the singular value decomposition algorithm [START_REF] Scharf | Statistical Signal Processing: Detection, Estimation, and Time Series Analysis[END_REF]. Otherwise, the subspace considered is constructed using principle component analysis. The threshold λ V is defined in Reference [START_REF] Nascimento | Vertex Component Analysis: A Fast Algorithm to Unmix Hyperspectral Data[END_REF] as follows:

λ V = 15 + 10 log 10 (K). (7)
Let us use S k to denote the subspace available at iteration k. A vector v k orthonormal to S k is calculated as follows:

v k = r k -S k S + k r k r k -S k S + k r k , ( 8
)
where r z is a zero-mean random vector, and S + z is the pseudo inverse of S k . The observed spectra A are then projected onto direction v k :

f k = v k A. (9)
The largest entry of f k in absolute values allows the designation of the spectrum in A to be considered as a new endmember. This endmember is then added to the set of endmembers that have already been extracted, to define the subspace S k+1 that is considered at the next iteration. The procedure is stopped when the number K of desired endmembers is reached. They are stored in matrix E V , and the abundance matrix U V is calculated by projecting A onto E V . For a fair comparison with other algorithms, we estimated the abundances with the FCLS algorithm from the endmember spectra E V estimated with VCA.

N-FINDR algorithm

N-FINDR also assumes that a pure pixel for each chromophore to unmix is present in the dataset [START_REF] Winter | An Algorithm for Fast Autonomous Spectral End-Member Determination in Hyperspectral Data[END_REF][START_REF] Plaza | An Improved N-FINDR Algorithm in Implementation[END_REF]. The first step is the random generation of a set of K endmembers to produce the matrix E 0 of endmembers that are available at iteration 0. At each iteration k, the following volume is calculated:

V k = | det(E k )| (K -1)! , ( 10
)
with

E k = 1 . . . 1 E (1) k . . . E (K) k
, where

E (p)
k is the p-th endmember in the matrix of endmembers E k available at the k-th iteration. All of the pixels A i in A take place successively in E k , where they are successively substituted to all possible columns p. Volume V k+1 is updated as follows, for instance:

V k+1 = det 1 . . . 1 . . . 1 E (1) k . . . A i . . . E (K) k (K -1)! . ( 11
)
If the new volume is greater than the previous one, V k and E k are updated with the new values. The procedure is stopped when the N pixels of the dataset have been tested. The resulting endmembers are stored in matrix E N . No abundance coefficients are calculated with this method, and the FCLS can be used to this end.

SSM-S Algorithm

Spatio-spectral mean-shift is a clustering method that was introduced by Reference [22] for chromophore discrimination in multispectral photoacoustic images. It is based on a spatio-spectral regularization, to cluster the pixels that are spatially and spectrally close.

Consider a pixel x i . Its neighbors in spatial dimensions at a radial distance less than R S are first considered. These pixels x j = s j , A j have to satisfy the following:

1 R S s i -s j 2 ≤ 1. ( 12
)
Assuming that spectra of the same chromophore are close, we also consider that two spectra can be from the same chromophore if their distance is less than

R λ ; namely, 1 R λ A i -A j ∞ ≤ 1. ( 13
)
Note that we found it interesting in Equation ( 13) to consider the infinite norm, to discriminate more effectively spectra that might differ only in narrow frequency bands.

The SSM-S only applies to pixels that satisfy both constraints of Equations ( 12) and ( 13). It consists of updating their locations and spectra as follows:

x [t+1] i = N ∑ j=1 g S (s [t] i ; s [t] j ) • g λ (A [t] i ; A [t] j ) • x [t] j N ∑ j=1 g S (s [t] i ; s [t] j ) • g λ (A [t] i ; A [t] j ) , ( 14
)
where g S (s 12) and ( 13) hold, respectively, and 0 otherwise. This iterative procedure is applied to all of the pixels x i until convergence; see Reference [START_REF] Dolet | Spatial and Spectral Regularization to Discriminate Tissues Using Multispectral Photoacoustic Imaging[END_REF] for details. Unlike the GLUP, VCA, and N-FINDR algorithms, note that SSM-S uses spatial information in addition to spectral information. SSM-S is studied here because of this particular capacity which is not often taken into account to process multispectral photoacoustic data. Updating of Equation ( 14) leads to clusters of pixels, with each defined by a centroid. Each centroid can be considered as an endmember and stored in a matrix, which is denoted by E S . As no abundance coefficients are calculated with this method, the FCLS can be used to this end.

[t] i ; s [t] j ) = 1 and g λ (A [t] i ; A [t] j ) = 1, if Equations (
Before concluding this section, we want to draw attention to the following: that, unlike the N-FINDR and SSM-S algorithms, GLUP and VCA converge to stationary points that can vary significantly depending on their initialization.

Abundances Estimation

The FCLS algorithm is widely used in remote sensing to estimate abundances, based on endmember knowledge. It considers the constraints of Equation ( 3) when solving the following optimization problem:

min U 1 2 N ∑ j=1 A -UE 2 F , (15)
with

     u ji ≥ 0, ∀ i, j N ∑ j=1 u ji = 1, ∀ i .

Materials

Acquisition System

A commercial multispectral photoacoustic system (Vevo LAZR; Visualsonics, Fujifilm) was used for this study; see Figure 2. The optical source of this system is a Nd:YAG pulsed laser with a pulse duration of 5 ns and a 20 Hz repetition rate, which is coupled to an optical parametric oscillator to select the transmission wavelengths [START_REF] Arthuis | Real-Time Monitoring of Placental Oxygenation during Maternal Hypoxia and Hyperoxygenation Using Photoacoustic Imaging[END_REF]. The wavelength can be set from 680 nm to 970 nm. The photoacoustic probe that was used (LZ400) is composed of 256 elements, which can acquire ultrasound signals in the frequency range from 18 MHz to 38 MHz. The imaging depth is approximately 1.5 cm and the imaged region of interest is around 1 cm. The images are reconstructed with a delay-and-sum algorithm. For every acquisition, the corresponding fluence is saved. For the acquisitions presented in this study, the fluence is varying of only 2% between the minimal and maximal fluences. As the variation is low and the proposed methods are based on unsupervised endmembers extraction with relative comparison between spectra, the images were not normalized by the fluence. The data were acquired using the full range of available wavelengths (680 nm, 970 nm) with 1 nm steps, which required less than 1 min of acquisition time. The time to switch from one wavelength to the other and to acquire the corresponding image is around 0.2 s. Usually, in multispectral photoacoustic imaging, from 5 to 10 wavelengths are used because the information collected is then sufficient to discriminate or quantify biological chromophores with an acceptable processing time. To get closer to this framework, only eight of the acquired wavelengths were selected to construct our dataset, from 680 nm to 820 nm, with 20 nm steps. These wavelengths were chosen because the used optical absorbers (blue and green inks) can be discriminated in this range (see Figure 3).

The algorithms were computed with MATLAB (The MathWorks Inc., Natick, MA, USA) on a computer with a 16-GB memory and an Intel Core i7-5500U CPU at 2.40 GHz.

Figure 3. Spectra, measured with a Perkin Elmer Lambda900 spectrometer for the study [START_REF] Dolet | Spectrophotometry and Photoacoustic Imaging: A Comparative Study[END_REF], of the blue and green inks used to make the phantoms (a.u., arbitrary units).

Imaged Phantoms

The homemade colored phantoms used to construct the dataset were composed of three different 4% colored agar parts. The ones on the left side in Figure 4a were prepared using only blue ink (400 g water, 16 g agar, 380 µL ink), and the ones on the right with only green ink (400 g water, 16 g agar, 950 µL ink). Both of these are considered as relative ink concentrations equal to 1 (i.e., pure chromophores). The blue and green ink quantities were not equal because both of these inks do not show the same maximum photoacoustic signal amplitude. Based on a previous study [START_REF] Dolet | Spectrophotometry and Photoacoustic Imaging: A Comparative Study[END_REF], the ink quantities were chosen to obtain the same maximum photoacoustic signal amplitude, on the range from 400 nm to 1200 nm, for both of these pure chromophores. The used blue and green inks spectra, measured with a Perkin Elmer Lambda900 spectrometer for the study [START_REF] Dolet | Spectrophotometry and Photoacoustic Imaging: A Comparative Study[END_REF], are shown in Figure 3, on the wavelength range of interest for the present study.

In Figure 4b, the sample of shown spectra seems to highlight that the fully concentrated blue and green parts do not exhibit the same, even approximately, spectral intensities for both phantoms. However, this is only due to the spectra shown here. Nevertheless, the Figure shows that the shape of blue or green spectra are quite similar which is the property used by the proposed methods in this study.

Chromophore Dilution

The central part of the left phantom in Figure 4a is composed of a 0.53 dilution of the blue ink relative to the pure concentration indicated above. Therefore, here 100 g water, 4 g agar, and 50 µL blue ink were used. This phantom is referred to as B-Bdil-G(Vevo) in the following.

Mixing of the Chromophores

The central part of the second phantom, shown in Figure 4a-right, is a mix of the blue and green inks (50 g water, 2 g agar, 20 µL blue ink, 80 µL green ink). This means that the mix corresponds to a blue concentration of 0.42 and a green concentration of 0.67, relative to the pure concentrations indicated above. This phantom will be referred to as B-mix-G(Vevo) in the following.

Performance Evaluation on Phantoms

The performances of the GLUP/FCLS, VCA/FCLS, N-FINDR/FCLS, and SSM-S/FCLS strategies were evaluated on the phantoms described in Section 4.2. This means that FCLS was used to estimate the concentrations of the endmembers extracted beforehand with the other methods presented in Section 3.2. The resulting abundance maps were compared to the ground truth.

For each abundance map of each endmember, a mean concentration was calculated for each part of the imaged region, as illustrated in Figure 5. The illustrative phantom considered in Figure 5 is composed of a pure blue chromophore, a dilution of this chromophore, and a pure green chromophore (from left to right). Considering the abundances related to the blue endmember, three average values were calculated: one for each part, as circled by the yellow dotted line in Figure 5. As we assume that the pure blue chromophore, which corresponds to the endmember of interest, is present in the image, the mean values were normalized by their maximum, as circled by the orange dashed line in Figure 5. Where the data would be perfectly unmixed, the maximum value should correspond to the pure blue part of the phantom (i.e., the left part), which is then set to 1 after normalization. The other normalized values should correspond to the relative concentrations of the pure chromophore, i.e., the relative concentration of blue in the illustrative example in Figure 5.

This calculation was performed to estimate the abundance maps provided by the different unmixing strategies. Indeed, by considering the light absorption and the diffusion in tissue, each pixel of similar region exhibit similar spectra shape tendency (see Figure 4b).

For each pixel, a particular concentration is then evaluated. The normalized mean values were compared with the ground truth to evaluate the performances of the algorithms on each whole part.

In Vivo Data Acquisitions

Images were acquired with the Vevo LAZR system. For imaging the mice tumors, the BALBc MMTV-Neu transgenic mouse model was used, which expresses the Her-2/neu V664E oncogene [START_REF] Boggio | Interleukin 12-mediated prevention of spontaneous mammary adenocarcinomas in two lines of Her-2/neu transgenic mice[END_REF]. This model develops spontaneous mammary tumors that are visible with the photoacoustic system after 17 weeks (between 5 mm 3 and 10 mm 3 in size). Palpable tumors are detected around 20 weeks old; see Figure 6, where the white arrows indicate some of the tumors, and the green square denote the imaged one. A 20-week-old mouse was used for the imaging, where the region of interest was shaved the day before using a commercial hair remover cream (Veet, Cream hair remover; Reckitt Benckiser, UK) to avoid any imaging noise coming from the mouse hair. The mouse was anesthetized with a mix of 3.5% isoflurane and 96.5% oxygen before the photoacoustic acquisitions, using an anesthesia mask. During the imaging procedure, the mouse was maintained asleep through inhalation of a mix that contained 1.5% isoflurane (Figure 6a), to minimize movements during the acquisitions at the different wavelengths. Heartbeat and breathing were monitored by several devices connected to the mouse paws (Figure 6b, red arrows) and the anesthesia mask, respectively. This allowed monitoring of the health of the mouse and the acquisition of images at the same position. In this way, heartbeat and breathing motions were avoided as much as possible. The mice were maintained in the specific pathogen free animal facility at the Cancer Research Center of Lyon, at Center Léon Bérard. All of the mice experiments were performed in accordance with the animal care guidelines of the European Union and were validated by the local Animal Ethics Evaluation Committee and the French Ministry of Higher Education and Research (C2EA-15 and 02296.02).

The dataset used to test this strategy was composed of 15 images that were acquired from 680 nm to 960 nm, in 20 nm steps. The results were compared to the Vevo LAZR sO 2 map. To this end, acquisitions on the same tumor were also performed using the Vevo LAZR Oxy-Hemo mode.

sO 2 Calculation with Vevo LAZR Oxy-Hemo Mode

The Oxy-Hemo mode acquires two photoacoustic images, at 750 nm and 850 nm. With both of these images, the concentrations of deoxygenated hemoglobin (Hb) and oxygenated hemoglobin (HbO 2 ) were measured, which are expressed as [Hb] and [HbO 2 ], respectively [START_REF] Visualsonics | Oxy-Hemo Mode[END_REF]. The oxygen saturation was then calculated as follows:

sO 2 = [HbO 2 ] [Hb] + [HbO 2 ] . ( 16
)
This ratio is classically used in biomedical applications to provide oxygenation information [START_REF] Li | Photoacoustic Tomography and Sensing in Biomedicine[END_REF][START_REF] Deán-Ben | Functional Optoacoustic Imaging of Moving Objects Using Microsecond-Delay Acquisition of Multispectral Three-Dimensional Tomographic Data[END_REF]. Both of the concentration maps were saved, as well as the Vevo LAZR sO 2 displayed map, for further comparisons. The sO 2 acquisition is conducted just after the acquisitions on the whole range of wavelength, allowing then a fair comparison of the estimated concentrations.

Results

As indicated in Section 3.2, the N-FINDR and SSM-S algorithms converge to a single set of estimated endmembers, regardless of their initialization. On the contrary, the GLUP and VCA methods can converge to different sets of endmembers, which depend on the initialization. For the acquisitions with the phantoms, the results presented for GLUP and VCA are the best ones that were achieved. For the preliminary in vivo results, the sets of endmembers obtained with the different techniques are presented and discussed.

Chromophore Dilution

The endmembers estimated with GLUP, VCA, N-FINDR, and SSM-S are presented in Figure 7a-d, respectively. They allow us to discriminate the pure blue and green chromophores even if they do not exactly fit the reference spectra measured with a spectrometer (Figure 3). Indeed, as already mentioned, light absorption and diffusion in tissue impact the pure chromophore spectrum. The shape tendency is, however, kept. The shadow endmember is also plotted in each panel.

The abundance maps and the estimated mean concentrations are presented in Figure 8, for both the endmembers. All of these methods localize both of the pure chromophores well. The green abundance maps clearly highlight the green part, i.e., that on the right. The SSM-S achieves the best performance, as it leads to small green abundances for both parts of the blue phantom. Indeed, the mean concentrations are 0.05 and 0.07 for the blue and the diluted blue parts, respectively. The blue abundance maps are, nevertheless, also characterized by relatively large values in the green part. The SSM-S blue abundance map shows the smallest value of 0.26, and estimates a dilution value of 0.51, as compared to 0.53 for the ground truth. In conclusion, for this dataset, SSM-S/FCLS achieves the best performance.

Mixing of the Chromophores

The unmixing methods were also tested with mixed chromophores using the B-mix-G(Vevo) dataset. The endmembers extracted by the different methods are presented in 9. The extracted spectra differ from one strategy to another, but all blue spectra have similar evolutions, which matches the ground truth multispectral characteristics, in term of curve evolution, of this blue ink (Figure 3). The shadow endmembers are also plotted.

The abundance maps and the estimated mean concentrations are presented in Figure 10, for both of the endmembers. The blue abundance maps calculated with the endmembers extracted by GLUP and VCA localize the pure blue chromophore part (left) well, while the pure green part is more difficult to delimit. Indeed, the abundance values are similar in the central and right parts. For the N-FINDR and SSM-S blue abundance maps, the three parts are more visible. On the SSM-S abundance maps, the green part is characterized by a low concentration of 0.10, which is close to the ground truth value. For the green abundance maps, all of the methods allow discrimination of the blue part (left) and the mixed part (center). The VCA abundance map does not allow for clear discrimination of the mixed (central) and green (right) parts, while this is possible with the other methods. For the estimated mean concentrations in Figure 10, the first result to be noted is the blue concentration of the mixed chromophore estimated using VCA, which is equal to the ground truth of 0.42, as circled by the red dashed line. Nevertheless, the other concentrations estimated with VCA are larger than the ground truth concentrations: 0.89 for the green in the mixed chromophore instead of 0.67, 0.27 for the blue in the pure green chromophore instead of 0, and 0.20 for the green in the pure blue chromophore instead of 0. Consider now the performance of the SSM-S/FCLS, as circled with the orange dashed line in Figure 10. The estimated concentrations in the mix of chromophores are 0.35 for the blue and 0.73 for the green, instead of 0.42 and 0.67 for the ground truth, respectively. This corresponds to an error of less than 7%, which is already an interesting result in terms of the targeting of biological applications. Note that the estimated abundances where the ground truth concentration equals 0 were estimated at values less than 0.1 by SSM-S, which is a satisfying performance compared to the other methods. As a conclusion, for this dataset, SSM-S/FCLS achieves the best performance.

Preliminary In Vivo Results

The analyzed strategies were used on biological tissues to calculate the oxygen saturation rate (sO 2 ). This measurement is of great interest for various medical applications, such as the follow-up of tumors and the evaluation of tissue aging [START_REF] Su | Optoacoustic 3D visualization of changes in physiological properties of mouse tissues from live to postmortem[END_REF].

The pre-processing steps described in Section 3.1 were followed. Then, the endmember estimation methods were applied to extract two endmembers, which corresponded to the Hb and HbO 2 chromophores (the theoretical Hb and HbO 2 spectra are shown Figure 11). Finally, the abundance maps of both of these endogeneous chromophores were estimated using the FCLS algorithm. The sO 2 at each pixel was determined using the resulting abundance maps and Equation ( 16). These results were compared to the one provided by the Vevo LAZR. Indeed, the Vevo LAZR can provide abundance maps of the Hb and HbO 2 chromophores. We used these abundance maps with Equation ( 16) to calculate the oxygen saturation rates; see Figure 12. The authors are aware that Vevo LAZR Oxy-Hemo mode is not the real ground truth. However, this scanner represents a significant reference since it is typically used for sO 2 measurements in photoacoustic imaging. Figure 13, left column, shows the endmember spectra estimated with GLUP, VCA, N-FINDR, and SSM-S. Recall that the GLUP and VCA algorithms can converge to different sets of endmembers, depending on their initialization. The GLUP and VCA spectra represented in Figure 13 are those that then lead to the best sO 2 maps compared to the Vevo LAZR map, which are considered as the relative ground truth. FCLS was used with all of these endmember spectra, to estimate the abundance maps of Hb and HbO 2 . The sO 2 maps were then calculated using Equation ( 16); see Figure 13, center column. With this strategy, it appears that many pixels have strong sO 2 values, when VCA and N-FINDR algorithms are used. In addition, the endmembers extracted with both these algorithms are really similar, but not exactly the same, even if it is quite hard to see their differences on Figure 13. However, the sO 2 maps calculated with this strategy are different, which highlights the endmembers differences.

To characterize the performances of these methods, the relative deviations between each estimated sO 2 map (see Figure 13, center column) and the Vevo LAZR map calculated with Equation ( 16) (see Figure 12, right column) were computed. This relative deviation defines the absolute difference between the Vevo LAZR sO 2 values and each of the estimated sO 2 values, for each pixel. This is provided in %, as the sO 2 values are also in %. The relative deviation maps are presented in Figure 13, right column. Note that skin, tumor, and other biological tissues are present in the imaged area. The mean relative deviation was first calculated for all of these tissues. However, as the study was focused on the tumor only, the mean relative deviation was also calculated for the restricted regions delimited by the green rectangles in Figure 13. These relative deviations for all of these tissues and just for the tumor are in the black and green characters on the right in Figure 13, respectively. For both regions, the best estimations were reached with GLUP and SSM-S, which give similar relative deviations: 24.04% and 15.86% for GLUP, and 24.07% and 15.29% for SSM-S. The mean values of the endmembers are: 0.17, 0.39, 0.41, and 0.61 for Hb (in blue) and 0.16, 0.24, 0.24, and 0.48 for HbO 2 (in red), respectively for GLUP, VCA, N-FINDR, and SSM-S. The sO 2 maps estimated using the FCLS with the extracted endmembers and Equation ( 16) are highlighted (center). The relative deviation maps between these sO 2 maps and the Vevo LAZR map in Figure 12 are presented (right). The mean relative deviations for the whole image (black) and limited to the tumor (green) are given on the far right. The minimal relative deviations are circled in red. The image axes are in millimeters.

Discussion

This study experimentally evaluates the interest of using multispectral photoacoustic imaging to estimate chromophore concentration. To tackle such objectives, several methods coming from the remote sensing field have been tested to evaluate their potential performances in the photoacoustic domain. Principally, the constraints of Equation (3) can benefit also for chromophore concentration evaluation purpose. The interests and disadvantages of each method, highlighted by the in vivo proof of concept to evaluate sO 2 , are here discussed.

First, the variabilities of GLUP and SSM-S have to be analyzed. As already mentioned, GLUP and VCA can converge to different sets of endmembers, which depend on their initialization. Here, we focus only on GLUP as VCA, regardless, gives inaccurate sO 2 values. Figures 14 and15 illustrate this phenomenon with the results from 200 runs for the GLUP method. If GLUP can lead to small relative deviations, it should be noted that only 17 runs provided mean relative deviations less than 25%. In contrast, 59 runs provided mean errors greater than 50%. It is, however, interesting to note that the endmembers extracted with VCA and N-FINDR have the closest spectra to the theoretical ones [START_REF] Beard | Biomedical photoacoustic imaging[END_REF]; see Figure 16a. This observation should be considered with caution though, as the theoretical spectra shown do not take into account the characteristics of the acquisition system, which might modify the acquired spectra. As well as for the inks study, light absorption and diffusion in tissue also impact the pure chromophore spectrum. This could explain the bad results given by these two methods even if the extracted spectra seem the similar ones. Finally, considering its robustness, SSM-S is an interesting solution, as it converges to a unique solution while achieving a small relative deviation. However, it is important to note that the endmembers extracted with SSM-S do not look like the theoretical ones; see Figure 16a. As mentioned earlier, this could be due to the characteristics of the acquisition system and the light absorption and diffusion in tissue. In addition, in tumor tissues, pure Hb and HbO 2 are not the most frequently present chromophores. The Vevo LAZR sO 2 values (Figure 12) confirm this observation as they are mainly comprised between 25% and 75%, and then correspond to a mix of pure chromophores. This means that the endmembers that were actually extracted probably represent mixed chromophores instead of the pure Hb and HbO 2 . Regardless, the results obtained with the SSM-S/FCLS strategy are very encouraging in allowing accurate measurements of sO 2 with photoacoustic imaging.

Table 1 summarizes all of the results, the abundances are calculated with the FCLS algorithm, and the mentioned methods are the ones used to extract the endmembers. This table shows that the minimum mean relative deviation was obtained with the theoretical spectra, but the mean relative deviations remain large, as these spectra do not take the characteristics of the acquisition system into account as well as the light absorption and diffusion in tissue. Nevertheless, this result is interesting, as it is independent of the Vevo LAZR Oxy-Hemo mode, and it provides a convenient basis for comparing results. In this sense, given the proximity to the performance obtained with SSM-S/FCLS, these complementary experiments confirm the potential of this strategy to quantify photoacoustic in vivo data. GLUP/FCLS also have similar performance but only for some runs, as mentioned before. This strategy could then be more studied and compared to SSM-S/FCLS on various cases. However, the proposed strategies of hyperspectral imaging suffer from spectral coulouring effect [START_REF] Gehrung | Development of a blood oxygenation phantom for photoacoustic tomography combined with online pO 2 detection and flow spectrometry[END_REF][START_REF] Gröhl | Learned spectral decoloring enables photoacoustic oximetry[END_REF]. Indeed, in the proposed strategies, the endmembers are evaluated without taking into account any fluence evolution through the sample. Such effect has demonstrated to bias the sO 2 estimation and would required deeper study to be conducted. However, in the context of this work and its comparison to the sO 2 evaluation of the Vevo LAZR, we did not evaluate its impact.

These results are encouraging to be used for various clinical applications requiring chromophore quantification. Indeed, the evaluation of chromophore concentrations is of major interest for various applications, like concentration of contrast agent in the body [START_REF] Mienkina | Evaluation of Ferucarbotran (Resovist) as a Photoacoustic Contrast Agent[END_REF], calculation of the concentration of oxygenated and deoxygenated blood for oxygenation rate evaluation [START_REF] Mercep | Combined pulse-echo ultrasound and multispectral optoacoustic tomography with a multi-segment detector array[END_REF][START_REF] Singh | In vivo demonstration of real-time oxygen saturation imaging using a portable and affordable LED-based multispectral photoacoustic and ultrasound imaging system[END_REF], or cancer tissue evaluation [START_REF] Jnawali | Automatic cancer tissue detection using multispectral photoacoustic imaging[END_REF]. The endmembers extraction could be either automatic or manually selected by the user coupled with a fast inverse problem (FCLS).

Conclusions and Perspectives

Several methods for characterization of multispectral photoacoustic data have been considered, with the aim to estimate chromophore concentrations. These strategies mainly consist of two steps: [START_REF] Beard | Biomedical photoacoustic imaging[END_REF] an endmember extraction step, to estimate the spectrum of each pure chromophore to be quantified; and (2) an abundance calculation step, to evaluate their concentrations. For the first step, different algorithms were compared (GLUP, VCA, N-FINDR, and SSM-S). The second step was performed using the FCLS algorithm.

For nonbiological tissues and using the Vevo LAZR acquisition system, we demonstrated that the SSM-S/FCLS procedure provides accurate estimates of chromophore con-centrations for dilutions and mixtures. The positivity and sum-to-one constraints imposed by the FCLS, and coming from the remote sensing field, are interesting to take into account for chromophore concentration evalution. In some cases, the sum-to-one constraint is, however, questionable. It is a strong assumption, particularly when dilutions are imaged. To relax this constraint, a shadow endmember was introduced in this study, among the other possible strategies.

The methods were tested on a mouse tumor, where the sO 2 calculation was performed and compared to the sO 2 values calculated by the Vevo LAZR Oxy-Hemo mode. GLUP and SSM-S provided equivalent optimal performance to extract endmembers. GLUP algorithm should, however, be used with caution because of its variability. Nonetheless, it again highlights the interest of the constraints of Equation (3) in the photoacoustic field for chromophore concentration evaluation purpose. The SSM-S algorithm, which converges to a unique solution, appears to be more appropriate for endmember extraction, but further study on various cases should be done to quantify the performances of SSM-S/FCLS and GLUP/FCLS on various cases.

The endmember extraction and sum-to-one constraint relaxation still needs to be studied in more detail, but these primary results with biological tissues are encouraging. Moreover, the ground truth of the oxygenation rate in the imaged tumor was not known. It is, thus, difficult to know exactly which method is the most accurate to estimate the sO 2 values. In addition, considering the presence of only Hb and HbO 2 in the imaged area is a strong hypothesis. The water and other tissue absorbances should probably be taken into account. Further, other wavelengths can be tested to characterize their influence on the results. In the future, these experiments should be correlated with blood extraction and sO 2 evaluation using blood measurements to have the real ground truth.

Figure 1 .

 1 Figure 1. Multispectral photoacoustic image. (a) The signals are acquired at several wavelengths. Two pixels are identified: x 1 in green, and x 2 in purple. (b) Spectra of the pixels x 1 and x 2 : A 1 (in green) is a background pixel with low and flat photoacoustic amplitudes at all wavelengths, and A 2 (in purple) is a pixel with significant photoacoustic amplitudes (a.u., arbitrary units).

Figure 2 .

 2 Figure 2. Commercial multispectral photoacoustic system: Vevo LAZR.

Figure 4 .

 4 Figure 4. (a) Left: Blue/0.53 blue dilution/green colored phantom. Right: Blue/0.42 blue & 0.67 green mix/green colored phantom. (b) Spectra of each region of interest. Each spectrum corresponds to a specific pixel. These pixels are highlighted on the images (c) at λ = 680 nm. (a.u., arbitrary units) (c) Left: The B-Bdil-G(Vevo) dataset (mean value over all 8 wavelengths: 12.05 a.u.). Right: The B-mix-G(Vevo) dataset (mean value over all 8 wavelengths: 26.72 a.u.). The image axes are in millimeters.

Figure 5 .

 5 Figure 5. Performance validation procedure on an illustrative synthetic phantom. The validation here is for the estimation of the blue concentration.

Figure 6 .

 6 Figure 6. Mouse with breast tumors in the Vevo LAZR system. White and red arrows indicated the tumors and heartbeat detectors, respectively. In (a) the anesthesia mask can be seen (white block). The green square in (b) highlights the imaged region-of-interest.

Figure 7 .

 7 Figure 7. Endmembers estimated with (a) GLUP, (b) VCA, (c) N-FINDR, and (d) SSM-S (a.u., arbitrary units).

Figure 8 .

 8 Figure 8. Unmixing method results for the B-Bdil-G(Vevo) dataset with the algorithms indicated, combined with FCLS. Left columns: green endmember and unmixing results summary. Right columns: blue endmember and unmixing results summary. The values in black are the average normalized values as explain in Figure 5. The ground truth is shown in both cases, in the first row. The image axes are in millimeters.

Figure 9 .

 9 Figure 9. Endmembers extracted with (a) GLUP, (b) VCA, (c) N-FINDR, and (d) SSM-S (a.u., arbitrary units).

Figure 10 .

 10 Figure 10. Unmixing method results for the B-mix-G(Vevo) dataset with the algorithms indicated, combined with FCLS. Left columns: green endmember and unmixing results summary. Right columns: blue endmember and unmixing results summary. The values in black are the average normalized values as explain in Figure 5. The ground truth is shown in both cases, in the first row. The image axes are in millimeters.

Figure 11 .

 11 Figure 11. Theoretical Hb and HbO 2 spectra (mean values are: 0.43 and 0.35, respectively).

Figure 12 .

 12 Figure 12. Oxygenated hemoglobin (HbO 2 ) and deoxygenated hemoglobin (Hb) concentration maps and sO 2 map (in %) calculated with the Vevo LAZR concentration maps and Equation (16), from left to right, respectively. The image axes are in millimeters.

Figure 13 .

 13 Figure 13. Results for (a) GLUP, (b) VCA, (c) N-FINDR, and (d) SSM-S. The extracted endmembers are shown (left).The mean values of the endmembers are: 0.17, 0.39, 0.41, and 0.61 for Hb (in blue) and 0.16, 0.24, 0.24, and 0.48 for HbO 2 (in red), respectively for GLUP, VCA, N-FINDR, and SSM-S. The sO 2 maps estimated using the FCLS with the extracted endmembers and Equation (16) are highlighted (center). The relative deviation maps between these sO 2 maps and the Vevo LAZR map in Figure12are presented (right). The mean relative deviations for the whole image (black) and limited to the tumor (green) are given on the far right. The minimal relative deviations are circled in red. The image axes are in millimeters.

Figure 14 .

 14 Figure 14. Endmembers extracted with GLUP. Left column: sets of extracted endmembers. Right column: mean spectra and standard deviations.

Figure 15 .

 15 Figure 15. Mean relative deviations on sO 2 maps estimated with GLUP over 200 runs.

Figure 16 .

 16 Figure 16. Results using the theoretical absorbance spectra of Hb and HbO 2 . (a) The theoretical spectra (mean values of Hb and HbO 2 are: 0.43 and 0.35, respectively), (b) sO 2 map estimated using the FCLS with the theoretical spectra and Equation (16), and (c) the relative deviation maps between the sO 2 map and the Vevo LAZR map in Figure12. The image axes are in mm.

Table 1 .

 1 Summary of the average sO 2 relative deviation values. The abundances are estimated using the FCLS algorithm.

	Endmember Extraction	Whole Image	Limited to Tumor
	GLUP	24.04%	15.86%
	VCA	58.79%	38.96%
	N-FINDR	76.81%	51.63%
	SSM-S	24.07%	15.29%
	Theoretical spectra	21.52%	13.33%

Sensors 2021, 21, 3366. https://doi.org/10.3390/s21103366 https://www.mdpi.com/journal/sensors

Acknowledgments:

The authors want to thank PHENOCAN ("Investissements d'Avenir") for imaging devices (ANR-11-EQPX-0035 PHENOCAN) and Federica Cavallo of Università degli Studi di Torino for the mouse.

Data Availability Statement:

The data presented in this study are available on request from the corresponding author.

Institutional Review Board Statement:

The study was conducted according to the animal care guidelines of the European Union, and were validated by the local Animal Ethics Evaluation Committee and the French Ministry of Higher Education and Research (C2EA-15 and 02296.02).

Informed Consent Statement: Not applicable.

Conflicts of Interest:

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

FCLS

Fully Constrained Least-Square GLUP Group Lasso with Unit sum and Positivity constraints PCA Principal Component Analysis SSM-S Spatio-Spectral Mean-Shift VCA Vertex Component Analysis


Amaury Nègre

James L Crowley

Christian Laugier

Scale Invariant Detection and Tracking of Elongated Structures

This paper describes a method for the detection of tracking of elongated structures that is robust under changes of scale and orientation. This method is based on extending the concept of scale invariant natural interest points to include elongated ridge structures. An operator is proposed that directly detects ridge points and provides an estimation of their elongation and orientation. A tracking process is used to follow elongated features over time and to robustly observe changes in scale and orientation. Changes in scale are used to directly estimate time to contact. Experimental results demonstrate that the method works well in cluttered scenes that are typical of urban environments.

Introduction

Scale invariant image descriptions have been studied since the late 1970s. Burt [START_REF] Burt | The laplacian pyramid as a compact image code[END_REF] proposed a rapid algorithm to compute a multi-scale Laplacian pyramid. Crowley [START_REF] Crowley | A representation for shape based on peaks and ridges in the difference of low-pass transform[END_REF] proposed a scale invariant algorithm for the Laplacian pyramid and showed that such representations could be used to compute scale invariant image descriptions composed of peaks (scale invariant interest points) and ridges. Lowe adopted scale invariant peaks in the Laplacian Pyramid as the basis for a Scale Invariant Feature Transform (SIFT) [START_REF] Lowe | Object recognition from local scale-invariant feature[END_REF]. This operation is now widely used in image matching, tracking and object categorization. Despite the widespread success of the SIFT operator, problems remain. For many real world objects, the detected feature points are not centered on the "physical" object they represent. More importantly in the presence of elongated shapes, multiple interest points will appear along the central axis of the objects and the position of such points along the object principle axis can be highly unstable.

In this article, we generalize the notion of scale invariant interest points to include elongated "ridge" lines. We develop a multi-scale version of the ridge-line detector proposed by [START_REF] Tran | A method for ridge extraction[END_REF], and show that such structures may be tracked with a particle filter to estimate the motion of an object. We demonstrate that the rate of change of scale of a ridge may be used to directly estimate the time to contact for an object in the scene. Such a method may be used for visual motion estimation and obstacle detection in cluttered urban environments.

2 Scale invariant ridge detector

Algorithm description

The SIFT operator uses the scale at which a peak in a Laplacian is a local maximum to determine the most appropriate scale for local image description. This maximum is referred to as the characteristic scale. We note in passing that such a characteristic scale is not limited to peaks in the Laplacian but can be determined at most points in an image [START_REF] Crowley | View invariant object recognition using coloured receptive fields[END_REF]. A precise estimate of the characteristic scale can be obtained at any image point by interpolating the Laplacian values over a range of scale using a 3rd order interpolation function.

Two problems can arise with the use of scale invariant feature points 1. The Laplacian exhibits local extrema in the center of object but also on the edges, where intrinsic scale is not meaningful; 2. In the case of elongated object such as pedestrians or light poles, the position of detected points may be unstable along the axis of the object.

These phenomena can give rise to spurious or unstable feature points that degrade the reliability of a tracking or image matching.

To solve these problems, we extend the concept of natural interest points to include ridge lines. Such ridge lines naturally occur along the central axis of any elongated object, at a scale that corresponds to the visual width. Just as natural image points can be detected and tracked, so can natural ridge lines. Unlike interest points, ridge lines can be used to directly provide local orientation without the need for local histograms of gradients.

The detection algorithm operates as follows:

1. For each pixel of the Laplacian pyramid we first remove the edge response by eliminating pixels where the ratio between the Laplacian and gradient values is greater than a threshold. 2. We then determine the ridge direction using a local Hessian [START_REF] Tran | A method for ridge extraction[END_REF]. This direction is normal to the principal curvature direction given by the eigen-vector associated to the biggest eigen-value of the Hessian matrix H =

∂ 2 f ∂x 2 ∂ 2 f ∂xy ∂ 2 f ∂xy ∂ 2 f ∂y 2 where ∂ 2 f ∂x 2 , ∂ 2 f
∂xy , ∂ 2 f ∂y 2 are the second derivatives of the Image. 3. For each pixel we search the length l that maximizes the following score function:

S(X, l, u) = k=0..l (|L(X + k • u X )| + |L(X -k • u X )| -2 • |L(X + k • u X ) -L(X -k • u X )| -2 • L min ) (1)
Where L(X) is the Laplacian value at pixel X and u is the previously computed direction. This score is high when many pixels along the line have a high

Laplacian (|L(X + k • u X )| + |L(X -k • u X )|),
and is maximal at the center of the object (as the term

|L(X + k • u X ) -L(X -k • u X )| increases near contrast boundaries).
The term L min makes the score decrease when k is greater than the object's characteristic length. It should be noted that the cost of the score function is only proportional to the maximal searched length. 4. We then search local maxima of the previous score in the pyramid to obtain best segments.

When the feature is not a perfect line, the local direction of the principal curvature is not strictly aligned with the feature. The score can then be improved by scanning in neighboring directions.

A typical example of ridge segments detected in an image of an urban environment is shown in Figure 2.1. The detected segments are represented by red ellipses, the main axis of the ellipses represents the segment's half-edges, and the width of the ellipses represents the object scales. The interesting point is that all detected segments are centered on the object they represent and the detected scales fit very well object's size.

Performance evaluation

We used the repeatability criterion as described in [START_REF] Schmid | Comparing and evaluating interest points[END_REF] to evaluate the stability of this method for ridge segments detection to various image transformations. The repeatability criterion is determined as the proportion of image features that can be matched when an image is subjected to a known transformation. For ridge segments, correspondence is determined using the the center position, orientation and scale of the ridge segment. The maximal accepted matching distance (the validation gate) is proportional to the scale of the ridge segment.

We compared the the repeatability of ridge segments to interest points detected using a Harris detector [START_REF] Harris | A combined corner and edge detector[END_REF] computed at different scales as well as the SIFT detector [START_REF] Lowe | Object recognition from local scale-invariant feature[END_REF]. The results are shown in Figure 2. In this experiment we observe that the ridge segments are not as much robust as rotation change, but actually the low performance comes from the localization error that can be high along the ridge direction. As expected, the change of scale doesn't affect the detector's performances. The time to contact "τ " (also called time-to-collision or time-to-crash) can provide important information for the detection and avoidance of obstacles. Time to contact can be seen as directly determining the distance between two objects in a termporal space. It has been shown in [START_REF] Negre | Real-time time-to-collision from variation of intrinsic scale[END_REF] that the time to contact between the camera and a visible obstacle can be computed in the image space by measuring the variation of the characteristic scale.

Let "s" represent the characteristic scale of an image feature. In this case, τ can be approximated by :

τ = s ∂s ∂t (2)
The difficulty with this method is to identify and measure the change of scale of an image feature in the video sequence. The characteristic scale of a ridge segment provides exactly such a measure. Tracking ridge segments makes it possible to estimate both the visual motion of a feature and the time to contact.

Tracking Ridge Segments

We use the score function, described in 1, as a measure function for tracking ridge segments. Because this function is non-Gaussian and non linear, we have employed a particle filter for the tracking process [START_REF] Doucet | Sequential Monte Carlo methods in practice[END_REF]. The particle filter uses a set of samples (particles) to represent a probability distribution resulting for a Bayesian filtering process. Formally, if we let X t be the state of the target at time t, Z t be the observation of the target at time t, then the goal of the particle filter is to estimate :

P (X t |Z t1 ...Z tn )
To this end, the filter relies on two models : (1) the observation model P (Z t | X t ) that predicts a target observation, and (2) a displacement model P (X t k | X t k-1 ) that predicts the motion of the particles.

Target tracking is implemented as follows: we use a set of particles to characterize the target's position and the target's speed in the image coded by three vectors, each vith 3 dimensions :

c : the segment's center position in the 3D Scale-Space; v : the segment's center speed ; r : the 3D vector between the center and an extremity (called "half-edge" in the following).

X =   c v r  
When a target is identified in an image for the first time, all particles are initialized around the detected position and with a zero speed. Next, each camera's image (Z t ) is used as an observation to update the particle filter. The observation model takes into account the score function described in Equation 1:

P (Z t |X t )) ∝ S(c t , ||r t ||, r ||r t || )
In between two images, we consider that the target's center is subjected to a Gaussian acceleration a and that the half-edge can also be affected by a Gaussian noise n. Each particle is then updated using the following model :

X t+∆t =   c t+∆t v t+∆t r t+∆t   =   c t + (v t + a • ∆t) • ∆t v t + a • ∆t r t + n • ∆t  
The output of the particle filter is a probabilistic estimation of the target state. In practice, to obtain a single state of the target, we compute the average pose of all particles weighted by their probabilities. It can be noted that the tracking is automatically ended either when the target get out the image frame or when the observation probability will decrease. As we track many segment at same time, we also add a mechanism to fuse differents segments that converge.

The global schema of the detection and tracking of scale invariant segments is shown on Figure 3.

Experimental results

Simple tracking

To demonstrate the utility of this method for tracking, we consider a simple video sequence with a black rectangle printed on a white sheet. The tracking result is shown on Figure 4(a). We can see that the estimated segment (red ellipse) is centered on the black rectangle, but that the length does not precisely fit the rectangle's length. None the less, the rectangle scale (height) is precisely approximated. To evaluate this precision, the height of the rectangle has been measured by hand in each of these images and used to estimate the actual time-to-contact. The results are plotted in Figure 4(b)(c). As expected the two curves fit very well.

Tracking Ridge Segments in Urban Scenes

To demonstrate the efficiency of the scale invariant segment detector and tracking, we have tested our algorithm on real world sequences. To obtain real-time performance, we implemented the algorithm on a Graphic Processing Unit (GPU), which made it possible to detect and track simultaneously 64 segments in 640x480 images at 17 frames per second. Some results are shown in Figure 5, tracked segments are represented by red ellipse and each segment's trajectory in the image is drawn with a black curve. In this sequence, we noted that the segment's tracker locks well on to road lines, components of the cars, trees and poles. For a more detailed analysis, we plotted the evolution of scale of the pole located in the center right part of the image (Figure 6(a)). As the camera moves forward and backward, the object's size increases and decreases which is approved by the tracker. The approximated time-to-contact (Figure 6(b)) is also appreciated as it begins positive and decreases (which mean the obstacle is approaching) and next it swaps negative, which mean the obstacle is getting away).

Conclusion

Ridge segments are a natural complement to interest points. Ridge segments are particularly useful in scenes, such as urban environments that contain many elongated objects. In this paper, we have have described a method for detecting ridge segments and shown that such a method can provide good repeatability under changes in orientation and scale. We have demonstrated a tracking algorithm using a particle filter for tracking ridge segments, and have demonstrated that such tracking can be used to directly estimate time to contact for elongated objects observed from a moving camera. Because this detector uses the center of an object, (as opposed to an edge or corner) it can provide improved stability over large changes in viewing angle for elongated objects.

Fig. 1 .

 1 Fig. 1. Scale invariant segment detector in an urban environment, the detected segments are represented by red ellipses.

  (a) Original test image and detected segments. (b) Ridges in rotated test image. (c) Ridges in Scaled test image. Repeatability with respect to scale factor.

Fig. 2 .

 2 Fig.2. Repeatability of Harris, SIFT and Ridge-Segment detector with respect to rotation and scale change. We can see that the segment detector obtain moderate score for rotation changes but gives similar results than the SIFT detector with change of scale.

3

  Estimating time-to-contact using the change of scale in ridge segments 3.1 Introduction of the Time-to-conctact

Fig. 3 .Fig. 4 .

 34 Fig. 3. Global schema of scale invariant segment detection and tracking.

Fig. 5 .

 5 Fig. 5. Segment detection (right) and tracking (left) in an urban context. Red ellipses with blue axes represent the tracked segments, black curves represents segment's trajectories in the image.

  Evolution of scale for the segments localized in the white pole.

  Estimated Time-To-Contact of the white pole.

Fig. 6 .

 6 Fig. 6. (a) Evolution of the scale of the tracked segments localized on the white pole (center right of images). As the car moves forward and backward, the scale increases and decreased. The estimated Time-To-Contact (b) is well approximated : at the beginning it is positive and it decreases (which mean the obstacle is approaching) and next the TTC swaps negative (which means the obstacle is getting away).