Jakub Ciecierski

Bartomiej Dybisz

Agnieszka Jastrzebska

Witold Pedrycz

A Geometrical Approach to Rejecting Option in Pattern Recognition Problem

Keywords: pattern recognition, rejecting option, native and foreign elements

Frequently it happens that during symbols recognition, not all of them are the proper ones. This may cause deterioration of a classifying process. In this paper we present a way to "separate the wheat from the chaff", by constructing a rejector, based on geometrical figures enclosing "wheat" and excluding "chaff". We assume that entities of wheat, called native elements, are structured in some way and that there is no a priori knowledge about chaff, named foreign symbols. For the purpose of this study we present simple geometrical figures to generalize the distribution of symbols and to govern the rejection process.

Introduction

In a standard attempt to pattern recognition an object is classified into one of given classes. However, in practice, this assumption is often too optimistic. Unfortunately, in important practical applications we do not only have elements belonging to one of proper classes, but also elements that do not belong to any of these classes, cf. [START_REF] Bromley | Improving rejection performance on handwritten digits by training with rubbish[END_REF][START_REF] Homenda | Optical Music Recognition: the Case Study of Pattern Recognition[END_REF]. Such elements, which we shall not classify into any of proper classes, could appear for example due to an error in input signal segmentation.

Surprisingly, the issue of rejecting foreign symbols is rarely considered in otherwise very impressive research volume on pattern recognition, despite its obvious importance. We may refer to [START_REF] Bromley | Improving rejection performance on handwritten digits by training with rubbish[END_REF] and [START_REF] Chow | On optimum recognition error and reject tradeoff[END_REF], where the issue of contaminated datasets is recalled. Among few recent papers, dealing with the issue of contaminated datasets we can list [START_REF] Bertolami | Rejection strategies for offline handwritten text line recognition[END_REF] and [START_REF] Elad | Pattern detection using maximal rejection classiffier[END_REF]. The literature study reveals that dominating technique for classification with rejection proposes to issue "flexible" decisions regarding class belongingness. This does not necessarily entail application of fuzzy sets, though in [START_REF] Meel | Disturbance classification and rejection using pattern recognition methods[END_REF] authors point out that this is a viable method. Among studies oriented on rejection in pattern recognition problems we may refer to [START_REF] Koerich | Rejection strategies for handwritten word recognition[END_REF][START_REF] Homenda | Optical Music Recognition: a Case of Pattern recognition with Undesirable and Garbage Symbols[END_REF].

Since elements not belonging to any of proper classes are usually not known a priori, i.e. at the stage of recognizer construction. Moreover, we cannot assume that they form their own class(es) and we cannot use them at this stage. To distinguish between these two types of elements the following terms are used:

-native elements for elements of recognized classes and -foreign elements for ones not belonging to any given class.

Objectives of this study are aimed on inventing and exploring geometrical methods that could be employed for rejecting foreign elements, to test proposed methods on synthetic data and validate them on real data. In addition, the study focuses on comparison of different geometrical figures used for symbol rejection. In particular, two figures are considered: n-dimensional ellipsoid and hyperrectangle. As a working hypothesis we assume that the latter one will perform worse.

Let us emphasise that it is not our goal to elaborate and work on native symbols classification, i.e. determining their class belongingness. Such techniques are already well wrought. Nonetheless, when the need occurs we allow ourselves to use some well-known methods of this area.

The paper is structured as follows. In Section 2 mathematical rudiments of presented work are presented. Section 3 contains information about conducted experiments for both synthetic (3.1) as well real data and with their transformations (3.1). Eventually, Section 4 concludes the work.

Preliminaries

Basic Concepts

Conventionally, a pattern recognition problem is an action of dividing a set of objects S = {s 1 , s 2 . . . , s |S| } into subsets, which include similar objects. Let us assume that S = {S 1 , S 2 , . . . , S |C| } such that (∀i = j)(S i ∩ S j = ∅). Mapping σ : S → C, where C = {1, 2, . . . , |C|}, constitutes a basic model for a task of splitting.

In order to compare and assign class labels to each object we identify them with a vector of some measurable characteristics, which are called features, and then perform a classification. Such approach can be split into two mappings: φ : S → X and ω : X → C, where the first one goes from the space of object to the space of features and the second one from the space of features to the space of classes. It can be easily seen that σ = ω • φ and from now on σ can be referred to as a classifier. In addition, throughout this paper we assume that X

= X 1 × X 2 × . . . × X n = R × R × . . . × R i.e
. space of features is the real coordinate space of n-dimensions -R n and we use lower-case and upper-case letters to distinguish vectors and sets, respectively.

Since set S is commonly not available as a whole, we will construct the rejector only on some part of considered space. Let

S ⊃ L = L 1 ∪ L 2 ∪ . . . ∪ L |C| such that (∀i ∈ 1, |C| )(L i ⊂ S i ).
We call L a learning set. Furthermore we split the learning set into training set(Tr) and testing set(Ts) in the following fashion: L = (T r 1 ∪ T s 1 ) ∪ (T r 2 ∪ T s 2 ) ∪ . . . ∪ (T r |C| ∪ T s |C| ) such that every class of the learning set is split into a training set and a testing set, namely:

(∀i, j ∈ 1, |C| )(T r i ∪ T s i = L i ∧ T r i ∩ T s j = ∅) and Tr = T r 1 ∪ T r 2 ∪ . . . ∪ T r |C| along with Ts = T s 1 ∪ T s 2 ∪ . . . ∪ T s |C| .
It is also worthy to note that the word element is used all over this paper to indicate examined object, as well as its features vector. Usually difference can be distinguished from the context, but wherever it may cause any doubts we explicitly state which is which.

Rejector construction

Since we are interested only in evaluating whether given element belongs or not to any of the proper classes (rather than specifying to which one exactly), each element will be classified either as a native or a foreign one. In this case we assume that the rejector is represented by following mapping: ρ : S → X → C = {native, f oreign}.

Let Z = Z 1 ∪ Z 2 ∪ . . . ∪ Z |C| represent a geometrical region in R n being a union of some figures Z i . This region is assumed to enclose all points in training sets, such that Z i is created based on training set T i , for all i ∈ 1, |C| . For convenience we assume that figures Z i are convex and compact. Such a region will be used to determine if a feature vector from X represents native or foreign element. We will write x ∈ Z to indicate the fact that some x from R n belongs to one or more figures from Z. Since precise information about geometrical objects are defined in Section (2.3) -we are not covering any further details about Z.

Clearly, we are interested in a mapping:

λ : X → C such that (∀x ∈ X) λ(x) = native if x ∈ Z f oreign otherwise
In other words if a vector of features belongs to the region Z, we classify it as native. Otherwise it will be treated as a foreign element. To sum up, we complete ρ with a new mapping λ: ρ : S φ -→ X λ -→ C, where λ is the object of our study.

Geometry

As mentioned in the previous Section, there is still a need to precise, which geometrical figures we will use to construct the rejector, i.e. the region Z. We have chosen n-dimensional ellipsoids and hyperrectangles and in what follows we define how to establish a membership to each of them. Choice of figures has been motivated by their simplicity (hyperrectangles are really straight-forward to compute) and intuitiveness (using ellipsoidal figures rely on the expectation that features vectors will be normally distributed among each class).

Ellipsoids An n-dimensional ellipsoid E can be defined as follows:

E = {x ∈ R n : (x -x 0 ) T A(x -x 0 ) ≤ 1}, (1)
where x 0 ∈ R n is the center of the ellipsoid E and A is a positive definite matrix, cf. [START_REF] Kumar | Minimum-Volume Enclosing Ellipsoids and Core Sets[END_REF][START_REF] Todda | On Khachiyan's algorithm for the computation of minimum-volume enclosing ellipsoids[END_REF]. Having this in mind, it is obvious that any ellipsoid is uniquely represented by its center x 0 and matrix A.

Given volume of the unit n-dimensional hypersphere V 0 , we can compute volume of the ellipsoid E in the following way:

V ol(E) = V 0 det(A) , (2)
Referring to Section 2.2, we define a region Z as union of ellipsoids:

