A Novel Robust Approach for Computing DE-9IM Matrices Based on Space Partition and Integer Coordinates

: A novel approach for a robust computation of positional relations of two-dimensional geometric features is presented which guarantees reliable results, provided that the initial data is valid. The method is based on the use of integer coordinates and a method to generate a complete, gap-less and non-overlapping spatial decomposition. The spatial relationships of two geometric features are then represented using DE-9IM matrices. These allow the spatial relationships to be represented compactly. The DE-9IM matrices are based on the spatial decomposition using explicit neighborhood relations. No further geometric calculations are required for their computation. Based on comparative tests, it could be proven that this approach, up to a predictable limit, provides correct results and thus offers advantages over classical methods for the calculation of spatial relationships. This novel method can be used in all ﬁelds, especially where guaranteed reliable results are required.


Introduction
An important functionality in GIS or CAD is the calculation of spatial relations of geographic or geometric objects (features). There are many publications dealing with the efficient and correct representation of these relations. This article is not about the different types of representations and their advantages and disadvantages, it is about the correct calculation of spatial relations.
On the one hand, a correct calculation of spatial relations is based on exact geometric computations and, on the other hand, on consistent topological relations. Exact geometric computations cannot be guaranteed in every case if floating point numbers are used. Nevertheless, floating point numbers are mostly used for computation (see e.g., Ref. [1]). This problem can be easily solved by using integer numbers. Also, in any case it must be guaranteed that the topology is consistent with the geometry, this is especially a problem when the features to be compared come from different data sources.
We developed an algorithm that overcomes these problems. The method is based on three fundamental principles: 1.
The use of integer coordinates and the representation of intersections as positions on edges using rational numbers.

2.
An algorithm for generating a complete, gapless and non-overlapping spatial decomposition. 3.
The generation of DE-9IM matrices without further geometric calculations by set operations of special types.
The main goal of our algorithm is to guarantee consistency between geometrical and topological relations of 2D features. In the following we will prove that the presented approach can achieve this for geospatial data, when determine topological relations  or checking geometric features for topological inconsistencies. To illustrate the practical relevance, we will show that conventional geo tools fail to achieve the main goal under certain conditions. In our analysis, we compare our approach with the NetTopologySuite library and the SpatiaLite database. Although this is only a small selection, it should only serve to highlight the basic problem.
As an illustration, point P(A) (a point of feature A) and edge E(3) (an edge of feature 3) are highlighted in Figure 1a. A question about point P(A) would be: is it in feature 1 or 2, or even directly on the edge between the two features? Similarly for edge E(3): does this cross feature B or is it on an edge? Topology and geometry should be consistent for clear resolutions of these questions. In step one, the vertex coordinates of the original geometries are converted into integer coordinates. This step guarantees that after the conversion the geometric calculations are unambiguous. For this purpose, the coordinates are scaled so that they can be represented as integers without changing their relative position to each other (Figure 1b). In this step it is possible to check the validity of the given features (now represented with integer coordinates), based on their properties such as length, area and orientation.
In the second step, a spatial decomposition consisting of vertices, edges and faces is generated (Figure 1d). As a first sub-step, a triangulation of all points is generated (Figure 1c). This triangulation forms the base mesh and guarantees that the spatial decom-position is gap-less. Then the edges of the features are reconstructed in this base mesh. This special algorithm guarantees that the model represents the correct topological relations of the initial geometries. The generated spatial decomposition is then complete, gap-less and non-overlapping. To achieve this, the edge intersections are represented as ordered positions on the edges using rational numbers. Since their number range is predictable, no data type for arbitrarily large integers is needed.
The generated spatial decomposition contains all geometrical and topological information of the source geometries; thus it is possible to determine the spatial relations without further geometrical calculations. To facilitate the queries, special 0, 1 and 2 dimensional objects are created, each representing a source feature. This is the last instance where incorrect source geometries can be detected and filtered out. Then, the DE-9IM matrices can be generated efficiently without geometric calculations, purely by set operations.

