Ahmet Iscen

Giorgos Tolias

Philippe-Henri Gosselin

Hervé Jégou

A comparison of dense region detectors for image search and fine-grained classification

Keywords: dense keypoints, image retrieval, fine-grained classification, Zernike polynomials

We consider a pipeline for image classification or search based on coding approaches like Bag of Words or Fisher vectors. In this context, the most common approach is to extract the image patches regularly in a dense manner on several scales. This paper proposes and evaluates alternative choices to extract patches densely. Beyond simple strategies derived from regular interest region detectors, we propose approaches based on superpixels, edges, and a bank of Zernike filters used as detectors.

The different approaches are evaluated on recent image retrieval and fine-grain classification benchmarks. Our results show that the regular dense detector is outperformed by other methods in most situations, leading us to improve the state of the art in comparable setups on standard retrieval and finedgrain benchmarks. As a byproduct of our study, we show that existing methods for blob and super-pixel extraction achieve high accuracy if the patches are extracted along the edges and not around the detected regions.

I. INTRODUCTION

L OCAL image description is a popular research topic in computer vision, as it is involved in many applications such as image classification and particular object detection. Extracting local descriptors from an image consists of two steps. The detection step selects regions of interest, which are normalized into fixed-size patches. The description step produces a vector representation for each of the detected patches. The SIFT descriptor [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] and its RootSIFT extension [START_REF] Arandjelovic | Three things everyone should know to improve object retrieval[END_REF] have been shown to perform very well for most applications. Many descriptors have been introduced in the last years to improve the description speed or descriptor compactness, such as SURF [START_REF] Bay | SURF: Speeded up robust features[END_REF], CHOG [START_REF] Chandrasekhar | CHoG: compressed histogram of gradients[END_REF] or BRIEF [START_REF] Calonder | Brief: Binary robust independent elementary features[END_REF]. The matching accuracy is improved by learning the descriptor design [START_REF] Winder | Learning local image descriptors[END_REF][59] [START_REF] Simonyan | Learning local feature descriptors using convex optimisation[END_REF]. In this paper, we are solely interested in the detection stage and therefore adopt the gold-standard SIFT and RootSIFT.

The early works in this line of research have focused on the detection of sparse interest points and regions, typically producing a few thousand descriptors per image. These approaches aim at extracting distinctive and repeatable image parts, such as blobs and corners [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] [START_REF] Mikolajczyk | A comparison of affine region detectors[END_REF][31] [START_REF] Bay | SURF: Speeded up robust features[END_REF], offering covariance properties: the same regions should be detected under some geometrical transformations. From a historical perspective, the choice of sparse representations was arguably motivated by the lack computational and memory resources. Although such methods perform well and are still widely used for image matching, they are not competitive in other application scenarios such as image classification.

INRIA, Campus Universitaire de Beaulieu 35042 Rennes Cedex FRANCE Local feature detection has recently shifted towards denser techniques. Dense sampling is an easy way to provide a large number of patches and a better coverage of the objects of interest. Fei-Fei and Perona [START_REF] Fei-Fei | A bayesian hierarchical model for learning natural scene categories[END_REF] were the first to show that dense patches leads to better classification accuracy. Nowak et al. [START_REF] Nowak | Sampling strategies for bag-offeatures image classification[END_REF] argue that the key parameter for classification is the number of extracted patches. Similar conclusions hold for other tasks, like fine-grained classification [START_REF] Gosselin | Revisiting the Fisher vector for fine-grained classification[END_REF] or action recognition in videos [START_REF] Wang | Action recognition by dense trajectories[END_REF] [START_REF] Jain | Better motion for better action recognition[END_REF]. Recent works [START_REF] Zhao | Oriented pooling for dense and non-dense rotation-invariant features[END_REF] [START_REF] Delhumeau | Revisiting the VLAD image representation[END_REF][51] [START_REF] Tolias | Visual query expansion with or without geometry: refining local descriptors by feature aggregation[END_REF] also evidence that methods for image and particular object retrieval, which traditionally rely on sparse regions typically extracted with the Hessian-Affine detector [START_REF] Mikolajczyk | Scale and affine invariant interest point detectors[END_REF], are improved when using a larger set of descriptors.

Nevertheless, regular dense sampling has serious limitations. Uniform sampling of patches ignores the image structure and extracts many uniform and uninformative patches. Additionally, the position of the features is less or not repeatable. This prevents the image engine from employing a spatial verification method, such as RANSAC [START_REF] Fischler | Random sample consensus[END_REF] [START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF], which typically filters out many outliers by enforcing the spatial consistency of the detected regions.

In this paper, our goal is to develop and evaluate dense detection strategies for image retrieval and fine-grained image classification. Our motivation is similar to that of Tuytelaars when she introduced "dense interest points" [START_REF] Tuytelaars | Dense interest points[END_REF]: we consider solutions in between localized sparse interest points and dense strategies, in order to produce a large number of localized regions. We depart from using traditional evaluation metrics for detector evaluation, which do not reflect the final goal. For instance, the repeatability score reflects the effectiveness in detecting inliers, but is not directly related to the determination of class membership. Instead, we evaluate the performance with the metrics employed for the target application scenario: mean average precision for image retrieval and accuracy for fine-grained classification. We stress that better localized dense patches (see Figure 1) are more important for these tasks than in traditional image classification: the objects are more repeatable and distinguishing between two classes often rely on tiny details that suffer from being loosely localized.

We make the following contributions in this context of dense detection for image retrieval and fine-grained classification:

1) We propose strategies derived from standard interest point detectors (Harris, Hessian and DoG) to extract patches densely. In particular, we modify the detection process by relaxing the standard local maxima criterion so that it focuses on edges and not only corners. This, jointly with the optimization of the scaling factor, is shown to be a key to achieve higher performance.

2) We depart from the typical choice of fitting an ellipse to describe a region of interest extracted by blob detectors (MSER or super-pixels). Instead, we sample patches at several scales along the region's borders. This increases the performance significantly when considering a large number of patches per image. We also show that sampling on the edges produced by a state-of-the-art edge detector [START_REF] Dollár | Structured forests for fast edge detection[END_REF] offers competitive performance. 3) Finally, we propose two novel response filters to select the patch locations. First, we propose to use a bank of Zernike polynomials [START_REF] Zernike | Beugungstheorie des schneidenver-fahrens und seiner verbesserten form, der phasenkontrastmethode[END_REF] as detectors. These filters have been proposed to construct descriptors [START_REF] Revaud | Improving Zernike moments comparison for optimal similarity and rotation angle retrieval[END_REF], but to our knowledge not as detectors. Our second strategy is descriptor-oriented: the response for each pixel is simply the norm of local descriptor associated with the patch centered at this location. These two new approaches appear to perform best in most cases of our experiments. This paper is organized as follows. Section II discusses previous work related to the purpose of our study. Section III describes existing approaches which are part of our evaluation, while in section IV we present modified schemes of existing detectors and the use of pseudo-Zernike polynomials for dense feature detection. Finally in section V we present the outcome of this study on 2 datasets on image retrieval and 4 datasets on fine-grained classification.

II. RELATED WORK

A great deal of works has focused on the detection of local interest points [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF][6][34] [START_REF] Matas | Robust wide baseline stereo from maximally stable extremal regions[END_REF]. They are rather designed to be appropriate for image matching applications. The standard evaluation metrics for such detectors [START_REF] Mikolajczyk | A comparison of affine region detectors[END_REF], e.g. repeatability, are designed to reflect sufficiency for matching applications. However, this could be far from the final application tasks we focus in this paper, i.e. fine-grained classification.