Z = E 1 ∪ E 2 ∪ . . . ∪ E |C|
such that the ellipsoid E i encloses all elements of the training set T r i and has minimal volume, for all classes in C. Let us notice that according to Formula (2) volume minimization is equivalent to maximization of the determinant. The task of finding minimal volume ellipsoids enclosing given sets of elements was solved by implementation of methods given in [START_REF] Todda | On Khachiyan's algorithm for the computation of minimum-volume enclosing ellipsoids[END_REF].

Finally, let us state that for any x ∈ X, x ∈ Z if and only if it is enclosed by at least one ellipsoid from Z.

Hyperrectangles Alike in the case of ellipsoids, we now construct a rejection region Z as union of n-dimensional hyperrectangles enclosing training sets of native elements:

Z = H 1 ∪ H 2 ∪ . . . ∪ H |C|
such that the hyperrectangle H i encloses all elements of the training set T r i , its edges are parallel to axes and it has minimal volume, for all classes in C. More precisely, the hyperrectangle H i enclosing the training set T r i is a Cartesian product of the following finite intervals:

H i = I 1 i × I 2 i × • • • × I n i
where, for all classes i ∈ C, the interval I k i is the minimal interval containing values of k-th parameter for all points of the training set T r i . Assuming that f k (x) is the value of k-th feature of an elements x, we get:

I k i = min f k (x)) : x ∈ T r i , max f k (x) :
x ∈ T r i Again, let us state that for any x ∈ X, x ∈ Z if and only if it is enclosed by at least one hyperrectangle from Z.

Evaluation

For a better understanding of how quality of classification with rejection should be measured we adopt parameters and quality measures used in signal detection theory. Since these parameters are widely utilized, we do not refer to their original sources here. The following parameters were used in defining several measures outlining classification's quality. These parameters create so called confusion matrix, which is given in Table 1. The parameters given in the matrix are numbers of elements of a testing set, which have the following meaning:

-T P -the number of elements of the considered class correctly classified to this class, -F N -the number of elements of the considered class incorrectly classified to other classes, -F P -the number of elements of other classes incorrectly classified to the considered class, -T N -the number of elements of other classes correctly classified to other classes (no matter, if correctly, or not). In this study, original pattern recognition problem is a multiclass problem. However, from the point of view of rejection task, it is important whether a native element is classified to any native class, even to an incorrect native class. Therefore, the question if a native element is classified to the correct or any other native class is not considered here. Hence, a multiclass pattern recognition with rejection is turned to a two-class problem. Finally, the following measures were used assess the quality of the classifier:

Accuracy = TP+TN TP+FN+FP+TN Sensitivity = TP TP+FN Precision = TP TP+FP F-measure = 2 • Precision • Sensitivity Precision + Sensitivity (3)
The following comments help to understand these characteristics:

-Accuracy is a measure of a rejecting performance. This measure describes the ability to distinguish between native and foreign elements. Of course, the higher the value of this measure, the better the identification. -Precision is the ratio of the number of not rejected native elements to the number of all elements (native and foreign) not rejected. Precision evaluates the rejection ability to separate foreign elements from native ones. The higher the value of this measure, the better ability to distinguish foreign elements from native ones. -Sensitivity is the ratio of the number of not rejected native elements to the number of all native ones. This measure evaluates the rejection ability to identify native elements. The higher the value of Sensitivity, the more effective identification of native elements. Unlike the Precision, this measure does not evaluate the effectiveness of separation between native and foreign elements. The higher the value of this measure, the better ability to identify native elements. -Precision and Sensitivity are complementary: increasing sensitivity can cause a drop in precision since, along with increasing the number of not rejected native elements, there might be more not rejected foreign ones. It is there to express the balance between precision and sensitivity since, in practice, these two affect each other. There exists yet another characteristic that combines them: the F-measure. The higher the value of this measure, the better balance between identification of native elements along with separation native elements from foreign ones.

Experiment

Experiments are performed for two kinds of data: synthetic and real, both described in Section 3.1. Elements of both kinds are characterized by 24 features. In both, ellipsoids and hyperrectangles were built on chosen set of elements, cf. Sections 2.2 and 2.3, and effectiveness of rejectors has been checked. Both kinds of datasets have been split into learning and testing set.

In addition, we observe how rejection behaves when we start to remove some of points from training classes and reconstruct geometric figures according to shrunken training sets. We shrink number of elements 4 times; at each stage we subtract 5% of points from the current cardinality of training set.

Figure 1 depicts the process of ellipsoids and hyperrectangles shrinking. In this Figure, elements of two classes are presented in the 2-dimensional space of features in order to keep explanation as lucid as possible.

Datasets

Synthetic data For this kind of data randomizing symbols demands several stages. First of all we are generating 24 intervals (using uniform distribution), where each interval can take values from 0 to 20 inclusively. These intervals serve as a scope of possible values for each of features. Next, we are creating 10 different classes of symbols by assigning to them random vectors of features, as well drawn from a uniform distribution. At this stage we have 10 points in our features space. Next, we distort them by creating clouds of points around each symbol. Process of spreading points is done as follows: take features vector of some symbol and for every value of this vector add a number drawn from N (0, 1). This is done 1500 times per each previously generated symbol. It gives us 10 clouds of points, each having 1500 elements. We treat clouds as learning sets, where 1000 points belong to training set and remaining 500 to testing set.

As for foreign symbols, two approaches were applied: -homogeneous distribution in the space of features, which spreads elements uniformly around the Cartesian product of features intervals. Elements accidentally overlapping with natives were redrawn. -non-homogeneous distribution -created by applying following process: for each segment between a pair of native clouds centers find its middle point.

Treat the points as a centre of new cloud of 1000 foreign points, created using the normal distribution.

In this way, in both manners, we get 10000 foreign points.

Real and semi-synthetic data The datasets were constructed based on the MNIST database of handwritten digits [START_REF] Lecun | The MNIST database of handwritten digits[END_REF]. The following figure [START_REF] Bromley | Improving rejection performance on handwritten digits by training with rubbish[END_REF], represents a sample of natives elements used in the experiments: Then, based on this database, two sets of foreign elements were constructed. Sample of such set is represented by Figure 3. On the other hand, we refer to collection of foreign symbols as 'semi-synthetic', because they are neither pure real, nor pure synthetic objects. They have been created by modifying in some specific fashion (therefore synthetic) the native set of symbols (therefore real). In this paper we incorporate two such sets: A and B. First one comprises of native symbols rotated by 90 degrees, while the second one of uniformly distorted native elements, cf. Figure 3.

Results

This Section is built mostly upon tables, which present results of the experiments. Although we separate synthetic data measurements from real/semisynthetic ones, at the end there is an array concluding both of them.

Synthetic data

First of all we present Tables 2 and3. Each of them is divided into considered geometric figure (upon which rejector has been based on), type of foreign symbols distribution and consecutive steps of shrinking figures numbered as 1000, 950, 903, 858, 816 (number of points in each figure at each stage). The only difference between them is that both training and testing sets are used in the first case, while test set in the second case. One can observe natural behaviour of hyperrectangles, which in the end of the volume minimizing, enclose more native points than ellipsoids. What is more, ellipsoids have better foreign recognition properties, as mostly they are smaller (in terms of volumes) than hypperrectangles. In addition, both figures work better with homogeneous distributed foreign, which is caused by an assumption that feature's vectors of foreign symbols are equally likely to appear within the considered space.

Real and semi-synthetic data As for this type, we provided only one Table [START_REF] Elad | Pattern detection using maximal rejection classiffier[END_REF], which concerns both sets: A and B. 'h-r' is an abbreviation for hyperrectangle, while 'ell' means ellipsoid. In this case each stage of shrinking is marked as '-0%', '-5%', '-10%', '-15%' or '-20%', because our samples contain approximately 1000 features' vectors per symbol class, which causes slight differences in number of points enclosed by each of geometrical figure, hence it can not be generalized. Of course '-0%' symbolizes situation, in which geometric figures are based on all points in training sets. As in the previous Section, data are presented in confusion matrices built with testing set. First remark is that in both foreign sets hyperrectangle approach is not always worse than the ellipsoid one -we have got 100% rejected foreign symbols in case of hyperrectangles in B, while ellipsoids rejected 98.20%. Second observation is that for two sets ellipsoid-based rejectors works well, while hyperrectangles Summary As stated in Section 1, one of the main aims of this study is to measure differences between approach based on ellipsoids and hyperretangles.

