junctions and each local line segment of the CCDZ can participate in several groups, the exact arms of a junction could therefore be greater than the line segments forming this junction (Figure 2.9). Given a local topology represented by $n$ straight line segments $\{P_iQ_i\}$ with $1 \le i \le n$, the process of determination of the exact arms of each junction is as follows:

* Let $O_J$ be a set of arms of junction $J$ where $O_J \leftarrow \emptyset$ at the beginning for every junction.

* If the line segment $P_iQ_i$ is clustered in a group whose an optimal junction $J$ is then constructed, the line $P_iQ_i$ is considered as one of the arms of the junction $J$: $O_J \leftarrow O_J \cup \{P_iQ_i\}$.

* For each line segment $P_iQ_i$ that is clustered in several groups, the corresponding junctions involved in $P_iQ_i$ are sorted in the order of increasing distance to $P_i$. Then, for each junction $J$ except the last one in the list, the corresponding set $O_J$ is updated as: $O_J \leftarrow O_J \cup \{JG\}$, where $JG$ is a straight line segment constructed at $J$ with the same length as $P_iQ_i$ but the point $G$ lies in the opposite direction of vector $P_i\tilde{Q}_i$.

<figure><img src="image_48.png" /><figcaption>Figure 2.9: (a) A CCDZ cropped from Figure 2.7 with the superposition of three clusters; (b) detected junctions $\{J_1, J_2, J_3\}$; (c) the detected junctions are classified as a 3-junction even though $J_2$ and $J_3$ are constructed from two clusters, $\{1, 2\}$ and $\{1, 5\}$, respectively, each of which only contains two local line segments.</figcaption></figure>

Once we have correctly determined the arms of each junction, the junction characterization could easily be accomplished as follows. Given a junction $J$ associated with a set of $m$ arms $\{U_iV_i\}_{i=0,...,m-1}$, the characterization of this junction is described as $\{p, s_p, \{\theta_i^p\}_{i=0}^{m-1}\}$, where:

* $p$ is the location of $J$;

* $s_p$ is the local scale computed as the mean length of the arms of $J$:

$$s_p = \frac{1}{m} \sum_{i=0}^{m-1} |U_iV_i|$$

* $\theta_i^p$ is the difference in degrees between two consecutive arms $U_iV_i$ and $U_{i+1}V_{i+1}$. These parameters $\{\theta_i^p\}_{i=0}^{m-1}$ are tracked in the counterclockwise direction and the $\theta_{m-1}^p$ is the difference in degrees between the arms $U_{m-1}V_{m-1}$ and $U_0V_0$.
It is noted that a similar way of junction characterization has been also exploited in the CV field, taking the work of [Xia, 2011] for example. In their work, the junction arms $\{\theta_i\}_{i=0}^{m-1}$ are described using the absolute angles in a 2D plane. Here, we use the relative difference in degrees between two successive arms to make the junction characterization invariant to the image plane and make it easier for further processes of junction matching.

The description of each junction point derived in this way is rather compact, distinctive, and general. The dimension of this descriptor is variable but limited to the number of arms of each junction point, and in practice, this value is quite small (e.g., 3 for a T-junction, 4 for an X-junction). This point constitutes a great advantage of the detected junctions that provides a very efficient approach to the subsequent task of junction matching. In addition, the junction descriptor is distinctive and general, such that we can describe any junction points appearing in a variety of complex and heterogeneous documents. After this step, junction matching can be performed by simply comparing the descriptors of two junctions. Figure 3.5 shows the corresponding matches of the junctions detected in a query symbol (left) and those of an image cropped from a large document (right). For simplicity, the matches are shown after performing geometry checking using the Generalized Hough Transform [Ballard, 1981].

<figure><img src="image_49.png" /><figcaption>Figure 2.10: Corresponding junction matches between a query symbol (left) and a cropped document (right).</figcaption></figure>

## 2.3 Complexity evaluation

In this section, we provide a detailed analysis of the complexity of the proposed method given an image *I* of the size *M* × *N*. In the pre-processing stage, before applying the (3,4)-distance transform skeletonization algorithm, several basic pre-processing steps, such as hole filling, small contour removing, and image dilation, are performed, as discussed in the original work of Di Baja [di Baja, 1994]. Such steps can be processed in parallel using one scan over the image. The skeletonization step is then applied, which requires two scans of the image to calculate the (3,4)-chamfer distance. In summary, the computation complexity for the pre-processing stage is basically linear (i.e., $O(MN)$).
In the next stage of scale selection and dominant point detection, the ROS determination step is applied for each skeleton point, thus using a single loop of length S, where S is the number of skeleton points. The technique of least squares line fitting is a second-order linear computation of the length that it traverses. In practice, it is not necessary to traverse a full skeleton branch; rather, a short path of the branch with a length $k_ρ = 50$ (pixels) may be sufficient, for example. As the technique of least squares line fitting is performed in both directions at each point, it is equal to a complexity of $O(2Sk_ρ^2)$ in total for this step. The 2-junctions are then detected as dominant points by applying Teh-Chin's algorithm, which is a sequential 4-pass process, where the first pass is performed on the full length of the median lines to detect a list of H candidate dominant points and the other passes are conducted on one of these candidate points, where H is much smaller than S. Furthermore, the crossing-points could be detected in parallel with a cost of first-order linear polynomial time $O(S)$. The overall computation complexity for these processes is essentially linear to the length of the median points (i.e., $O(Sk_ρ^2)$).

For the last stage of junction reconstruction, let $K$ be the number of candidate junctions comprised of 2-junction points and crossing-points. The distorted zone $Z_i$ defined at each candidate junction $p_i$ ($1 \le i \le K$) has an area of $\pi r_i^2/4$ where $r_i$ is the line thickness at $p_i$. Given a distorted zone $Z_i$, the maximum complexity of the computation to find an optimal junction in $Z_i$ is $O(T\pi r_i^2/4)$ where the first factor, $T$, is the number of times that steps $\{1, 2, 3\}$ are repeated and the second factor, $\pi r_i^2/4$, is the number of foreground pixels in $Z_i$ (i.e., the local searching neighborhood). Our investigation has shown that the number of iterations, $T$, is very small and is typically less than 10. As $Z_i$ can contain multiple junctions, say $L$ junctions, it implies that the junction optimization process applied to $Z_i$ will be terminated after running $L$ iterations of the steps $\{1, 2, 3, 4, 5\}$. The value of $L$ is also very small in practice, often 2; thus, for a wide range of situations, we have set $L = 5$ in our implementation. Overall, the maximum complexity of computation for this stage, applied to $K$ distorted zones, is $O(KLTr^2)$ where $r$ is the average line thickness of the image $I$. In other words, this stage is linear time complexity for the areas of the distorted zones. Note that the distorted zones, in practice, could intersect, resulting in connected component distorted zones and making the searching areas much smaller.

# 2.4 Experimental results

## 2.4.1 Evaluation metric and protocol

We use *repeatability* criterion to evaluate the performance of our junction detector because this criterion is standard for the performance characterization of local keypoint detectors in CV [Tuytelaars and Mikolajczyk, 2008]. This criterion works as follows. Given a reference image $I_{ref}$ and a test image $I_{test}$ taken under different transformations (e.g., noise, rotation, scaling) from $I_{ref}$, the repeatability criterion signifies that the local features detected in $I_{ref}$ should be repeated in $I_{test}$ with some small error $\epsilon$ in location. We denote $D(I_{ref}, I_{test}, \epsilon)$ as the set of points in $I_{ref}$ that are successfully detected in $I_{test}$ in the sense that for each point $p \in D(I_{ref}, I_{test}, \epsilon)$, there exists at least one corresponding point $q \in I_{test}$ such that *distance*($p$, $q$) $\le \epsilon$. Let $n_r$ and $n_t$ be the number of keypoints detected by one detector from $I_{ref}$ and $I_{test}$, respectively. The repeatability score of this detector
applied to the pair $(I_{ref}, I_{test})$ is computed as follows:

$$r(I_{ref}, I_{test}, \epsilon) = \frac{|D(I_{ref}, I_{test}, \epsilon)|}{Mean(n_r, n_t)} \qquad (2.11)$$

<figure><img src="image_50.png" /><figcaption>Figure 2.11: The evaluation strategy applied to each detector.</figcaption></figure>

Our strategy to perform evaluation in the experiments is described in Figure 2.11. This strategy follows the general characterization protocol for keypoint detection in CV. Particularly, we first apply each detector to the reference images and the test images in each dataset to obtain the reference junctions ($S_r$) and the detected junctions ($S_t$), respectively. Then, we use groundtruth information to compute repeatability scores of this detector from two sets of junctions: $S_r$ and $S_t$. The overall repeatability score of each detector in each experiment is computed as an average score from the repeatability scores of the detector applied for all model symbols and test symbols in each dataset. We vary the value of parameter $\epsilon$ in the range of [1, 8] to obtain a ROC-like curve of the repeatability score.

### 2.4.2 Baseline methods

To compare the proposed system with other methods, we have selected two baseline systems [Liu et al., 1999, Hilaire and Tombre, 2006] dedicated to junction and fork-point detection. The work of [Liu et al., 1999] is dedicated to fork-point detection in handwritten Chinese characters, whereas the work of [Hilaire and Tombre, 2006] is a vectorization-based system for line-drawings. We wish to highlight that although the later work is designed for vectorization, the major contribution in this work is the process of skeleton optimization to correct skeletons and reconstruct junctions. As the implementations of these works are not publicly available, we have developed our own implementation for these two systems<fnref n="1" />. For each system, several running trials have been performed to select the best parameter settings, and only junctions or fork-points are compared with our detected junctions. It is worth noting that we have applied the same pre-processing steps for all three systems and used the same parameter settings in all the experiments.

### 2.4.3 Datasets

The datasets used in the experiments are summarized in Table 2.1, including the final datasets from Symbol Recognition Contest in GREC2011 (SymRecGREC11)<fnref n="2" />, the UMD

<footnote n="1">¹The source codes for the three systems and the demonstration of our junction detector and symbol spotting are publicly available at https://sites.google.com/site/ourjunctiondemo/</footnote>

<footnote n="2">²http://iapr-tc10.univ-lr.fr/index.php/symbol-contest-2011</footnote>
<table><caption>Table 2.1: Datasets used in our experiments</caption><thead><tr><th rowspan="2">No.</th><th rowspan="2">Dataset</th><th rowspan="2">Type</th><th rowspan="2">Noise</th><th colspan="2">#Images</th></tr><tr><th>#References</th><th>#Tests</th></tr></thead><tbody><tr><td>#1</td><td>GREC11</td><td>Line-drawing</td><td>Rotation</td><td>150</td><td>1339</td></tr><tr><td>#2</td><td>GREC11</td><td>Line-drawing</td><td>Scaling</td><td>150</td><td>1200</td></tr><tr><td>#3</td><td>GREC11</td><td>Line-drawing</td><td>Kanungo+Rotation+Scaling</td><td>150</td><td>15000</td></tr><tr><td>#4</td><td>GREC11</td><td>Line-drawing</td><td>Context</td><td>18</td><td>1800</td></tr><tr><td>#5</td><td>SESYD</td><td>Line-drawing</td><td>Low Resolution</td><td>100</td><td>936</td></tr><tr><td>#6</td><td>UMD Logos</td><td>Filled-shape</td><td>Kanungo+Rotation+Scaling</td><td>104</td><td>1272</td></tr></tbody></table>

Logo Database of University of Maryland, Laboratory for Language and Media Processing (LAMP)<fnref n="3" />, and low resolution diagram dataset from SESYD<fnref n="4" />. The SymRecGREC11 dataset is composed of 4 folders, namely setA, setB, setC, and setD with respect to 2500, 5000, 7500, and 1800 test images, respectively. The first three folders are distorted by a mixture of Kanungo noise and geometric transformations (i.e., scaling and rotation), whereas the last dataset is disturbed by context noise (i.e., symbols cropped from full line-drawing images). The UMD Logo Database consists of 104 model logos, which have been used to generate 1272 test images by applying a combination of Kanungo noise and geometric transformations. The low resolution diagram SESYD dataset contains 100 reference images and 936 test images by applying 4 levels of low resolution, corresponding to the scaling factors {1/2, 1/4, 1/8, 1/16}. Consequently, the image resolution is varied from 1700 × 1700 to 100 × 100, and the line thickness is also varied in the range of [2, 18]. These test images are then exported in PNG format in which some blurring effect will be automatically incorporated to these images. In addition, for the evaluation of single parameter changes (i.e., rotation and scaling), we have used 150 model symbols from GREC2011 to generate 1339 test images taken under different levels of rotation (e.g., from 10⁰ to 90⁰) and 1200 test images taken under different scaling factors (i.e., {1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0}).

### 2.4.4 Comparative results

#### 2.4.4.1 Evaluation of rotation and scaling change

In this experiment, the repeatability scores are computed over single parameter changes while the location error is fixed at 4 (pixels). Figure 2.12(a) presents the effect of rotation change for three systems. It can be noticed that the proposed approach far outperforms (almost 25%) the two other systems and that our repeatability scores tend to be quite stable and almost over 88% when varying the rotation parameter. These results confirm that the proposed system is very robust to rotation change. The systems of Liu et al. and Hilaire et al. are theoretically rotation invariant; however, the results reported here show that these systems are less adaptive to rotation change under real-world conditions.

In the context of scaling change, as shown in Figure 2.12(b), the same situation is repeated for the baseline methods, whereas the proposed system still strongly outperforms

<footnote n="3">³http://lampsrv02.umiacs.umd.edu/projdb/project.php?id=47</footnote>

<footnote n="4">⁴http://mathieu.delalandre.free.fr/projects/sesyd/</footnote>
<figure><img src="image_51.png" /><figcaption>Figure 2.12: Repeatability scores of three systems on rotation change (a), and scaling change (b), where location error is set at 4 pixels.</figcaption></figure>

the others. It is noticed, in particular, that the results obtained by the system of Hilaire et al. are significantly degraded when increasing the scaling factor. This degradation could be due to the two baseline systems being quite sensitive to the digitization effect caused by rotating and scaling the input images. In fact, the results reported in the work of [Hilaire and Tombre, 2006] are applied to several line-drawing images that are typically used in the context of vectorization contests, whereas the results reported by [Liu et al., 1999] are applied to handwritten Chinese characters, which are not taken under extreme rotation/scaling changes.

### 2.4.4.2 Evaluation of a mixture of Kanungo noise and rotation/scaling change

We have selected the first three sets, *setA*, *setB*, and *setC*, from the final recognition datasets in GREC2011 to evaluate the performance of three systems under different combination of binary noise and geometric transformations. Some examples of such degradation are shown on Figure 2.14 (a, b). The purpose of this experiment is to justify how well each system can work under different levels of degradation. The results are presented on Figure 2.13 where the proposed system achieves much better results on all three *setA*, *setB*, and *setC*, compared to the systems of Hilaire et al. and Liu et al. On average, the repeatability scores obtained by the proposed system are 15% higher than those of the other systems, especially for the first small range of location errors (e.g., see the first part of the score curve of the proposed system). If we have a more detailed analysis at our results by fixing the location error at 4 (pixels), it can be noticed that the repeatability scores of the proposed system are almost 80% on all three *setA*, *setB*, and *setC*. These obtained results are quite interesting considering a severe degree of degradation applied on these datasets. Under the same conditions, we can see that the scores of two baseline systems are approximately 60%, which is much less than that of the proposed system. These results suggest that the
proposed system can resist to a satisfactory level of degradation composing of common binary noise and geometric transformations.

Two main factors explain these results. First, the polygonization process in the system of Liu et al. and the skeleton segmentation step in the system of Hilaire et al. are rather sensitive to the distortion of contours. Second, the post-process of junction merging using *Criterion A* is quite sensitive to the variation and distortion of the line thickness of foreground objects. The results of our system suggest that the proposed system can satisfactorily resist degradation including common binary noise and geometric transformations. A last noticeable point exhibiting on Figure 2.13 is that there is a little difference of performance obtained on the *setA*, *setB*, and *setC* for all three systems because there is, in fact, no significant difference of degradation among the images of these three test sets.

<figure><img src="image_52.png" /><figcaption>Figure 2.13: Repeatability scores of three systems for *setA*, *setB*, *setC*, and *setD* of GREC2011.</figcaption></figure>
### 2.4.4.3 Evaluation of context noise

Although Kanungo noise and geometric transformations are very common degradation models in DIA, a more realistic type of degradation is known as *context* noise. By definition, context noise is concerned with a type of disturbance caused by background or context information. For this purpose, we have selected *setD* from the final recognition dataset in GREC2011. The test images in this set have been cropped from full line-drawing documents where each reference image could be touched with other context information. The repeatability scores of three systems are reported in Figure 2.13(d). Although the proposed system still outperforms the others, the repeatability scores of all systems are quite low. This finding is attributed to the fact that the images in *setD* are embedded into other context information, resulting in many false positives being detected in this set, as shown in Figure 2.14 (c, d). However, without any prior knowledge about groundtruth information, these false alarms correspond to the mismatches of correct detected junctions missing in the reference images.

<figure><img src="image_53.png" /><figcaption>Figure 2.14: Junctions detected (red dots) by the proposed system for *setC* ((a) and (b)) and *setD* ((c) and (d)) of GREC11. False positives caused by context noise in *setD* are marked by dashed-line boxes.</figcaption></figure>

### 2.4.4.4 Evaluation of the low resolution dataset

In this experiment, we wish to assess the performance of the three systems for very severely low resolution images. We have selected the low resolution diagram SESYD dataset, in which the test images have been generated from the reference images by applying four exponential levels of low resolution corresponding to the scaling factors {1/2, 1/4, 1/8, 1/16} and incorporating a compression scheme (JPG) from gray-scale images. The results of the three systems are shown in Figure 2.15(a), where the proposed
system again performs better than the baseline systems. On average, the proposed system provides 10% and 5% better results than those obtained from the systems of Hilaire et al. and Liu et al., respectively. All three systems perform quite well on the first levels of low resolution, but their performance rapidly degrades for the later levels of low resolution. This behavior is mainly due to the loss of much of the original information, especially finer features, when reducing the resolution.

<figure><img src="image_54.png" /><figcaption>Figure 2.15: Repeatability scores of three systems for SESYD low resolution dataset (a), and the UMD logo dataset (b).</figcaption></figure>

### 2.4.4.5 Evaluation of the filled-shape and non-uniform stroke dataset

In this experiment, we want to justify how different kind of images, such as filled-shape and non-uniform stroke images, can impact the performance of the three systems. For this purpose, we have selected the UMD logo dataset, which typically composes of filled-shape and non-uniform stroke objects at a hard level. The repeatability scores are presented in Figure 2.15(b). Even though the skeleton-based representation for such kind of images is not perfect, the obtained results are encouraging. The system of Hilaire et al. achieves the lowest scores because of a lower number of outputted junctions. As the line thickness of the filled shapes is greatly varied compared to the typical line-drawings, many short skeleton segments are produced. Consequently, few long skeleton segments are retained, and thus the number of detected junctions is rather limited in the system of Hilaire et al. Our system also produces a limited number of junctions, even less than that of the Hilaire's system, but it still noticeably outperforms the system of Hilaire et al. and almost gives the same results as those of the system of Liu et al. These results confirm the accuracy of the detected junctions of the proposed system. Some visual results of the proposed system, applied to the logo images and Chinese characters, are shown in Figure 2.16.
<figure><img src="image_55.png" /><figcaption>Figure 2.16: Junctions detected (red dots) by the proposed system for few Chinese characters and logo images.</figcaption></figure>

### 2.4.4.6 Evaluation of the built-in aspects

In addition to the evaluations discussed above, we have investigated several additional trials to understand the behavior of the proposed approach at the system level. In particular, we wish to present a detailed analysis of the impact of the stage of determination of ROS and the computation time of our junction detector. Regarding the first aspect, we have computed the repeatability scores for the setC up to the stage of dominant point detection over three different scenarios: the use of ROS based on local line thickness (i.e., the ROS at given a point *p* is set as the local line thickness at *p*), the use of ROS proposed by [Teh and Chin, 1989], and the use of ROS proposed by our approach. The results are presented in Table 2.2. It can be seen that our method achieves much better results than the others (by almost 23%). The results linked to the ROS proposed by Teh-Chin are quite low because, as we have discussed in Section 2.1.2, the Teh-Chin's ROS determination step is sensitive to digitization effects, whereas in this dataset, the noise applied to these images is quite severe, distorting their shapes. The results shown in Table 2.2 also reveal that line thickness could be a good feature to estimate local scales.

<table><caption>Table 2.2: Comparison of the dominant point detection rates for three scenarios.</caption><thead><tr><td>Dominant point detection mode</td><td>Repeatability Score</td></tr></thead><tbody><tr><td>With ROS of the proposed method</td><td>67.5 %</td></tr><tr><td>With ROS based on line thickness</td><td>44.3 %</td></tr><tr><td>With ROS proposed by Teh-Chin</td><td>21.3 %</td></tr></tbody></table>

We also performed an additional experiment to study the impact of the two parameters
$E_{min}$ and $E_{max}$ in the stage of ROS determination. For this purpose, we vary the values of $E_{min}$ and $E_{max}$, and compute the repeatability score of the proposed system for the setC up to the stage of dominant point detection. The obtained results are presented in Figure 2.17. These results show that the detection rate is quite stable (e.g., varying in the range of [63, 68]) given various settings of $E_{min}$ and $E_{max}$. This is expected as It is noted that the proposed system achieves the repeatability score of 67.5% in Table 2.2 with respect to the following setting: $E_{min} = 1.3$ and $E_{max} = 1.8$.

<figure><img src="image_56.png" /><figcaption>Figure 2.17: Impact of the parameters $E_{min}$ and $E_{max}$: the repeatability score is computed on the setC.</figcaption></figure>

For time complexity evaluation, we provide in Table 2.3 some information about the processing time (excluding the pre-process step) of three systems applied on several images with different sizes. The processing time has been recorded on our specific computer configuration: Intel(R) Core(TM) i5 CPU 2.4 GHz, RAM 2.4 GB, Windows 8.

<table><caption>Table 2.3: Report of the processing time (ms) and the number of detected junctions (in brackets).</caption><thead><tr><th rowspan="2">System</th><th colspan="3">Image size (Width × Height)</th></tr><tr><th>900 × 984</th><th>3600 × 3938</th><th>2100 × 4433</th></tr></thead><tbody><tr><td>Our system</td><td><strong>16.0 (67)</strong></td><td><strong>187.0 (79)</strong></td><td><strong>140.0 (105)</strong></td></tr><tr><td>Hilaire et al.</td><td>110.0 (89)</td><td>297.0 (147)</td><td>265.0 (146)</td></tr><tr><td>Liu et al.</td><td>563.0 (82)</td><td>14953.0 (157)</td><td>4078.0 (174)</td></tr></tbody></table>