Related Research
The problem area of exact results of geometrical calculations is considered in many publications (e.g., Refs. [2][3][4][5]). It is called "exact geometry", "exact computation", "geometric predicates" or "exact geometric computation". An overview of the problem area is provided in [6]. The problem is usually considered from the computational side. One can say that integer coordinates represent the best solution, but at the expense of computational performance. This then leads, for example, to the use of floating point numbers and special computational techniques [4]. We commit ourselves to the use of integers, since only these make it possible to represent the source data exactly (Section 3.2).
The problem actually arises when calculating intersections [7]. Coordinates of intersections can be represented lossless by integers or floating point numbers only in exceptional cases. As a general rule, rational numbers are required for this purpose [2]. Thus, topological inconsistencies may occur after the computation. Therefore, we use a special method for space decomposition, presented in [8] for the context of digital building models. This guarantees that the topology of the created model is free of inconsistencies. The problem that coordinates of intersections are not exact is not avoided by this. They have to be calculated at least for representation, but this does not change the previously correctly generated topology. One could, for example, flag these intersections so that they are easy to identify.
In contrast to the approach in [9], our goal is not to validate data for topological consistency, but to ensure that the topology between the processed features is consistent for determining spatial relationships. This is similar, for example, to the approach in [10], whose focus is on topologically consistent data models. Our algorithm cannot guarantee that the model is still topologically consistent after the coordinates of the intersections have been computed. However, these coordinates are not needed for the determination of the spatial relations.
The topologically correct model we generate is then used to generate the DE-9IM matrices. How to map topological relations is explained in [11]. In [12] the advantages of 9IM matrices are shown. The extension to the output of the dimension of the intersection can be found in [13]. Later, the output of the dimension has extended to the 9IM matrices, leading to dimensionally extend 9IM matrices (DE-9IM). Their definition is represented, for example, in [14]. Because of standardization ( [14], pp. 34-35), we use the DE-9IM matrices, although there have been further developments [15].
Some first ideas of the presented algorithm have already been presented in [16] in the context of object recognition from terrestrial laser scanning. This shows that our method is very broadly applicable. However, the focus of this article is to present a fundamentally new concept for the exact calculation of topological predicates.

Materials and Methods
In Figure 2 the preparation steps of the geometric features are depicted, in order to use them for the generation of DE-9IM matrices. These are explained in more detail below. All input features are reconstructed in a space partitioning as described below.

Geometric Features
The Open Geospatial Consortium (OGC) "Simple Features Specification" [14] is a wellestablished conceptual standard for various (mostly two-dimensional) geometric features. The standard also includes the specifications for the DE-9IM matrices and for the Well Known Text (WKT) representation used in our algorithm and test scenario. Table 1 shows the used Simple Features and their WKT representation. The third column shows the geometric dimension of the corresponding types. This is important for the subsequent object creation for DE-9IM matrix determination. Our test contains only basic OGC WKT types, however all basic geometric operations, that are of our interest, are investigable with this limited set of types.  The representation of geometric features as WKT is a widely used standard, so data exchange is straightforward.

Integer Coordinates
The algorithm is using integer coordinates. This is the only way to guarantee that all calculations are exact and that there are no inconsistencies between topology and geometry. To do this, the vertex coordinates of the geometric features must first be converted into integers. The coordinates are scaled and translated in such a way that the original relative spatial relations are preserved.
The coordinate systems used are not described in detail, Cartesian coordinates are assumed. For projected coordinates, a conformal coordinate system should be used to avoid distorting the spatial relationships. In the numerical examples it is assumed that meter is used as the unit of length.
The source data, the set of geometric features (G), are represented as WKT simple features. Thus, the coordinates are presented as string. These strings have to be converted lossless into a number format for calculations. There are different options for this depending on the computer platform used. The important issue is to ensure that there is no loss of quality during this step. For example, it would be sufficient to use the decimal128 type [17] or an equivalent. It is essential not to use binary floating point numbers, since in this case the coordinates can already be corrupted when they are imported. This can be seen easily in the examples in Table 2. Caused by the base change it comes with rational numbers inevitably to conversion losses, since many decimal finite rational numbers become infinitely periodic in the binary system.

