Determining similarity in histological images using graph-theoretic description and matching methods for content-based image retrieval in medical diagnostics

Overview of Content Based Image Retrieval

In Information Retrieval (IR) systems, the user specifies a query either in the form of text, documents, images, or sounds and the system is expected to return the items that are semantically similar to the query in some sense. CBIR is an information retrieval system that includes techniques for retrieving digital images by their visual content. The horizon of CBIR includes methods ranging from image similarity functions to highly complex image annotation systems [9].

At present, CBIR is an extremely active area of research. Descriptions of a variety of CBIR approaches implemented in the past are given in reviews [10] and [11]. CBIR has been applied to medical domain and a comprehensive review on medical CBIR systems is given by Müller et al. [12]. However, most of the recently developed retrieval methods are dedicated to radiological images [13]. Specifically for histological images, research in this field has been comparatively less. An application with histopathological images is described in [14], using a property concept frame representation for morphological characteristics based on fuzzy logic. However, it does not emphasize on the spatial relationships between the various tissue components, which are considered an important aspect in our work, in order to describe the overall topology of the breast tissues.

In Diamond Project [15], the interactive search in large distributed data repositories was addressed. Particularly, the most relevant to medical domain are MassFind [16], FatFind [17] and PathFind [18]. MassFind is an application for diagnosing lesions in mammograms, which focuses on performance of different distance metrics to define similarity between ROI images. FatFind exploits the property of perfect round shape of adipocytes in cell microscopy images for their automatic counting, by making use of low-level shape features of cells. PathFind is a tool employing “discard-based search” for content-based retrieval of WSIs. However, in all the applications, less attention is given to high-level structural representation and retrieval algorithms specific to histological images; but more emphasis is on development of design and implementation strategies of the search methods, to handle huge data collections in large-scale efficient network-distributed frameworks.

Graph Theory

A graph is a set containing a finite number of points, called nodes (or vertices), which are connected by lines called edges (or arcs). In this paper, a graph is considered as a 4-tuple G=(V,E,α,β), where

Figure 1 gives an example of a basic graph.

Attributed Relational Graph

A graph G is said to be an Attributed Relational Graph (ARG) when both the nodes and the edges are represented with attributes. The node attributes for node n _i are denoted as a vector $a_i=\left[a_i^{\left(k\right)}\right],(k=1,2,3,\dots ,K)$, where K is the number of node attributes in the vector a _i , and the edge attributes (or weights) for edge e _j by the vector denoted as $b_j=\left[b_j^{\left(m\right)}\right],(m=1,2,3,\dots ,M)$, where M is the number of edge attributes in the vector b _j . In Figure 2, ARG with node attribute vectors a _i , i=1,2,3 and edge attribute vectors b_j,j=1,2,3 is shown. Node attributes represent quantities such as size, position, shape and colour of an object whereas edge attributes define relationships between nodes like the distance between two points or dissimilarity between objects. ARGs act as convenient structures for physical representation and are frequently used in applications ranging from computer-aided design to machine vision [21]. ARGs have been employed in this work for representing the information content of histological images.

Regional Adjacency Graph

Regional Adjacency Graph (RAG) is an ARG whose vertices represent regions and edges represent connections between adjacent regions. Node attributes are assigned according to characteristics of the region corresponding to each node and edge attributes (or weights) describe the adjacency relationships. RAGs give a spatial view of the images and are effective in applications for representing image information where neighbourhood relationships can be taken into account. RAGs have been used in [22] for segmentation of colour images. In this work, segmented histological images are represented as RAGs. Figure 3 shows a simple RAG example.

Graph Matching

Graph Matching is the process of comparing two graphs to find an appropriate correspondence between their nodes and edges. It refers to the process of finding a mapping F from the nodes of one graph G to the nodes of another graph G ^′ that satisfies some constraints or optimality criteria, ensuring that similar substructures in one graph are mapped to similar substructures in the other. Standard structural matching concepts include the following:

1. 1.

   Graph Isomorphism: It finds an exact structural correspondence between two graphs. It is a bijective mapping that preserves the number of nodes and edges. It is illustrated in Figure 4.

    

1. 2.

   Subgraph Isomorphism: If nodes along with their corresponding edges are deleted from a graph G, a subgraph G ^′ denoted by G ^′ ⊆ G is obtained. A subgraph isomorphism from G to G” is an isomorphism from a graph G” to a subgraph G ^′ of G. It is shown in Figure 5.

    

