Antoine Tran
email: antoine.[tran@ensta-paristech.fr

Antoine Manzanera
email: manzanera]@ensta-paristech.fr

Fast Growing Hough Forest as a Stable Model for Object Detection

Hough Forest is a framework combining Hough Transform and Random Forest for object detection. The purpose of the present paper is to improve the efficiency and reliability of the original framework by the mean of two contributions. First, instead of generating the image samples by drawing patches randomly from the training set, we bias this step toward the most relevant image content by selecting a proportion of patches from a geometrical criterion. Second, during the creation of non-leafnodes of the trees, instead of sampling uniformly the parameter space for choosing the binary tests aimed at splitting the set of image samples, we choose them according to a probability map constructed from the sample set. We aim to drastically reduce the training time without impacting the accuracy, and decreasing the variability of the produced detectors. The interest of this improved model is shown in the context of car and pedestrian detection by evaluating it on academic datasets.

I. INTRODUCTION

Object detection is a fundamental and difficult problem of computer vision. In our context it consists in localising objects of a certain class in an input image. Typical applications of object detection are video-surveillance and image retrieval. For one given object class (such as cars or humans), the variability of all instances, due to pose, shape, and appearance changes within this class, makes the object representation very challenging. Moreover, even for a single instance, occlusions, illumination and point of view changes also make the object localisation a difficult task.

The principle of the Hough Transform is to consider a set of elements from an image (such as pixels, patches or interest points), and to make them vote for a set of parameters (related to position, scale or orientation). Historically, Hough [START_REF] Hough | Method and means for recognizing complex patterns[END_REF] used contour points to detect lines, circles or ellipses. In object detection context, Leibe and al. [START_REF] Leibe | Interleaving object categorization and segmentation[END_REF] used a keypoint representation in their Implicit Shape Model (ISM) algorithm. During the training step, keypoints from example objects are collected and quantised into a codebook. Every keypoint is also associated to a displacement vector corresponding to its position with respect to the object centroid. During the detection step, all keypoints found in the image are matched with an entry from the codebook, and then vote for all displacements stored in the entry. Finally, a peak detection in the vote collecting map (the Hough transform) is used to detect objects. This detector has proven effective against occlusion and illumination change; it can cope with scale and appearance variations. However, the object representation is sparse, and to be effective, a codebook requires a large number of entries, which implies a time-consuming matching process.

On the other hand, Random Forest classifier [START_REF] Breiman | Random forests[END_REF] is popular in object detection. Moosmann [START_REF] Moosmann | Randomized clustering forests for image classification[END_REF] used it for multi-class detection, using visual descriptor vectors to train the random trees. Then, he used it to generate class-specific saliency maps on unknown images, obtained from sampling patches from the image.

Gall [START_REF] Gall | Class-specific hough forests for object detection[END_REF] proposed a patch-based appearance method, providing a denser level of representation. By using a Random Forest [START_REF] Breiman | Random forests[END_REF] framework associated to a sequence of binary tests related to local geometry, his algorithm is capable to detect different kinds of objects (pedestrian, cars, horse) with a low time-consuming detection step. His training step consists in generating decision trees, which recursively partition an input set of positive and negative patches. Each non-leaf node contains a binary test, while each leaf node contains displacement vectors built according to positive patches positions. During the detection step, each pixel of the image is described by its neighbouring patch, which goes through every decision tree, until it reaches a leaf. The patch then casts votes according to the displacements stored in the leaf. Hough Forest has proven a solid algorithm, giving good results on several academic datasets. It is also a strong basis for diverse enhancements: Ciolini [START_REF] Ciolini | Efficient hough forest object detection for low-power devices[END_REF] proposed different simplifications during the training and detection steps to be able to process the algorithm in embedded systems, while maintaining accuracy performances. Murai [START_REF] Murai | Weighted hough forest for object detection[END_REF] proposed a method to add a weight on patches during training step, reducing the weights of positive patches which are too similar to negative ones, during the training step, in order to improve the accuracy of detection. Gall [START_REF] Gall | An introduction to random forests for multi-class object detection[END_REF] also extended his original Hough Forest algorithm to multi-class and multi-view object detection.

Even though Hough Forest is a strong and versatile algorithm, training a whole forest currently lasts a couple of days following the recommended parameters in a midrange desktop machine, which makes the parametric study and evaluation of the algorithm a complicated task. Simply decreasing some parameters may reduce the training time, however, due to the high level of randomness, the variance of the resulting performance can increase dramatically. In this paper, we propose a one-class object detector based on Hough Forest, providing the following contributions:

• During the sampling step, to generate the positive and negative patch set for the training step, we use a geometric criterion to draw some patches in a deterministic way.

• For the construction of the Hough trees, we propose a way to drastically reduce the number of potential nodes by considering all patches used to train the non-leaf node, and generate probability maps from them.

The main objective of our contribution is to reduce the variance of the object models, while keeping accuracy close to the baseline Hough Forest, and with low time-consuming operations. This paper is organised in three sections: after recalling the baseline Hough Forest algorithm, we present our contributions. Finally, experimentations on academic datasets will be presented and discussed.

II. HOUGH FOREST

Random Forests [START_REF] Breiman | Random forests[END_REF] belong to the random decision tree classification techniques. For classic random decision trees, bagging (generating a subset of training data by randomly sampling with replacement from the full training set) is applied to increase the representativeness of the model. The difference between classic random decision trees and random forest is the use of attribute bagging (randomly choosing the feature from the full feature set).

The Hough Forest is made of binary decision trees, built from a set P of training patches of a constant size W × H, sampled from the learning (positive and negative) images using a uniform random process. P is recursively partitioned during the training phase, and its resulting subsets, corresponding to the leaves of the tree, form the input of the voting process. Given a patch set P, let P + and P -denote the subsets of P coming from the positive and negative images respectively.

Each patch is sampled from a colour image, described with these feature channels:

• 3 components from the Lab colour space • 2 absolute first derivatives:

| ∂L ∂x | and | ∂L ∂y | • 2 absolute second derivatives: | ∂ 2 L ∂x 2 |, | ∂ 2 L ∂y 2 | • 9
bins from the HOG [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF] descriptor. Every channel is additionally spatially enhanced by a min (erosion) and a max (dilation) filter, which finally provides 2 • (3 + 2 + 2 + 9) = 32 features. Each positive patch π is also described with the displacement d π from its centre to the object's centroid.

Each non-leaf node is identified to a binary test T using as input a patch π, and defined from these four parameters: with the following expression:

• 2 pixels (p l , p r ) ∈ π 2 • 1 feature channel c • 1 threshold value τ
T (π) = 0 if c(p l ) > c(p r ) + τ 1 otherwise (1)
A test T will be specified as {p l , p r , c, τ } for the remaining of this paper. The creation of a test (i.e. a non-leaf node) T works as follows. Let P T denote the subset of P, input of the node T . The purpose is to find the test T that realises an optimal partition of P T into two subsets. The training step starts by generating N binary tests {T i } i∈{1,...,N } , with, for each i, p l i , p r i , c i drawn randomly (uniform process), and τ i drawn using a random uniform process in the range formed by the min and the max of {c i (p l i ) -c i (p r i )} i∈{1,...,N } . Each potential binary test T i generates two disjoint subsets

P t Ti = {π ∈ P T |T i (π) = t} with t ∈ {0, 1}. A score
evaluating the quality of the test is then calculated, using either of the two following criteria, chosen randomly:

• One related to the class (positive or negative) repartition of the patches in the two subsets, based on the Shannon entropy H:

U e (T i ) = |P 0 Ti | • H(P 0 Ti ) + |P 1 Ti | • H(P 1 Ti ) (2)
where

H(P) = -|P + | |P| log |P + | |P| -|P -| |P| log |P -|
|P| , and |S| denotes the cardinal of set S.

• One related to the variance of displacements associated to positive patches within the two subsets:

U d (T i ) = t∈{0,1} π∈(P t T i ) + (d π -d (P t T i ) + ) 2 (3)
where d S is the average displacement on set S. Finally, the best test T opt is the one with the lowest score, whatever the chosen criterion.