Decimal
Binary Decimal The imported OGC Simple Features are stored in an equivalent data type (Point, Linestring, Multilinestring, Polygon). The coordinates are imported lossless (e.g., with decimal128). In addition, the minimum bounding box of all points from the features in G and their maximum scaling factor s max of the decimal places are determined (e.g., The factor s for the number 12.345 is 1000). The bounding box directly defines the smallest coordinate value (x min , y min ) and the largest extent r max (Equation (1)).
To store the coordinates as integers, it is necessary to specify the integer type to be used. One could simply use an arbitrary large integer type for the integer coordinates and have no problems converting the coordinates [6] (pp. 79-81). However, this would not be practical because it would lead to an unpredictable increase in the amount of memory needed which also effects the computation time. So it is reasonable to use native integer types. Consequently, it must first be checked whether this is possible with the usual 32 bit and 64 bit integer types.
For this, first consider the essential formula for the following algorithm. This is the determinant calculation of a matrix of three points (Equation (2)). It is used to determine whether three points are collinear, or if they are not, to determine the orientation of a triangle made of the three points [3,4].
x a x b x c y a y b y c 1 1 1 The geometric origin can be used to determine the maximum coordinate value for the determinants calculation. Figure 3 shows triangles which represent the set of all possible largest triangles within the positive integer coordinate range. The formula for calculating the triangle area (Equation (3)) makes it easy to deduce that all triangles have the same area.
Thus, the determinant can result in a maximum of twice the value (positive and negative) of the triangle area (Equation (4)). From this, Equation (5) can be derived, by which the maximum possible coordinate value c max can be determined from the integer type used for the determinant calculation. Now we have to consider how integer coordinates are calculated from the given coordinates (Equation (6)).
Equation (7) illustrates that the size of the integer coordinates depends on the maximum extent r max and the maximum scale s max . If their product is larger than the maximum available coordinate value c max , data loss occurs because the coordinates are truncated by their remaining decimal places after scaling. The conversion is no longer bijective in this case.
If we use the 64 bit integer type for the determinant calculation, with a number range from [−2 63 to 2 63 − 1], we get for the largest possible coordinate value. This number is greater than the maximum positive numerical value of 32 bit integers. Thus, they can be used for the coordinates. Assuming a necessary computational accuracy in the millimeter range (s max = 1000), the extent of an region must not be larger than 2,147,483.645 m, i.e., about 2000 km, in order to still calculate correctly with 32 bit integers. Thus, this range is only recommended for small-scale applications. However, in order to be able to use the algorithm in all geospatial cases, it is advisable to work with 64 bit integers for the coordinates. This would even allow an extent of approx. 9,000,000,000,000 km with millimeter accuracy. Since the largest distance possible on Earth is approx. 40,000 km, 64 bit integers can map all coordinates up to a computational accuracy of 11 decimal places (10 pm). Thus, this number range is sufficient for all use cases in the GIS area. The disadvantage is that for the determinants 128 bit integers are needed, which are not always available natively in all programming languages or platforms.
After the conversion, the geometric features are stored with their integer coordinates, and a point set of all integer coordinates P is formed (Figure 1b). During the generation of integer objects it must be checked that none of the features in G is degenerated (Figure 4a) or if the orientation has changed (Figure 4b). A degeneration of a feature can be recognized by the fact that the end points of a line become congruent or it gets an area of zero, where a change of orientation is associated with a sign change of the area. In this case the feature must be discarded and cannot be analyzed with this method. For this reason, we use 64 bit integers. Alternatively, it is possible to decide based on the spatial extent (Equation (7)). Due to the previously explained usability of 64 bit integers, it can be assumed that degenerated features were already degenerated before the conversion.

Intersections as Rational Positions
Integer numbers are used for the coordinates of P in the basic mesh M, but integer numbers cannot be used for intersection points, when reconstructing edges E G in M, since intersection points are integer numbers only in exceptional cases. Even floating point numbers could not represent their positions correctly. Rational numbers are required.
We have decided not to calculate the intersection points as coordinates. We represent them as ordered positions on an edge ( Figure 5) of the base mesh M or a feature G i . However, this also requires rational numbers which are represented as a fraction of two integers. The given end points of the edges have integer values as coordinates. As a consequence, the position of the intersection point on these edges can be exactly expressed by a fraction of two integer values. This is true to all intersection points on an edge so that the order of these intersection points on an edge can be determined exactly as well. Equations (9) and (10) Equation (11) shows an example for the position calculation, which is also depicted in Figure 6. Figure 6 also shows the geometric meaning of pos ab and pos cd , whose values were derived from the ratio of the partial distances to the total distances. Since determinant calculations (similar to Equation (2)) are needed, the integer type resulting from this calculation (64 bit for 32 bit coordinates resp. 128 bit for 64 bit coordinates) must be used.
It is important to note that the coordinates of the intersection points are not required for our algorithm. The presented algorithm uses only relative positions on edges. However, absolute xy-coordinates can be calculated, for example, for display purposes. But the calculated xy-coordinates of the intersection points should not be used for further calculations, due to potential computational inconsistencies between geometry and topology.

