An Integrated Fast Hough Transform for Multidimensional Data

Line, plane and hyperplane detection in multidimensional data has many applications in computer vision and artificial intelligence. We propose Integrated Fast Hough Transform (IFHT), a highly-efficient multidimensional Hough transform algorithm based on a new mathematical model. The parameter space of IFHT can be represented with a single k-tree to support hierarchical storage and “coarse-to-fine” search strategy. IFHT essentially changes the least square data-fitting in Li's Fast Hough transform (FHT) to the total least squares data-fitting, in which observational errors across all dimensions are taken into account, thus more practical and more resistant to data noise. It has practically resolved the problem of decreased precision of FHT for target objects mapped to boundaries between accumulators in the parameter space. In addition, it enables a straightforward visualization of the parameter space which not only provides intuitive insight on the number of objects in the data, but also helps with tuning the parameters and combining multiple instances if needed. In all simulated data with different levels of noise and parameters, IFHT surpasses Li's Fast Hough transform in terms of robustness and precision significantly.

deployment of Light Detection and Ranging (LiDAR) equipment in vehicles and satellites [11], [12], Hough transform for 3D point clouds attracts more and more attention in auto-piloting and remote sensing research [13], [14].
The computational cost of Hough transform usually increases dramatically when it is extended to high dimensional data due to the "curse of dimensionality" [15]. Among many variants of Hough transform, Li's Fast Hough transform (FHT) is one of the best that can efficiently deal with high dimensional data [16]. We noticed literatures in the field [17] sometimes refer Fast Hough transform to Brady's radon transform-based algorithm [18]. To avoid confusion, we hereafter refer FHT to Li's algorithm unless stated otherwise. Li's FHT algorithm is highly efficient in terms of computation and storage in high-dimensional data. This is achieved mainly through a "coarse-to-fine" strategy. For n-dimension data, FHT transforms input data points into target objects (lines in 2D, planes in 3D, hyperplanes in higher dimension) in the parameter space and searches for intersection. It subsequently divides the parameter space into n subspaces, each represented by a k-tree. Each subspace is divided recursively into hypercubes from low to high resolution and further analysis is performed only on the hypercubes with votes exceeding a certain threshold. A hypercube receives a vote only when a hyperplane, corresponding to an input data point in the data space, intersects it. This hierarchical approach leads to a significant reduction on computational load and storage. The algorithm has wide applications [19] and has been implemented by some open source image analysis libraries, such as GANDALF [20].
In an n-dimensional dataset, FHT divides the parameter space into n independent subspaces. For example, a straight line in a 2D space can be modelled as y = mx + c, with m → ∞ to describe vertical lines. A direct search of parameter space [m, c] would be difficult if not impossible. FHT models straight lines with two equations y = mx + c, |m| ≤ 1 and x = m y + c , |m | < 1. This effectively divides the search space into two parameter subspaces (m, c) and (m', c'), with |m| ≤ 1 , |m | < 1. Since m, c, m' and c' are all limited, efficient strategies can be designed to exploit all search space. The same model has been used in a plenary of other HT variants such as adaptive Hough transform [21].
Searching in multiple independent parameter subspaces poses several serious challenges. Though in theory any line would not appear in both parameter subspaces, e.g., y = mx + c, |m| ≤ 1 and x = m y + c , |m | < 1 , in practice lines with the slope close to 1 need to be pursuit in both parameter subspaces when noise is considered. Thus, the n parameter subspaces in FHT must overlap marginally with each other to avoid potential gaps. With the increase of the data dimension, significant computational and storage costs are needed for the overlapped regions in the search subspaces.
To make matters worse, FHT operates in n independent parameter subspaces and each parameter subspace uses a unique transformation. As a consequence, testing whether objects detected in different parameter subspaces are from the same object in the original data space is a necessary and very challenging task and needs complicated algorithms. This was indeed noticed in both FHT patent and the original FHT publication and is listed as future work. However, to our best knowledge, there is still a lack of reliable method to resolve it after the initial publication of FHT decades ago.
In the paper, we propose an integrated fast Hough transform (IFHT) parameter space to model lines, planes and hyperplanes. IFHT allows us to represent the parameter space with a single k-tree. Apart from the obvious advantage on computation and storage due to the elimination of redundancy in search space, the model also had two advantages. First, it essentially changes the least square data-fitting in FHT to the total least squares data-fitting, in which observational errors across all dimensions are taken into accounts, thus more resistant to data noise. Furthermore, it has practically resolved the problem of decreased precision of FHT for target objects mapped to boundaries between accumulators in the parameter space. The second advantage stems from the fact that the integrated parameter space allows an intuitive presentation of vote distribution in the parameter space for 2D and 3D data, which provide an intuitive insight on tuning the parameters and combining multiple instances from the same object.

