Ewa Skubalska-Rafaj Lowicz

Aneta Górniak
email: aneta.gorniak@pwr.edu.pl

Kohonen SOM for Image Slides Sequencing

Keywords: image sequencing, SOM Kohonen map, Euclidean travelling salesman problem, image processing, image registration, 3-D image reconstruction

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In this paper we present a new approach to sequencing of 2-D slide images with the goal of 3-D object reconstructions. We focus our attention on microscope images obtained from serial sections of biological tissues.It is assumed that subsequent section images may be located inadequately in the image series. Translations, rotations and even small non-rigid transformations of the object of interest can occur.

Image registration is a very important step in processing sequences of images [START_REF] Flusser | Image registration methods: a survey[END_REF]. A large number of algorithms have been developed to perform registration of medical and biological images. Usually the image registration is the process of aligning multiple images representing the same scene that were captured at a different time, at or by a different set of modalities, or by aligment of successive anatomical slices [START_REF] Flusser | Image registration methods: a survey[END_REF]. We concentrate on the image registration in the case of 3-D tissue reconstruction [START_REF] Lippolis | Automatic registration of multi-modal microscopy images for integrative analysis of prostate tissue sections[END_REF], [START_REF] Ourselin | Reconstructing a 3D structure from serial histological sections[END_REF], [START_REF] Schwier | Registration of histological whole slide images guided by vessel structures[END_REF], [START_REF] Stille | 3D reconstruction of 2D fluorescence histology images and registration with in vivo MR images: Application in a rodent stroke model[END_REF]. The problem becomes remarkably difficult if the order of the slides is not known. In such a case the process involves not only the transformation of the coordinate systems of the reference image and the input image into the joint coordinate system, but the most important problem is determine the reference images, i.e. the nearest neighbors (the preceding and the subsequent slide image) for each section image. It is clear that aligning every pair of slides is possible only when the number of images is small enough as it was assumed in our previous paper [START_REF] Górniak | Registration and sequencing of vessels section images at macroscopic levels[END_REF]. Otherwise, the problem of common aligning and sequencing of slide images becomes extremely expensive from computational complexity problem of view, especially when separate slide image consists of over a million pixels.

If the correct order of slides is known the reconstruction of a 3-D object body from a series of 2-D sections consists in the co-registration of any two consecutive slices. Nevertheless, the choice of the reference image has a significant impact on registration quality and 3-D reconstruction, even when the correct slide order is known, because the image chosen as a reference slide may contain distorted parts. Thus, it is advisable to repeat the registration in the reverse order and finally identify and correct possibly erroneous slides. When we do not know which 2-D slide should be the first (the last) in the reconstructed sequence, then the performed registration may lead to undesirable slides' deformations. Thus, it is important to find possibly correct image sequence before performing precise registration. It is clear that image matching could be realized on different levels of image precision. In the paper [START_REF] Arganda-Carreras | 3D Reconstruction of Histological Sections: Application to Mammary Gland Tissue[END_REF] the imaging system for simultaneous morphological and molecular analysis of section images is presented, but in that approach it is also assumed that correct order of slides is given at the start of the procedure. Another worth of attention method of the automatic registration of histological slices is given in [START_REF] Bagci | Fully Automatic 3D Reconstruction of Histological Images[END_REF]. The standardized slices are divided into small groups. The best reference slice selection in every group of slides is performed. The method of selection is developed based on an iterative assessment of image entropy and mean square error of the registration process. Nevertheless, in the paper it is also assumed that the order of slides presentation have not been destroyed.

Here we concentrate on providing a method of sequencing slides which is based on 2-D image features, independent of image translation and rotation (i.e., invariant to rigid transformations), but sensitive to a scale of an object on the slide. First, in order to avoid starting a work on a cellular level, we standardize images removing intensity variations between slices and after that binarize these images providing only contours of sections at the morphological level. We assume that microscopic slides are obtained from the geometrically continuous anatomical structure [START_REF] Ourselin | Reconstructing a 3D structure from serial histological sections[END_REF]. Thus, simple image measurements, such as an object diameter, object area, eccentricity of the best-fit ellipse containing the object, and the mean distance of all object's pixels from the center of gravity position of the binarized image, occur highly informative and allow us not only to distinguish one slide from another but also to determine geometrical similarity between them. We propose to consider these image measurements as coordinates in an abstract Euclidean space where each slide is a point with coordinates given by mentioned image measurements. In such a way the problem of slide sequencing can be formalized as a problem of finding the shortest open path between the points. We decided to use the sum of squared Euclidean distances as an optimization criterion. However, other goal functions can also be adequate.

The main reason that we prefer the sum of squared Euclidean distances criterion to the simple sum of Euclidean distances is that: the squared Euclidean distances favor the paths with more uniformly divided distances between separate points than in the case of the sum of strict Euclidean distances. Formally, we are looking for a path that visits every point once, and have separate (not-joined) endpoints, in contrast to the Hamiltonian cycle [START_REF] Lawler | The Traveling Salesman Problem: A Guided Tour of Combinatorial Optimization[END_REF], where the path is closed, i.e., forms a cycle. The numer of slides may vary between let say 10 to 200 (see for example [START_REF] Ourselin | Reconstructing a 3D structure from serial histological sections[END_REF], [START_REF] Stille | 3D reconstruction of 2D fluorescence histology images and registration with in vivo MR images: Application in a rodent stroke model[END_REF]). Thus, it is important to propose any efficient method of solving this in fact combinatorial problem which is known to be NP-hard [START_REF] Lawler | The Traveling Salesman Problem: A Guided Tour of Combinatorial Optimization[END_REF], [START_REF] Papadimitriou | The Euclidean traveling salesman problem is NP-complete[END_REF]. All method of solving the Travelling Salesmen Problem with Euclidean metric can be easily adapted to solving the open path problem. We decided to use modified version of Self-organizing Map (SOM) network with open chain topology [START_REF] Kohonen | Self-organizing map[END_REF], since it is fast, stable and provide very accurate solution when the number of points is not very huge. Furthermore, this method is more appropriate when the squared Euclidean distance instead of the Euclidean metric is used (see for example [START_REF] Budinich | A Self-Organising Neural Network for the Travelling Salesman Problem that is Competitive with Simulated Annealing[END_REF]).

Nevertheless, other general global optimization methods as simulated annealing, particle swarm optimization or genetic algorithms [START_REF] Kennedy | Particle Swarm Optimization Proceedings of IEEE International Conference on Neural Networks[END_REF], [START_REF] Kirkpatrick | Optimization by Simulated Annealing[END_REF], [START_REF] Saadatmand-Tarzjan | A Novel Constructive-Optimizer Neural Network for the Traveling Salesman Problem[END_REF] can be also applied for solving our sequencing problem.

It should be noted that the approach for slide images sequencing proposed here will be also helpful for detection of possibly missing slides, errors in the image order or even destroyed sections.

It is clear that the method proposed here generates image sequences optimal with respect to general shape and size parameters. When sections are very similar with respect to these parameters, other, more precise and taking into account internal section details should be used. Furthermore, in the case of simple artificial 3-D bodies, for example a ball, it is obvious that the ordering may be inadequate. Fortunately, biological objects usually characterize large shape variability, allowing the method to generate sequences which with high probability are closed to the true ones. Thus, a sequence obtained on the basis of morphological parameters should be additionally inspected and corrected using internal, possibly cellular structure of subsequent sections, but this problem is outside the scope of our paper. In this paper we do not concentrate on image processing tasks. We have used image processing procedures provided by Wolfram Mathematica 10.1. These procedures occurred to be adequate to our purposes.