Space Partition
The integer coordinates P now form the basis of a Delaunay triangulation (Figure 1c). Theoretically, a simple triangulation would suffice, but it was proved in [18] that a Delaunay triangulation reduces the number of intersections. This triangulation provides the base mesh and remains structural unchanged in the next upcoming steps.
Then, successively the edges E G of the geometric features G can be inserted into the base mesh M (Figure 7). This step is subdivided into the following substeps: 1.
Since all points are already contained in the base mesh M, the starting point of the edge has to be searched first (Figure 7a).

2.
Then the base mesh edge E Mi to the right of or on the feature edge E G i is searched for (Figure 7b). Once all sector triangle edges are inserted, the intersection points of these edges within the individual triangles of the basic mesh can be determined. The positions of the intersection points are stored as rational numbers as explained in Section 3.3.
In Figure 8 the individual steps for identifying edge intersections E G with E G within the triangles of M are represented. For this purpose, an iteration is performed over the triangles of the base mesh. For each triangle the following substeps are processed:

1.
Determine all intersection points of feature edges on the triangle edge and sort counterclockwise around the triangle. (Figure 8a).

2.
Determine the inner intersection points according to the scheme in (Figure 8b). Now each point (P), each half-edge (HE ) and each facet (F ) is checked to which geometric feature it belongs. Thereby an element can belong to none, one or more geometric features. Figure 9 shows how to find all elements of an areal geometric feature. First, all half-edges that lie on the boundary of the feature are determined (Figure 9a). The adjacent facets in the interior (Figure 9b) belong to the geometric feature. Then, the set of facets is extended by their neighboring facets until no new ones are found (Figure 9c). Finally, all points, half-edges and facets that were detected or are part of the detected elements are assigned to the feature. At this stage we have a complete, gap-less and non-overlapping two-dimensional decomposition, generated fully automatically from the original data as seen in Figure 1.

Point, Line and Region Types
To increase the efficiency in using the space decomposition to determine the DE-9IM matrices, a preparation step is required.
The comparison to determine the DE-9IM matrices operates on three different types: Point, Line, and Region. They contain fields representing the set of Interior and Boundary elements ( Table 3). The types indicate the geometric dimension of the geometric features to be represented.

Name Mesh Element Point Line Region
Interior0 Vertex HalfEdge HalfEdge x Since in the last step of the spatial decomposition the elements of the base mesh were assigned to belong to a geometric feature, now the Point, Line and Region objects can be easily formed, depending on the OGC type of the associated features.
At this stage, the preparation for generating the DE-9IM matrices is complete. Only the generated Point, Line and Region objects with their associated elements are needed.

DE-9IM Matrices
With the now existing Point, Line and Region objects, one can easily determine the DE-9IM matrices. They have 0, 1, and 2 dimensional Interior (I) fields, and 0 and 1 dimensional Boundary (B) fields ( Table 3). The Exterior (E) is made up of those elements that do not belong to the object.
In Equation (12) it can be seen how DE-9IM matrices are defined. The set operations of the definition can be applied directly to the previously created Point, Line and Region objects.
The dim operator can be implemented in such a way that it always returns the highest dimension that has a non-empty intersection between the comparison objects. The tests in the last column and row of the matrix have to be adapted, because the exterior is not stored directly. However, the tests can be easily adapted. Since we are working in a bounded region, the exterior is simply that part of the overall set which is not part of the interior and the boundary. Therefore, the types always have combined InteriorBoundary fields to allow easy querying (Table 3).
With this, the DE-9IM matrix can be easily created. If a simple binary 9IM matrix is needed, it can be easily converted. At all positions occupied by a 0, 1 or 2, a T (true) is set. Figure 10 demonstrates how these set operations work. Figure 10a shows the region object of feature 2 in Figures 1a and 10b shows the region object of feature A. Figure 10c then depicts the determined intersection. The DE-9IM matrix for this example is presented in Equation (13). In words, this is an overlaps relation. The individual cells of the matrix are created as shown in Table 4.   Figure 10 (Equation (13)).