1. 3.
   Monomorphism: It is a more relaxed matching than subgraph isomorphism as extra edges are also allowed between nodes in the larger graph. Figure 6 illustrates monomorphism from graph G to graph G ^′ . Formally, it can be stated as: Let G and G’ be graphs. A graph monomorphism between G(V,E,α,β) and G ^′ (V ^′ ,E, ^′ α ^′ ,β ^′ ) is an injective mapping F _mono :V→V ^′ such that:

   $\alpha \left(v\right)={\alpha }^{\prime }\left(F_{\mathrm{mono}}\right(v\left)\right)\forall v\in \mathrm{V.}$

   (1)
    

For any edge e=(u,v) ∈ E, there is an edge e ^′ =(F _mono (u),F _mono (v)) ∈ E ^′ such that β(e)=β ^′ (e ^′ ).

1. 4.

   Maximum Common Subgraph (MCS): An MCS of two graphs, G and G ^′ , is a graph G” that is a subgraph of both G and G ^′ , such that it has the maximum number of nodes among all possible subgraphs of G and G ^′ . MCS of two graphs is usually not unique. It can be used to measure the similarity of objects as the larger the MCS, higher will be the similarity. It is shown in Figure 7. A common subgraph of G and G ^′ , CS(G,G ^′ ), is a graph G” such that there exist subgraph isomorphisms from G to G” and from G ^′ to G” or vice-versa. G” is a MCS of G and G ^′ , MCS(G,G ^′ ), if it is the common subgraph with maximum nodes.

    

The basic A* Search Algorithm

Dijkstra’s algorithm [25] starts with the source node and traverses the nodes in a graph such that shortest path from the source, found so far, is prolongated first. Thus, by reaching the goal node the shortest path is guaranteed to be found. On the other hand, the Greedy Best-First-Search algorithm [26] selects the node closest to the goal by using a heuristical estimate of the distance of a node from the goal node irrespective to distance to source and thus finds a path to the goal in shortest time, which is not necessarily the shortest path. A* algorithm [27] was developed to combine formal approaches like Dijkstra’s algorithm and heuristic approaches like Greedy Best-First-Search algorithm.

Implementation

A* algorithm can be implemented by maintaining of two lists: the Open List and the Closed List. The Open List contains nodes that are candidates for examining. It is generally maintained as a priority queue, as the node with highest priority is the one with least f(n) cost. Initially, it contains just one element: the source node. The Closed List contains those nodes that have already been traversed and form the optimal path. At each step of the algorithm, the node n with the lowest f(n) value is examined from the Open List. If n is the goal, then algorithm stops. Otherwise, it is removed from Open List and added to Closed List, and the f(n), g(n) and h(n) values of its child nodes are updated accordingly. These nodes are then inserted into the Open List queue and are synchronised according to priorities of other existing elements. The process continues till the goal node is reached, or no more nodes are available in Open List (a case of no solution). The algorithm is shown in stepwise manner in Figure 8.

Image segmentation

The maximum size of the pixel area for decision making is tissue type related. In these experiments 16× 16 pixels for epithelial, 32×32 pixels for connective, 48×48 for lobular and 64×64 for fat tissues were used. Here, the effective pixel size for these images (i.e. how large the physical area of the tissue which corresponds to one pixel) was roughly 1.0 micrometer ×1.0 micrometer. The segmentation algorithm is not described in further details here and is a subject to a separate publication. Segmentation results were just provided for this study. The segmentation is done into four tissue types: lobules, fibrous connective tissue, epithelial lining cells, lumens-&-fat (centres of ducts and adipose tissue).

The multilabel (L=4 here) segmented image is decomposed into binary images, one image for each label. Then morphological operations closing and opening are performed twice each on each binary image. These operations aim to remove small artefacts, fill in the potential gaps between tissue fragments and smooth the contours of the shapes. The size of structuring element chosen depends on the magnification of the WSIs used in the study. Then connected components are identified in each binary image. A connected component analysis ensures that only connected pixels are assigned the same label and form a region. It is required for distinguishing the regions within the image.

Graph-theoretic Representation

Each of the images in the database as well as the query image have been described by corresponding graphs. Namely, ARGs have been constructed, where each node corresponds to one connected region in the image and edge is obtained between neighbouring regions which share a common boundary. The procedure involves describing nodes and edges with attributes explained below.