The paper is organized as follow. In the next Section the problem of image representation that is invariant to location and rotation on the shape and the size level is addressed. Section 3 provides a concise description of image sequencing by determining the shortest open path through a given set of points in Euclidean feature space. Section 4 presents the details of the proposed modification of the Kohonen SOM algorithm for sequencing image feature points. Section 5 describes results of two computational experiments of sequencing a set of artificial section images of a mechanical device and a set of microscopic histological slices. Finally, in Section 6 some brief conclusions are presented.

Image Representation on the Shape and the Size Level

We have assumed that each object image is represented roughly by chosen image measurements which are invariant to rigid object image transformation, i.e., are invariant to translations and/or rotations. Thus, in order to obtain exact measurements, it is important to separate precisely an object image from its background. The methods used for this task depend on the slide quality and visual inspection of results is very important [START_REF] Gonzalez | Digital Image Processing[END_REF]. In our experiments we have used algorithms supplied by Mathematica 10.1. We skip details because in every case the segmentation procedure was different. For example, in our histological sections it was important to remove images of patches used for fixing every section. After image segmentation process the background on each slide was deleted. The whole image was binarized, i.e., every object pixel intensity was set to 1 whereas the rest of pixels was set to 0. On the basis of the object silhouette we can to compute simple object features such as:

-The object area A (areas of the two largest image segments are summarized), i.e., number of object's pixels normalized by the total number of pixels forming the whole image; -The mean distance of to the center of gravity M CD (mean distance of all object's elements from the center of gravity position) divided by number of pixels in the image column; -Maximal perimeter distance of the object M P D ( i.e., object diameter measured in pixels and divided by number of pixels in the image column); -The eccentricity of the object E, i.e., the eccentricity of the best-fit ellipse containing the object.

E = √ 1 -(b/a) 2
, where a and b are major and minor axes of the ellipse.

We assume that subsequent slides should be characterized by similar feature values. Each image is represented by a vector of its features. We expect that good image sequence can be obtained by minimizing the sum of squared errors (distances) between image feature vectors. The good sequence means the sequence close to the true one. It is clear that in general the problem of choosing image features providing such good sequences depends on properties of object and on a quality of images.

Image Sequencing by Determining the Shortest Open Path through a Given Set of Points in the Euclidean Space

Similarly as each city in the travelling-salesman problem with Euclidean distances between cities (ETSP) [START_REF] Lawler | The Traveling Salesman Problem: A Guided Tour of Combinatorial Optimization[END_REF] is represented by a point on the map (usually 2-D space) also in our problem every image is represented by a point in some compact subset of R d . We assume that this set is a unit cube

I d = [0, 1] × . . . × [0, 1].
The problem of determining the shortest open path through a given set of N points is very close to the well known problem of finding the shortest closed path through a given set of N points, i.e., ETSP. Let X = {x 1 , . . . x N } be a given set of point in I d . Our goal is to find a sequence of points X, i.e., a permutation of points indexes Π = (π(1), π(2) . . . , π(N )), such that

Q P (Π, X) = N ∑ i=2 ||x π(i) -x π(i-1) || 2 , ( 1
)
is minimized over all possible permutations. ||.|| denotes the Euclidean norm. Due to a symmetry property of the Euclidean distance between a pair of points (and obviously also the squared Euclidean metric), there exist at most 0.5N ! different solutions of the problem. Testing all 0.5N ! possible image sequences one can always find the shortest path connecting all N points, but the computational effort for this exhaustive search strategy, rises exponentially with N and rapidly becomes unmanageable.

For N = 10 it is 1814400 different solutions, for N = 20 the number of solutions rises to 1216451004088320000, and for N = 30 it is definitely unmanageable to perform the brute force search.

Usually in the case of ETSP only the sum of Euclidean distances, and not squared Euclidean distances, is minimized:

C T SP (Π, X) = N ∑ i=2 ||x π(i) -x π(i-1) || + ||x π(N ) -x π(1) ||.
(

Furthermore, it can be assumed that π(1) = 1, so that it is possible to find only 0.5(N -1)! different closed paths between points from set X.

It is known that the Euclidean TSP, similarly as a more general TSP, is NP-hard [START_REF] Papadimitriou | The Euclidean traveling salesman problem is NP-complete[END_REF]. There are known many exact and approximate methods of solving TSP and ETSP (see for example [START_REF] Saadatmand-Tarzjan | A Novel Constructive-Optimizer Neural Network for the Traveling Salesman Problem[END_REF] and the literature cited in there). Most of them can be adopt to solving min Π Q P (Π, X). Here we concentrate on using Kohonen SOM network [START_REF] Kohonen | Self-organizing map[END_REF] with open chain topology since similar type SOM networks with closed chain topology used for solving ETSP [START_REF] Angeniol | Selforganizing feature maps and the traveling salesman problem[END_REF], [START_REF] Favata | A study of the application of Kohonen-type neural networks to the traveling salesman problem[END_REF], [START_REF] Budinich | A Self-Organising Neural Network for the Travelling Salesman Problem that is Competitive with Simulated Annealing[END_REF], [START_REF] Bai | An Efficient Growing Ring SOM and Its Application to TSP[END_REF], [START_REF] Bai | An new self-organizing maps strategy for solving the traveling salesman problem[END_REF], [START_REF] Maire | Comparison of Neural Networks for Solving the Travelling Salesman Problem[END_REF], [START_REF] Xu | A Self-organizing Map Algorithm for the Traveling Salesman Problem[END_REF] are known to minimize the sum of squared Euclidean distances (SSD) rather than the sum of Euclidean distances between TSP points [START_REF] Budinich | A Neural Network for the Travelling Salesman Problem with a Well Behaved Energy Function[END_REF].

In our work we have used the algorithm presented in detail in the next Section.

Kohonen SOM Algorithm for Sequencing Data Points

Let M denote the number of neurons. It is assumed that M ≥ N [START_REF] Angeniol | Selforganizing feature maps and the traveling salesman problem[END_REF], [START_REF] Bai | An new self-organizing maps strategy for solving the traveling salesman problem[END_REF]. Recall that N is the number of data points in set X. Each neuron is represented by a point in R d labeled by a number j = 1, 2, . . . , M . Coordinates of such a point, let say labeled by j, form a vector W j ∈ I d , i.e., W j is the j-th neuron's weight. Each data point x i ∈ X can be uniquely represented by neuron j ⋆ such that

i ⋆ (x i ) = arg min j=1,...,M ||x i -W j || 2 . ( 3
)
Such neuron is often called "a winner". In the case of ties we will randomly pick exactly one winning neuron. As a consequence, each neuron j ∈ {1, . . . , M } may be "empty", may be closest to exactly one data point from X or it may happen that two or even more than two data points are connected to the same neuron [START_REF] Favata | A study of the application of Kohonen-type neural networks to the traveling salesman problem[END_REF]. Let A(j) denote the set of data points connected to neuron j. It is clear that

M ∪ j=1 A(j) = {1, 2, . . . , N }. ( 4
)
Furthermore, A(j), j = 1, . . . , M define at least partial order in the set of data points. In the case when there exist A(j) such that |A(j)| > 1, strict order can be obtain by random choice or using the lexicographical order defined in X. We will be represent A j by lists rather than sets, which allows us to retain ordering. Thus, we obtain a linear order of data points from X uniquely providing sequence Π. Computational complexity of that ordering is O(N M log M ) (we skip the dependence on d).

SOM Learning

Kohonen learning rule is applied for updating neuron weights (their positions in the data space). Denote by W j (n) the weight of j-th neuron in n-iteration of learning process. The process of learning is some kind of stochastic approximation process [START_REF] Ritter | Neural Computation and Self Organizing Maps: An Introduction[END_REF]. In every iteration only one data point from X is chosen at random and used for neuron position updating. Kohonen's method of learning is very efficient since not only a winning neuron selected according to (3) is updated but its neighbors in neuron's chain topology are also updated according to:

W j (n) = W j (n) + α(n)h(x i , j, n)(x i -W j (n)), j = 1, . . . , M (5)
where h(x i , j, n) is a neighborhood function taking values in [0, 1] and α(n) is a learning rate. α(n) is usually very small and decreases to zero with n. The common choice of the neighborhood function is the Gaussian function

h(x i , j, n) = exp (- |i ⋆ (x i ) -j| 2σ(n) 2 ), (6)
were σ(n) → 0 when n → ∞. There could be different methods of changing σ and α. Here we will assume that the learning is performed according to the following rules of learning parametres updating:

α(n) = 0.5 exp (-n/1000) (7)
and σ(n) = σ 0 0.01 n/1000 , σ 0 = 20, 10, 5.

However, it should be noted that in the case of ETSP also other formulas are used with success [START_REF] Bai | An Efficient Growing Ring SOM and Its Application to TSP[END_REF], [START_REF] Budinich | A Self-Organising Neural Network for the Travelling Salesman Problem that is Competitive with Simulated Annealing[END_REF], [START_REF] Xu | A Self-organizing Map Algorithm for the Traveling Salesman Problem[END_REF].

Algorithm

Step1. Initialization: Let n be a discrete learning time starting from 0. The algorithm inputs are points X = {x 1 , . . . x N } ⊂ I d . Choose M ≥ N , (for example, let M = 1.5N ). We initialize neuron's weights W at random using the uniform distribution on cube [w 1,min , w The main loop of the algorithm consists of the following steps:

Step 2. If n > n max go to Step End. Otherwise, compute parameters α(n) = 0.5 exp (-n/1000) and σ(n) = σ 0 0.01 n/1000 . Randomly draw point x from set of points X. Let say x = x i . Step 3. Find neuron i ⋆ (x i ) closest to x i ( 3). Update weights according to [START_REF] Bai | An new self-organizing maps strategy for solving the traveling salesman problem[END_REF].

Step 4. Set A j = {} j = 1, . . . , M . For i = 1, 2, . . . , N find j = i ⋆ (x i ) and add index i to the list A j .

Step 5. Obtain new sequence Π by concatenating lists A j , j = 1, . . . , M (empty lists are rejected).

Step 6. Compute Q(Π, X) according to (1). If Q(Π, X) ≤ Q opt set: Q opt = Q(Π, X), Π opt = Π and W opt = W . Set n = n + 1 and go to Step 2.
Step End Return the best results: Π opt and Q opt .

Experiments

In this section we demonstrate results of computational experiments with sequencing both artificial and natural 2-D image slides. The artificial set of slides is a part of provided by Mathematica image example "ExampleData/CTengine.tiff" (20 subsequent image sections from the set of 110 sections). The histological set of section images consists of 10 slides. In this case we do not know how the proper sequence should look like. The sample images of each set are shown in Fig. 1 and Fig. 6.

Sequencing 2-D Section Images of a Technical Device

We have 20 subsequent image sections at our disposal. Each section image is of the size 256 × 256 pixels. In this example the true sequence of slides is known. Thus, we have randomly changed the numbering of slides. The true sequence (after random permutations of section numbers) is [START_REF] Papadimitriou | The Euclidean traveling salesman problem is NP-complete[END_REF][START_REF] Arganda-Carreras | 3D Reconstruction of Histological Sections: Application to Mammary Gland Tissue[END_REF][START_REF] Bagci | Fully Automatic 3D Reconstruction of Histological Images[END_REF][START_REF] Flusser | Image registration methods: a survey[END_REF][START_REF] Favata | A study of the application of Kohonen-type neural networks to the traveling salesman problem[END_REF][START_REF] Maire | Comparison of Neural Networks for Solving the Travelling Salesman Problem[END_REF][START_REF] Kennedy | Particle Swarm Optimization Proceedings of IEEE International Conference on Neural Networks[END_REF][START_REF] Kohonen | Self-organizing map[END_REF][START_REF] Bai | An new self-organizing maps strategy for solving the traveling salesman problem[END_REF][START_REF] Angeniol | Selforganizing feature maps and the traveling salesman problem[END_REF][START_REF] Lippolis | Automatic registration of multi-modal microscopy images for integrative analysis of prostate tissue sections[END_REF][START_REF] Ourselin | Reconstructing a 3D structure from serial histological sections[END_REF][START_REF] Bai | An Efficient Growing Ring SOM and Its Application to TSP[END_REF][START_REF] Górniak | Registration and sequencing of vessels section images at macroscopic levels[END_REF][START_REF] Budinich | Sorting with Self-Organising Maps[END_REF][START_REF] Budinich | A Self-Organising Neural Network for the Travelling Salesman Problem that is Competitive with Simulated Annealing[END_REF][START_REF] Gonzalez | Digital Image Processing[END_REF][START_REF] Budinich | A Neural Network for the Travelling Salesman Problem with a Well Behaved Energy Function[END_REF][START_REF] Lawler | The Traveling Salesman Problem: A Guided Tour of Combinatorial Optimization[END_REF][START_REF] Kirkpatrick | Optimization by Simulated Annealing[END_REF]. In Figure 1 four subsequent slides with new numbers 13, 15, 5, 1, respectively, have been shown.

We have computed three image features: the area of the binarized object A, the maximal diameter distance of the object M P D and the mean distance to the center of gravity M CD (see Figure 2 -left panel). The optimal sequence Fig. 2. Image coordinates of the randomly rearranged sequence of the mechanical device slides (see Fig. 1): area of the binarized object A (thick), maximal diameter distance of the object M P D (dashed), mean distance to the center of gravity M CD (dotted) (left panel). Image coordinates of the original (untidy) sequence of microscopy slides (see Fig. 6): area of the binarized object (thick), mean distance to the center of gravity (dashed), eccentricity of the best-fit ellipse containing the object (dotted) (right panel)

was obtained using the proposed here algorithm (see section 4.2) as: [START_REF] Papadimitriou | The Euclidean traveling salesman problem is NP-complete[END_REF][START_REF] Arganda-Carreras | 3D Reconstruction of Histological Sections: Application to Mammary Gland Tissue[END_REF][START_REF] Bagci | Fully Automatic 3D Reconstruction of Histological Images[END_REF][START_REF] Flusser | Image registration methods: a survey[END_REF][START_REF] Favata | A study of the application of Kohonen-type neural networks to the traveling salesman problem[END_REF][START_REF] Maire | Comparison of Neural Networks for Solving the Travelling Salesman Problem[END_REF][START_REF] Kennedy | Particle Swarm Optimization Proceedings of IEEE International Conference on Neural Networks[END_REF][START_REF] Bai | An new self-organizing maps strategy for solving the traveling salesman problem[END_REF][START_REF] Angeniol | Selforganizing feature maps and the traveling salesman problem[END_REF][START_REF] Kohonen | Self-organizing map[END_REF][START_REF] Lippolis | Automatic registration of multi-modal microscopy images for integrative analysis of prostate tissue sections[END_REF][START_REF] Ourselin | Reconstructing a 3D structure from serial histological sections[END_REF][START_REF] Bai | An Efficient Growing Ring SOM and Its Application to TSP[END_REF][START_REF] Górniak | Registration and sequencing of vessels section images at macroscopic levels[END_REF][START_REF] Budinich | Sorting with Self-Organising Maps[END_REF][START_REF] Budinich | A Self-Organising Neural Network for the Travelling Salesman Problem that is Competitive with Simulated Annealing[END_REF][START_REF] Gonzalez | Digital Image Processing[END_REF][START_REF] Budinich | A Neural Network for the Travelling Salesman Problem with a Well Behaved Energy Function[END_REF][START_REF] Lawler | The Traveling Salesman Problem: A Guided Tour of Combinatorial Optimization[END_REF][START_REF] Kirkpatrick | Optimization by Simulated Annealing[END_REF]. This sequence differs from the true one in the position of slide 15. Its correct location is just before slides 5 and 1 (not after them). Process of Kohonen network learning was performed ten times starting from randomly chosen neuron's positions. The number of neurons was M = 20, 30, 40 and σ 0 = 20. In the case of M = 30 the optimal solution was obtained in every simulation. Otherwise, when the number of neurons was smaller or larger (i.e., M = 20 or M = 40), the network has produced also suboptimal paths. Trajectories of the sum of the squared Euclidean distances (SSD) obtained by the SOM algorithm in two exemplary runs are depicted on Figure 3. Figure 4 shows the sketch of the neuronal Fig. 3. The results of the section images sequencing using the proposed method. SSD obtained in the subsequent iterations of the algorithm for 20 neuron (left panel) and for 30 neuron (right panel). In the both cases the same -optimal sequence -was found paths corresponding to the obtained optimal sequence. The network equipped with M = 30 neurons visibly converges to the optimal solution, while in the case M = 20 the best sequence is computed and memorized in the Step 5 of the algorithm, but the network further converges to the slightly worse sequence. Finally, optimal with respect to SSD section sequence. The differences are not visible, although small inconsistencies can be found inside the object, since the location of section 15 should be corrected in the precise matching of neighboring slides.

Sequencing Microscopy Section Images

The analyzed set of microscopy section images contains 10 slides of the size 184 × 184 pixels (see Fig 6). This time also three image features have been computed: the area of the binarized object A, the eccentricity of the best-fit ellipse containing the 2-D object E, and the mean distance to the center of gravity M CD. Due to the vary small differences between measurements and relatively large differences between mean values of the separate features, values of A and M CD features rescaled by 0.15 and 0.2, respectively. The resacaled feature values are shown on the right panel of Fig. 2. The optimal with respect to the SSD section sequence was: (5, 6, 10, 4, 8, 7, 2, 3, 9, 1) with SSD = 0.0075 and the second best was (5, 6, 4, 10, 8, 7, 2, 3, 9, 1) with SSD = 0.0084. Proces of Kohonen network learning was performed ten times starting from randomly chosen neuron's positions. The number of neurons was M = 10 and M = 15 and σ 0 = 5 and σ 0 = 20, respectively. Trajectories of the sum of the SSD obtained by the SOM algorithm in two exemplary runs are depicted on Figure 7. In the both cases the same -optimal sequence was found, but only for 15 neurons the networks has converged to the optimal solution. The outline of the neuronal paths corresponding to the obtained optimal sequence is shown in Fig. 8 The network equipped with M = 15 neurons visibly converges to the optimal solution, while in the case M = 10 the best sequence is computed and memorized in the Step 5 of the algorithm, but the network further converges to the slightly worse sequence. 3-D reconstruction of an biological object obtained on the basis of the optimal (with respect to SSE) section sequence is presented in Fig 9.

Conclusion

The paper presents a problem of 2-D section images ordering for 3-D body reconstruction. We propose a slightly modified SOM Kohonen networks with the open chain topology for adequate sequence specification based on 2-D image features invariant to section image translation and rotation, but sensitive to a scale of an object on the slide. This means that location and orientation of objects on slides can be coincidental and sequencing of slides at the shape and the size level is performed without precise image registration. It should be noted that the approach proposed here for slide images sequencing will be also helpful for detection of possibly missing slides, errors in the image order or even destroyed sections. The method of slide sequencing proposed here provides image orders optimal with respect to its general shape and size parameters. When sections are very similar with respect to these parameters, other, more precise and taking into account internal section details should be used. We stress that a sequence obtained on the basis of morphological parameters should be additionally inspected and corrected using internal, possibly cellular features of the subsequent sections. This problem is outside the scope of the paper and will be developed later taking into account the non-rigid registration problems.

The presented SOM algorithm is shown to perform well on artificial and on natural series section images when the number of neurons is greater than the number of slides, but not too large. 1.5N neurons, where N is the number of section images was a best choice. It is less than in the case of ETSP, where 2N -3N neurons are recommended. Inspecting sequence path picture (Fig. 4,8) one can observe that the points representing images are located rather regularly, also forming almost regular routes. This can be reason that a relatively small number of neuron is here needed. Due to the lack of space we have shown only the small part of experiments we have performed, but it is worth to stress that the number of iterations needed for obtaining good solution was relatively stable, independent of the number of slides. This observation agrees with experiments performed by Budinich [START_REF] Budinich | Sorting with Self-Organising Maps[END_REF] (with respect to one dimensional ordering) and indicates that properties of SOM Kohonen networks seem to be worth of further investigations.

We also hope that the proposed approach can be also usable in others sequencing problems where object of complicated structure should be ordered.

Fig. 1 .

 1 Fig. 1. The part of the image set consisting of 20 section images of a mechanical device; the size of image 256 × 256 pixels

Fig. 4 .

 4 Fig. 4. The path of 20 neurons obtained for image sequence optimal with respect to SSD (left panel) and the path of 30 neurons obtained for image sequence optimal with respect to SSD (right panel) . Large dots indicate position of images

Figure 5

 5 Figure 5 presents comparison of the 3-D image of a part of the mechanical device (on the left image) and and its 3-D reconstruction obtained on the basis of the

Fig. 5 .

 5 Fig. 5. 3-D image of a part of an mechanical device (left panel) 3-D reconstruction of an mechanical device object obtained on the basis of the optimal with respect to SSD section images sequence (right panel)

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

 6789 Fig. 6. The image set consisting of 10 microscopy section images of the size 184 × 184 pixels

  1,max ] × . . . × [w d,min , w d,max ], where wk,min = M in 1≤i≤N [x i,k ], w k,max = M ax 1≤i≤N [x i,k ](k = 1, . . . , d), and (x i,1 , . . . , x i,k , . . . , x i,d ) T = x i , i = 1, . . . , N . Set W opt = W .Denote the value of the best sequence according to (1) by Q opt and the optimal sequence consisting of N element list by Π opt . Start with Q opt = 100 and Π opt = (1, 2, . . . , N ) Let the initial value of neighborhood function parametr be σ 0 = 20. And let n max be the value of the maximum number of iterations.