II. METHOD
Hough transform for line, plane or hyperplane detection in general comprised two steps: input data transformation into parameter space and voting for the accumulators. The workflow of HT is illustrated in Fig. 1.

A. The Introduction of FHT
In 2D data, a point (x, y) can lie on any line y = mx + c, |m| ≤ 1 or x = m y + c , |m | < 1. The FHT parameter space comprises two sub-spaces (m, c) and (m , c ). Any point (x, y) can be mapped into the first parameter subspace (m, c) into a 0 + a 1 m + a 2 c = 0 (|m| ≤ 1). The transformation is the transformation to map the point into the second parameter subspace (m , c ) as line a 0 + a 1 m + a 2 c = 0 (|m | < 1) is It can be re-formulated as testing the intersection of a square m ∈ [m 1 , m 2 ], c ∈ [c 1 , c 2 ] and line a 0 + a 1 m + a 2 c = 0 with known values of a 0 , a 1, a 2 . If [m 1 , m 2 ] and [c 1 , c 2 ] have the same length, which can be achieved through data normalization, the test can be relaxed to check whether the line intersects with the circumscribing circle (or hypersphere in high dimension), resulting in dramatically reduction in the computational cost. For example, we only need to test whether |a 0 + a 1 m * + a 2 c * | < r, where m * and c * are the coordinates of the square's center and r is the radius of the circle (Fig. 2). This approximation could lead to over-counting since a line may intersect the circumscribing circle but not the square. Interestingly, experiments indicate it can significantly reduce the computational complexity in high dimensional data and has trivial effect with the precision of the algorithm.
In general, FHT has the following workflow. Let the input be a set of n-dimensional feature points in data space The framework of FHT in a high dimensional data space can be formulated as follows: 1) Constructing a k-dimensional parameter space denoted by {X 1 , X 2 , …, X k } with n subspaces and n transformation functions that maps each F(j) into a hyperplane in the parameter space which results in that any hyperplane is represented by where a i (j) is a function of F(j) and is normalized such 2) Defining a threshold for number of vote "T" needed to justify a hyperplane and the desired finest resolution "q". 3) Dividing the parameter space recursively into hypercubes from low to high resolution, representing with n k-trees, and performing the further analysis only on the hypercubes with votes exceeding a selected threshold.

B. The Integrated Mathematical Model of Hyperplanes by IFHT
IFHT models hyperplanes in the data space in a very intuitive way. Let the input be a set of n-dimensional feature points in data space Here τ i is a scale factor for normalization purpose, and the way to select it will be discussed later. For simplicity, we here relax the restriction of the model by dropping the requirement β 1 ≥ 0. Suppose we know the initial ranges of In this new parameter space, each dimension X i is in [-1, 1] and the parameter space can be quantized into hypercubes for voting. In addition, the requirement of β 1 ≥ 0 can be easily satisfied by dropping specific hypercubes with X 1 < 0 as X 1 = β 1 L 1 . Eq. (2) is the full model of IFHT, which allows targeted searching, i.e., searching objects within a specific part of the space. For the general application aiming for all objects in the space, (2) can be simplified into Here β n+1 is limited and can be quantified into small intervals for investigation. It can be easily proved that the β n+1 is smaller than or equal to the maximum distance D of data points {F 1 , F 2 , …, F n } to the origin in the space coordinate system. For simplicity, we can rewrite the above equation as where d equals to the maximum distance D of data points {F 1 , F 2 , …, F n } to the origin in the space coordinate system. In practice, d could be a bit smaller, for example 0.8D, for the better resolution. Here has the same length along each dimension and can be quantized into cubes for voting.
In the above explanation, the dimension of the parameter space is n+1 for the n-dimensional input data. This is not necessarily so. For example, if we are interested in target objects (lines, planes or hyperplanes) with zero intercept, i.e., objects that pass through the origin, the IFHT parameter space will be n-dimensional. We hereafter denote k as the dimension of the IFHT parameter space for generality.