In this Section we present Table 5, which is the resume of performed tests, based on measuring sensitivity, accuracy, precision and F-measure features, as described in 2.4. This array helps us observing behaviour of both considered geometrical figures.

We have divided the array into datasets types, foreign types (homogeneous, non-homogeneous; set A, set B), rejector kind which has been used ('h-r' means hyperrectangle; 'ell' describes ellipsoid) and results for consecutive shrinking of geometrical figures.

First, we observe that in ellipsoid approach sensitivity decreases quicker than in hyperrectangle one. Also, interesting is fact that for non-homogeneous foreign symbols, F-measure of hyperrectangle technique actually elevates, while in all others is getting smaller (especially in the ellipsoid type of rejectors). In all cases during shrinking process precision rises, but some of the accuracy is lost.

Conclusion and Future Research Directions

We have proposed an approach to native elements recognition/foreign elements rejection. Native elements belong to proper classes, which we plan to process, while foreign elements may appear in a data set due to an unexpected error. Because classification itself has not been in the scope of this study, presented methods may play crucial role for example in data preprocessing, where one can firstly dispose of foreign elements and then focus on proper classification. Foreign elements removed early from a dataset would not be subjected to a classification process and therefore the proposed technique could improve the overall pattern recognition process. Two methods have been analysed: based on hyperrectangles and on ellipsoids. Despite the fact that ellipsoid-based rejector exhibits high precision in detecting foreign elements, it forfeits significant amount of native elements (especially in shrinking process). On the other hand, hyperrectangular rejector detects native elements on acceptable level but its ability to discern foreign ones is closely related to the learning set layout. As for real life application, it is beneficial to use ellipsoid-based rejector, without shrinking technique, but one must keep in mind that with this method approximately 14% of native points may be wrongly rejected.

As to future directions, we plan to employ models of bipolar uncertainty to features and their aggregation [START_REF] Homenda | Processing of uncertain information in linear space of fuzzy sets[END_REF][START_REF] Homenda | Balanced Fuzzy Sets[END_REF]. Since final decision is crisp (a pattern is either native, or foreign) preferred are univariate bipolar models based on the uncertainty interval [-1,1] rather than bivariate models based on two unit intervals representing negative/positive uncertainties.

Fig. 1 :

 1 Fig. 1: Consecutive shrinking process applied for: ellipsoids (upper part) and hyperrectangles (bottom part).

Fig. 2 :

 2 Fig. 2: Samples of elements of native classes.

Fig. 3 :

 3 Fig. 3: Transformed native elements: rotated by 90 degrees (upper part) and randomly distorted (bottom part).

Table 1 :

 1 Confusion matrix for rejecting in pattern recognition problem

		Classification to the class	Classification to other classes
	The class	True Positives (TP)	False Negatives (FN)
	Other classes	False Positives (FP)	True Negatives (TN)

Table 2 :

 2 Confusion matrices for consecutive shrinking of hyperrectangle (h-r) and ellipsoid (ell) at training set (in %).

	training set			Rejector built for synthetic data sets
	data set		homogeneous		non-homogeneous
	figure	h-r		ell		h-r		ell
	confusion	TP	FN	TP	FN	TP	FN	TP	FN
	matrix	FP	TN	FP	TN	FP	TN	FP	TN
		100.00	0.00	100.00	0.00	100.00	0.00	100.00	0.00
	-0%	0.67 99.33	0.14 99.86	25.39 74.61	14.96 85.04
		98.46	1.54	94.84	5.16	98.46	1.54	94.84	5.16
	-5%	0.48 99.52	0.07 99.93	21.67 78.33	10.54 89.46
		96.98	3.02	90.10	9.90	96.98	3.02	90.10	9.90
	-10%	0.38 99.62	0.03 99.97	17.62 82.38	7.91 92.09
		95.53	4.47	85.64 14.36	95.53	4.47	85.64 14.36
	-15%	0.31 99.96	0.02 99.98	17.13 82.87	5.38 94.62
		94.23	5.77	81.56 18.44	94.23	5.77	81.56 18.44
	-20%	0.29 99.71	0.01 99.99	16.99 83.01	4.83 95.17

Table 3 :

 3 Confusion matrices for consecutive shrinking of hyperrectangle (h-r) and ellipsoid (ell) at testing set (in %).

	testing set			Rejector built for synthetic data sets
	data set		homogeneous		non-homogeneous
	figure	h-r		ell		h-r		ell
	confusion	TP	FN	TP	FN	TP	FN	TP	FN
	matrix	FP	TN	FP	TN	FP	TN	FP	TN
		95.10	4.90	87.96	12.04	95.10	4.90	87.97	12.03
	-0%	0.71	99.29	0.03	99.97	31.03	68.97	8.33	91.67
		93.36	6.64	80.66	19.34	93.36	6.64	80.66	19.34
	-5%	0.47	99.53	0.02	99.98	27.64	72.36	6.31	93.69
		92.12	7.88	76.16	23.84	92.12	7.88	76.16	23.84
	-10%	0.39	99.61	0.00	100.0	24.32	75.68	5.63	94.37
		91.02	8.98	71.68	28.32	91.02	8.98	71.68	28.32
	-15%	0.33	99.67	0.00	100.0	22.55	77.45	4.55	95.45
		89.52	10.48	66.74	33.26	89.52	10.48	66.74	33.26
	-20%	0.26	99.74	0.00	100.0	21.43	78.57	3.84	96.16

Table 4 :

 4 Confusion matrices for consecutive shrinking in ellipsoid (ell) and hyperrectangle (h-r) approach. Considered are sets A and B of foreign elements.

	testing set			Rejector built for real data set
	data set		A				B
	geometry	h-r		ell		h-r		ell
	confusion	TP	FN	TP	FN	TP	FN	TP	FN
	matrix	FP	TN	FP	TN	FP	TN	FP	TN
		98.42	1.58	85.06	14.93	98.45	1.55	85.98	14.02
	-0%	62.41	37.59	12.10	87.90	3.18	96.82	0.11	99.89
		96.45	3.55	80.82	19.18	96.74	3.26	81.83	18.17
	-5%	54.16	45.84	7.53	92.47	2.57	97.43	0.10	99.90
		95.60	4.40	77.25	22.75	95.45	4.55	76.93	23.07
	-10%	49.96	50.04	5.51	94.49	1.95	98.05	0.05	99.95
		93.92	6.08	73.61	26.39	93.88	6.12	73.37	26.63
	-15%	46.41	53.59	4.27	95.73	1.87	98.13	0.02	99.98
		92.77	7.23	69.53	30.47	66.74	33.26	93.15	6.85
	-20%	43.20	56.80	3.18	96.82	0.00	100.0	1.80	98.20
	acquire rather average results in A. Hyperrectangle approach in sets A and B
	manifests difference in points classified as true positive. In the first one, at stage
	'-20%', it is 92.77% while in the second 66.74%. This reveals difference in native
	elements distribution i.e. in set B they are more 'compressed' than in A, hence
	shrinking process has created hyperectangles of a bigger volume in the case of
	set A than in B. In addition in set A, ellipsoid technique, we achieved high
	foreign symbols rejection rate, but at the same time the rate of native elements'
	classification decreased.

Table 5 :

 5 Measuring rejecting results for synthetic and real data for shrunken figures ('-0%', '-5%', '-10%', '-15%' or '-20%'). h-r -hyperrectangle, ell -ellipsoid.

	testing sets		synthetic		real/semi-synthetic
	data	homogeneous	non-homogen.	90ccl		AND
	figure	h-r	ell	h-r	ell	h-r	ell	h-r	ell
	Sensitivity	95.10 87.96	95.10 87.96	98.41 85.06	98.45 85.98
	Accuracy	97.89 95.97	77.68 90.43	57.99 86.95	97.37 95.23
	Precision	98.52 99.93	60.51 84.07	44.32 78.01	93.97 99.75
	F-measure	96.78 93.56	73.96 85.97	61.11 81.38	96.16 92.35
	Sensitivity	93.36 80.66	93.36 80.66	96.45 80.82	96.74 81.83
	Accuracy	97.47 93.54	79.36 89.35	62.82 88.56	97.20 93.85
	Precision	99.00 99.95	62.81 86.47	47.33 84.42	94.99 99.76
	F-measure	96.10 89.27	75.10 83.46	63.50 82.58	95.86 89.90
	Sensitivity	92.12 76.16	92.12 76.16	95.60 77.25	95.45 76.93
	Accuracy	97.11 92.05	81.16 88.30	65.32 88.70	97.18 92.24
	Precision	99.16 100.0	65.44 87.12	49.13 87.62	96.10 99.87
	F-measure	95.51 86.47	76.52 81.27	64.90 82.11	95.78 86.91
	Sensitivity	91.02 71.68	91.02 71.68	99.91 73.61	93.88 73.37
	Accuracy	96.79 90.56	81.97 87.53	67.12 88.31	96.71 91.07
	Precision	99.28 100.0	66.87 88.73	50.53 89.69	96.19 99.94
	F-measure	94.97 83.50	77.10 79.30	65.70 80.86	95.02 84.62
	Sensitivity	89.52 66.74	89.52 66.74	92.77 69.52	93.15 69.53
	Accuracy	96.33 88.91	82.22 86.35	68.86 87.66	96.50 89.79
	Precision	99.42 100.0	67.62 89.68	52.01 91.69	96.30 99.97
	F-measure	94.21 80.05	77.05 76.53	66.65 79.08	94.70 82.02