Acknowledgments. The authors express their special thanks to Dr. Agnieszka Malińska and to professor Maciej Zabel from University School of Medicine, Poznań, Poland for the microscopy section images.

This research was supported by S50242 grant at the Faculty of Electronics, Wroc law University of Science and Technology.


Alexandros Agapitos
email: alexandros.agapitos@ucd.ie

Andreas Konstantinidis

Haris Haralambous

Harris Papadopoulos
email: h.papadopoulos@frederick.ac.cy

Evolutionary Prediction of Total Electron Content over Cyprus

Keywords: Evolutionary Algorithms, Genetic Programming, Global Positioning System, Total Electron Content

Total Electron Content (TEC) is an ionospheric characteristic used to derive the signal delay imposed by the ionosphere on trans-ionospheric links and subsequently overwhelm its negative impact in accurate position determination. In this paper, an Evolutionary Algorithm (EA), and particularly a Genetic Programming (GP) based model is designed. The proposed model is based on the main factors that influence the variability of the predicted parameter on a diurnal, seasonal and long-term time-scale. Experimental results show that the GP-model, which is based on TEC measurements obtained over a period of 11 years, has produced a good approximation of the modeled parameter and can be implemented as a local model to account for the ionospheric imposed error in positioning. The GP-based approach performs better than the existing Neural Network-based approach in several cases.

Introduction