A node is declared as a leaf if it reaches the maximal tree depth (D max = 15) or if the number of positive input patches is too small (N min = 20). A leaf contains all the displacements from the positive patches, and a weight given by ω L =

|P + L | |P L | , i.
e. the proportion of positive patches in the leaf L.

During the detection step, given an image test I, patches are densely extracted from I. Each patch goes through all trees of the forest, and when it reaches a leaf node, votes for all displacements stored in the leaf, with a vote weight equal to

ω L /|P + L | = 1/|P L |
, such that the cumulative weight for all the votes casted from a leaf L is finally ω L . Finally, the detection output is provided by the set of maxima into the vote map smoothed by a Gaussian filter as in [START_REF] Gall | Class-specific hough forests for object detection[END_REF], or by using Mean-Shift [START_REF] Cheng | Pattern Analysis and Machine Intelligence[END_REF] algorithm, as in [START_REF] Leibe | Robust object detection with interleaved categorization and segmentation[END_REF]. To handle scale variations, Gall also proposed [START_REF] Gall | Class-specific hough forests for object detection[END_REF] to rescale I to different sizes, and then perform the maxima search in the scale-space dimension. III. OUR CONTRIBUTION Hough Forest has proven an effective and versatile framework for object detection. However, considering parameters in [START_REF] Gall | Class-specific hough forests for object detection[END_REF] (forest composed of 15 trees, and 20 000 candidate binary tests to create a non-leaf node), the training step is extremely time consuming. Our goal is to evaluate the performance of the algorithm with those critical parameters significantly lower (9 trees and 500 binary tests), and to propose modifications in the original Hough Forest to reduce variance in these conditions. Our two contributions provide improvements at the following levels:

• During the patch sampling step, instead of using a spatially uniform random process, a proportion of the patches is chosen in a deterministic way, using a multiscale saliency map [START_REF] Lindeberg | Feature detection with automatic scale selection[END_REF]. This guarantees the presence of relevant samples, whatever the content of the training image set. • During the node training step, we improve the efficiency of the patch partition by drawing the binary tests according to a probability map constructed from the current patch collection, instead of choosing them by uniform random selection.

A. Patch sampling

In the baseline method [START_REF] Gall | Class-specific hough forests for object detection[END_REF], the patches are randomly drawn using a uniform probability on the training set. The weakness of this method is the possibility to draw patches in irrelevant (e.g. homogeneous) areas. Such patches are hard to classify using binary tests defined in (Eq. 1), and are sensitive to noise.

Based on this remark, we propose a method to guarantee a minimal number of interesting patches, chosen in a deterministic way. As a saliency measure, we use the multiscale junction-ness function [START_REF] Lindeberg | Feature detection with automatic scale selection[END_REF], defined as follows: for an image I, let I σ x (resp. I σ y ) denote the estimated first derivatives at scale σ, i.e. the image I convolved with the derivative of the Gaussian kernel with variance σ 2 , with respect to x (resp. y). The second derivatives I σ xx , I σ xy , and I σ yy are defined accordingly. The junction-ness map of image I at scale σ is:

J σ I = I σ xx (I σ y ) 2 -2I σ x I σ y I σ xy + I σ yy (I σ x ) 2 (4)
Let N be the number of patches to be selected to train one tree. We set the level of determinism to a rate α ∈ [0, 1], as the proportion of patches to be selected with respect to their saliency, the remaining part being chosen randomly, as in the original method. In order to get structured patches at different sizes, we choose to work with N s different scales. Then, for a given image, and for each scale, starting from the smallest one, we take the N p = α•N Ns pixels with the highest absolute junction-ness as patch centres. To avoid patches spatially too close, we perform non-maximum suppression for each scale, by setting to 0 all pixels at a distance from the selected centre lower than the square of the current scale.

To promote the variety in the different trees of the Hough forest, the non-maximum suppression operations done for one tree remain for the next trees.

The rationale for using the multiscale junction-ness as a saliency measure for choosing the patches is that, as a second order geometric measure, it makes all the involved feature channels significant and it allows to focus on the most discriminant points while representing a negligible time overhead, once the channels calculated. Fig. 2 shows examples of junctions.

B. Node training

In this part, we focus on the non-leaf node training process, which consists in choosing the optimal parameters {p l , p r , c, τ } of the binary test (Eq. 1). As mentioned earlier, Gall's method [START_REF] Gall | Class-specific hough forests for object detection[END_REF] consists in drawing uniformly a high number (Gall mentioned 20 000) of tests, and choosing the best one. However, the number of candidate tests can only be negligible compared to the number of potential tests. Indeed, by using a unique patch size of 16 × 16 pixels, and considering 32 feature channels, there are 16 2 2 × 32 = 2 097 152 potential binary tests, without mentioning the threshold, which is also randomly chosen.

Our method consists in constructing, for each node level, a particular probability law for drawing the binary tests, using the current patch collection as an input. Our goal is to dramatically reduce the number of potential binary tests required to create a node, while improving the quality of the partition with respect to the two criteria (Eq. 2 and 3).

Following the notations of Sec. II, T denotes a binary test (or a non-leaf node), which has to be created, using a patch set P T . Let us consider a deterministic partition of P T into two subsets P T and P T . The binary label denoted by the white or black square will be instanciated later: its nature depends on the criterion to optimise (entropy vs spatial deviation).

The first step is to create global patches (let us call them superpatches), summarising information about P T and P T .

For any patch π let π denote the binarised version -valued in {-1, +1} -of π using Otsu's adaptive threshold method [START_REF] Otsu | A threshold selection method from gray-level histograms[END_REF]. The role of π is to provide contrast invariance, by merely indicating, for a given feature and patch, which pixels belong to the "upper" or "lower" part of the patch. Then, for a ∈ { , }, we define:

Π a = π∈P a T π (5)
For a given label a, Π a can be interpreted as a probability map indicating, for each feature, which pixel is more likely to be a "upper" (resp. "lower") pixel for all patches from P a T . However, directly using as probability density the normalised Π a , whose dynamics can be high, may turn the process too deterministic. We chose instead to drastically quantise the density by using the binary superpatch Π a , resulting from the threshold operation of Π a by 0 and valued in {-1, +1}.

Finally, for a binary test T , and a prior partition P T and P T , we define the two superpatches Π l T and Π r T as follows:

• Π l T = 1 + 1 2 (1 + Π ) + 1 2 (1 -Π ) • Π r T = 1 + 1 2 (1 + Π ) + 1 2 (1 -Π ) Π l
T and Π r T are superpatches valued in {1, 2, 3}, whose purpose is twofold:

• For a given feature channel, normalised Π l T (resp. Π r T ) is used as a probability map to draw the left (resp. right) pixel member p l (resp. p r ) used in Eq. 1.

• To select the most relevant feature channel, the (normalised) sum of the standard deviations of Π l T and Π r T is used as a probability map to draw the feature channel. To conclude with the random selection process of the binary test, the parameters {p l , p r , c} are drawn as explained above, and only the threshold τ is set the same way as in the original Gall's method. Now the missing step consists in determining the two subsets P T and P T , depending on the criterion to minimise:

• When the entropy criterion is chosen, { , } = {-, +}, meaning that we are using directly the labels from negative and positive example patches of P T . • With the spatial deviation criterion, considering the equation (3) to minimise, we apply a k-means (with k = 2) clustering of the patch set P + T with respect to the displacement vector d π . Indeed, the k-means algorithm is obviously a better way to minimise the spatial deviation criterion than Gall's random method. In this case, P T and P T are the two sets resulting of the 2-means of P + T . IV. EXPERIMENTS This section deals with the impact of our contributions by evaluating the improved models in two academic datasets:

• The UIUC-Cars dataset [START_REF] Agarwal | Learning to detect objects in images via a sparse, part-based representation[END_REF]. The training dataset is composed of 550 images of cars, approximatively of the same size (about 100×40 pixels), and 440 negative images. For the detection step, it proposes two datasets: one with 210 single scale cars in 170 images, and one with 139 cars of different scales in 108 images. All cars are taken from side viewpoints. Difficulties are due to partial occlusion, low contrast and complex backgrounds. Images are in grayscale. This dataset is only used for evaluating the impact of the patch selection; we reduce the feature channels to absolute first and second derivatives only. • The TUD-Pedestrians dataset [START_REF] Andriluka | People-trackingby-detection and people-detection-by-tracking[END_REF]. Its training set is composed of 400 images of people walking in different urban backgrounds, with silhouette based ground truth. Following the protocol used in [START_REF] Gall | Class-specific hough forests for object detection[END_REF] to compensate for the low variability of the background (4 different backgrounds), we picked negative images from the INRIA-Person dataset [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF] to train our forest. We also added pedestrians images generated by y-axis symmetry from all positive images of the training set. For the detection step, the test data set is made of 250 images containing 311 pedestrians. Difficulties are due to partial occlusions, scale, pose, aspect and illumination changes. Our objective is to show how our improved models perform with a much smaller number of potential binary tests generated for a node creation (compared to the 20 000 used for the original version). Such drastic reduction of the node parameter space sampling does not affect much the results for the UIUC-Cars dataset (all objects being similar in shape and pose). However for the TUD-Pedestrian dataset, performances can vary a lot from one experiment to the other. In this section, experiments are made following the protocols proposed by the owners of the datasets, and the results are displayed using ROC curves, with error percentage (100 -precision) in abscissa and recall percentage in ordinate.

A. Results on UIUC-Cars

For this dataset, we use a forest of 9 trees, and select 50 patches per image. Three versions of the Hough Forest were evaluated: Fig. 3. ROC-curves for the UIUC-multiscale dataset. In black, the Equal-Error Rate line (EER, Precision = Recall). Green: V1, Blue: V2, Red: V3.

• V1: feature channels reduced to absolute first derivative (as Gall suggested [5]) • V2: feature channels reduced to absolute first and second derivative features • V3: same as previous, plus our deterministic patch selection process (α = 0.1) Fig. 3 shows the resulting ROC-curves for the three versions on the multi-scale test dataset. It can be seen that adding the absolute second order derivative slightly improves the detector. Our model (V3) obtains the best results. However, comparing to results provided by Gall [START_REF] Gall | Class-specific hough forests for object detection[END_REF], particularly EER = 98.5% with 20 000 binary tests, while in our case, EER = 96% with 2 000 tests, and considering that results are very close to 100%, it turns out that this dataset is not relevant enough to evaluate our improved model.

B. Results on TUD-Pedestrians

The TUD-Pedestrians data set is much more challenging than the previous one, as it contains deformable objects (pedestrians) with variable aspects, sizes, and illuminations. Considering our set of parameters and the variability of pedestrian's aspects, performances can vary a lot. Three versions of the Hough Forest were evaluated:

• Va: Original Gall's method • Vb: Improved model with superpatches approach (Part. III-B) but without sampling enhancement (Part. III-A) • Vc: Improved model with both superpatches and sampling enhancement (α = 0.10) For each version, 9 trees per forest were constructed (instead of 15 for [START_REF] Gall | Class-specific hough forests for object detection[END_REF]), 25 000 negative and positive patches, and 500 generated binary tests for each non-leaf node (instead of 20 000 for [START_REF] Gall | Class-specific hough forests for object detection[END_REF]). Junction-ness defined by Eq. 4 is measured at three scales: 1.0, 2.0 and 4.0. For each version, we ran 10 simulations, corresponding to 10 different Hough forests. Fig. 5 and Fig. 6 summarise those simulations.

Even though our results are slightly less accurate than the one obtained by the original Hough Forest with the same parameters, the main interest of our model lies on the stability of the produced detectors in degraded parameter conditions. Indeed, considering standard deviation in Fig. 5, it can be noticed that the red curve, corresponding to the original Hough Forest is slightly the highest one from the three curves, but it is also the one with the highest variance, all along the ROC curve. Additionally, Fig. 6 displays the best and the worst ROC curves obtained within the 10 simulations. While the best results are obtained with the original Hough Forest, this method is also the one with the largest margin between the best and worst results.

If we focus on our methods, Vc provides better results than Vb, but both are more stable than the original Hough Forest Va, and the loss of performance is still very low (the average EER passes from 84.5% in Va to 84.0% in Vc). Finally, Vc results are more stable than Vb, and its accuracy is also closer to Va.

The benefits on the training time, and then on the tractability of the parameter setting, is easy to quantify, as the complexity of the learning phase is linear with respect to the number of sampled tests. Besides, the computational overheads of our two improvements are negligible. Thus, using 500 tests instead of 20 000 typically decreases the training phase of one tree on the TUD-Pedestrians dataset from 8 hours to 15 minutes on a standard desktop PC. With our set of parameters, our training step is slightly slower than the original one. It is mainly due to the superpatches calcutation at each node, requiring to sum every patches from one set. This operation is as slow as the number of patch is high, which is the case for nodes at low depth. Finally, the impact is limited, the time required to train the forest in our case and Gall's case are in the same order of magnitude

In term of memory consumption, we need to store one scalar saliency map for one scale and one image. Each scalar saliency map needs as much memory as one feature channel. Consequently, for one image, we pass from 32 maps for the original algorithm to 35 in ours (3 scales), which is still low.

V. CONCLUSION

In this paper, two contributions to improve the Hough Forest models were presented. Our goal was to drastically reduce the training time, while keeping an accuracy equivalent to the original algorithm, and reducing the variance of the model. Compared to the original Hough Forest, with the same parameters, our algorithm is slightly slower, and need a bit more memory, but is more stable. Our contributions focused on enhancing patch sampling and non-leaf node training, by reducing the randomness to the benefit of the most relevant image content. First, we proposed a method to select some patches based on the multiscale junction-ness measure [START_REF] Lindeberg | Feature detection with automatic scale selection[END_REF]. Second, for a non-leaf node creation, we generated a probability map whose purpose is to summarise the landscapes corresponding to the patch subset used to create this node, and optimise the node quality criterion. These two contributions are simple and computationnally efficient. They allow to reduce the variance of the models in degraded parameters conditions, and then provide a faster and easier way to evaluate the impact of the different parameters on the detection.

Compared to the original parameters [START_REF] Gall | Class-specific hough forests for object detection[END_REF], we dramatically reduced the size of the forest and the number of generated tests used for creating each node. At short term, our objective is to see how our algorithm behaves when reducing another parameter which has a strong impact on the training time: the size of the patch set. We will also work on methods to enhance the way to draw the threshold τ used in Eq. 1. Then, we plan to look for the parameters providing the best results and testing our algorithm on different academic datasets, compared to the state-of-the-art. Finally, as Gall suggested [START_REF] Gall | Hough forests for object detection, tracking, and action recognition[END_REF], we will extend our contribution to different applications, such as object tracking or action recognition.

Fig. 1 .

 1 Fig.1. Vote on one image from TUD-Pedestrian dataset[START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF]. The two peaks correspond to the two pedestrians.

Fig. 2 .

 2 Fig. 2. Junctions shown by circles for one landscape image and two pedestrian images. The radius corresponds to the scale, and the colour to the absolute value of junction-ness (the red corresponding to the highest values).

Fig. 4 .

 4 Fig. 4. Example of false and correct detections in UIUC-Cars and TUD-Pedestrian data sets. In the top-left image, the pedestrian on the left is detected after the false detection, in the middle. In the lop-right one, the scale of the pedestrian detected in the middle is wrong. Below, examples of correct detection.

Fig. 5 .

 5 Fig. 5. Stability of Hough Forest with and without our contributions on the TUD-Pedestrian dataset. The three versions were run 10 times. Each ROC curve corresponds to the average values, while the crosses indicate the standard deviation in precision recall for each threshold level. Va (Gall): Red, Vb (Partially improved): Green, Vc (Fully improved): Blue. Black: EER line.

Fig. 6 .

 6 Fig. 6. Best and worst ROC-curves for each version. Colour code is the same as in Fig. 5.

ACKNOWLEDGMENT

The research was supported by a DGA-MRIS scholarship.


Asma Kerkeni
email: asma.kerkeni@gmail.com

Asma Ben Abdallah
email: asma.babdallah@cristal.rnu.tn

Antoine Manzanera
email: antoine.manzanera@ensta-paristech.fr

Ibtihel Nouira
email: ibtihelnouira@gmail.com

Mohamed Hedi Bedoui
email: medhedi.bedoui@fmm.rnu.tn

Automatic Bifurcation Detection in Coronary X-Ray Angiographies

Keywords: hessian, eigenvector, vesselness, histogram, bifurcationness

come

I. INTRODUCTION

Cardiovascular diseases are the leading cause of deaths all over the world and therefore, researches in modern medical image processing aim at developing reliable and robust medical tools to assist clinicians in data analysis and assistance during diagnostic and interventional procedures. Nowadays, X-ray Coronary Angiography (XCA) is the gold standard for the assessment of cardiovascular diseases (CVD). The angiograms obtained by the XCA enable to reveal the initial CVD symptoms by the morphological features of the coronary arteries such as diameter, length, branching angle, and tortuosity.

The detection of vascular bifurcations in X-ray images is particularly important since these branching points can be used as landmarks for image registration and vessel tracking [START_REF] Feng | Junction-aware shape descriptor for 3d articulated models using local shape-radius variation[END_REF]. Although several methods deal with vascular bifurcation detection [START_REF] Calvo | Automatic detection and characterisation of retinal vessel tree bifurcations and crossovers in eye fundus images[END_REF]- [START_REF] Núñez | Growing algorithm for intersection detection (graid) in branching patterns[END_REF], little work has been devoted to coronary bifurcations in X-ray images. In many blood vessel images such as retinal images, branching point detection generally focuses on binary images. Thus, the original structures are reduced to one-pixel-wide vascular trees through segmentation and thinning. The bifurcations and crossovers are then obtained by counting the pixels local connectivity in the skeleton image [START_REF] Calvo | Automatic detection and characterisation of retinal vessel tree bifurcations and crossovers in eye fundus images[END_REF], [START_REF] Qureshi | A probabilistic model for the optimal configuration of retinal junctions using theoretically proven features[END_REF], [START_REF] Nguyen | Automatic detection of retinal vascular landmark features for colour fundus image matching and patient longitudinal study[END_REF]. Such approaches are usually affected by the thinning artifacts like tailing, necking, spurs or staircase artifact [START_REF] Núñez | Growing algorithm for intersection detection (graid) in branching patterns[END_REF]. Other methods estimate initial bifurcation candidates by the classic Harris corner detector [START_REF] Harris | A combined corner and edge detector[END_REF] based on the first image derivatives. However, in coronary angiographic images, the first-derivative information is usually too noisy and many false detections occur. Hence many detectors rely on the second image derivative, presented by the Hessian matrix, which have strong response on blobs and ridges [START_REF] Mikolajczyk | A comparison of affine region detectors[END_REF]. General Hessian-based features detectors such as difference of Gaussian (DOG) [START_REF] Lowe | Object recognition from local scale-invariant features[END_REF], Hessianaffine [START_REF] Mikolajczyk | A comparison of affine region detectors[END_REF], speeded up robust features (SURF) [START_REF] Bay | Surf: Speeded up robust features[END_REF], smallest Hessian eigenvalue [START_REF] Baboiu | Vascular bifurcation detection in scale-space[END_REF], Determinant of Hessian [START_REF] Baboiu | Vascular bifurcation detection in scale-space[END_REF] and the Laplacian were proposed for blob detection in the computer vision field. The major limit of such detectors is the large rate of false positives since blob points, high curved vessel location, and even some noise artifacts are often confused with branching points.

In [START_REF] Zhao | Bifurcation detection in 3d vascular images using novel features and random forest[END_REF] a novel bifurcation detector named Histogram of Eigenvectors (HOE) is proposed by analogy to the Histogram of Gradient (HOG) feature [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF]. The HOE is combined with other features in order to define a random forest classifier for bifurcation detection in 3D vascular images [START_REF] Zhao | Bifurcation detection in 3d vascular images using novel features and random forest[END_REF]. The essential thought behind the HOE is that bifurcation can be detected thanks to the observation of local vessel orientation provided by the eigenvector corresponding to the smaller Hessian eigenvalues. Thus, for each pixel, an HOE is computed in a local neighbourhood and pixels with three peaks are most likely to be branching points. However, it is unclear how the considered local neighbourhood is selected, though its size is critical since the images often present a wide range of vessel width and a fixed window may lead to underestimate a large bifurcation or combine two nearby vessels in the same window.

In this work, we propose a novel scale specific Histogram of Eigenvector detector, where the local neighbourhood of each pixel is taken as the optimal vesselness scale, yielding a more accurate branching detection and a better bifurcation size estimation. We also propose a novel bifurcationness measure on the basis of this descriptor.

The remaining part of the paper is organized as follows: in section 2 we present the proposed method by describing the different Hessian features and the HOE calculation steps. The experiments and results are presented in section 3. Section 4 summarizes our paper and discusses the possible ways for future works.

II. METHOD

The overview of our proposed vessel bifurcation detector is shown in Fig. 1. First, the eigenvalues and eigenvectors of the Hessian matrix are obtained in multiple scales to provide the structural and directional information. Then, the maximum response with the corresponding optimal scale is retained for each pixel of the image using the Frangi vesselness filter which is built on the Hessian eigenvalues. Finally, for each pixel, the Histogram of Eigenvectors is computed in an appropriate local neighborhood given by the optimal scale of the pixel.

A. Multiscale Hessian Feature extraction

The idea behind multiscale image analysis is to add a new dimension to the analysis which is the image scale. The image is transformed into a set of blurred images, each representing the original image, but at a different scale [START_REF] Lindeberg | Principles for automatic scale selection[END_REF]. These blurred images are obtained by convolving the initial image I 0 (p) = I 0 (x, y) with a Gaussian kernel to represent the information at a certain scale.

I σ (p) = I σ (x, y) = I 0 (x, y) ⊗ G σ (x, y) (1)
where I σ is an image of the scale space, p = (x, y) is a pixel location, ⊗ represents the convolution operation and G σ (x, y) is the 2D Gaussian kernel with standard deviation σ defined as:

G σ (p) = G σ (x, y) = 1 √ 2πσ 2 exp - (x 2 + y 2 ) 2σ 2 (2)
where σ ∈ Σ = {σ min , .., σ max }, σ min and σ max are set according to the approximate width of the smallest and largest vessel to be detected [START_REF] Lindeberg | Principles for automatic scale selection[END_REF]. The 2-by-2 Hessian matrix at a given pixel is composed of the second-order derivatives of an image intensity function Thus, in the scale space, it can be obtained at each point by Eq.3:

H σ (p) = ∂ 2 Iσ(x,y) ∂x 2 ∂ 2 Iσ(x,y) ∂y∂x ∂ 2 Iσ(x,y) ∂x∂y ∂ 2 Iσ(x,y) ∂y 2 . ( 3
)
The Hessian matrix is symmetric, therefore has two real eigenvalues λ 1 and λ 2 and two associated eigenvectors e 1 and e 2 . The signs as well as the relative and absolute magnitudes of eigenvalues characterize the local shape of the intensity in the image. If we assume that we seek to characterize black vessels on a white background, which is the case of the angiographic images, a pixel belonging to a vessel region will be given by small (ideally 0) λ 1 and a high positive value of λ 2 (|λ

1 | < |λ 2 |)

B. Multiscale Vesselness Filter

The filter of Frangi [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF] is often considered as the current gold-standard due to its simplicity, intuitive formulation and good vascular structure enhancement of medical angiographic images [START_REF] Kerkeni | Coronary artery multiscale enhancement methods: A comparative study[END_REF], [START_REF] Kerkeni | A coronary artery segmentation method based on multiscale analysis region growing[END_REF]. The thought behind this filter is that the Hessian eigenvalues are geometrically interpreted as principal vessel curvatures. Generally, vessels have strong variation perpendicular to the vessel centreline (high λ 2 ) and weak variation along it (weak λ 1 ). On the basis of these observations, Frangi et al. introduce two measures to describe structures in images: R b = λ1 λ2 is the bloblike structure measure and S = H F = λ 2 1 + λ 2 2 , the Frobenius norm of the Hessian matrix, is the second-order structureness measure. These measures were combined in a vesselness function as given by Eq.4.

V σ (p) = 0 if λ 2 < 0 exp(- R 2 b 2β 2 )(1 -exp(-S 2 2c 2 )) otherwise. (4)
Where parameters β and c are thresholds which control the filter's sensitivity to R b and S respectively. According to scale space theory, V σ (p) will be maximum up to a normalisation factor when the width of the vessel in pixel p matches a suitable scale factor σ. Thus, In the multi-scale vessel enhancement algorithm, V σ is computed for different scales, then the largest one is taken as the final filter output.

V opt (p) = max σ∈Σ V σ (p) (5)
The optimal scale σ * corresponding to the maximum of vesselness is given by Eq.6.

σ * (p) = arg max σ∈Σ V σ (p) (6)

C. Vessel Direction Information

The eigenvectors e 1 and e 2 of the Hessian matrix characterize the local shape of the intensity. e 1 gives the direction of lower second derivative, i.e. the direction of a potential linear structure, while e 2 gives the direction of the stronger second derivative, i.e. the normal direction to a potential linear structure. Thus the proposed direction information D σ can be simply given by the following expression:

D σ (p) = 0 if λ 2 < 0 e 1 otherwise. (7)
The final vessel direction is given by the direction D σ at the scale σ that provided the optimal response of the filter V opt (p) as given by Eq.8.

D opt (p) = D σ * (p) (8)

D. Histogram of Eigenvectors

Once the vesselness and local orientation features are computed, the Histogram of Eigenvectors (HoE) is obtained through the following steps:

1) For each pixel p, estimate the optimal scale σ * (p) using the equation Eq.6. 2) Construct a circular histogram by casting a vote for each direction vector on the circle centred on p with radius σ * (p) and weight it by the corresponding optimal vesselness (Eq.5). 3) Normalize each bin count on the circle so that the total count adds up to 1; 4) Sort the histogram bins and save the three largest peaks and their corresponding angles. Fig. 3 illustrates circular histograms of a tubular region in (a) and junction location in (b). One can see that the eigenvectors at neighbouring points inside the vessel are oriented in the same direction. Thus, in this region, the histogram of neighbourhood orientation has only one peak. Meanwhile, at the bifurcation, where the vessel splits into two smaller vessels, there are three principal directions and the histogram has three peaks.

E. Bifurcationness measure

According to the HoE feature, a junction candidate has at least three peaks in the histogram. Thus, the proposed bifurcationness measure is set as the height of the third peak in this histogram: it will have a significant value at bifurcations only. The higher the third peak is and the stronger the bifurcationnes will be. Fig. 4 shows the bifurcationness of two different images.

III. EXPERIMENTAL RESULTS

The method has been tested on five real clinical X-ray angiograms acquired during routine cardiac catheter examinations performed in the Cardiology Department of The University Hospital Fattouma Bourguiba, Monastir, Tunisia. The images have different levels of difficulty in terms of vessel illumination and noise. All the results presented in this work use β = 0.5 and half the value of the maximum Hessian norm for c (c = 1 2 max H(x) F ) as proposed in [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF]. Scales were set according to the approximate width of the smallest and largest vessel to be detected. The largest radii (σ max ) was measured and provided by the user whereas σ min was set to 1 since vessels with smaller radii aren't of interest. Four scales between σ min and σ max are generated using a logarithmic progressing.

For the HoE method, only candidates with a bifurcationness measure greater than 0.1 are retained. In terms of absolute time consuming, the Matlab code of the whole procedure on a Core i7 computer of 2.6GHz and 8GB spends an average of 5 seconds per image We have compared the HoE detector with the classic determinant of Hessian detector (DoH). Only pixels with DoH ratio higher than 40 th percentile of the accumulated of the DOH image were retained. Both HoE and DoH candidates are grouped into different connected components using a labelling algorithm assuming the 8-connectivity. Then, the centre of each component is identified as the most suitable point to be centre of the junction. The set of junction candidates is further refined through a non maximal-suppression (NMS) step based on a set of hypothesis related to the nature of images observed on our dataset. For example, the minimum distance between two junction was set to 15 pixels, thus a window of 30-by-30 pixel was used for NMS and the number of points in each connected component was used as the score of the corresponding junction point. The minimum size of a retained component has been set to 5 pixels.The yellow circles in Fig. 5

IV. CONCLUSION

In this paper, we presented a novel method for the detection of junctions in coronary X-ray angiograms. This method is based on the multiscale Hessian analysis. The Hessian eigenvectors, which are a good marker of local vessel orientation have been used to design an Histogram of Eigenvectors detector whereas its eigenvalues were used to weight this histogram and estimate the size of the junction. This leads to an adaptive HoE which behaviour is better than the classic fixed-size version. Based on this detector, the bifurcationness measure is set to the height of the third peak in the histogram. The method has been tested on real coronary artery angiographies and shows better results compared to the classic DoH detector. Overall, the proposed landmark detector approach has potential for application with other 2D and 3D vascular images. In the future work, we would like to refine the bifurcation location and accurately detect its branches. The second image (Image 3) in the table is more difficult than the first due to the low contrast and non-uniform illumination of the background.

Figure 2 .

 2 Figure 2. Hessian-based features: (a) original image; (b) Frangi vesselness response; (c) a small part of the original image; (d) corresponding first eigenvector field superimposed to the vesselness response.

Fig. 2

 2 Fig.2 shows the vesselness and direction features of an X-ray angiogram where Fig.2(b) is the vesselness Map obtained by the Frangi filter and Fig.2(d) is the final eigenvector field of the region of interest shown in Fig.2(c).

Figure 3 .Figure 4 .

 34 Figure 3. Illustration of Histogram of eigenvectors. (a) Eigenvector field of a tubular part and its corresponding circular Histogram. (a) Eigenvector field of a junction location and its corresponding circular Histogram.

  (a) and Fig.5(b) show the final bifurcation points obtained by DoH while the HoE candidates are represented by the red color in Fig.5(c) and Fig.5(d). The quantitative results of the five images for the efficiency of the junction detection process are shown in Table.I. Both DoH and HoE methods can tackle with principal bifurcations however the rate of false positives produced by DoH is very high (10.8 against 4.4 for HoE). The wrong detection in our method is due to background artifacts and noise. The HoE detector also gives less false negatives. The missed junctions are generally located on the crossing between very thin and noisy vessels.

Figure 5 .

 5 Figure 5. Comparison between DoH (yellow) and HoE (red) on two different images. The first image is Image 1 in the table.The second image (Image 3) in the table is more difficult than the first due to the low contrast and non-uniform illumination of the background.

Multiscale Hessian Features extraction Input Image Hessian Eigenvectors Hessian Eigenvalues Multiscale Vesselness Vesselness Map Optimal scales Histogram of Eigenvectors Bifurcationness Map

Figure 1. Overview of the multiscale bifurcation detection method.

Table I COMPARISON

 I BETWEEN DOH AND HOE IN TERMS OF NUMBER OF FALSE POSITIVES AND NUMBER OF FALSE NEGATIVES.


S Graillat
email: stef.graillat@upmc.fr

C Jeangoudoux
email: clothilde.jeangoudoux@lip6.fr

Ch Lauter

MPDI: A Decimal Multiple-Precision Interval Arithmetic Library

Keywords: interval arithmetic, multi-precision, IEEE 754-2008, decimal arithmetic. AMS subject classifications: 65-G30 ://gmplib.org

We describe the mechanisms and implementation of a library that defines a decimal multiple-precision interval arithmetic. The aim of such a library is to provide guaranteed and accurate results in decimal. This library contains correctly rounded (for decimal arbitrary precision), fast and reliable basic operations and some elementary functions. Furthermore the decimal representation is IEEE 754-2008 compatible, and the interval arithmetic is compliant with the new IEEE 1788-2015 Standard for Interval Arithmetic [2].

Introduction

In many situations, binary computer arithmetic is based on a decimal designed system, for the very simple reason that we learn to count in decimal. The conversion between decimal and binary introduces a number of well-known issues: a finite decimal number is not always representable in binary, and the error introduced in this conversion is often amplified along with the combination of operations.

The natural question that comes to mind is then why use decimal instead of binary in the first place. Indeed, the easiest way to avoid such conversion errors is simply to stop designing systems intended to be run on a computer in decimal and do it all in binary. But in many cases this is not possible, due to the difference of cultural 1 background and history of the field, as in the aeronautic industry, or for some financial applications.

Secondly we might want to know why use interval arithmetic. It is very common that applications developed in those fields have critical features. Interval arithmetic consists in representing a number by a certified enclosure instead of an approximated floating-point number. Hence every operation is done with intervals enclosing the real numbers, and computes a result enclosed in an interval too. In the case of critical applications, interval arithmetic can be very useful during the validation and testing process, to ensure a bound on the error of the software.

And thirdly we might want an explanation to why use multiple precision with interval arithmetic. Multiple-precision enables the user to choose the precision of each variable. This precision can be arbitrarily large, and allows us to compute an operation with more digits than in double precision for example. In order to compute the bounds of an interval, we use directed rounding modes. Using multiple-precision arithmetic can in most cases minimize the error made by the computation on those bounds and so to have an interval as narrow as possible.

The application fields of such a library are wide, going from financial calculation to aeronautic engineering. In those fields, sophisticated algorithms are generally designed in high level with decimal numbers. Thus there is a need for correctly rounded decimal arithmetic multiple-precision interval library.

A lot of work has been done in IEEE 754-2008 decimal compliant format [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Std[END_REF]. Our contribution is to propose a correctly rounded, fast and reliable decimal multipleprecision interval decimal library. In the first part we will describe a short overview of the state of the art in decimal scientific libraries (Section 2). Then we will talk about the design of this library (Section 3), and next focus on the operations such as the conversion between decimal and binary (Section 4). Finally, after describing our testing methods and results (Section 5), we will finish by a conclusion and an outlook of the work that can still be done (Section 6).

State of the Art

There are numerous arithmetic libraries available that allows us to perform computation more accurately than the standard floating point arithmetic. Two main types of arithmetic are used to perform more accurate operations: multiple-precision arithmetic and interval arithmetic. Each one has its advantages and drawbacks, and they can very well be combined together as multiple-precision interval arithmetic.

In this section we will shortly describe some of the interesting arithmetic libraries for our problem, and we will develop our motivations for the choices in the implementation.

Short Taxonomy of Arithmetic Libraries

The GNU MP, MPFR and MPFI libraries

When we talk about multiple-precision libraries, a lot of high level Bignum libraries in all kinds of languages like Python or Java come to mind. However, most of them rely on the C GNU Multiple Precision Arithmetic Library, or GMP 1 [START_REF] Granlund | GNU MP: The GNU Multiple Precision Arithmetic Library[END_REF], which provides a fast and reliable arbitrary precision arithmetic on signed integers, rational numbers and floating point numbers. GMP has a rich set of functions, and the functions have a regular interface. This means that each GMP function is defined the same way, as mpz function name(output, inputs). The only practical limit to the precision that can be employed is the memory available in the machine GMP runs on.

The MPFR library [START_REF] Fousse | MPFR: A multiple-precision binary floating-point library with correct rounding[END_REF] is the multiple-precision floating-point extension of GMP and compliant with the IEEE 754 standard. It offers a larger set of functions, and ensures correctly rounded floating-point arithmetic. Each MPFR numbers is represented by a multiple-precision binary significand, an integer exponent, the sign and the expected output binary precision.

The MPFI library [START_REF] Revol | Motivations for an arbitrary precision interval arithmetic and the MPFI library[END_REF] is based on MPFR and provides a multiple-precision floatingpoint arithmetic on intervals. In MPFI, the intervals are represented by two MPFR numbers as the lower and upper bound. The computation of those bounds is made possible by using direct rounding while calling the MPFR functions. The correctness of those bounds is ensured by the correct rounding property of MPFR.

MPFR and MPFI follow the same regular interface and well defined semantics as GMP. They are low level libraries, written in C designed to be fast and reliable.

INTLAB: MATLAB interval library with multiple precision

INTLAB [START_REF] Rump | INTLAB -INTerval LABoratory[END_REF] is a multiple-precision interval toolbox for MATLAB. INTLAB supports real and complex intervals, but also vectors, full and sparse matrices over those.

INTLAB is designed to be very fast, especially on optimized algorithms for matrices with intervals. It also uses the MATLAB interactive programming interface.

The main propose of INTLAB is to produce reliable results fast while keeping the interface with the user at a high level. The user does not have to pay attention to rounding modes or the type of the variables. The interval arithmetic is anabled by the conversion of floating-point numbers with the conversion routine str2intval.

Jinterval for decimal intervals

When we think of fast computing, we mostly think about low level languages. However Java has its advantages, and as shown in [START_REF] Yu | Jinterval library: Principles, development, and perspectives[END_REF] with the Jinterval library there is a way to compute interval arithmetic in an efficient way.

Jinterval performs interval arithmetic on rationals, with the ExtendedRational class. This choice of representation is motivated by the desire to build an algebraic system closed under arithmetic operations. Jinterval is a fast and efficient library, which ensures compliance to the IEEE 1788-2015 Standard for Interval Arithmetic [START_REF]IEEE Standard for Interval Arithmetic. IEEE Std 1788-2015[END_REF].

Jinterval also provides several interval representations, such as classic intervals (defined in Section 3.3), Kaucher intervals [START_REF] Goldsztejn | Modal Intervals Revisited Part 1: A Generalized Interval Natural Extension[END_REF], or complex rectangular or circular intervals. With this flexibility in the choice of interval arithmetic for computations, the user can switch from one arithmetic to another if they are compatible.

Java BigInteger and BigDecimal

The IEEE 754-2008 [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Std[END_REF] revision introduced the definition of decimal floating-point format, with the purpose of providing a basis for a robust and correctly rounded decimal arithmetic. Although there are a lot of libraries implementing decimal arithmetic with words of 32 bits or 64 bits, there are very few in multiple precision arithmetic, and none in interval arithmetic.

Two of those Java classes that provide multiple-precision integer and decimal arith-metic are BigInteger2 and BigDecimal3 . The latter relies on the former, as the BigDecimal type is defined by a BigInterger significand, an integer exponent and a second integer scale, that gives the number of digits to the right of the decimal point. The precision of the computation and the rounding mode are provided by a MathContext constructor. This representation has interesting attributes, but the arithmetic functions implemented in those classes are very limited. They implement the basic operations of addition, multiplication and square root, yet lacking the exponential, logarithm and trigonometric functions. That is why we did not direct our choice of implementation toward those classes.

Choice of GNU MP and MPFR

Studies of implementation of decimal arithmetic libraries have shown that one of the simplest and efficient way to implement decimal functions is to use their binary counterparts [START_REF] Harrison | Decimal transcendentals via binary[END_REF]. This method makes it possible to rely on the efficiency of the binary arithmetic and focus more effort only on some critical points.

Several criteria are taken into account for the choice of implementation of our decimal multiple-precision interval library. They depend on the two goals we want to achieve with the development of this library.

First we want our library to be correctly rounded in decimal and reliable. But to do so we do not want to put too much effort in the implementation of new decimal functions. That is why we choose to rely on the GMP and MPFR libraries to provide binary correctly rounded arithmetic.

Second we want our library to be fast. As we choose to rely on binary functions, we want them to be as fast as possible since we will add some conversion operations. Those operations are not inexpensive, hence the wish for a fast multiple-precision integer and floating-point arithmetic to rely on.

Design of a Decimal Multiple-Precision Interval Library

In this section, we will first set our notations for the rest of the article. Then we will define our decimal format, followed by a short overview of the interval arithmetic chosen for our implementation. To finish, we will talk about the high level design of our library.

Some definitions

In the rest of the article, we will use the following notations to describe the different sets we work with. Let us define:

• R the set of real numbers,

• F 2 p = m • 2 E |E ∈ Z, m ∈ Z, Emin ≤ E ≤ Emax, 2 p-1 ≤ |m| ≤ 2 p -1 ∪ {0} the set of binary representable numbers of precision p, • F 10 k = n • 10 F |F ∈ Z, n ∈ Z, Fmin ≤ F ≤ Fmax, |n| ≤ 10 k -1 the set of deci- mal representable numbers of precision k.
We also choose specific notation for the numbers in those sets, such that:

x, y, z ∈ R are used for real numbers, a, b, c ∈ F 2 p are used for binary numbers, α, β, γ ∈ F 10 k are used for decimal numbers, ε, η, θ will be used to describe errors.

Decimal Number Representation

Up to our knowledge, there is no low level library for decimal multiple-precision arithmetic. So we need to begin our work by creating such a library. To do this we need to create a multiple precision representation of a decimal number.

α := n • 10 F , α ∈ F 10 k . (1)
The decimal number α is defined by:

• n the significand is a GMP signed integer,

• F the exponent is a 64 bit signed integer,

• an additional information about the class of x (NaN, ±∞, ±0 or decimal number)

• k the decimal precision in digits, the library will ensure that 1 ≤ |n| ≤ 10 k -1 at all times.

The IEEE 754-2008 standard [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Std[END_REF] defines a decimal format which guarantees certain properties. For example the use of a quantum guaranties the unique representation of a decimal number.

This property is not provided by our library, hence the same decimal number can have several representations. We seek to have conformity with the decimal format described in IEEE754-2008, not compliance.

Interval Arithmetic

As presented by Moore [START_REF] Moore | Introduction to Interval Analysis[END_REF], an interval can be defined as follows:

[x, x] := {x ∈ R such that x ≤ x ≤ x for some x, x ∈ R, x ≤ x} . (2)
Interval arithmetic provide certified lower and upper bounds on the result of an operation. An operation on this interval is defined as:

[x, x] • [y, y] := x • y such that x ∈ [x, x], y ∈ [y, y] . (3)
In our library, we represent the interval as two multiple precision decimal numbers:

[α, α] := x ∈ R such that α ≤ x ≤ α for some α, α ∈ F 10 k , α ≤ α . (4)
Intervals enable us to represent every number exactly in machine. It is very useful to compute certified algorithm with numbers that are not representable in binary such as 0.1 or even transcendental numbers such as π. The error introduced by the finite binary representation is enclosed.

MPD

High Level Design

In order to implement our decimal multiple-precision interval library, we rely on the GMP library for the definition of our decimal format. Then, assuming a correctly rounded conversion function, we want to use binary functions to implement our library. That is why we use MPFR functions to compute some operations in binary.

Then the MPDI library is built on MPD the same way MPFI relies on MPFR, as we can see from Figure 1. An interval is represented by the two decimal bounds, which are MPD numbers. The computation of those bounds is realized by the computation of the decimal MPD function with directed rounding modes.

Figure 2 illustrates the supported decimal operations, and their dependencies to GMP and MPFR. In each decimal operation is called at least one binary operation. The correctness of the multiple-precision decimal arithmetic relies on a correctly rounded conversion algorithm.

Supported Operations

In this section, we go a little bit further into the description of our library, in particular the basic operations, as the multiplication and addition, and the conversion algorithm. Then we will talk about elementary functions such as the logarithm or the exponential.

Basic Operations

It may seem strange to begin by describing the multiplication algorithm. But the reader will see that with our decimal format, the most simple operation is the multiplication.

Multiplication

The multiplication algorithm is clearly the simplest one, and getting simpler by the use of the decimal representation. It consists of three steps:

α × β = nα • 10 Fα × n β • 10 F β , α, β, γ ∈ F 10 k . (5)
• Compute the result's exponent by adding the two exponents Fα + F β ,

• Compute a temporary significand by the multiplication of the two signifiands nα × n β , with the GMP multiplication function, such as nα × n β = gmp mul(ntemp, nα, n β )

• Round the intermediary result ntemp at the precision k required and adjust the exponent accordingly.

Step three is made possible by our conversion algorithm that also serves the purpose of a rounding algorithm when the size of the number exceeds the precision required by the result.

Addition and Subtraction

The addition and subtraction algorithms are a bit more tricky than the multiplication. Indeed, we cannot follow the basic binary floating-point algorithm, that aligns the numbers according to their exponent by shifting the bits of their significand. In binary the shifting operation is inexpensive, as it comes back to memory manipulation.

To perform the same kind of alignment in decimal, we would need to use decimal division. However this operation is not cheap at all, that is why we must find a better way to do it.

As we describe it in section 4.1.1, the multiplication algorithm is the simplest we can have in decimal. Furthermore in the worst case it only requires one rounding operation. This rounding is performed by our conversion algorithm described below.

Once the numbers are aligned by the multiplication by a power of ten, we can perform the addition or subtraction, without any cancelation. Then we round the result to the precision k required.

Conversion Algorithm

Let us note, with i ∈ N * : i the rounding mode towards + ∞ with i bits of precision, i the rounding mode towards -∞ with i bits of precision, i the rounding mode to the nearest with i bits of precision.

We now describe the guidelines of the conversion algorithm. It is presented as a conversion from binary to decimal, but the same principles are applied for the conversion from decimal to binary.

This algorithm corresponds to the mpd set fr function, or the function that takes a MPFR binary floating point number and converts it into a MPD decimal number.

We want to convert the binary a into the decimal α such that:

a = m • 10 E convert into α = n • 10 F , a ∈ F 2 p , α ∈ F 10 k . (6)
Algorithm 1 Binary to decimal conversion from MPFR to MPD Require: a: the MPFR such that a = m • 10 E k: the decimal precision expected for the output Ensure:

α: the MPD such that α = n • 10 F 1: procedure mpd set fr(a, α)

2:

Compute the exponent F as described in Equation [START_REF] Moore | Introduction to Interval Analysis[END_REF] 3:

Compute the intermediary precision p as described in Equation (8)

4:

Compute r = p (a × p (10 -F )) with MPFR functions if n low = n high then 11:

One of the bound is the signed multiple-precision integer n 12:

else 13:

Increase the precision p and goto line 4

The computation of the exponent can be done with a fixed binary precision because we defined it as an integer. Hence we have

F = t66 -(k -1) with t66 = 66(log 10 (|x|)) ∈ F 2 66 . (7)
All those operations are easily done with MPFR. We use 66 bits of precision because we only need two more bits to ensure that F , which is a 64 bits signed integer, is correct.

The first step to compute the significand is to compute the binary precision p needed in the general case to compute a decimal number of precision k, such that:

p = k • log 2 (10) + c, ( 8
)
where c is a constant allowing p to be a slight overestimation of the binary equivalent of the output precision k, to perform internal calculation. In practice, we choose c = 4.

In this algorithm, we use what is called Ziv's loop [START_REF] Ziv | Fast evaluation of elementary mathematical functions with correctly rounded last bit[END_REF]. It consists of computing a result with more precision than the required precision for the output and decide whether to round this result toward on or the other closest decimal number. To make that decision, the result is represented by an interval and compared with the midpoint between the two numbers. If the interval does not contain the midpoint, then the rounding is easy, else the computation of the interval must be done with a higher precision.

We can prove the termination and the correctness of this algorithm. For the termination, the computation of Equation Line 4 in the algorithm with a greater precision p either provides new information on the computed digits, or is exact, and in that case the computation of the interval is not necessary. For the correctness, we can easily prove that there is a sufficiently large interval to manage the inherent error of the number we try to round. The trick is to take the interval large enough to ensure the correct rounding and the smallest possible to minimize the computation time. We have shown that this corresponds to an interval of three binary shifts in both directions around the number we want to round.

The decimal to binary conversion is essentially the same, except the computation of the result is not separated in exponent and significand. We only compute r = p (α × p (10 F )).

(

) 9
Then we perform the same Ziv's loop as for the binary to decimal conversion, to round the result according to the table maker's dilemma.

Elementary Functions

Some Notations and Discussion about Error Propagation

To go further in building this arithmetic library, we need to take a moment to analyze the propagation of the error. Indeed, the computation of transcendental functions such as the exponential or the logarithm requires to perform several conversions between decimal and binary, and thus introducing an error that will propagate through the algorithm.

There are three kinds of error that can occur during the general case of the computation of a decimal function D : F 10 k → F 10 k approximating a real function f : R → R. First the conversion error ε from a decimal α ∈ F 10 k to a binary a ∈ F 2 p :

•2(α) = a • (1 + ε) with |ε| ≤ 2 -p . (10)
Second the approximation error θ of the binary function B such that:

B : F 2 p → F 2 p , a → f (a)(1 + θ). (11)
If B provides correct rounding, then we have |θ| ≤ 2 -p .

And third the conversion error η from binary to decimal:

•10(a) = α • (1 + η) with |η| ≤ 10 -k . (12)
To sum-up, we can describe the decimal function D as follow:

D(α) = •10(B(•2(α))) = B(a • (1 + ε))(1 + η) D(α) = f (a • (1 + ε))(1 + η)(1 + θ) (13)
with α ∈ F 10 k , a ∈ F 2 p , and ε, θ, η ∈ R.

The errors η and θ are easy to handle. It is the first conversion error ε that gets amplified during the computation of D. Hence we have

D(α) = f (a • (1 + ε)).
We can approximate the function by f (a•(1+ε)) = f (a+aε). Let n be an integer, if f is n times differentiable at the point a ∈ F 2 p ⊂ R, then applying the Taylor-Lagrange formula, there exists a real number ξ between a and a(1 + ε) such that:

f (a + aε) = f (a) + f (a) aε 1! + • • • + f (n) (a) (aε) n n! + f (n+1) (ξ) (n + 1)! (x -a) n+1 = f (a) 1 + af (a) f (a) ε + • • • + a n f (n) (a) n!f (a) ε n + (aε -a) n+1 f (n+1) (ξ) (n + 1)!f (a) . (14)
We can easily neglect the higher order terms and focus on the first. In that case the stability of the error is linked to the condition number of the function. We can then say that if the condition number is bounded, the propagation of the error is not important. However in case it is not, we will have to find some methods to compute the function in another way.

Exponential

Let us study the case of the exponential. We want to bound the conversion error by following the protocol described in Section 4.2.1:

exp(α) = exp(a • (1 + ε)), = exp(a) • exp(aε), = exp(a) • (1 + exp(aε) -1) ε . (15)
We have then the error such as:

ε = exp(aε) -1, 1 + ε = exp(aε), log(1 + ε ) = aε, ε = log(1 + ε ) log(a) . (16)
This equations shows us that we can bound ε such that |ε| ≤ 2 -p for every value α ∈ F 10 k . This implies that using only the conversion, computation and convert back method, we can implement the exponential function in decimal. If the error increases with the exponential, it is always negligible towards the precision wanted for the result.

Logarithm

The case of the logarithm is a bit trickier. To study it, we can write the logarithm as follows:

log(α) = log(a • (1 + ε)) log(α) = log(a) • 1 + 1 log(a) • log(1 + ε) ε • ε (17)
In this equation, we can see that the error ε is related to 1 log(a) . This means that whenever this quantity is unbounded, the error grows tremendously.

The point α = 1 is manually defined by log(α) = 0. But near that point, the error produced by the conversion is amplified by the computation of the logarithm. This means that hence having a correctly rounded result with our conversion algorithm will take exponential time.

But for the logarithm case, we can avoid this problem quite easily by calling the MPFR function mpfr log1p near 1. To do so, we need first to compute β = α -1 in decimal. Then use mpfr log1p on the binary counterpart b of β. The function will compute •2(log(1 + b)) = •2(log(1 + a -1)) = •2(log(a)). In the end we will have the result we are looking for. This solution is not always possible, and in most cases a function does not have only one point of amplification of the conversion error. Like the trigonometric function sin and cos, which require some more work to compute in our decimal multiple-precision interval library.

Experimental Results

In this section we will discuss our experimental approach and results by first talking a bit more about our implementation. In a second part the problem of correctness of the result will be addressed, and finally the timings of our libraries.

Correction Testing

To perform the validation of the correctness of our interval library, we first need to certify the correct rounding of our decimal multiple-precision library. That ultimately boils down to put a decimal number as entry of the test, and verify the correctness of the computed result.

This approach raises two issues: what points in the functions to test, and how to verify the correctness of the multiple-precision decimal result. To do so we need to determine which input we want to test, with which precision, and also the expected output.

The MPFR test frame is very detailed and optimized by bug reports and user experiences. Hence we used it as inspiration for our MPD test frame, in particular the input values used to test specific sensible points of the functions to validate our algorithms.

In MPFR, the results given for the verification of the validation of the tests are decimals converted in binary with the mpfr set str conversion function. To test our MPD library, we cannot put the same output as MPFR for a given entry, as MPFR ensures a binary number of bits correctly rounded, and MPD ensures it for a decimal The tests of MPDI are essentially the same as the MPFI test framework, as it mimics its behavior. The same consideration than with MPD and MPFR tests are taken into account for the process of validation of the correctness of the result.

Performance Testing

The tests are performed on a Intel R Core TM i5-3320M CPU @ 2.60GHz 4, with the operating system Ubuntu 16.04.1 LTS 64 bits. The version of the libraries used for our implementation are GMP 6.1.1 and MPFR 3.1.5. The comparison of the timing results is made with the libraries MPFR 3.1.5 and MPFI 1.5.1.

Implementation Details

Apart from the arithmetical functions described in Section 4, other basic functions have been implemented in our MPD and MPDI libraries. Among them there is the comparison, implemented for both the decimal library and the interval library. These functions also use the conversion algorithm.

So basically, in almost all of our functions, the conversion algorithm is called at least once. This algorithm calls both MPFR and GMP functions, and when the internal precision is not accurate enough to compute the expected precision, the computation is done again with more digits. This implies that we expect the MPD library to more than two times slower than the MPFR library.

The MPDI library behaves as MPFI by computing two times the same algorithm with the bounds of the interval. It is then expected to be as much as twice slower than the MPD library.

Computation Time

To illustrate the computation time considerations of our MPD and MPDI libraries, we choose to compare the result of the exponential and the logarithm algorithms with their binary counterparts in MPFR and MPFI. As expected, and for the reasons highlighted in Section 5.2.1, we can see in Figure 3 that the decimal algorithms are slower than the binary ones. This Figure illustrates the computation time in seconds for the exponential and the logarithm for a set of entries, with a required decimal output precision k = 26 digits, and the equivalent precision p for the binary algorithms. On this Figure, we can see that the MPD functions are in general four times slower than the binary ones, and as expected the MPDI algorithms are two times slower than MPD. We can see peaks in the graphs, that can be explained by the computation of an output number which needed more intermediary precision to achieve the expected precision k = 26.

We tried a different approach to time the algorithms, where we compare the computation time between decimal and binary algorithm in function of the output precision required, as illustrated in Figure 4. The overall conclusion is the same as before, the decimal algorithm is in general four times slower, and the interval one eight times.We separated the timing of the computation of the logarithm with an input close to 1 and afar from 1, but we do not see any difference in the computation time.

However, the increase in the computation time is balanced by the fact that the decimal accuracy of the result is in general higher than with the binary counterparts.

Conclusion and Perspectives

In this paper we have presented MPDI, a decimal multiple-precision ¡interval library written in C. This library is a correctly-rounded arithmetic for decimal arbitrary precision. It contains fast and reliable basic operations and some elementary functions. Furthermore the decimal representation is IEEE 754-2008 compatible, and the interval arithmetic is compliant with the new IEEE 1788-2015 Standard for Interval Arithmetic.

MPDI is build upon a decimal-multiple precision library MPD, that rely on GMP and MPFR. This library is still under development, and new operations such as the trigonometric functions sin and cos will be implemented soon. Another goal is to increase the performance of the different algorithms, by narrowing the use of the conversion algorithm and perform other optimizations.

Figure 2 :

 2 Figure 2: Global architecture of MPD library

Figure 3 :

 3 Figure 3: Timings of exp and log for a precision k = 26 digits

Figure 4 :

 4 Figure 4: Timings of exp and log in function of the precision k for a given value

https://docs.oracle.com/javase/7/docs/api/java/math/BigInteger.html

https://docs.oracle.com/javase/7/docs/api/java/math/BigDecimal.html