Acknowledgement

The research is supported by the National Science Center, grant No 2012/07/B/ST6/01501, decision no UMO172012/07/B/ST6/01501.


Jaroslav Pokorný
email: pokorny@ksi.mff.cuni.cz

Graph Databases: Their Power and Limitations

Keywords: graph database, graph storage, graph querying, graph scalability, Big Graphs

Real world data offers a lot of possibilities to be represented as graphs. As a result we obtain undirected or directed graphs, multigraphs and hypergraphs, labelled or weighted graphs and their variants. A development of graph modelling brings also new approaches, e.g., considering constraints. Processing graphs in a database way can be done in many different ways. Some graphs can be represented as JSON or XML structures and processed by their native database tools. More generally, a graph database is specified as any storage system that provides index-free adjacency, i.e. an explicit graph structure. Graph database technology contains some technological features inherent to traditional databases, e.g. ACID properties and availability. Use cases of graph databases like Neo4j, OrientDB, InfiniteGraph, FlockDB, AllegroGraph, and others, document that graph databases are becoming a common means for any connected data. In Big Data era, important questions are connected with scalability for large graphs as well as scaling for read/write operations. For example, scaling graph data by distributing it in a network is much more difficult than scaling simpler data models and is still a work in progress. Still a challenge is pattern matching in graphs providing, in principle, an arbitrarily complex identity function. Mining complete frequent patterns from graph databases is also challenging since supporting operations are computationally costly. In this paper, we discuss recent advances and limitations in these areas as well as future directions.

Introduction

A graph database is any storage system that uses graph structures with nodes and edges, to represent and store data. The most commonly used model of graphs in the context of graph databases is called a (labelled) property graph model [START_REF] Robinson | Graph Databases[END_REF]. The property graph contains connected entities (the nodes) which can hold any number of properties (attributes) expressed as key-value pairs. Nodes and edges can be tagged with labels representing their different roles in application domain. Some approaches refer to the label as the type. Labels may also serve to attach metadata-index or constraint information-to certain nodes.

Relationships provide directed, semantically relevant connections (edges) between two nodes. A relationship always has a direction, a start node, and an end node. Like nodes, relationships can have any properties. Often, relationships have quantitative properties, such as weight, cost, distance, ratings or time interval. Properties make the nodes and edges more descriptive and practical in use. Both nodes and edges are defined by a unique identifier.

As relationships are stored efficiently, two nodes can share any number or relationships of different types without sacrificing performance. Note that although they are directed, relationships can always be navigated regardless of direction. In fact, the property graph model concerns data structure called in graph theory labelled and directed attributed multigraphs.

Sometimes we can meet hypergraphs in graph database software. A hypergraph is a generalization of the concept of a graph, in which the edges are substituted by hyperedges. If a regular edge connects two nodes of a graph, then a hyperedge connects an arbitrary set of nodes.

Considering graphs as a special structured data, an immediate idea which arises is, how to store and process graph data in a database way. For example, we can represent a graph by tables in a relational DBMS (RDBMS) and use sophisticated constructs of SQL or Datalog to express some graph queries. Some graphs can be represented as JSON or XML structures and processed by their native database tools. A more general native solution is offered by graph databases.

One of the more interesting upcoming growth areas is the use of graph databases and graph-based analytics on large, unstructured datasets. A special attention is devoted to so-called Big Graphs, e.g. Facebook with 1 Billion nodes and 140 Billion edges, requiring special storage and processing algorithms [START_REF] Pallavi | Review: Graph Databases[END_REF].

Graph databases are focused on:

 processing highly connected data,  be flexible in usage data models behind graphs used,  exceptional performances for local reads, by traversing the graph.

Graph databases are often included among NoSQL databases 1 . We should also mention lower tools for dealing with graphs. They include frameworks, such as Google's Pregel [START_REF] Malewicz | Pregel: a system for large-scale graph processing[END_REF] -a system for large-scale graph processing on distributed cluster of commodity machines, and its more advanced variant Giraph2 suitable for analytical purposes. They do not use a graph database for storage. These systems are particularly suitable for OLAP and offline graph analytics, i.e. they are optimized for scanning and processing Big Graphs in batch mode. Also the notion of a Big Analytics occurs in this context.

In traditional database terminology, we should distinguish a Graph Database Management Systems (GDBMS) and a graph database. Unfortunately, the latter substitutes often the former in practice. We will also follow this imprecise terminology.

There are a lot of papers about graph models, graph databases, e.g. [START_REF] Larriba-Pey | Introduction to Graph Databases[END_REF], [START_REF] Pallavi | Review: Graph Databases[END_REF], [START_REF] Shimpi | An overview of Graph Databases[END_REF], and theory and practise of graph queries, e.g. [START_REF] Holzschuher | Performance of graph query languages: comparison of cypher, gremlin and native access in Neo4j[END_REF]. Now the most popular book is rather practically oriented work [START_REF] Robinson | Graph Databases[END_REF]. A performance comparison of some graph databases is presented, e.g., in [START_REF] Kolomičenko | Experimental Comparison of Graph Databases[END_REF], [START_REF] Mccoll | A Performance Evaluation of Open Source Graph Databases[END_REF].

In this paper, a lot of examples from the graph database technology will be documented on the most popular graph database Neo4j3 , particularly in its version 2.2. In Section 2 we describe some basic technological features of graph databases. Section 3 presents an overview of graph databases categories as well as some their representatives, i.e., some commercial products. Section 4 presents some facts concerning the paper title and offers some research challenges. Finally, Section 5 concludes the paper.

2

Graph database technology

According to other DBMS, we can distinguish a number of basic components of graph database technology. They include graph storage, graph querying, scalability, and transaction processing. We will discuss them in the following subsections.

Graph storage

An important feature of graph databases is that provide native processing capabilities, at least a property called index-free adjacency, meaning that every node is directly linked to its neighbour node. A database engine that utilizes index-free adjacency is one in which each node maintains direct references to its adjacent nodes; each node, therefore acts as an index of other nearby nodes, which is much cheaper than using global indexes. This is appropriate for local graph queries where we need one index lookup for starting node, and then we will traverse relationships by dereferencing physical pointers directly. In RDBMS we would probably need joining more tables trough foreign keys and, possibly, additional index lookups. Obviously, more advanced indexes are used. For example, it is desirable to retrieve graphs quickly from a large database via graph-based indices, e.g. path-based methods. The approach used in [START_REF] Yan | Graph Indexing: A Frequent Structure -based Approach[END_REF] introduces so called gIndex using frequent substructures as the basic indexing features. Unfortunately, most of these techniques are usable only for small graphs.

Some graph stores offer a graph interface over non-native graph storage, such as a column store in the Virtuoso Universal Server4 in application for RDF data. Often other DBMS is used as back-end storage. For example, the graph database FlockDB5 stores graph data, but it is not optimized for graph-traversal operations. Instead, it is optimized for very large adjacency lists. FlockDB uses MySQL as the basic database storage system just for storing adjacency lists.

2.2

Graph querying

