Kohonen SOM for Image Slides Sequencing

Abstract

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 mostly on microscope images obtained from serial sections of biological tissues. For adequate 3-D body reconstruction the correct sequence of section should be available. We propose a modified SOM Kohonen networks with open chain topology for true sequence estimation based on based on 2-D image features invariant to section image translation and rotation, but sensitive to a scale of an object on the slide. We provide an experimental validation of the proposed method using artificial section images of a mechanical device and real microscopic histological section images.

1 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 [10]. 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 alignment of successive anatomical slices [10]. We concentrate on the image registration in the case of 3-D tissue reconstruction [17, 19, 23, 24]. 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 [12]. 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 [2] 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 [3]. 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 [19]. 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 [16], where the path is closed, i.e., forms a cycle. The number of slides may vary between let say 10 to 200 (see for example [19, 24]). Thus, it is important to propose any efficient method of solving this in fact combinatorial problem which is known to be NP-hard [16, 20]. 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 [15], 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 [6]).

Nevertheless, other general global optimization methods as simulated annealing, particle swarm optimization or genetic algorithms [13, 14, 22] 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 Sect. 6 some brief conclusions are presented.

2 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 [11]. 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 MCD (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 MPD (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= \sqrt{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.

3 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) [16] 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] \times \ldots \times [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}, \ldots 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 $\varPi =(\pi (1), \pi (2)\ldots ,\pi (N))$, such that

$$\begin{aligned} Q_{P}(\varPi ,X)= \sum _{i=2}^{N} ||x_{\pi (i)}-x_{\pi (i-1)}||^2, \end{aligned}$$

(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 1216451 004088 320000, 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:

$$\begin{aligned} C_{TSP}(\varPi , X)=\sum _{i=2}^{N} ||x_{\pi (i)}-x_{\pi (i-1)}|| + ||x_{\pi (N)}-x_{\pi (1)}||. \end{aligned}$$

(2)

Furthermore, it can be assumed that $\pi (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 [20]. There are known many exact and approximate methods of solving TSP and ETSP (see for example [22] and the literature cited in there). Most of them can be adopt to solving $\min _{\varPi } Q_{P}(\varPi ,X)$. Here we concentrate on using Kohonen SOM network [15] with open chain topology since similar type SOM networks with closed chain topology used for solving ETSP [1, 4–6, 9, 18, 25] are known to minimize the sum of squared Euclidean distances (SSD) rather than the sum of Euclidean distances between TSP points [8].

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

4 Kohonen SOM Algorithm for Sequencing Data Points

Let M denote the number of neurons. It is assumed that $M \ge N$ [1, 5]. 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, \ldots , M$. Coordinates of such a point, let say labeled by j, form a vector $W_{j} \in I^{d}$, i.e., $W_{j}$ is the j-th neuron’s weight. Each data point $x_{i} \in X$ can be uniquely represented by neuron $j^{\star }$ such that

$$\begin{aligned} i^{\star }(x_{i}) = arg \min _{j=1,\ldots ,M} ||x_{i}-W_{j}||^2. \end{aligned}$$

(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 \in \{1, \dots , 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 [9]. Let A(j) denote the set of data points connected to neuron j. It is clear that

$$\begin{aligned} \bigcup _{j=1}^{M} A(j) = \{1,2, \ldots , N\}. \end{aligned}$$

(4)

Furthermore, $A(j), \; j=1, \dots , 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 $\varPi$. Computational complexity of that ordering is $O(NM\log {M})$ (we skip the dependence on d).

4.1 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 [21]. 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:

$$\begin{aligned} W_{j}(n) = W_{j}(n) + \alpha (n) h(x_{i}, j, n) (x_{i} - W_{j}(n)),\; j=1, \ldots , M \end{aligned}$$

(5)

where $h(x_{i}, j, n)$ is a neighborhood function taking values in [0, 1] and $\alpha (n)$ is a learning rate. $\alpha (n)$ is usually very small and decreases to zero with n. The common choice of the neighborhood function is the Gaussian function

$$\begin{aligned} h(x_{i}, j, n)= \exp {(- \frac{|i^{\star }(x_{i})-j|}{2 \sigma (n)^2})}, \end{aligned}$$

(6)

were $\sigma (n) \rightarrow 0$ when $n \rightarrow \infty$. There could be different methods of changing $\sigma$ and $\alpha$. Here we will assume that the learning is performed according to the following rules of learning parameters updating:

$$\begin{aligned} \alpha (n)= 0.5 \exp {(-n/1000)} \end{aligned}$$

(7)

and

$$\begin{aligned} \sigma (n)= \sigma _{0} 0.01^{n/1000},\; \sigma _{0}=20, \; 10, \; 5. \end{aligned}$$

(8)

However, it should be noted that in the case of ETSP also other formulas are used with success [4, 6, 25].

4.2 Algorithm

 

Step 1. :

Initialization: Let n be a discrete learning time starting from 0. The algorithm inputs are points $X=\{x_{1}, \ldots x_{N}\} \subset I^{d}$. Choose $M \ge 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_{1,max}]\times \ldots \times [w_{d,min},w_{d,max}]$, where $w_{k,min}=Min_{1 \le i \le N}[x_{i,k}]$, $w_{k,max}=Max_{1 \le i \le N}[x_{i,k}](k=1, \ldots ,d)$, and $(x_{i,1}, \ldots , x_{i,k}, \ldots , x_{i,d})^{T}= x_{i}, \;i=1, \ldots , 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 $\varPi _{opt}$. Start with $Q_{opt}=100$ and $\varPi _{opt}=(1,2, \ldots , N)$. Let the initial value of neighborhood function parameter be $\sigma _{0}=20$. And let $n_{max}$ be the value of the maximum number of iterations.

  The main loop of the algorithm consists of the following steps:  

Step 2. :

If $n > n_{max}$ go to Step End. Otherwise, compute parameters $\alpha (n)= 0.5 \exp {(-n/1000)}$ and $\sigma (n)= \sigma _{0} 0.01^{n/1000}$. Randomly draw point x from set of points X. Let say $x=x_{i}$.

Step 3. :

Find neuron $i^{\star }(x_{i})$ closest to $x_{i}$ (3). Update weights according to (5).

Step 4. :

Set $A_{j}= \{ \} \; j=1, \ldots , M$. For $i=1,2, \ldots , N$ find $j=i^{\star }(x_{i})$ and add index i to the list $A_{j}$.

Step 5. :

Obtain new sequence $\varPi$ by concatenating lists $A_{j}, \; j=1, \ldots , M$ (empty lists are rejected).

Step 6. :

Compute $Q(\varPi ,X)$ according to (1). If $Q(\varPi ,X) \le Q_{opt}$ set: $Q_{opt}=Q(\varPi ,X)$, $\varPi _{opt}=\varPi$ and $W_{opt}=W$. Set $n=n+1$ and go to Step 2.

Step End :

Return the best results: $\varPi _{opt}$ and $Q_{opt}$.

 

Fig. 1.

The part of the image set consisting of 20 section images of a mechanical device; the size of image $256 \times 256$ pixels

5 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 Figs. 1 and 6.

5.1 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 \times 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

$$\begin{aligned} {(20, 2, 3, 10, 9, 18, 13, 15, 5, 1, 17, 19, 4, 12, 7, 6, 11, 8, 16, 14)}. \end{aligned}$$

In Fig. 1 four subsequent slides with new numbers 13, 15, 5, 1, respectively, have been shown.

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 MPD (dashed), mean distance to the center of gravity MCD (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)

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

We have computed three image features: the area of the binarized object A, the maximal diameter distance of the object MPD and the mean distance to the center of gravity MCD (see Fig. 2 – left panel). The optimal sequence was obtained using the proposed here algorithm (see Sect. 4.2) as:

$$\begin{aligned} {(20, 2, 3, 10, 9, 18, 13, 5, 1, 15, 17, 19, 4, 12, 7, 6, 11, 8, 16, 14)}. \end{aligned}$$

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

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)

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 $\sigma _{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 Fig. 3. Figure 4 shows the sketch of the neuronal 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, Fig. 5 presents comparison of the 3-D image of a part of the mechanical device (on the left image) and its 3-D reconstruction obtained on the basis of the 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.

Fig. 6.

The image set consisting of 10 microscopy section images of the size $184 \times 184$ pixels

5.2 Sequencing Microscopy Section Images

The analyzed set of microscopy section images contains 10 slides of the size $184 \times 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 MCD. Due to the vary small differences between measurements and relatively large differences between mean values of the separate features, values of A and MCD features rescaled by 0.15 and 0.2, respectively. The rescaled 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$. Process 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 $\sigma _{0}=5$ and $\sigma _{0}=20$, respectively. Trajectories of the sum of the SSD obtained by the SOM algorithm in two exemplary runs are depicted on Fig. 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.

Fig. 7.

The results of the section images sequencing using the proposed method. SSD obtained in the subsequent iterations of the algorithm for 15 neurons net (left panel) and for 10 neurons net (right panel)

Fig. 8.

The 15 neurons path obtained for optimal image sequence for starting value of $\sigma =20$ (left panel) and $\sigma =5$ (right panel). Large dots indicate position of images

Fig. 9.

Two views of 3-D reconstruction of an biological object obtained on the basis of the optimal with respect to SSD section sequence after additional alignment

6 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 (Figs. 4 and 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 [7] (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.