An iterative approach to Hough transform without re-voting

The Hough transform is a well known pattern recognition technique to detect features in images when specified in some parametric form [5]. The classical Hough transform is most commonly used for the detection of simple curves such as straight lines, circles, ellipses [5, 3]. However, in recent papers its effectiveness has been proved in the case of more complex classes of curves, specifically, families of irreducible algebraic plane curves like elliptic curves, conchoids of Slüse, oblique bifolium, curves with 3 convexities, Wassenaar curves, quartic curves with a tacnode, and piriform curves [1, 6]. Some of these curves have been successfully applied for identifying profiles of interest on CT slices of the human skeleton, including lumbar vertebrae and the spinal canal at different heights of the vertebral column [6].

The basic idea to perform this identification is that points on a curve in the image space are mapped into curves (more generally hypersurfaces) in a parameter space which intersect at a uniquely determined point. The coordinates of this intersection point identify the curve in the image space. The procedure that computationally realizes this identification is based on a discretization of the parameter space and on the definition of an accumulator matrix counting for the number of points mapped into curves passing through (i.e., voting) its cells. The set of parameters corresponding to the maximum in the accumulator allows the recognition of the desired curve in the image space.

From a mathematical viewpoint, the use of the Hough transform to recognize curves of more complexity in an image requires to prove the duality type result that associates a unique point in the parameter space to many different points on a curve in the image space [1]. From a computational viewpoint, the number of parameters in the curve equation defines the dimension of the parameter space. Thus, the computational effort of the voting procedure accordingly increases as the dimension of the accumulator matrix. Further, more practically, many recognition tasks in real applications require the identification of profiles that cannot be associated to a single curve but are more realistically represented by piecewise defined curves. For example, the identification of bone profiles in the human skeleton requires a method to assign an appropriate piece of curve to a specific portion of the real skeleton profile. Therefore, an approximation of the whole district shape can be obtained by using an iterative approach to Hough transform in order to iteratively approximate different parts of the same bone profile with different curves.

In order to realize such a procedure, one can eliminate those points in the profile that have been recognized by using the Hough transform according to a specific family of curves, and re-applying the method by searching for a different curve from the family using the remaining set of points. However, this procedure is computationally expensive and therefore in this paper we introduced a layered accumulator matrix that allowed us to discard the contribution of some points in the voting procedure without re-voting. This way, we reached the twofold objective of iteratively approximating the entire profile with a piecewise defined curve without dramatically increasing the computational cost of the algorithm.

The paper is organized as follows. In Section 2 we recall some theoretical results [1] on the Hough transform of algebraic curves. We also prove that, for the family of curves in the form used in the presented application, the requested conditions for the uniqueness of the intersection point in the parameter space are verified. In Section 3 we define the layered accumulator and we present the recognition algorithm according to this new accumulator definition. In Section 4 we present some numerical results considering the recognition of a lumbar vertebra in a X-ray CT image. Finally, in Section 5 our conclusions are offered.

This section surveys some definitions and basic facts introduced in [1] and [6], to which refer for more details and examples. Although these results hold over an infinite integral ring $K$, having in mind applications to imaging, we assume $K=\bbbr$. We introduce a family of irreducible polynomials

where the coefficients $g_{ij}(\lambda)$ are polynomials in the independent parameters $\lambda=(\lambda_{1},\ldots,\lambda_{t})$ varying in an euclidean open set $\mathcal{U}\subseteq\bbbr^{t}$, and with the degree, $d$, of $f_{\lambda}(x,y)$ not depending on $\lambda$.

Let ${\mathcal{F}}$ be the corresponding family of zero loci ${\mathcal{C}}_{\lambda}$ of $f_{\lambda}(x,y)$ and assume that each ${\mathcal{C}}_{\lambda}$ is a curve in the affine plane ${\mathbb{A}}_{(x,y)}^{2}(\bbbr)$. So we want a family $\mathcal{F}=\{{\mathcal{C}}_{\lambda}\}$ of irreducible curves which share the degree. If $P=(x_{P},y_{P})$ is a point in the image space ${\mathbb{A}}_{(x,y)}^{2}(\bbbr)$, the polynomial

defines a hypersurface in the parameter space ${\mathbb{A}}_{(\Lambda_{1},\ldots,\Lambda_{t})}^{t}(\bbbr)$, where $\Lambda_{1},\ldots,\Lambda_{t}$ are affine indeterminates. We say that this hypersurface is the Hough transform $\Gamma_{P}(\mathcal{F})$ of the point $P$ with respect to the family $\mathcal{F}$.