C. Hierarchical Vote Counting Scheme
Every point in the data space corresponds to one hyperplane in the IFHT parameter space. To achieve a robust detection, small perturbation in the parameter should be allowed to deal with the noise inside the data. This can be easily implemented in IFHT by evaluating a small interval along each dimension in the parameter space rather than a specific value. The length of interval then corresponds to the required resolution of the target hyperplane (straight line for 2D space). To test whether a hyperplane consisted of more than T points exists in the data space, IFHT quantified the parameter space into small accumulators and let each data point to vote for them. Since the IFHT parameter space has the same length along each dimension, the parameter space can be quantized into small accumulators (squares, cubes, or hypercubes) with the same length along each dimension too.
IFHT parameter space is linear. If an accumulator collects T votes, no any sub-accumulator can surpass T votes. This allows the IFHT parameter space to be quantized hierarchically. We can divide IFHT parameter space into several large accumulators for voting and then divide only those with enough votes into smaller accumulators for further testing until the required resolution is achieved. This coarse-to-fine strategy pose great advantage compared to other polar-coordinate based Hough transform methods.
As the accumulator has the same length along each dimension, we can use the same approximation method shown in Fig. 2 for vote testing. The test equation is where [C 1 , C 2 , …, C k ] is the center of accumulator and r is the radius of the hypersphere circumscribing the accumulator (hypercube).

D. A Single K-Tree Representation
In IFHT, the parameter space is a hollow hyper-semi-cylinder (Fig. 3). We represent the parameter space with a nested hierarchy hypercube and associate a k-tree with it. The root node of the tree corresponds to a hypercube centered at a vector C 0 with side length S 0 . Each node of the tree has 2 k children arising when that node's hypercube is halved along each of its k dimensions.
where each b i is −1 or 1 is used as an index for a child of a node. If a node at level l has center C l , then the center of its child with index where S l+1 is the length of child at level l+1 and S l+1 = S l /2. IFHT borrows the coarse-to-fine mechanism from FHT. In detail, for a hypercube at level l centered at [C l1 , C l2 , …, C lk ], the distance D l to a hyperplane defined in (3) can be estimated as D l = k i=1 a i (j)C li . By normalizing distance with sidelength S l , we have We can iteratively calculate this normalized distance. Suppose we known the initial ranges of {X 1 , X 2 , …, X k ) in (3) are centered The test criterion in 4) can be simplified as In (2), we use constraint k−1 i=1 (τ i · β i ) 2 = 1, k = n+1 to guarantee that the hyperplane in the data space is unique. When neither prior information nor preference about where the target hyperplane exists, the initial ranges of In the parameter space, this is equivalent to test whether an accumulator intersect with (k-1)-dimensional unit hypersphere. Following the same principle of the simplification used for (4), for an accumulator centered at [C 1 , C 2 , …, C k ] at level l, with r the radius of the hypersphere circumscribing the accumulator (hypercube), the test equation is

E. Workflow for IFHT
In data space {F 1 , F 2 , …, F n }, when giving a set of data point F(j) = [F 1 (j), F 2 (j), …, F n (j)], IFHT parameters that need to be predetermined are as follows: 1) The minimum votes count "T" as threshold and the desired finest resolution "q". 2) A transformation that maps each F(j) into the parameter space {X 1 , X 2 , …, X k } which results in that any hyperplane is represented by where a i (j) is a function of F(j) and is normalized so that k i=1 a 2 i (j) = 1. For a typical case where no preference or prior information exists, the target hyperplanes can be modelled as where d is a scale smaller than or equal to the maximum distance of data points to the origin. The mapping function for (10) is k = n + 1, a i (j) = F i (j) W (j)d , j = 1, . . . , n, and a n+1 (j) = 1

W (j)
where W (j) is a weighted scale used to insure n+1 i=1 a 2 i (j) = 1. IFHT can be summarized as follows: 1) compute the initial normalized distance from all hyperplane to the root node and count its votes. The root node is usually centered at (0, …, 0) with half of the side-length as (1, …, 1). 2) for each child node, calculate whether it intersects with k-dimensional unit hypersphere with (9). If not, stop processing that child node. If yes, calculate the normalized distance for that child node with (6) and vote for that child node with (7). If it has vote count larger than T, subdivide it into k-tree child nodes. For the child nodes of the root node, only those with first index b 1 = 1 are kept for further analysis. 3) When there is no node with the resolution equal to q and the vote count larger than T, record the node with the highest resolution. Otherwise, output all nodes with the resolution equal to q and the vote count larger than T.