In general, the system of Liu et al. is subjected to a high computation load because distorted skeleton correction using *Criterion A* is very time-consuming. The system of Hilaire et al. seems to provide a reasonable level of processing time because the criterion
to merge two discrete primitives is somewhat similar in spirit to *Criterion A* but with the elimination of much of the redundant computation. Our system works most efficiently, not only for the cases of the several images reported in Table 2.3 but also throughout the extensive experiments we have performed. It is also noted that the number of detected junctions (in brackets) provided by our system is much smaller than those outputted by the other systems. We provided few illustrative examples of the detected junctions of our approach applied to different kinds of images as shown in Figures {2.18, 2.19, 2.20, 2.21}.

## 2.5 Discussion

This chapter presents a new approach for junction detection and characterization in line-drawing images. The main contribution of this work is three-fold. First, a new algorithm for the determination of the region of support is presented using the linear least squares technique. The crossing-points, in combination with the dominant points detected from median lines, are treated as candidate junctions. Next, using these candidate junctions, an efficient algorithm is proposed to detect and conceptually remove all distorted zones, retaining reliable median line segments only. These line segments are then locally characterized to construct the topological representations of the crossing zones. Finally, a novel junction optimization algorithm is presented, yielding accurate junction localization and characterization. The proposed approach is extremely robust to common geometry transformations and can resist a satisfactory level of noise/degradation. Furthermore, it works very efficiently in terms of time complexity and requires no prior knowledge of the document content. The proposed method is also independent on any vectorization systems. All of these prominent features of the proposed approach have been validated relative to other baseline methods by our extensive experiments.

In addition to these advantages, the proposed approach has several shortcomings. First, as this approach is dedicated to working with line-like primitives, its performance would be degraded if applied to filled-shape objects, such as logo images. In addition, the junction optimization process could lead to some difficulties in correctly interpreting the junction position as originally produced by craftsmen. However, although this point is valid for some specific domains of exact line-drawing representation, such as vectorization, we are interested in detecting local features that would be useful to addressing the problem of large-scale document indexing and retrieval. In this sense, a low rate of false positives in the final results is not problematic. A last noticeable point is that the detected junctions could be used in combination with additional features (e.g., end-points, *isolated* straight lines, arcs, and circles) to obtain the complete representation of a graphical document image. The detected junctions can be also used to address the problem of vectorizing the line-drawings.
<figure><img src="image_57.png" /><figcaption>Figure 2.18: Detected junctions for a synthetic symbol with different levels of noise.</figcaption></figure>

<figure><img src="image_58.png" /><figcaption>Figure 2.19: Detected junctions for a real musical score image.</figcaption></figure>
<figure><img src="image_59.png" /><figcaption>Figure 2.20: Detected junctions for part of an electronical image.</figcaption></figure>
<figure><img src="image_60.png" /><figcaption>Figure 2.21: Detected junctions for a mechanical-text image.</figcaption></figure>
# Chapter 3

# Application to symbol localization

In this chapter, an application of symbol localization in line-drawing images is developed to demonstrate that our junction detector is robust and discriminated enough to be used in the context of object localization. The experimental results, applied to several public datasets, show that our system is very time- and memory-efficient. Our precision and recall results in terms of symbol localization highlight that we outperform other methods in the literature on this problem.

## 3.1 Introduction

A common problem of any symbol processing systems, recognition or spotting, is localization or detection of the symbols. Symbol localization can be defined as the ability of a system to localize the symbol entities in the complete documents. It could be embedded in the recognition/spotting method or works as a separated stage in a two-step system [Qureshi et al., 2008]. The approaches used for localization are similar for recognition and spotting. All systems rely first on a primitive extraction step (e.g., connected components, loops, key-points, lines, etc.). These systems differ mainly in the way that the detected primitives are processed, using machine learning or retrieval and indexing techniques. Different approaches have been investigated in the literature to deal with the localization problem.

One of the earliest approach employed in many systems is subgraph matching. Graph is a very effective tool to represent line drawings. Attributed Relational Graphs (ARGs) can be used to describe the primitives, their associated attributes and interconnections. However, subgraph isomorphism is known to be a NP-hard problem, making it difficult
to use the graph for large images and document collections, despite the approximate solutions of subgraph isomorphism developed in the literature [Messmer and Bunke, 1996, Bodic et al., 2009]. In addition, subgraph isomorphism remains very sensitive to the robustness of the feature extraction step, as any wrong detection can result in strong distortions in the ARGs.

An alternative approach to subgraph matching is the use of a "triggering" mechanism. Such a system looks for some specific primitives in line drawing images and triggers a matching process at the symbol level within the Regions of Interest (ROIs) around these primitives. The system in [Nguyen et al., 2009] is a typical example. In this work, given a query symbol, the keypoints (i.e., Difference of Gaussian features) and its corresponding vocabularies are computed and used to find the matched keypoints from the database documents. For each pair of two matched keypoints, the local scale and orientation extracted at the keypoint in the query symbol are used to generate the ROI in the document that probably contains the instance of the symbol. Because the number of detected keypoints would be very large and the local scale computed at each keypoint could be far from satisfaction, the ROI extraction step is thus fragile and time-consuming. Triggering mechanisms have been also developed from graph-based representations, as in [Rusinol and LLados, 2006, Qureshi et al., 2008]. These proposed systems work from the ARGs, where the structures and attributes of the graphs are exploited to identify the ROIs without recognizing the symbols. Triggering-based localization is very sensitive to robustness of the mechanism in that any missed detection at the triggering level will result in the failure of symbol localization.

In the other place, a different approach to deal with object localization is framing [Dosch and LLados, 2004, Kong et al., 2011, Dutta et al., 2011]. These techniques involve the decomposition of the image into frames (i.e., tiles, buckets, windows) in which the frames could be overlapped [Kong et al., 2011] or disjointed [Dosch and LLados, 2004]. Local signatures are computed from the primitives contained in the frames and matched to identify the candidate symbols. The size of the frames can be determined based on the symbol models [Dosch and LLados, 2004, Kong et al., 2011] or set at different resolutions [Dutta et al., 2011]. In this way, framing is not scale invariant as the size of the frames cannot be dynamically adapted. The position of the frames can be set with a grid [Dosch and LLados, 2004, Dutta et al., 2011] or by sliding [Kong et al., 2011]. Sliding could be performed by steps to reduce the entire processing time [Kong et al., 2011], as any computations with overlapping would be subjected to a polynomial complexity.

A part from the before-mentioned approaches, a recent framework for symbol localization is the use of geometry consistency checking as presented in [Nayef and Breuel, 2011, Jain and Doermann, 2012, Rusinol et al., 2013]. Such a method concerns a pipeline of decomposing a graphical document into a set of primitives, matching the primitives of the model against the test image, and checking the geometry constraints among the matches. Different techniques for geometric verification have been applied in these works. [Nayef and Breuel, 2011] proposed to use a branch and bound algorithm to search for a transformation that maps a maximal subset of the matches between the primitives of the model and the test images. [Rusinol et al., 2013] employed a classical technique RANSAC [Fischler and Bolles, 1981] to achieve this goal. [Jain and Doermann, 2012] incorporated the use of orientation information of the detected features to prune the matches, following
the verification step of pair-wise angles between every two triangles in the model and test images.

The common drawback of all these methods is the computation complexity of the geometric verification process. In this application, we mainly aim at demonstrating that our contribution on junction detection can be used for symbol localization by incorporating some addition processing steps. In that concern, we show that (1) the detected junctions are useful to deal with the problem of symbol localization, and (2) these junctions support the process of geometry verification in a very efficient way. Particularly, our system is composed of four main stages, each of which is briefly described below with respect to the Figure 3.1.

* In the first stage, the junction points are detected and characterized into different types such as T-, L-, and X-junctions.

* The second stage decomposes a document image into a set of smooth primitives, composing of isolated shapes (e.g., isolated circles and straight lines) and curve segments bounded between either two junctions or a junction and an end-point. These primitives are then associated with a new set of keypoints including Line-, Arc-, and Circle-keypoints. The obtained keypoints, in combination with the junction points and end-points, form a complete and compact representation of graphical documents.

* In the third stage, keypoint matching is performed to find the correspondences among the keypoints of the query and those of database documents.

* Finally, geometry consistency checking is applied to the obtained matches using a new and efficient algorithm, which is designed to work on our specific keypoints.

<figure><img src="image_61.png" /><figcaption>Figure 3.1: Overview of our symbol localization system.</figcaption></figure>

As the first stage is accomplished by simply applying our contribution as presented in chapter 2, it will not be detailed in the following sections. Instead, we directly describe hereafter the last three stages.

## 3.2 Document decomposition

We use the detected junctions to decompose a document image into a set of smooth primitives. Here, we define the smooth primitives as those composing of isolated shapes (e.g., isolated circles and straight lines) and curve segments bounded between either two junctions or a junction and an end-point. This definition is derived based on the fact that after the process of junction detection, every median line segment bounded between two
junctions are sufficiently smooth. Otherwise, some new junction points are likely to be detected on this segment. In this work, we restrict the smooth primitives to three kinds of segment: straight line segment, arc segment, and circle. These basic-shape primitives could be derived using the linear least squares (LLS) fitting technique as follows (see Figure 3.2):

* For each smooth primitive $P$, we try first to fit $P$ to a straight line segment by comparing the average distance error to a fixed threshold (e.g., 1.5 pixels in our implementation). The distance error is simply computed as the Euclidean distance from each point of $P$ to the fitted straight line.

* If $P$ is not fitted to a straight line, circle fitting will be performed next. This is also accomplished by comparing the average distance error to a threshold but the distance error is now computed as the Euclidean distance from each point of $P$ to the fitted circle.

* If $P$ is fitted to a circle, it will be further classified as an arc primitive (e.g., opened curve) or a circle (e.g., closed curve).

<figure><img src="image_62.png" /><figcaption>Figure 3.2: Basic primitive decomposition and description.</figcaption></figure>

Next, each type of these primitives is characterized as a specific structural keypoint as follows:

* A straight line primitive is represented by a triple {$p_L, p_{L1}, p_{L2}$} corresponding to the middle point and two extremity points, respectively. A triple {$p_L, p_{L1}, p_{L2}$} is regarded as a Line-type keypoint or *L-keypoint*.

* An arc primitive is represented by {$p_A, p_{A1}, p_{A2}$} with the same meaning as that of a straight line primitive. A triple {$p_A, p_{A1}, p_{A2}$} is regarded as an Arc-type keypoint or *A-keypoint*. It is noted that the characterization of an A-keypoint is proceeded in the same spirit as that of a junction whose two arms are $p_A p_{A1}$ and $p_A p_{A2}$.
* A circle primitive is represented by {$p_C, r_C$} corresponding to its centroid and radius. A couple {$p_C, r_C$} is regarded as a Circle-type keypoint or C-keypoint.

For completeness, we call the junction points as J-keypoints and end-points as E-keypoints. The description of a J-keypoint is performed using the same process of junction characterization, while the description of the E-keypoints is no needed. To this end, a document image is completely represented by a set of structural keypoints, composing of L-keypoints, A-keypoints, C-keypoints, J-keypoints, and E-keypoints. Figure 3.3 shows a decomposition of a document image into a set of structural keypoints.

<figure><img src="image_63.png" /><figcaption>Figure 3.3: Keypoint-based representation of a simple document image.</figcaption></figure>

## 3.3 Keypoint matching

In the previous stage of document decomposition, a document image is completely represented by a set of structural keypoints. In this stage, keypoint matching will be performed to establish the correspondences between the keypoints of the query *Q* and those of the database document *D*. Keypoint matching is independently processed for each type of keypoint. Particularly, matching of the L-, A-, C-, and E-keypoints is proceeded very simply as outlined below:

* An E-keypoint (*resp.* C-keypoint) is always matched with any other E-keypoints (*resp.* C-keypoints).

* A L-keypoint is matched with another L-keypoint if the two extremity points of one L-keypoint are matched with those of the other L-keypoint (Figure 3.4 (a)). It is noted that the extremity points of one L-keypoint could be the E-keypoints and J-keypoints. Therefore, matching of two extremity points is also performed in the same manner as keypoint matching.

* An A-keypoint is matched with another A-keypoint if the difference in direction of the two keypoints is not significant and the two extremity points of one A-keypoint are matched with those of the other A-keypoint (Figure 3.4 (b)).
The matching process of the J-keypoints is often performed by pair-wise matching between the angles of the two corresponding junction points [Xia, 2011]. Particularly, given two J-keypoints characterized as {$p, s_p, \{\theta_i^p\}_{i=0}^{m_p-1}$} and {$q, s_q, \{\theta_j^q\}_{j=0}^{m_q-1}$}, the information of junction location and junction scale is used to quickly refine the matches to be described later, and the rest is used to compute a similarity score, $C(p, q)$, of matching two junctions $p$ and $q$ as follows:

$$C(p, q) = \max_{i,j} \left\{ \frac{1}{H} \sum_{k=0}^{h-1} D(\theta_{(i+k) \bmod m_p}^p, \theta_{(j+k) \bmod m_q}^q) \right\} \quad (3.1)$$

where $h = min(m_p, m_q)$, $H = Max(m_p, m_q)$, and

$$D(\theta_i^p, \theta_j^q) = \begin{cases} 1, & \text{if } |\theta_i^p - \theta_j^q| \le \theta_{\text{thres}} \\ 0, & \text{for otherwise.} \end{cases} \quad (3.2) \quad (3.3)$$

<figure><img src="image_64.png" /><figcaption>Figure 3.4: The matching process of L-keypoints (a) and A-keypoints (b).</figcaption></figure>

<figure><img src="image_65.png" /><figcaption>Figure 3.5: An example of context distortion of the detected junctions: a 3-junction p in (a) is distorted as a 5-junction q in (b).</figcaption></figure>

The similarity score $C(p, q)$ is in the range [0, 1] and $\theta_{\text{thres}}$ is an angle difference tolerance. Most of the time, two J-keypoints are matched if their similar score is higher than
a threshold. However, in some specific domains taking object localization for example, a query object or symbol is often connected to different context information appearing in a document. Figure 3.5 provides an example where the two instances of the query symbol are touched resulting in some context distortion of the detected junctions (e.g., a 3-junction **p** is distorted as a 5-junction **q**). In such cases, using the similarity score $C(p, q)$ could be too restricted to find corresponding junctions.

We therefore relax the junction matching step by introducing a new constraint as follows. Two J-keypoints $p$ and $q$ are matched if an *inclusion test* is hold for these two junctions. Here, we consider that $p$ is included in $q$ if there are exact $m_p - 1$ angle matches between the angles of $p$ and $q$. This implies: $C(p, q) * \text{Max}(m_p, m_q) = m_p - 1$.

## 3.4 Geometry consistency checking

Given a query symbol $Q$ and a database document $D$, their keypoints are first detected and matched as described in the previous sections. The obtained matches are finally verified by checking geometry consistency. This step will remove false matches and cluster the remaining matches into different clusters, each of which indicates an instance of the query symbol. Concerning this problem of geometry consistency checking, two main strategies are often exploited in the literature. A brief review of these two strategies is given in the following.

The first strategy treats data (i.e., the matches) in a top-down way. One typical technique belonging to this strategy is known as RANSAC (RANdom Sample Consensus) [Fischler and Bolles, 1981]. The key idea of RANSAC is to randomly select $k$ matches for estimating a transformation model (typically an affine transformation and thus $k = 2$ or 3). The model is then assigned with a confidence factor, which is calculated as the number of matches fitting well to this model. Next, these steps are repeated a number of times to find the model with the highest confidence. RANSAC is often used to find a *single transformation model* between two images with a high degree of accuracy of the derived parameters, provided that the ratio of inliers and outliers occurring in the data is sufficiently high ($\ge 50\%$). However, when this is not the case, it is difficult to use RANSAC. In addition, RANSAC would be time-consuming because the number of iterations is often large to ensure that an optimal solution could be found. Precisely, the computation complexity is $O(MN)$ where $M$ is the number of iterations and $N$ is the size of the data. A family of advanced algorithms based on the RANSAC technique is reported in [Choi et al., 2009] where the accuracy and robustness are thoroughly investigated. The computation complexity of these algorithms is also evaluated as the trade-off between accuracy and robustness.

The second strategy treats data in a bottom-up manner by performing a voting process starting from all data points, and then finding the parameters (typically composing of 4 parameters: orientation, scaling, and $x, y$-translation) corresponding to the dense density areas of support. One typical technique falling this strategy is known as Generalized Hough Transform (GHT) [Ballard, 1981]. GHT is most commonly used for the cases where *multiple transformation models* are presence in the data. It is less accurate than RANSAC in terms of parameter estimation but very robust to noise even if a large number of outliers are presence. However, because GHT requires a process of parameter quantization, it is
subjected to very high cost of memory space $O(M^4)$ and processing time $O(N^2)^1$, and
is sensitive to the quantization of the parameters. A comparative evaluation of the GHT
techniques can be found in [Kassim et al., 1999]. In this work, different extensions of
the GHT method are presented to reduce the memory requirement and computational
complexity.

In our case, as each keypoint of the query symbol is often occurred in a database
document with a high frequency, the outliers are thus significantly higher than the in-
lier matches. In addition, as multiple instances of a query symbol can be appeared in a
database document, it is therefore not a good idea to use some techniques like RANSAC
or GHT because of the aforementioned weaknesses. We therefore present, below, an effi-
cient algorithm to deal with the problem of geometry consistency checking. The proposed
algorithm incorporates the advantages of both RANSAC and GHT while avoiding their
weaknesses. It exploits the information extracted from the matches of L- and A-keypoints
to speed up the process of estimating the affine transformation models.

<caption>**Algorithm 1** Estimation of the affine transformation models</caption>

**Input:** Two set of structural keypoints of a query $Q$ and a document image $D$, and a match list $T$ between the keypoints of $Q$ and $D$
**Output:** A set of affine transformation models ($F_{out}$)
$K \leftarrow$ sort the L- and A-keypoints of $Q$ in descending order of the primitive's length
$m \leftarrow |K|$
$i \leftarrow 1$
$F_{out} \leftarrow \emptyset$
**while** $i \le m/2$ **do**
  $p \leftarrow K_i$
  **for** each match in $T$ between $p \in Q$ and $q \in D$ **do**
    $F \leftarrow$ solve the two linear equations formed by the extremity points of $p$ and $q$
    $n_F \leftarrow$ count the matches fitting to $F$
    **if** $n_F > n_{thres}$ **and** $F \notin F_{out}$ **then**
      $F_{out} \leftarrow F_{out} + \{F\}$
    **end if**
  **end for**
  $i \leftarrow i + 1$
**end while**

Our method of geometry consistency checking is outlined in Algorithm 1. The basic idea is to directly estimate a geometry model $F$ (i.e., an affine transformation) based on every match formed by a pair of either two L-keypoints or two A-keypoints. This idea is inspired by the fact that a pair of two matched lines (or arcs) provides us with 4 parameters (i.e., orientation, scaling, and $x, y$-translation) of an affine transformation $F$. As a result, we need only one match to estimate a model $F$ other than two matches as the cases of GHT and RANSAC. Next, we apply the transformation $F$ to all the keypoints detected on the query $Q$, resulting in a new set of projected points on the database document $D$. If there is a sufficiently large overlap among the projected points and the keypoints of $D$ matched

<footnote n="1">¹N and M, are the number of matches and sampling bins, respectively.</footnote>
with those of $Q$, the transformation $F$ is accepted and used to localize the position of the
corresponding instance of $Q$ on $D$. Here, we consider two points are overlapped if the
distance between them is less than a threshold. Alternatively, we can consider two points
are overlapped if they are positioned in a local window of the size $\epsilon_{dist} \times \epsilon_{dist}$. This can be
efficiently done by using a 2D lookup table. The memory complexity is thus linear to the
image size. It was found empirically that the parameter $\epsilon_{dist} \in [10, 20]$ is a common setting.
It is also noted that we can set $\epsilon_{dist}$ to the minimum distance between two keypoints of
the query document $Q$ projected on the database document $D$ to ensure no mismatch in
our geometry consistency checking step.

Regard to the computation complexity, as the number of the L- and A-keypoints of $Q$
is quite small, and few real computations are needed, the proposed method is very time-
efficient. Particularly, the computation complexity of the proposed method is limited up
to a linear order $O(\epsilon_{dist}^2 N_1 N_2)$, where $N_1$ is the number of keypoints of $Q$, and $N_2$ is the
number of matches corresponding to the L- and A-keypoints of $Q$. In addition, with a bit
prior knowledge of the dataset, we can quickly prune a large number of matches by setting
the lower and upper scales for the query symbol. In this way, the local scales associated
to the keypoints are used to prune the matches. It is also noted that there is no need to
process all the L- and A-keypoints of $Q$. In our experiments, we first sort $m$ L- and A-
keypoints of $Q$ in a decreasing order of their length and then choose the first $m/2$ keypoints
to be processed. Since the primitives having higher length would be less distorted by the
transformation, the estimation of the parameters is thus more robust. For each accepted
model $F$, we obtain an instance of the query. Therefore, multiple instances of the query
are thus successfully detected with respect to the number of accepted models. Figure 3.6
demonstrates the result of applying this step of geometry consistency checking.

<figure><img src="image_66.png" /><figcaption>Figure 3.6: Geometry consistency checking: (a) the matches before checking; (b) the matches after checking.</figcaption></figure>

## 3.5 Experimental Results

For performance evaluation of the proposed approach, we selected the latest dataset
for symbol spotting in GREC2011² [Valveny et al., 2011]. The detail of this dataset is

<footnote n="2">²http://iapr-tc10.univ-lr.fr/index.php/final-test-description</footnote>
described on Table 3.1. For performance evaluation, we selected the evaluation metric in [Rusinol and Llados, 2009b] in order to make the comparison with other methods [Dutta et al., 2013a, Dutta et al., 2013b]. This metric includes precision (P), recall (R) and Fscore computed as follows.

$$P = \frac{S_{Int}}{S_{Ret}}, R = \frac{S_{Int}}{S_{GT}}, Fscore = 2 \cdot \frac{P \cdot R}{P + R}$$