The ionosphere is defined as a region of the earth's upper atmosphere where sufficient ionisation can exist to affect the propagation of radio waves. It ranges in height above the surface of the earth from approximately 50 km to 1000 km. The influence of this region on radio waves is accredited to the presence of free electrons. The impact of the ionosphere on communication, navigation, positioning and surveillance systems is determined by variations in its electron density profile and total electron content along the signal propagation path [START_REF] Goodman | HF Communications, Science and Technology[END_REF][START_REF] Maslin | The HF Communications, a Systems Approach[END_REF]. As a result satellite systems for communication and navigation, surveillance and control that are based on transonospheric propagation may be affected by complex variations in the ionospheric structure in space and time leading to degradation of the accuracy, reliability and availability of their service. Total Electron Content (TEC) is an important parameter in trans-ionospheric links since when multiplied by a factor which is a function of the signal frequency, it yields an estimate of the delay imposed on the signal by the ionosphere due to its dispersive nature.

This paper describes an attempt to develop a model to predict TEC above Cyprus to encapsulate its variability on a diurnal, seasonal and long-term scale. The model development is based on around 60000 hourly TEC measurements recorded above Cyprus from 1998 to 2009. The practical application of this model lies in its possible use as an alternative candidate local model to the existing Klobuchar global model [START_REF] Klobuchar | Ionospheric Time-Delay Algorithm for Single-Frequency GPS Users[END_REF] that is currently being used in single frequency GPS navigation system receivers to improve positioning accuracy.

Metaheuristics and more specifically Evolutionary Algorithms were proven efficient and effective in dealing with difficult-to-solve real-life problems [START_REF] Reeves | Genetic algorithms, Handbook of Metaheuristics[END_REF]. Particularly, Genetic Programming (GP) based approaches performed well in evolving computer programs, controllers and models [START_REF] Koza | Genetic Programming: on the programming of computers by means of natural selection[END_REF] in the past. In this paper, we have adopted a Pareto-based Genetic Programming approach for dealing with the TEC problem. The proposed GP is a panmictic, generational, elitist genetic algorithm with an expression-tree representation [START_REF] Koza | Genetic Programming: on the programming of computers by means of natural selection[END_REF]. To the best of our knowledge this is the first time that a GP-based approach is applied to the proposed problem. The main contribution of our paper is: a GP-based prediction model designed for the TEC over Cyprus that outperforms the previously proposed Neural Network-based model [START_REF] Haralambous | A local Total Electron Content Neural Network model over Cyprus[END_REF]. The rest of this paper is structured as follows. Section 2 describes the measurements and characteristics of TEC values used for the model and section 3 discusses the model parameters involved in the model development. Section 4 gives a description of the proposed evolutionary algorithm and section 5 outlines the experimental results. Finally, section 6 gives the concluding remarks of the paper.

Measurements and characteristics of TEC

Dual-frequency GPS data recorded by GPS receivers enable an estimation of the TEC measured in total electron content units, (1 TECU = 10 16 elecrons m -2 ). This is the total amount of electrons along a line of sight between the receiver and a GPS satellite in a column of 1 m 2 cross-sectional area and represents a parameter of interest to GPS users (see Figure 1a). TEC is the integral of the electron density (Ne see Figure 1b) from the ground to an infinite height (the height of the satellite). The electron density of free electrons within the ionosphere and therefore TEC depend upon the strength of the solar ionizing radiation which is a function of time of day, season, geographical location and solar activity [START_REF] Goodman | HF Communications, Science and Technology[END_REF], [START_REF] Maslin | The HF Communications, a Systems Approach[END_REF]. Since solar activity has an impact on ionospheric dynamics which in turn influence the electron density of the ionosphere, TEC also exhibits variability on daily, seasonal and long-term time scales in response to the effect of solar radiation. It is also subject to abrupt variations due to enhancements of geomagnetic activity following extreme manifestations of solar activity disturbing the ionosphere from minutes to days on a local or global scale. The most profound solar effect on TEC is reflected on its daily variation as shown in Figure 2. As it is clearly depicted, there is a strong dependency of TEC on local time which follows a sharp increase of TEC around sunrise and gradual decrease around sunset. This is attributed to the rapid increase in the production of electrons due to the photo-ionization process during the day and a more gradual decrease due to the recombination of ions and electrons during the night. The long-term effect of solar activity on TEC is also clearly shown as we can observe a correlation of the mean level of TEC and sunspot number which is an established index of solar activity.

There is also a seasonal component in the variability of TEC which can be attributed to the seasonal change in extreme ultraviolet (EUV) radiation from the Sun. This can be clearly identified in Figure 3 for noon values of TEC for high and low solar activity periods (years 2001 and 2008).

Model Parameters

The diurnal variation of TEC is clearly evident by observing

       
(1) and cos 2 24

hour coshour      (2)
A seasonal variation is also an underlying characteristic of TEC as shown in fig. 3 and is described by day number daynum in the range 1 ≤ daynum ≤ 365. Again to avoid unrealistic discontinuity between December 31st and January 1st daynum is converted into its quadrature components according to: sin 2 365

daynum sinday         (3)
and cos 2 365

daynum cosday         (4)
Long-term solar activity has a prominent effect on TEC. To include this effect in the model specification we need to incorporate an index, which represents a good indicator of solar activity. In ionospheric work the 12-month smoothed sunspot number is usually used, yet this has the disadvantage that the most recent value available corresponds to TEC measurements made six months ago. To enable TEC data to be modelled as soon as they are measured, and for future predictions of TEC to be made, the monthly mean sunspot number values were modeled using a smooth curve defined by a summation of sinusoids.

The proposed Evolutionary Algorithm

Genetic Programming (GP) is an Evolutionary Computation (EC) technique that evolves populations of computer programs as solutions to problems. It is a general purpose evolutionary search technique that can be applied in both regression and classification problems. In contrast to linear and non-linear regression, the technique does not presuppose a functional form as it generally proceeds for a definition of lower-level building blocks (the function set). Therefore Koza [START_REF] Koza | Genetic Programming: on the programming of computers by means of natural selection[END_REF] speaks of symbolic regression or system identification as the object of search is both the functional form and the optimal coefficients.

The term evolutionary algorithm describes a class of stochastic search procedure inspired by principles of natural genetics and survival of the fittest. They operate through a simulated evolution process on a population of solution structures that represent candidate solutions in the search space. Evolution occurs through (1) a selection mechanism that implements a survival of the fittest strategy, and (2) diversification of the selected solutions to produce offspring for the next generation. In GP, programs are usually expressed using hierarchical representations taking the form of syntax-trees. It is common to evolve programs into a constrained, and often problem-specific user-defined language. The variables and constants in the program are leaves in the tree (collectively named as terminal set), whilst arithmetic operators are internal nodes (collectively named as function set). It is common in the GP literature to represent expressions in the prefix notation similar to that used in LISP or Scheme. For example, x+3*x becomes (+ x (* 3 y))). This representation eases the expression-tree data structure formation, and its manipulation during the application of variation operators.

GP finds out how well a program works by running it, and then comparing its behaviour to some ideal. We might be interested, for example, in how well a program predicts a time series or controls an industrial process. This comparison is quantified to give a numeric value called fitness. Those programs that do well are chosen to breed, and produce new programs for the new generation. The primary variation operators to perform transitions within the space of computer programs are crossover and mutation. Once a stopping criterion has been met the algorithm terminates and the best program is designated as the output of the run.

The most commonly used form of crossover is subtree crossover [START_REF] Koza | Genetic Programming: on the programming of computers by means of natural selection[END_REF]. Given two parents, subtree crossover randomly (and independently) selects a crossover point (a node) in each parent tree. Then, it creates the offspring by replacing the subtree rooted at the crossover point in a copy of the first parent with a copy of the subtree rooted at the crossover point in the second parent. Copies are used to avoid disrupting the original individuals. This way, if selected multiple times, they can take part in the creation of multiple offspring programs. Note that it is also possible to define a version of crossover that returns two offspring, but this is not commonly used. Often crossover points are not selected with uniform probability. Typical GP primitive sets lead to trees with an average branching factor of at least two, so the majority of the nodes will be leaves. Consequently the uniform selection of crossover points leads to crossover operations frequently exchanging only very small amounts of genetic material (i.e., small subtrees); many crossovers may in fact reduce to simply swapping two leaves. To counter this, Koza [START_REF] Koza | Genetic Programming: on the programming of computers by means of natural selection[END_REF] suggested the widely used approach of choosing functions 90% of the time and leaves 10% of the time.