In [1, Theorem 2.2] it is proved that, given a family ${\cal{F}}$ of irreducible curves defined as above, the Hough transforms $\Gamma_{P}(\mathcal{F})$, when $P$ varies on a curve ${\mathcal{C}}_{\lambda}$, all pass through the same point $\lambda$ and, if there exists another point $\lambda^{\prime}=(\lambda_{1}^{\prime},\ldots,\lambda_{t}^{\prime})$ where the Hough transforms $\Gamma_{P}(\mathcal{F})$ intersect, then ${\mathcal{C}}_{\lambda}={\mathcal{C}}_{\lambda^{\prime}}$. Moreover, the following two conditions are equivalent.

1. 1.

   For any pair of curves ${\mathcal{C}}_{\lambda}$ and ${\mathcal{C}}_{\lambda^{\prime}}$ in ${\mathcal{F}}$, the equality ${\mathcal{C}}_{\lambda}={\mathcal{C}}_{\lambda^{\prime}}$ implies $\lambda=\lambda^{\prime}$.

2. 2.

   For each curve ${\mathcal{C}}_{\lambda}$ in ${\mathcal{F}}$, one has $\cap_{P\in{\mathcal{C}}_{\lambda}}\Gamma_{P}(\mathcal{F})=\lambda$.

A family ${\cal{F}}$ of curves which meets one of the above equivalent conditions is said Hough regular. The curves of those families share the same property of the straight lines as in the original definition of the Hough transform (see [3]). Therefore Hough regularity has a crucial computational role in the solution of the pattern recognition problem since for a Hough regular family of curves it is possible to formulate a computationally effective generalization of the algorithm for line detection, able to automatically extract in the image a profile described by the curves from the family. Clearly, the above equivalent properties are not true in general. For example, for the family of conics of equation $a^{2}x^{2}+b^{2}y-1=0$ we have that ${\mathcal{C}}_{a,b}={\mathcal{C}}_{-a,-b}$. It is therefore of special interest to find Hough regular families of curves. To this aim, we observe that, if

are the equations of ${\mathcal{C}}_{\lambda}$ and ${\mathcal{C}}_{\lambda^{\prime}}$ respectively, then ${\mathcal{C}}_{\lambda}={\mathcal{C}}_{\lambda^{\prime}}$ reads $f_{\lambda}(x,y)=kf_{\lambda^{\prime}}(x,y)$ for some $k\in\bbbr^{*}:=\bbbr-\{0\}$, that is,

for each pair of indices $i$, $j$ (see [4, Theorem 1, p. 365]). Relations (3) can be used to check the Hough regularity of the family, i.e., that $\mathcal{C}_{\lambda}=\mathcal{C}_{\lambda^{\prime}}$ implies $\lambda=\lambda^{\prime}$.

In the following example we describe the family of curves we deal with throughout the paper.

Let $P_{j}$, $j=1,\ldots,\nu$, be the set of points of interest in the image space, and let $\Gamma_{P_{j}}(\mathcal{F})$ be their Hough transforms in the parameter space ${\mathbb{A}}_{(\Lambda_{1},\ldots,\Lambda_{t})}^{t}(\bbbr)$. The first step of the standard recognition algorithm is based on a discretization of a subset $\mathcal{T}$ of the parameter space, i.e., on the definition of a set of cells

where $\lambda=(\lambda_{1},\lambda_{2},\ldots,\lambda_{t})$ is a point in $\mathcal{T}$, $\lambda_{\bf n}=(\lambda_{1,n_{1}},\lambda_{2,n_{2}},\ldots,\lambda_{t,n_{t}})$ is a sampling point, ${\bf n}=(n_{1},n_{2},\ldots,n_{t})$ is a multi-index identifying both a specific cell and the corresponding sampling point, $d_{k}$ and $N_{k}$ are the sampling distance and the number of samples considered along the $k$-component, respectively. Following the standard approach, the accumulator matrix $H=(H_{\bf n})$ is then defined as

and it is usually computed by means of the iterative process

with initialization $H^{(0)}_{\bf n}=0$ and final result $H_{\bf n}=H^{(\nu)}_{\bf n}$, for each ${\bf n}$. The search for the maximum of $H$ at position ${\bf n^{*}}=(n_{1}^{*},n_{2}^{*},\ldots,n_{t}^{*})$ allows us to identify the set of parameters $\lambda^{*}=\lambda_{{\bf n}^{*}}=(\lambda_{1,n_{1}^{*}},\lambda_{2,n_{2}^{*}},\ldots,\lambda_{t,n_{t}^{*}})$ and finally the curve ${\mathcal{C}}_{\lambda^{*}}$ best approximating the set of points $P_{j}$ in the image.

This approach can be effectively applied to the recognition of significant profiles in medical and astronomical images [1, 6]. However, in practical applications, a better accuracy can be achieved by approximating the profile of interest using more than one curve, each one fitting a specific portion of the profile.

A possible approach is based on the search of successive maxima in the accumulator matrix, starting from the absolute one and moving on in descending order, up to a fixed threshold. This scheme has the main drawback that local maxima often cluster in the same region of the parameter space, thus identifying curves characterized by very similar sets of parameters, well approximating the same subsets of points in the same portion of the profile.

