Edge Detection in Biomedical Images Using Self-Organizing Maps

The application of self-organizing maps (SOMs) to the edge detection in biomedical images is discussed. The SOM algorithm has been implemented in MATLAB program suite with various optional parameters enabling the adjustment of the model according to the user’s requirements. For easier application of SOM the graphical user interface has been developed. The edge detection procedure is a critical step in the analysis of biomedical images, enabling for instance the detection of the abnormal structure or the recognition of different types of tissue. The self-organizing map provides a quick and easy approach for edge detection tasks with satisfying quality of outputs, which has been verified using the high-resolution computed tomography images capturing the expressions of the Granulomatosis with polyangiitis. The obtained results have been discussed with an expert as well.


Introduction
The application of self-organizing maps (SOMs) to the edge detection in biomedical images is discussed. The SOM algorithm has been implemented in MATLAB program suite with various optional parameters enabling the adjustment of the model according to the user's requirements. For easier application of SOM the graphical user interface has been developed. The edge detection procedure is a critical step in the analysis of biomedical images, enabling for instance the detection of the abnormal structure or the recognition of different types of tissue. The self-organizing map provides a quick and easy approach for edge detection tasks with satisfying quality of outputs, which has been verified using the high-resolution computed tomography images capturing the expressions of the Granulomatosis with polyangiitis. The obtained results have been discussed with an expert as well.

Self-organizing map 2.1. Self-organizing map in edge detection performance
The self-organizing map (SOM) [5,9] is widely applied approach for clustering and pattern recognition that can be used in many stages of the image processing, e. g. in color image segmentation [18], generation of a global ordering of spectral vectors [26], image compression [25], binarisation document [4] etc.
The edge detection approaches based on SOMs are not extensively used. Nevertheless there are some examples of SOM utilization in edge detection, e. g. texture edge detection [27], edge detection by contours [13] or edge detection performed in combination with conventional edge detector [24] and methods of image de-noising [8] .
In our case, the SOM has been utilized in edge detection process in order to reduce the image intensity levels.

Structure of self-organizing map
A SOM has two layers of neurons, see Figure 1. The input layer (size Nx1) represents input data x 1 , x 2 , . . . , x M (M inputs, each input is N dimensional). The output layer (size KxL), that may have a linear or 2D arrangement, represents clusters in which the input data will be grouped. Each neuron of the input layer is connected with all neurons of the output layer through the weights W (size of the weight matrix is KxLxN).

Training of self-organizing map
A SOM is neural network with unsupervised type of learning, i. e. no cluster values denoting an a priori grouping of the data instances are provided.
The learning process is divided in epochs, during which the entire batch of input vectors is processed. The epoch involves the following steps: 1. Consecutive submission of an input data vector to the network.
2. Calculation of a distances between the input vector and the weight vectors of the neurons of the output layer.
3. Selection of the nearest (the most similar) neuron of the output layer to the presented input data vector.
4. An adjustment of the weights.
SOM can be trained in either recursive or batch mode. In recursive mode, the weights of the winning neurons are updated after each submission of an input vector, whereas in batch mode, the weight adjustment for each neuron is made after the entire batch of inputs has been processed, i. e. at the end of an epoch.
The weights adapt during the learning process based on a competition, i. e. the nearest (the most similar) neuron of the output layer to the submitted input vector becomes a winner and its weight vector and the weight vectors of its neighbouring neurons are adjusted according to where W is the weight matrix, x i the submitted input vector, λ the learning parameter determining the strength of the learning and φ s the neighbourhood strength parameter determining how the weight adjustment decays with distance from the winner neuron (it depends on s, the value of the neighbourhood size parameter).
The learning process can be divided into two phases: ordering and convergence. In the ordering phase, the topological ordering of the weight vectors is established using reduction of learning rate and neighbourhood size with iterations. In the convergence phase, the SOM is fine tuned with the shrunk neighbourhood and constant learning rate. [23] 2.3.1. Learning Parameter The learning parameter, corresponding to the strength of the learning, is usually reduced during the learning process. It decays from the initial value to the final value, which can be reached already during the learning process, not only at the end of the learning. There are several common forms of the decay function (see Figure 2):