The most commonly used form of mutation in GP (which we will call subtree mutation) randomly selects a mutation point in a tree and substitutes the subtree rooted there with a randomly generated subtree. Another common form of mutation is point mutation, which is GP's rough equivalent of the bit-flip mutation used in genetic algorithms. In point mutation, a random node is selected and the primitive stored there is replaced with a different random primitive of the same arity taken from the primitive set. When subtree mutation is applied, this involves the modification of exactly one subtree. Point mutation, on the other hand, is typically applied on a pernode basis. That is, each node is considered in turn and, with a certain probability, it is altered as explained above. This allows multiple nodes to be mutated independently in one application of point mutation.

Like in other evolutionary algorithms, in GP the individuals in the initial population are typically randomly generated. Two dominant methods are the full and grow and the widely used combination of the two known as Ramped half-and-half [START_REF] Koza | Genetic Programming: on the programming of computers by means of natural selection[END_REF]. In both the full and grow methods, the initial individuals are generated so that they do not exceed a user-specified maximum depth. The depth of a node is the number of edges that need to be traversed to reach the node starting from the tree's root node (the depth of the tree is the depth of its deepest leaf). The full method generates full tree-structures where all the leaves are at the same depth, whereas the grow method allows for the creation of trees of more varied sizes and shapes.

The evolutionary algorithm employed is a standard elitist (i.e. the best is always preserved), generational (i.e. populations are arranged in generations, not steadystate), panmictic (i.e. no program mating restrictions) genetic algorithm [START_REF] Goldberg | Genetic Algorithms in Search Optimisation and Machine Learning[END_REF]. The algorithm uses tournament selection with a tournament size of 7. Evolution proceeds for 50 generations, and the population size is set to 1000 individuals. Ramped-halfand-half tree creation with a maximum depth of 6 is used to perform a random sampling of rules during run initialisation. Throughout evolution, expression-trees are allowed to grow up to depth of 12. The evolutionary search employs a mixture of mutation-based variation operators, where subtree mutation is combined with pointmutation; a probability governing the application of each, set to 0.6 in favour of subtree mutation. Neither recombination, nor reproduction was used. The primitive language consisted of the basic arithmetic operators (+, -, *, /) serving as the function set, whereas the terminal set consisted of the five independent variables sinhour, coshour, sinday, cosday and modeled monthly mean sunspot number.

Experimental Results and Discussion

The primary goal of our experimental studies is to investigate the performance of our GP-based approach in designing a prediction model for the TEC over Cyprus with a good approximation to the measured values, compared to the previously developed Neural Network based model. The NN used in [START_REF] Haralambous | A local Total Electron Content Neural Network model over Cyprus[END_REF] had a fully connected two-layer structure, with 5 input, 10 hidden and 1 output neurons. Both its hidden and output neurons consisted of hyperbolic tangent sigmoid activation functions. The number of hidden neurons was determined by trial and error. The training algorithm used was the Levenberg-Marquardt back propagation algorithm.

The data-set was segmented in 10 continuous folds similarly to the case of NNs [START_REF] Haralambous | A local Total Electron Content Neural Network model over Cyprus[END_REF]. In each cross-validation cycle, 9 folds were used as the training set, whereas the evolved model was tested on the remaining 10 th fold. The training set was further randomly divided into two data-sets (with no overlapping): the fitness evaluation data-set, with 67% of the training data, and the validation data-set with the remaining 33%. The fitness measure consists of minimising the RMSE on the fitness evaluation data-set. At each generation, a two-objective sort is conducted in order to extract a set of non-dominated individuals [START_REF] Deb | Multi-Objective Optimization Using Evolutionary Algorithms[END_REF] (Pareto front) with regards to the lowest fitness evaluation data-set RMSE, and the smallest model complexity in terms of expressiontree size, as measured by the number of tree-nodes [START_REF] Gagne | Genetic Programming, Validation Sets, and Parsimony Pressure[END_REF]. The rationale behind this is to create a selection pressure towards simpler prediction models that have the potential to generalise better. These non-dominated individuals are then evaluated on the validation data-set, with the best-of-generation prediction model selected as the one of these with the smallest RMSE. During tournament selection based on the fitness evaluation data-set performance, we used the model complexity as a second point of comparison in cases of identical error rates. The two approaches were coded in Java and run on an Intel® Pentium 4 3.2 GHz Windows XP server with 1.5 GB RAM. We performed 50 independent evolutionary runs for each test fold, in order to account for the stochastic nature of the adaptive search algorithm, and obtain statistically meaningful results. Figure 4 shows a comparison between the GP and the NN approaches in terms of min RMSE per fold and average RMSE, respectively. The results in Figure 4 show that the GP model provides a lower RMSE compared to the NN-based model in 6 out of 10 folds, giving a better prediction and consequently a better approximation in around 60% of the TEC measured values. However, the NN-based model outperforms the GP-model in folds 2, 3, 4 and 5. The average RMSE obtained by the GP approach is around 4% less than the average RMSE obtained by the NN approach. The two approaches have a similar standard deviation of around 4 RMSE. Some examples of measured and predicted TEC values are given in Figure 5. These demonstrate both the good performance of the developed GP and its superiority over the NN model.

Conclusions and Future Work

In this paper, a Genetic Programming based approach is used to design a prediction model for the Total Electron Content over Cyprus. Particularly, a panmictic, generational, elitist genetic algorithm with an expression-tree representation is used. A prediction model is developed based on a data set obtained during a period of eleven years covering a full sunspot cycle. The GP-model has shown a good approximation of the different time-scales in the variability of the modelled parameter and it has outperformed the existing Neural Network based model. The proposed model can therefore be used in single frequency GPS navigation system receivers to account for the ionospheric imposed error in positioning.

There are a number of avenues for future research. For example, it will be interesting to investigate different genetic operators and primitive languages to further improve the performance of the GP approach. Moreover, the hybridization of the GP with NNs and the design of a more robust approach is also a future possibility.

Fig. 1 .

 1 Fig. 1. TEC representation and typical electron density profile of the ionosphere over Cyprus.

Fig. 2 .

 2 Fig. 2. Diurnal variability of TEC for low, medium and high solar activity.

Fig. 3 .

 3 Fig. 3. Seasonal variation of TEC at 12:00.

Figure 2 .

 2 We therefore include hour number as an input to the model. The hour number, hour, is an integer in the range 0 ≤ hour ≤ 23. In order to avoid unrealistic discontinuity at the midnight boundary, hour is converted into its quadrature components according to: sin 2 24 hour sinhour

Fig. 4 .

 4 Fig. 4. GP versus NN in terms of Root Mean Square Error (RMSE) per fold and average Root Mean Square Error (RMSE).

Fig. 5 .

 5 Fig. 5. Examples of measured and predicted TEC values.


Tomáš Kocyan
email: tomas.kocyan@vsb.cz

Kateřina Slaninová
email: katerina.slaninova@vsb.cz

Jan Martinovič
email: jan.martinovic@vsb.cz

Flexible Global Constraint Extension for Dynamic Time Warping

Keywords: dynamic time warping, flexible global constraint, longest common subsequence, comparison of sequences