Where $S_{Int}$ is the sum of intersection areas between the bounding boxes retrieved by the spotting system and ground-truth, $S_{Ret}$ is the sum of areas of the bounding boxes retrieved by the spotting system, and $S_{GT}$ is the sum of areas of the bounding boxes in ground-truth. It is worth mentioning that each bounding box ($B_{gt}$) in ground-truth is counted at most once. Typically, $B_{gt}$ will be marked as already considered if there exists a bounding box ($B_{ret}$) retrieved by the spotting system such that the ratio of the intersection of their areas to the union of their areas exceeds a given threshold (e.g., 75% in our experiments). This is needed to avoid biased scores caused by multiple detections of a same symbol at a same location. Furthermore, we have set a strict constraint for our system in that the ratio of overlapping area to union area of any two retrieved bounding boxes is always less than a threshold (10% in our experiments).

<table><caption>Table 3.1: Dataset used for symbol spotting in GREC2011.</caption><thead><tr><th>Test Set</th><th>Models</th><th>Images</th><th>Queries</th><th>Symbols</th><th>Noise</th><th>Image Size</th></tr></thead><tbody><tr><td>Elec1.</td><td>21</td><td>20</td><td>118</td><td>246</td><td>Ideal</td><td rowspan="2">Min:1700x1600</td></tr><tr><td>Elec2.</td><td>21</td><td>20</td><td>127</td><td>274</td><td>Level 1</td></tr><tr><td>Elec3.</td><td>21</td><td>20</td><td>114</td><td>237</td><td>Level 2</td><td rowspan="2">Max:4400x2100</td></tr><tr><td>Elec4.</td><td>21</td><td>20</td><td>156</td><td>322</td><td>Level 3</td></tr><tr><td>Archi1.</td><td>16</td><td>20</td><td>247</td><td>633</td><td>Ideal</td><td rowspan="2">Min:2300x2500</td></tr><tr><td>Archi2.</td><td>16</td><td>20</td><td>245</td><td>597</td><td>Level 1</td></tr><tr><td>Archi3.</td><td>16</td><td>20</td><td>245</td><td>561</td><td>Level 2</td><td rowspan="2">Max:5400x2900</td></tr><tr><td>Archi4.</td><td>16</td><td>20</td><td>249</td><td>593</td><td>Level 3</td></tr></tbody></table>

<table><caption>Table 3.2: Experimental results of our system (%).</caption><thead><tr><th>Test Set</th><th>Precision</th><th>Recall</th><th>F-Score</th><th>Max<sub>Overlap Union</sub></th><th>Mean Time (ms)</th></tr></thead><tbody><tr><td>Elec1.</td><td>0.85</td><td>0.84</td><td>0.84</td><td>0.00 %</td><td>723.38</td></tr><tr><td>Elec2.</td><td>0.86</td><td>0.79</td><td>0.82</td><td>0.00 %</td><td>677.83</td></tr><tr><td>Elec3.</td><td>0.92</td><td>0.81</td><td>0.86</td><td>5.76 %</td><td>655.39</td></tr><tr><td>Elec4.</td><td>0.80</td><td>0.81</td><td>0.80</td><td>5.76 %</td><td>1097.05</td></tr><tr><td>Archi1.</td><td>0.91</td><td>0.96</td><td>0.93</td><td>6.69 %</td><td>1521.47</td></tr><tr><td>Archi2.</td><td>0.91</td><td>0.92</td><td>0.92</td><td>9.92 %</td><td>1341.90</td></tr><tr><td>Archi3.</td><td>0.94</td><td>0.89</td><td>0.92</td><td>9.92 %</td><td>1605.80</td></tr><tr><td>Archi4.</td><td>0.91</td><td>0.90</td><td>0.90</td><td>10.63 %</td><td>1409.88</td></tr></tbody></table>

The detailed results of our system are reported on Table 3.2 including precision, recall, Fscore, maximum ratio of overlapping area to union area of any two bounding boxes re-
trieved by our system, and mean processing time of a complete query. General speaking, the proposed system achieves quite good results for both detection and accuracy scores. On average, the Fscore(s) of the proposed system are 0.83 and 0.92 for electrical and architectural datasets, respectively. In addition, these scores are obtained under a very small overlap of the detected bounding boxes. It is noted that the results obtained on the architectural dataset are much better than those on the electrical dataset. The reason lies in the fact that in the electrical dataset, more query symbols are used and many query symbols looks very similar making them difficult to be correctly distinguished. The processing time is calculated as the mean time of the whole process (i.e., both the online and offline phases), subjected to our specific computer configuration: Intel(R) Core(TM) i5 CPU 2.4GHz, RAM 2.4 GB, Windows XP.

The last experiment is performed using the SESYD database<fnref n="3" /> [Delalandre et al., 2010]. The detail of this dataset is summarized on Table 3.3. This time, we measure the processing time for the online phase only. That is, we compute the query time for returning the full list of the detected symbol entities given a query symbol. Table 3.4 reports the obtained results of our system. It can be seen that we obtain very good results in terms of both detection rate and recall. On average, the proposed system achieves the precision of 0.92 and the recall of 0.95, resulting in the F-score of 0.93. The mean query time is just 0.3 (s).

In order to make some comparative results, we provide on Table 3.5 the results of some recent symbol localization systems. It is clear that our system much outperforms all these baseline methods. It is worth mentioning that the system of [Dutta et al., 2013b] gives the lowest processing time (i.e., 0.7 (s)) because it was integrated with an indexing hashing-based scheme. Our system requires 0.3 (s), on average, to perform a query without using any indexing methods.

<table><caption>Table 3.3: The detail of the SESYD (floorplans) dataset.</caption><thead><tr><td>Test Set</td><td>Images</td><td>Models</td><td>Symbols</td><td>Noise</td><td>Image Size</td></tr></thead><tbody><tr><td>floorplans16-01</td><td>100</td><td>16</td><td>2671</td><td>None</td><td>6775 × 2858</td></tr><tr><td>floorplans16-02</td><td>100</td><td>16</td><td>2488</td><td>None</td><td>3059 × 3341</td></tr><tr><td>floorplans16-03</td><td>100</td><td>16</td><td>2661</td><td>None</td><td>2218 × 2475</td></tr><tr><td>floorplans16-04</td><td>100</td><td>16</td><td>3251</td><td>None</td><td>2056 × 1837</td></tr><tr><td>floorplans16-05</td><td>100</td><td>16</td><td>2148</td><td>None</td><td>2596 × 2313</td></tr><tr><td>floorplans16-06</td><td>100</td><td>16</td><td>2068</td><td>None</td><td>2352 × 2507</td></tr><tr><td>floorplans16-07</td><td>100</td><td>16</td><td>3898</td><td>None</td><td>5498 × 2961</td></tr><tr><td>floorplans16-08</td><td>100</td><td>16</td><td>2260</td><td>None</td><td>3026 × 2967</td></tr><tr><td>floorplans16-09</td><td>100</td><td>16</td><td>3948</td><td>None</td><td>4307 × 1893</td></tr><tr><td>floorplans16-10</td><td>100</td><td>16</td><td>2653</td><td>None</td><td>4349 × 2227</td></tr></tbody></table>

Figure 3.7 shows an example of our symbol localization system where the query is perfectly localized even though it is embedded into a complicated database document. There are, however, some queries as shown in Figure 3.8 that the proposed system fails to detect the symbol "outlet" because of a very limited number of keypoints detected on the instance of the symbol on the database document. Besides, we showed in Figure 3.9

<footnote n="3">³http://mathieu.delalandre.free.fr/projects/sesyd/symbols/floorplans.html</footnote>
<table><caption>Table 3.4: The results of our system for the SESYD (floorplans) dataset.</caption><thead><tr><th>Test Set</th><th>Precision</th><th>Recall</th><th>F-score</th><th>Mean time (ms)</th></tr></thead><tbody><tr><td>floorplans16-01</td><td>0.97</td><td>0.95</td><td>0.96</td><td>339.5</td></tr><tr><td>floorplans16-02</td><td>0.89</td><td>0.98</td><td>0.93</td><td>284.2</td></tr><tr><td>floorplans16-03</td><td>0.87</td><td>0.91</td><td>0.89</td><td>242.3</td></tr><tr><td>floorplans16-04</td><td>0.93</td><td>0.95</td><td>0.94</td><td>418.6</td></tr><tr><td>floorplans16-05</td><td>0.92</td><td>0.98</td><td>0.95</td><td>112.3</td></tr><tr><td>floorplans16-06</td><td>0.90</td><td>0.97</td><td>0.93</td><td>161.1</td></tr><tr><td>floorplans16-07</td><td>0.96</td><td>0.97</td><td>0.96</td><td>621.2</td></tr><tr><td>floorplans16-08</td><td>0.94</td><td>0.95</td><td>0.94</td><td>226.8</td></tr><tr><td>floorplans16-09</td><td>0.90</td><td>0.98</td><td>0.94</td><td>406.8</td></tr><tr><td>floorplans16-10</td><td>0.88</td><td>0.81</td><td>0.84</td><td>252.3</td></tr><tr><td><strong>Average</strong></td><td><strong>0.92</strong></td><td><strong>0.95</strong></td><td><strong>0.93</strong></td><td><strong>306.5</strong></td></tr></tbody></table>

<table><caption>Table 3.5: Comparison of recent methods for symbol localization on the SESYD (floorplans-01) dataset.</caption><thead><tr><th>System</th><th>Precision</th><th>Recall</th><th>F-score</th><th>Mean time (s)</th></tr></thead><tbody><tr><td>Our system</td><td>0.97</td><td>0.95</td><td>0.96</td><td>0.34</td></tr><tr><td>[Dutta et al., 2013a]</td><td>0.62</td><td>0.95</td><td>0.74</td><td>0.57</td></tr><tr><td>[Dutta et al., 2013b]</td><td>0.41</td><td>0.82</td><td>0.52</td><td>0.07</td></tr><tr><td>[Nguyen et al., 2009]</td><td>NA</td><td>NA</td><td>0.82</td><td>NA</td></tr></tbody></table>

another example that the system correctly detects multiple instances of the symbol "sofa1" even these detections take part of a different symbol (i.e., "table2"). This further confirms the interesting results of our system.

## 3.6 Discussion

We have presented an application to symbol localization in line-drawing images using junction feature and geometry consistency checking. This system proves that our contribution on junction detection can be used in the context of object detection and localization by incorporating some feature extraction steps at primitive level. As the junction detectors is robust and accurate, the obtained primitives are stable under the different contexts of the documents. These primitives are used to support object matching, geometry consistency checking, and object localization. In that sense, this highlights that our detector can support and to be combined in an efficient way within a vectorization process. The experimental results in terms of symbol localization confirm the advantages of the system for both efficiency and accuracy.
<figure><img src="image_67.png" /><figcaption>Figure 3.7: An example of symbol localization: a query symbol (left) and the detected instances of the query (red bounding boxes).</figcaption></figure>

<figure><img src="image_68.png" /><figcaption>Figure 3.8: The system fails to detect the symbol "outlet" (left) due to the omission of the E-keypoint and the displacement of the L-keypoint on the instance of the symbol (right).</figcaption></figure>
<figure><img src="image_69.png" /><figcaption>Figure 3.9: Few examples of the detected symbols of our system.</figcaption></figure>
# Part II

## Feature indexing in high-dimensional vector space
<figure><img src="image_70.png" /></figure>
# Chapter 4

## State-of-the-art in feature indexing

In the previous chapter, we have provided a local detector of junction points. As it is often the case that the detected keypoints shall be described by local descriptors followed by a further step of descriptor indexing. Hence, the next two chapters are attributed to the problem of feature indexing. This chapter reviews the state-of-the-art in feature indexing in high-dimensional feature vector space. The main ideas, favourable features and shortcomings of each method are carefully discussed. We also provide our subjective remarks for these methods and highlight the need of an advanced contribution for an efficient indexing technique.

### 4.1 Introduction

As we have discussed earlier, robust feature extraction is of central importance for an image processing system. Furthermore, feature indexing is of a crucial need for all the real-time image processing applications. For these reasons, the two next chapters attempt to deal with the problem of feature indexing. At first, an overview of the existing techniques is provided. Next, an attempt is made to give a new contribution in feature indexing for quickly answering the queries of proximity search. Let us first describe the general context of the fast proximity search problem. Considering a scenery where the objects are represented by real feature vectors in a feature vector space $S$, the problem of finding the nearest neighbor of a given query object $q$ over a dataset $X$ has been well-established in the literature. Usually, two problematic factors make the difficulty of this problem. First, the dataset $X$ is often composed of very large data points (e.g., millions of feature vectors). Second, each data point is of a high-dimensional feature space (e.g., $> 100$). A
conventional solution is to sequentially search for each object $p \in X$, to find the closest one of $q$, based on some similarity distance function $d: S \times S \rightarrow \mathbb{R}$. While this brute-force search is practicable in a low-dimensional and small-scale space, it is subjected to the problem of *curse of dimensionality* when working on high-dimensional and large-scale space [Yamamoto et al., 1999].

To deal with this problem, many approaches have been proposed in the literature supporting the tasks of *exact/approximate* nearest neighbor (ENN/ANN) search. These approaches are mainly concerned with designing an efficient indexing scheme. Basically, an indexing algorithm is formulated as the task of reorganizing the data so that it is able to answer quickly the query of proximity search [Beis and Lowe, 1997, Chávez et al., 2001]. For completeness, we introduce hereafter popular formulations of the ENN and ANN search in the literature.

* Exact nearest neighbor search [Gionis et al., 1999]: given a query object $q$ and a database $X$, return an object $p^* \in X$ such that $\forall p \in X: d(p^*, q) \le d(p, q)$.

* Exact K-nearest neighbor search [Gionis et al., 1999]: given a query object $q$ and a database $X$, return at most $K$ objects $\{p_k^*\}_{k=1}^K$ in $X$ such that $\forall p \in X: d(p_1^*, q) \le d(p_2^*, q) \le \dots \le d(p_K^*, q) \le d(p, q)$.

* Range search [Böhm et al., 2001]: given a query object $q$ and a database $X$, and a distance $\epsilon > 0$, return all objects $p^* \in X$ such that $d(p^*, q) \le \epsilon$.

* *c-approximate* nearest neighbor search (*c-ANN*) [Gionis et al., 1999]: given a query object $q$ and a database $X$, and a parameter $c \ge 1$ and let $p_{exact}$ be the closest object to the query $q$, return an object $p^* \in X$ such that $d(p^*, q) \le c \cdot d(p_{exact}, q)$. The parameter $c$ is treated as *approximate factor* or *approximate tolerance*. If $c = 1$, we obtain the exact nearest neighbor of $q$.

Several surveys of indexing algorithms in vector space are presented in [Böhm et al., 2001] and [Liu et al., 2004]. Such methods are often categorized into four classes as shown in Figure 4.1. These approaches are detailed in the following.

<figure><img src="image_71.png" /><figcaption>Figure 4.1: Different approaches for feature indexing in vector space.</figcaption></figure>
<table><caption>Table 4.1: Indexing methods based on space partitioning in $\mathbb{R}^k$ vector space</caption><thead><tr><th>Method</th><th>Worst-case search complexity</th><th>Memory cost</th><th>Dataset & Feature dimension (k)</th><th>Main features</th></tr></thead><tbody><tr><td>Friedman et al., 1977 (KD-tree)</td><td>O(kn<sup>1-1/k</sup>)</td><td>O(n)</td><td>Uniform dataset, k ≤ 20</td><td>Static binary tree, fast tree building, ANN search, good for low-dimensional data space</td></tr><tr><td>Silpa-Anan and Hartley, 2008 (NKD-trees and PKD-trees)</td><td>O(n)</td><td>O(nm)</td><td>SIFT dataset, 5 · 10<sup>5</sup> records, k = 128</td><td>Static binary trees, building m randomized KD-trees, ANN search</td></tr><tr><td>McNames, 2001 (PA-tree)</td><td>O(kn<sup>1-1/k</sup>)</td><td>O(n)</td><td>Uniform dataset, 10<sup>5</sup> records, k ≤ 20</td><td>Static m-ary tree, data alignment with PCA analysis at each partitioning level, high computation cost of lower bounds and tree building, ANN search</td></tr><tr><td>Guttman, 1984 (R-tree)</td><td>O(n)</td><td>O(n)</td><td>Rectangle dataset, 5000 records, k = 2</td><td>Dynamic and balanced m-ary tree, supporting both data points and extended spatial data, difficulties of maintaining minimized overlap and empty space, range and ANN search</td></tr><tr><td>Sellis et al., 1987 (R<sup>+</sup>-tree)</td><td>O(n)</td><td>O(n log n)</td><td>Point dataset, 10<sup>5</sup> records, k = 2</td><td>Dynamic and balanced m-ary tree, supporting both data points and extended spatial data, difficulties of maintaining minimized overlap and empty space, range and ANN search</td></tr><tr><td>Beckmann et al., 1990 (R<sup>*</sup>-tree)</td><td>O(n)</td><td>O(n)</td><td>Rectangle dataset, 10<sup>5</sup> records, k = 2</td><td>Dynamic and balanced m-ary tree, supporting both data points and extended spatial data, difficulties of maintaining minimized overlap and empty space, range and ANN search</td></tr><tr><td>White and Jain, 1996 (SS-tree)</td><td>O(n)</td><td>O(n)</td><td>Uniform dataset, 10<sup>5</sup> records, k ≤ 11</td><td>Dynamic and balanced m-ary tree, supporting both data points and extended spatial data, difficulties of maintaining minimized overlap and empty space, range and ANN search</td></tr><tr><td>Katayama and Satoh, 1997 (SR-tree)</td><td>O(n)</td><td>O(n)</td><td>EigenFace dataset, 10<sup>4</sup> records, k ≤ 100</td><td>Dynamic and balanced m-ary tree, high memory space for storing both the sets of hyper-spheres and the hyper-rectangles, range and ANN search</td></tr><tr><td>Berchtold et al., 1996 (X-tree)</td><td>O(n)</td><td>O(n)</td><td>Uniform dataset and Fourier points, k ≤ 16</td><td>Dynamic and balanced m-ary tree, high memory space for storing both the sets of hyper-spheres and the hyper-rectangles, range and ANN search</td></tr></tbody></table>
### 4.2 Space-partitioning-based methods

The main spirit of the space-partitioning-based indexing methods is that they hierarchically divide the underlying data space into sub-spaces such that the cardinality of every sub-space at the same level is roughly equivalent. The partitioning process is repeated until the cardinality of every new sub-space at the highest level is small enough. As the resulting sub-spaces form a hierarchical representation of the original data space, a tree-based structure is often chosen to represent these sub-spaces. Fast proximity search is then accomplished using various tree traversing strategies including back tracking and priority search. Table 4.1 briefly outlines the main characteristics of these methods. For each method, the main ideas, favourable features, and shortcomings are discussed. The detail of each method is given in the following.

**KD-tree:** KD-tree [Friedman et al., 1977] is probably argued as one of the most popular techniques for feature indexing. Its basic idea is to iteratively perform a partition of X into two roughly equal-sized subsets $X_l$ and $X_r$ by using a hyperplane perpendicular to one of the axes in $R^k$. Particularly, a partition is performed at the $i^{th}$ dimension ($1 \le i \le k$) using a pivot value $s_{med}$ to split X into two subsets $X_l$ and $X_r$. The $X_l$ contains the points whose values at the $i^{th}$ dimension are smaller than $s_{med}$, while the $X_r$ contains the rest of X. A new node is then constructed to record the split information such as the split axis and the pivot value. The pivot value $s_{med}$ is often chosen as the median value at the $i^{th}$ dimension of all the points contained in the underlying dataset. This process is then repeated for the two new subsets $X_l$ and $X_r$ until the size of each subset falls below a pre-defined threshold. The resulting tree is known as a balanced KD-tree whose leaf nodes form a full splitting space of the original data set X (Figure 4.2).

<figure><img src="image_72.png" /><figcaption>Figure 4.2: Illustration of data partitioning of the KD-tree in an 2D space: the resulting KD-tree (left) and the corresponding partitioned space (right).</figcaption></figure>

Searching for a nearest neighbor of a given query point $q$ in the KD-tree is proceeded using a branch-and-bound technique whose pruning rule works as follows: a node $u$ is selected to explore if its hyper-rectangle *does* intersect the hyper-sphere centered at $q$ with a radius equal to the distance of $q$ to the nearest neighbor found so far (Figure 4.3). The most efficient technique for constructing such a KD-tree has been reported by
[Wald and Havran, 2006] with the time complexity $O(N\log(N))$ and the worst case search time was reported as $O(kN^{1-1/k})$ in [Lee and Wong, 1977], where N is the number of nodes in the tree. The KD-tree has been shown to work efficiently when the dimensionality is low (e.g., most often $k \le 20$), otherwise its performance is almost close to brute-force search.

<figure><img src="image_73.png" /><figcaption>Figure 4.3: Illustration of the searching process in the KD-tree: only the regions whose bounding boxes do intersect the circle are inspected.</figcaption></figure>

The KD-tree has been extended to deal with the problem of approximate nearest neighbor search by several researchers. [Beis and Lowe, 1997] introduced an extension of the KD-tree, known as *Best-Bin-First (BBF)* or *priority search*, to quickly find approximate nearest neighbors of a given query point. The improvement of the BBF technique is twofold: (1) it limits the maximum number of leaf nodes to be visited; and (2) it visits the points in the order of increasing distance from the query point to the bins of splitting space. The first improvement indicates that the search procedure will be terminated early once the number of visited nodes in the tree is greater than a specific value. The second improvement guides the search procedure to the branches (or bins) which are "closer" to the query. This is accomplished by defining the distance from the query to a bin as the minimum distance to the points in the bin boundary. Their experimental results showed the efficiency of the BBF algorithm for a synthetic dataset with the dimensionality $k \in \{8, 15\}$.

**NKD-trees and PKD-trees:** The use of priority search was further improved in the work of [Silpa-Anan and Hartley, 2008], where the author constructed multiple KD-trees, called *NKD-trees*, with different settings of the orientation parameter. That is, the data are rotated by different angles before constructing the KD-tree. The obtained results are quite interesting. Based on that, the idea of using multiple KD-trees has been developed in two new different ways. In the first place, multiple *randomized KD-trees* (*RKD-trees*) are constructed, where each tree is built by selecting randomly, at each node, a split axis from few dimensions having the highest variance. In the second place, multiply randomized KD-trees are constructed in the same manner but the data are initially aligned to the principal axes obtained from PCA analysis. The later indexing scheme is thus called principal component KD-trees (*PKD-trees*). Experimental results showed significantly out-
standing performance compared to the use of a single KD-tree because the use of multiple randomized trees gives a better chance of reaching the true answers than using a single tree.