F. The Complexity of IFHT
In FHT, it has been proved that the number of active accumulators does not increase exponentially with the quantization level and the dimensionality of the parameter space. This is fully guaranteed by the "thin tree" theorem. In IFHT, the "thin-tree" theorem is as follows: Thin tree Theorem: For a hyperplane with M data points, each data point is transformed with (2) into a hyperplane in the parameter space and M hyperplanes in the parameter space intersect at a single point P. Assume they are uniformly distributed in orientation at the point P. Given a minimum vote threshold T, the number of active accumulators of size q that receives T or more votes is less than a constant number K that does not depend on q.
Proof: For IFHT, if we drop the constraint n i=1 (τ i · β i ) 2 = 1 and β 1 ≥ 0, the transform can be regarded as a special instance of FHT in (n+1)-dimension. If M hyperplanes (corresponding to the M data points in the data space) in the parameter space intersect at point P with a uniformly distribution in orientation, it satisfies FHT "thin tree" theorem. As the constraint only decreases the number of active accumulators to be tested, so the number of active accumulators of size q that receive T or more votes is less than a constant number K that does not depend on q.
In Fig. 5, we collected the overall number of active accumulators by IFHT for two 3D datasets. It is clear that in both examples, the number of active accumulators does not increase exponentially with the quantization level.

G. Two Possible Variants of IFHT
In IFHT described above, a hyperplane is modelled as It is possible to modify the equation slightly in certain scenarios. One modification could be to put the intercept β n+1 of the line into the normalization part. This would relieve us from specifying the range parameter for β n+1 which often requires calculating the maximum distance of data points to the origin. The corresponding mathematical model is as follows: In addition, the implementation of IFHT algorithm could be modified slightly. In the IFHT mathematic model, n+1 variables have been used. However, when {β 2 , …, β n+1 } is determined, β 1 can be uniquely decided and serves as a dummy variable. It is possible to calculate β 1 straightforwardly without putting it into the parameter space in the IFHT implementation.