Dynamic Time Warping algorithm (DTW) is an effective tool for comparing two sequences which are subject to some kind of distortion. Unlike the standard methods for comparison, it is able to deal with a different length of compared sequences or with reasonable amount of inaccuracy. For this reason, DTW has become very popular and it is widely used in many domains. One of its the biggest advantages is a possibility to specify definable amount of benevolence while evaluating similarity of two sequences. It enables to percept similarity through the eyes of domain expert, in contrast with a strict sequential comparison of opposite sequence elements. Unfortunately, such commonly used definition of benevolence cannot be applied on DTW modifications, which were created for solving specific tasks (e.g. searching the longest common subsequence). The main goal of this paper is to eliminate weaknesses of commonly used approach and to propose a new flexible mechanism for definition of benevolence applicable to modifications of original DTW.

Introduction

Nowadays, searching and comparing time series databases generated by computers, which consist of accurate time cycles and which achieve a determined finite number of value levels, is a trivial problem. Main attention is focused rather on optimization of the searching speed. A non-trivial task occurs while comparing or searching signals with different length, which are not strictly defined and have various distortions in time and amplitude. As a typical example, we can mention the measurement of functionality of human body (ECG, EEG) or the elements (precipitation, flow rates in riverbeds), that does not contain any accurate timing for signal generation. Therefore, comparison of such sequences is significantly difficult, and almost impossible while using standard functions for similarity (distance) computation [START_REF] Ding | Querying and mining of time series data: Experimental comparison of representations and distance measures[END_REF], such as Euclidean distance [START_REF] Elmore | Euclidean distance as a similarity metric for principal component analysis[END_REF], cosine measure [START_REF] Lee | Document ranking and the vector-space model[END_REF], Mean Estimate Error [START_REF] Singh | Hydrological modeling of the iroquois river watershed using hspf and swat[END_REF], etc. Examples of such signals are presented in Figure 1. A problem of standard functions for similarity (distance) computation consists in sequential comparison of the opposite elements in the both sequences (comparison of elements with the identical indices). Fortunately, such lack of commonly used approach can be easily eliminated by the Dynamic Time Warping algorithm, which is able to percept similarity through the eyes of a domain expert, in contrast with a strict sequential comparison. However, such commonly used definition of benevolence cannot be applied on DTW modifications, which were created for solving specific tasks (e.g. searching the longest common subsequence).

The main goal of this paper is to eliminate weaknesses of commonly used approach and to propose a new flexible mechanism for definition of benevolence applicable to modifications of the original DTW. It is organized as follows: First, the DTW algorithm for comparing two distorted sequences and its several modifications will be described in Section 2. In Section 3, commonly used approaches for definition of benevolence will be introduced. It will be followed by a proposal of a new Flexible Global Constraint. Finally, an effect of the algorithm's settings will be visualized and the proposed solution will be discussed.

Dynamic Time Warping

Dynamic Time Warping (DTW) is a technique for finding the optimal matching of two warped sequences using pre-defined rules [START_REF] Müller | Information Retrieval for Music and Motion[END_REF]. Essentially, it is a nonlinear mapping of particular elements to match them in the most appropriate way. The output of such DTW mapping of sequences from Figure 1 can be seen in Figure 2. At first, this approach was used for comparison of two voice patterns during an automatic recognition of voice commands [START_REF] Rabiner | Fundamentals of Speech Recognition[END_REF]. Since this time, it was widely used in many domains, e.g. for efficient satellite image analysis [START_REF] Petitjean | Efficient satellite image time series analysis under time warping[END_REF], in analysis of student behavioral patterns [START_REF] Slaninová | Dynamic time warping in analysis of student behavioral patterns[END_REF] or in protein fold recognition [START_REF] Lyons | Protein fold recognition by alignment of amino acid residues using kernelized dynamic time warping[END_REF]. As it is correctly noted in [START_REF] Keogh | Scaling up dynamic time warping for datamining applications[END_REF], a common problem of many DTW applications lies in the fact, that the DTW is too computationally expensive. In order to speed up the algorithm run, several lower bounding methods [START_REF] Keogh | Exact indexing of dynamic time warping[END_REF] or parallelization techniques were created [START_REF] Rakthanmanon | Addressing big data time series: Mining trillions of time series subsequences under dynamic time warping[END_REF][START_REF] Sart | Accelerating dynamic time warping subsequence search with gpus and fpgas[END_REF]. Moreover, the DTW was modified many times for solving specific tasks (e.g. searching the longest common subsequence [START_REF] Kocyan | Searching the longest common subsequences in distorted data[END_REF]) or for better algorithm behavior (e.g. Derivative Dynamic Time Warping [START_REF] Keogh | Derivative dynamic time warping[END_REF]). Since the proposed approach is also an extension of this algorithm, the original DTW algorithm will be described in more detail for better understanding.

Formally, the main goal of DTW method is a comparison of two time dependent sequences x and y, where x = (x 1 , x 2 , . . . , x n ) and y = (y 1 , y 2 , . . . , y m ), and finding an optimal mapping of their elements. To compare partial elements of sequences x i , y j ∈ R, it is necessary to define a local cost measure c : R×R → R ≥0 , where c is small if x and y is similar to each other, and otherwise it is large. Computation of the local cost measure for each pair of elements of sequences Then the goal is to find an alignment between x and y with a minimal overall cost. Such optimal alignment leads through the black valleys of the cost matrix C, trying to avoid the white areas with a high cost. Such alignment is demonstrated in Figure 3(b). Basically, the alignment (called warping path) p = (p 1 , . . . , p q ) is a sequence of q pairs (warping path points) p k = (p kx , p ky ) ∈ {1, . . . , n} × {1, . . . , m}. Each of such pairs (i, j) indicates an alignment between the ith element of the sequence x and jth element of the sequence y.

Retrieval of optimal path p * by evaluating all possible warping paths between sequences x and y leads to an exponential computational complexity. Fortunately, there exists a better way with O(n • m) complexity based on dynamic programming. It involves the use of an accumulated cost matrix D ∈ R n×m described in [START_REF] Müller | Information Retrieval for Music and Motion[END_REF].

Accumulated cost matrix computed for the cost matrix from Figure 3(a) can be seen in Figure 4(a). It is evident that the accumulation highlights only a single black valley. The optimal path p * = (p 1 , . . . , p q ) is then computed in a reverse order starting with p q = (n, m) and finishing in p 1 = (1, 1). An example of such found warping path can be seen in Figure 4 The final DTW cost can be understood as a quantified effort for the alignment of the two sequences (see Equation 1).

DT W (x, y) = q k=1 C(x p kx , y p ky ) = D(n, m) (1)

Subsequence DTW

In some cases, it is not necessary to compare or align the whole sequences. A usual goal is to find an optimal alignment of a sample (a relatively short time series) within the signal database (a very long time series). This is very usual in situations, in which one manages with a signal database and wants to find the best occurrence(s) of a sample (query). Using the slight modification [START_REF] Müller | Information Retrieval for Music and Motion[END_REF], the DTW has the ability to search such queries in a much longer sequence. The basic idea is not to penalize the omission in the alignment between x and y that appears at the beginning and at the end of the sequence y. Suppose we have two sequences x = (x 1 , x 2 , . . . , x n ) of the length n ∈ N and y = (y 1 , y 2 , . . . , y m ) of the much larger length m ∈ N. The goal is to find a subsequence y a:b = (y a , y a+1 , . . . , y b ) where 1 ≤ a ≤ b ≤ m that minimizes the DTW cost to x over the all possible subsequences of y. An example of such searching the best subsequence alignment can be seen in Figure 5. Both constructed matrices including the found warping path are then shown in Figure 6. Despite the fact that the DTW has its own modification for searching subsequences, it works perfectly only in case of searching an exact pattern in a signal database. However, in real situations, exact patterns are not available because they are surrounded by additional values, or even repeated several times in a sequence (see Figure 7). Unfortunately, the basic DTW is not able to handle these situations and it fails or returns only a single occurrence of the pattern. To deal with this type of situations, several DTW modifications were created and described for example in [START_REF] Kocyan | Searching the longest common subsequences in distorted data[END_REF] or [START_REF] Movchan | Time series subsequence similarity search under dynamic time warping distance on the intel many-core accelerators[END_REF] in detail. The biggest difference is in the approach for searching the warping path. In simple terms, the algorithm does not search the warping path from the upper right corner to the bottom left one (shown in the case of classical DTW in Figure 8(a)) and also it does not connect the opposite sides of the matrix (shown in the case of subsequence DTW in Figure 8(b)). The main idea is to find warping paths as long as possible from any element to another one, parallel to a diagonal, as it is outlined in Figure 8(c). An example of such found common subsequences can be seen in Figure 10. The corresponding warping paths are also visualized in the cost matrix in Figure 9.