Val. Definition Calculation
is always 2, because the exteriors of objects would always overlap if the space is larger than the objects.

Results
The examples presented in this section investigate the application of the algorithm. They follow a specific test scenario. This test scenario is also conducted in established GIS tools. Then, the results are compared to demonstrate the capabilities of the presented approach. The test scenario was developed and tested in the context of a master's thesis [19]. The rough test flow is represented in Figure 11 and will be presented in the following.  Figure 11. Test procedure in rough diagram form.

Examined Tools
A well-known opensource library for computing spatial predicates is the JTS Topology Suite (JTS) [1]. It implements the Simple Features Specification for SQL of OGC [20]. Since JTS can only work within a Java environment, there are also ports in other programming languages. A port to C++ is the GeometryEngine-OpenSource (GEOS) [21]. This port makes the functionality also usable by C/C++ tools. GEOS is for example used in the popular opensource spatial database PostGIS as well as in the opensource GIS tool QGIS. The NetTopologySuite (NTS) [22] is a direct port of the JTS targeting the .Net environment.
Since JTS and its ports are used so extensively, our algorithm is tested directly against this library respectively its port. Because of the ease of implementation in our used environment (Windows, .Net Framework) we have chosen NTS [22] and SpatiaLite (SpL.) [23]. The latter is a spatial database that is directly based on SQLite and can therefore run without a server installation.
These two tools are compared with the implementation of our algorithm. We use our algorithm in a version with 32 bit coordinates and a version with 64 bit coordinates. The 32 bit version is used to show possible problems with spatial tests in areas with large extent or high accuracy.
It is important that the tested tools meet the following conditions: 1. Data transfer by using OGC WKT types.

Test Design
The goal of the tests is to determine if any invalid results will be generated by the exanimated software tools. If so, it is further of interest for which geometric constellations (intersections) and (homeomorphic) variants the tests fail.

Tested Geometries
The three types Point, Line and Region have different properties (interior, exterior, boundary), so there are many possible combinations. For our tests, we decided to use the selection of combinations shown in Figure 12  The configuration of the tests was chosen to produce as many variations of DE-9IM matrices as possible. This is represented in Table 5.
It is important to mention, that not all possible relationships of features are covered in this document. Only cases in which different DE-9IM matrices occur. Tests with the reversed order of features are not represented in Table 5, but they are calculated and checked with every test. This is easily done because the DE-9IM matrix of the inverted comparison corresponds to the transposed original matrix. Line Line FF1FF0102 7 Line Line 1FFF0FFF2 8 Line Line 101FF0FF2 9 Line Line 101F00FF2 10 Line Line FF10F0102 The actual geometries of the objects are set to have integer coordinates. This results in a s max of 0 for the untransformed features. It guarantees that the base tests match the DE-9IM matrices represented in Table 5. The features themselves are stored as OGC WKT (see Table 6). The extent (r max ) of all test areas is not greater than 1056 so that all tested tools compute these base tests correctly. Table 6. WKT Geometries of the first 10 tests. Since the basic tests are designed in such a way that they should produce correct results with every tool to be tested, they cannot be used directly to point out problems with the tools. In the following it is shown how the features are changed by a homeomorphic mapping for the tests in such a way that the topology remains, but other coordinates are used in the computation. When creating the WKT features to be tested, care must be taken that no errors are produced just by the type of calculation. For this reason also decimal128 numbers are used. This avoids invalid test geometries due to number conversion when applying the homeomorphic mapping for the test cases.

Test # Geometry
After creating the transformed WKT features, all tools to be tested get them passed as a string.

Translation
A simple way to perform homeomorphic mapping is the translation. It transforms the coordinates of the base features as shown in Equation (14).
The values for T can be seen in Table 7. There are in total 60 different values. These are chosen so that, firstly, the number of decimal places increases from 1-6. Secondly, the integer part increases from 0 to 1,000,000 by a power of ten correspondingly. In this way, coordinates up to an accuracy in the nanometer range are obtained during translation. Also region extents are obtained which are larger than the maximum possible distance on the earth's surface (40,000 km). Thus, the tests should cover all GIS application areas.
An example of a translation of features shows Listing 1. The example uses test number 4.