**PA-trees:** The PA-tree (Principal Axis tree) [McNames, 2001] extends the KD-tree at twofold. First, it constructs a bigger fanout tree by partitioning dataset at each step into $n_c$ subsets ($n_c \ge 2$). Second, PCA is applied to the underlying dataset at each partition to select the split axis having the highest variance. Therefore, the obtained regions are treated as hyper-polygons rather than hyper-rectangles as in the KD-tree. Consequently, this approach complicates the process of computing the lower bound of the distance from the query to a node of the tree. To facilitate this matter, the author proposed a convenient solution for computing the lower bound using the cosine rule as illustrated in Figure 4.4. Let $x$ be any point in some node of the tree, say the Region 2, and $b_2$ be the projection of a query point $q$ on the common split hyper-plane between the Region 2 and the Region 1. By applying the cosine rule, the distance $d(q, x)$ is computed as follows:

$$d(q, x)^2 = d(q, b_2)^2 + d(b_2, x)^2 - 2d(q, b_2)d(b_2, x) \cos \angle qb_2x \quad (4.1)$$

As the angle $\angle qb_2x \ge 90^\circ$, Equation 4.1 implies: $d(q, x)^2 \ge d(q, b_2)^2 + d(b_2, x)^2$. Hence, the lower bound $d(q, x)^2$ can be recursively computed via $d(b_2, x)^2$ as follows:

$$d(q, x)^2 \leftarrow d(q, b_2)^2 + d(b_2, x)^2 \quad (4.2)$$

It is noted that each time of computing an intermediate distance, it requires $O(k)$ real computations to compute $d(q, b_2)^2$ as the PCA is applied to every level of partitioning. Despite of the heavy computations on the lower bounds, the PA-tree was showed to outperform many other indexing schemes for many datasets.

<figure><img src="image_74.png" /><figcaption>Figure 4.4: Illustration of computing the lower bound using the cosine rule of the PA-tree (reprinted from [McNames, 2001]).</figcaption></figure>
**R-tree:** The R-tree (Rectangle tree) [Guttman, 1984] is *dynamically* built based on an idea that nearby points are grouped together to form an internal node represented by a minimum bounding hyper-rectangle (Figure 4.5). Each internal node of the R-tree has at most *M* children and the root contains at least two children. Each leaf node contains also at most *M* entries corresponding to data points. Inserting a new data object, *p*<sub>n</sub>, into the R-tree is performed by two processes: target node searching and node splitting. The former process traverses down the tree and considers all the subtrees whose hyper-rectangles contain *p*<sub>n</sub>. Tie<sup>1</sup> resolving is treated by choosing the child node whose bounding hyper-rectangle needs least area enlargement to include *p*<sub>n</sub>, and further selecting the one with the smallest volume if necessary. This continues when reaching a leaf node. If the number of entries at this leaf node does not exceed *M*, the new point *p*<sub>n</sub> is simply inserted to the leaf node. Otherwise (e.g., this case is called *overflow*), node splitting is performed to split (M + 1) points into two new nodes, subjected to an optimal criterion that the total volume of two new bounding hyper-rectangles is minimized. Then, upward change propagation is performed to update the bounding hyper-rectangles, and to split higher-level nodes, if necessary. Once the R-tree has been constructed, the *range* search algorithm starts from the root and visits every subtree whose the bounding hyper-rectangle does intersect the range query. This continues when reaching the leaf nodes, where only the entries lying inside the query range are returned as final answers. In this way, many branches of the tree are pruned, providing rapid answers of proximity queries. The main challenge of this tree is realized on the process of constructing a balanced tree, while maintaining at the same time the minimal overlap among the bounding hyper-rectangles.

<figure><img src="image_75.png" /><figcaption>Figure 4.5: The construction of the R-tree (right) for a dataset (left); (reprinted from [Guttman, 1984]).</figcaption></figure>

**R<sup>+</sup>-tree:** The R<sup>+</sup>-tree [Sellis et al., 1987] improves the R-tree in a different way. It optimizes the overlapping of internal nodes by allowing a new data object *p*<sub>n</sub> to be inserted into several leaf nodes. To insert *p*<sub>n</sub> into the tree, the target leaf node is first detected in a similar manner as in the R-tree. If overflow happens at the leaf node, node splitting will be performed. In this case, the leaf node is partitioned into two subsets by finding a

<footnote n="1">¹This means that more than a single solution is happened.</footnote>
cutting hyper-plane perpendicular to one of the axes, subjected to several optimal criteria, including minimal coverage of dead space (i.e., empty space of each new subset) and minimal number of hyper-rectangle splits. In contrast to the splitting process of the R-tree, node splitting in $R^+$-tree may be propagated to both upward and downward nodes. Range searching in the $R^+$-tree is more efficient than in the R-tree as the overlap space is taken place at the leaf nodes only. However, $R^+$-tree requires a logarithmically increase of memory space as the points are duplicated at the leaf nodes. The overlap is therefore shifted to the leaves, making the inspection of multiple paths at the leaf nodes necessary.

**$R^*$-tree:** The $R^*$-tree [Beckmann et al., 1990] is another development of the R-tree, where the main improvements are attributed to two process of searching for a target leaf node and inserting a new data object. Instead of using the sole criterion of area enlargement to determine the best path, the $R^*$-tree makes use of both area enlargement and overlap enlargement when searching for the best path. Particularly, given a new data object $p_n$ to be inserted, the process of searching for the target leaf node is basically relying on that in the R-tree, but with the following modifications at each step:

* If the children of the current node, $p_{curr}$, point to the leaves, then choose the entry whose hyper-rectangle needs least overlap enlargement to include $p_n$. Tie resolving is handled by selecting the one, which needs least area enlargement, followed by the selection of the entry with the smallest hyper-rectangle volume, if needed.

* If the children of $p_{curr}$ do not point to the leaves, then choose the entry whose hyper-rectangle needs least area enlargement to include $p_n$, followed by selecting the one having the smallest hyper-rectangle volume.

Then, $p_n$ is added as a new entry of the obtained target node if its capacity is not full. If it is not the case, overflow is handled by either reinsertion or node splitting. Reinsertion is applied if the children of the target node have not been reinserted previously, and is accomplished by reinserting the $k$ entries farthest away from the current node. If node splitting is applied, several optimal criteria are employed, including minimization of total volume of two new nodes, minimization of overlap volume of two new nodes, minimization of total surface of two new nodes, minimization of storage utilization. Upward propagation for handling overflow might be needed. Experiments showed a superiority search performance, compared to the original R-tree. However, as mentioned in [Berchtold et al., 1996], the $R^*$-tree performs poorly on high-dimensional space (> 5) because the overlap at the internal nodes increases rapidly with respect to the increase of data dimensionality.

**SS-tree:** The SS-tree (similarity search tree) [White and Jain, 1996] is a new variation of the $R^*$-tree where its nodes are represented by hyper-spheres rather than hyper-rectangles. A hyper-sphere in the SS-tree is defined by a centroid and a radius computed as the biggest distance from the centroid to all elements contained in the hyper-sphere. To insert a new object $p_n$ into the tree, a target leaf node is first detected to contain $p_n$ by descending the tree and choosing at each step the subtree whose centroid is closest to $p_n$. Next, the process of node insertion and reinsertion is quite similar to that in the $R^*$-tree. The only difference is in the process of node splitting, where the dimension with the highest variance is selected as the split axis, and the split location is selected such that it minimizes the sum of variances on each side of the split plane. Exper-
iments shows better results of the SS-tree in comparison with the *R*\*-tree. However, the problem with the SS-tree is the difficulty in yielding overlap-free of the splitting process [Böhm et al., 2001]. In addition, the use of the hyper-spheres, in general, occupies higher space than the hyper-rectangles in high-dimensional space, making the similarity search less efficient [Katayama and Satoh, 1997, Böhm et al., 2001].

<figure><img src="image_76.png" /><figcaption>Figure 4.6: The construction of the SR-tree (reprinted from [Katayama and Satoh, 1997]).</figcaption></figure>

**SR-tree:** The SR-tree (Sphere/Rectangle-tree) [Katayama and Satoh, 1997] overcomes the weakness of the SS-tree by assembling the spirit of the *R*\*-tree and SS-tree into an unified scheme. A node of the SR-tree is represented by the common space (R) of the hyper-sphere and the hyper-rectangle as shown in Figure 4.6. Unfortunately, the space *R* is not explicitly computed due to the complicated computation of the intersection between the hyper-rectangle and hyper-sphere. Instead, each node of the SR-tree records the information of both the hyper-sphere and hyper-rectangle. Given a data object *q*, an estimation of the distance from *q* to *R* is provided as the biggest value of the minimum distances from *q* to the corresponding hyper-sphere and hyper-rectangle. Next, the insertion algorithm is performed essentially similar to the SS-tree. Experimental results reported better performance of the SR-tree, compared to the SS-tree and *R*\*-tree.

**X-tree:** One common weakness of the R-tree, *R*\+-tree, and *R*\*-tree is to maintain a minimized overlap volume of the bounding hyper-rectangles when proceeding on high-dimensional space. All the heuristic solutions introduced in these works address this problem to some extent but not completely resolved [Berchtold et al., 1996, Böhm et al., 2001]. The X-tree (eXtended node tree) [Berchtold et al., 1996] gives more investigation for the optimization of overlap in the nodes of R-tree-based structures. The new investigation of the X-tree is two-fold. First, it introduces a new kind of internal node, so-called super-node. A super-node is similar to an internal node used in the R-tree-based structures except with a big capacity for containing its entries. Second, it introduces a new split procedure for optimizing the overlap of the hyper-rectangles. The split procedure first tries to find an optimal split of the overflow node by using the same heuristic rules as presented in the R\+-tree and R\*-tree. If the obtained overlap is still high enough, the split procedure tries to find an overlap-free split relying on split history obtained previously. If the obtained split results in unbalance nodes (i.e. the difference in cardinality of each new node is high enough), the
split procedure terminates without any available splits. In this case, the current node is extended to a super-node. The super-nodes again can be extended by one addition block if no available splits are found. Experiments showed the efficiency of the X-tree compared to R\*-tree and TV-tree by up to two orders of magnitude in higher-dimensional space.

## 4.3 Clustering-based methods

The clustering-based indexing methods differ from the space-partitioning-based methods mainly in the step of tree construction. Instead of dividing the data using a hyper-plane, these methods employ a clustering method such as K-means and K-medoids to iteratively partition the underlying data into sub-clusters. The partitioning process is repeated until the size of every sub-cluster falls below a threshold. A tree-based structure is constructed to hierarchically represent the resulting sub-clusters at all levels of decomposition. Proximity search is often handled using a branch-and-bound algorithm. Table 4.2 briefly outlines the main characteristics of these methods.

**K-means clustering tree:** One of the first clustering-based trees was reported by [Fukunaga and Narendra, 1975]. The proposed algorithm recursively divides all points in the dataset into smaller regions using the K-means clustering technique, and constructs a corresponding clustering tree. Each node *p* of the tree has the following parameters: {$S_p, M_p, N_p, r_p$} corresponding to the set of data points contained in *p*, the cluster center, the number of data points, and the farthest distance from $M_p$ to an $X_i \in S_p$, respectively. The iterative clustering process terminates when the size of each obtained region falls below a threshold. Searching for k-nearest neighbors of a given query *q* is then proceeded by a branch-and-bound algorithm. Let *Y* be the current nearest neighbor of *q*, two following pruning rules are used to eliminate the branches too far from the query:

* Rule 1: A node *p* will be not searched if $d(q, Y) + r_p < d(q, M_p)$ as illustrated in Figure 4.7 (a).

* Rule 2: A point $X_i \in S_p$ will not be the nearest neighbor of *q* if $d(q, Y) + d(X_i, M_p) < d(q, M_p)$ as illustrated in Figure 4.7 (b).

Experiment results demonstrate the efficiency of the proposed algorithm for a small dataset (1000 data points).

<figure><img src="image_77.png" /><figcaption>Figure 4.7: Illustration of the rule 1 (a) and rule 2 (b) for tree pruning; (reproduced from [Fukunaga and Narendra, 1975])</figcaption></figure>
<table><caption>Table 4.2: Indexing methods based on clustering in $\mathbb{R}^k$ vector space.</caption><thead><tr><th>Method</th><th>Worst-case search complexity</th><th>Memory cost</th><th>Dataset & Feature dimension (k)</th><th>Main features</th></tr></thead><tbody><tr><td>Fukunaga and Narendra, 1975</td><td>O(n)</td><td>O(n)</td><td>Uniform distribution data, three thousand records, k = 8</td><td>Single K-means-based clustering tree, high computation cost of tree building, ANN search</td></tr><tr><td>Nister and Stewenius, 2006</td><td>O(n)</td><td>O(n)</td><td>SIFT features, one million images, k = 128</td><td>Single hierarchical vocabulary tree based on K-means, high computation cost of tree construction, ANN search</td></tr><tr><td>Leibe et al., 2006</td><td>O(n)</td><td>O(n)</td><td>SIFT features, one million records, k = 128</td><td>Combined partitional-agglomerative clustering algorithm, fast matching with the ball tree, fast clustering and tree building, range and ANN search</td></tr><tr><td>Muja and Lowe, 2009</td><td>O(n)</td><td>O(n)</td><td>SIFT features, 31 billions records, k = 128</td><td>Single K-means-based clustering tree, priority search, ANN search</td></tr><tr><td>Muja and Lowe, 2012</td><td>O(n)</td><td>O(n + m)</td><td>SIFT and SURF features, k ∈ {128, 64}</td><td>Multiple randomized K-medoids-based clustering trees, random selection of centers, modest computation cost of tree building, priority search, ANN search</td></tr></tbody></table>
**Vocabulary K-means tree:** [Nister and Stewenius, 2006] proposed a hierarchical vocabulary tree for representation of feature vectors. At first, the K-means algorithm is employed to partition the feature vectors into K groups, where the parameter K is treated as the branching factor of the tree rather than the number of cluster centers as usual. Each group contains the feature vectors closest to a particular center. This process is recursively repeated for each of the obtained groups until the height of the tree is higher than a pre-defined level L. The nodes at the highest level are considered as leaf nodes, containing the actual feature vectors. In this way, the tree is hierarchically constructed and defines the quantized cells of feature vectors treated as the visual vocabularies. Figure 4.8 illustrates this process. Next, each node i of the tree is assigned with a weight $w_i$ computed as follows:

$$w_i = \ln \frac{N}{N_i} \qquad (4.3)$$

where N is the number of images in the database and $N_i$ is the number of images in the database such that at least one feature vector contained in each image passes through node i. Given a query feature vector, a path is defined by traversing down the tree until a leaf node is reached. At each intermediate level, the node with the center closest to the feature vector is selected to be further explored. Image retrieval is performed by computing a relevant score between a query image and each database image. This relevant score is computed as follows. A global descriptor is created for each image by accumulating the weighted frequencies of the feature vectors of the image whose corresponding paths pass through the nodes of the tree. The relevance score of two images is then defined as the normalized difference between the global descriptors.

<figure><img src="image_78.png" /><figcaption>Figure 4.8: Construction of the vocabulary tree with the branching factor K = 3 (reprinted from [Nister and Stewenius, 2006]).</figcaption></figure>

**Agglomerative clustering tree:** [Leibe et al., 2006] proposed an efficient clustering and fast matching methods by combining the advantages of the K-means and agglomerative techniques. In contrast to the traditional K-means technique, the agglomerative algorithm builds a hierarchical clustering tree in a bottom-up manner. Particularly, starting from the data points, it iteratively selects and merges the pairs of two clusters up to the root node
of the tree. At the beginning, each data point is a cluster. The criterion for selecting two clusters to be merged is based on the average similarity which is computed as the average pairwise similarity between the points of the two clusters [Guha et al., 1998]. Therefore, the proposed algorithm can be considered as a fusion of both top-down and bottom-up clustering strategies. First, an initial partitioning of the data points is performed by dividing them into two subsets corresponding to Laplacian maxima and minima. Each subset is then further partitioned using the K-means algorithm where K is the number of clusters in each subset and is set to be small for gaining efficiency. Next, the agglomerative clustering technique is independently applied to each resulting partition. Finally, the agglomerative algorithm is run once more on all the clusters obtained previously. For fast matching and searching, a binary ball tree is constructed in which two children of a node correspond to two clusters merged during the agglomerative clustering process. Fast searching is accomplished with the use of the triangle inequality rule. Experimental results showed the efficiency of the proposed clustering algorithm and significant speedup is achieved for feature matching and searching.

**Priority K-means clustering tree:** [Muja and Lowe, 2009] extended the work of [Fukunaga and Narendra, 1975] by incorporating the use of priority queue or best-first [Beis and Lowe, 1997] to the hierarchical clustering tree. In their work, approximate nearest search is proceeded by traversing down the tree and always choosing the node whose cluster center is closest to the query point. Each time when a node is selected for further exploration, the other sibling nodes are inserted into a priority queue, which contains a sequence of nodes stored in the increasing order of the distances to the query. This continues when reaching a leaf node, followed by a sequence search for the points contained in this node. Backtracking is then invoked starting from the top node stored in the priority queue. During the search process, the algorithm maintains adding new candidate nodes to the priority queue along the search path. The search process terminates early when a given number of leaf nodes have been visited. The authors reported in their experiments that the use of priority search along with the hierarchical clustering tree gives significant improvements, compared to the locality sensitive hashing algorithm and the KD-tree for many datasets.

**Multiple K-medoids clustering trees:** The hierarchical clustering tree has been further extended by [Muja and Lowe, 2012] to build up multiple hierarchical clustering trees. Single tree construction is essentially similar to that in [Muja and Lowe, 2009] except the use of a new K-medoids clustering algorithm. Particularly, the cluster centers are randomly selected from the input points. The clusters are then built by assigning each input point to the closest center. To achieve fast construction of the tree, the step of cluster center optimization by using squared error minimization is discarded. If the size of a cluster is sufficiently small, a leaf node is constructed containing all points in this cluster. Otherwise, an internal node is constructed associated to the cluster center. These steps are recursively applied to the new clusters corresponding to the internal nodes. Following the advanced work of [Silpa-Anan and Hartley, 2008] concerned with the advantage of using multiple randomized trees for nearest neighbor search, multiple hierarchical clustering trees are thus constructed in this work. Approximate nearest search is proceeded in parallel from multiple hierarchical clustering trees. Single search for each tree is accomplished in a same manner to that in [Muja and Lowe, 2009]. Once the single search is finished for
every tree, backtracking is applied starting from the top node stored in the priority queue. Experimental results showed that significant improvements are achieved by using parallel search in multiple hierarchical clustering trees.

## 4.4 Hashing-based methods

Apart from the space-partitioning-based and clustering-based approaches, there is also a wide range of hashing-based techniques dealing with the problem of ANN search. Such methods are often applied to binary features employing the Hamming distance as the similarity score of two feature vectors. The basic idea is to hash data points into buckets in which similar points might likely be hashed in the same or adjacent buckets, and dissimilar points are likely to be hashed into different buckets. Multiple hash tables are often used to ensure a good chance of collision among the similar points. ANN search can be efficiently accomplished with a sub-linear time complexity. Table 4.2 summarizes the representative methods in this domain, and the details are discussed in the following paragraphs.

**Locality Sensitive Hashing (LHS):** The LSH-based indexing scheme has been known as one of the most popular methods, which can perform proximity search with a sub-linear time for high-dimensional data [Gionis et al., 1999, Indyk and Motwani, 1998]. The key idea of LHS is to design the hash functions that the similar points may be hashed with a high probability of collision, while the dissimilar points may likely be hashed with different keys. Particularly, let $D(p, q) \in [0, 1]$ be a similarity function of two given points $p$ and $q$ in a dataset $X$, a LSH function family $H = \{h : X \to U^1\}$ is defined as follows:

$$ Pr_{h \in H}[h(p) = h(q)] = D(p, q) \qquad (4.4) $$