Query capabilities are fundamental for each DBMS. Those used in graph databases, of course, come from the associated graph model [START_REF] Angeles | Survey of Graph Database Models[END_REF]. The simplest type of a query preferably uses the index-free adjacency. A node v k є V is said to be at a k-hop distance from another node v 0 є V, if there exists a shortest path from v 0 to v k comprising of k edges. In practice, the basic queries are the most frequent. They include look for a node, look for the neighbours (1-hop), scan edges in several hops (layers), retrieve an attribute values, etc. Looking for a node based on its properties or through its identifier is called point querying.

Retrieving an edge by id, may not be a constant time operation. For example, Titan 6 will retrieve an adjacent node of the edge to be retrieved and then execute a node query to identify the edge. The former is constant time but the latter is potentially linear in the number of edges incident on the node with the same edge label.

As more complex queries we meet very often subgraph and supergraph queries. They belong to rather traditional queries based on exact matching. Other typical queries include breadth-first/depth-first search, path and shortest path finding, finding cliques or dense subgraphs, finding strong connected components, etc. Algorithms used for such complex queries need often iterative computation. This is not easy, e.g., with the MapReduce (MR) framework used usually in NoSQL databases for BigData processing. But the authors of [START_REF] Rastogi | Finding Connected Components on Map-reduce in Logarithmic Rounds[END_REF] show for finding connected components that some efficient MR algorithms exist. In Big Graphs often approximate matching is needed. Allowing structural relaxation, then we talk about structural similarity queries.

Inspired by the SQL language, graph databases are often equipped by a declarative query language. Today, the most known graph declarative query language is Cypher working with Neo4j database. Cypher commands are loosely based on SQL syntax and are targeted at ad hoc queries of the graph data. A rather procedural graph language is the traversal language Gremlin7 .

The most distinctive output for a graph query is another graph, which is ordinarily a transformation, a selection or a projection of the original graph stored in the database. This implies that graph visualization is strongly tied to the graph querying [START_REF] Pokorny | Big Graph Storage, Processing and Visualization[END_REF].

Scalability

Sharding (or graph partitioning) is crucial to making graphs scale. Scaling graph data by distributing it across multiple machines is much more difficult than scaling the simpler data in other NoSQL databases, but it is possible. The reason is the very nature way the graph data is connected. When distributing a graph, we want to avoid having relationships that span machines as much as possible; this is called the minimum point-cut problem. But what looks like a good distribution one moment may no longer be optimal a few seconds later. Typically, graph partition problems fall under the category of NP-hard problems. Scaling is usually connected with three things:

 scaling for large datasets,  scaling for read performance,  and scaling for write performance.

In practice, the former is most often discussed. Today, it is not problem in graph databases area. For example, Neo4j currently has an arbitrary upper limit on the size of the graph on the order of 10 10 . This is enough to support most of real-world graphs, including a Neo4j deployment that has now more than half of Facebook's social graph in one Neo4j cluster.

Scaling for reads usually presents no problem. For example, Neo4j has historically focused on read performance. In master-slave regime read operations can be done locally on each slave. To improve scalability in highly concurrent workloads, Neo4j uses two levels of caching.

Scaling for writes can be accomplished by scaling vertically, but at some point, for very heavy write loads, it requires the ability to distribute the data across multiple machines. This is the real challenge. For example, Titan is a highly scalable OLTP graph database system optimized for thousands of users concurrently accessing and updating one Big Graph.

Transaction processing

As in any other DBMS, there are three generic use cases for graphs:

 CRUD (create, read, update, delete) applications,  query processing -reporting, data warehousing, and real-time analytics,  batch mode analytics or data discovery.

Graph databases are often optimized and focused on one or more of these uses. Particularly, the first two uses are focused on transactions processing, i.e. OLTP databases. When dealing with many concurrent transactions, the nature of the graph data structure helps spread the transactional overhead across the graph. As the graph grows transactional conflicts typically falls away, i.e. extending the graph tends to the more throughputs. But not all graph databases are fully ACID. However, the variant based on the BASE properties often considered in the context of NoSQL databases is not too appropriate for graphs.

In general, distributed graph processing requires the application of appropriate partitioning and replication strategies such as to maximise the locality of the processing, i.e., minimise the need to ship data between different network nodes.

For example, Neo4j uses master-slave replication, i.e. one machine is designated as the master and the others as slaves. In Neo4j, all writes directed towards any machine are passed through the master, which in turn ships updates to the slaves when polled. If the master fails, the cluster automatically elects a new master.

Neo4j requires a quorum in order to serve write load. It means that a strict majority of the servers in the cluster need to be online in order for the cluster to accept write operations. Othervise, the cluster will degrade into read-only operation until a quorum can be established. Emphasize, that today's graph databases do not have the same level of write throughput as other types of NoSQL databases. This is a consequence of master-slave clustering and proper ACID transactions.

Some more complex architectures occur in the world of graph databases. Typically, a simple database is used to absorb load, and then feed the data into a graph database for refinement and analysis. The architecture Neo4j 2.2 contains even a bulk loader which operates at throughput of million records per second.

Categories of graph databases

There is a lot of graph databases. The well-maintained and structured Web site 8 included 20 products belonging among GDBMSs in 2011. The development of graph databases until 2011 is described in [START_REF] Angeles | A Comparison of Current Graph Database Models[END_REF]. Wikipedia 9 describes 45 such tools. One half of them are ACID compliant. We distinguish general purpose GDBMs, like Neo4j, InfiniteGraph 10 , Sparksee 11 , Titan, GraphBase 12 , and Trinity 13 , and special ones, e.g. the Web graph database In-foGrid 14 and FlockDB, or multimodel databases such as document-oriented databases enabling traversing between documents. For example, OrientDB 15 brings together the power of graphs and the flexibility of documents into one scalable database even with an SQL layer. HyperGraphDB 16 stores not only graphs but also hypergraph structures. All the graph information is stored in the form of key-value pairs. An interesting question is which graph databases are most popular today. In June 2015, the web page DB-Engines Ranking of GDBMS 17 considering 17 graph products presented Neo4j, OrientDB, and Titan on the first three places. GDBMS Sparksee is on the 6th place.

In Section 3.1 we present two typical representatives of the general purpose category. From those special ones, more attention will be devoted to RDF triplestores in Section 3.2. The framework Pregel is explained in Section 3.3.

General graph purpose databases -examples

We describe shortly two successful graph GDBMSs -Neo4j and Sparksee -in some detail. In both GDBMSs a graph is a labelled directed attributed multigraph, where edges can be either directed or undirected.

Example 1: Neo4j Neo4j (now in version 2.2) is the world's leading GDBMS.

It is an open-source, highly scalable, robust (fully ACID compliant) native graph database.

Neo4j stores data as nodes and relationships. Both nodes and relationships can hold properties in a key-value form. Values can be either a primitive or an array of one primitive type. Nodes are often used to represent entities, but depending on the domain the relationships may be used for that purpose as well. The nodes and edges have internal unique identifiers that can be used for the data search. Nodes cannot reference themselves directly [START_REF] Hurwitz | Big Data for Dummies[END_REF]. The semantics can be expressed by adding directed relationships between nodes Graph processing in Neo4j entails mostly random data access which can be unsuitable for Big Graphs. Graphs that cannot fit into main memory may require more disk accesses, which significantly influences graph processing. Big Graphs similarly to other Big Data collections must be partitioned over multiple machines to achieve scalable processing (see Section 2.3).

Example 2: Sparksee

In addition to the basic graph model, Sparksee also introduces the notion of a virtual edge that connects nodes having the same value for a given attribute. These edges are not materialized. A Sparksee graph is stored in a single file; values and identifiers are mapped by mapping functions into B+-trees. Bitmaps are used to store nodes and edges of a certain type.

The architecture of Sparksee includes the core, that manages and queries the graph structures, then an API layer to provide an application programming interface, and the higher layer applications, to extend the core capabilities and to visualize and browse the results. To speed up the different graph queries and other graph operations, Sparksee offers these index types:  attributes,  unique attributes,  edges to index their neighbours, and  indices on neighbours.

Sparksee implements a number of graph algorithms, e.g. shortest path, depth-first searching, finding strong connected components.

Triplestores

Some graph-oriented products are intended for special graph applications, mostly RDF data expressed in the form of subject (S) -predicate (P)object (O). RDF graphs can be viewed as a special kind of a property graph. At the logical level, an RDF graph is then represented as one table. For example, AllegroGraph 18 works with RDF graphs. BrightStarDB 19 , Bigdata 20 and SparkleDB 21 (formerly known as Meronymy) serve for similar purposes. These triple stores employ intelligent data management solutions which combine full text search with graph analytics and logical reasoning to produce deeper results. Sometimes, quad stores are used holding a fourth attribute -the graph name (N) corresponding normally with the namespace of the ontology. AllegroGraph deals even with quints (S, P, O, N, ID), the ID can be used to attach metadata to a triple. Now, GraphDB™ 22 is the world's leading RDF triple store that can perform semantic inferring at scale allowing users to create new semantic facts from existing facts. GraphDB™ is built on OWL (Ontology Web Language). It uses ontologies that allow the repository to automatically reason about the data. AlegroGraph also supports reasoning and ontology modelling.

However, existing triple store technologies are not yet suitable for storing truly large data sets efficiently. According to the W3C Wiki, AllegroGraph leads the largest deployment with loading and querying 1 Trillion triples. The load operation alone took about 338 hours.

We remind also that triple stores create only a subcategory of graph databases. Rather a hybrid solution is represented by Virtuoso Universal Server 23 . Its functionality covers not only processing RDF data, but also relations, XML, text, and others.

A list of requirements often required by customers considering a triple store is introduced in [START_REF]Ontotext: The Truth about Triplestores[END_REF]:

 inferring,  integration with text mining pipelines,  scalability,  extensibility,  enterprise resilience,  data integration and identity resolution,  semantics in the cloud,  semantic expertise.

Pregel and Giraph

Pregel and Giraph are systems for large-scale graph processing. They provide a faulttolerant framework for the execution of graph algorithms in parallel over many machines. Giraph utilizes Apache MR framework implementation to process graphs. A significant approach to the design, analysis and implementation of parallel algorithms, hardware and software in Pregel is now the Bulk Synchronous Processing (BSP) model. BSP offers architecture independence and very high performance of parallel algorithms on top of multiple computers connected by a communication network. BSP is a powerful generalization of MR. A subclass of BSP algorithms can be efficiently implemented in MR [START_REF] Pace | BSP vs MapReduce[END_REF]. BSP is superfast on standard commodity hardware, orders of magnitude faster than the MR alone. It is an easy parallel programming model to learn, it has a cost model that makes it simple to design, analyse, and optimize massively parallel algorithms. It can be considered as a strong candidate to be the programming model for parallel computing and Big Data in the next years. For example, Google is already moving in its internal infrastructure from MR to BSP/Pregel.

Limitations of graphs databases

Despite of the long-term research and practice in this area, there are many important and hard problems that remain open in graph data management. They have influence on functionality restrictions of graph databases (Section 4.1). Others are specifically related to Big Analytics (Section 4.2). Challenges concerning some specific problems of graph database technology are summarized in Section 4.3.

Functionality restrictions

Declarative querying: Most commercial graph databases cannot be queried using a declarative language. Only few vendors offer a declarative query interface. This implies also a lack of query optimization abilities.

Data partitioning: Most graph databases do not include the functionality to partition and distribute data in a computer network. This is essential for supporting horizontal scalability, too. It is difficult to partition a graph in a way that would not result in most queries having to access multiple partitions.

Vectored operations: They support a procedure which sequentially writes data from multiple buffers to a single data stream or reads data from a data stream to multiple buffers. Horizontally scaled NoSQL databases support this type of data access. It seems that it is not the case in graph databases today.

Model restrictions: Possibilities of data schema and constraints definitions are restricted in graph databases. Therefore, data inconsistencies can quickly reduce their usefulness. Often the graph model itself is restricted. Let us recall, e.g., Neo4j nodes cannot reference themselves directly. There might be real world cases where selfreference is required.

Querying restrictions: For example, FlockDB overcomes the difficulty of horizontal scaling the graph by limiting the complexity of graph traversal. In particular, FlockDB does not allow multi-hop graph walks, so it cannot do a full "transitive closure". However, FlockDB enables very fast and scalable processing of 1-hop queries.

Big Analytics requirements

Graph extraction: A question is how to efficiently extract a graph, or a collection of graphs, from non-graph data stores. Most graph analytics systems assume that the graph is provided explicitly. However, in many cases, the graph may have to be constructed by joining and combining information from different resources which are not necessarily graphical. Even if the data is stored in a graph database, often we only need to load a set of subgraphs of that database graph for further analysis.

High cost of some queries: Most real-world graphs are highly dynamic and often generate large volumes of data at a very rapid rate. One challenge here is how to store the historical trace compactly while still enabling efficient execution of point queries and global or neighbourhood-centric analysis tasks. Key differences from temporal DBMSs developed in the past are the scale of data, focus on distributed and inmemory environments, and the need to support global analytics. The last task usually requires loading entire historical snapshots into memory.

Real time processing: As noted, graph data discovery takes place essentially in batch environments, e.g., in Giraph. Some products aimed at data discovery and complex analytics that will operate in real-time. An example is uRIKA24 a Big Data Appliance for Graph Analytics. It uses in-memory technology and multithreaded processor to support non-batch operations on RDF triples.

Graph algorithms: More complex graph algorithms are needed in practice. The ideal graph database should understand analytic queries that go beyond k-hop queries for small k. Authors of [START_REF] Mccoll | A Performance Evaluation of Open Source Graph Databases[END_REF] did a performance comparison of 12 open source graph databases using four fundamental graph algorithms (e.g. simple source shortest path problem and Page Rank) on networks containing up to 256 million edges. Surprisingly, the most popular graph databases have reached the worst results in these tests. Current graph databases (like relational databases) tend to prioritize low latency query execution over high-throughput data analytics.

Parallelisation: In the context of Big Graphs there is a need for parallelisation of graph data processing algorithms when the data is too big to handle on one server. There is a need to understand the performance impact on graph data processing algorithms when the data does not all fit into the memory available and to design algorithms explicitly for these scenarios.

Heterogeneous and uncertain graph data:

There is a need to find automated methods of handling the heterogeneity, incompleteness and inconsistency between different Big Graph data sets that need to be semantically integrated in order to be effectively queried or analysed.

4.3

Other challenges

Other challenges in the development of graph databases include:

Design of graph databases: Similarly to traditional databases, some attempts to develop design models and tools occur in last time. In [START_REF] De Virgilio | Model-Driven Design of Graph Databases[END_REF], the authors propose a modeldriven, system-independent methodology for the design of graph databases starting from ER-model conceptual schema.

Need for a benchmark: Querying graph data can significantly depend on graph properties. The benchmarks built, e.g., for RDF data are mostly focused on scaling and not on querying. Also benchmarks covering a variety of graph analysis tasks would help towards evaluating and comparing the expressive power and the performance of different graph databases and frameworks.

Developing heuristics for some hard graph problems: For example, partitioning of large-scale dynamic graph data for efficient distributed processing belongs among these problems, given that the classical graph partitioning problem is NP-hard.

Graph pattern matching: New semantics and algorithms for graph pattern matching over distributed graphs are in development, given that the classical subgraph isomorphism problem is NP-complete.

Compressing graphs: Compressing graphs for matching without decompression is possible. Combining parallelism with compressing or partitioning is also very interesting.

Integration of graph data: In the context of Big Data, query formulation and evaluation techniques to assist users querying heterogeneous graph data are needed.

Visualization: Improvement of human-data interaction is fundamental, particularly a visualization of large-scale graph data, and of query and analysis results.

Graph streams processing: Developing algorithms for processing Big Graph data streams with goal to compute properties of a graph without storing the entire graph.

Conclusions

Graph databases are becoming mainstream. As data becomes connected in a more complicated way and as the technology of graph databases matures, their use will increase. New application areas occur, e.g. the Internet of Things, or rather Internet of Connected Things. Comparing to traditional RDBMS, there is a difficulty for potential users to identify the particular types of use case for which each product is most suitable. Performance varies greatly across different GDBMSs depending upon the size of the graph and how well-optimized a given tool is for a particular task. It seems that especially for Big Graphs and Big Analytics a lot of previous results and designs will have to be re-considered and re-thought in next research and development.

  20 http://www.systap.com/ (retrieved on 9.3.2015) 21 https://www.sparkledb.net/ (retrieved on 9.3.2015) 22 http://www.ontotext.com/products/ontotext-graphdb/ (retrieved on 9.3.2015) 23 http://virtuoso.openlinksw.com/ (retrieved on 9.3.2015)

http://nosql-database.org/

http://giraph.apache.org/

http://www.neo4j.org/ (retrieved on9.3.2015)

http://virtuoso.openlinksw.com/ (retrieved on9.3.2015)

https://github.com/twitter/flockdb (retrieved on9.3.2015)

http://thinkaurelius.github.io/titan/ (retrieved on9.3.2015)

http://gremlindocs.com/

http://franz.com/agraph/ (retrieved on9.3.2015)

http://brightstardb.com/ (retrieved on9.3.2015)

http://www.cray.com/products/analytics/urika-gd

Acknowledgments

This paper was supported by the GAČR grant No. P103/13/08195S.


Andrea Mendizabal

Rémi Bessard Duparc

Huu Phuoc Bui

Christoph J Paulus

Igor Peterlik

Stéphane Cotin

Face-based Smoothed Finite Element Method for Real-time Simulation of soft tissue

Keywords: smoothed finite element method, strain smoothing, real-time simulation, soft tissue modeling 1

come   L'archive ouverte pluridisciplinaire

In medical simulations, real-time performance is a strong requirement for clinical use. The FS-FEM was proposed by Liu et al 1 in order to improve the accuracy of the FE-based modeling. Recent works have adapted the FS-FEM to biological research. [START_REF] Duong | A face-based smoothed finite element method for hyperelastic models and tissue growth[END_REF] However, to our best knowledge, the method has not been studied or adapted in the context of real-time simulations and no implementation of the S-FEM approach allowing for real-time computation exists.

In this paper, we employ the FS-FEM approach in the context of navigation during the brain tumor surgery [START_REF] Clatz | Robust nonrigid registration to capture brain shift from intraoperative mri[END_REF][START_REF] Miga | In vivo quantification of a homogeneous brain deformation model for updating preoperative images during surgery[END_REF][START_REF] Bucki | Framework for a low-cost intra-operative image-guided neuronavigator including brain shift compensation[END_REF] and partial nephrectomy. During the intervention, the position of the tumor has to be accurately estimated from the pre-operative CT or MR data. In the case of brain surgery, a burr hole is drilled in the skull in order to insert surgical instruments and the cerebrospinal fluid (CSF) leaks resulting in the motion and deformation of the brain known as the brain shift. As it is difficult to predict the exact amount of the brain shift before the leak of the CSF, the pre-operative planning of the intervention becomes inaccurate. [START_REF] Nabavi | Serial intraoperative magnetic resonance imaging of brain shift[END_REF][START_REF] Hastreiter | Strategies for brain shift evaluation[END_REF] Hence, it is necessary to design a novel technique combining the pre-operative data with intra-operative tracking which would be used to drive a simulation based on an accurate biomechanical model. Such technique would dramatically improve the accuracy of navigation, thus reducing the risk to the patient. The partial removal of a kidney is a complex procedure that requires a highly experienced surgeon to be successful. A consequence of this is that several surgeons prefer doing a total nephrectomy in order to reduce the risk of failure. Hence, it is important to design Corresponding author: Stéphane Cotin, e-mail: stephane.cotin@inria.fr a technique that predicts accurately the tumor location intra-operatively in order to remove as less as possible healthy kidney tissue and preserve its function.

In this work, we show that FS-FEM has similar accuracy as the standard FEM in the simulations of the brain-shift and of the kidney's deformation, while it outperforms the standard FEM in a bending beam scenario.

METHOD

A three-dimensional computational domain Ω is divided into n e tetrahedral elements. Over each element the displacement field is defined as follows

u(x) = N d I=1 N I u I (1)
and the compatible strain is given by

(x) = N d I=1 B I u I . ( 2
)
N I is the shape function at node I, N d is the total number of independent nodal basis shape functions, and B I is the strain displacement matrix at node I which size is (6 × 3). u and are the displacement and the strain computed using the standard FEM. The S-FEM smooths the strain over the boundaries of the elements using the so-called smoothing domains. For the face-based smoothed finite element method, these domains are associated to the faces of the elements and are obtained by connecting portions of elements sharing a face k (see Fig. 1 left). They are constructed on top of the element mesh, as a non overlapping and no gap subdivision of Ω such that :

Ω = ∪ ns k=1 Ω s k
where n s is the total number of smoothing domains. In our case, the number of smoothing domains is equal to the number of faces i.e. n s = n f where n f is the number of faces. The smoothed strain ¯ is obtained through a spatial average of the standard strain over the smoothing domains :

¯ (x k ) = Ω s k (x) Ŵ (x k -x) dΩ. (3)
Ŵ (x k -x) is the smoothing (or weight) function, that is compactly supported, positive, centered at x k and verifying Ω s k Ŵ (x k -x) dΩ = 1. The most common expression for it is given by :

Ŵ (x k -x) = 1 V s k if x ∈ Ω s k ∪ Γ s k 0 elsewhere. (4)
V s k is the volume of the smoothing domain Ω s k that is obtained by

V s k = 1 4 n e k m=1
V em . In this formula, V em is the volume of the m-th element associated to face k and n e k is the number of elements associated to face k (e.g. n e k = 1 for a boundary face and n e k = 2 for an inner face). Plugging ( 4) into (3), the smoothed strain field obtained is :

¯ (x k ) = ¯ = 1 V s k Ω s k (x) dΩ. (5)
In other words, the smoothed strain displacement matrix linked to face k weights and averages the strain displacement matrix B em of the two elements sharing the face k (e.g. m = 1 and m = 2). In the case of a tetrahedral mesh, the smoothing is performed during the local assembling process (see Fig. 1 right for details) which is expressed as

Bk = 1 V s k n e k m=1 1 4 V em B em . ( 6
)
As a result of the assembling process, Bk is a (6 × 12) matrix for a boundary face and a (6 × 15) matrix for an inner face. The entries of Bk are constant over the smoothing domain Ω s k . For linear tetrahedral elements, the strain is constant over the smoothing domains Ω s k and the local smoothed stiffness matrix Kk is given by : Kk =

Ω s k BT k C Bk dΩ = BT k C Bk V s k . (7)
with the standard constitutive matrix C for linear elasticity. For an inner face, Kk is a (15 × 15) matrix, while it is a (12 × 12) matrix (same as in standard FEM) for a boundary face.

Once the local matrix is computed, the global assembly differs from that of the standard FEM since on the local level, the matrix sizes are different which results in an increased density of the system matrix on the global level. The global smoothed stiffness matrix K is computed by an assembling process as :

K = ns k=1 Kk . ( 8
)
Since we are in the linear case, this assembling can be done at the first time step once for all because it does not change over time. K is a (n d × n d ) matrix where n d is the number of degrees of freedom of the problem (e.g. n d = 3 × n n where n n is the number of nodes). Having the global matrix, the use of a direct solver such as the LDL factorization is possible. This will further push the simulation closer to the real-time performance as it will allow for pre-computing of the factorization of the matrix before the simulation. Another effect is the significant reduction of the computational load required directly during the simulation.

RESULTS

To verify the implementation of the FS-FEM performed in SOFA * , a simulation of bending is performed on a cantilever beam with different mesh resolutions (1183, 2210, 4326, 8743 and 12410 elements) and compared to the FEM solution computed using a fine mesh with 22700 elements. The beam is statically loaded by a force of 0.0005 N applied vertically to the right face of the beam having dimensions 0.178 × 0.01 × 0.01 m. The Young's modulus and Poisson's ratio are set to 100 kPa and 0.3, respectively.

Fig. 2 left shows the convergence of both FEM and FS-FEM for different densities of the mesh. The relative error e is computed as