No decay
2. Linear decay 3. Gaussian decay 4. Exponential decay where T is the total number of iterations, λ 0 and λ t are the initial learning rate and that at iteration t, respectively. The learning parameter should be in the interval 0.01, 1 .

Neighbourhood
In the learning process not only the winner but also the neighbouring neurons of the winner neuron learn, i. e. adjust their weights. All neighbour weight vectors are shifted towards the submitted input vector, however, the winning neuron update is the most pronounced and the farther away the neighbouring neuron is, the less its weight is updated. This procedure of the weight adjustment produces topology preservation.
There are several ways how to define a neighbourhood (some of them are depicted in Figure 3).
The initial value of the neighbourhood size can be up to the size of the output layer, the final value of the neighbourhood size must not be less than 1. The neighbourhood strength parameter, determining how the weight adjustment of the neighbouring neurons decays with distance from the winner, is usually reduced during the learning process (as well as the learning parameter, analogue of Equations 2-5 and Figure 2). It decays from the initial value to the final value, which can be reached already during the learning process, not only at the end of the learning process. The neighbourhood strength parameter should be in the interval 0.01, 1 . The Figure 4 depicts one of the possible development of neighbourhood size and strength parameters during the learning process.

Weights
The resulting weight vectors of the neurons of the output layer, obtained at the end of the learning process, represent the centers of the clusters. The resulting patterns of the weight vectors may depend on the type of the weights initialization. There are several ways how to initialize the weight vector, some of them are depicted in Figure 5.

Distance Measures
The criterion for victory in the competition of the neurons of the output layer, i. e. the measure of the distance between the presented input vector and its weight vectors, may have many forms. The most commonly used are:

Block distance
where x i is i-th component of the input vector, w ji i-th component of the j-th weight vector, N dimension of the input and weight vectors, x mean value of the input vector x, w j mean value of the weight vector w j , σ x standard deviation of the input vector x, σ w j standard deviation of the weight vector w j , x i length of the input vector x and w ji length of the weight vector w j .

Learning Progress Criterion
The learning progress criterion, minimized over the learning process, is the sum of distances between all input vectors and their respective winning neuron weights, calculated after the end of each epoch, according to where x n is the n-th input vector belonging to cluster c i whose center is represented by w i (e. i. the weight vector of the winning neuron representing cluster c i ).
The weight adjustment corresponding to the smallest learning progress criterion is the result of the SOM learning process, see Figure 6. These weights represent the cluster centers.
For the best result, the SOM should be run several times with various settings of SOM parameters to avoid detection of local minima and to find the global optimum on the error surface plot.

Errors
The errors of trained SOM can be evaluated according to 1. Learning progress criterion (see Equation 10), 2. Normalized learning progress criterion 3. Normalized error in the cluster 4. Error in the i-th cluster 5. Normalized error in the i-th cluster where x n is the n-th input vector belonging to cluster c i whose center is represented by w i (e. i. the weight vector of the winning neuron representing cluster c i ), M is number of input vectors, M i is number of input vectors belonging to i-th cluster and k is number of clusters.
For more information about the trained SOM, the distribution of the input vectors in the clusters and the errors of the clusters can be visualized, see Figure 7.

Forming Clusters
The U-matrix (the matrix of average distance between weights vectors of neighbouring neurons) can be used for finding of realistic and distinct clusters. The other approach for forming clusters on the map can be to utilize any established clustering method (e. g. K-means clustering). [23,28] 2.3.8. Validation of Trained Self-organizing Map The validation of the trained SOM can be done so, a portion of the input data is used for map training and another portion for validation (e. g. in proportion 70:30). Different approach for SOM validation is n-fold cross validation with the leave-one out method. [16,23]