We propose here an alternative approach based on the principle one man one vote, i.e., each point in the dataset is allowed voting for just one curve. Once a curve has been recognized, all the points voting for it have to be removed from the dataset and the computation of the accumulator matrix has to be repeated without their contribution. In order to realize this idea without increasing the computational effort associated to the computation of the process defined by (8), we introduce a layered accumulator matrix $\mathcal{A}=(\mathcal{A}_{{\bf n},j})$, where the $j$-th layer is a binary matrix containing the votes of the $j$-th point $P_{j}$ of the original dataset to the whole set of cells. That is,

By definition $\mathcal{A}_{{\bf n},j}\in\{0,1\}$ is then set to $1$ if and only if the Hough transform $\Gamma_{P_{j}}(\mathcal{F})$ of the point $P_{j}$ passes across the cell ${\bf C}_{\bf n}$. The usual accumulator $H$ can be easily computed from $\mathcal{A}$ summing over all the layers:

Then a first curve ${\mathcal{C}}_{\lambda^{*}}$ can be recognized in the usual way. However, the main advantage of this re-formulation is that one can easily identify the points which voted for the cell ${\bf C}_{{\bf n}^{*}}$ by looking for the non-zero entries $\mathcal{A}_{{\bf n}^{*},j}$, and then remove the corresponding layers from $\mathcal{A}$. Calling $\mathcal{S}$ the whole set of values the index $j$ can assume in the original dataset, and $\mathcal{G}\subset\mathcal{S}$ the subset containing the indices of the removed layers, we can re-compute the usual accumulator $H$ without the contribution of points $P_{j}$, $j\in\mathcal{G}$, and without re-voting, by simply computing

A new curve can now be recognized through optimization of the updated accumulator matrix, and so on. The iterative approach inspired by previous considerations is summarized in Algorithm 1. A pictorial representation of the layered accumulator $\mathcal{A}$ is shown in Fig. 1 (where we refer to the case $t=2$, so that $\mathcal{A}$ is a $3$-dimensional matrix).

We applied the piecewise Hough transform algorithm described in the previous section to recognize a lumbar vertebra profile by using the family $\mathcal{F}$ of elliptic curves of equation (4). In particular, we have considered an axial X-ray CT slice of a human lumbar vertebra (see Fig. 2(a)) and we have applied a standard Canny [2] edge detection algorithm in order to select the set of edge points shown in Fig. 2(b). From this initial dataset we then extracted the edges corresponding to the outer profile of the vertebra by applying an ad hoc $\alpha$-concave hull algorithm. This algorithm looks for the subset of points that uniquely defines a polygon enclosing the whole dataset, without having any angle between two neighboring points exceeding $\pi+\alpha$. In our application we set $\alpha=\pi/4$ and we extracted the subset $\mathcal{S}$ of $\nu=959$ blue points, as seen in Fig. 2(d). Each point $P_{j}\in\mathcal{S}$ is then Hough transformed in a quartic hypersurface of equation (5), in the parameter space ${\mathbb{A}}_{(A,B,M,N)}^{4}(\bbbr)$. This parameter space is then discretized in order to have a $4$-dimensional accumulator matrix $H$ made of $N_{H}=324135$ cells. As a consequence, the layered accumulator $\mathcal{A}$ results in a $5$-dimensional matrix with $\nu$ layers of $N_{H}$ entries each. The iterative scheme described in the previous section (with threshold_value to stop the iterations set to $35\%$ of $\nu$) gives us a set of $7$ different elliptic curves from the fixed family, best approximating different portions of the outer profile of the vertebra. The whole set of curves is superimposed to the profile edge points in Fig. 2(d), together with green points highlighting the boundaries of each different portion. As for the inner profile, the same process gives us a set of $9$ different curves from the family $\mathcal{F}$. The piecewise defined curves providing the best approximations for both internal and external vertebra profiles, are then shown in Fig. 2(e) and Fig. 2(c) on the whole set of edge points and on the original X-ray CT image, respectively.

We describe an iterative approach of the Hough transform technique to provide an efficient procedure for approximating the profile of a bone in the human skeleton as a collection of different but continuously attached curves. The core of the proposed procedure is the definition of a layered accumulator matrix. From its knowledge, it is possible to compute the usual accumulator matrix at each iteration and to recognize and remove from the dataset the points lying on portions of the profile, already well approximated. More importantly, the layered accumulator allows us to iterate the algorithm on an input dataset progressively smaller, without the need of repeating the computationally demanding task of voting. We tested this iterative approach to identify bone profiles in a clinical X-ray tomography image, showing that it accurately recognizes the inner and outer profiles of a human lumbar vertebra. Clinical applications of this recognition method may include either the characterization of the compact bone asset for the prognosis of chronic leukemia or the estimate of the bone marrow asset from integrated CT and nuclear medicine data [7].