PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Int J Rob Res. Author manuscript; available in PMC 2010 July 2.
Published in final edited form as:
Int J Rob Res. 2010 June 1; 29(7): 801–812.
doi:  10.1177/0278364909352700
PMCID: PMC2896220
NIHMSID: NIHMS152034

Generating Uniform Incremental Grids on SO(3) Using the Hopf Fibration

Anna Yershova
Duke University, Durham, NC 27707, USA, ude.ekud.sc@avohsrey
Swati Jain
Duke University, Durham, NC 27707, USA, ude.ekud@niaj.itaws
Steven M. LaValle
University of Illinois, Urbana, IL 61801 USA, ude.cuiu.sc@ellaval

Abstract

The problem of generating uniform deterministic samples over the rotation group, SO(3), is fundamental to computational biology, chemistry, physics, and numerous branches of computer science. We present the best-known method to date for constructing incremental, deterministic grids on SO(3); it provides: 1) the lowest metric distortion for grid neighbor edges, 2) optimal dispersion-reduction with each additional sample, 3) explicit neighborhood structure, and 4) equivolumetric partition of SO(3) by the grid cells. We also demonstrate the use of the sequence on motion planning problems.

1 Introduction

Numerical computations on continuous spaces often require generation of a representative set of samples. The performance of various methods in engineering and scientific fields, such as numerical optimization and integration as well as collisionfree path generation in robot motion planing, rely heavily on the quality of the sampling technique. Hence, it is important that the underlying samples are as good as possible.

A particular problem of discretization of SO(3), the space of 3D rotations, arises in applications, such as biological protein docking problems, robot motion planning, aerospace trajectory design, and quantum computations. Typical operations on this space include numerical optimization, searching, integration, sampling, and path generation. Multiresolution grids are widely used for many of these operations on other spaces which are nicely behaved, such as rectangular subsets of R2 or R3.

It would be wonderful to achieve the same for SO(3); however, the space of 3D rotations is substantially more complicated. In its basic form, SO(3) is defined as a set of matrices that satisfy orthogonality and orientation constraints. It is an implicitly defined, three-dimensional surface embedded in R9. One approach is to place a coordinate system on the surface, causing it to behave like a patch in R3. However, many of such coordinates cause metric distortions in comparison to distances on the original surface. Only few representations of SO(3), such as quaternions and Hopf coordinates, preserve distances and volumes. They treat SO(3) as a unit sphere S3 [subset or is implied by] R4 with antipodal points identified. The volumes of surface patches on S3 correspond to the unique Haar measure for SO(3), which is the only way to obtain distortion-free notions of distance and volume. This implies that if we want to make multiresolution grids on SO(3), we are faced with warping them onto S3. It may seem that such curvature prohibits the introduction of distortion-free grids, similar to the problem of making distance-preserving maps of the world (e.g., Greenland usually looks too big on a flat map). In addition, the identification of antipodal points causes a minor complication in which only half of S3 is used, with unusual connectivity in the equatorial three-plane. However, in this paper we use intrinsic properties unique to SO(3) (first described in [17]) to build almost distortion-free grids and avoid the issue of having to identify the antipodal points on S3.

Due to widespread interest in discretizing SO(3) in numerous fields, there have been considerable efforts in the past. The problem of generating point sets on spheres minimizing various criteria, such as energy functions, discrepancy, dispersion, and mutual distances, has been extensively studied in mathematics and statistics [8, 15, 22, 25, 28, 29]. Random sampling methods were also developed in [2, 24, 27, 34]. Problems of sampling rotational groups and spheres have been studied and applied in the context of computational structural biology, physics, chemistry, computer graphics and robotics [4, 7, 16, 19, 21, 26, 30, 31, 32].

In this paper, we introduce the best-known deterministic method to date for SO(3) in terms of providing:

  1. incremental generation,
  2. optimal dispersion-reduction with each additional sample,
  3. explicit neighborhood structure,
  4. the lowest metric distortion for grid neighbor edges,
  5. equivolumetric partition of SO(3) into grid regions.

The rest of the paper is organized around the presentation of the method. Section 2 defines the topological properties of SO(3) together with its representations that are crucial for presenting our method. Section 3 overviews sampling requirements for the sequence. We discuss the relevant sampling methods that influenced our work in Section 4. Finally, we present our method in Section 5; experimental results and its application to motion planning problems in Section 6. We conclude our work in Section 7.

2 Properties and Representations of SO(3)

The special orthogonal group, SO(3), arises from rotations around the origin in R3. Each rotation, by definition, is a linear transformation that preserves the length of vectors and orientation of space. The elements of SO(3) form a group, with the group action being the composition of rotations. SO(3) is not only a group, but also a manifold, which makes it a Lie group.

To sample SO(3) uniformly, it is necessary to understand its topology. Any method known to date that produces uniform rotations relies on topology and Haar measure of SO(3) (see Section 4).

Topology of SO(3)

SO(3) is diffeomorphic to the real projective space, RP3. It is hard to visualize the real projective space, because it cannot be embedded in R3. Fortunately, it can be represented as RP3=S3(x~x), the more familiar 3-sphere, S3, embedded in R4, with antipodal points identified. Topologists say that the 3-sphere is a double cover of RP3, since one point of the projective space has two corresponding points on the 3-sphere.

Haar Measure on SO(3)

Up to a scalar multiple, there exists a unique measure on SO(3) that is invariant with respect to the group action. This is called the Haar measure. That is, the Haar measure of a set is equal to the Haar measure of all of the rotations of the set. In our particular situation, we can think of the Haar measure as being invariant under all orthogonal coordinate changes. The Haar measure is an intrinsic property of SO(3) which comes from the group structure, and is independent of its topological structure.

We have not used any coordinate system or parametrization of SO(3) yet, since the notion of Haar measure is abstracted from representations of SO(3). One has to use extreme caution when expressing the measure in terms of any of the representations we describe next. Not all of these naturally preserve the Haar measure.

Orthogonal Matrices

The elements of SO(3) are defined as 3 × 3 orthogonal matrices with determinant +1. The group operation is multiplication of matrices. Because rotation matrices are less efficient and less numerically stable than quaternions, they are generally used less often than quaternions.

Quaternions

One of the most useful representations of the projective space is the set of quaternions. Let x=(x1,x2,x3,x4)R4 be a unit quaternion, x1+x2i+x3j+x4k,||x||=1, representing a 3D rotation. Because of the topological relationship between the projective space and the 3-sphere, once the identifications of the antipodal points on the 3-sphere are taken into account, metrics similar to those defined for the 3-sphere can be used for the projective space. Moreover, such metrics will respect the Haar measure on SO(3).

The most natural way to define a metric for any two points x,y [set membership] SO(3) is as the length of the shortest arc between x and y on the 3-sphere, which quaternions conveniently allow to do:

ρSO(3)(x,y)=cos1(xy),
(1)

in which (x · y) denotes the dot product for vectors in R4, and the absolute value, |.|, guarantees that the shortest arc is chosen among the identifications of the two quaternions [11].

Quaternion representation is also convenient for calculating the composition of rotations, which is expressed as multiplication of quaternions. Any rotation invariant surface measure on S3 naturally preserves the Haar measure for SO(3) and can be used for quaternions. However, the surface measure is not straightforwardly expressed using quaternions. Other representations, such as spherical or Hopf coordinates, are more convenient for measuring the volume of surface regions.

Spherical Coordinates for SO(3)

Because of the topological relationship between the 3-sphere and SO(3), hyperspherical coordinates can be used for SO(3). Consider a point (θ, ϕ, ψ) [set membership] S3, in which ψ [set membership] [0,π/2] (to compensate for the identifications, we consider only one hemisphere of S3), θ [set membership] [0,π], and ϕ [set membership] [0,2π). For each ψ, the full ranges of θ and ϕ define a 2-sphere of radius sin(ψ). The quaternion x = (x1,x2,x3,x4) corresponding to the rotation (θ, ϕ, ψ) can be obtained using the formula:

x1=cosΨx2=sinΨcosθx3=sinΨsinθcosϕx4=sinΨsinθsinϕ.
(2)

The volume element on SO(3) defines the Haar measure and has the following expression in spherical coordinates:

dV=sin2ΨsinθdθdϕdΨ.
(3)

This representation is not as convenient for integration as the Hopf coordinates, which have a simpler expression for the Jacobian. Spherical coordinates are also cumbersome for computing compositions of rotations.

Hopf Coordinates for SO(3)

As opposed to spherical coordinates for hyperspheres, the Hopf coordinates are unique for both SO(3) and the 3-sphere. They naturally describe the intrinsic structure of both of these spaces and provide a natural tool for obtaining uniform distributions on these spaces.

The Hopf fibration describes SO(3) in terms of the circle S1 and the ordinary 2-sphere S2. Intuitively, SO(3) is composed of non-intersecting fibers, such that each fiber is a circle S1 corresponding to a point on the 2-sphere. This fiber bundle structure is denoted as SO(3)S1~S2. The Hopf fibration has the important property of locally being a Cartesian product space. The space SO(3), however, is not (globally) the Cartesian product of S2 and S1. Intuitively, SO(3) is the product of S2 and S1 similarly to the way the Möbius band is locally the Cartesian product of an interval and a circle S1. That is, locally, a sequence of coordinates from each subspace results in a global parametrization of the space, whereas the global embedding into the Euclidean space introduces a twist, and does not have the Cartesian product structure. The Hopf coordinates can also be used for the 3-sphere, because of the topological relationship between the 3-sphere and SO(3).

Each rotation in Hopf coordinates can be written as (θ, ϕ, ψ), in which ψ [set membership] [0,2π) parametrizes the circle S1, and θ [set membership] [0,π] and ϕ [set membership] [0,2π) represent spherical coordinates on S2. The transformation to a quaternion x = (x1,x2,x3,x4) can be expressed using the formula:

x1=cosθ2cosΨ2x2=cosθ2sinΨ2x3=sinθ2cos(ϕ+Ψ2)x4=sinθ2sin(ϕ+Ψ2).
(4)

A detailed derivation of the Hopf Coordinates is shown in Appendix 8.1. Briefly, Eq. (4) represents each rotation from SO(3) as a rotation by angle ψ [set membership] S1 around the z axis, followed by the rotation, which places z in a position (θ, ϕ) [set membership] S2. Eq. (4) is obtained after the composition of these two rotations. The Hopf coordinates define exactly half of S3, since the coordinate x2 never takes negative values. The Hopf coordinates can be extended to the entire S3 by increasing the range of ψ to be [0, 4π).

The volume element on SO(3), which is also the surface volume element on S3, can be computed from Eq. (4) (see Appendix 8.2 for a detailed derivation), and has the following form:

dV=18sinθdθdϕdΨ.
(5)

Note that sin θ dθ dϕ represents the surface area on the 2-sphere, and dψ is the length element on S1. This formula additionally demonstrates that the length of a portion of S1 is multiplied by the surface area of the base space, S2, to obtain the volume on SO(3). The coefficient 1/8 results from the fact that neither fibers S1 nor the base space S2 are unit. In fact, in Appendix 8.3 we compute the lengths of the fibers and the surface area of the base space, which is used later for determining the grid cell sizes for our sequence.

As we have shown, the Hopf coordinates preserve the fiber structure of SO(3) and are convenient for integration on SO(3). However, composition of rotations is best expressed using quaternions.

Angle-Axis Representation for SO(3)

One of the most intuitive ways to represent rotations is by using Euler's theorem, which states that every 3D rotation is a rotation by some angle θ around a unit axis n = (n1,n2,n3), ||n|| = 1. The transformation from the angle-axis representation to quaternions is achieved by:

x=(cosθ2,sinθ2n1,sinθ2n2,sinθ2n3).
(6)

The angle-axis representation is useful for visualizing the projective space in R3. Each rotation is drawn as a vector with direction n and a magnitude corresponding to θ (a multiple or a function of θ can be used; see Section 5.6, and [3]). Figure 1 shows the visualization of the spherical and Hopf coordinates on SO(3) using the angle-axis representation. From this visualization one can immediately notice the singularities introduced by the spherical coordinates. It is also possible to see the advantage of using Hopf coordinates from this visualization. Hopf coordinates do not introduce singularities. The circles represented by the range of the variable ψ are all of equal length (see App. 8.3) and non-intersecting; they uniformly cover SO(3). The fiber structure formed by these circles is also seen in the figure.

Fig. 1
Visualization of the spherical and Hopf coordinates on SO(3) using angle and axis representation. This representation corresponds to a projection of the S3 onto the equatorial solid sphere which we draw in R3. (a) The full range of the spherical coordinate ...

Euler Angles Representation

Euler angles are often used in robotics to represent rotations. Each rotation is then a vector (x1,x2,x3),xi [set membership] [−π,π]/−π ~ π. The topology of the resulting space is S1 × S1 × S1, and, therefore, Euler angles do not correctly capture the structure of SO(3). There are many detrimental consequences of this. Special tricks (see [11]) are needed to implement metric and measure that preserve Haar measure. Moreover, Euler angles are hard to compose, and present problems of singularities and the gimbal lock [23].

Understanding and respecting the global topology of SO(3) is crucial for performing other numerical computations on the space. For example, sampling and interpolation of SO(3) using Euler angles in [33] led to failure in producing motion planning path on a relatively simple problem, which was solved in seconds using the correct parametrizations. In the rest of the paper we use Hopf coordinates and quaternions to represent rotations.

3 Sampling Terminology and Problem Formulation

In applications such as motion planning, the algorithms are often terminated early and the particular order in which samples are chosen becomes crucial. Sampling literature distinguishes between a sample set and a sample sequence. For a sample set, the number of points, n, is specified in advance, and a set of n points is then chosen to satisfy the requirements of the method. The notion of ordering between points is not defined for a sample set but becomes important for sequences. Successive points in a sequence should be chosen carefully so that the resulting sample sets are all of good quality. Sequences are particularly suitable for motion planning algorithms, in which the number of points needed to solve the problem is not known in advance.

Now that the background definitions for SO(3) have been presented in Section 2, to generate samples over SO(3) we need to formulate the desirable properties for the samples. The first requirement is that samples form a sequence. We also require that samples get arbitrarily close to every point in SO(3), i.e. that the sequence of samples is dense in SO(3). Next, we formulate several requirements on the uniformity properties of samples.

3.1 Discrepancy and Dispersion

Additional requirements that the sequence needs to satisfy are described by the uniformity measures, discrepancy and dispersion.

Intuitively, discrepancy can be thought of as enforcing two criteria: first, that no region of the space is left uncovered; and second, that no region is left too full. Dispersion eliminates the second criterion, requiring only the first. It can be shown that low discrepancy implies low dispersion [18].

To define discrepancy, choose a range space, R, as a collection of subsets of SO(3). Let R [set membership] R denote one such subset. Range spaces that are usually considered on spheres are the set of spherical caps (intersections of the 3-sphere with half-spaces) or the set of spherical slices (intersections of two 3-hemispheres) [20], which can be used on SO(3) once the identifications of the 3-sphere are taken into account.

Let μ (R) denote the Haar measure of the subset R. If the samples in the set P are uniform in some ideal sense, then it seems reasonable that the fraction of these samples that lie in any subset R should be roughly μ (R) divided by μ (SO(3)) (which is simply π2, see App. 8.3). We define the discrepancy [18] to measure how far from ideal the sample set P is:

D(P,𝓡)=supR𝓡PRNμ(R)μ(SO(3)),
(7)

in which | · | applied to a finite set denotes its cardinality. Figure 2 (a) demonstrates the notion on the 2-sphere.

Fig. 2
An illustration of the notions of dispersion and discrepancy for a set of points on a 2-sphere. (a) The discrepancy searches for the subset R for which the deviation from the measure of R to the number of samples placed inside R is the largest. (b) The ...

While discrepancy is based on measure, a metric-based criterion, dispersion, can be introduced:

δ(P,ρ)=maxqSO(3)minpP(q,p).
(8)

Above, ρ denotes any metric on SO(3) that agrees with the Haar measure, such as (1). Intuitively, this corresponds to the spherical radius of the largest empty ball that fits in between the samples (assuming all ball centers lie on SO(3)). See Figure 2(b) for an illustration.

3.2 Problem Formulation

In summary, the goal of this paper is to define a sequence of elements from SO(3) which:

  • is incremental,
  • is deterministic,
  • minimizes the dispersion (8) and discrepancy (7) on SO(3),
  • has grid structure with respect to the metric (1) on SO(3).

4 Sampling Methods Overview

Our work was influenced by many successful sampling methods developed recently for spheres and SO(3). As demonstrated in Table 1, several of them are highly related to the problem formulated in Section 3.2. However, none of the methods known to date has all of the desired properties.

Table 1
The comparison of different sampling methods related to the problem of Section 3.2. The rows correspond to the desired properties of these methods.

Random Sequences of Elements from Sd

To generate uniformly distributed random points on a hypersphere Sd, spherical symmetry of the multidimensional Gaussian density function can be exploited [6]. For each of the i = 0…d +1 coordinates use a zero-mean Gaussian distribution with the same variance to generate xi. This is done approximately by generating k uniformly distributed values from the interval [−1,1] and adding them following the Central Limit Theorem. In practice, any k ≥ 12 is a reasonable choice. Then the normalized vector (xi/||xi||) is uniformly distributed over the hypersphere Sd.

Random Sequence of SO(3) Rotations

There are several ways of sampling the space of rotations uniformly at random [2, 24, 27, 34]. The main difficulty in doing so is the choice of a convenient parametrization of SO(3). If a parameter space is sampled uniformly, the resulting samples on SO(3) are not necessarily uniform. As was shown in Section 2, not all of the parametrizations of SO(3) are natural representations of rotations, and some of them lead to measure distortions, and even singularities. Only few parametrizations, such as the Hopf coordinates, result in a local isometry to SO(3).

It is easy to make the mistake of sampling rotations using a wrong parametrization [1]. The subgroup algorithm [5] for selecting random elements for SO(3) is the correct and most popular method for uniform random sampling of SO(3). It uses the fact that any Lie group can be uniformly sampled, by combining elements from a subgroup (in case of SO(3) it is S1), and the quotient, or coset space (S2) at random. Essentially, this method utilizes the Hopf coordinates. Random sequences of rotations are used in many applications. However, they lack deterministic uniformity guarantees, and the explicit neighborhood structure.

Successive Orthogonal Images on SO(n)

Related to the subgroup method for generating random rotations, is the deterministic method of Successive Orthogonal Images [17], which generates lattice-like sets with a specified length step based on uniform deterministic samples from the subgroup, S1, and the coset space, S2. The method utilized Hopf coordinates, and is also generalized to arbitrary SO(n).

The deterministic point sets from [17] can be applied to the problems in which the number of the desired samples is specified in advance. If the set of samples on S2 is chosen so that it has a grid structure, the resulting set of samples on SO(3) has the explicit neighborhood structure. Part of the current work will be in applying this method in a way that provides the incremental quality necessary for our motion planning applications.

Layered Sukharev Grid Sequence for Spheres and SO(3)

Uniform, deterministic sequences were first designed for the unit cube [0,1]d [13]. To minimize dispersion, the method places one resolution of grid at a time inside of the unit cube. A discrepancy-optimal ordering is then generated for each resolution. The sequence can be extended to spheres and SO(3) [33] using the projection from faces of an inscribed cube. For SO(3), though, the distortions produced by the method result in some grid cells being four times the volume of others.

The general method for designing Layered Sukharev Grid sequences inside Cartesian products was later presented in [14]. Our current paper builds on top of these works by combining the method in [14], with the Successive Orthogonal Images [17] generation of rotations, Hopf coordinates, and the HEALPix spherical sampling method [7] described next.

HEALPix Multiresolution Grids on S2

The HEALPix package [7] was designed for efficient and incremental discretization of full-sky maps in application to the satellite missions to measure the cosmic microwave background in astrophysics. It provides a deterministic, uniform, and multiresolution sampling method for the 2-sphere. Moreover, it possesses additional qualities, such as equal area partitioning of the 2-sphere, and isolatitude sampling on the 2-sphere, which make computations of the spherical harmonics integrals even more efficient.

The method takes advantage of the measure preserving cylindrical projection of the 2-sphere. This intrinsic property of the 2-sphere cannot be generalized directly to higher dimensional spheres. However, this work shows that an extremely uniform grid can be constructed on such a non-trivial curvature space as the 2-sphere. It is also not difficult to make this grid incremental using the method from [14]. We have done this as part of our implementation for the current work, and the code can be found at [9].

5 Our Approach

In this section we present our approach of sampling SO(3) which satisfies all of the requirements of Section 3.2, and Table 1. The fiber bundle structure of SO(3) locally behaves similarly to the Cartesian product of two spaces, S1 and S2. Therefore, the method presented in [14] for constructing multiresolution grid sequences for Cartesian products of spaces, can be used for constructing a grid sequence on SO(3). The resulting rotations are computed using the Hopf coordinates, as was first described in [17]. It is a much simpler problem to construct nicely behaved grids on the 1-sphere and 2-sphere. Hopf coordinates allow the two grids to be lifted to the space of rotations without loss of uniformity. Next, we outline the details of this construction.

5.1 Description of the Grid Structure on SO(3)

Let ψ be the angle parametrizing the circle, S1, and (θ, ϕ) be the spherical coordinates parametrizing the sphere, S2. Using these coordinates, define T1 to be the multiresolution grid over the circle and T2 to be the multiresolution grid over the 2-sphere. Let m1 and m2 be the number of points at the base resolution 0 of the grids T1 and T2 respectively.

There are numerous grids that can be defined on S2 (see Figure 3 for an illustration of some). In this work we have selected the HEALPix grid [7] on S2, and the ordinary grid for S1. Both of these grids are uniform, have simple neighborhood structure, and can have multiple resolutions. Moreover, HEALPix divides the surface of the 2-sphere into subregions of equal area. After multiplying these by equal length fibers of S1, this results in equivolumetric partition of SO(3) into the grid regions.

Fig. 3
Different sampling methods on S2. (a) 200 random samples (b) 192 Sukharev grid samples [33] (c) icosahedron samples (d) 216 HEALPix samples [7]

Next, consider the space S2~S1. The multiresolution grid sequence that we define for SO(3) has m1 ·m2 · 23l points at the resolution level l, in which every 23 points falling into a single grid cell comprise a cube in Hopf coordinates. Each element of the sequence is obtained by combining the corresponding coordinates in the subspaces, S1 and S2, using the Eq. (4). The dispersion and discrepancy of the resulting sequence can then be computed using the representation for the metric and volume element from Eq. (1) and (5).

5.2 Choosing the Base Resolution onS1 andS2

One of the issues arising when combining the two grids from S1 and S2 is the length of a grid cell edge along each of the coordinates. For this, we have to match the number of cells in each base grids on both of the subspaces, so that they have cell sides of equal lengths [17]. That is, the following equation should hold for m1 and m2:

μ(S1)m1=μ(S2)m2,
(9)

in which μ (S1) is the circumference of the circle S1 and μ (S2) is the surface area of S2. In Appendix 8.3 we explicitly show that both of these values are equal to π.

In our particular case, the base HEALPix grid consists of m2 = 12 cells (Figure 4). Therefore, the number of points in the base resolution of the grid on S1 is m1 = 6. The base grid of the sequence for SO(3) then consists of m1 ·m2 = 6·12 = 72 points (the projections of the grid regions on the Hopf coordinates are shown in Figure 5).

Fig. 4
The base grid of the HEALPix sequence consists of 12 points. The cylindrical projection of the grid cells from S2 to (cos(θ),ϕ) coordinates is shown. Each next resolution subdivides each of the spherical squares into 4 squares of equal ...
Fig. 5
The base grid of the proposed SO(3) sequence consists of 72 points. For the Hopf coordinates (θ,ϕ,ψ) the projections of the grid cells on each of the coordinates are shown. Grid cells for ψ are chosen according to the ordinary ...

5.3 Choosing the Base Ordering

The next step is to choose the ordering of the m = m1m2 points within the base resolution on SO(3). In general, the initial ordering will influence the quality of the resulting sequence, and a method similar to [14] can be used for deciding the ordering of the base sequence.

In our case, we have to define the ordering on the first 72 points of the sequence (see Figure 5 for the illustration of the associated grid regions). In our implementation [9] we have manually selected such an ordering. However, it is possible to design a program that would run through the orderings and select the one that minimizes the discrepancy or any other desired property. For the purpose of further analysis we assume that such an optimal ordering function fbase:N[1,72] is available.

In our implementation [9] we first have selected an ordering of the 12 base points on S2 and 6 base points on S1 (these orderings are shown in Figure 5). For each point on S2 we then generated the 6 points on S1 according to the S1 ordering. The points on S2 are chosen according to the S2 ordering.

5.4 The Sequence

The sequence for SO(3) is constructed one resolution level at a time. The order in which the points from each resolution level are placed in the sequence can be described as follows. The ordering fbase:N[1,,m] of the first m points in the base resolution determines the order of the grid regions within SO(3) and is taken from the previous section. The points in other resolutions fall into the base resolution grids according to the function fbase(i) = fbase(imodm). Every successive m points in the sequence should be placed in these grid regions in the same order. Each of the grid regions is isomorphic to the [0,1]3, and is subdivided into 8 grid regions in each successive resolution. Where exactly each point should be placed within each of the grid regions is determined by the ordering fcube:N[1,8] and recursion procedure defined for the cube [0,1]3 in [13].

The resulting procedure for obtaining the coordinates of the ith element in the sequence is the following:

  1. Assign fbase(i) to be the index of the base grid region that the i-th element has to be placed within.
  2. Assign the floor of the division, icube = [left floor]i/m[right floor], to be the number of subregions already generated in the base grid region. This index then determines the subregion of the region fbase(i) that the i-th element has to be placed within.
  3. Call the recursive procedure from [13] to determine the coordinates of the subregion of the cube [0,1]3 determined by the index icube and the ordering fcube. The i-th element is then placed within this subregion of the fbase(i) region.

5.5 Analysis

Several claims, similar to those obtained in [13], can be made for the new approach. The most important distinction is that the new sequence provides equal volume partition of SO(3) which results in a strong dispersion guarantee.

Proposition 0.1. The dispersion of the sequence T at the resolution level l satisfies:

δ(T)2sin1(12δ2(T2)+(πm12l)2),

in which δ (T2) is the dispersion of the sequence T2 defined over S2.

Proof: It was shown in [17] that the fibers S1 are locally orthogonal to the base space S2 in the sense that an equivalent of the Pythagorean theorem holds for the Hopf coordinates. The bound follows directly from the Pythagorean theorem, and the dispersion bound on the ordinary grid on S1 at the resolution level l. □

Proposition 0.2. The sequence T has the following properties:

  • It is discrepancy-optimal with respect to the set of axis-aligned grid regions defined over S1 and S2.
  • The position of the i-th element of T can be generated in O(logi) time.
  • For any i-th sample, any of the 2d nearest grid neighbors from the same layer can be found in O((logi)/d) time.

Proof: The proof closely follows similar considerations in [13]. □

5.6 Visualization of the Results

To visualize our sequence and compare it with other sequences designed for SO(3), we use the angle-axis, (θ,n), representation from Section 2. It can be shown that if the rotations are uniformly distributed, then the distribution of the angle θ is (sin(θ)−θ)/π. This allows us to draw the elements of SO(3) as the points inside a ball in such a way that every radial line has uniform distribution of elements. This provides a more intuitive visualization, which partially preserves the uniformity. See Figure 6 for visualization of several of the methods of sampling over SO(3), compared to the proposed approach. Specifically, the images show points in the direction of the axis of rotation and with distance to the origin equal to (sin(θ)−θ)/π. Using this representation, the distribution of points increases linearly as a function of distance from the origin. In comparison, a set of points that was uniform with respect to the measure on R3 would have a distribution that varies as the cube of distance from the origin.

Fig. 6
Different sets of samples on SO(3) (a) 2000 random samples (b) 2048 Sukharev grid samples (c) 1944 icosahedral samples (d) 1944 HEALPix samples

6 Implementation and Application to Motion Planning

We have implemented our sampling algorithm in C++ as part of the new library publicly available at [9]. The experiments reported here were performed on a 2.2 GHz Pentium IV running Linux and compiled under GNU C++.

We first compared the uniformity of the new sequence with the Layered Sukharev sequence and random sequence. To demonstrate the importance of understanding the topology of SO(3), we have included the evaluation of uniformity for random Euler angles in this experiment. For each of the deterministic sequences, we generated a fixed set of points, in the range 50 to 100,000, which is shown on the x axis of the graph in Figure 7. We then calculated the distance from a randomly generated point on SO(3) to the nearest neighbor in each of the sets, and selected the largest such distance among 10,000 random points. We have averaged the same computations over 10 runs for each of the random sequences. The obtained value approximates the dispersion. The results are shown for each of the sequences as a separate curve in Figure 7. In all of the cases, the smallest obtained value was the one generated with our new method, which demonstrates that the resulting samples are more uniformly distributed compared to other sequences known for SO(3). Even though it might appear that the actual difference in dispersion is not significant for data sets of a particular size, there is another interpretation of the results on the graph. Consider a particular value of dispersion on the graph, for example, d = 0.06. If a sample set has this dispersion then no ball of radius r > 0.06 can be placed between the samples. To achieve such dispersion, the Hopf sequence required around 50,000 samples, whereas the sample set generated using random Euler angles with twice as many points does not reach the same resolution (see Fig. 7(b)).

Fig. 7
(a) For each of the deterministic sequences, we generated a fixed set of points, in the range 50 to 100,000, which is shown on the x axis of the graph. We then calculated the distance from a randomly generated point on SO(3) to the nearest neighbor in ...

We also used our library as the sampling method in the implementation of PRM-based planner [10] in the Motion Strategy Library [12]. It is important to note that the experiments we present here are just one of possible applications of the developed sequences to motion planning problems. Alternate applications may exist in other areas of computer science, or related fields.

In our experimental setup we consider the rotation-only models for which the configuration space is SO(3). For the two problems shown in Figure 8 we have compared the number of nodes generated by the basic PRM planner using the pseudorandom sequence (with quaternion components [24]), the layered Sukharev grid sequence, and the new sequence. For the first problem the results are: 258, 250, and 248 nodes, respectively. To solve the second problem the PRM planner needed 429, 446 and 410 nodes, respectively. In each trial a fixed, random quaternion rotation was premultiplied to each deterministic sample, to displace the entire sequence. The results obtained were averaged over 50 trials.

Fig. 8
Motion planning problems involving: a) moving a robot (black) from the north pole to the south pole. Multiple views of the geometry of the problem are shown (obstacles are drawn in lighter shades); and b) moving a robot along the corridor.

Based on our results we have observed that the performance of our method is equivalent or better than the performance of the previously known sequences for the basic PRM-based planner. This makes our approach an alternative approach for use in motion planning. It is important to note, however, that for some applications, such as verification problems, only strong resolution guarantees are acceptable.

7 Conclusions

In conclusion, we have developed and implemented a deterministic incremental grid sequence on SO(3) that is highly uniform, can be efficiently generated, and divides the surface of SO(3) into regions of equal volume. Sequences that minimize uniformity criteria, such as dispersion and discrepancy, at each step of generation are especially useful in applications in which the required number of samples is not known in advance.

In the paper we report the performance of the deterministic sequence on motion planning examples for demonstration purposes only. We believe that the value of this method is in providing strong provable guarantees (the bound on dispersion, neighborhood structure, and deterministic generation). These guarantees are not always required in motion planning. Therefore, we do not perform extensive experimental evaluation. However, for the motion planning examples in Section 6, the provable guarantees of the method come at no additional cost. Moreover, it is consistent with the performance of other deterministic sampling methods on Euclidean spaces we have observed on numerous motion planning examples in our previous works. Therefore, we conclude that the sequence can be applied in the context of motion planning, in case deterministic guarantees are required, or in any other applications with such guarantees.

There are a number of ways to improve the current work which we consider as future directions. It is an interesting problem to determine the criteria for an optimal selection of the base sequence on SO(3) to improve the performance of the sequence. It is also tempting to assess the general rate of convergence for motion planning solutions using different sampling sequences.

There are many general open problems related to the presented work. Nicely distributed grids are not yet developed for general n-spheres, n > 3. Implicitly defined manifolds, such as the ones arising from motion planning for closed linkages, are very hard to efficiently and uniformly sample. Such manifolds also arise as the conformation spaces of protein loops. In such cases, efficient parametrization is the bottleneck for developing sampling schemes.

Acknowledgements

Primary funding for this project was provided by the National Science Foundation (NSF CISE-0535007) in support of Steve LaValle and Anna Yershova at UIUC, along with the Alfred P Sloan Foundation and US Department of Energy (DE-FG02-04ER25627) in support of Julie Mitchell. In addition, Anna Yershova's work while at Duke University was supported by the following grants to Bruce R. Donald: NIH grants R01 GM-65982 and R01 GM-078031.

8 Appendix

8.1 Derivation of the Hopf Coordinates

The group of all 3D rotations, SO(3), has a subgroup of all planar rotations, S1. Therefore, SO(3) can be represented as a disjoint union of all of the cosets of this subgroup. Consider a subgroup S1 of all rotations around z axis by angle ψ [set membership] [0,2π). Such rotations can be represented by a unit quaternion of the form:

h1=(cosΨ2,0,0,sinΨ2).

To obtain the left coset of this subgroup, the z axis should be placed in an arbitrary position on the sphere S2. If the sphere is parametrized using the spherical coordinates θ [set membership] [0,π] and ϕ [set membership] [0,2π), then the rotation of the z axis to the (θ, ϕ) position on the sphere corresponds to the quaternion of the form:

h2=(cosθ2,sinθ2sinϕ,sinθ2cosϕ,0).

To obtain the Hopf coordinates for an element x [set membership] SO(3), we compose the two rotations x = h1 * h2 and obtain the expression in Eq. (4):

x=(cosθ2cosΨ2,sinθ2sin(ϕ+Ψ2),cos(ϕ+Ψ2)sinθ2,cosθ2sinΨ2).

8.2 Derivation of the Volume Element Using the Hopf Coordinates

The Jacobian, J = ([partial differential]x/[partial differential]θ,[partial differential]x/[partial differential]ψ,[partial differential]x/[partial differential]ϕ), of the coordinate transformation in Eq. (4) is the following:

J=(12sinθ2cosΨ212cosθ2sinΨ2012sinθ2sinΨ212cosθ2sinΨ2012cosθ2sin(ϕ+Ψ2)12sinθ2cos(ϕ+Ψ2)sinθ2cos(ϕ+Ψ2)12cosθ2cos(ϕ+Ψ2)12sin(ϕ+Ψ2)sinθ2sinθ2sin(ϕ+Ψ2)).

We next compute the volume element of the transformation by taking a square root of the following determinant:

dV=det(JTJ)dθdϕdΨ=18sinθdθdϕdΨ.

8.3 Lengths of the S1 Fibers and the Area of the S2 Base Space for the Hopf Coordinates

By fixing θ and ϕ in Eq. (4) we obtain half circles on the hypersphere S3 (which are S1 fibers on SO(3), after the identifications of the antipodal points are taken into account on S3). To compute the length of a fiber, we follow similar derivations to the Appendix 8.2. For the fixed values of θ and ϕ, the Jacobian of the resulting transformation is the matrix:

J=(xΨ)(12cosθ2sinΨ212cosθ2cosΨ212sinθ2cos(ϕ+Ψ2)12sin(ϕ+Ψ2)sinθ2).

Then the length element for S1 is:

det(JTJ)dΨ=12dΨ.

Therefore, the length of each fiber is:

1202πdΨ=π.

The area of the corresponding base sphere, S2, is then obtained from the volume V(SO(3)) = π2 as the remaining contribution. Therefore, it is π. Note, that these derivations demonstrate that the base space S2 and the fibers S1 are not unit for the Hopf fibration.

References

1. Arvo J. Random rotation matrices. In: Arvo J, editor. Graphics Gems II. Academic Press; Boston: 1991. pp. 355–356.
2. Arvo J. Fast random rotation matrices. Academic Press; Boston: 1992. pp. 117–120.
3. Chirikjian GS, Kyatkin AB. Engineering Applications of Noncommutative Harmonic Analysis. CRC Press; Boca Raton: 2001.
4. Crowther RA. The fast rotation function in the molecular replacement method. Int. Sci. Rev. Ser. 1972;12:173–178.
5. Diaconis P, Shahshahani M. The subgroup algorithm for generating uniform random variables. Prob. in Eng. and Info. Sci. 1987;1:15–32.
6. Fishman GF. Monte Carlo: Concepts, Algorithms, and Applications. Springer-Verlag; Berlin: 1996.
7. Górski KM, Hivon E, Banday AJ, Wandelt BD, Hansen FK, Reinecke M, Bartelmann M. HEALPix: a framework for high-resolution discretization and fast analysis of data distributed on the sphere. arXiv:astro-ph/0409513. 2005 Apr.622:759–771.
8. Hardin DP, Saff EB. Discretizing manifolds via minimum energy points. Notices of the American Mathematical Society. 2004;51(10):1186–1194.
9. Jain S. Library for uniform deterministic sequences/sets of samples over 2-sphere and SO(3) 2009. http://www.cs.duke.edu/~yershova/so3sampling/so3sampling.htm.
10. Kavraki LE, Svestka P, Latombe J, Overmars MH. Probabilistic roadmaps for path planning in high-dimensional configuration spaces. IEEE Trans. Robot. & Autom. 1996 June;12(4):566–580.
11. Kuffner J. IEEE Int. Conf. Robot. & Autom. IEEE; May, 2004. Effective sampling and distance metrics for 3D rigid body path planning.
12. LaValle SM. MSL: Motion strategy library. http://msl.cs.uiuc.edu/msl/
13. Lindemann SR, LaValle SM. Incremental Low-Discrepancy lattice methods for motion planning. IEEE Int'l Conf. on Robotics and Automation.2003. pp. 2920–2927.
14. Lindemann SR, Yershova A, LaValle SM. Workshop on the Algorithmic Foundations of Robotics. 2004. Incremental grid sampling strategies in robotics.
15. Lubotsky A, Phillips R, Sarnak P. Hecke operators and distributing points on the sphere I. Comm. Pure Appl. Math. 1986;XXXIX:S149–S138.
16. Mandell JG, Roberts VA, Pique ME, Kotlovyi V, Mitchell JC, Nelson E, Tsigelny I, Ten Eyck LF. Protein docking using continuum electrostatics and geometric fit. Protein Engineering. 2001;14(2):105–113. [PubMed]
17. Mitchell JC. Sampling rotation groups by successive orthogonal images. SIAM J. Sci. Comput. 2007;30(1):525–547.
18. Niederreiter H. Random Number Generation and Quasi-Monte-Carlo Methods. Society for Industrial and Applied Mathematics; Philadelphia, USA: 1992.
19. Ramamoorthy S, Rajagopal R, Ruan Q, Wenzel L. Proc. of the Workshop on Algorithmic Foundations of Robotics. 2006. Low-discrepancy curves and efficient coverage of space.
20. Rote G, Tichy RF. Spherical dispersion with applications to polygonal approximation of the curves. Anz. Österreich. Akad. Wiss. Math.-Natur, Kl. Abt. II. 1995;132:3–10.
21. Rovira J, Wonka P, Castro F, Sbert M. Point sampling with uniformly distributed lines. Point-Based Graphics, 2005. Eurographics/IEEE VGTC Symposium Proc.Jun, 2005.
22. Saff EB, Kuijlaars ABJ. Distributing many points on a sphere. Mathematical Intelligencer. 1997;19(1):5–14.
23. Shoemake K. Animating rotation with quaternion curves. Proc. of the 12th Annual Conference on Computer Graphics and Interactive Techniques; ACM Press; 1985.
24. Shoemake K. Uniform random rotations. In: Kirk D, editor. Graphics Gems III. Academic Press; 1992. pp. 124–132.
25. Sloane NJA, Hardin TSDRH, Conway JH. Minimal-energy clusters of hard spheres. Disc Comp Geom. 1995;14:237–259.
26. Sternberg MJE, Moont G. Modelling protein-protein and protein-DNA docking. In: Lengauer T, editor. Bioinformatics – From Genomes to Drugs. Wiley-VCH; Weinheim: 2002. pp. 361–404.
27. Stewart GW. The efficient generation of random orthogonal matrices with an application to condition estimation. SIAM J Numer Anal. 1980;17(3):403409.
28. Sun X, Chen Z. Spherical basis functions and uniform distribution of points on spheres. J. Approx. Theory. 2008;151(2):186–207.
29. Wagner G. On a new method for constructing good point sets on spheres. Journal of Discrete and Computational Geometry. 1993;9(1):119–129.
30. Wales D, Doye J. Global optimization by Basin-Hopping and the lowest energy structures of Lennard-Jones clusters containing up to 110 atoms. J Phys Chem A. 1997;101
31. Yan AK, Langmead CJ, Donald BR. A Probability-Based similarity measure for saupe alignment tensors with applications to residual dipolar couplings in NMR structural biology. Int. J. Robot. Res. 2005;24(2–3):162–182.
32. Yang G, Chen I-M. Equivolumetric partition of solid spheres with applications to orientation workspace analysis of robot manipulators. IEEE Transactions on Robotics. 2006 Oct.22(5):869–879.
33. Yershova A, LaValle SM. Deterministic sampling methods for spheres and SO(3); Proc. IEEE International Conference on Robotics and Automation; 2004.
34. Zyczkowski K, Kus M. Random unitary matrices. J. Phys. A. 1994;27(12):4235–4245.