Edge Description

The process of describing edges involves identifying the edges and assigning weights to them. The edge information (weights) is obtained as:

• Distance between centroids: It is taken as the Euclidean distance between the centroids of two regions.

• Common boundary length: It is the number of pixels lying on the common border between two neighbouring regions. It has been calculated by considering the 4-connectivity of each pixel. The algorithm counts those 4-connected neighbours of a pixel which have a different label than the pixel itself.

Same as with nodes, other characteristics can also be included in edge attributes. However, it was found that the properties given above are suitable for representing a histological image. Area is important for matching of similarly-sized regions, whereas perimeter conveys approximate shape information. The distance between centers determines how far the nodes are placed with respect to each other, and common boundary length denotes upto what extent two regions are adjacent to each other. Thus, these properties are most useful for determining the structure and neighbourhood relationships for histological image analysis. One example of such a graph-representation is given in Figure 10.

A*-based graph matching method

Given a query image, its ARG is matched to each ARG in the database and the cost of matching is assigned to each pair of graphs. The graph matching problem has been formulated as an A* based tree search problem. Functions for the cost g(n), heuristic h(n) and total cost f(n) have been designed using the information present in the corresponding image ARGs. The heuristic h(n) is designed to be a consistent lower bound estimate of the exact cost, hence, admissibility criterion is satisfied that leads to the optimal solution.

For the proposed method, K and M both have value 2, where k=1 for area and k=2 for perimeter in node attributes, and m=1 for distance between centroids and m=2 for common boundary length in edge attributes. Note that for node attributes, area has been assigned double the weight of perimeter, as it is considered more important feature during the matching of nodes.

where, δ _i in equation 10 describes the distance between the attributes of already matched nodes and edges, and $a_i^{\left(1\right)}$ in equation 11 refers to areas of the nodes in G not matched yet. Note that only the area attribute is used at this point as computing δ _i will not be possible for unmatched nodes and the most important attribute that needs to be considered in the heuristic function is area.

where n _p ∈ G and $n_q^{\prime }\in G^{\prime }$ are matching nodes and a _p ⁽¹⁾ and ${a_q^{\prime }}^{\left(1\right)}$ denote the first attributes corresponding to feature vectors a _p and ${\mathbf{a}}_{\mathbf{q}}^{\mathbf{\prime }}$. They describe the area of the two regions being matched.

The graph matching is implemented through a tree-based search using A* algorithm which always extends the partial mapping of nodes towards an optimum. The tree is the representation of Open List realized using a priority queue, containing partial mappings in increasing order of their costs. It is constructed by first allowing each test node to be mapped to each available model node, if the match is permissible, the pairs forming the first level of the tree. The cost of each pair is computed, and the pair with the lowest total cost f(n) is expanded. Each leaf of the tree now represents a combination of matched nodes or partial mapping from nodes of test graph to those of model graph. The Closed List consists of the latest and most favourable partial mapping constructed. The tree is expanded until best optimum mapping of maximum nodes is found. An example of the tree-based search method employed is illustrated in Figure 11.

The main problem with optimal graph matching is its high computational complexity. The complexity of the described search is exponential in the worst case, however, practically, it depends on the data to be handled as only nodes of same label can be matched. It considerably reduces the search space and complexity is scaled down.

Dataset

The data used for this work consists of histological images provided by The Charité Hospital, Berlin. These are biopsy images of the breast tissue. The samples have been stained with the H&E dye. The WSI images are produced by a Zeiss MIRAX SCAN WSI scanner. We used selected archived slides from daily clinical workload that were not older than 6 months at the time of digitalization. The glass slides have been produced in AP-laboratory of the Institute of Pathology at Charité hospital. They have not been modified in any way. We have evaluated the method on 3 WSI images of FEA-suspected breast biopsies divided into sub-images representing possible retrieval results. Our aim was to demonstrate the potential of the graph-based approach, leaving the in-depth performance evaluation for future research. One of the reasons for this is the relatively high computational complexities of the segmentation, the description and the retrieval algorithm, which are subject for future change and improvements too.

The images have been pre-segmented to four categories describing different types of tissue. They are then divided for one approach (described in Section Experimental Approaches and Results) into smaller sub-images to obtain the database for different image sizes of 64×64, 128×128, 256×256 and 512×512. Query image is selected by giving a choice of four sizes, and the selection is resized to the size selected by the user. The number of images used in the database, depending on the size of query image is:

64×64: 70869 images

128×128: 27596 images

256×256: 9132 images

512×512: 2485 images

Graph representations of the images stated above were obtained and stored for future reference.

Experimental Approaches and Results

Two types of configurations were used for experiments. These are as follows:

Subgraph isomorphism approach

The approach aims to find the subgraph isomorphism between a smaller query graph and the graphs obtained for the whole size images. The user submits a query image of any size by selecting a rectangular section of the whole histological image presented to him. The group of pre-segmented regions present in the selection is identified. Each of the regions is then extended to its original size and shape. An RAG is then constructed whose attributes are computed from the properties and spatial information of the extended regions. The matching process returns those subgraphs of the graphs obtained from the database of entire WSIs, which are closest to the graph obtained for the query. An example graph, generated for an entire histological image, is shown in Figure 12. The selection of the query image, the obtained RAG and three nearest matching results are shown in Figure 13. Here, a query consisting of lobular cells surrounded by epithelium and adjacent to a layer of fat tissue is selected. The closest matches show similar structural groups of lobules, fat cells and epithelial cells.

Inexact graph matching approach

In this approach, an inexact matching between a pair of graphs generated for a query image and a database image of the same size is determined. For this, first the database images are divided into smaller sub-images, query image is selected of a predefined size, and then the graphs of each sub-image is compared with the graph of the query image. The sub-images with closest similarity are retrieved. Note that in contrast to the previous method the regions of the query image are not extended to their original size. An example of this approach with results shown for a query size of 512×512 is shown in Figure 14. The query image has a duct with a lumen (center) with an outer lining of epithelial tissues and having lobules and connective tissue in the background. The retrieved results show similar regions. First result is not exact but depicts similar physiological structures and spatial relationships between them. Moreover, although the first and the third match basically show the same histological structure, due to the splitting of the whole histological image into sub-images they appear as different images in our database. They both have been selected due to the high similarity to the query image.

Observations

It can be observed in both approaches that:

Performance evaluation

where P _s is precision (in %), s is scope length and score is a value from {0,0.25,0.5,0.75,1}. The score values express the similarity of the retrieved images to the query image in terms of structure and spatial relationships. Higher score is assigned for higher resemblance. The evaluation is subjective and coarse, so quantitative results (precision values) have been rounded down to integers. Resulted plots are plotted between the precision vs. different scope lengths.

The proposed method has been compared with a common, histogram-based retrieval system. The histograms for segmented sub-images have been found using 4 bins. Then the distances between the histograms of query image and sub-images have been calculated. Similar as for the graph-based approach, exponential distance has been used. Finally, the results were compared for both methods.

For the histogram-based technique, the table of precision for different scope lengths at different window sizes (or query sizes) is given in Table 1, with corresponding line plots illustrated in Figure 15. The table of precision for scope lengths at the same query sizes and the line plots for graph-theoretic method are shown in Table 2 and Figure 16, respectively. Next, a comparison is established between the average precision of both methods, calculated across all window sizes. The relative improvement of graph-based technique over histogram-based technique is given in Table 3, with corresponding line plots visualised in Figure 17.

The tables and graphs obtained justify our choices of the parameters used and methods employed for the system proposed. It can be concluded that:

The results obtained using graph-theoretic technique are better than simple histogram based method, as it takes into account the structural characteristics of the image and neighbourhood relationships between regions, which are completely neglected in the histogram-based method.

As scope length increases, the precision declines for proposed method, which shows that it gives the most relevant results earlier in the list of retrieved results. This is a desirable property of any CBIR system that the results initially obtained are the most useful. However, it is evident that it does not hold for all cases of histogram-based method.

The results so obtained by the proposed method are not as high as reported for general CBIR applications (about 90% precision or more). The highest precision reported for image size 64×64 and scope length 10 is 80%. The reason behind this is the complexity and subjectivity associated with histological images. The evaluation was biased strongly with the subjective scoring as even when the histological image shows the same tissue composition, several factors have to be kept in mind before assigning a score. The closeness to query image depends on the type of tissue regions, size and shape of regions as well as the neighbourhood relationships between them. Due to this relative scores have been used, however, more thorough evaluation should be performed especially by employing medical professionals.

The performance depends on the characteristics of the query image, i.e. the number of tiled images available in the database that lie close to the position of the query window.

Execution Time Requirements