Typically, local interest points are tuned to detect image structures, such as corners and blobs, highly distinctive and repeatable. One of the exceptions is the work of Mikolajczyk et al. [START_REF] Mikolajczyk | Shape recognition with edge-based features[END_REF] where they focus on edges to represent objects. Similarly, in this study we argue that such a choice is beneficial for image retrieval and fine-grained classification.

Tuytelaars's work on dense interest points [START_REF] Tuytelaars | Dense interest points[END_REF] has a motivation close to ours. Initially, a dense grid of patches is considered, and then for each point, local refinement of its location and scale parameters is performed. In particular, the point of maximum interestingness is selected within a bounded area. This choice suffers from quantization artifacts, since a true local maximum might not exist in this search area. Apparently, a large number of patches are located on smooth and uninformative regions.

A previous survey closely related to this study is the work of Nowak et al. [START_REF] Nowak | Sampling strategies for bag-offeatures image classification[END_REF]. They are also not interested in the repeatability of patches, but directly measure classification accuracy. According to them, detectors designed to obtain high repeatability perform the same as randomly selected patches. They observe that performance is, up to some extent, an increasing function of the number of points per image. Inspired by their finding, we conduct a similar evaluation and try to study different ways to control the average number of patches. Similarly, Avrithis and Rapantzikos [START_REF] Avrithis | The medial feature detector: Stable regions from image boundaries[END_REF] do not restrict comparison to repeatability, but further compared image retrieval performance, while considering the average number of points as a crucial parameter.

One of the parameters we consider in our study is the size of the measurement region for local descriptor extraction. This has been given particular attention by the recent work of Simonyan et al. [START_REF] Simonyan | Learning local feature descriptors using convex optimisation[END_REF], who shows that introducing a scaling factor increases retrieval performance. Similarly, the improved Hessian-Affine detector by Perdoch et al. [START_REF] Perdoch | Efficient representation of local geometry for large scale object retrieval[END_REF] uses a larger measurement region to improve the performance.

In regular sampling local descriptors are agnostically extracted on a dense grid, and as a consequence many uninformative descriptors are derived from smooth regions. A significant improvement is achieved by filtering out descriptors with low 2 -norm [START_REF] Gosselin | Revisiting the Fisher vector for fine-grained classification[END_REF]. We combine this approach with examined detectors and study its effect on different setups.

Convolutional neural networks are very effective for image classification [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF][14] and detection [START_REF] Girshick | Rich feature hierarchies for accurate object detection and semantic segmentation[END_REF]. However, for the tasks addressed in this paper, namely fine-grained classification, image and particular object, this is not (yet) the case. Fisher vectors based on a dense representation and employing medium-sized codebooks and spatial coding achieved much better results in the fine-grained challenge [START_REF] Gosselin | Revisiting the Fisher vector for fine-grained classification[END_REF]. For image retrieval, recent work demonstrates that CNNs [START_REF] Babenko | Neural codes for image retrieval[END_REF] achieve excellent performance for aggressive operating points for compact representations. However, their best reported performance is lower than that of the state-of-the-art Fisher vectors, and outperformed by a large margin by state-of-the-art methods like the selective match kernels [START_REF] Tolias | To aggregate or not to aggregate: Selective match kernels for image search[END_REF].

III. BACKGROUND

In this section we briefly describe existing methods for interest point detection or dense patch sampling, that are part of our study. In all those methods, a local descriptor is extracted from each region of interest and an image is represented by a set of such descriptors. We visualize detected features for all the detectors examined or proposed in this work on Figure 2.

A. Interest points

Harris-Laplace detector. Harris detector localizes corners, based on the fact that gradient values will change in multiple directions around a corner. It uses a scale adapted version of the second moment matrix, known as Harris matrix [START_REF] Harris | A combined corner and edge detector[END_REF]:

M " σ 2 D ¨gpσ I q ˚" L 2 x px, σ D q L x L y px, σ D q L x L y px, σ D q L 2 y px, σ D q  , (1)
where σ D is the differentiation scale, σ I is the integration scale and L z is the derivative computed in z direction. Differentiation scale σ D is used to compute the local derivatives with Gaussian kernels, and a Gaussian window with a size σ I is used to smooth and average the neighborhood around the point. The eigenvalues of this matrix represent the gradient changes in two directions. Consequently, if one eigenvalue is large while the other is small, there exists an edge, whereas if both eigenvalues are large, there exists a corner. The interestingness of a point is captured by the cornerness function, defined as cornerness " detpM q ´αtrace 2 pM q, where α is usually set to 0.05. Extending Harris detector to be scale invariant, Mikolajczyk et al. [START_REF] Mikolajczyk | Scale and affine invariant interest point detectors[END_REF] use Laplacian-of-Gaussian (LoG) response and detect local extrema over multiple scales to perform scale selection [START_REF] Lindeberg | Feature detection with automatic scale selection[END_REF]. This is the well known Harris-Laplace detector. Only points with cornerness value higher than threshold τ are retained and final point locations are chosen by a local maxima search procedure.

Hessian-Affine detector [START_REF] Mikolajczyk | A comparison of affine region detectors[END_REF], similar to Harris detector, detects image locations that have large derivatives in both directions. Point locations are selected as local maxima of the Hessian matrix determinant, which now constitutes the interestingness measure. The Hessian matrix is defined as

H " " L xx px, σ D q L xy px, σ D q L xy px, σ D q L yy px, σ D q  , (2)
where L zz is the second order partial derivative. All points with interestingness below threshold τ are discarded. Although Hessian and Harris detectors are quite similar, the detected points may be slightly different. In particular, unlike Harris, Hessian detector tends to select locations with texture variations, in addition to corners.

The affine shape of the point is estimated by the eigenvalues of the second moment matrix M . An iterative procedure modifies the point's location, scale and shape until the estimated affine transform is able to map the detected region into one that has equal eigenvalues of its second moment matrix.

MSER Maximally Stable Extremal Regions were introduced by Matas et al. [START_REF] Matas | Robust wide baseline stereo from maximally stable extremal regions[END_REF] as a feature detector that tends to localize blobs. An extremal region has all its intensity values greater (or less) than the outer region boundary pixels. Then, a sequence of nested extremal regions is considered. Along this sequence scale change between neighboring regions is estimated. Maximally stable are the local minima of this quantity. In this fashion, nested regions are also likely to appear. The detected regions can have very irregular shapes, therefore an ellipse is fitted to detected regions in order to extract local descriptors. A parameter ∆ controls the locality of the scale change computation.

Difference of Gaussians (DoG) was originally used for local feature detection by Lowe [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF]. DoG is an efficient way to approximate the Laplacian of Gaussian and to detect edges at various image scales. A Gaussian kernel is used to create multiple blurred versions of the image per octave. Simple subtraction of two consecutive blurred images produces the DoG response. Interest points are located in the scale-space as local-maxima in a 3D search area of size 3.

The original SIFT detection algorithm [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] employs DoG, but further applies additional steps to filter out points belonging to edges or low-contrast regions. However, in our experiments with DoG detector, we keep all the points for a denser representation.

B. Dense patch sampling

Regular grid dense sampling. In contrast to interest points, dense sampling methods give less importance to high repeatability and try to provide a dense coverage of the depicted objects. The most popular method is to sample points on a regular grid, every δ xy pixels. Depending on the application δ xy can be really small, such as 3, or quite larger, such as 16.

In order to provide some scale tolerance, different scales are considered by following the same procedure at n σ multiple scales of the image. All the patches from different scales are pooled together in the end. A typical value for n σ is 5, with 2 scales per octave. We adopt this choice in our experiments.