Using of trained self-organizing map
A trained SOM can be used according the following steps for clustering: 1. Consecutive submission of an input data vector to the network.
2. Calculation of a distances between the input vector and the weight vectors of the neurons of the output layer.
3. Selection of the nearest (the most similar) neuron of the output layer (e. i. the cluster) to the presented input data vector.

Implementation of self-organizing map
The SOM algorithm has been implemented in MATLAB program suite [17] with various optional parameters enabling the adjustment of the model according to the user's requirements. For easier application of SOM the graphical user interface has been developed facilitating above all the setting of the neural network, see Figure 8.

Edge detection
Edge detection techniques are commonly used in image processing, above all for feature detection, feature extraction and segmentation.
The aim of the edge detection process is to detect the object boundaries based on the abrupt changes in the image tones, i. e. to detect discontinuities in either the image intensity or the derivatives of the image intensity.

Conventional edge detector
The image edge is a property attached to a single pixel. However it is calculated from the image intensity function of the adjacent pixels.
Many commonly used edge detection methods (Roberts [22], Prewitt [20], Sobel [19], Canny The essential step in edge detection process is thresholding, i. e. determination of the threshold limit corresponding to a dividing value for the evaluation of the edge detector response either as the edge or non-edge. Due to the thresholding, the result image of the edge detection process is comprised only of the edge (white) and non-edge (black) pixels. The quality of thresholding setting has an impact on the quality of the whole edge detection process, i. e. exceedingly small value of the threshold leads to assignment of the noise as the edge, on the other hand exceedingly large value of the threshold leads to omission of some significant edges.

Self-organizing map
A SOM may facilitate the edge detection task, for instance by reducing the dimensionality of an input data or by segmentation of an input data.
In our case, the SOM has been utilized for the reduction of image intensity levels, from 256 to 2 levels. Each image has been transformed to mask (3x3), see Table 1, that forms the set of input vectors (9-dimensional). The input vectors have been then classified into 2 classes according to the weights of the beforehand trained SOM. The output set of the classified vectors has been reversely transformed into the binary (black and white) image.
Due to this image preprocessing using SOM, the following edge detection process has been strongly simplified, i. e. only the numerical gradient has been calculated where G is the edge gradient, G x and G y are values of the first derivative in the horizontal and in the vertical direction, respectively.
The ability of SOM for edge detection in biomedical images has been tested using the high-resolution computed tomography (CT) images capturing the expressions of the Granulomatosis with polyangiitis disease.

Diagnosis of granulomatosis with polyangiitis
Granulomatosis with polyangiitis (GPA), in the past also known as Wegener's granulomatosis is a disease belonging to the group of vasculitic diseases affecting mainly small caliber blood vessels [7].
They can be distinguished from other vasculitides by the presence of ANCA antibodies (ANCA -Anti-Neutrophil Cytoplasmic Antibodies). GPA is quite a rare disease, its yearly incidence is around 10 cases per million inhabitants. The course of the disease is extremely variable on the one hand there are cases of organ limited disease that affects only single organ, on the other hand there is a possibility of generalized disease affecting multiple organs and threatening the patients life. Diagnostics of this disease has been improved by the discovery of ANCA antibodies and their routine investigation since the nineties of the past century. The onset of GPA may occur at any age, although patients typically present at age 35-55 years [11].
A classic form of GPA is characterized by necrotizing granulomatous vasculitis of the upper and lower respiratory tract, glomerulonefritis, and small-vessel vasculitis of variable degree. Because of the respiratory tract involvement nearly all the patients have some respiratory symptoms including cough, dyspnea or hemoptysis. Due to this, nearly all of them have a chest X-ray at the admission to the hospital, usually followed by a high-resolution computed CT scan. The major value of CT scanning is in the further characterization of lesions found on chest radiography.
The spectrum of high-resolution CT findings of GPA is broad, ranging from nodules and masses to ground glass opacity and lung consolidation. All of the findings may mimic other conditions such as pneumonia, neoplasm, and noninfectious inflammatory diseases [2].
The prognosis of GPA depends on the activity of the disease and disease-caused damage and response to therapy. There are several drugs used to induce remission, including cyclophosphamide, glucocorticoids or monoclonal anti CD 20 antibody rituximab. Once the remission is induced, mainly azathioprin or mycophenolate mofetil are used for its maintenance.
Unfortunately, relapses are common in GPA. Typically, up to half of patients experience relapse within 5 years [21].

Analysis of granulomatosis with polyangiitis
Up to now, all analysis of CT images with GPA expressions have been done using manual measurements and subsequent statistical evaluation [1,3,10,12,14]. It has been based on the experience of a radiologist who can find abnormalities in the CT scan and who is able to classify types of the finding. The standard software 1 used for CT image visualization usually includes some interactive tools (zoom, rotation, distance or angle measurements, etc.) but no analysis is done (e. g. detection of the abnormal structure or recognition of different types of tissue). Moreover, there is no software customized specifically for requirements of GPA analysis.
Therefore, there is a place to introduce a new approach for analysis based on SOM. CT finding analysis is then less time consuming and more precise.

Results and discussion
The aim of the work has been to detect all three expression forms of the GPA disease in high-resolution CT images (provided by Department of Nephrology, First Faculty of Medicine and General Faculty Hospital, Prague, Czech Republic) using the SOM, i. e. to detect granulomatosin, mass and ground-glass, see Figure 10a, 11a, 12a. The particular expression forms occur often together, therefore there has been the requirement to distinguish the particular expression forms from each other.
Firstly, the SOM has been trained using the test image, see Figure 9 (For detailed information about the training process, please see Table 2).
Secondly, the edge detection of the validating images has been done using the result weights from the SOM training process, see Figure 10b, 11b, 12b (For detailed information about the edge detection process, please see Table 3).

Setting description Value
Size of the output layer    In the first case, a high-resolution CT image capturing a granuloma in the left lung has been processed, see Figure 10a. The expert has been satisfied with the quality of the GPA detection (see Figure 10b) provided by the SOM, since the granuloma has been detected without any artifacts.
In the second case, the expert has pointed out that the high-resolution CT image, capturing both masses and ground-glass (see Figure 11a, 12a), was problematic to be evaluated. The possibility of the detection of the masses was aggravated by the ground-glass surrounding the lower part of the first mass and the upper part of the second mass. The artifacts, which originated in the coughing movements of the patient ('wavy' lower part of the image), made the detection process of the masses and the ground-glass difficult as well. Despite these inconveniences, the expert has confirmed, that the presented software has been able to distinguish between the particular expression forms of GPA with satisfying accuracy and it has detected the mass and ground-glass correctly (see Figure 11b, 12b, 12c).
In conclusion, the presented SOM approach represents a new helpful approach for GPA disease diagnostics.

Conclusions
The edge detection procedure is a'critical step in the analysis of biomedical images, enabling above all the detection of the abnormal structure or the'recognition of different types of tissue.
The application of SOM for edge detection in biomedical images has been discussed and its contribution to the solution of the edge detection task has been confirmed. The'ability of SOM has been verified using the'high-resolution CT images capturing all three forms of the expressions of the GPA disease (granulomatosin, mass, ground-glass). Using SOM, particular expression forms of the GPA disease have been detected and distinguished from each other. The obtained results have been discussed by the expert who has confirmed that the SOM provides a quick and easy approach for edge detection tasks with satisfying quality of output.
Future plans are based on the problem extension to three-dimensional space to enable CT image analysis involving (i) pathological finding 3D visualization and (ii) 3D reconstruction of the whole region (using the whole set of CT images).