Flexible Global Constraints

In the practical applications [START_REF] Toyoda | Discovery of cross-similarity in data streams[END_REF][START_REF] Cheng | An image-to-class dynamic time warping approach for both 3d static and trajectory hand gesture recognition[END_REF][START_REF] Xu | Automated detection of burned-out luminaries using indoor positioning[END_REF][START_REF] Zhao | Adaptive fuzzy clustering based anomaly data detection in energy system of steel industry[END_REF], the construction of a warping path has to be controlled. The reason is possible uncontrolled high number of warpings, i.e. alignment of a single element to a high number of the elements in the opposite sequence [START_REF] Müller | Information Retrieval for Music and Motion[END_REF]. In this manner, dissimilar sequences can get low DT W Cost and they can be evaluated as similar. This situation is demonstrated on sequences in Figure 11, and on appropriate cost matrix in Figure 12.

Generally, this can be easily fixed by definition of a global constraint region R ⊆ D. This region then determines the elements of the cost matrix, which can be used for searching the warping path. In the original paper about DTW [START_REF] Müller | Information Retrieval for Music and Motion[END_REF], there are two global constraints for warping path mentioned -Itakura parallelogram (Figure 13 However, for purpose of searching subsequences and other DTW modifications, the Itakura parallelogram seems to be inappropriate, because it was designed to limit warpings at the start and end of the classical DTW warping path, where the first and last warping points are exactly known. Fortunately, the Saoke-Chiba band looks more preferable. The warping path respecting this band for sequences from Figure 11 is visible in Figure 14. However, one may ask what width of band to choose. The width essentially defines the maximal number of warpings in a found sequence. For this reason, it is almost impossible to define a universal number applicable both on shorter and longer sequences. It is evident that allowing five warpings on a path comparing sequences of the length ten or hundred has absolutely different meaning. In this example, the results look satisfactorily, but this belt was also designed for searching the warping path through the whole sequences. This inaccuracy is evident in the following example:

Lets have two sequences x = (x 1 , x 2 , . . . , x n ) and y = (y 1 , y 2 , . . . , y 2n ), where y is created by stretching x into the double length (i.e. ∀i ∈ {1, . . . , 2n} : y i = x i/2 ). The matrix will stretch in one dimension and the line of minima will slightly bend (see Figure 15(a)). It causes some warpings, but it is still acceptable. Using the standard Sakoe-Chiba band, the warping path cannot follow the minima trajectory and have to continue in straight direction, as shown in Figure 15(b).

More elegant solution is to allow a band to bend itself and provide a warping path with reasonable freedom. For this purpose, we designed a flexible band allowing configurable bending. The band is based on Saoke-Chiba band, but it

Proposed modifications to Saoke-Chiba band

In our modification, the center of the band varies and passes through one of the previous points of the currently constructed warping path, called control point. Such control point is always located in the fixed distance from the currently processed point. This distance is called control point distance and it is defined as a number of warping path points preceding the currently processed point. The center of constructed band always moves to a newly established control point.

Formally, suppose we have a currently constructed warping path p defined as p = (p 1 , . . . , p q ) consisting of a sequence of q path points p k = (p kx , p ky ) ∈ {1, . . . , n} × {1, . . . , m}, p Figure 16 demonstrates how the increasing distance of the control point d causes higher toughness of the band, and how the ability to bend loses. The shorter distance makes the band more flexible, the higher distance causes inflexibility. It is especially evident from Figure 16(d) (with d = 4), where the band became too much tough to follow the black valleys.

An effect of predefined toughness can be also easily quantified by the received DT W Cost defined in Equation 1. With an original Saoke-Chiba Band (see Figure 13), received DT W Cost = 3.6433. On the other hand, with using the proposed flexible constraint (distances of the control point d) and appropriately adjusted benevolence, the sequences can be evaluated as almost equal (DT W Cost = 0, 0182). Table 1 illustrates how the received DT W Cost reflects the adjusted amount of benevolence (various distances of the control point d).

In order to set the control point distance up correctly, it is necessary to have some domain knowledge. At this point, the domain expert has to define the benevolence for the evaluation.

Conclusion

The Dynamic Time Warping algorithm has become widely used technique for comparing two sequences and evaluating their mutual similarity. Its many modifications, created for solving specific tasks, subsequently requested additional adjustments of partial steps of this algorithm. As a typical example, the DTW approach for searching the longest common subsequence can be mentioned. In this type of modification, none of commonly used constraints for construction of the warping path can be used. Therefore, the mail goal of this paper was to provide a solution for such situations and to propose a new flexible mechanism for definition of the constraint applicable to the modifications of the original DTW. The proposed solution consists in a new flexible constraint, which is based on the original Saoke-Chiba band. The constraint enables the control over the process of warping path construction and it generally offers more flexibility and predictable behaviour. Moreover, definition of its conduct (i.e. rigidity of the band) can be defined by a single number, which is not dependent on the length of the processed sequences. The use of the proposed solution is not limited only for searching the common subsequences, but it can be utilized in all DTW modifications, whose constructed warping paths are not defined by exactly beginnings and ends.

Fig. 1 .

 1 Fig. 1. Standard Metrics Comparison

Fig. 3 .

 3 Fig. 3. DTW Cost Matrices

  (b).

Fig. 4 .

 4 Fig. 4. DTW Accumulated Cost Matrices

Fig. 5 .Fig. 6 .

 56 Fig. 5. Found DTW Subsequence

Fig. 7 .Fig. 8 .

 78 Fig. 7. Basic DTW Subsequence Inaccuracies

Fig. 9 .Fig. 10 .

 910 Fig. 9. Cost Matrix with Found Warping Paths

  (a)) and Sakoe-Chiba band (Figure 13(b)).

Fig. 11 .Fig. 12 .Fig. 13 .

 111213 Fig. 11. Mapping of Dissimilar Sequences

Fig. 14 .

 14 Fig. 14. Examples of Applied Saoke-Chiba Bands

Fig. 15 .

 15 Fig. 15. Cost Matrices for Stretched Sequence

  1 = (n, m). Each such pair (p kx , p ky ) indicates an alignment between the ith element of the sequence x and jth element of the sequence y. The path point (p kx , p ky ) lies in the Saoke-Chiba band of a width w, if |p kx -p ky | < w. With the flexible band of the width w and with a control point distance d, the path point (p kx , p ky ) lies in the band if |(p kx -p (k-d)x ) -(p ky -p (k-d)y )| < w. The distance d of such control point from the end of the warping path defines a rigidity of the band.

6 Fig. 16 .

 616 Fig. 16. Various Distances of Control Point

Table 1 .

 1 DTW Costs for Various Control Point Distances

	Control point distance DT W Cost
	1	0.0182
	2	0.0182
	3	0.0182
	4	1.0329
	5	2.3110
	6	2.9819
	7	3.1504
	8	3.9598
	9	3.6433
	10	3.6433
	Without Flexible Band	3.9634

Acknowledgment

This work was supported by The Ministry of Education, Youth and Sports from the National Programme of Sustainability (NPU II) project 'IT4Innovations excellence in science -LQ1602'.