e 2 = Ω (u SF EM -u F EM ref ) T (u SF EM -u F EM ref ) dΩ Ω u F EM ref T u F EM ref dΩ . (9)
Fig. 2 right shows the beam deflection obtained from FEM (in green) and FS-FEM (in blue) using the mesh with 2210 elements compared to the reference FEM solution (in red) employing the mesh with 22710 elements. It is observed that using the same mesh resolutions, the solutions obtained by FS-FEM are closer to the reference solution when compared to FEM. This confirms our hypothesis that a significantly coarser mesh is sufficient for FS-FEM to achieve the same error as the FEM using finer mesh. This can be in turn exploited in design of a real-time simulation. Further, when using the iterative solver, the rate of convergence is higher when using FS-FEM. It should be emphasized that the amount by which the error decreases is proportional to the difference in number of DOFs required by FEM and FS-FEM. We performed the same simulation on the bending beam scenario, but this time employing a FEM with a tetrahedral mesh (FEM-T4), a FEM with an hexahedral mesh (FEM-H8) and a FS-FEM using a tetrahedral mesh (FS-FEM-T4). We run the simulation with 5 different meshes as in the previous experiment. For the tetrahedral meshes we used 373, 693, 1311, 2536 and 3484 nodes, and for the hexahedral ones 76, 333, 845, 1775 and 3852 nodes. The reference solution in this experiment is the result of FEM-H8 having a mesh of 28 000 elements. In Fig. 3 are shown the relative errors computed using equation (9). It is observed that the curve of FS-FEM-T4 is between that of FEM-T4 and FEM-H8. Hence, FS-FEM using T4 elements is closer to a simulation having an hexahedral topology than standard FEM-T4. An hexahedral topology is known for its accuracy but also hexahedral meshes are difficult to generate for complex geometries. Hence, a tetrahedral topology is preferred but a high resolution is needed for a high accuracy of the solution. The FS-FEM allows for reducing the mesh resolution while preserving a good accuracy. Overall, the FS-FEM performs better than standard FEM-T4 in a bending scenario. In the next experiment we analyze the behavior of the FS-FEM in a traction and compression scenario that are common situations in medical contexts. To that purpose we simulated a 3D cube undergoing a uniform load at its upper face with fixed boundary conditions at its lower face as depicted in Fig. 4 left. The Young's modulus and Poisson's ratio of the cube are respectively 2 kPa and 0.45. We used FEM-T4 and FS-FEM-T4 in meshes having 1423, 2146, 4536, 9553 and 14553 elements. To analyze this compression scenario, we mapped a grid to the upper part of the cube where the deformation is the biggest (cf Fig. 4 right). We computed the relative errors (cf equation ( 9)) comparing to a reference obtained using a FEM-T4 having 31536 tetrahedral elements. The values are shown in Fig. 5. It is observed that the FS-FEM-T4 has a slightly lower relative error than the FEM-T4 but globally we can say that there is no difference between the two methods in a compression scenario. We find similar results when traction is applied on that cube. Since the strain field is homogeneous over the cube, it was expected to obtain the same result no matter the method used.

Because of its properties, the FS-FEM seems to be a good choice for simulation which requires both high accuracy and good performance. Since this is typically the case of the brain-shift simulation and of the kidney simulation during a partial nephrectomy, we employ the FS-FEM in this two contexts. Fig. 6 shows an image of a human head which was previously used in work of A.Bilger et al. [START_REF] Bilger | Biomechanical Simulation of Electrode Migration for Deep Brain Stimulation[END_REF] The anatomy of the brain comes from a generic atlas and its deformations are modeled with the FS-FEM using a tetrahedral mesh. The right image of Fig. 6, shows the result of the simulation, e.g., the completed brain shift predicted by the model. The parameters Young's modulus, Poisson's ratio and mass are set to 3 kPa, 0.4 and 0.808 respectively. Fig. 7 left shows the position of several tumors inside the brain at rest position. Fig. 7 right shows the positions after the deformation due to brain-shift, predicted by a simulation using a linear FEM (blue points), a FS-FEM (red points) and a non-linear FEM (green points). The mesh used for the three methods is a mesh having 7924 tetrahedral elements per hemisphere. Unexpectedly, it is observed that the blue and the red points overlap in most cases. This means that the linear FEM and the FS-FEM do not perform in a significantly different way in this scenario. Furthermore, the red points are in general far from the two other colors. It is very likely that the elements rotate during the deformation and this leads to geometrical non-linearities that are treated better by a non-linear FEM than by the linear FEM or FS-FEM. Overall, the FS-FEM performs qualitatively in the same manner as standard linear FEM in this scenario. We now perform a similar simulation and analysis over a kidney under gravity. In Fig. 8, a tetrahedral mesh of a kidney having 727 elements is shown. Its Young's modulus and Poisson's ratio are set to 3 kPa and 0.3, respectively. We used a direct solver (LDL factorization) in order to be compatible with real-time simulations.

We want to verify whether the FS-FEM behaves in a similar manner as in the brain. To that purpose, we perform a simulation of the kidney's deformation using the FS-FEM, the linear FEM and a non-linear FEM over a mesh having 727 elements. We placed 8 tumors inside the organ. The position of the tumors after the deformation are shown in Fig. 9 left. The red spheres stand for the FS-FEM, the green ones for the linear FEM and the blue ones for the non-linear FEM. It is observed that all the tumors are very close to each other and as in the brain, in some cases the green spheres are between the red and the blue ones. This means that for the kidney scenario, the FS-FEM performs in a similar manner to the linear FEM, just as in the brain-shit scenario. In Fig. 9 right are shown the results of a simulation comparing the linear FEM and the FS-FEM employing a mesh of 727 elements to a reference solution computed with linear FEM using 3224 elements. The red spheres stand for the FS-FEM, the green ones for the linear FEM and the blue ones for the linear FEM with a very fine mesh. It is observed that the green and the red spheres are overlapping. Hence, it is difficult to draw a conclusion on the most appropriate method for this context. Overall, there is no difference between the FS-FEM and the linear FEM in this scenario. Figure 9. Left : Tumor positions after deforming the kidney under gravity using a FS-FEM having a mesh of 727 elements (red spheres), linear FEM with the same mesh (green spheres) and non-linear FEM with the same mesh (blue spheres). Right : Tumor positions after deforming the kidney under gravity using a FS-FEM having a mesh of 727 elements (red spheres), linear FEM with the same mesh (green spheres) and linear FEM as a reference with a mesh of 3224 elements (blue spheres).

CONCLUSIONS AND PERSPECTIVES

In a bending scenario, the FS-FEM-T4 is proved to better approximate the displacement solution than the standard FEM-T4. Moreover, the FS-FEM having a tetrahedral mesh is closer to a simulation having an hexahedral topology than standard FEM-T4. Hence, the FS-FEM allows for reducing the mesh resolution while preserving a good accuracy of the solution of a bending scenario.

Not surprisingly, our implementation of the FS-FEM presents no advantage on compression and traction scenarios. We have seen that the reason for this is the homogeneous strain field. Indeed, the smoothing has no impact in such cases. Contrary to our expectations the results on the brain-shift and on the kidney where rather unexpected. We did not find a significant improvement when using the FS-FEM in this contexts.

Hence, further work will concentrate on simulating other organs with different elasticities and shapes such as the liver or the eye in the context of liver tumor removal and retinal surgery. To the best of our knowledge, the co-rotational FEM is the fastest method accounting for geometrical non-linearities and this is a very common situation in medical contexts. As a consequence, our next step is the development of the co-rotational FS-FEM that has not been presented in the literature before.

Figure 1 .

 1 Figure 1. Left: Two adjacent tetrahedral elements and the smoothing domain (shaded) formed based on their interface k in FS-FEM. Right: Example of the assembling process of the smoothed strain displacement matrix.

Figure 2 .

 2 Figure 2. Left : Relative errors of FEM and FS-FEM solution compared to the reference FEM solution (with very fine mesh composed of T4 elements). Right : Zoom-in at the tip of the beams.

Figure 3 .

 3 Figure 3. Relative errors in log scale of FEM-T4, FEM-H8 and FS-FEM-T4 solution compared to the reference FEM-H8 solution (with very fine mesh composed of H8 elements).

Figure 4 .

 4 Figure 4. Left : 3D cube under compression. Right : Tetrahedral mesh of the 3D cube and mapped grid in red for the computation of the relative errors.

Figure 5 .

 5 Figure 5. Relative errors in log scale of FEM-T4 and FS-FEM-T4 in cube compression

Figure 6 .

 6 Figure 6. Brain at rest position and shifted brain simulated with FS-FEM.

Figure 7 .

 7 Figure 7. Grid representing tumor positions using a mesh of 7924 elements for linear FEM in blue, FS-FEM in red and non-linear FEM in green after the brain-shift. Left : rest position. Right : position after the deformation due to brain-shift

Figure 8 .

 8 Figure 8. 3D mesh of a kidney with fixed constraints (pink points).