Scaling
Another simple homeomorphic mapping to perform is scaling. Point coordinates are scaled by the application of Equation (15) with a value S = 0.
The values for S in the tested examples are identical to the values for translation and can be seen in Table 7.
In the case of scaling, as well as in the case of translation, coordinates are obtained down to an accuracy in the nanometer range. In the case of scaling, the resulting range also covers all possible GIS application.
In Listing 2 the effect of scaling is represented by the example of the test number 4.

Test Results
This section shows the results of the tests. Listing 3 shows the total number of tests which results from the number of comparison tests (33) and the number of translation and scaling values (60 each). Table 8 demonstrates for which tests and for how many translation or scaling values the tested tools fail or calculate with data loss. Failure means that the test of the two features from Figure 12 resp. We can see that the translation has no influence on the result of our algorithm, because the translation does not change the maximum extent (r max ) of the test region. However, this is not the case for the other tools. They both fail with the same translation values and test cases, which refers to the identical computational kernel. It is significant that always such tests fail, which are based on a test "point on straight line" resp. "straight line collinear to straight line".    Similarly, it can be seen that our algorithm in the 32-bit version fails exactly the same tests as the other tools, because of the data loss in the integer conversion of the coordinates when the scaling factor is too large (s max r max > c max ). However, unlike the other tools, this could be accurately predicted beforehand with Equations (7). Similarly, our test in the 64 bit version exhibits the predicted behavior and does not fail in any of the transformation tests.
The behavior across all tests can be seen in Figures 13 and 14. It can be clearly seen, that only our algorithm (64 bit) produces correct results. Other tools always produce unreliable results (at least as soon as "point to straight line" resp. "straight line collinear to straight line" tests are involved).

Discussion
As the presented tests proved, our novel approach offers a better alternative in determining spatial relations of geometric features.
The other two tested tools always failed when the tests involve "point on line" respectively "line collinear to line" tests. Both tools are based on the JTS. Therefore, it is not possible to determine whether the errors are due to the implementation or the use of floating point numbers. But as the examples in Table 2 show, this problem is fundamental since binary floating point numbers are not suitable for lossless representation of the digits after the decimal point.
In this context, the performance of our algorithm was not tested in relation to the other tools for a large amount of data since this is only a proof of concept. However, our algorithm should be at least used when the reliability of the results is important, regardless of the performance.
The initial claim was that our algorithm always produces correct results when comparing geometric features. Based on the tests, this was clearly proven for the algorithm with 64 bit integers. Although these tests do not cover all possible combinations and region extents, the proof in Section 3.2 shows that the possible region extent is perfectly sufficient for all GIS use cases. For smaller regions, the 32 bit implementation can also be used. This has the advantage that no 128 bit integers are required during the computation. In addition, it can be automatically identified whether the computation must be better accomplished with the 64 bit version, if a data loss is expected.
An important part of the presented approach is also the space decomposition from Section 3.4 because only this guarantees the correct functionality of the algorithm. In addition, this method of space decomposition has the advantage that, if it is used consistently for model generation, the consistency of topology and geometry can be guaranteed.
Another advantage of the presented approach is that the used Point, Line and Region types can also be used, for example, to easily generate the results of Boolean operations of two geometries.
All in all, our algorithm offers a real alternative to established methods in the comparison of geometric features and additionally has a great potential to extend the range of applications.

Outlook
The approach presented in this paper is restricted to the 2-dimensional space. It shows that the application of a space partitioning concept works and give results without topological inconsistencies. In a space partitioning concept, objects are not tested against each other with all well know problems of the fact that the computer cannot represent all real numbers.
It is a challenge to apply the concept of space partitioning to the 3-dimensional space. Two different ways can be applied: development of a theory and a tool which offers functionality to construct 3-dimensional space partitions from scratch and the transformation of existing boundary representations into a space partition. The first approach is under development [24]. The second approach is the basis of the research presented in this paper. Investigations of the application of this approach to the 3-dimensional space are also under development [25].
There are still challenges in this field. Research in the 3-dimensional space has the potential to close the gap between applications of GIS and digital buildings models in the field of Building Information Modeling. We accept these challenges and focus in our future collaboration on a next step to improve digital geometrical modeling.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: