SPReAD: On Spherical Part Recognition
by Axial Discretization in 4D Hough Space
Radhika Mittal and Partha Bhowmick
Department of Computer Science and Engineering
Indian Institute of Technology, Kharagpur, India
{radhikamittal.iitkgp@gmail.com, bhowmick@gmail.com}
Abstract. A novel algorithm is proposed to locate the sets of adjacent co-spherical triangles for a given object, which enables us to detect
spheres and spherical parts constituting the object. An extension of the
idea of Hough transform has been used, aided by axial discretization
and restricted searching, along with the geometric data structure of doubly connected edge list. The algorithm has been analyzed and shown
to achieve significant efficiency in space and run-time. On testing the
algorithm with various 3D objects, it is found to produce the desired
result. Effects of different input parameters have been explained and the
robustness of the algorithm has been shown for rough/noisy surfaces.
1
Introduction
Recognition and analysis of 3D shapes in general, and spherical components in
particular, have gained a significant research perspective for several important
applications, e.g., detection of pulmonary nodules [19], medical analysis of brains
[15], etc. Among different geometric approaches for 3D shape analysis [9,11], the
most commonly used is Generalized Hough Transform (GHT) [1, 6], as it can
easily be adapted based on specific needs of the application, e.g., [3, 7, 16–18].
We propose here an extension of Hough transform to detect spherical components, using a more space- and time-efficient technique of restricted discretization
and localized searching instead of processing a 4D array containing the 3D coordinates of the center and the radius. To work with a significantly large range
with a high accuracy, we use the novel idea of axial discretization. We take discrete points along the circum-normal (axis) of a face (triangle) and use restricted
Hough transform to locate the parameters of the sphere it potentially belongs to,
with respect to the given object. The restricted HT is based on the observation
that if a face belongs to a particular sphere in the given object, then with a
very high probability, its adjacent faces would also belong to the same sphere.
The geometric data structure of DCEL has been used extensively to achieve an
attractive run-time, as shown in this paper.
2
Proposed Algorithm
One of the novel principles while designing the algorithm SPReAD lies in the
concept of axial discretization. Such a discretization has been resorted to in
2
R. Mittal and P. Bhowmick
order to reduce the computational burden of O(n4 ) for the conventional Hough
transform (HT) [10, 14]. This technique has not been used so far for recognizing
spherical parts in 3D objects, although proposed to detect circular arcs in 2D
images for several applications [4, 5, 8, 12].
2.1
Principle of SPReAD
We use the following observation: If a (triangular) face f belongs to a particular spherical part in the given 3D object, then with a very high probability,
the three faces adjacent to f would also be co-spheric with f . Based on this
observation, we estimate the parameters, i.e., center and radius of the sphere,
which all the adjacent faces potentially belong to. The algorithm SPReAD considers the discrete points along the circum-normals of these faces as input, and
finds the approximate center. (The circum-normal of a face f is the normal to
f that passes through its circum-center.) The first face f is selected at random.
Although there exists a unique circle passing through the three vertices of any
triangle, there exists an infinitude of spheres passing through them. The centers
of all these spheres are collinear in 3D space. Hence, for a given face f with
a center point taken on its circum-normal, the sphere parameters are first obtained, and then we build the initial spherical part on it by considering its three
adjacent triangles to obtain the common solution space of these four triangles.
The larger spherical part is built from this by doing a breadth-first search iteratively on the set of faces adjacent to each other whose corresponding ε-balls in
the parameter space have all pair-wise intersection, as explained next.
The parameter space has four dimensions, the first three being corresponding
to the three coordinates of the object space and the fourth to the radius. The
distance metric used by us in this parameter space is based on L2 -norm. If pi
be the point in the 4D parameter space corresponding to a sphere Si passing
through the vertices of a face f , then the ε-ball corresponding to Si is denoted
by φ(pi , ε) and defined as the 4D hyper-sphere with center at pi and radius ε.
Note that, if Si has center ci = (xi , yi , zi ) and radius ri , then it is uniquely
represented in the parameter space by the point pi = (ci , ri ) = (xi , yi , zi , ri ).
Hence, two faces fi (with circum-normal Ni ) and fj (with circum-normal Nj )
are said to be co-spheric with each other if and only if there exist two points
pi (ai , ri ) and pj (aj , rj ) in the parameter space such that ai ∈ Ni , aj ∈ Nj , and
there exists a non-empty intersection between φ(pi , ε) and φ(pj , ε). Evidently,
two ε-balls are said to intersect each other if and only if the distance between
their centers does not exceed 2ε. In general, a set of faces are mutually co-spheric
if the two ε-balls corresponding to some axial point of one and some axial point
of the other in each pair from the set have a non-empty intersection in the 4D
parameter space.
The proposed algorithm SPReAD is made flexible to suit the needs of different objects by providing the hyper-ball radius, ε. The value of ε is specified
w.r.t. the normalized scale, and varies from 0.01 to 0.04, as shown in Sec. 3.
Smaller the value of ε, tighter is the detection of spherical parts, and larger the
value of ε, looser is the detection. The speed and precision of the algorithm is
SPReAD: Spherical Part Recognition
3
decided by another parameter, namely N , which indicates the number of discrete
points to be considered along the circum-normal of each triangle while detecting
a spherical part. And lastly, the minimum number of constituent faces of a cospheric part, namely F , is considered to decide whether the concerned part can
be reported as a practically spherical part. The exact usage of these parameters
is described in detail later in this section and also in Sec. 3.
To maintain uniformity across different sizes of objects, we apply an isotropic
scaling on the object by normalizing the coordinates of all the vertices of the
triangulated object surface such that the maximum over the distances for all
pairs of vertices is a constant (taken as 10 in our experimentation). We maintain
a Doubly Connected Edge List (DCEL) [2] which helps in finding the adjacent
faces. In this list, we store, for each edge, its source vertex, destination vertex,
incident face, and twin edge. The twin twin[e] of an edge e is the edge whose
incident face is adjacent to the face incident on e. Note that, for each face, each
of its incident (directed) edges is defined in a manner so that the face lies left of
the edge. Thus, the source vertex of e coincides with the destination vertex of
twin[e], and vice versa (Fig. 1).
For each triangle, we find its circum-center coordinates (xc , yc , zc ) as follows.
xc =
2
∑
i=0
))
( (
2
2
2
xi li2 l(i+1)mod
−
l
+
l
i
(i+2)mod 3
3
2
2 (
∑
i=0
) ∑
2
2
li2 l(i+1)mod
li4
3 −
i=0
where li (i = 0, 1, 2) denotes the length of the side opposite to the ith vertex
having coordinates (xi , yi , zi ). The corresponding formulas for finding yc or zc
are similar, and obtained simply by replacing xi with yi or zi , as applicable. The
unit circum-normal of each triangle is found by normalizing the cross product
of two of its edges, and then by translating it to the circum-center.
We start the process of detecting co-spherical sets of faces with all faces as
initially ‘unmarked’. Then a face f is picked at random. We find the three faces
adjacent to f using the DCEL in constant time (Fig. 1). For each such adjacent
face fi , we uniformly discretize a finite segment of the circum-normal N passing
through its circum-center. The discretization procedure samples 2N + 1 points
on each circum-normal: N points on either side of the plane of the face and one
lying on it. They are sampled in a manner such that the distance between the
farthest pair of sample points along the concerned circum-normal is a constant
(taken to be 100 times the size of the normalized object, in our experimentation).
If there exists a sample point a ∈ N for f and a sample point ai ∈ Ni for a face
fi with p = (a, r) and pi = (ai , ri ), where r is the distance of the point a from a
vertex of f (all vertices of f being equidistant from a) and ri is that of ai from
a vertex of fi , such that φ(p, ε) ∩ φ(pi , ε) ̸= ∅, then there exists a solution for
the common sphere of f and fi (see Fig. 1). The parameters corresponding to
this solution are obtained from φ(p, ε) ∩ φ(pi , ε).
We maintain a queue Q to store all the faces that are co-spheric with f . Thus,
all the adjacent faces which have been taken into consideration are enqueued in
4
R. Mittal and P. Bhowmick
Fig. 1. Left: Information in DCEL. The red triangle denotes a face f having three
(directed) edges, shown in blue. Each of of three edges is defined by a pair of vertices
(blue), and has a twin edge in its adjacent face (orange). The unit normals are computed
from the edges corresponding to the faces. Right: Finding the center of the sphere built
by a face f (shown in red) and its three adjacent faces {f1 , f2 , f3 } (orange), using axial
discretization of their circum-normals (N, N1 , N2 , N3 ).
Q along with the face f ; these faces are marked on before being enqueued. We
calculate the mean of all satisfying points on the discretized circum-normals to
estimate the effective center and radius of the sphere S built by the face f and
its surrounding faces. If the estimated radius is larger than a sufficiently large
constant (taken to be twice the normalized size of the object in our experimentation), then Q is emptied, indicating that no sphere is formed by f . This is
done because a very large radius implies a planar segment instead of a spherical
patch.
We build a spherical part S using breadth-first-search. For each face in Q,
we find the adjacent faces using the DCEL. If an adjacent face fi , thus found, is
not marked on, then we check whether there exists a point ai in the discretized
circum-normal of fi such that φ(pi , ε) has a non-empty intersection with the
existing parameter space corresponding to S. If it has a non-empty intersection,
then the face fi is marked on and enqueued in Q. We continue expansion on
each face in Q, until Q stops growing or all faces in Q have been considered.
Speeding up the Algorithm To speed up the procedure of searching the
matching points of discretization, if any, while verifying the intersection of one
ε-ball with the solution space of the current sphere S, we use the range of
indices of the discrete points that correspond to the current solution. Note that,
the index varies from −N to N , 0 inclusive. Let [smin , smax ] be the range of
indices that have produced the current solution. Here, smin is the lowest index
of one or more discrete points of some of the circum-normals corresponding to
S, and smax is the highest. If a new face fi belongs to S, then it is expected
SPReAD: Spherical Part Recognition
5
that the index of the discrete point ai from the circum-normal of fi would lie
in or close to the range [smin , smax ]. In fact, if the face fi is comparable in
size with the ‘marked faces’ corresponding to S, then the index of ai would
lie in [smin , smax ]; otherwise, ai would lie outside but close to [smin , smax ]. The
length of [smin , smax ] is also quite small if the sizes of the ‘marked faces’ have a
small variation. Hence, the search for the matching discrete point for the new
face fi is guided by [smin , smax ]. The first point picked up from the discretized
circum-normal of fi has index s0 = 12 ⌊smin + smax + 1⌋. If it matches, then no
more search is necessary; otherwise the indices of two subsequent points under
consideration are s0 ±1, next are s0 ±2, and so on, until the match is found or the
face fi is rejected for a diverging match with the current solution. A diverging
match is reported when for two points in succession, the ε-ball is found to be
digressing away from the current solution.
Post-processing The co-spherical components are colored by the algorithm
SPReAD using distinct colors. The coloring is constrained by the condition
that each spherical part must have a minimum number of faces. This is to avoid
detection of small undesired spherical components and also to avoid detection
of cylindrical or similar shapes whose faces may be misinterpreted to be forming
a spherical part. Hence, an important restriction is that a spherical component
will be colored only if its area is greater than F · Amax , where F is the tuning parameter lying in (0, 1) and Amax denotes the area of the largest spherical
component found. Post-processing is done by SPReAD using the color information in order to fill up the gaps created by the undetected faces in the detected
spherical parts. If a set A of adjacent faces is uncolored and all the faces adjacent
to the boundary faces of this set have the same color (i.e., belong to the same
sphere, S), then we check whether the area formed by the set A is less than a
threshold (taken as 0.1 times the area of the corresponding sphere). If so, the
set A is considered to be ‘falsely undetected’, an hence the color of each face of
A is set to the color of the sphere S.
2.2
Complexity Analysis
Let the number of vertices be v, the number of faces be n, and the number of
edges be e. Then the DCEL takes O(v + n + e) = O(n) space, since v = O(n)
and e = O(n). Apart from the conventional information as mentioned earlier,
the DCEL also stores the details of faces, vertices, face areas, normal equations,
and circum-centers, taking O(n) space. The queue Q stores all the faces that are
currently co-spherical. The best case occurs when there is no overlap between two
co-spherical sets, resulting in a space complexity of O(n). In the worst case, there
can be many overlaps among these co-spherical sets. An instance of the queue
may contain as many as O(n) faces and there can be O(n) such instances of the
queue. However, when a particular instance of Q is processed, other instances
do not occupy any space, since they either have been processed earlier or would
be processed later. This gives us a space complexity of O(n) for Q.
6
R. Mittal and P. Bhowmick
It takes O(v+n) time to read the object, O(v 2 ) time to normalize the vertices
by computing the maximum distance, and O(n2 ) time to prepare the DCEL.
The algorithm checks no more faces adjacent to a face fki when fki is found
to be not matching with the spherical part Sk currently under consideration.
For each spherical part, the first two faces are obtained in O(N ) time by index
searching. Let the number of faces of Sk be nk . All these nk faces can have
O(nk ) adjacent faces which are not part of Sk ; and each of these faces will be
tested, resulting to no positive result. Hence, the time complexity to report Sk
is given by the sum of time complexity to report the faces of Sk and that to
check the mismatching faces. The best-case time complexity is, therefore, given
by O(N )+(nk −2)O(1)+O(nk )O(1) = O(nk +N ); summing up over all spherical
parts and observing that nk is bounded by O(n) for this sum, we get O(n+mN ),
where m is the number of spherical parts reported by SPReAD. For the worst
case, we sum up O(N ) + (nk − 2)O(N ) + O(nk )O(N ) = O(nk N ) over all parts,
hence getting O(nN ).
3
Experiment and Results
We have implemented the algorithm SPReAD in C using OpenGL on an Intel(R) Core(TM)2 Duo CPU E4500 2.20 GHz machine, the OS being Ubuntu
Release 10.10. Table 1 shows its results on some objects having one or more
spherical parts. The object Sphere is completely covered with ε = 0.03; note
that 759 out of 840 triangles are initially detected, and the remaining 81 triangles get included during post-processing (Sec. 2.1).
Table 1. Summary of results for running SPReAD on different objects (ε − 0.03, F =
0.10, N = 1000). n = number of faces; m = number of spheres detected; nk = number(s)
of faces for each sphere (shown comma-separated).
Parameter
n
m
nk
CPU time (sec.)
Sphere
840
1
759
0.596
Globe
7618
2
3968, 927
2.610
Balloon
2080
2
1384, 134
0.759
Elk
6518
4
513, 492, 483, 476
3.542
We take N in the range of 1000–2000; for a low value, the points on the
normal are largely separated from each other, thereby giving a bad estimate
(Fig. 2). The specified value of ε determines the accuracy with which spheres
SPReAD: Spherical Part Recognition
7
(a) (0.04, 0.10, 200) (b) (0.04, 0.08, 2000) (c) (0.06, 0.08, 2000) (d) (0.04, 0.3, 1000).
Fig. 2. Output instances showing effect of user paramester (ε, F, N ): (a) too small
value of N ; (b, c) change in ε; (d) large value of F . (For proper results, see Table 1.)
are detected (Fig. 2). Changing ε by even 0.02 units (in normalized object scale
of unity) significantly affects the number of co-spherical triangles. A high value of
ε would encompass unwanted triangles in a spherical component, whereas a low
value would discard many desired triangles, thereby hindering the formation of
a larger spherical part. For F , a high value may not correctly report a spherical
part (Fig. 2); when F = 0.3 (see correct result in Table 1), base part of the
globe is not reported as spherical as its area is appreciably less than the red
one.
Post-processing takes care of missing triangles in a spherical part, possibly due to staggering circum-normals. Figure 3(a) shows the result of running
SPReAD on elk, without post-processing. For testing SPReAD against noisy
or rough surface, we have prepared some rough surfaces by introducing Gaussian noise. Figure 3(b) shows that on introducing Gaussian noise with a peak
of 1% of the maximum distance (0.05), fairly good results are obtained without
post-processing. However, for a Gaussian noise of 2%, we get better results on
doing post-processing (Fig. 3(c, d)).
Fig. 3. Role of post-processing: (a) Result on elk without post-processing (see Table 1
for result after post-processing. (b-d) Results for rough surfaces, with and without 2nd
stage processing. (b) 1% noise, no post-processing; (c) 2% noise, no post-processing;
(d) 2% noise, after post-processing (N = 1000, ε = 0.06, F = 0.1 for all three cases).
8
R. Mittal and P. Bhowmick
4
Conclusion
The proposed algorithm SPReAD can detect the adjacent co-spherical triangles
for a given object using axial discretization based on the contemporary idea of
restricted Hough transform, resulting in usage of less space as compared to the
generalized Hough transform. We have tested the algorithm on various 3D objects, have analyzed the effects of different input parameters, and have tested its
robustness rough or noisy surfaces. All these results exhibit that the algorithm
is quite resistant to noise and can generate fairly accurate results with a spaceand time-efficient manner. As a future work, the algorithm can be extended to
detect other 3D shapes like ellipsoids and paraboloids by properly applying the
concept of axial discretization.
Acknowledgement
A part of this work has been sponsored by the NSF project: NEBULA (CNS-1038695).
References
1. Ballard, D.H.: Readings in computer vision: Issues, problems, principles, and
paradigms, 1987.
2. Berg, M.D. et al.: Computational Geometry Algorithms and Applications, 2000.
3. Cao, M.Y. et al.: Spherical parameter detection based on hierarchical Hough
transform. PRL, 27:980–986, 2006.
4. Chiu, S.H. and Liaw, J.J.: An effective voting method for circle detection. PRL,
26(2):121–133, 2005.
5. Davies, E.R.: A high speed algorithm for circular object detection. PRL, 6:323–
333, 1987.
6. Duda, R.O. and Hart, P.E.: Use of the Hough transformation to detect lines and
curves in pictures. Commun. ACM, 15:11–15, 1972.
7. Ecabert, O. and Thiran, J.P.: Adaptive Hough transform for the detection of
natural shapes under weak affine transformations. PRL, 25(12):1411–1419, 2004.
8. Foresti, G.L. et al.: Circular arc extraction by direct clustering in a 3D Hough
parameter space. Signal Processing, 41:203–224, 1995.
9. Gelfand, N. and Guibas, L.J.: Shape segmentation using local slippage analysis.
In Proc. SGP ’04, pages 214–223, 2004.
10. Gonzalez, R.C. and Woods, R.E.: Digital Image Processing. Addison-Wesley,
California, 1993.
11. Hofer, M.: et al. 3d shape recognition and reconstruction based on line element
geometry. In Proc. ICCV ’05, pages 1532–1538, 2005.
12. Ioannoua, D. et al.: Circle recognition through a 2D Hough transform and radius
histogramming. IVC, 17:15–26, 1999.
13. Kim, H.S. and Kim, J.H.: A two-step circle detection algorithm from the intersecting chords. PRL, 22(6-7):787–798, 2001.
14. Leavers, V.: Survey: Which Hough transform? CVGIP, 58:250–264, 1993.
15. Nain, D. et al.: Statistical shape analysis of brain structures using spherical
wavelets. In Proc. ISBI ’07, pages 209–212, 2007.
16. Ogundana, O.O. et al.: Fast Hough transform for automated detection of spheres
in three-dimensional point clouds. Opt. Eng. 46, 2007.
SPReAD: Spherical Part Recognition
9
17. Rong, F. et al.: A novel Hough transform algorithm for multi-objective detection.
In Proc. IITA ’09, pages 705–708, 2009.
18. Yang, M.C.K. et al.: Hough transform modified by line connectivity and line
thickness. IEEE Trans. PAMI, 19:905–910, 1997.
19. Zhang, X. et al.: A new method for spherical object detection and its application to computer aided detection of pulmonary nodules in CT images. In Proc.
MICCAI ’07, pages 842–849, 2007.