III. RESULTS
In all variants of Hough transform, to determine the desired the quantization level q is not a trivial task. When a very coarse quantization level with a small value of q is chosen, too many lines, planes or hyperplanes could be reported which cannot satisfy the required precision. When a large value of q is chosen, the final results could be very sensitive to statistical noise or quantization noise in the data. In FHT, it is very frequently observed the same lines, planes or hyperplanes are reported as multiple instances by several accumulators at a specific level l, but no lines or hyperplanes can be found at level l+1. As a consequence, to automatically report how many lines or hyperplanes inside the high-dimensional data with HT or FHT is still challenging. In IFHT, the integrated parameter space allows an intuitive presentation of vote distribution in the parameter space, which provides an intuitive insight on combining detected lines or tuning the parameters as shown in Fig. 7 for 2D data and Fig. 8 for 3D data. It is very clear that when a coarse resolution is chosen, the multiple instances from the same line or hyperplane always cluster together. This observation is a manifest of IFHT thin-tree theorem in the parameter space. It provides a very robust and intuitive way to determine the number of objects to be reported, and it is also instructive in parameter tuning and merging of multiple instances if necessary.
The proposed IFHT method has some interesting properties. In FHT, the noise in different dimension is not treated equally. For example, in 2D data, y = mx + c, |m| ≤ 1 and x = m y + c , |m | < 1 are used for modeling hyperplanes. When fitting the model to available data points, FHT minimizing the error parallel to the y axis, as shown in Fig. 8(a). IFHT treats each dimension of data equally, and essentially changes the least square data-fitting in FHT to the total least-square data-fitting [22], in which observational errors in both dependent and independent variables were taken into account. In IFHT, we test whether it is possible that β 1 x + β 2 y + β 3 = 0 where β 2 1 + β 2 2 = 1 for a specific pair of β 1 , β 2 . This is achieved through calculating the direct distance a square to a line in 2D case. The data-fitting can be illustrated in Fig. 8(b). We compared the performances of our IFHT method and HFT with simulated data using Monte Carlo method. In the first experiment, a specific number N of data points from a line in 2D space were generated, with x ∼ U (−2, 2), y = 0.8x + 0.24 + g.
Here g is Gaussian white noise with standard deviation 0.04 and 0.08. Then the same number of random data points with x, y ∼ N (0, 1), without those in the proximity of the line, were generated for the background data points. For each set of parameters, 100 images were generated. IHFT and HFT with the same parameters (T = 0.75N and l = 5) were used to detect the line in each image. The number of correctly reported points on the detected line were recorded and the ratio to the overall number of true points of the line is reported as detection rate. Fig. 9(a) shows that in both scenarios, the performance of IHFT is significantly better than that of FHT. It has been noted that when the initial parameters have been chosen, the positioning and size of accumulators in FHT will be fixed. In some cases, peaks may straddle the boundaries between accumulators [21]. Due to the different parameter spaces have been used, only very few lines cause this type of problem to both FHT and IFHT methods. Line y = 0.5 is one of them when parameters (l = 5 and d = 2) are used in both methods. Fig. 9(b) shows the detection rate of both methods when the standard deviation of the Gaussian white noise is 0.04 and 0.08, respectively. The performance of IFHT was not significantly affected, while the performance of FHT was indeed seriously undermined in this scenario.
In above simulations, the white noise has only been added into the second dimension of the observed data. This would be rare in the real data since the noise usually distributed across both dimensions of observed data. We repeated the first experiment with noise presented across both dimensions. The observed data point for the line is t ∼ U (−2, 2), x = t + g1, y = 0.8t + 0.24 + g2. Here g1 and g2 are Gaussian white noise with the same standard deviation. We tested IFHT and FHT using with the same parameters (T = N/2 and l = 5). Their performances are shown in Fig. 10. Apparently, IFHT has significantly better performances than FHT.
We further tested the performance of IFHT for the real 3D data. An RGB-D image from a publicly available image database [23] was chosen. The RGB image and the corresponding depth image are shown in Fig. 11(a) and (b) respectively. The edge   of objects was detected using the Canny operator implemented in OpenCV Fig. 11(c), and then transformed into 3D data points using the transformation suggested by Kinect. The input was centered and scaled with standard deviation along each dimension (by R function "scale"). IFHT was then applied with parameters (T = 600, l = 7 and the range of intercept = 4). IFHT successfully detected the plane containing the edge of the table and the inner surface of the sofa as shown in Fig. 11(d).

IV. CONCLUSION
The detection of line, plane or hyperplane in multidimensional data has many applications in computer vision, data clustering and artificial intelligence. Though there are many 2D algorithms, few of them are practical for high-dimensional data due to the increased computational and storage demands. The Fast Hough transform approach suggested by Li et al. is so far one of the most efficient with relative space efficiency R s = (α/β) k for a k dimensional case. However, the mathematical model used by FHT imposes a restriction that multiple independent parameter subspaces need to be constructed and searched. The burden of dealing with multiple parameter subspaces not only increases the complexity of algorithm design and implementation, but also makes the intuitive visualization of the parameter space almost impossible. As a consequence, the parameter tuning of FHT is difficult and there is lack of an efficient method to merge multiple instances from the same object (line, plane or hyperplane). In addition, the sensitivity and precision of FHT vary with the parameters of target objects. Lines, planes and hyperplanes which are mapped to the boundaries of accumulators in the parameter space could have very low precision. These disadvantages have seriously limited the application of FHT algorithm.
We have proposed an integrated mathematical model for Hough transform, corresponding to an integrated fast Hough transform parameter space. The parameter space of IFHT can be represented with a single k-tree to support hierarchical searching. Apart from the obvious advantage on simplicity of algorithm design and reduces redundancy in computation and storage, the model allows the intuitive presentation of the parameter space. This new insight can be used to decide the number of lines, planes or hyperplanes in the data, and to help with tuning the parameters and combining detected instances if necessary. Furthermore, IFHT essentially changes the least square data-fitting in FHT to the total least squares data-fitting, in which observational errors across all dimensions are taken into accounts, thus more resistant to data noise. In our test, IFHT successfully avoided the drawback of varied sensitivity and precision embed in FHT algorithm, with very high precision for almost all lines or planes in the data space. In almost all simulated images with different levels of noise and parameters, IFHT surpassed FHT in terms of robustness and precision significantly. IFHT also shows good performance for the real 3D data.