Dense interest points were introduced by Tuytelaars [START_REF] Tuytelaars | Dense interest points[END_REF] as a hybrid solution to trade-off between sparse interest point detection and dense sampling. Instead of selecting the center pixel of grid cell, as in regular dense sampling, they conduct local search inside each cell over spatial and scale space. The point with maximum response is kept per cell. Selected points are not necessarily local maxima. Since a single interest point is selected from each cell, patches are very likely to be localized on smooth regions. Moreover, quite a few patches are localized on the cell borders. In the original work, a large local search area of 16 ˆ16 pixels ˆ8 scale levels is used, while in our study we experiment with smaller search areas in order to supply more points.

IV. PROPOSED METHODS

Our purpose is to provide a dense set of patches without targeting high repeatability. We rather focus on accuracy in retrieval and fine-grained classification. The proposed methods try to provide a nice coverage of the depicted objects, while focusing on edges in addition to corners and blobs.

In this section we present the proposed approaches for dense patch representation. We relax the Harris-Laplace detector to detect mostly edges. We propose to use convolution by Zernike polynomials [START_REF] Zernike | Beugungstheorie des schneidenver-fahrens und seiner verbesserten form, der phasenkontrastmethode[END_REF] as a response function to select point locations. Moreover, we modify the well-known region detectors, such as MSER, to sample points on edge maps, or the edges of the detected regions. Finally, we propose a detector that exploits the 2 -norm of already extracted local descriptors to select the patches. In the following, whenever we mention that local maxima search is performed, a 3 ˆ3 neighborhood is used.

A. Relaxing Harris-Laplace detector

Harris-Laplace detector is designed to detect corners. We propose to replace Harris-Laplace cornerness function with Frobenius norm in order to select other points in addition to corners. We also propose to adjust the local extrema selection criterion to sample denser points for classification context.

The first modification is to adjust the standard cornerness function to sample denser points. We estimate the energy of the Harris matrix by its Frobenius norm. It is defined as the square root of the sum of the absolute squares of its elements. This offers the ability to produce high response for structures with high energy other than corners.

As our second modification, we propose to modify the local maxima criterion. Regularly, keypoints are selected as local maxima of some interestingness measure. Instead of using the 8-neighborhood of a pixel as local maxima criterion as the standard procedure, we propose to use multiple 2neighborhoods in different directions, as shown in Figure 4. We select the points that are local maxima in any of those neighborhoods. In this fashion, we retrieve more points along the edges of objects or object parts.

In order to investigate the contribution of different modifications, we end up with four different methods. These are, the regular Harris-Laplace (Harris), Frobenius norm based detector (Frobenius), Harris-Laplace with relaxed local maxima selection (relaxed-Harris) and Frobenius norm based detector with relaxed local maxima selection (relaxed-Frobenius).

In the example of Figure 3, we observe that points are located on edges as well as corners using the Frobenius norm. The relaxed local maxima criterion also provides more points outlining the shapes on the image.

B. Zernike

Zernike polynomials were originally introduced by Zernike [START_REF] Zernike | Beugungstheorie des schneidenver-fahrens und seiner verbesserten form, der phasenkontrastmethode[END_REF] and are traditionally used in optics and also as shape descriptors [START_REF] Khotanzad | Invariant image recognition by Zernike moments[END_REF] [START_REF] Revaud | Improving Zernike moments comparison for optimal similarity and rotation angle retrieval[END_REF]. They form a set of orthogonal basis functions defined on the unit disc. We focus on the use of pseudo-Zernike polynomials [START_REF] Bhatia | On the circle polynomials of Zernike and related orthogonal sets[END_REF] which have proven to be more robust to image noise [START_REF] Teh | On image analysis by the methods of moments[END_REF]. They are defined as

v l n " R nl prqe ilθ , (3)
with