To design such a kind of hash functions, the authors proposed to conceptually transform the $d$-dimensional feature vectors into $d'$-dimensional unary representations, where the $L_1$ distance is preserved ($d' > d$). Next, the hash functions are constructed by selecting $l$ subsets, $\{g_i\}_{i=1}^l$. Each is composed of $k$ elements uniformly sampled with replacement from the set $R = \{1, 2, \dots, d'\}$ (i.e., the axes in the unary space). Each subset $g_i$ can be regarded as a hash function composing of $k$ random lines:

$$ g_i : X \to U^k \qquad (4.5) $$

Equivalently, a hash function $g_i$ is composed of $k$ LSH functions:

$$ g_i(p) = \{h_{i1}(p), h_{i2}(p), \dots, h_{ik}(p)\} \qquad (4.6) $$

where $h_{it} \in H$. As there are $l$ hash functions, $l$ hash tables are created to store all the projected feature vectors. More precisely, given any data point $x \in X'$ in $d'$-dimensional space, the hash table $T_i$ ($1 \le i \le l$) is constructed as follows:

$$ T_i = g_i(x) = \{h_{i1}(x), h_{i2}(x), \dots, h_{ik}(x)\} \qquad (4.7) $$
<table><caption>Table 4.3: Hashing-based indexing methods in $\mathbb{R}^k$ vector space.</caption><thead><tr><td>Method</td><td>Worst-case search complexity</td><td>Memory cost</td><td>Dataset &amp; Feature dimension (k)</td><td>Main features</td></tr></thead><tbody><tr><td>Indyk and Motwani, 1998</td><td>O(n<sup>1/c</sup>)</td><td>O(nl)</td><td>Histograms of color, 20000 records, k = 64</td><td>Locality Sensitive Hashing (LHS), randomized projection, using l hash tables, (c, r)-ANN search, expensive space cost</td></tr><tr><td>Kulis and Grauman, 2009</td><td>O(n<sup>1/c</sup>)</td><td>O(nl)</td><td>80 million tiny images, GIST features, k = 384</td><td>Kernelized LSH, good for kernelized data, (c, r)-ANN search, expensive space cost</td></tr><tr><td>Panigrahy, 2006</td><td>O(n<sup>2.09/c</sup>)</td><td>O(n)</td><td>None</td><td>Entropy LSH, less cost of memory space, (c, r)-ANN search</td></tr><tr><td>Lv et al., 2007</td><td>O(n<sup>1/c</sup>)</td><td>O(nl/m)</td><td>Audio dataset, 2.7 million records, k = 192</td><td>Multi-probe LSH, search for multiple adjacent buckets, less cost of memory space (m ∈ [14, 18]), (c, r)-ANN search</td></tr><tr><td>Auclair, 2009</td><td>O(n)</td><td>O(nl)</td><td>512000 images, 256 SIFT descriptors/image, K = 128</td><td>Making use of the highest distinctive dimensions as projection axes, poor collision of similar data points, ANN search</td></tr><tr><td>Jain and Doermann, 2012</td><td>O(n)</td><td>O(nl)</td><td>Logo dataset, SURF features, 664 million records, k = 32</td><td>Making use of the highest and lowest distinctive dimensions as projection axes, poor collision of similar data points, ANN search</td></tr></tbody></table>
For ANN search, a query $q$ is first projected using the designed hash functions resulting in $l$ indices $\{g_1(q), g_2(q), \dots, g_l(q)\}$. Each is a $k$-dimensional vector. Using these indices, the search algorithm looks up in the hash tables to obtain a set of candidate points. Next, sequence search is applied to return the $K$-nearest neighbor points. This method has been shown to provide a sub-linear complexity of computation, $O(n^{1/c})$, for answering the $(c, r)$-ANN query. However, the main drawback is the utilization of a huge amount of memory, $O(nl)$, to store $l$ hash tables, raising another problem known as bottleneck as argued in [Auclair, 2009].

**Kernelized LSH:** [Kulis and Grauman, 2009] extended the LSH to the case, when the similarity function is an arbitrary kernel function $\kappa$: $D(p, q) = \kappa(p, q) = \phi(p)^T \phi(q)$. Given an input feature vector $x$, the problem is then to design a specific LSH function over the feature vector $\phi(x)$ rather than $x$ itself, where $\phi(x)$ is some unknown embedding function. In their work, the LSH hash function is constructed as: $h(\phi(x)) = \text{sign}(r^T \phi(x))$, where $r$ is a random hyperplane drawn from $N(0, I)$ and is computed as a weighted sum of a subset of the database feature vectors. Since the newly derived $h(\phi(x))$ satisfies the LSH property (i.e., $\text{Pr}_{h \in H}[h(p) = h(q)] = D(p, q)$), the new indexing scheme is thus capable of performing similarity searching in a sub-linear time complexity, while being useful to the cases of kernelized data. However, the problem of expensive cost of memory space is still unsolved.

**Entropy-based LSH:** To reduce the hash tables used in the original LSH indexing scheme, [Panigrahy, 2006] introduced an entropy-based LSH indexing technique. The basic idea of the entropy-based LSH method is as follows: given a query $q$ and a parameter $r$ of the distance from $q$ to its nearest neighbor, the synthetic neighbors of $q$ within a distance $r$ are randomly generated to be hashed using the LSH functions. The obtained hash keys are then used to access the bucket candidates in which we expect that the true nearest neighbor of $q$ may be present. The rest is then shifted to determine a sufficient large number of the synthetic neighbors of $q$ to ensure that the c-approximate nearest neighbor search (c-ANN) can be resolved under a sub-linear time. For this purpose, the author proposed to compute the entropy $I(h(p)|q, h)$, given a query $q$ and the LSH function $h$. This entropy is then used to estimate the required search time of $O(n^{(1+o(1))1.47/c})$. In other words, it is required to generate $n^{(1+o(1))1.47/c}$ synthetic nearest neighbors of $q$, each of which is used to access the corresponding bucket. Sequence search is applied to the data points contained in the obtained buckets. A detailed analysis of the entropy-based LSH algorithm was presented in [Lv et al., 2007] showing that this method is subjected to several weaknesses. First, it is not trivial to determine a good setting for the parameter $r$ since this value is highly dependent on the dataset. Second, without prior knowledge about the data, the estimation of $r$ is not reliable. Third, the process of generating the synthetic neighbors of $q$ by sampling the feature space from a normal distribution is inefficient and sensitive to the parameter quantization. Finally, experiment results performed by [Lv et al., 2007] show that the entropy-based LSH algorithm does not give noticeable search improvement, compared to the original LSH method.

**Multi-probe LSH:** [Lv et al., 2007] proposed another approach, known as *multi-probe* LSH, to reduce the hash tables used in the original LSH indexing scheme. This method requires less memory space, while retaining the same search precision as in the original LSH algorithm. The basic idea is to search multiple buckets, which probably contain the
<figure><img src="image_79.png" /><figcaption>Figure 4.9: Illustration of searching process in the multi-probe LSH algorithm: $g_i(q)$ is the hash values of the query $q$ (green buckets), and $g_i(q) + \Delta_i$ are the adjacent buckets for probing (reprinted from [Lv et al., 2007]).</figcaption></figure>

potential nearest neighbors of the query. The rationale of looking at multiple buckets is realized based on the fact that if a data point $p$ is close to another data point $q$ but they are not hashed into the same bucket, then there is a high chance that they are hashed into two "close" buckets. Therefore, the multi-probe LSH method works by constructing a sequence of hash perturbation vectors $\{\Delta_1, \Delta_2, \dots, \Delta_l\}$. Next, let $\{g_1(q), g_2(q), \dots, g_l(q)\}$ be the indices of a given query $q$ obtained from the original LSH algorithm, the new indices $\{g_i(q) + \Delta_i\}$ are used to look up in the hash tables for obtaining the nearest neighbor candidates (Figure 4.9). Finally, sequence search is applied to return the $K$-nearest neighbor points. In this way, the proposed method reduces the space requirement and increases the chance of finding the true answers. Experimental results showed a significant improvement in space efficiency of the multi-probe LSH method, compared with the original LSH scheme, when working on a wide range of search precision. However, our experiments, using the implementation provided from the open source FLANN library<fnref n="2" />, showed that the multi-probe LSH algorithm performs poorly, compared to the randomized KD-trees and hierarchical K-means tree, applied to the SIFT and GIST features.

**Distinctive-dimension-based hashing:** [Auclair, 2009] proposed a new indexing algorithm for near duplicate images based on local descriptors (typically SIFT descriptor). Once the feature vectors are computed, the value in each dimension of this vector is normalized by using the mean and standard deviation in the same dimension, computed from all the feature vectors. These values are then sorted to select the $K$ strongest scores whose corresponding dimensions are considered as the distinctive positions. The hashing function is then defined by selecting the first $K$ distinctive dimensions of a given feature vector. Given a query feature vector, its hash keys are first computed. Next, the hash keys are used to access the appropriate buckets in the hash tables to obtain a shortlist of nearest candidates. Sequence search is finally applied to return the approximate nearest answers.

<footnote n="2">²http://www.cs.ubc.ca/~mariusm/index.php/FLANN/FLANN</footnote>
For this approach, it can be noted that the use of distinctive dimensions is too sensitive to discriminate the objects since very similar points could be hashed with different keys (Figure 4.10). In addition, the task of selecting a good setting for $K$ needs to be more thoroughly studied.

<figure><img src="image_80.png" /><figcaption>Figure 4.10: Example of the hashing algorithm in [Auclair, 2009]: Two points $x^1$ and $x^2$ are quite similar in a 8-dimensional space; $D(x^1)$ and $D(x^2)$ are the sorted vectors with respect to the distinctive dimensions; $Key(x^1)$ and $Key(x^2)$ are the obtained hash keys where $K = 2$ in this example; (reprinted from [Auclair, 2009]).</figcaption></figure>

[Jain and Doermann, 2012] improved the work of [Auclair, 2009] in a different way. In their work, the hash functions are defined by selecting 6 dimensions with the highest distinctive scores and 6 dimensions with the lowest distinctive scores. In the retrieval stage, the hash keys of a given descriptor are used to find out all the keypoints having collision with these hash keys to form the candidate list. All these candidate keypoints are finally filtered using a geometric consistency process. Experimental results showed quite good performance under the context of logo retrieval. The main issue of this method remains the same as in the [Auclair, 2009] which concerns the sensitiveness of the distinctive dimensions.

## 4.5 Other methods

All the methods discuss so far can be regarded as direct-based indexing algorithms since they directly address the problem of feature indexing by performing an offline phase of reorganization of the data. In addition to these direct-based indexing method, indirect-based techniques for feature indexing have also been studied in the literature. Such techniques often address the problem of fast proximity search by improving the feature-extraction step [Mikolajczyk and Matas, 2007, Roy et al., 2012] and/or the matching step [Cheng et al., 1984, Mori et al., 2001]. Table 4.4 reviews some typical methods in this domain, and the details are discussed in the following.

[Mikolajczyk and Matas, 2007] employed linear projection to improve the descriptors
<table><caption>Table 4.4: Other indexing methods in $\mathbb{R}^k$ vector space.</caption><thead><tr><th>Method</th><th>New descriptors</th><th>New matching</th><th>Search performance</th><th>Main features</th></tr></thead><tbody><tr><td>Mikolajczyk and Matas, 2007</td><td>×</td><td></td><td>+++</td><td>Feature improvement by de-correlation and dimension reduction, limited to Gaussian distribution of feature variation</td></tr><tr><td>Roy et al., 2012</td><td>×</td><td>×</td><td>++</td><td>Coarse-to-fine matching, dedicated to historical documents, sensitive to glyphs segmentation</td></tr><tr><td>De-Yuan Cheng, 1984</td><td></td><td>×</td><td>+++</td><td>Partial distance search (PDS), loop unrolling</td></tr><tr><td>Mori et al., 2001</td><td>×</td><td>×</td><td>+++</td><td>The use of few representative shape contexts to filter the matches, feature improvement with *shapeme*, and vector quantization of shape contexts</td></tr></tbody></table>

for fast tree matching. In their work, the descriptor space are transformed into a new space by applying a whitening linear transformation (i.e., for feature de-correlation), followed by a PCA transformation (i.e., for dimension reduction). The experiments, applied to the SIFT descriptors, demonstrated that the nearest searching procedure based on the transformed descriptors provides high accuracy and low running time on various tree-based indexing structures. It is noticed that even though the proposed linear transformation can be applied for any descriptors, it is restricted to one constraint that the variation of two descriptors of the same scene or surface has a Gaussian distribution.

[Roy et al., 2012] presented a coarse-to-fine indexing technique for fast text retrieval in historical documents. In their work, the input text lines are indexed at two levels as illustrated in Figure 4.11:

* Coarse level: the text lines are segmented and described with respect to the connected components (CCs).

* Finer level: the text lines are segmented and described with respect to the glyphs.

For each level of features, a training step is applied to generate a codebook. For each codebook, the representative elements are learnt using an unsupervised algorithm. After that, the text lines of the books are indexed at two levels corresponding to two codebooks. Given a query text line, a retrieval step is first applied to the CC-based codebook and then to the glyph-based codebook, if necessary. The approximate string matching algorithm is finally applied to retrieve the answers. The authors demonstrated a good performance (i.e., high accuracy and low cost of processing time) of the proposed system for the historical books. The advantage of this system is the avoidance of the word segmentation step which is very sensitive to noise. However, the robustness of the primitive segmentation step (i.e.,
CCs and glyphs) needs to be further studied as this step is sensitive to the noise and the variation of characters. For example, the primitive segment step is sensitive to the text characters having multiple components, such as 'i', 'à', 'é'.

<figure><img src="image_81.png" /><figcaption>Figure 4.11: Illustration of multi-level indexing in [Roy et al., 2012].</figcaption></figure>

[Cheng et al., 1984] presented a simple and efficient improvement of the matching step for brute-force search, so-called partial distance search (PDS). Given two feature vectors $u$ and $v$ in a $k$-dimensional space and a reference value $D_{ref}$ computed as the distance from the query point to the nearest neighbor found so far, let $sum_i(u, v)$ be the partial sum of square differences up to the $i^{th}$ coordinate ($1 \le i \le k$):

$$sum_i(u, v) = \sum_{t=1}^{i} (u_t - v_t)^2 \qquad (4.8)$$

The basic idea of the PDS technique is that the calculation of the distance between $u$ and $v$ is terminated early at the dimension $i^{th}$ ($1 \le i \le k$) if $sum_i(u, v) \ge D_{ref}$. This partial distance search has been widely used in the literature [McNames, 2001, Muja and Lowe, 2009, Muja and Lowe, 2012] as an alternative solution for brute-force search.

[Mori et al., 2001] introduced two new matching methods for efficient retrieval of similar shapes based on the shape context descriptor. The first method is so-called representative shape contexts, which concerns the use of few shape contexts of a given query shape to match against the full set of shape contexts of every database image. The rationale of this idea is that when trying to match two sufficiently different shapes, few or none of the shape contexts from the first shape have good correspondences on the other. Therefore, if there is a lack of good correspondences between two shapes, we can assume that the two shapes are different. In this method, only five representative shape contexts are randomly computed for the query. Then, a matching step is performed using these representative shape contexts to obtain a shortlist of candidate shapes.

The second matching method is so-called shapemes, which concerns the use of vector quantization to construct a discrete set of the shape context labels. Its basic idea is to quantize the shape context features into different classes using an unsupervised algorithm. To be more specific, in the training stage, a clustering technique such as K-means, can be used to cluster a set of shape context feature vectors into discrete groups. Each group is then quantized by an integer number called shape context label. Once all the context labels are derived, each shape is then represented by a global descriptor, as illustrated in Figure 4.12, constructed by counting the frequency of each label from the shape contexts
of the shape. Both these shape context-based matching methods have been verified on different shape datasets, demonstrating quite good results.

<figure><img src="image_82.png" /><figcaption>Figure 4.12: Shapeme descriptor: (a) an input image; (b) shapeme representation of (a); (reproduced from [Mori et al., 2001]).</figcaption></figure>

4.6 Discussion

Table 4.5 presents some key conclusions about all the indexing methods discussed pre-
viously. The hashing-based approaches, particularly the LHS-based indexing methods,
provided a solid basis dealing with the problem of ANN search. These methods have
been theoretically shown to perform approximate nearest neighbor search in a sub-linear
time, given a sufficient space requirement. The main drawback is the high memory cost
to construct the hash tables. Several works [Panigrahy, 2006, Lv et al., 2007] have been
investigated to deal with this weakness but the obtained results are not satisfactory. In
addition, the search precision could be a problem because the true nearest neighbors could
be hashed into many adjacent buckets, making the access to a single hash bin insufficient
to recover the true answers. Furthermore, these hashing-based methods are known to
work well on binary features, their behaviors on real feature vectors would need further
investigation.

The clustering-based approaches have shown quite good performance in a wide range
of feature types, and scaled well to the increase of the dataset size. Extensive experiments
have been performed in [Muja and Lowe, 2009, 2012] showing that these approaches are
able to give a speedup of hundreds of times, compared to the brute-force search, while
maintaining a reasonably high search precision (i.e., > 80%). They are also known to
work well on both binary features and real feature vectors. The main disadvantage is the
expensive time of the process of tree construction. However, this matter would make sense
only to several application domains where a regular update of data and tree construction
is needed.

The space-partitioning-based approaches, particularly the KD-tree-based algorithms,
seem to be the most appropriate solutions for all aspects of search precision, search
speedup, and tree construction time. Extensive experiments have been performed by
[McNames, 2001, Silpa-Anan and Hartley, 2008, Muja and Lowe, 2009], demonstrating the
same search performance as the clustering-based approaches. These experiments also reveal
that the space-partitioning- and clustering-based approaches significantly outperform the
<table><caption>Table 4.5: Summary of different approaches for feature indexing in vector space</caption><thead><tr><td>Approach</td><td>Pros.</td><td>Cons.</td><td>Search performance</td></tr></thead><tbody><tr><td>Space partitioning</td><td>Fast tree building, applicable to both real and binary features, good generality</td><td>Poor ENN search, difficulties of dynamic update of the tree, tree-traversing-inefficient</td><td>+++</td></tr><tr><td>Clustering</td><td>Applicable to both real and binary features, good generality, memory space-efficient</td><td>Offline computation complexity, poor ENN search, difficulties in dynamic update of the tree</td><td>+++</td></tr><tr><td>Hashing</td><td>Applicable to binary features only, easy to index new data points</td><td>Memory space complexity, poor ENN search</td><td>+</td></tr><tr><td>Others</td><td>Good retrieval performance</td><td>Poor generality, domain-dependent, poor ENN search, poor support of dynamic update</td><td>++</td></tr></tbody></table>

LHS-based indexing schemes. All in all, while the aforementioned approaches are shown to give satisfactory results for approximate nearest neighbor search, their performance for *exact* nearest neighbor search is not significantly improved.

Finally, some other indexing methods, which rely on improving the features and/or matching step, might work well under some particular domains but they face several difficulties. At first, most of these methods are dedicated to specific applications, where prior knowledge about the domain is often known, and thus it is difficult to apply them to other tasks. Furthermore, since these methods rely on feature handling, the real problem of feature indexing still remains. As a result, such methods tend to perform poorly with respect to the increase of dataset size and data dimensionality.

In the next chapter, we attempt to bring a new contribution for feature indexing in feature vector space. We also empirically demonstrate that the proposed indexing algorithm performs much better for both the tasks of ENN and ANN search, compared to the state-of-the-art KD-tree-based and clustering-based methods.
# Chapter 5

## An efficient indexing scheme based on linked-node m-ary tree (LM-tree)

Following the conclusions given in the previous chapter, it is highly agreed that feature indexing is very important for a real-time image processing system. This chapter presents a new contribution for feature indexing in high-dimensional feature vector space. The proposed algorithm, called linked-node m-ary tree (LM-tree), has many favourable properties that make it different from all the existing methods. Extensive experiments, applied to a wide corpus set of features, have been performed to validate the proposed algorithm in comparison with the state-of-the-art indexing schemes.

### 5.1 Introduction

As we argued earlier, despite the fact that a large number of indexing algorithms have been proposed in the literature, few of them (e.g., LHS-based schemes [Lv et al., 2007], randomized KD-trees [Silpa-Anan and Hartley, 2008], randomized K-medoids clustering trees [Muja and Lowe, 2012], and hierarchical K-means tree [Muja and Lowe, 2009]) have been well validated on extensive experiments to give satisfactory performance on specific benchmarks. Although the mentioned algorithms can produce quickly the answers of approximate nearest neighbors of a given query, their search performances are still limited in the case where a pretty high search precision is desired (e.g., > 90%). Especially, in some applications where exact search is required, these algorithms give little or even no better search performance compared to the brute-force search. These arguments leave a room for advanced indexing algorithms.
In this chapter, an attempt of bringing a new and efficient indexing algorithm in feature vector space is made. Particularly, a linked-node m-ary tree (LM-tree) indexing algorithm is presented, which works really well for both exact and approximate nearest neighbor search. The proposed indexing algorithm consists of three main parts, each of which is described as follows. First, a new polar-space-based method of data decomposition is presented to construct the LM-tree. The new decomposition method employs randomly two axes from a few axes having the highest variance of the underlying dataset to iteratively partition the data into $m$ ($m > 2$) roughly equally sized subsets. This spirit is in contrast to many existing tree-based indexing algorithms, where only one axis is employed for the same purpose. Second, a novel pruning rule is proposed to efficiently narrow down the search space. Furthermore, the computation of the lower bounds is very simple, avoiding the overhead of complicated computation as often seen in many existing approaches. Finally, a bandwidth search method is introduced to explore the nodes of the LM-tree. Experimental results, applied to one million 128-dimensional SIFT features and 250000 960-dimensional GIST features, show that the proposed algorithm gives a significant improvement of search performance, compared to many state-of-the-art indexing algorithms. An additional application to image retrieval is also investigated to further demonstrate the efficiency of the proposed indexing scheme.

## 5.2 The proposed algorithm

The proposed indexing scheme is a tree-based structure, composing of three main components: constructing the LM-tree, doing exact nearest search with the LM-tree, and doing approximate nearest search with the LM-trees. The detail of each component is described in the following.

### 5.2.1 Construction of the LM-tree

Given a dataset $X$ that is composed of $N$ feature vectors in a $D$-dimensional space $\mathbb{R}^D$, we present, in this section, an indexing structure to index the dataset $X$ while supporting an efficient proximity search. For a better presentation of our approach, we use the notation $\mathbf{p}$ as a point in the $\mathbb{R}^D$ feature vector space, and $p_i$ as the $i^{th}$ component of $\mathbf{p}$ ($1 \le i \le D$). We also denote $p = (p_{i_1}, p_{i_2})$ as a point in a $2D$ space. We adopted here the conclusion made in [Silpa-Anan and Hartley, 2008] about the use of PCA for aligning the data before constructing the LM-tree. This approach enables us to partition the data via the narrowest directions. In particular, the dataset $X$ is translated to its centroid following a step of data rotation to make the coordinate axes aligned with the principal axes. Note that no dimension reduction is performed in this step. In fact, PCA analysis is used only to align the data. Next, the LM-tree is constructed by recursively partitioning the dataset $X$ into $m$ roughly equal-sized subsets. Next, we present the main steps of the LM-tree's construction process with respect to the outline in Algorithm 2:

* Sort the axes in decreasing order of variance, and choose randomly two axes, $i_1$ and $i_2$, from the first $L$ highest variance axes ($L < D$).
* Project every point $\mathbf{p} \in X$ into the plane $i_1\mathbf{c}i_2$, where $\mathbf{c}$ is the centroid of the set $X$, and then compute the corresponding angle: $\phi = \arctan(p_{i_1} - c_{i_1}, p_{i_2} - c_{i_2})$.

* Sort the angles $\{\phi_t\}_{t=1}^n$ in increasing order ($n = |X|$), and then divide the angles into $m$ disjointed sub-partitions: $(0, \phi_{t_1}] \cup (\phi_{t_1}, \phi_{t_2}] \cup \dots \cup (\phi_{t_m}, 360]$, each of which contains roughly the same number of elements (e.g., the data points projected into the plane $i_1\mathbf{c}i_2$).

* Partition the set $X$ into $m$ subsets $\{X_k\}_{k=1}^m$ corresponding to $m$ angle sub-partitions obtained in the previous step.

<figure><img src="image_83.png" /><figcaption>Figure 5.1: Illustration of the iterative process of data partitioning in an 2D space: the 1st partitioning is applied to the dataset $X$, and the 2nd partitioning is applied to the subset $X_6$ (the branching factor $m = 6$).</figcaption></figure>

For each subset $X_k$, a new node $T_k$ is constructed and then attached to its parent node, where we also store the following information: the split axes (i.e., $i_1$ and $i_2$), the split centroid ($c_{i_1}, c_{i_2}$), the split angles $\{\phi_{t_k}\}_{k=1}^m$, and the split projected points $\{(p_{i_1}^k, p_{i_2}^k)\}_{k=1}^m$, where the point $(p_{i_1}^k, p_{i_2}^k)$ corresponds to the split angle $\phi_{t_k}$. For efficient access across these child nodes, a direct link is established between two adjacent nodes $T_k$ and $T_{k+1}$ ($1 \le k < m$), and the last one $T_m$ is linked to the first one $T_1$. Next, we repeat this partitioning process for each subset $X_k$ that is associated with the child node $T_k$ until the number of data points in each node falls below a pre-defined threshold $L_{max}$. Figure 5.1 illustrates the first and second levels of the LM-tree construction with a branching factor of $m = 6$.

It is worthwhile pointing out that each time that a partition proceeds, two axes are employed for dividing the data. This approach is in contrast to many existing tree-based indexing algorithms, where only one axis is employed to partition the data. Consequently, as argued in [Silpa-Anan and Hartley, 2008], considering a high-dimensional feature space, such as 128-dimensional SIFT features, the total number of axes that are involved in the tree construction is rather limited, making any pruning rules less-efficient, and the tree is less discriminative for later usage of searching. Naturally, the number of principal axes involved in partitioning the data is proportional to both the search efficiency and precision.
<caption>**Algorithm 2** LMTreeBuilding($X, m, L, L_{max}$)</caption>

1: **Input:** A dataset $X \in R^D$, the branching factor ($m$), the number of highest variance axes to be selected ($L$), and the maximum number of data points in a leaf node ($L_{max}$).
2: **Output:** The LM-tree of representing all the data points.
3: $A \leftarrow$ sort the axes in the decreasing order of variance
4: $t_1 \leftarrow random(L)$ {select randomly a number in $[1, L]$}
5: $t_2 \leftarrow random(L)$
{make sure that $t_1 \neq t_2$}
6: $i_1 \leftarrow A[t_1]$, $i_2 \leftarrow A[t_2]$ {select two axes $i_1$ and $i_2$}
7: $c \leftarrow$ compute the centroid of $X$
8: $i \leftarrow 1$, $n \leftarrow |X|$
{project every point $\mathbf{p} \in R^D$ into the plane $i_1c_{i2}$ to compute $\phi_i$}
9: **while** $i \le n$ **do**
10: &nbsp;&nbsp;&nbsp;&nbsp;**p** $\leftarrow X[i]$ {for each point $\mathbf{p} \in X$}
11: &nbsp;&nbsp;&nbsp;&nbsp;$\phi_i$ = arctan($p_{i1} - c_{i1}, p_{i2} - c_{i2}$) {compute the corresponding angle w.r.t the selected axes $i_1$ and $i_2$}
12: &nbsp;&nbsp;&nbsp;&nbsp;$i \leftarrow i + 1$
13: **end while**
14: $B \leftarrow$ sort the angles $\phi_i$ in the increasing order
15: $\Delta \leftarrow n/m$
16: $X_k \leftarrow \emptyset$ {reset each subset $X_k$ to empty ($1 \le k \le m$)}
{start partitioning $X$ into $m$ equal-sized subsets $\{X_k\}$}
17: **for** each $\mathbf{p} \in X$ **do**
18: &nbsp;&nbsp;&nbsp;&nbsp; $\phi_{\mathbf{p}}$ = arctan($p_{i1} - c_{i1}, p_{i2} - c_{i2}$)
19: &nbsp;&nbsp;&nbsp;&nbsp; $k \leftarrow 1$
20: **while** $k \le m$ **do**
21: &nbsp;&nbsp;&nbsp;&nbsp; **if** $\phi_{\mathbf{p}} \le B[k\Delta]$ **then**
22: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $X_k \leftarrow X_k \cup \{p\}$ {put $\mathbf{p}$ into a proper subset $X_k$}
23: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $k \leftarrow m + 1$ {break the loop}
24: **end if**
25: $k \leftarrow k + 1$
26: **end while**
27: **end for**
{Recursively building the tree for each subset $\{X_k\}$}
28: $k \leftarrow 1$
29: **while** $k \le m$ **do**
30: &nbsp;&nbsp;&nbsp;&nbsp; $T_k \leftarrow$ construct a new node corresponding to $X_k$
31: &nbsp;&nbsp;&nbsp;&nbsp; **if** $|X_k| < L_{max}$ **then**
32: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; LMTreeBuilding($X_k, m, L, L_{max}$) {recursive tree building for each subset $X_k$}
33: **end if**
34: **end while**
### 5.2.2 Exact nearest neighbor search in the LM-tree

Exact nearest neighbor search in the LM-tree is proceeded by using a branch-and-bound algorithm. Given a query point **q**, we first project **q** into a new space by using the principal axes, which is similar to how we processed the LM-tree construction. Next, starting from the root, we traverse down the tree and we use the split information stored at each node to choose the best child node for further exploration. Particularly, given an internal node **u** along with the corresponding split information {$i_1, i_2, c_{i1}, c_{i2}, \{\phi_{t_k}\}_{k=1}^m, \{(p_{i1}^k, p_{i2}^k)\}_{k=1}^m$} which is already stored at **u**, we first compute an angle: $\phi_{q_u} = \arctan(q_{i1} - c_{i1}, q_{i2} - c_{i2})$. Next, a binary search is applied to the query angle $\phi_{q_u}$ over the sequence {$\phi_{t_k}\}_{k=1}^m$ to choose the child node of **u** that is closest to the query **q** for further exploration. This process continues until a leaf node is reached, followed by the partial distance search (PDS) [Cheng et al., 1984, McNames, 2001] to the points contained in the leaf. Backtracking is then invoked to explore the rest of the tree. Algorithm 3 outlines the main steps of this process.

<caption>**Algorithm 3 ENNSearch(u, q, dist<sub>lb</sub>)**</caption>

1: **Input:** A pointer to a current node of the LM-tree (u), a query (q) and the current lower bound distance (dist<sub>lb</sub>)
2: **Output:** The exact nearest neighbor of q
3: **if** u is a leaf node **then**
4:     {p<sub>best</sub>, dist<sub>best</sub>} ← apply sequence search for the points contained in u
5:     **else**
6:         {i<sub>1</sub>, i<sub>2</sub>, c<sub>1</sub>, c<sub>2</sub>, {φ<sub>t<sub>k</sub></sub>}<sub>k=1</sub><sup>m</sup>} ← access split information stored at u
7:         φ<sub>q<sub>u</sub></sub> = arctan(q<sub>i<sub>1</sub></sub> - c<sub>i<sub>1</sub></sub>, q<sub>i<sub>2</sub></sub> - c<sub>i<sub>2</sub></sub>)
8:         s<sub>k</sub> ← find a son s<sub>k</sub> of u where φ<sub>q<sub>u</sub></sub> is contained in (φ<sub>t<sub>k-1</sub></sub>, φ<sub>t<sub>k</sub></sub>]
9:         ENNSearch(s<sub>k</sub>, q, dist<sub>lb</sub>) {explore s<sub>k</sub> first}
10:      m ← the number of sons of u
11:     L<sub>ord</sub> ← ∅ {construct an ordered list of visiting the nodes}
12:     s<sub>left</sub> ← s<sub>k</sub>, s<sub>right</sub> ← s<sub>k</sub>, i ← 1
13:     **while** i ≤ m/2 **do**
14:         L<sub>ord</sub> ← L<sub>ord</sub> ∪ move2left(s<sub>left</sub>) {get the adjacent node in the left}
15:         L<sub>ord</sub> ← L<sub>ord</sub> ∪ move2right(s<sub>right</sub>) {get the adjacent node in the right}
16:         i ← i + 1
17:     **end while**
18:     **for** each node s in L<sub>ord</sub> **do**
19:         {dist<sub>nlb</sub>, q<sub>n</sub>} ← ComputeLB(s, q, dist<sub>lb</sub>) {update new lower bound and query}
20:     **if** dist<sub>nlb</sub> < dist<sub>best</sub> **then**
21:         ENNSearch(s, q<sub>n</sub>, dist<sub>nlb</sub>)
22:     **end if**
23:     **end for**
24: **end if**
25: **return** p<sub>best</sub> {the exact nearest neighbor of q}

Each time when we are positioned at a node **u**, the lower bound is computed as the distance from the query **q** to the node **u**. If the lower bound is higher than the distance
from $\mathbf{q}$ to the nearest point found so far, we can safely avoid exploring this node and proceed with other nodes. In this section, we present a novel rule for efficiently computing the lower bound. Our pruning rule was inspired by the work presented in the principal axis tree (PAT) [McNames, 2001]. PAT is a generalization of the KD-tree, where the page regions are hyper-polygons rather than hyper-rectangles, and the pruning rule is recursively computed based on the law of cosines.

The drawbacks of the PAT's pruning rule are the computational cost of the complexity (i.e., $O(D)$) and the inefficiency when working on a high-dimensional space because only one axis is employed at each partition. As our method of data decomposition (i.e., LM-tree construction) is quite different from that of the KD-tree-based structures, we have developed a significant improvement of the pruning rule used in PAT. Particularly, we have incorporated two following major advantages for the proposed pruning rule:

* The lower bound is computed as simple as in a 2D space, regardless of how large the dimensionality $D$ is. Therefore, the time complexity is just $O(2)$ instead of $O(D)$ as in the case of PAT.

* The magnitude of the proposed lower bound is significantly higher than that in PAT. This enables the proposed pruning rule to work efficiently.

<figure><img src="image_84.png" /><figcaption>Figure 5.2: Illustration of the lower bound computation.</figcaption></figure>

We now come back to the description of computing the lower bound. Let $u$ be the node in the LM-tree at which we are positioned, let $T_k$ be one of the children of $u$ that is going to be searched, and let $p_k = (p_{i_1}^k, p_{i_2}^k)$ be the $k^{th}$ split point, which corresponds to the child node $T_k$ (see Figure 5.2). The lower bound $LB(\mathbf{q}, T_k)$, from $\mathbf{q}$ to $T_k$, is recursively computed from $LB(\mathbf{q}, u)$. The main steps of this process are as follows (i.e., Algorithm 4):

* Compute the angles: $\alpha_1 = \angle qcp_k$ and $\alpha_2 = \angle qcp_{k+1}$, where $q = (q_{i_1}, q_{i_2})$ and $c = (c_{i_1}, c_{i_2})$.
* If one of two angles $α_1$ and $α_2$ is smaller than 90°, then we have the following fact due to the rule of cosines [McNames, 2001]:

$$d(q, x)^2 \geq d(q, h)^2 + d(h, x)^2 \qquad (5.1)$$

where $x$ is any point in the region of $T_k$, and $h = (h_{i_1}, h_{i_2})$ is the projection of $q$ on the line $cp_k$ or $cp_{k+1}$, while accounting for $α_1 ≤ α_2$ or $α_1 > α_2$, and $d(h, x)$ is the Euclidean distance between two points. Then, we apply the rule of computing the lower bound in PAT in an 2D space as follows:

$$LB^2(\mathbf{q}, T_k) \leftarrow LB^2(\mathbf{q}, u) + d(q, h)^2 \qquad (5.2)$$

Next, we treat the point $h = (q_1, q_2, \dots, h_{i_1}, \dots, h_{i_2}, \dots, q_{D-1}, q_D)$ in place of $\mathbf{q}$ by the means of lower bound computation for the descendant of $T_k$.

* If both angles, $α_1$ and $α_2$, are higher than 90° (e.g., the point $q_2$ in Figure 5.2), we have a more restricted rule as follows:

$$d(q, x)^2 \geq d(q, c)^2 + d(c, x)^2 \qquad (5.3)$$

Therefore, the lower bound is easily computed as:

$$LB^2(\mathbf{q}, T_k) \leftarrow LB^2(\mathbf{q}, u) + d(q, c)^2 \qquad (5.4)$$

Again, we treat the point $c = (q_1, q_2, \dots, c_{i_1}, \dots, c_{i_2}, \dots, q_{D-1}, q_D)$ in place of $\mathbf{q}$ by the means of lower bound computation for the descendant of $T_k$.

As the lower bound $LB(\mathbf{q}, T_k)$ is recursively computed from $LB(\mathbf{q}, u)$, an initial value must be set for the lower bound at the root node. Obviously, we set $LB(\mathbf{q}, root) = 0$. It is also noted that when the point $\mathbf{q}$ is fully contained in the region of $T_k$, no computation of the lower bound is required. Therefore: $LB(\mathbf{q}, T_k) \leftarrow LB(\mathbf{q}, u)$.

### 5.2.3 Approximate nearest neighbor search in the LM-tree

Approximate nearest neighbor search is proceeded by constructing multiple randomized LM-trees to account for different viewpoints of the data. To achieve the optimized utilization of memory space, the data points are stored in a common buffer, while the leaf nodes of the trees contain the indexes to these actual data points. The idea of using multiple randomized trees for ANN search was originally presented in [Silpa-Anan and Hartley, 2008], where the authors proposed to construct multiple randomized KD-trees. This technique was then incorporated with the priority search and successfully used in many other tree-based structures [Muja and Lowe, 2009], [Muja and Lowe, 2012]. Although the priority search was shown to give better search performance, it is certainly subjected to a high computational cost because the process of maintaining a priority queue during the online search is rather expensive.

Here, we exploit the advantages of using multiple randomized LM-trees but without using the priority queue. The basic idea is to restrict the search space to the branches
<caption>**Algorithm 4 ComputeLB(u, q, dist<sub>lb</sub>)**</caption>

1: **Input:** A pointer to a node of the LM-tree (u), a query (q) and the current lower bound distance (dist<sub>lb</sub>)
2: **Output:** New lower bound (dist<sub>nlb</sub>) and new query point (q<sub>n</sub>)
3: p<sub>u</sub> ← get the parent node of u
4: {i<sub>1</sub>, i<sub>2</sub>, c<sub>1</sub>, c<sub>2</sub>} ← access split information stored at p<sub>u</sub>
5: α<sub>min</sub> ← min{α<sub>1</sub>, α<sub>2</sub>} {see explanation of α<sub>1</sub> and α<sub>2</sub> in the text}
6: q<sub>n</sub> ← q
7: **if** α<sub>min</sub> ≤ 90<sup>0</sup> **then**
8:     **h** ← projection of q on the split axis {see explanation in the text}
9:     q<sub>n</sub>[i<sub>1</sub>] ← h<sub>i<sub>1</sub></sub> {update the new query at the coordinates i<sub>1</sub> and i<sub>2</sub>}
10:    q<sub>n</sub>[i<sub>2</sub>] ← h<sub>i<sub>2</sub></sub>
11:    dist<sub>nlb</sub> ← dist<sub>lb</sub> + (q<sub>i<sub>1</sub></sub> − h<sub>i<sub>1</sub></sub>)<sup>2</sup> + (q<sub>i<sub>2</sub></sub> − h<sub>i<sub>2</sub></sub>)<sup>2</sup> {update the new lower bound distance}
12: **else**
13:    q<sub>n</sub>[i<sub>1</sub>] ← c<sub>i<sub>1</sub></sub> {update the new query based on the centroid}
14:    q<sub>n</sub>[i<sub>2</sub>] ← c<sub>i<sub>2</sub></sub>
15:    dist<sub>nlb</sub> ← dist<sub>lb</sub> + (q<sub>i<sub>1</sub></sub> − c<sub>i<sub>1</sub></sub>)<sup>2</sup> + (q<sub>i<sub>2</sub></sub> − c<sub>i<sub>2</sub></sub>)<sup>2</sup> {update the new lower bound distance}
16: **end if**
17: **return** {dist<sub>nlb</sub>, q<sub>n</sub>}

that are not very far from the considering path. In this way, we introduce a specific search procedure, so-called *bandwidth* search, which proceeds by setting a search bandwidth to every intermediate node of the ongoing path. Particularly, let $P = \{u_1, u_2, \dots, u_r\}$ be a considering path obtained by traversing down a single LM-tree, where $u_1$ is the root node and $u_r$ is the current node of the path. The proposed bandwidth search indicates that for each intermediate node $u_i$ of $P$ ($1 \le i \le r$), every sibling node of $u_i$ at a distance of $b+1$ nodes ($1 \le b < m/2$) on both sides from $u_i$ does not need to be searched. The value $b$ is called search bandwidth. Taking one example as shown in Figure 5.3, where $X_6$ is an intermediate node on the considering path $P$, then only $X_1$ and $X_5$ are candidates for further inspection given a search bandwidth of $b=1$.

There is a notable point that when the projected query $q$ is too close to the projected centroid $c$, all of the sibling nodes of $u_i$ should be inspected as in the case of an ENN search. Particularly, this scenario occurs at a certain node $u_i$ if $d(q, c) \le \epsilon D_{med}$, where $D_{med}$ is the median value of the distances between $c$ and all of the projected data points that are associated with $u_i$, and $\epsilon$ is a tolerance radius parameter. In addition, in order to obtain a varying range of the search precision, we would need a parameter $E_{max}$ for the maximum data points to be searched on a single LM-tree. Algorithm 5 outlines the flowchart of our bandwidth search process.

As we are designing an efficient solution dedicated to an ANN search, it would make sense to use an approximate pruning rule rather than an exact rule. Particularly, we have used only formula (5.4) as an approximate pruning rule. This adaption offers two favourable features. First, it reduces much of the computational cost. Recall that rule (5.4) requires the computation of $d(q, c)$ in a 2D space, where point $c$ has been already computed during the offline tree construction. The computation of this rule is thus very
<figure><img src="image_85.png" /><figcaption>Figure 5.3: Illustration of our bandwidth search with $b = 1$: $X_6$ is an intermediate node of the considering path, and its adjacent sibling nodes, $X_1$ and $X_5$, will also be searched; if $\mathbf{q}$ is too close to the centroid (e.g., inside the circle for the case of $X_6$), then all of the sibling nodes of $X_6$ will be searched.</figcaption></figure>

efficient. Second, it also ensures that a larger fraction of nodes will be inspected but few of them would be actually searched after checking the lower bound. In this way, it increases the chance of reaching the true nodes that are closest to the query. More generally, we have adapted formula (5.4) as follows:

$$ LB^2(\mathbf{q}, T_k) \leftarrow \kappa \cdot (LB^2(\mathbf{q}, u) + d(q, c)^2) \qquad (5.5) $$

where $\kappa \ge 1$ is the pruning factor that controls the rate of pruning the branches in the trees. This factor can be adaptively estimated during the tree construction given a specific precision and a specific dataset. However, we have set this value as $\kappa = 2.5$ and we show, in our experiments, that it is possible to achieve satisfactory search performance on many datasets. Algorithm 6 sketches the basic steps of computing the approximate pruning rule.

## 5.3 Experimental results

We evaluated our system versus several representative fast proximity search systems in the literature, including randomized KD-trees (RKD-trees) [Silpa-Anan and Hartley, 2008] that use the priority search implementation in [Muja and Lowe, 2009], hierarchical K-means tree (K-means tree) [Muja and Lowe, 2009], randomized K-medoids clustering trees (RC-trees) [Muja and Lowe, 2012], and the multi-probe LSH algorithm [Lv et al., 2007]. These indexing systems were well-implemented and widely used in the literature thanks to the open source FLANN library<fnref n="1" />. The source code of our system is also publicly available at this address<fnref n="2" />. Note that the partial distance search was implemented in these systems

<footnote n="1">¹http://www.cs.ubc.ca/~mariusm/index.php/FLANN/FLANN</footnote>

<footnote n="2">²https://sites.google.com/site/ptalmtree/</footnote>
**Algorithm 5 ANNSearch(u, **q**, dist<sub>lb</sub>, E<sub>max</sub>, ε, κ, b)**

1: **Input:** A pointer to a current node of the LM-tree (u), a query (**q**), the current lower bound distance (dist<sub>lb</sub>), and the 4 parameters E<sub>max</sub>, ε, κ, and b.

2: **Output:** The approximate nearest neighbor of **q**

3: **if** u is a leaf node **then**

4: {**p**<sub>best</sub>, dist<sub>best</sub>} ← apply sequence search for the points contained in u

5: **else**

6: {i<sub>1</sub>, i<sub>2</sub>, c<sub>1</sub>, c<sub>2</sub>, D<sub>med</sub>, {φ<sub>t<sub>k</sub></sub>}<sup>m</sup><sub>k=1</sub>} ← access split information stored at u

7: φ<sub>q<sub>u</sub></sub> = arctan(q<sub>i<sub>1</sub></sub> - c<sub>1</sub>, q<sub>i<sub>2</sub></sub> - c<sub>2</sub>)

8: s<sub>k</sub> ← find a son s<sub>k</sub> of u where φ<sub>q<sub>u</sub></sub> is contained in (φ<sub>t<sub>k-1</sub></sub>, φ<sub>t<sub>k</sub></sub>]

9: ANNSearch(s<sub>k</sub>, **q**, dist<sub>lb</sub>, E<sub>max</sub>, ε, κ, b) {explore s<sub>k</sub> first}

10: m ← the number of sons of u

11: dist ← (q<sub>i<sub>1</sub></sub> - c<sub>1</sub>)<sup>2</sup> + (q<sub>i<sub>2</sub></sub> - c<sub>2</sub>)<sup>2</sup>

12: m<sub>bd</sub> ← b {compute the search bandwidth}

13: **if** dist < (εD<sub>med</sub>)<sup>2</sup> **then**

14: m<sub>bd</sub> ← m/2

15: **end if**

16: L<sub>ord</sub> ← ∅ {construct an ordered list of visiting the nodes}

17: s<sub>left</sub> ← s<sub>k</sub>, s<sub>right</sub> ← s<sub>k</sub>

18: i ← 1

19: **while** i ≤ m<sub>bd</sub> **do**

20: L<sub>ord</sub> ← L<sub>ord</sub> ∪ move2left(s<sub>left</sub>) {get the adjacent node in the left}

21: L<sub>ord</sub> ← L<sub>ord</sub> ∪ move2right(s<sub>right</sub>) {get the adjacent node in the right}

22: i ← i + 1

23: **end while**

24: **for** each node s in L<sub>ord</sub> **do**

25: {dist<sub>nlb</sub>, **q**<sub>n</sub>} ← ComputeLB_ANN(s, **q**, dist<sub>lb</sub>, κ) {compute the lower bound}

26: **if** dist<sub>nlb</sub> < dist<sub>best</sub> **then**

27: ANNSearch(s, **q**<sub>n</sub>, dist<sub>nlb</sub>, E<sub>max</sub>, ε, κ, b)

28: **end if**

29: **end for**

30: **end if**

31: **return** **p**<sub>best</sub> {the approximate nearest neighbor of **q**}

**Algorithm 6 ComputeLB_ANN(u, **q**, dist<sub>lb</sub>, κ)**

1: **Input:** A pointer to a node of the LM-tree (u), a query (**q**), the current lower bound distance (dist<sub>lb</sub>) and the pruning factor (κ)

2: **Output:** New (approximate) lower bound (dist<sub>nlb</sub>) and new query point (**q**<sub>n</sub>)

3: p<sub>u</sub> ← get the parent node of u

4: {i<sub>1</sub>, i<sub>2</sub>, c<sub>1</sub>, c<sub>2</sub>} ← access split information stored at p<sub>u</sub>

5: **q**<sub>n</sub> ← **q**

6: **q**<sub>n</sub>[i<sub>1</sub>] ← c<sub>i<sub>1</sub></sub> {update the new query at the coordinates i<sub>1</sub> and i<sub>2</sub>

7: **q**<sub>n</sub>[i<sub>2</sub>] ← c<sub>i<sub>2</sub></sub>

8: dist<sub>nlb</sub> ← dist<sub>lb</sub> + κ · ((q<sub>i<sub>1</sub></sub> - c<sub>i<sub>1</sub></sub>)<sup>2</sup> + (q<sub>i<sub>2</sub></sub> - c<sub>i<sub>2</sub></sub>)<sup>2</sup>) {update the new lower bound distance}

9: **return** {dist<sub>nlb</sub>, **q**<sub>n</sub>}
in order to improve the efficiency of the sequence search at the leaf nodes. Two datasets, ANN\_SIFT1M and ANN\_GIST1M [Jégou et al., 2011], were used for all of the experiments. The ANN\_SIFT1M dataset contains a database of one million 128-dimensional SIFT features and a test set of 5000 SIFT features, while the ANN\_GIST1M dataset is composed of a database of one million 960-dimensional GIST features and a test set of 1000 GIST features. Because the dimensionality of the GIST feature is very high and our computer configuration is limited (i.e., Windows XP, 2.4G RAM), we were not able to load the full ANN\_GIST1M dataset into memory. Consequently, we have used 250000 GIST features for search evaluation. Following the convention of the evaluation protocol used in the literature [Beis and Lowe, 1997], [Muja and Lowe, 2009], [Muja and Lowe, 2012], we computed the search precision and search time as the average measures obtained by running 1000 queries taken from the test sets of the two datasets. To make the results independent on the machine and software configuration, the speedup factor is computed relative to the brute-force search. The details of our experiments are presented in this section.

### 5.3.1 ENN search evaluation

For ENN search, we use a single LM-tree and set the parameters that are involved in the LM-tree as follows: $L_{max} = 10$, $L = 2$, and $m \in \{6, 7\}$ with respect to the GIST and SIFT datasets (see section 5.2.1). By setting $L = 2$, we choose exactly the two highest variance axes at each level of the tree for data partitioning. We compared the performance of the ENN search of the following systems: the proposed LM-tree, the KD-tree, and the hierarchical K-means tree. Figure 5.4(a) shows the speedup over the brute-force search for the three systems, when applied to the SIFT datasets with different sizes. We can note that the LM-tree outperforms the other two systems on all of the tests. Figure 5.4(b) presents the search performance for the three systems for the GIST features. The proposed LM-tree again outperforms the others and even performs far better than the SIFT features. Taking the test where #Points = 150000 on Figure 5.4(b), for example, the LM-tree gives a speedup of 17.2, the KD-tree gives a speedup of 3.53, and the K-means tree gives a speedup of 1.42 over the brute-force search. These results confirm the efficiency of the LM-tree for the ENN search relative to the two baseline systems.

### 5.3.2 ANN search evaluation

For the ANN search, we adopted here the benefit of using multiple randomized trees, which were successfully used in the previous works of [Silpa-Anan and Hartley, 2008] and [Muja and Lowe, 2012]. Hence, we set the required parameters as follows: $L_{max} = 10$, $L = 8$, $b = 1$, and $m \in \{6, 7\}$. By setting $L = 8$, we choose randomly two axes from the eight highest variance axes at each level of the tree for data partitioning. In addition, the parameter $b = 1$ indicates that our bandwidth search process visits three adjacent nodes (including the node in question) at each level of the LM-tree. Four systems participated in this evaluation, including the proposed LM-trees, RKD-trees, RC-trees, and K-means tree. Following the conclusion made in [Muja and Lowe, 2012] about the optimal number of parallel trees for ANN search, we used 8 parallel trees in the first three systems as the search precision is expected to be high (> 70%) in our experiments. It is worth explaining
<figure><img src="image_86.png" /><figcaption>Figure 5.4: Exact search performance for the SIFT (a) and GIST (b) features.</figcaption></figure>

a bit the meaning of the search precision in our context. For 1-NN ANN search, the search precision is computed as the ratio of the number of exact answers to the number of queries. In our experiments, we performed 1000 queries for each test and if one system achieves a precision of 90% implying that it produces 900 exact answers out of 1000 queries.

It is also noted that we used a single tree for the K-means tree indexing algorithm because it was shown in [Muja and Lowe, 2009] that the use of multiple K-means trees does not give better search performance. For all of the systems, the parameters $E_{max}$ and $\epsilon$ (i.e., the LM-tree) are varied to obtain a wide range of search precision.

<figure><img src="image_87.png" /><figcaption>Figure 5.5: Approximate search performance for the SIFT (a) and GIST (b) features.</figcaption></figure>

Figure 5.5(a) shows the search speedup versus the search precision of the four systems
for 1 million SIFT features. As can be seen, the proposed LM-trees algorithm gives significantly better search performance everywhere compared with the other systems. When considering the search precision of 95%, for example, the speedups over a brute-force search of the LM-trees, RKD-trees, RC-trees, and K-means tree are 167.7, 108.4, 122.4, and 114.5, respectively. To make it comparable with the multi-probe LSH indexing algorithm, we converted the real SIFT features to the binary vectors and tried several parameter settings (i.e., the number of hash tables, the number of multi-probe levels, and the length of the hash key) to obtain the best search performance. However, the result obtained on one million SIFT vectors is rather limited. Taking the search precision of 74.7%, for example, the speedup over a brute-force search (using Hamming distance) is only 1.5.

Figure 5.5(b) shows the search performance of all of the systems for 200000 GIST features. Again, the LM-trees algorithm clearly outperforms the others and tends to perform much better than the SIFT features. The RC-trees algorithm also works reasonably well, while the RKD-trees and K-means tree work poorly for this dataset. Considering the search precision of 90%, for example, the speedups over a brute-force search of the LM-trees, RKD-trees, RC-trees, and K-means tree are 113.5, 15.0, 45.2, and 21.2, respectively.

<figure><img src="image_88.png" /><figcaption>Figure 5.6: ANN search performance as a function of the dataset size: (a) the SIFT datasets (search precision = 96%); (b) the GIST datasets (search precision = 95%).</figcaption></figure>

In Figure 5.6, we present the ANN search performance of the four systems as a function of the dataset size. For this purpose, the search precision is set to a rather high degree, especially at 96% and 95% for the SIFT and GIST features, respectively. This time, the LM-trees algorithm still gives a substantially better search performance than the others and tends to perform quite well, considering the increase in the dataset size. For the SIFT features, the RC-trees algorithm works reasonably well, except for the point where #Points = 800K, at which its search performance is noticeably degraded. It is also noted that the speedups of the RKD-trees and K-means tree for the GIST features are quite low, even lower than the speedup of the LM-tree for the ENN search (see Figure 5.4(b)).

Three crucial factors explain these outstanding results of the LM-trees. First, the
use of the two highest variance axes for data partitioning in the LM-tree gives a more discriminative representation of the data in comparison to the common use of the sole highest variance axis as in the literature. Second, by using the approximate pruning rule, a larger fraction of nodes will be inspected, but many of them would be eliminated after checking the lower bound. In this way, the number of data points that will be actually searched, is retained under the pre-defined threshold $E_{max}$, while covering a larger number of inspected nodes, and thus increasing the chance of reaching the true nodes that are closest to the query. Finally, the use of bandwidth search gives much benefit in terms of the computational cost, compared to the priority search that is used in the baseline indexing systems.

The last experiment is presented on Figure 5.7 to evaluate the distance error (i.e., error ratio) of the approximate answers to the exact ones. By approximate nearest neighbor search, it means that not all the answers are the exact ones. More precisely, the quality of approximate nearest neighbor search is evaluated by two metrics: the search precision and the distance error ratio. Taking the search precision of 95% provided by one indexing system, for instance, implies that the system produces 950 exact answers out of 1000 queries (e.g., considering a test of 1000 queries). It also tells us that there are 50 answers which are not the exact ones but the approximate nearest neighbors. The quality of these approximate answers are evaluated by the distance error ratio to the exact ones. [Gionis et al., 1999] defined the distance error ratio as follows:

$$ dist_{err} = \frac{1}{Q} \sum_{i=1}^{Q} \frac{dist(p_i, q_i)}{dist(p_i^*, q_i)} \qquad (5.6) $$

where $p_i$ is the approximate answer of the query $q_i$, whereas $p_i^*$ is the exact nearest neighbor of $q_i$, and $dist(\cdot, \cdot)$ is the Euclidean distance between two points.

Figure 5.7 shows the distance error ratio with respect to the increase of search precision of all the systems for the SIFT and GIST features. For both cases, although the proposed system does not always perform best, the distance error ratio is quite low, especially for the GIST features. On average, the distance error ratios of the proposed system are 1.0414 and 1.0193 for the SIFT and GIST features, respectively. In other words, the approximate answers given by the proposed system are, on average, located on the hyper-spheres with the radii of 4.1% and 1.9% higher than that of the exact nearest neighbors for the SIFT and GIST features, respectively.

### 5.3.3 Parameter tuning

In this section, we study the effects of the parameters that are involved in the LM-tree on the search performance. Two types of parameters are concerned in the LM-tree: static and dynamic parameters. The static parameters are those that are involved in the offline phase of building the LM-tree, including the maximum number of data points at a leaf node ($L_{max}$), the number of axes having the highest variance ($L$), and the branching factor ($m$). The dynamic parameters are used in the online phase of searching, including the search bandwidth ($b$), the maximum number of data points to be searched in a single LM-tree ($E_{max}$), the pruning factor ($\kappa = 2.5$), and the tolerance radius $\epsilon$. For the static
<figure><img src="image_89.png" /><figcaption>Figure 5.7: Evaluation of distance error ratio for ANN search applied to SIFT (a) and GIST (b) features.</figcaption></figure>

parameters, our investigation, obtained on the experiments performed so far, has shown that they produce a negligible effect on the search performance provided an appropriate setting of the dynamic parameters. Specifically, it was found that the following settings for the static parameters often lead to relatively good search performance: $L_{max}$ = 10, L = 8, and m ∈ {6, 7}. As a rule of thumb, the branching factor m should be high for a large-scale dataset and vice versa. The dynamic parameters ($E_{max}$, ϵ, b) are used to achieve a given specific search precision. They are thus treated as precision-driven parameters. This implies that given a specific precision P (0% ≤ P ≤ 100%), we can design a method to determine automatically the optimized settings for these parameters to achieve the best search time. The basic idea is that given a specific setting of b and ϵ, the parameter $E_{max}$ is estimated by using a binary search. This approach enables the method to work very efficiently. In our case, we have set the parameter b to 1 and designed the following procedure to estimate the optimized setting for $E_{max}$ and ϵ:

* Step 1: Sample the parameter ϵ into discrete values: {$ϵ_0, ϵ_0 + Δ, ..., ϵ_0 + lΔ$}. In our implementation, we set: ϵ₀ = 0, Δ = 0.04, l = 20.

* Step 2: For each value $ϵ_i = ϵ_0 + iΔ$ ($0 ≤ i ≤ l$):

- Step 2(a): Estimate an initial value for $E_{max}$ by running the approximate search procedure without consideration of the parameter $E_{max}$. In this way, the search procedure terminates early with respect to the current settings of $ϵ_i$ and b. Assume Q be the number of searched points during this process.

- Step 2(b): Compute the precision $P_i$ and speedup $S_i$ by using the groundtruth information. If $P_i < P$, then proceed with a new value of $ϵ_{i+1}$ and go to Step 2(a). Otherwise, go to Step 2(c).
- Step 2(c): If $P_i \ge P$, we apply a binary search to find the optimized value of $E_{max}$ in the range of $[0, Q]$. Particularly, we first set $E_{max} = Q/2$ and run the approximate search procedure. Next, the search range is updated as either $[0, \frac{Q}{2}]$ or $[\frac{Q}{2}, Q]$ taking into account that $P_i > P$ or $P_i \le P$. This process continues until the search range is unit-length.

- Step 2(d): Update the best speedup, the parameter $\epsilon_i$, and the optimized parameter $E_{max}$ obtained from Step 2(c). Proceed with a new value of $\epsilon_{i+1}$ and go to Step 2(a) to find the better solution.

* Step 3: Return the parameters $\epsilon_i$ and $E_{max}$, corresponding to the best speedup found so far.

The rationale of setting the parameter $b$ to 1 was inspired from the work of multi-probe LSH [Lv et al., 2007]. When using the two highest variance axes for partitioning the data, two close points can be divided into two adjacent bins. It is thus necessary to have a look at the adjacent bins while exploring the tree. Our experiments have revealed that the parameter $b = 1$ is often a good setting for obtaining satisfactory search performance.

<figure><img src="image_90.png" /><figcaption>Figure 5.8: Search performance as a function of the pruning factor: precision is set to 95%, 90%, and 80% for both the SIFT and GIST features.</figcaption></figure>

Figure 5.8 shows the search performance as a function of the pruning factor $\kappa$. In this experiment, we study the effect of $\kappa$ for a wide range of search precision. For this purpose, the precision is set to 95%, 90%, and 80% for both the SIFT and GIST features. It can be seen that for both types of features, the speedup increases strongly with respect to the increase in $\kappa$ to a certain extent. For example, taking the curve that corresponds to the precision $P = 90\%$ for the SIFT features, the speedup starts to decrease for $\kappa \ge 3.5$. Figure 5.8 also reveals that when increasing $\kappa$, we achieve much of the speedup for low precision
(P = 80%) compared with that for higher precision (P ≥ 90%). This relationship is re-
ferred to as the problem of over-pruning because a large number of branches has been elimi-
nated as κ increases. Hence, the search process would miss some true answers. Fortunately,
this matter can be resolved by using a cross-validation process [Muja and Lowe, 2009] dur-
ing the tree construction to adaptively select an appropriate value for κ, when provided a
specific precision and dataset. However, we have not employed this approach here, and we
set κ = 2.5 for all of our experiments.

**5.4 Application to image retrieval**

In this section, an application to image retrieval is investigated using the proposed
LM-tree indexing scheme. We carry out here a classical approach composing of three
main steps: feature extraction, indexing, and retrieval. For this purpose, we selected a
wide corpus set<fnref n="3" /> of historical books containing ornamental graphics. The dataset, called
"Lettrine", is composed of 23062 isolated graphical images that were extracted from old
documents. Some examples of the drop caps in this dataset are shown in Figure 5.9.

<figure><img src="image_91.png" /><figcaption>Figure 5.9: Several example images of ornamental drop caps in the dataset.</figcaption></figure>

Our challenge here is to prove that the proposed LM-tree indexing algorithm can be
used in the context of a CBIR system. Hence, for the first step of feature extraction, we
selected the GIST feature as it is a commonly used descriptor in the literature for image
retrieval [Oliva and Torralba, 2001]. The GIST descriptor comprises a set of perceptual
properties that represent the dominant spatial structure of a scene such as naturalness,
openness, roughness, expansion, and ruggedness. The GIST descriptor has been widely
used in the literature [Kulis and Grauman, 2009, Jégou et al., 2011] on image retrieval,
scene recognition and classification because of its distinctiveness and efficiency. We used
the original implementation of the 512-dimensional GIST descriptor that was provided by
the authors at this address<fnref n="4" />.

In order to use the GIST descriptor for image matching and retrieval, it is necessary
to normalize the image size. This action was already addressed in the GIST descriptor's
implementation by centering, cropping and resizing the image such that the normalized
image preserves the size ratio of the original image. In our experiments, the common size

<footnote n="3">³http://www.bvh.univ-tours.fr/</footnote>

<footnote n="4">⁴http://people.csail.mit.edu/torralba/code/spatialenvelope/</footnote>
for normalization is set to 256 × 256. Among the 23062 ornamental graphics, 500 images
were included in the query set, and the remaining served as the database set. Next, the
GIST features are computed once for all of the images in an offline phase.

In the second step, all of the database GIST features are indexed by using our LM-tree
algorithm. As this application is dedicated to the image retrieval domain, the main goal
here is to illustrate that the system can produce sufficiently relevant results with a critical
constraint of fast processing time. Therefore, it makes sense to apply an ANN indexing
algorithm rather than the ENN method. Here, we keep the same parameter settings for
the LM-trees as before. In other words, we used 8 parallel LM-trees, and each of which is
associated with the same parameter configuration, as follows: Lmax = 10, L = 8, b = 1,
and m = 6.

For performance evaluation, it would be interesting to use standard evaluation metrics such as precision and recall. However, because there is no groundtruth information that is included in this dataset, the use of precision and recall is not possible. Instead, we used an alternative metric to quantify the retrieval performance of our system. This metric was introduced in [Kulis and Grauman, 2009] for the same purpose as ours. Its basic idea is to measure how well a retrieval system can approximate an ideal linear scan (i.e., a brute-force search). Specifically, we computed the fraction of the common answers between our retrieval system and the ideal linear scan to the number of answers of our system. Figure 5.10 (a) quantitatively shows how well our retrieval system approaches the ideal linear scan. These quantitative results are computed using the top 1, 5, 10, 15, 20, and 25 ranked nearest neighbors (NNs) of our system. To have a more detailed analysis of the results, we can derive hereafter several key remarks:

* For the top 5-NNs of the LM-trees, 93.1% of the time, the retrieved answers are covered by that of the top 15-NNs of the ideal linear scan.

* For the top 15-NNs of the LM-trees, 80.6% of the time, the retrieved answers are covered by that of the top 30-NNs of the ideal linear scan.

* For the top 25-NNs of the LM-trees, 74.4% of the time, the retrieved answers are covered by that of the top 50-NNs of the ideal linear scan.

The results presented in Figure 5.10 (a) quantitatively show the quality of our retrieval
system, which can be regarded as a function of search precision over recall. Furthermore,
we want to measure how fast the system is with respect to these results. For this purpose,
Figure 5.10 (b) shows the speedup of our system relative to the ideal linear scan. In this
test, both of the systems are evaluated by using the same parameter *k* for the number of
nearest neighbors. For more detail, we extract hereafter several key results of our system:

* The 5-NNs LM-trees achieved a speedup factor of 113.8 relative to the 5-NNs ideal linear scan, given a precision of 88.1%.

* The 15-NNs LM-trees achieved a speedup factor of 92.8 relative to the 15-NNs ideal linear scan, given a precision of 80.1%.

* The 50-NNs LM-trees achieved a speedup factor of 78.2 relative to the 50-NNs ideal linear scan, given a precision of 67.1%.
<figure><img src="image_92.png" /><figcaption>Figure 5.10: Search quality of the k-NN LM-trees (a), and search quality versus speedup of the k-NN LM-trees (b).</figcaption></figure>

With regard to the results presented in Figure 5.10 (b), we report in Table 5.1 the absolute search time (ms) and the fraction of searched points for the case of the 5-NNs LM-trees, which were averaged over the 500 queries. Taking a search precision of 73.9% for example, our system must explore 0.78% of the whole database and takes only 0.3 (ms) to return the top 5-NN answers.

<table><caption>Table 5.1: Report of search time and fraction of searched points for the 5-NN LM-trees.</caption><thead><tr><th rowspan="2">Performance</th><th colspan="7">Search precision (%)</th></tr><tr><th>59.9</th><th>61.8</th><th>68.1</th><th>73.9</th><th>78.6</th><th>88.1</th><th>96.5</th></tr></thead><tbody><tr><th>Search fraction (%)</th><td>0.27</td><td>0.32</td><td>0.52</td><td>0.78</td><td>1.17</td><td>2.45</td><td>6.21</td></tr><tr><th>Mean search time (ms)</th><td>0.24</td><td>0.25</td><td>0.26</td><td>0.30</td><td>0.36</td><td>0.38</td><td>0.60</td></tr></tbody></table>

Figure 5.11 shows some examples of our retrieval results in comparison with the ideal linear scan. These retrieval results are obtained using the 5-NN LM-trees that are associated with the precision of 78.6% (see Table 5.1).

## 5.5 Discussion

In this chapter, a novel and efficient indexing algorithm in feature vector space has been presented. Three main features are attributed to the proposed LM-tree. First, a new polar-space-based method of data decomposition has been presented to construct the LM-tree. This new decomposition method differs from the existing methods in the literature in that the two highest variance axes of the underlying dataset are employed to iteratively partition the data. Second, a novel pruning rule is proposed to quickly eliminate the search paths that are unlikely to contain good candidates of the nearest
<figure><img src="image_93.png" /><figcaption>Figure 5.11: Few examples of the retrieval results: for each query in the left, the top ranked 5-NNs are showed for the ideal linear scan (the top row) and our LM-trees (the bottom row).</figcaption></figure>

neighbors. Furthermore, the lower bounds are easily computed, as if the data were in 2D
space, regardless of how high the dimensionality is. Finally, a bandwidth search method is
presented to explore the nodes of the LM-tree. Its basic idea is inspired by the fact that
only a limited number of relevant bins are searched to avoid the overhead of complexity.
The proposed LM-tree has been validated on one million SIFT features and 250000 GIST
features, demonstrating that it works very well for both ENN and ANN search, compared
to the state-of-the-art indexing algorithms, including randomized KD-trees, hierarchical
K-means tree, randomized clustering trees, and the multi-probe LHS scheme.

For further improvements to this work, more experiments on binary features (e.g., BRIEF [Calonder et al., 2010], ORB [Rublee et al., 2011]) would be interesting to evaluate the proposed LM-tree. Dynamic insertion and deletion of the data points in the LM-tree would be also investigated to make the system adaptive to dynamic changes in the data. Automatic parameter tuning using cross-validation approach could be also integrated to make the system more robust and well-fitting to specific datasets. The use of many high variance axes (e.g., more than two) could be considered to study the new behaviors of the system. In addition to these open works, the study of designing an indexing system, which can work on the data stored on an external disk, will be investigated to deal with extremely large datasets in which they are often not fully loaded into the main memory.
# Conclusions

In this chapter, we summarize the main contributions of this dissertation for junction detection in line-drawing images and feature indexing in high-dimensional feature vector space. The favourable features and the shortcomings of the two contributions are carefully discussed. Possible lines of future research for each of these contributions are also given.

Beginning with the chapter 1, we have provided the description of the state-of-the-art methods for junction detection. The merits and limitations of each method are discussed in detail. The connections among these methods are studied to establish potential links to our contribution on junction detection.

Based on the discussion of the existing approaches for junction detection, a novel approach is presented in chapter 2 for robust and accurate junction detection and characterization in line-drawing images. The proposed approach has been evaluated through extensive experiments in which we achieved better results compared to other baseline methods. Besides, the computation complexity of the proposed system has been analyzed from the theoretical point of view, showing that its complexity is essentially linear to the image size. An application of symbol localization has been also investigated, confirming very good results in terms of both detection rate and efficiency. In short, the proposed approach has the following major features:

* Junction distortion avoidance: The problem of junction distortion is avoided by completely removing all distorted zones from foreground. For this purpose, an efficient algorithm is presented to detect and remove crossing regions. In this way, the proposed approach works on the remaining line segments only.

* Accurate junction detection: A new algorithm based on linear least square technique is presented to correctly determine local scales (or regions of support) of every median point. Junction localization is achieved by a novel optimization algorithm which analyzes the local structure relations among the characterized line segments to produce precise junction detection.
* Multiple junction detection: The proposed approach can deal with the problem of multiple detection of junction points for a given complicated crossing zone. This was accomplished by iteratively clustering the incident branches into different groups, each of which forms an optimized location of the junction.

* Efficiency: Computation complexity for the whole process of junction detection is $O(MN + k^2S)$ where M, N are the dimensions of input image, S is the number of median points, and k is a small constant value. In addition, the detected junctions are characterized, classified, and matched in very efficient way that makes them suitable to deal with high complexity problems such as image indexing, symbol recognition, and symbol spotting.

* Robustness: The junction detector is stable to common transformations such as rotation, scaling, translation, and can resist a satisfactory level of degradation/noise. Furthermore, the proposed approach requires no prior knowledge about document content and is independent on any vectorization systems.

* Usefulness: The detected junctions can be used to address different applications of symbol localization and spotting, vectorization, and engineering document retrieval. As the junction detector is robust and accurate, the information of the detected junctions can be used to support different tasks including junction matching, geometry consistency checking, and primitive detection. Such favourable features are highly expected to address the time-critical applications.

In addition to these positive characteristics, we are also aware of several shortcomings of the proposed approach. First, as this approach is dedicated to working with line-like primitives, its performance would be degraded if applying to filled-shape objects, such as logo images. The reason is due to the fact that median lines are not an appropriate means for representing these filled shapes. For the filled shapes, as the local line-thickness at each crossing zone is often high, a large part of median lines is eliminated around the crossing zones. For this reason, we have not enough information to reconstruct the junctions. Second, although we have improved the stage of region of support determination to make it more robust to the digitalization effect, the step of dominant point detection is still dependent on the threshold of deciding a low curvature point. Furthermore, the junction optimization process could lead to some difficulties in correctly interpreting the junction position as originally produced by craftsmen. However, although this point is valid for some specific domains of exact line-drawing representation, such as vectorization, we are interested in detecting local features that would be useful to addressing the problem of large-scale document indexing and retrieval. In this sense, a low rate of false positives in the final results is not problematic. To promote the evaluation of this work, the source code and demonstration of the symbol localization application are publicly available at this address<fnref n="5" />. Furthermore, several potential directions of research have planned for improvement of this work:

* It can be noted that the junction optimization algorithm is dependent on the detection of junction candidates. Any false detections of the junction candidates could

<footnote n="5">⁵https://sites.google.com/site/ourjunctiondemo/</footnote>
result in the false acceptance of the final junctions. A potential solution for this matter would probably rely on the context information. By incorporating some null hypothesis $H_0$ of the underlying context information, the powerful *a contrario* detection theory [Desolneux et al., 2000] can be applied to justify the meaningfulness of a given detected junction. That means, a junction is meaningful if it is unlikely to occur at random under the hypothesis $H_0$. For instance, a successful application of this theory for detecting the junctions in natural images was presented in [Xia, 2011].

* To this end, the evaluation protocol applied to the junction detectors is detector-dependent. Different evaluation metrics and protocols shall be used to study the behavior of the proposed junction detector. It would be also interesting to evaluate the junction detectors using some datasets provided with semantic groundtruth such as the BSDS benchmark presented in [Maire et al., 2008].

* Junction characterization has been addressed in our work as a natural benefit of the junction optimization process. It has several favourable properties of simplicity, distinctiveness, scale invariant, and low dimensionality. However, this characterization is not shift invariant (i.e., it depends on the selection of the start junction arm). Therefore, further works to address these points would make the whole system completely robust. Besides, the idea of employing an off-the-shelf local descriptor in the CV field would be also interesting to characterize the junctions.

Regarding the second contribution of this dissertation about feature indexing, we first provided, in chapter 4, a deep review of the state-of-the-art methods for feature indexing in high-dimensional feature vector space. The main ideas, favourable features, and shortcomings of each method are thoroughly studied. We also give our subjection remarks for these methods and highlight the need of an advanced contribution for efficiently indexing the feature vectors.

Following the conclusions derived from the discussion of the existing indexing methods, a new contribution for feature indexing in high-dimensional feature vector space has been presented in chapter 5. The proposed algorithm, called linked-node m-ary tree (LM-tree), has many desirable features that make it different from all the existing methods. Followings are three main advancements attributed to the proposed LM-tree.

* A new polar-space-based method for data decomposition has been presented to construct the LM-tree. This new decomposition method differs from the existing methods in that two axes are randomly selected from a few axes that have the highest variance of the underlying dataset to iteratively partition the data.

* An efficient pruning rule is proposed to eliminate the search paths that are unlikely to contain the true answers. Furthermore, the lower bounds are easily computed, as if the data were in 2D space, regardless of how high the dimensionality is.

* A bandwidth search method is introduced to explore the nodes of the LM-tree. Its basic idea is inspired by the fact that searching at multiple adjacent bins gives a good chance of reaching the true answers, while reducing the computational overhead. Like this, it makes unnecessary the expensive computations of finding the best bins that is a matter of the priority search technique [Beis and Lowe, 1997].
The proposed LM-tree has been validated on a wide corpus dataset composing of one mil-
lion SIFT features and 250000 GIST features, demonstrating that the proposed algorithm
gives a significant improvement of search performance, compared to the state-of-the-art
indexing algorithms, including randomized KD-trees, hierarchical K-means tree, random-
ized clustering trees, and multi-probe LHS scheme. To further confirm the efficiency the
proposed indexing algorithm, an additional application to image retrieval was developed
using a wide corpus set of historical books containing ornamental graphics. Performance
evaluation has been reported, showing that the LM-tree indexing algorithm is very time-
efficient. At last, the source code of this contribution is also made available for the interest
of researchers at this address<fnref n="6" />.

Notwithstanding the obtained results are interesting, we realize that many different
lines of researches for this work are still possible. Followings are several improvements
planned in the future work.

* To this end, the proposed LM-tree is a balance but static tree. That said, our indexing algorithm would not work for the cases where the data are dynamically changed over time. Accordingly, dynamic insertion and deletion of the data points in the LM-tree would be also investigated in the future to make the system adaptive to the dynamic change of the data.

* As the study of the binary features is more and more becoming a great interest of researches, different tests on binary features would be interesting to evaluate the proposed LM-tree.

* It is agreed that performance of an indexing algorithm is highly dependent on partic-
ular dataset. It is therefore a good idea to train the system on a particular dataset
provided some prior knowledge about the desired search precision. As already re-
ported in [Muja and Lowe, 2009], automatic parameter tuning using cross-validation
approach would be integrated to make the system more robust and well-fitting to
specific datasets.

* The use of many high variance axes for data partitioning (e.g., more than two) would be considered to study the new behavior of the system.

* In addition to these open works, the study of designing an indexing system, which can work on the data stored on external disk, will be investigated to deal with extremely large datasets that are not able to be fully loaded into the main memory.

As the last conclusion, we expected that the two contributions in this dissertation shall
draw much of interest and attention from the research community.

<footnote n="6">⁶https://sites.google.com/site/ptalmtree/</footnote>
# Liste des publications et prix

## Revues internationales avec comité de lecture

* The Anh Pham, Mathieu Delalandre, Sabine Barrat and Jean-Yves Ramel, "Accurate junction detection and characterization in line-drawing images". Pattern Recognition (2013), In Press. DOI: http://dx.doi.org/10.1016/j.patcog.2013.06.027

* The Anh Pham, Sabine Barrat, Mathieu Delalandre and Jean-Yves Ramel, "An efficient tree structure for indexing feature vectors". Invited submission to a special section of Pattern Recognition Letters (2013).

## Conférences internationales avec comité de lecture et actes

* The Anh Pham, Sabine Barrat, Mathieu Delalandre and Jean-Yves Ramel, "An efficient indexing scheme based on linked-node m-ary tree structure". In proceedings of the 17th International Conference on Image Analysis and Processing (ICIAP 2013), A. Petrosino (Ed.): Part I, LNCS 8156, pp. 752-762, 2013.

* The Anh Pham, Mathieu Delalandre, Sabine Barrat and Jean-Yves Ramel, "Robust symbol localization based on junction features and efficient geometry consistency checking". In proceedings of the 12th International Conference on Document Analysis and Recognition (ICDAR 2013), pp. 1083-1087, 2013.

* The Anh Pham, Mathieu Delalandre, Sabine Barrat and Jean-Yves Ramel, "Accurate Junction Detection and Reconstruction in Line-Drawing Images", In proceedings of the 21st International Conference on Pattern Recognition (ICPR 2012), pp. 693-696, 2012.

* The Anh Pham, Mathieu Delalandre, Sabine Barrat and Jean-Yves Ramel, "A robust approach for local interest point detection in line-drawing images", In proceedings of the 10th IAPR International Workshop on Document Analysis Systems (DAS 2012), pp. 79-84, 2012.

* The Anh Pham, Mathieu Delalandre and Sabine Barrat, "A contour-based method for logo detection", In proceedings of the 11th International Conference on Document Analysis and Recognition (ICDAR 2011), pp. 718-722, 2011.
Conférence nationale avec comité de lecture et actes

* The Anh Pham, Sabine Barrat, Mathieu Delalandre and Jean-Yves Ramel, "Une approche robuste pour la détection de points d'intérêts dans les images de documents techniques", Colloque International Francophone sur l'Écrit et le Document 2012 (CIFED), pp. 445-460, March 21th-23th, 2012, Bordeaux, France.

Prix

* "IAPR Best Student Paper Award" for the paper: The Anh Pham, Sabine Barrat, Mathieu Delalandre and Jean-Yves Ramel, "An efficient indexing scheme based on linked-node m-ary tree structure". In proceedings of the 17th International Conference on Image Analysis and Processing (ICIAP 2013), A. Petrosino (Ed.): Part I, LNCS 8156, pp. 752–762, 2013
Bibliography

[Aslan et al., 2008] Aslan, C., Erdem, A., Erdem, E., and Tari, S. (2008). Disconnected skeleton: Shape at its absolute scale. *IEEE Trans. Pattern Anal. Mach. Intell.*, 30(12):2188-2203.

[Auclair, 2009] Auclair, A. Vincent, N. C. L. (2009). Hash functions for near duplicate image retrieval. In *Workshop on Applications of Computer Vision (WACV2009)*, pages 1-6.

[Awrangjeb and Lu, 2008] Awrangjeb, M. and Lu, G. (2008). An improved curvature scale-space corner detector and a robust corner matching approach for transformed image identification. *IEEE Trans. Image Process.*, 17(12):2425-2441.

[Bai et al., 2007] Bai, X., Latecki, L. J., and Liu, W.-Y. (2007). Skeleton pruning by contour partitioning with discrete curve evolution. *IEEE Trans. Pattern Anal. Mach. Intell.*, 29(3):449-462.

[Ballard, 1981] Ballard, D. (1981). Generalizing the hough transform to detect arbitrary patterns. *Communications of the ACM*, 13(2):111-122.

[Bay et al., 2008] Bay, H., Ess, A., Tuytelaars, T., and Van Gool, L. (2008). Speeded-up robust features (surf). *Comput. Vis. Image Underst.*, 110(3):346-359.

[Beckmann et al., 1990] Beckmann, N., Kriegel, H.-P., Schneider, R., and Seeger, B. (1990). The r*-tree: an efficient and robust access method for points and rectangles. *ACM SIGMOD Record*, 19(2):322-331.

[Beis and Lowe, 1997] Beis, J. S. and Lowe, D. G. (1997). Shape indexing using approximate nearest-neighbour search in high-dimensional spaces. In *Proceedings of the 1997 Conference on Computer Vision and Pattern Recognition*, CVPR'97, pages 1000-1006.

[Belongie et al., 2002] Belongie, S., Malik, J., and Puzicha, J. (2002). Shape matching and object recognition using shape contexts. *IEEE Trans. Pattern Anal. Mach. Intell.*, 24(4):509-522.

[Berchtold et al., 1996] Berchtold, S., Keim, D. A., and Kriegel, H.-P. (1996). The x-tree: An index structure for high-dimensional data. In *Proceedings of the 22nd International Conference on Very Large Data Bases*, VLDB'96, pages 28-39.
[Bergevin and Bubel, 2004] Bergevin, R. and Bubel, A. (2004). Detection and characterization of junctions in a 2d image. *Comput. Vis. Image Underst.*, 93(3):288-309.

[Biederman, 1986] Biederman, I. (1986). Human image understanding: recent research and a theory. In *The second workshop on Human and Machine Vision II*, volume 13, pages 13–57.

[Bodic et al., 2009] Bodic, P. L., Locteau, H., Adam, S., Heroux, P., Lecourtier, Y., and Knippel, A. (2009). Symbol detection using region adjacency graphs and integer linear programming. In *International Conference on Document Analysis and Recognition (ICDAR)*, pages 1320–1324.

[Böhm et al., 2001] Böhm, C., Berchtold, S., and Keim, D. A. (2001). Searching in high-dimensional spaces: Index structures for improving the performance of multimedia databases. *ACM Comput. Surv.*, 33(3):322–373.

[Calonder et al., 2010] Calonder, M., Lepetit, V., Strecha, C., and Fua, P. (2010). Brief: binary robust independent elementary features. In *Proceedings of the 11th European conference on Computer Vision: Part IV*, ECCV’10, pages 778–792.

[Carmona-Poyato et al., 2005] Carmona-Poyato, A., Fernández-García, N. L., Medina-Carnicer, R., and Madrid-Cuevas, F. J. (2005). Dominant point detection: A new proposal. *Image Vision Comput.*, 23(13):1226–1236.

[Chávez et al., 2001] Chávez, E., Navarro, G., Baeza-Yates, R., and Marroquín, J. L. (2001). Searching in metric spaces. *ACM Comput. Surv.*, 33(3):273–321.

[Cheng et al., 1984] Cheng, D.-Y., Gersho, A., Ramamurthi, B., and Shoham, Y. (1984). Fast search algorithms for vector quantization and pattern matching. In *The IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP'84)*, volume 9, pages 372–375.

[Chiang et al., 1998] Chiang, J. Y., Tue, S., and Leu, Y. C. (1998). A new algorithm for line image vectorization. *Pattern Recognition*, 31(10):1541–1549.

[Choi et al., 2009] Choi, S., Kim, T., and Yu, W. (2009). Performance evaluation of ransac family. In *Proceedings of British Machine Vision Conference*, BMVC'09, pages 1–12.

[Damiand et al., 2009] Damiand, G., Higuera, C., Janodet, J.-C., Samuel, E., and Solnon, C. (2009). A polynomial algorithm for submap isomorphism. In *Graph-Based Representations in Pattern Recognition*, volume 5534 of *Lecture Notes in Computer Science*, pages 102–112.

[Delalandre et al., 2010] Delalandre, M., Ramel, J.-Y., Valveny, E., and Luqman, M. (2010). A performance characterization algorithm for symbol localization. In *Graphics Recognition. Achievements, Challenges, and Evolution*, volume 6020 of *LNCS*, pages 260–271.

[Delalandre et al., 2008] Delalandre, M., Valveny, E., and Llados, J. (2008). Performance evaluation of symbol recognition and spotting systems: An overview. In *Workshop on Document Analysis Systems (DAS)*, pages 497–505.
[Deschênes and Ziou, 2000] Deschênes, F. and Ziou, D. (2000). Detection of line junctions and line terminations using curvilinear features. *Pattern Recogn. Lett.*, 21(6-7):637-649.

[Deseilligny et al., 1998] Deseilligny, M. P., Stamon, G., and Suen, C. Y. (1998). Veinerization: A new shape description for flexible skeletonization. *IEEE Trans. Pattern Anal. Mach. Intell.*, 20(5):505-521.

[Desolneux et al., 2000] Desolneux, A., Moisan, L., and Morel, J.-M. (2000). Meaningful alignments. *Int. J. Comput. Vision*, 40(1):7-23.

[di Baja, 1994] di Baja, G. S. (1994). Well-shaped, stable, and reversible skeletons from the (3,4)-distance transform. *J. Vis. Commun. Image R.*, 5(1):107-115.

[Dori and Liu, 1999] Dori, D. and Liu, W. (1999). Sparse pixel vectorization: An algorithm and its performance evaluation. *IEEE Trans. Pattern Anal. Mach. Intell.*, 21(3):202-215.

[Dosch and LLados, 2004] Dosch, P. and LLados, J. (2004). Vectorial signatures for symbol discrimination. In *Workshop on Graphics Recognition (GREC), LNCS*, volume 3088, pages 154-165.

[Douglas and Peucker, 1973] Douglas, D. H. and Peucker, T. K. (1973). Algorithms for the reduction of the number of points required to represent a digitized line or its caricature. *Cartographica: The International Journal for Geographic Information and Geovisualization*, 10(2):112-122.

[Dutta et al., 2013a] Dutta, A., LladóS, J., Bunke, H., and Pal, U. (2013a). Near convex region adjacency graph and approximate neighborhood string matching for symbol spotting in graphical documents. In *Proceedings of 12th International Conference on Document Analysis and Recognition (ICDAR13)*, pages 1078-1082.

[Dutta et al., 2011] Dutta, A., Llados, J., and Pal, U. (2011). Symbol spotting in line drawings through graph paths hashing. In *International Conference on Document Analysis and Recognition (ICDAR)*.

[Dutta et al., 2013b] Dutta, A., LladóS, J., and Pal, U. (2013b). A symbol spotting approach in graphical documents by hashing serialized graphs. *Pattern Recognition*, 46(3):752-768.

[Embley et al., 2011] Embley, D. W., Machado, S., Packer, T., Park, J., Zitzelberger, A., Liddle, S. W., Tate, N., and Lonsdale, D. W. (2011). Enabling search for facts and implied facts in historical documents. In *Proceedings of the 2011 Workshop on Historical Document Imaging and Processing*, HIP '11, pages 59-66.

[Faas and van Vliet, 2007] Faas, F. G. A. and van Vliet, L. J. (2007). Junction detection and multi-orientation analysis using streamlines. In *Proceedings of the 12th international conference on Computer analysis of images and patterns*, CAIP'07, pages 718-725.

[Fan et al., 1998] Fan, K.-C., Chen, D.-F., and Wen, M.-G. (1998). Skeletonization of binary images with nonuniform width via block decomposition and contour vector matching. *Pattern Recognition*, 31(7):823-838.
[Fischler and Bolles, 1981] Fischler, M. A. and Bolles, R. C. (1981). Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography. *Communications of the ACM*, 24(6):381-395.

[Förstner, 1994] Förstner, W. (1994). A framework for low level feature extraction. In *Proceedings of the third European conference on Computer Vision (Vol. II)*, ECCV'94, pages 383-394.

[Förstner and Gülch, 1987] Förstner, W. and Gülch, E. (1987). A fast operator for detection and precise location of distinct points, corners and centres of circular features. In *Intercomission Conference on Fast Proc. of Photogrammetric Data*, pages 281-305.

[Friedman et al., 1977] Friedman, J. H., Bentley, J. L., and Finkel, R. A. (1977). An algorithm for finding best matches in logarithmic expected time. *ACM Trans. Math. Softw.*, 3(3):209-226.

[Fukunaga and Narendra, 1975] Fukunaga, K. and Narendra, M. (1975). A branch and bound algorithm for computing k-nearest neighbors. *IEEE Trans. Comput.*, 24(7):750-753.

[Gionis et al., 1999] Gionis, A., Indyk, P., and Motwani, R. (1999). Similarity search in high dimensions via hashing. In *Proceedings of the 25th International Conference on Very Large DataBases*, VLDB'99, pages 518-529.

[Graves et al., 2009] Graves, A., Liwicki, M., Fernández, S., Bertolami, R., Bunke, H., and Schmidhuber, J. (2009). A novel connectionist system for unconstrained handwriting recognition. *IEEE Trans. Pattern Anal. Mach. Intell.*, 31(5):855-868.

[Guha et al., 1998] Guha, S., Rastogi, R., and Shim, K. (1998). Cure: an efficient clustering algorithm for large databases. In *Proceedings of the 1998 ACM SIGMOD international conference on Management of data*, SIGMOD '98, pages 73-84.

[Guttman, 1984] Guttman, A. (1984). R-trees: a dynamic index structure for spatial searching. *ACM SIGMOD Record*, 14(2):47-57.

[Han and Fan, 1994] Han, C.-C. and Fan, K.-C. (1994). Skeleton generation of engineering drawings via contour matching. *Pattern Recognition*, 27(2):261-275.

[Hansen and Neumann, 2004] Hansen, T. and Neumann, H. (2004). Neural mechanisms for the robust representation of junctions. *Neural Comput.*, 16(5):1013-1037.

[Harris and Stephens., 1988] Harris, C. and Stephens., M. (1988). A combined corner and edge detector. *In Alvey Vision Conference*, pages 147-151.

[Heitger, 1995] Heitger, F. (1995). Feature detection using suppression and enhancement. Technical Report TR-163, Image Science Lab.

[Hilaire and Tombre, 2001] Hilaire, X. and Tombre, K. (2001). Improving the accuracy of skeleton-based vectorization. In *Proceedings of the 4th International Workshop GREC'01, LNCS*, volume 2390, pages 273-288.