R nl prq " n´|l| ÿ s"0 p´1q s p2n `1 ´sq! s!pn ´|l| ´sq!pn `|l| `1 ´sq! r n´s , (4
) where n is the order and l is the repetition. In Figure 5 we present pseudo-Zernike polynomials of order 1 up to order 4. Starting from low frequency patterns they end up being more complex patterns and more localized.

We normalize pseudo-Zernike functions to a 2D rectangular patch (patch width equal to 11 in all our experiments) and convolve the input image using them. Each pseudo-Zernike function is used as a filter and this response constitutes the interestingness measure for detection of point locations. Local maxima and local minima are detected for each filter independently. We now define the maximum number of patches to be extracted per image as N z . This is used to define the available capacity per filter, in order to retain only the strongest detections. Capacity is uniformly shared among different filters, and also among maxima and minima. We apply the same procedure in n σ " 5 scales, similarly to regular sampling. The ratio of capacities of two consecutive scales is 2. That is equal to the ratio of the down-sampled image areas. Filter responses are ranked per filter, and point locations are selected until each filter capacity is filled or until there are no more local extrema left.

We employ N f filters, which detect complementary structures as being orthogonal and produce high responses for different structures. In contrast to previous feature detectors, we do not focus particularly on detecting corners, blobs or edges and claim that all such structures are useful for image representation. For given capacity, the more filters we use the stronger features we select per filter. In Figure 6, we show the responses produced by various filters. Observe, from the responses and from the filters of Figure 5, that lower order polynomials will detect vertical edges, horizontal edges and

v ´1 1 v 0 1 v 1 1 v ´2 2 v ´1 2 v 0 2 v 1 2 v 2 2 v ´3 3 v ´2 3 v ´1 3 v 0 3 v 1 3 v 2 3 v 3 3 v ´4 4 v ´3 4 v ´2 4 v ´1 4 v 0 4 v 1 4 v 2

blobs (v ´1

1 ,v 1 1 and v 0 1 respectively.), while filters of higher orders will detect more complex image structures.

C. Edge maps

We stress that it is advantageous to describe an image by patches centered on its edges. We do not focus on solving the edge detection problem, where breakthrough research have been done, such as the recent fast edge detection method based on structured forests [START_REF] Dollár | Structured forests for fast edge detection[END_REF]. We use the existing techniques to dense points from images. We initially consider MSER arbitrarily shaped regions to create an edge map. This 2D map is equal 1 on the region borders and 0 everywhere else. We compute the gradient magnitude on edge pixels and use it as interestingness measure. A local maxima selection is performed to select keypoint locations. Note that only edge pixels are selected in this manner. We center n σ patches at each detected location, one at each of the multiple scales used. We typically use the same number of scales employed in regular dense sampling, as reported in Section III. We will refer to this detection approach as MSER-edge. Comparing MSER to MSER-edge is a way to compare performance based on region description by its interior or by its boundaries. A similar prior attempt is to describe an MSER region by a set of Local Affine Frames extracted with multiple affine-covariant procedures [START_REF] Obdržálek | Object recognition using local affine frames on maximally stable extremal regions[END_REF].

We now consider exactly the same procedure, but applied on the edge map derived from the selective search algorithm of Uijlings's et al. [START_REF] Uijlings | Selective search for object recognition[END_REF]. Their starting point is Felzenszwalb's segmentation algorithm [START_REF] Felzenszwalb | Efficient graph-based image segmentation[END_REF] applied on multiple color spaces.

(a) v ´1 1 (b) v 0 1 (c) v 1 1 (d) v ´2 2 (e) v 2 2 (f) v ´3 3 (g) v ´1 3 (h) v 2 4
Fig. 6. The convolution response for some of the pseudo-Zernike polynomials. The response map is visualized only for the first scale. We use green (blue) patches to visualize the local maxima (minima). Best viewed in color.

Similar segments are merged in a hierarchical manner, and the resulting segments form candidate object regions. We tune it to produce smaller regions, via parameter k, potentially capturing object parts. A super-pixel [START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF] like segmentation is produced in this fashion.

Except for centering patches on selective search region (SSR) borders, we also fit an upright ellipse to each one. We refer to the latter as SSR, and to the former as SSR-edge.

Finally, we consider the fast edge detection algorithm of Dollár and Zitnick [START_REF] Dollár | Structured forests for fast edge detection[END_REF] and apply the same process. The only difference is that the interestingness measure now comes directly from the edge detection process, that is the edge strength. We further eliminate all pixels with strength less than τ , in order to control the number of patches per image. We refer to this approach as fast-edge in our experiments.

D. Dense 2 -norm local maxima

Gosselin et al. [START_REF] Gosselin | Revisiting the Fisher vector for fine-grained classification[END_REF] show that the classification performance increases when SIFT vectors with low energy are filtered out. This reasoning is supported visually, as these features correspond to homogeneous regions with very little visual information. However, finding a threshold value is not a feasible operation as it requires cross validation. The optimal threshold value may vary depending on the image and the size of sampled patches. The effect of this filtering strategy on regular dense is shown in Figure 7, where the descriptors whose squared 2 -norm is lower than the threshold τ are removed.

It appears that SIFT vectors and in particular their 2 -norm provide an interestingness information subsequent at feature detection. In contrast to traditional detectors that localize interest points by low-level statistics, such as gradient changes, we propose to use the 2 -norm of each SIFT as an interestingness measure. Initially, we compute SIFT on a very dense uniform grid of step δ xy " 1 and at 5 different scales. The 2 -norm Detector Section Dense Section III-B Dense Interest Points (Dense-IP) [START_REF] Tuytelaars | Dense interest points[END_REF] Section III-B Dense 2 -norm Section IV-D Harris-Laplace [START_REF] Mikolajczyk | Scale and affine invariant interest point detectors[END_REF], +relaxed Sections III-A, IV-A Frobenius, +relaxed Section IV-A Hessian Affine (HesAff) [START_REF] Mikolajczyk | A comparison of affine region detectors[END_REF] Section III-A Difference of Gaussians (DoG) [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] Section III-A Zernike Section IV-B MSER [START_REF] Matas | Robust wide baseline stereo from maximally stable extremal regions[END_REF], +edge Sections III-A,IV-C Selective Search [START_REF] Uijlings | Selective search for object recognition[END_REF] (SSR), +edge Section IV-C Edge Detector (Fast-edge) [START_REF] Dollár | Structured forests for fast edge detection[END_REF] Section IV-C of SIFT descriptors constitutes a response map. That is the interestingness value per pixel. Following the standard detector procedure, we select the local maxima of this response map. Similar to interest point detectors, we introduce a threshold parameter τ , to only retain the stronger points.

As a consequence, patches from homogeneous regions are discarded and the ones with highest energy of SIFT descriptors and more structured regions are retained. As shown in Figure 2 the retained points are mostly localized along edges and corners. It is the only proposed method that is based on the local descriptor to compute the interestingness measure.

V. EXPERIMENTS

We compare different detection strategies in two different applications, image retrieval and fine-grained image classification. In our experiments, we analyze aspects of feature detection and extraction separately and report our findings. List of detectors used in our experiments are shown in Table I. Note that all the detectors we use are rotation variant (up-right).

A. Experimental Setup

Image Retrieval datasets. INRIA Holidays dataset [START_REF] Jégou | Hamming embedding and weak geometric consistency for large scale image search[END_REF] consists of 1,491 personal holiday photographs. There are 500 image groups, with each group having a single query image. We use the up-right version of this dataset, where each image is rotated according to the natural orientation of objects depicted. Flickr60k dataset [START_REF] Jégou | Hamming embedding and weak geometric consistency for large scale image search[END_REF] is used as a training set for this dataset. Oxford5K Buildings dataset [START_REF]Oxford5k dataset[END_REF] has 5,063 images of Oxford landmarks collected from Flickr. There are 55 query images. Training is performed using Paris 6k dataset [START_REF]Lost in quantization: Improving particular object retrieval in large scale image databases[END_REF].

Fine-grained classification datasets. We use the following evaluation benchmarks:

' FGVC-Aircraft dataset [START_REF] Maji | Finegrained visual classification of aircraft[END_REF] [START_REF] Wah | The Caltech-UCSD Birds-200-2011 Dataset[END_REF]. ' Oxford-IIIT Pet dataset [START_REF] Parkhi | Cats and dogs[END_REF] has 37 categories and 7,349 images of cats and dogs. We randomly select 50% of the training images as validation set [START_REF] Parkhi | Cats and dogs[END_REF].

' The 102 Category Flower dataset [START_REF] Nilsback | Automated flower classification over a large number of classes[END_REF] consists of 102 flower categories, with each category having between 40 and 258 images. We use the provided train, validation and test sets. We did not make any use of the provided objects regions or segments in any of the aforementioned datasets. Although using segmentation information is shown to improve results [START_REF] Wang | Collaborative linear coding for robust image classification[END_REF], we do not consider such an approach. Evaluating the performance on full images offers a direct and fair way to compare different detectors.

Image retrieval pipeline. All images are initially downsampled to 150k pixels. We extract 128-dimensional SIFT descriptors from the detected set of keypoints. We apply RootSIFT [START_REF] Arandjelovic | Three things everyone should know to improve object retrieval[END_REF], [START_REF] Jain | Hamming embedding similarity-based image classification[END_REF], center, rotate along the PCA-axis and 2 normalize the local descriptors. Codebook of 256 visual words is trained with k-means on an independent dataset. We encode the set of local descriptors with the VLAD [START_REF] Jégou | Aggregating local descriptors into a compact image representation[END_REF] representation. Finally, power-law [START_REF] Perronnin | Improving the Fisher kernel for large-scale image classification[END_REF] and 2 normalization is applied to the VLAD vectors. During evaluation, nearest neighbors are found for each query vector, and the performance is measured with mean Average Precision (mAP).

Fine-grained classification pipeline. The down-sampling of images and post-processing of local descriptors is the same as in the retrieval task, except to the fact that we reduce to 80 dimensions using PCA, instead of only rotating and keeping all dimensions. A GMM of 256 components is trained, and the Detector conversion p=11 p=21 p=41 p=61 Dense s " pp ´1q{20 0.5 1 2 3 Dense-IP [START_REF] Tuytelaars | Dense interest points[END_REF] HesAff [START_REF] Mikolajczyk | A comparison of affine region detectors[END_REF] s " pp ´1q{20 0.5 1 2 3 MSER [START_REF] Matas | Robust wide baseline stereo from maximally stable extremal regions[END_REF] SSR [START_REF] Uijlings | Selective search for object recognition[END_REF] Harris-Laplace [START_REF] Mikolajczyk | Scale and affine invariant interest point detectors[END_REF] s " p{2.88 Fisher vector [START_REF] Perronnin | Fisher kernels on visual vocabularies for image categorization[END_REF] is adopted to represent an image. Powerlaw [START_REF] Perronnin | Improving the Fisher kernel for large-scale image classification[END_REF] normalization is applied to the Fisher vectors, with its optimal parameter found through cross-validation. A linear classifier is trained using stochastic gradient descent [START_REF] Akata | Good practice in large-scale learning for image classification[END_REF], and the accuracy is reported.

B. Impact of the scaling factor

Traditionally, a detected region of interest is mapped to a rectangular patch of width equal to p pixels and a local descriptor is extracted from it. It is beneficial to enlarge this measurement region by a scaling factor [START_REF] Simonyan | Learning local feature descriptors using convex optimisation[END_REF], denoted by s.

The ratio between the descriptor measurement region (used to extract local descriptors) and the detection measurement region (used to compute the interestingness measure to detect keypoints) is exactly the scaling factor s. When s " 1, these two are identical.

We investigate the impact of the scaling factor for a variety of detection processes. We linearly relate the scaling factor to the patch size, such that larger regions will be normalized to patches of larger resolution. This exact relation is presented in Table II for all the examined detectors.

Affine regions are isotropically enlarged, and when no enlargement is involved, the patch size is equal to 21. For Zernike, no scaling factor corresponds to descriptor measurement region equal to 11, since this is the filter size we use (see Section IV). For regular dense sampling, where there is no detector measurement region, we have arbitrarily defined the scaling factor to be equal to 1 when the descriptor measurement region size is 21 pixels.

We perform experiments to investigate the impact of the scaling factor and present results in Figures 8 and9. It is observed that such kind of region enlargement is beneficial for all detectors. Performance increases by increasing the descriptor measurement regions, while for large increase it saturates or even drops. Our conclusions agree with the conclusions of Simonyan et al. [START_REF] Simonyan | Learning local feature descriptors using convex optimisation[END_REF], but we evaluate the scaling factor for a larger number of detectors and both classification and retrieval tasks. We adopt the scaling factor that corresponds to a patch size of 41 for all detectors, as it seems to be a good choice overall. The same stands for detectors not included in the corresponding figures. Note that we do not compare different detectors with each other in this experiment. We investigate the impact of scaling factor for each detector individually. [START_REF] Mikolajczyk | Scale and affine invariant interest point detectors[END_REF] 0,50,100,200,300,400 HesAff [START_REF] Mikolajczyk | A comparison of affine region detectors[END_REF] 0,50,100,200,300,400 DoG [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] 0.1,0.2,0.3 Zernike N number of features 1k,5k,10k,16k,20k MSER [START_REF] Matas | Robust wide baseline stereo from maximally stable extremal regions[END_REF] ∆ stable region parameter 3,5,10,15 MSER-edge SSR [START_REF] Uijlings | Selective search for object recognition[END_REF] k initial segment size 25,50,100 SSR-edge Fast-edge [START_REF] Dollár | Structured forests for fast edge detection[END_REF] τ edge strength threshold 0,0.1,0.2,0.3 Each detector produces different number of features per image in this experiment. We consider this a crucial parameter for performance and investigate it in Section V-D.

C. Impact of focusing on edges

In Section IV we proposed detectors that have different localization properties than the existing ones. In particular, we try to focus the local representation on image edges.

In Figure 10, we present the performance measured for Harris-Laplace detector and the relaxed modifications that we proposed. Relaxed-Harris produces more keypoints and consistently improves performance. The relaxation based on Frobenius norm also seems to improve. An exception is the Oxford5K Buildings dataset, where due to particular object matching, original Harris keypoints are localized to have higher repeatability and perform better.

We argue that the typical choice of fitting an ellipse and trying to describe a region by its interior is not necessarily the optimal option for the tasks we consider. Figure 11 shows the results of comparing such a traditional approach to our proposal of extracting local descriptors from patches centered on the region borders. We perform this comparison for MSER regions and selective search regions (SSR). Interestingly, the standard choice does not appear to be the optimal. A larger number of features is obtained when focused on edges, which generally results in better performance. Finally, the region detector based on the fast edge detection algorithm [START_REF] Dollár | Structured forests for fast edge detection[END_REF] (fastedge) performs rather well, especially for intermediate number of descriptors per image.

Additionally, we conduct localization experiments for the proposed Zernike detector by directly controlling the number of keypoints via parameter N z . We also capture more complementary information by increasing the number of filters N f . We evaluate performance using all polynomials up to order 2, 3 and 4. These correspond to 8, 15 and 24 filters respectively. We present the results in Figure 12. Less filters seem to perform better on retrieval datasets. On the contrary, more filters are needed for fine-grained classification, possibly to capture larger intra class variations that exist. We set N f equal to 8 for retrieval and 15 for fine-grained classification. 11. Performance comparison between approaches that describe a detected region by fitting a single ellipse and approaches that sample patches along the region edges. Fast-edge is another detector that samples patches on edges. We show performance versus average number of descriptors per image N .

D. Performance comparison versus number of descriptors

A single parameter for each detector controls the number of keypoints per image. These are shown in Table III. We evaluate performance for multiple values. This allows a fair comparison of all detectors with respect to the number of descriptors used to describe an image. The descriptor set cardinality is directly related to the complexity of encoding stage.

Results of performance comparison between different detectors versus the average number of features per image N are shown in Figures 13 and14, for retrieval and classification respectively. Only the best performing variants of the previous experiments are included in this comparison.

On retrieval it appears that the more descriptors the better, while for fine-grained classification this is not always the case. However, such a behavior seems to be consistent among all detectors for a particular dataset. Thus, it seems to be related with the nature of images and objects depicted in them.

The relaxed Frobenius detector achieves the best performance on Oxford5k, however Hessian-Affine also performs very well with less number of features, due to its better localization properties. However, we observe that the standard dense and relaxed-Frobenius detectors give the best performance using larger number of features on Holidays. This fact is related to the nature of the dataset; it contains more scenes and not particular objects. On fine-grained classification, the Zernike detector outperforms all other detectors by far, with an exception on Oxford-Flowers dataset. The latter is attributed to many points localized on background texture, as shown in the example of Figure 15. Interestingly, over both tasks our dense Finally, note that detectors which extract descriptors from patches of 5 constant scales perform better than the detectors that detect local maxima in the scale space. It appears that the former attains enough tolerance to scale changes, although it is not scale invariant.

E. Effect of 2 -norm based filtering

Gosselin et al. [START_REF] Gosselin | Revisiting the Fisher vector for fine-grained classification[END_REF] show that the accuracy on fine-grained classification based on regular dense sampling significantly increases when the SIFT vectors with low 2 -norm are filtered out. We apply the same strategy and show results on Figure 16. We compute squared 2 -norm to avoid the calculation of square root, and choose T " 5000. As expected, this strategy which aims at removing descriptors extracted from too uniform regions improves performance for regular dense sampling, but the improvement is not significant for detectors that are better localized.

F. State of the art in fine-grained classification

We compare accuracy of Zernike detector, which is best performing among our contributions, with the state of the art results on fine-grained classification. Fisher vectors are employed once more to encode an image based on Zernike patches. This time, we also apply horizontal mirroring to images, in order to obtain another Fisher vector representation per image. As recently shown, encoding spatial information can boost performance [START_REF] Gosselin | Revisiting the Fisher vector for fine-grained classification[END_REF]. We adopt the choice of spatial coordinate coding (SCC) [START_REF] Mccann | Spatially local coding for object recognition[END_REF][26] that allows the use of much larger codebooks, without increase of dimensionality, in contrast to spatial pyramid matching. We choose SCC instead of spatial pyramid matching following the work of Gosselin [START_REF] Gosselin | Revisiting the Fisher vector for fine-grained classification[END_REF], which shows that SCC performs better than spatial pyramid matching for fine-grained image classification. It allows the use of much larger codebooks. In particular, we append regularized 2D coordinates of patches to RootSIFT descriptors before quantizing them. Regularization weights are learned via cross validation. We have set N z " 10k for the Zernike detector. Several works focus on a single domain of fine-grained classification, e.g. birds, and optimize their methods only for that particular domain. Since this is not the case with our study, we compare with methods that are similar. We present our results in Table IV, and compare them with the work of Murray and Perronnin [START_REF] Murray | Generalized max pooling[END_REF]. Using the same codebook, with no additional features, we outperform their results in 2 out of 3 datasets. They also employ color descriptors [START_REF] Clinchant | XRCEs participation to imageval[END_REF] which further improves the accuracy. This can complement also for Zernike patches. We also increase the codebook size and combine with SCC that drastically boosts performance.

To our knowledge, there are not any published results for the FGVC-Aircraft dataset yet, except for the initial paper that achieved accuracy equal to 48.69% [START_REF] Maji | Finegrained visual classification of aircraft[END_REF]. For the other datasets, it is possible to improve the results using the provided annotations of bounding boxes or object parts. We do not consider such an option. However, the best reported result on CUB2011 is 85.40% [START_REF] Branson | Bird species categorization using pose normalized deep convolutional nets[END_REF], which trains convolutional neural networks (CNN) using such annotations. On the other hand, we are superior to methods which use segmentation information on Oxford-IIIT dataset, such as Wang et al. [START_REF] Wang | Collaborative linear coding for robust image classification[END_REF] that report 59.29%. In the case of Oxford-Flowers dataset, the highest reported score is 86.6% [START_REF] Razavian | CNN features off-the-shelf: an astounding baseline for recognition[END_REF], where the authors also use CNNs without the provided segmentation information.

VI. CONCLUSIONS

We investigated dense keypoint detection solutions that lie in between sparse interest points and dense sampling on a uniform grid. We propose to modify existing interest point detectors by relaxing the cornerness criterion and the local maxima selection. We introduce a new detection method using Zernike filters, which provides dense, yet localized image patches. We also show that sampling patches on the borders of a region of interest performs better than the standard choice of fitting an ellipse and describing it by a single descriptor. Finally, we propose to detect dense patches by using the 2norm of the descriptor instead of low-level pixel information. To our knowledge, this is the first detection strategy which focuses on descriptors, and the results seem very promising.

Interestingly, solutions employing patches of multiple fixed scales perform better than patches detected as local maxima in the scale space. Albeit not scale invariant, apparently this option provides enough scale tolerance for the tasks of image retrieval and fine-grained classification.

Compared with the existing studies, Zernike patches encoded with a standard technique, such as Fisher vectors, appear to outperform state of the art approaches for some of the finegrained classification datasets. An exception is the Oxford-Flowers dataset that Zernike seem to perform poorly.

Fig. 1 .

 1 Fig. 1. Locations of patches detected by regular dense sampling (left), by Zernike detector (middle), and by dense 2 -norm based detector (right). We visualize the patches only for the first scale.

Fig. 2 .

 2 Fig. 2. Visualization of regions detected by each method. For methods with predefined scales, we visualize the center of points detected on the first scale (a)-(g), while for the ones that perform scale selection we draw the corresponding ellipses or circles (h)-(l).

Fig. 4 .

 4 Fig. 4. The search neighborhood typically used to detect local maxima is shown on the left. We relax it by taking the union of extrema obtained from multiple neighborhoods shown on the right.

Fig. 3 .

 3 Fig. 3. The detected keypoints for our Harris Laplace modifications. We visualize the centers of detected points.

Fig. 5 .

 5 Fig. 5. Visualization of pseudo-Zernike polynomials for n " 1 . . . 4.

  (a) response map (b) τ " 0.1 ¨10 4 (c) τ " 0.5 ¨10 4 (d) τ " 10 4

Fig. 7 .Fig. 8 .Fig. 9 .

 789 Fig. 7. The response map of 2 -norm of SIFT descriptors is shown in (a). Dark pixels correspond to higher responses. Filtering with different threshold values are shown in (b)-(d). We use the patch size of 41 ˆ41 pixels.

Fig. 10 .

 10 Fig.10. Performance comparison between Harris-Laplace and its relaxed modifications that we propose. We show performance versus average number of descriptors per image N .

  Fig. 11. Performance comparison between approaches that describe a detected region by fitting a single ellipse and approaches that sample patches along the region edges. Fast-edge is another detector that samples patches on edges. We show performance versus average number of descriptors per image N .

Fig. 12 .

 12 Fig. 12. Performance versus average number of descriptors for Zernike detector with different number of filters.

Fig. 13 .

 13 Fig. 13. Performance comparison on image retrieval versus number of descriptors per image N . The parameter values shown in Table III are used.

Fig. 14 .

 14 Fig. 14. Performance comparison on fine-grained classification versus number of descriptors per image N . The parameter values shown in Table III are used.

Fig. 15 .

 15 Fig. 15. Detection example showing the failure of Zernike detector on images of the Oxford-Flowers dataset.

Fig. 16 .

 16 Fig. 16. Impact of 2 -norm based filtering when applied on top of a variety of detectors. We show performance versus average number of descriptors per image N . ˚corresponds to detection with 2 -norm filtering.

TABLE I LIST

 I OF DETECTORS USED IN EXPERIMENTS.

TABLE IV PERFORMANCE

 IV COMPARISON WITH THE STATE OF ART METHODS ON FINE-GRAINED IMAGE CLASSIFICATION. SCC DENOTES SPATIAL COORDINATE CODING AND #C IS THE CODEBOOK SIZE. ALL METHODS USE SIFT DESCRIPTORS, EXCEPT FOR THE LAST ONE WHICH ADDITIONALLY MAKES USE OF COLOR DESCRIPTORS.

-norm local maxima detector achieves best performance for an intermediate amount of points.

ACKNOWLEDGMENT

This work was supported by ERC grant VIAMASS no. 336054 and ANR project FIRE-ID.


Marco Pizzolato

Aurobrata Ghosh

Timothé Boutelier

Rachid Deriche

Perfusion MRI Deconvolution with Delay Estimation and Non-Negativity Constraints

Keywords: Perfusion, Deconvolution, Exponential Bases, Delay, Dispersion, Non-Negative, DSC-MRI

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

Perfusion imaging aims to recover parameters related to the passage of blood in the parenchyma (i.e. the functional part) of a tissue. The amount of perfusion is related to the functionality of the parenchyma and its grade of activity.

The perfusion can be characterized in vivo using Dynamic Susceptibility Contrast MRI (DSC-MRI) where a bolus of paramagnetic contrast agent (PA) is injected in the subject's vascular system. In each voxel, a signal related to the PA concentration is acquired on a time sampling grid to derive the correspondent tissue concentration curve C ts (t). According to the indicator-dilution theory [START_REF] Meier | On the theory of the indicator-dilution method for measurement of blood flow and volume[END_REF], the measured C ts (t) in a voxel is expressed as the convolution between the arterial input concentration to the voxel C a (t), and a curve R(t) which expresses at each time the residual amount of PA in the voxel

C ts (t) = t 0 C a (θ)R(t -θ)dθ (1)
The authors express their thanks to Olea Medical and to the Provence-Alpes-Côte d'Azur Regional Council for providing grant and support.

where R(t) is unknown and must be recovered from the observed C a (t) and C ts (t) by means of deconvolution. The quality of the reconstructed R(t) directly affects the estimation of the perfusion parameters such as the blood flow (BF), the blood volume (BV) and the mean transit time (MTT).

Since deconvolution is ill-posed [START_REF] Keeling | Deconvolution for dce-mri using an exponential approximation basis[END_REF] several regularization methods have been proposed to suppress noise, such as the truncated Singular Value Decomposition (tSVD) [START_REF] Østergaard | High resolution measurement of cerebral blood flow using intravascular tracer bolus passages. part i: Mathematical approach and statistical analysis[END_REF]. More recently a deconvolution method [START_REF] Keeling | Deconvolution for dce-mri using an exponential approximation basis[END_REF] based on exponential bases approximation (EA) showed to outperform tSVD by modeling R(t) as a weighted sum of exponentially decaying functions.

However those techniques show some limitations related to a technical issue inherent the perfusion data acquisition and the physiological nature of the residue function. In real practice, indeed, it is impossible to measure the actual arterial concentration in input to a specific voxel and normally C a (t) is selected as the tissue concentration of a voxel within a region containing only arterial blood. Therefore the concentrations C ts (t) of all the other voxels are physiologically delayed with respect to the selected arterial input C a (t), or can anticipate it, causing the failure of the deconvolution due to the lack of causality assumptions. In addition, the bolus of PA may undergo dispersion which causes the shape of C a (t) to modify, thus changing the natural decreasing shape of R(t) into a non-monotonic one [START_REF] Calamante | Delay and dispersion effects in dynamic susceptibility contrast mri: simulations using singular value decomposition[END_REF]. Such a shape is sometimes considered in myocardial perfusion applications [START_REF] Zarinabad | Voxelwise quantification of myocardial perfusion by cardiac magnetic resonance[END_REF]. Moreover, since R(t) expresses the time-dependent remaining quantity of PA in the voxel, negative values are non-physiological and need to be avoided. The EA does not take into account the time lag (delay) between C a (t) and C ts (t). In addition it does not consider dispersion effects since it cannot guarantee non-negativity unless the residue function is constrained to be monotonically decreasing [START_REF] Keeling | Deconvolution for dce-mri using an exponential approximation basis[END_REF].

We propose Modified Exponential Bases (MEB) allowing time delay estimation and a straight forward implementation of non-negativity constraints also in case of non-monotonic dispersed residue functions. We discuss an iterative linear implementation of the MEB and compare the results with the standard block-circulant formulation of tSVD [START_REF] Wu | Tracer arrival timinginsensitive technique for estimating flow in mr perfusionweighted imaging using singular value decomposition with a blockcirculant deconvolution matrix[END_REF], henceforth addressed as oSVD, which is insensitive to the delay and thus comparable with the proposed technique.

THE DECONVOLUTION PROBLEM

Assuming the tissue concentration C ts (t) and the arterial one C a (t) both measured on an equally spaced time grid t 1 , t 2 , . . . , t M of size M , with ∆t = t i+1 -t i , then Eq. 1 is discretized as

C ts (t j ) = ∆t j i=0 C a (t i )R(t j -t i ) (2)
which can be formulated in matrix form as c ts = Ar, where A is the M ×M convolution matrix containing the samples of the arterial input concentration, c ts contains the M samples of the tissue concentration and r contains the M unknown samples of the residue function. If a model R(t) = G(t, p) for the residue function is defined, then the deconvolution problem is finding the set of parameters p. In case of linearity in the parameters the residue function can be written as r = Gp, where G is the M × N design matrix and p the N × 1 vector of coefficients with N ≤ M . Therefore the convolution is formulated as c ts = AGp, where p is unknown.

In [START_REF] Keeling | Deconvolution for dce-mri using an exponential approximation basis[END_REF] it is proposed the exponential approximation (EA) to perform deconvolution, which searches a solution for R(t) in the span of decaying exponential functions, thus modeling the residue function as a sum truncated to an order N

R(t) = N n=1 k n e -λnt u(t) (3)
where u(t) is the Heaviside step function and the time rates λ n are chosen according to the harmonic distribution λ n = n/T , with T the observation time interval. The unknown parameters p = [k 1 , ..., k N ] are found with R(t) constrained to be at the same time non-negative R(t) ≥ 0 and non-increasing R (t) ≤ 0 for t > 0, where the non-negativity holds due to the non-increasing assumption [START_REF] Keeling | Deconvolution for dce-mri using an exponential approximation basis[END_REF]. Therefore EA, in order to guarantee non-negativity, neglects dispersion effects assuming a non-increasing shape for R(t). In addition it does not take into account the time delay between C a (t) and C ts (t). In the following section we propose a more general solution solving these limitations.

MODIFIED EXPONENTIAL BASIS

In order to estimate the time delay τ and guarantee the non-negativity for both monotonically decreasing and nonmonotonic residue functions, we propose Modified Exponential Bases (MEB). Each basis is constituted as the sum of an exponential decay and the opposite of its derivative with respect to the time rate, respectively weighted by two different constants, a n and b n . Each basis shape can thus vary between a pure exponential decay (b n = 0) and its derivative term (a n = 0). Therefore the MEB include the EA and generalize it by integrating the exponential derivative terms. Moreover the time delay τ is explicitly considered with a Heaviside step function and is included in the exponentials, giving

R(t) = u(t -τ ) N n=1 (a n + b n (t -τ ))e -αn(t-τ ) . (4)
Hence, we can guarantee non-negativity just by constraining all the a n and b n to be non-negative, and at the same time recover both monotonically decreasing and non-monotonic convolution kernels.

The deconvolution problem incorporating Eq. 4 has a nonlinear solution since the unknown α n and τ appear in the exponential. However the time rates α n can be preset without loss of generality [START_REF] Keeling | Deconvolution for dce-mri using an exponential approximation basis[END_REF], thus allowing to find an iterative linear procedure to estimate the unknowns.

In order to fix the rates α n we found that the harmonic rule proposed in [START_REF] Keeling | Deconvolution for dce-mri using an exponential approximation basis[END_REF] is too sensitive to the chosen maximum basis order N , as also reported in [START_REF] Batchelor | Arma regularization of cardiac perfusion modeling[END_REF]. Despite such a rule can be used, we decide to calculate the rates based on the signal rather than the observation time interval. We define M T T max as the maximum M T T expected. This can be fixed empirically, or set according to a conservative value calculated from the MTT found with mono-exponential approximation or via oSVD. Then the time rates α n of Eq. 4 are calculated as the reciprocal of the M T T values harmonically distributed in the range [M T T min , M T T max ] where M T T min = M T T max /N .

In order to consider the delay τ in the estimation we notice that, due to the linearity of the convolution, it can be incorporated in the residue function estimation. To account for this, the problem formulated in Sec. 2 must consider the circular formulation of the convolution [START_REF] Wu | Tracer arrival timinginsensitive technique for estimating flow in mr perfusionweighted imaging using singular value decomposition with a blockcirculant deconvolution matrix[END_REF]. More precisely C ts (t) and C a (t) of length M are extended by zero-padding up to a length L ≥ 2M to avoid aliasing. Then the new L×L squared circular convolution matrix A c has entries

A c i,j = C a (t i-j+1 ) if j ≤ i and A c i,j = C a (t L+i-j+1 ) if j > i.
In the rest we drop the superscript for clarity, A c i,j = A i,j . The convolution problem is then formulated as c ts = AG(τ )p, where G(τ ) is the L × 2N τ -dependent design matrix extended to the circular time sampling grid, and p is the vector containing the 2N coefficients a i and b i to be estimated. We perform the estimation of τ via grid search over a range [τ min , τ max ] seconds, with a time step τ s ≤ ∆t. The set of estimated parameters p is obtained as the Least Squares solution subject to a i , b i ≥ 0, ∀i = 1, ..., N when the estimated delay τ is that minimizing ||c ts -AG(τ )p|| 2 2 among all τ ∈ [τ min , τ max ].

EXPERIMENTAL DESIGN

We test the performance of the compared methods with two synthetic datasets, one generated with a decreasing residue function and the other with a non-monotonic (dispersed) one. In the first case we adopt the bi-exponential model

R bi-exp (t) = A • e -τ1t + (1 -A) • e -τ2t
, where τ 1 = 0.68, τ 2 = 0.05 and A = 0.95, according to the normal tissue values given in [START_REF] Mehndiratta | Modeling the residue function in dscmri simulations: Analytical approximation to in vivo data[END_REF]. In the second case, we use the Pharmacokinetic model as formulated in [START_REF] Batchelor | Arma regularization of cardiac perfusion modeling[END_REF], R pk (t) = (e -λ1t -e -λ2t )/(λ 2 -λ 1 ) -(e -λ1t -e -λ3t )/(λ 3 -λ 1 ) where we fix λ 2 = 0.21, λ 3 = 0.36 according to [START_REF] Zarinabad | Voxelwise quantification of myocardial perfusion by cardiac magnetic resonance[END_REF] and

λ 1 = (λ 3 -λ 2 )/(M T T • λ 2 • λ 3 ) with M T T ≈ 2.
2 s, as obtained from the bi-exponential model. We generate the actual residue function via multiplication of R bi-exp (t) or R pk (t) with BF = 30 ml/100g/min, giving peak values of BF for the bi-exponential model and ≈ 5.7 ml/100g/min for the pharmacokinetic one. In the latter case, the peak value will be considered as BF for the experiments. With the chosen parametrization, BV ≈ 1.1 ml/100g in all cases. The arterial concentration is generated according to [START_REF] Calamante | Delay and dispersion effects in dynamic susceptibility contrast mri: simulations using singular value decomposition[END_REF] as C a (t) = γ 0 (t-t 0 ) ν •e -(t-t0)/β with ν = 3, β = 1.5, γ 0 = 1, and t 0 = 30 s. The tissue concentration C ts (t) is then obtained via Eq. 2. The concentration curves are then converted to signal intensities according to

S(t) = S 0 • e -κ•C(t)•T E ,
where T E is the echo time and κ is the susceptibility constant fixed according to [START_REF] Østergaard | High resolution measurement of cerebral blood flow using intravascular tracer bolus passages. part i: Mathematical approach and statistical analysis[END_REF]. We choose S 0 and T E separately for the tissue concentration S t 0 = 200, T E t = 55 ms and for the arterial one S a 0 = 600, T E a = 13 ms, according to [START_REF] Wirestam | Wavelet-based noise reduction for improved deconvolution of time-series data in dynamic susceptibility-contrast mri[END_REF]. The signals, obtained with repetition ∆t = 1 s, are then corrupted with Gaussian noise with zero mean and standard deviation σ = S 0 /SN R, and finally reconverted in concentration curves.

We performed the experiments separately for C ts (t) obtained after convolution of C a (t) with both the bi-exponential and the pharmacokinetic kernel. In each case we delayed C a (t We notice that, for both MEB and oSVD, the quality of the estimated parameters depends on the chosen convolution kernel (Fig. 1). In the case of a smooth kernel, such as the pharmacokinetic one, the MEB show less effective results than oSVD in estimating the parameters (table. 1(c)). Indeed oSVD cuts the high frequency oscillations of C a (t) which, in such a case, mainly correspond to noise. A possible explanation of the MEB performance is that the delay estimation is harder on the dispersed kernel due to its smoothness. Indeed, we found that the delay estimation error for the pharmacokinetic kernel doubles that obtained for the bi-exponential one (≈ 42% against ≈ 19%, SN R = 80).

RESULTS AND DISCUSSION

On the other hand oSVD show major limitations in the case of bi-exponential residue function. For instance the estimation error of M T T almost doubles the magnitude of the correspondent ground truth value. This performance drop of oSVD is highly related to the marked underestimation of BF , with an error around the 57% of the actual blood flow as shown in Fig. 1(a) and in table 1(b).

The MEB outperform oSVD for the bi-exponential residue function as shown in Fig. 1(a) and in table 1(b). Indeed the MEB manage to reduce the underestimation of BF giving an average relative error in a range [15%, 21%] of the ground truth value. The MEB also allow to obtain an error for M T T always at least three times smaller than that of oSVD. The calculation of BV as the analytical integral of R(t) obtained via MEB is comparable to that obtained via oSVD. However it is normally preferable to estimate BV directly from the concentration curves C a (t) and C ts (t) to avoid introducing errors related to the deconvolution procedure.

The MEB show less residual error then oSVD for both of the tested kernels and at any SN R, as shown in Fig. 2(a) for the dispersed residue function and 2(c) for the bi-exponential one. We further notice that the MEB fitting is more robust to noise. Indeed Fig. 2(a) shows that with dispersed residue function, the difference between the oSVD residual error at SN R = 40 and SN R = 100 is almost four times higher than the correspondent difference for the MEB. The fitting capability of the MEB is more marked for the bi-exponential kernel, Fig. 2(c), where the amount of residual error stands between the 35% and the 60% of the oSVD error for any SN R. In general, the residue function recovered via MEB looks smoother and without ripples, which are a contributing cause to the oSVD residual error.

CONCLUSIONS

We introduced Modified Exponential Bases to perform perfusion deconvolution. The proposed MEB generalize the exponential approximation by including both exponentially decaying terms and their derivatives. Indeed the MEB allow to take into account residue functions which are either purely decreasing or dispersed, providing at the same time a straight forward implementation of non-negativity constraints. The MEB also explicitly take into account the time delay between the arterial and the tissue concentrations. The oSVD showed a lower error in the perfusion parameters estimation for the dispersed residue function, in agreement with the latter's smoothness. However, in the case of the bi-exponential kernel the MEB outperform oSVD. In addition, for any type of kernel, the fitting results of the MEB compare favorably with respect to those obtained via oSVD, and are more robust to noise. We believe that time delay estimation and non-negativity constraints help in recovering a physiologically meaningful residue function. Moreover we find that performing deconvolution with continuous basis functions improves the robustness to noise and we think that proposed bases should take into account the variability in the possible shapes of the residue function.

Fig. 1 .

 1 Fig. 1. Absolute relative error of estimates (mean and standard deviation).

  ) and C ts (t) with integer values in range [-5, 5] seconds. For each delay we calculated the residue function via oSVD (20%) and MEB (N = 30, harmonic α n , M T T max = 4 • M T T oSV D , τ min = -10 s, τ max = 15 s, τ s = 0.25 s) deconvolution, and obtained the estimates of BF , BV and M T T . We tested four SN R values [9], [40, 60, 80, 100], generating 100 noisy repetitions for each possible combination of kernel, delay and SN R.

Fig. 2 .

 2 Fig. 2. Residual error (a,c) and residue function fitting (b,d).

Figure 1 (

 1 Figure 1(a) illustrates the mean absolute relative error in the estimation of each parameter for the representative case of SN R = 80. Figures 1(b) and 1(c) show the tables giving the estimates mean errors and standard deviations for all the SN R, in the case of bi-exponential and pharmacokinetic kernel respectively. For both oSVD and MEB, we computed the mean residual error of the recovered residue function with respect to the ground truth, as shown in Fig. 2(a) and 2(c), for increasing noise levels, in case of dispersed and the bi-exponential kernel respectively; the std areas are shown. Fig. 2(b) and 2(d) show examples of the correspondent MEB and oSVD fittings against the ground truth.