1

Pseudometrically Constrained Centroidal Voronoi Tessellations: Generating
uniform antipodally symmetric points on the unit sphere with a novel acceleration
strategy and its applications to Diffusion and 3D radial MRI

Cheng Guan Koay
Department of Medical Physics
University of Wisconsin School of Medicine and Public Health
Madison, WI 53705

Corresponding author:
Cheng Guan Koay, PhD
Department of Medical Physics
University of Wisconsin School of Medicine and Public Health
1161 Wisconsin Institutes for Medical Research (WIMR)
1111 Highland Avenue
Madison, WI 53705
E-mail: cgkoay@wisc.edu

Word Counts: 4953
Keywords: antipodal symmetry, uniform point set, Centroidal Voronoi Tessellations,
diffusion MRI, 3D radial MRI, uniform points on the sphere.

2

ABSTRACT
Purpose:
The purpose of this work is to investigate the hypothesis that uniform sampling measurements
that are endowed with antipodal symmetry play an important role when the raw data and image
data are related through the Fourier relationship as in q-space diffusion MRI and 3D radial MRI.
Currently, it is extremely challenging to generate large uniform antipodally symmetric point sets
suitable for 3D radial MRI. A novel approach is proposed to solve this important and longstanding problem.
Methods: The proposed method is based upon constrained centroidal Voronoi tessellations of
the upper hemisphere with a novel pseudometric. Geometrically intuitive approach to
tessellating the upper hemisphere is also proposed.
Results: The average time complexity of the proposed centroidal tessellations was shown to be
effectively linear. For small sample size, the proposed method was comparable to the state-ofthe-art iterative method in terms of the uniformity. For large sample size, in which the state-ofthe-art method is infeasible, the reconstructed images from the proposed method has less
streak and ringing artifact as compared to those of the commonly used methods.
Conclusion: This work solved a long-standing problem on generating uniform sampling points
for 3D radial MRI.

3

INTRODUCTION
A centroidal Voronoi tessellation (2) is a Voronoi tessellation (3) in which the center of mass (or
centroid) of the Voronoi region is also its generator and it has been found useful in many
applications, from analysis of cellular pattern (4), neuronal density (5) and territorial behavior of
animals(6) in biological sciences to optimal data acquisition, data quantization, compression and
clustering (7) in the engineering and physical sciences. However, it remains a challenge to
generate these tessellations via Lloyd's algorithm even though Lloyd's algorithm has been found
to be very robust (8).
The main computational bottleneck in the computation of the centroidal Voronoi
tessellations has been the reconstruction step in which the Voronoi tessellations are
reconstructed anew from the newly computed centroids at each iteration.

From the

conventional point of view, if the time complexity of the Voronoi tessellations of a specific
manifold, e.g., plane, sphere or other esoteric surfaces, is O(n log n), then the time complexity
of the centroidal Voronoi tessellations on this manifold must be O(m n log n), i.e., invoking the
O(n log n) algorithm m number times. Note that m is the number of iterations needed to reach
convergence and n is the number of generators.
We have discovered strong heuristic strategies to reduce the average time complexity of
the problem from O(m n log n) to O(m n) by using less optimal, i.e., O(n2), but simpler algorithm
innovatively and will present these strategies in the context of generating uniform antipodally
symmetric points on the sphere. The approach used in this work for accelerating the centroidal
Voronoi tessellations is completely novel, representing a clear departure from the current
paradigms.

4

The problem of generating uniform points on the sphere was proposed by J.J. Thomson
(9) more than a century ago and has spurred many theoretical and computational investigations,
and applications (1,10-15). A variant of the Thomson problem is that of generating uniform
antipodally symmetric points on the sphere (16); the resultant point set plays a particular
important role in MRI, which will be discussed next. It is interesting to note that the problem of
generating nearly uniform antipodally symmetric points on the sphere via a deterministic method
was solved only recently (13,14).
Conceptually, diffusion MRI (17) and 3D radial MRI (18) share many similar features. The
most notable feature relevant to this work is the data redundancy associated with symmetries
inherent in each of these imaging techniques. Specifically, antipodal symmetry of the diffusion
propagator (19) and Hermitian symmetry due to the real-valuedness of the object in the image
domain influence how data are acquired in diffusion MRI and 3D radial MRI, respectively.
Hermitian symmetry and antipodal symmetry are closely related to each other and both
symmetries are defined through parity transformation or spatial inversion with Hermitian
symmetry having an additional operation, which is complex conjugation. Specifically, a realvalued function f

possesses antipodal symmetry if f ( x)  f (  x) and a complex-valued
______

function g possesses Hermitian symmetry if it is a Hermitian function, i.e., g (x)  g ( x) where
complex conjugation is denoted by the overhead bar. Due to these symmetries, it is sufficient to
sample only half of the space, e.g., q-space in diffusion MRI or k-space in MRI. However,
certain constraints on the imaging gradients (20) force the sampling trajectory in 3D radial MRI
to take the usual diametrical line or a curve consists of two straight radial lines adjoined with a
slight non-smooth bend at the center of k-space (21); the former sampling strategy cannot take
advantage of the data redundancy because sampling along a diametrical line through k-space
only is equivalent to acquiring repeated measurements and not new measurements and the

5

latter sampling strategy is more desirable with the potential of acquiring new and different points
in k-space.
The use of uniformity and antipodal symmetry in sampling scheme was advocated by
Jones (16) in diffusion tensor imaging on the basis that the diffusion tensor is symmetric. Here,
we argue that uniformity and antipodal symmetry should be incorporated into sampling design in
diffusion MRI (q-space formalism) and 3D radial MRI (k-space formalism) because both the qspace and k-space possess Hermitian symmetry. Uniform sampling measurements endowed
with antipodal symmetry can maximize sampling coverage and at the same time take advantage
of the inherent symmetry of the imaging systems. It is to the best of our knowledge that none of
the works in 3D radial MRI (18,21-28) used antipodally symmetric point sets to acquire k-space
data. In fact, many studies (22,23,25,27) used the method of Wong and Roos (15) or the
generalized spiral scheme (10,11). The point sets generated from any of the methods, i.e.,
generalized spiral scheme, the method of Wong and Roos and the method of Bauer (12), have
the spiral pattern although each of these methods was formulated and derived differently.
Antipodally symmetric point sets have not been used in any 3D radial MRI study may be due to
the following reasons: the lack of a robust iterative method or a deterministic method capable of
generating large number of antipodally symmetric points on the sphere, and the availability of
other deterministic methods (1,11,12,15,29) for generating non-antipodally symmetric point sets.
Ad hoc strategy of collecting half of the non-antipodally symmetric points on one of the
hemispheres is usually adopted to come up with an antipodally symmetric point set. We should
note that the diffusion MRI sampling method based on vertices of a multifold-tessellated
icosahedral hemisphere was used by Tuch et al. (30) but this method is inflexible and limited in
that it cannot generate point set of any size except certain "magic" numbers, e.g., 6, 26, 46, 91,
126, 196, 246, 341, 406, 526 and so on. Note that fivefold-tessellated icosahedral hemisphere

6

produces 126 points on the upper hemisphere, see Tuch (30). Due to the above mentioned
limitation, this sampling method used by Tuch et al. may not be very appealing to the
practitioners of 3D radial MRI.
Even with the advent of the deterministic method capable of generating nearly uniform
antipodally symmetric points on the sphere (13,14), it is still extremely challenging to generate
highly uniform antipodally symmetric point set for use in 3D radial MRI through iterative function
minimization even when the initial solution is taken from the deterministic methods. The number
of points (on the sphere) required in 3D radial MRI and diffusion MRI is in the thousands (4000
to 20000 or more) and in the hundreds (6 to 500), respectively. Therefore, the number of
parameters in spherical coordinate representation to be optimized are twice as many (8000 to
40000) and optimization methods that rely on approximate or exact Hessian or its inverse are
completely infeasible. While some of the iteratively optimized antipodally symmetric point sets
(with sample size in the hundreds) have been tabulated and available for public use (31), the
lack of highly uniform and antipodally symmetric point sets applicable to 3D radial MRI
motivates the present investigation.
Here, we propose a constrained centroidal Voronoi tessellation endowed with a
pseudometric to robustly and efficiently generating uniform antipodally symmetric point sets on
the unit sphere. The uniformity of the point sets generated from the proposed method is shown
to be comparable to the state-of-the-art iterative method (31) for small sample size relevant to
diffusion MRI but not to 3D radial MRI. For large point set with sample size in the range of
thousands, which is relevant to 3D radial MRI and other applications (32), the state-of-the-art
method is completely infeasible but the proposed approach remains feasible. Further, we used
our recently developed three-dimensional analytical MRI phantom in the Fourier domain (33),
which was based upon the three-dimensional version of the famous Shepp-Logan phantom

7

(34,35) in the image domain, to compare and contrast the qualitative features of the
reconstructed images obtained through the proposed method, the method of Wong and Roos,
the generalized scheme and the tessellated-icosahedral of the upper hemisphere (the method
used by Tuch et al.). We found that the reconstructed image from the proposed method has less
streak and ringing artifacts as compared to those of other methods.

8

METHODS
A centroidal Voronoi tessellation (2) is a Voronoi tessellation in which the center of mass (or
centroid) of the Voronoi region is also its generator and the Voronoi regions can be prescribed
with a density function. To the best of our knowledge, a centroidal Voronoi tessellation capable
of generating uniform antipodally symmetric points on the unit sphere has never been proposed.
Instead of dealing with a density function that is constant and invariant with respect to spatial
inversion or antipodally symmetry, it is more convenient and efficient to introduce a novel
pseudometric in the centroidal Voronoi tessellation so as to make generating uniform antipodally
symmetric points on the unit sphere possible.
The spherical Voronoi regions, { Vi }iN1 , on the upper hemisphere are characterized by
the following properties:


No two distinct regions share the same point. That is, the intersection between any two
distinct regions is empty, Vi  V j   with i  j . Points on the boundary between any
two regions belong to the closure of these regions, which is denoted by an overhead bar,
e.g., Vi .



The union of all the closures of the spherical Voronoi regions covers the upper
hemisphere, denoted by S 2 .



Vi  {x  S2 | d (x, gi )  d (x, g j ) for

j  1,..., N

and

j  i}

[1]

Each of the unit vectors, { g i }iN1 , on S 2 is called the generator of its respective Voronoi
region. The pseudometric, d , will be defined later.
The center of mass of each Voronoi region does not necessarily coincide with the
generator of that region. An iterative method such as the Lloyd algorithm (8,36) may be used to

9

make the generators from each successive iterations closer and closer to the centers of mass of their
respective regions. The resultant tessellations are known as the centroidal Voronoi tessellations (2).
The center of mass of a spherical Voronoi region,

ĝi 

 vd i
 vd i

Vi , can be expressed as:

.

[2]

Note that  i is the spherical surface of Vi , v represents unit vector normal to the spherical
surface element d i and . denotes the Euclidean norm, which is needed to ensure ĝ i is unit
length. In practice, the computation of ĝ i is based on the discretized version of Eq.[2].
Specifically, it is obtained through the sum of the products between the area of the spherical
triangle formed by the generator and each pair of consecutive vertices at the boundary
surrounding the generator in counterclockwise order. We further note that the density function,
which usually appears as a factor in the integrand of Eq.[2], is taken to be a unit constant
function in order to ensure uniformity of the generators.
The distance measure or metric, d (,) , used in this work is a novel extension of the
modified electrostatic potential energy term studied in Ref.(37). For completeness, we will
introduce the notion of real and virtual points for manipulating the antipodally symmetric point
set suggested in Ref.(37). Due to the constraint of antipodal symmetry, we classify points as
real and their corresponding antipodal points as virtual. If we have N real points (on the upper
hemisphere), denoted by unit vectors ri with i  1,  , N , then the total electrostatic energy for
the complete configuration (37) of 2 N points of both real and virtual points on the whole sphere
is given by:

10

 

N
2


 1
 2 
 r 
i  1 j  i  1  ij

N 1

N




4  rij2 

1

[3]

with rij  ri  r j . Note that Eq.[3] is expressed solely in terms of real points. If we define

S (ri , r j ) 

1
1
then S (ri , r j ) may be thought of as a reciprocal metric (or reciprocal of

2
rij
4  rij

the distance measure) between two real points.). Here, we exploit this reciprocal metric fully by
defining our pseudometric, d (,) , as

d (ri , r j )  1 / S (ri , r j ) .

[4]

Note that d (r j , r j )  0  d (r j ,r j ) , hence the term pseudometric. Please refer to Appendix B
for the proof that d (,) is indeed a pseudometric.
We should note that implementation of spherical Voronoi tessellations is not a trivial
computational task (38,39) and our approach is different from the existing methods,
Refs.(38,39). To compare and contrast our proposed approach, we will first present an outline of
the Lloyd's algorithm below:
Pseudometrically Constrained Centroidal Voronoi Tessellation:
Let n be the number of desired points on the upper hemisphere.
Step 1. Deterministically generate 2 n highly uniform points on the unit sphere
via Analytically Exact Spiral Scheme(1) and select those points that are
on the upper hemisphere as the generators.
Step 2. Construct the spherical Voronoi regions of the upper hemisphere with the
chosen generators.
Step 3. Compute the normalized centroids using the discretized version of Eq.[2].
Step 4. Adopt the normalized centroids as the generators.
Step 5. Iterate Steps 2, 3, 4 until convergence it reached.

11

We would like to point out that time complexity of Step 2 above is O(n log n), see Ref.(39).
Therefore, the time complexity of the centroidal Voronoi tessellation is O(m n log n) where m is
the number of iterations needed to reach convergence.
Our approach to generating the Voronoi tessellations of the upper hemisphere is based
on the following steps and observations:
1. The Voronoi region of a generator can be found by first identifying surrounding
generators (approximately neighbors of neighbors) within a spherical cap with a
prescribed radius, which depends on the total number of generators, away from
the generator of interest, see Figure 1A. This step can always be done because
the initial point set obtained from the Analytically Exact Spiral Scheme(1) is
already reasonably uniform. When the total number of generators is large, the
number of surrounding generators remain constant and is around 20. The time
complexity of this step for all the generators is O(n2).
2. The generator and its surrounding generators are rotated in such a way that the
generator is situated along the z-axis, see Figure 1A. The time complexity of this
step for all the generators is O(n).
3. To compute the vertices of the Voronoi region, the smallest convex region formed
by the surrounding generators and around the generator is needed, see Figure
1B. Finding this smallest convex region is equivalent to finding the boundary
points of the convex hull of the surrounding generators after these generators
have been stereographically projected (40) onto the x-y plane, see Figure 1C. The
validity of this step hinges on the angle preserving (or conformal) property of
steoreographic projection (41). The boundary points of the convex hull of a set of
points on the plane is then computed through the Graham's scan, (42). The time

12

complexity of Graham's scan is O(k log k) and k is the number of points in the
convex hull, i.e., the number of surrounding generators, which is about 20.
Therefore, the time complexity of this step for all the generators is O(n) because k
is a constant.
4. The vertices of the Voronoi region, the area of spherical triangles formed by the
Voronoi vertices and the generator can be computed in O(n) for all the generators.
The steps proposed above are intuitive and geometrically motivated but we should note that
Step 1 in our proposed approach is O(n2), which is the main bottleneck of our approach if this
step were to be invoked at every iteration. In what follows, we will describe a simple way to
avoid tessellating at every iteration and provide empirical evidence that our approach is
effectively O(n) for generating the centroidal Voronoi tessellations. We can compute the
distance between a generator at the current iteration and the same generator from the previous
iteration. This information has already been used to determine convergence, i.e., Ref.(8). The
novelty of our proposal is in using this information to compute the cumulative sum of Euclidean
distances made by each generators. If the maximum value of these cumulative sums is less
than some prescribed value, see Appendix A, then Step 1 will not be invoked. Otherwise, Step 1
will be invoked and the cumulative sum is reset to zero. When Step 1 is not invoked, the
connectivity network between a generator and its surrounding generators will not be altered but
the coordinates of the generator and its surrounding generators will likely be different from one
iteration to the next.
Here, we discuss other important implementation details. Note that some of the vertices
around the equator may be located at the lower hemisphere and the resultant centroids may
also be on the lower hemisphere. Therefore, the centroids or generators should be reoriented

13

onto S 2 at each iteration. Finally, we note that the convergence of the proposed algorithm is
based on that of local deviation as described in Ref.(8).

14

RESULTS
We implemented the proposed pseudometrically constrained centroidal Voronoi tessellation in
Java. We conducted four tests to illustrate the proposed method. The first two tests were run on
a machine with an Intel® Core™ i7 CPU at 1.73 GHz and with 8 GB of RAM and the last two
tests were run on a different machine, Four X 8-Core 2.3GHz AMD Opteron Processors (6134)
with 128GB of RAM.
The first test is to show that the performance of the point sets generated by the proposed
method is of comparable quality in terms of modified electrostatic energy to that of the state-ofthe-art iterative method and to show that the proposed centroidal Voronoi algorithm is O(m n)
through empirical analysis of the execution time per iteration as a function of number of
generators (or points on the upper hemisphere), which shows a linear trend, see Figure 2.
Figure 2A shows the percent of relative error in terms of the modified electrostatic energy
of initial deterministic generators, which were obtained our recent analytically exact spiral
scheme(1) and of the final generators obtained by the proposed method with respect to the
state-of-the-art method. At the level of percent of relative error achieved by the proposed
method, the uniformity of the points between the proposed method and the state-of-the-art
method is indistinguishable visually or quantitatively. However, we observed that the point sets
generated centroidal Voronoi tessellations has consistently shown to be higher in modified
electrostatic energy than those obtained through iterative scheme. We believe this observation
has to do with the intrinsic property of centroidal Voronoi tessellations, that is the fundamental
centroidal constraint itself. The centroidal constraint forces generators to be in geometrically
frustrated positions even though the effect of the centroidal constraint on geometric frustration is
very slight.

15

Figure 2B shows the performance of the proposed method in terms of execution time (in
seconds) per iteration of the proposed method as a function of the number of points on the
upper hemisphere. It is clear that the execution time per iteration has a linear trend with respect
to N.
Figure 2C shows the frequency of invocation of Step 1, which is an O(N2) algorithm, as a
function of N, number of points on the upper hemisphere and Figure 2D shows the total number
of iterations as a function of N and the instances of invocation of Step1 that occurred within
these iterations are shown in red. Note that the highest peak is at (N=324, 4161 iterations).
The second test is to give an example of a point set with sample size large enough that
the state-of-the-art method is currently not feasible. In this example, we used sample size of
888. Figure 3A shows the initial generators and their Voronoi regions as obtained the
analytically exact spiral scheme(1). Figures 3B and 3C shows the centroids and centroidal
Voronoi regions on the upper and the lower hemispheres, respectively. Figures 3D shows the
centroids on the whole sphere by combining Figures 3B and 3C.
The third test is intended to show the feasibility of the proposed method in generating point
set of sample size large enough for 3D MRI applications and beyond. Figure 4 shows only the
Voronoi regions on the sphere with 16000 generators on the upper hemisphere. It took 10.38
minutes to generate the centroidal Voronoi tessellations (at the tolerance level of 10-12) on the
machine with Four X 8-Core 2.3GHz AMD Opteron Processors (6134).
The final test is intended to show the qualitative features of the reconstructed images of
the three-dimensional analytical MRI phantom (33) obtained from the proposed method, the
method of Wong and Roos, the generalized spiral scheme and the tessellated icosahedral of
the upper hemisphere. Note that the reconstruction algorithm is based upon the following nonuniform discrete Fourier transform for the image domain signal at r :

16

m

I (r )   S (k i ) | k i |2 exp(2 i k i  r ) ,
i 1

where S (k i ) is the signal value of the three-dimensional Shepp-Logan phantom evaluated
analytically at

k i , see (33). The term, | k i |2 , came about from the spherical coordinate

transformation in the Fourier expansion. Two target matrix sizes of the reconstructed imaging
volume were chosen and they were 128  128  128 (standard resolution) and 256  256  256
(high resolution). We have m  526  64  33664 , kmax  32 (arbitrary units) and k  0.5
(arbitrary units) at the standard resolution and m  526  128  67328 , kmax  64 (arbitrary
units) and k  0.5 (arbitrary units) at the high resolution. Figure 5 shows the reconstructed
images generated from the different schemes studied in this work at the standard resolution (the
first two rows) and at the high resolution (the last two rows). At the standard resolution of
128x128x128 with z = -0.181, the images generated by the method of Wong and Roos and the
generalized spiral scheme have more noticeable ringing artifacts around the bright region than
that of the proposed method. The image generated by the tessellated icosahedral scheme has
more streak artifacts in the dark regions than that of the proposed method. At the same
resolution with z = 0.228, the images generated from these methods have more ringing artifacts
than that of the proposed method. Similar patterns of artifacts showed up more noticeably at a
higher resolution of 256x256x256. At this resolution, the images generated from the proposed
method have less ringing artifacts than those from other methods.

17

DISCUSSION
Even though our recent works (13,14) were the first to solve the problem of deterministically
distributing antipodally symmetric points on the unit sphere, the generated points are structurally
very regular and are prone to causing degeneracy in triangulation (43), i.e., when four points
rather than three sit on the same circumcircle, and hence, point sets from these works could not
be used as the initial generators for the proposed method. The point sets generated from the
analytically exact spiral scheme did not have the above mentioned defect and contribute
immensely to the feasibility of the present approach by providing highly efficient technique for
generating nearly uniform points on the sphere and half of which are then used as the initial
generators in the present approach.
The novelty of this work lies in the realization that the reciprocal of the reciprocal metric
can be treated as a pseudometric in the antipodally symmetric space, in the extension of the
centroidal Voronoi tessellations to the antipodally symmetric space, in the heuristic strategies
suggested above for accelerating centroidal Voronoi tessellations without tessellating at every
iteration, and in the utilization of the proposed constrained centroidal Voronoi tessellations to
generate uniform antipodally symmetric points on the sphere. It should be clear that the
pseudometric used in this work can also be used in any K-mean clustering algorithm of discrete
data that are endowed with antipodal symmetry. The results showed that antipodally symmetric
uniform sampling strategy does play a positive role in image quality and such a strategy should
be incorporated into sampling design considerations of any 3D radial MRI study.
One of the most exciting results of this work is the heuristic strategies used to accelerate
centroidal Voronoi tessellations. This result highlights the emergence of strong heuristics in
which traditional optimal algorithm that is invoked iteratively may not outperform less optimal
algorithm that is invoked sparingly in an innovative way to accomplish the same computational

18

task. For small sample size, the point sets generated from the proposed method are
comparable to those generated from the state-of-the-art method in all the cases tested. One of
the important highlights of this work is that the proposed method is capable of generating large
uniform point sets relevant to 3D radial MRI applications; the state-of-the-art iterative method is
completely infeasible in this case particularly when the desired number of points is in the
thousands.
We must note that the reconstruction algorithm adopted in this work is for testing
purposes as it is computational expensive but it serves the present purpose very well because
whatever the differences in the reconstructed images these differences have to come from
uniformity or lack thereof of each of the methods studied in this work. It is not fortuitous that the
reconstructed image from the proposed method has less streak and ringing artifacts as
compared to those of other methods studied in this work.
Uniform antipodally symmetric point sets have continued to play a major role in MRI
acquisition design (14,16,23,37,44-49) and are beginning to make an impact in other fields (32).
We hope this work contributes to further advances to these scientific endeavors beyond
diffusion and 3D radial MRI.

19

Acknowledgment
The author dedicates this work to Lean Choo Khoo. He would like to thank Drs. M. Elizabeth
Meyerand, Charles A. Mistretta and Peter J. Basser for the encouragement. Thanks go to Dr.
Orhan Unal, Director of Medical Physics Computing, for making available the computing
resources needed for this work and for providing the specifications of the Opteron Processors.
The author would like to thank Drs. Steven R. Kecskemeti, Kevin M Johnson, and Michael A.
Speidel for helpful discussion on image reconstruction and image quality assessment. This work
was supported in part by the National Institutes of Health IRCMH090912. The software related
to

this

work

will

be

made

freely

http://sites.google.com/site/hispeedpackets/ .

Appendix A

available

for

research

use

at

URL:

20

In this appendix, we present a method, which is an special case of our previous work(50), to
compute the radius of a spherical cap so as to cover sufficient number of surrounding
generators for the purposes of computing Voronoi region. We also mention the criterion adopted
to decide whether to invoke Step 1 or otherwise.
If the center of a circle of radius r on a plane that is tangent to a unit sphere at point
(0,0,1) or on the z-axis, then its inverse Gnomonic projection is a spherical cap, see Figure 2A
of Ref.(50). Since the area of the upper hemisphere of unit radius is 2 , the area associated
with a generator can be approximated by 2 / N

where N is the total number of generators

on the upper hemisphere. This approximate works well for large N , i.e., 15 and above. The
unnormalized areal measure, Eq.(33) of Ref.(50), of a spherical cap with area equals to 2 / N
is related to r by the following expression:

2 / N  2 (1 

Therefore, r 

1 ).
1 r 2

NN1 2  1 .

The angle subtended by an longest arc passing through the interior of the spherical cap is given
by:

  2 tan 1( r ) .
This angle is also the approximate spherical distance between two generators. Therefore, any
generators that are within 5 / 2 in spherical distance away for the generator of interest will be
classified as surrounding generators.

Previously, we mentioned that if the maximum value of the cumulative sums of distances
made by generators in successive iterations is less than some prescribed threshold value,

21

which is a small fraction of the prescribed radius in Step 1, then Step 1 will not be invoked. The
prescribed value is based on Euclidean distance for simplicity as the local convergence criterion
is also based on Euclidean distance. The prescribed threshold value is given

C 2(1  cos(tan 1( r ))) where C is a positive real number less than unity. The specific value of
C adopted in this study was C  0.15 . Note that if C assumes too high a value, one may run
into the problem in which the initial surrounding generators may not be the current surrounding
generators, which may cause the error in Voronoi tessellations. If C assumes too low a value,
Step1 may be invoked many more times than necessary.

Appendix B: The proposed pseudometric for the unit sphere

22

The distance measure, d (,) , proposed in this work is a novel extension of the modified
electrostatic potential energy term studied in (37). For the sake of completeness, the notion of
real and virtual points for manipulating the antipodally symmetric point set as suggested in (37)
will be reintroduced here.
Due to the constraint of antipodal symmetry, we classify points as real and their
corresponding antipodal points as virtual. If we have N real points (say on the upper
hemisphere), denoted by unit vectors ri with i  1, , N , then the total electrostatic energy for
the whole configuration (37) of 2 N points of both real and virtual points on the whole sphere is
given by:



N
2


1
1
2    
4  rij2
i 1 j i 1  rij

N 1 N



,



[B1]

with

 ri  r j
rij  
 ri  r j

if
if

arccos(ri  r j )   / 2
.
arccos(ri  r j )   / 2

[B2]

We should note that we could have defined rij to be simply ri  r j , which is in fact the actual
definition used in the implementation of the proposed algorithm, and the expression in Eq.[B1]
would still be valid but the definition of rij in Eq.[B2] would bring much clarity to the proof later.
Note that Eq.[B1] is expressed solely in terms of real points. The interaction term in Eq.[B1],

1
I (ri , r j )  1 
rij
4  rij2

[B3]

was first postulated in (37) as a reciprocal metric (or reciprocal of the distance measure)
between two real points.). Here, we prove that m(,) given by

23

m(ri , r j )  1 / I (ri , r j ) 

1
rij



1

[B4]

1

4  rij2

is indeed a pseudometric.
Before we begin the proof to show that m (,) is a pseudometric on the unit sphere, we
first list all the necessary information about (1) the pseudometric conditions on the unit sphere,
(2) concave functions, (3) sub-additive functions and (4) a well known property about the
vertices of any spherical triangle on the unit sphere that is endowed with antipodal symmetry.

Definition 2.1. Let S 2 be the three-dimensional unit sphere. A pseudometric on S 2 is a
mapping, d : S 2  S 2  [0,2] , such that
(a) d (ri , r j )  0 ,
(b) d (ri , r j )  d (r j , ri ) ,
(c) ri  r j implies d (ri , r j )  0 , and
(d) d (ri , r j )  d (ri , rk )  d (r j , rk ) .
Note that ri , r j , and rk are unit vectors. Note further that the d (,) is considered a metric if it
satisfies

the

above

conditions

and

also

satisfies

the

converse

of

(c),

i.e.,

d (ri , r j )  0 implies ri  r j . In brief, a metric is automatically a pseudometric but the converse
is not generally true. It is not hard to see that m (,) satisfies the first three conditions trivially by
virtue of Eq.[B2]. The proof will focus on showing that m (,) also satisfies the last condition, i.e.,
the triangle inequality.

24

Definition 2.2. A function g is concave on an interval [ a , b] if, for any x and y in [ a , b] and for
any   [0,1] , we have

g ( x  (1   ) y )  g ( x)  (1   ) g ( y ) .

[B5]

Lemma 2.1. If a function g is concave on an interval [0, b] for some b  0 and g (0)  0 , then

g is subadditive, i.e., g ( x  y )  g ( x )  g ( y ) for any x, y  0 .
Proof. By setting y  0 and invoking the fact that g is concave and g (0)  0 , Eq.[B5] reduces
to

g ( x )  g ( x ) .

[B6]

Therefore,
y
y
g ( x  y)  x g ( x  y) 
g ( x  y )  g ( x ( x  y ))  g (
( x  y ))  g ( x)  g ( y ).
x y

x y

x y

x y

The following lemma is well known and it is stated for completeness.
Lemma 2.2. Any three points, ri , r j , and rk , on S 2 that are endowed with antipodal symmetry,
two of the eight triplets, ( ri ,r j ,rk ) , have vertices bounded by a spherical octant.

ri , ~
rj , and ~
rk , satisfies the following
Consequently, any one of these two triplets, say, ~
properties:

arccos((~
rk  ~
ri )  (~
rj  ~
ri ))   / 2 ,
arccos((~
rk  ~
rj )  (~
ri  ~
r j ))   / 2 , and
arccos((~
rj  ~
rk )  (~
ri  ~
rk ))   / 2 .

25

If you define a new function based on Eq.[B4] as

f (r ) 

1

r

1

1

,

[B7]

4r 2

it should be clear that f is concave on [0,2] because f is continuous and twice-differentiable,
and

d2 f
dr 2

 0 for any r  (0,2) . Further, f (0)  0 . Therefore, f is sub-additive. We should

note that f is an increasing function on [0, 2 ] and a decreasing function on [ 2 ,2] .

ri , ~
rj ,
To prove the triangle inequality, we invoke Lemma 2.2 to focus on one of the triplets, i.e., ~

rk , and write
and ~
~ ~

~ ~

~ ~

~ ~

( r  r )( r  r )
( r  r )( r  r )
~
ri  ~
r j  a  b  2 , with a  i ~j ~ k j  0 , and b  j ~i ~ k i  0 .
r r
r r
i



j

i

j



By Lemma 2.1, f ~
ri  ~
r j  f (a  b)  f (a )  f (b) .
By Lemma 2.2 again, we have
ri  ~
r j )( ~
rk  ~
rj)
rj  ~
ri )( ~
rk  ~
ri )
(~
(~
,
and

1
 1 , and therefore,
~
~
ri  ~
rj ~
rk  ~
rj
ri  ~
rj ~
rk  ~
ri

a ~
rk  ~
r j  2 and b  ~
rk  ~
ri  2 .
Note that





f (a )  f (b)  f ~
rk  ~
rj  f~
rk  ~
ri



because f is an increasing function on [0, 2 ] . Therefore,



 



f ~
ri  ~
rj  f ~
rk  ~
rj  f~
rk  ~
ri  .
We conclude that m (,) satisfies the triangular inequality and therefore, it is a pseudometric.

26

References
1.
2.
3.

Koay CG. Analytically exact spiral scheme for generating uniformly distributed points on the unit sphere. Journal of Computational Science
2011;2(1):88-91.
Du Q, Faber V, Gunzburger M. Centroidal Voronoi Tessellations: Applications and Algorithms. SIAM Review 1999;41(4):637-676.
Okabe A, Boots B, Sugihara K, Chiu DSN. Spatial Tessellations: Concepts and Applications of Voronoi Diagrams: John Wiley & Sons;
2009.

27

4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
31.
32.
33.
34.
35.
36.
37.
38.
39.
40.

Honda H. Description of cellular patterns by Dirichlet domains: The two-dimensional case. Journal of Theoretical Biology 1978;72(3):523543.
Duyckaerts C, Godefroy G, Hauw J-J. Evaluation of neuronal numerical density by Dirichlet tessellation. Journal of Neuroscience Methods
1994;51(1):47-69.
Barlow GW. Hexagonal territories. Animal Behaviour 1974;22, Part 4(0):876-IN871.
Gray RM, Kieffer JC, Linde Y. Locally optimal block quantizer design. Information and Control 1980;45(2):178-198.
Du Q, Emelianenko M, Ju L. Convergence of the Lloyd Algorithm for Computing Centroidal Voronoi Tessellations. SIAM Journal on
Numerical Analysis 2006;44(1):102-119.
Thomson JJ. On the structure of the atom: an investigation of the stability and periods of oscillation of a number of corpuscles arranged at
equal intervals around the circumference of a circle; with application of the results to the theory of atomic structure. Philosophical Magazine
1904;7(39):237-265.
Rakhmanov E, Saff E, Zhou Y. Minimal discrete energy on the sphere. Mathematical Research Letters 1994;1:647-662.
Saff E, Kuijlaars A. Distributing many points on a sphere. The Mathematical Intelligencer 1997;19:5-11.
Bauer R. Distribution of points on a sphere with application to star catalogs. Journal of guidance control and dynamic 2000;23:130-137.
Koay CG. A simple scheme for generating nearly uniform distribution of antipodally symmetric points on the unit sphere. Journal of
Computational Science 2011(2):376-380.
Koay CG, Hurley SA, Meyerand ME. Extremely efficient and deterministic approach to generating optimal ordering of diffusion MRI
measurements. Med Phys 2011;38(8):4795-4801.
Wong S, Roos M. A strategy for sampling on a sphere applied to 3D selective RF pulse design. Magnetic Resonance in Medicine
1994;32:778-784.
Jones DK, Horsfield MA, Simmons A. Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging.
Magnetic Resonance in Medicine 1999;42(3):515-525.
Basser PJ, Mattiello J, Le Bihan D. MR diffusion tensor spectroscopy and imaging. Biophys J 1994;66(1):259-267.
Lai CM, Lauterbur PC. TRUE 3-DIMENSIONAL IMAGE-RECONSTRUCTION BY NUCLEAR MAGNETIC-RESONANCE
ZEUGMATOGRAPHY. Phys Med Biol 1981;26(5):851-856.
Callaghan PT. Principles of nuclear magnetic resonance microscopy. New York: Oxford University Press; 1991.
Lustig M, Seung-Jean K, Pauly JM. A Fast Method for Designing Time-Optimal Gradient Waveforms for Arbitrary <formula> <img
src="/images/tex/348.gif" alt="k"> </formula>-Space Trajectories. Medical Imaging, IEEE Transactions on 2008;27(6):866-873.
Liu J, Redmond MJ, Brodsky EK, Alexander AL, Lu A, Thornton FJ, Schulte MJ, Grist TM, Pipe JG, Block WF. Generation and
visualization of four-dimensional MR angiography data using an undersampled 3-D projection trajectory. Medical Imaging, IEEE
Transactions on 2006;25(2):148-157.
Mistretta CA, Wieben O, Velikina J, Block W, Perry J, Wu Y, Johnson K. Highly constrained backprojection for time-resolved MRI.
Magnetic Resonance in Medicine 2006;55(1):30-40.
Barger AV, Block WF, Toropov Y, Grist TM, Mistretta CA. Time-resolved contrast-enhanced imaging with isotropic resolution and broad
coverage using an undersampled 3D projection trajectory. Magnetic Resonance in Medicine 2002;48(2):297-305.
Stehning C, Börnert P, Nehrke K, Eggers H, Dössel O. Fast isotropic volumetric coronary MR angiography using free-breathing 3D radial
balanced FFE acquisition. Magnetic Resonance in Medicine 2004;52(1):197-203.
Nagel AM, Laun FB, Weber M-A, Matthies C, Semmler W, Schad LR. Sodium MRI using a density-adapted 3D radial acquisition
technique. Magnetic Resonance in Medicine 2009;62(6):1565-1573.
Rahmer J, Börnert P, Groen J, Bos C. Three-dimensional radial ultrashort echo-time imaging with T2 adapted sampling. Magnetic
Resonance in Medicine 2006;55(5):1075-1082.
Nielles-Vallespin S, Weber M-A, Bock M, Bongers A, Speier P, Combs SE, Wöhrle J, Lehmann-Horn F, Essig M, Schad LR. 3D radial
projection technique with ultrashort echo times for sodium MRI: Clinical applications in human brain and skeletal muscle. Magnetic
Resonance in Medicine 2007;57(1):74-81.
Piccini D, Littmann A, Nielles-Vallespin S, Zenge MO. Spiral phyllotaxis: The natural way to construct a 3D radial trajectory in MRI.
Magnetic Resonance in Medicine 2011;66(4):1049-1056.
Ahmad R, Deng Y, Vikram DS, Clymer B, Srinivasan P, Zweier JL, Kuppusamy P. Quasi Monte Carlo-based isotropic distribution of
gradient directions for improved reconstruction quality of 3D EPR imaging. Journal of Magnetic Resonance 2007;184(2):236-245.
Tuch DS, Reese TG, Wiegell MR, Makris N, Belliveau JW, Wedeen VJ. High angular resolution diffusion imaging reveals intravoxel white
matter fiber heterogeneity. Magnetic Resonance in Medicine 2002;48(4):577-582.
Cook PA, Bai Y, Nedjati-Gilani S, Seunarine KK, Hall MG, Parker GJ, Alexander DC. Camino: Open-source Diffusion-MRI
Reconstruction and Processing. 2006; 14th Scientific Meeting of the International Society for Magnetic Resonance in Medicine, Seattle,
WA, USA. p 2759.
King C, Brown WR, Geller MJ, Kenyon SJ. Identifying star streams in the milky way halo. Astrophysical Journal 2012;750(1).
Koay CG, Sarls JE, Ozarslan E. Three-dimensional analytical magnetic resonance Imaging phantom in the Fourier domain. Magnetic
Resonance in Medicine 2007;58(2):430-436.
Shepp LA, Logan BF. Reconstructing Interior Head Tissue from X-Ray Transmissions. Nuclear Science, IEEE Transactions on
1974;21(1):228-236.
Shepp LA. Computerized Tomography and Nuclear Magnetic Resonance. Journal of Computer Assisted Tomography 1980;4(1):94-107.
Lloyd S. Least squares quantization in PCM. Information Theory, IEEE Transactions on 1982;28(2):129-137.
Koay CG, Özarslan E, Johnson KM, Meyerand ME. Sparse and optimal acquisition design for diffusion MRI and beyond. Med Phys
2012;39(5):2499-2511.
Augenbaum JM, Peskin CS. On the construction of the Voronoi mesh on a sphere. Journal of Computational Physics 1985;59(2):177-192.
Renka RJ. Algorithm 772: STRIPACK: Delaunay triangulation and Voronoi diagram on the surface of a sphere. ACM Trans Math Softw
1997;23(3):416-434.
Coxeter HSM. Introduction to geometry. New York: Wiley; 1989.

28

41.
42.
43.
44.
45.
46.
47.
48.
49.
50.

Pólya G, Latta G. Complex variables. New York: Wiley; 1974.
Graham RL. An efficient algorith for determining the convex hull of a finite planar set. Information Processing Letters 1972;1(4):132-133.
Beichl I. Dealing with degeneracy in triangulation. Computing in Science & Engineering 2002;4(6):70-74.
Dubois J, Poupon C, Lethimonnier F, Le Bihan D. Optimized diffusion gradient orientation schemes for corrupted clinical DTI data sets.
Magnetic Resonance Materials in Physics, Biology and Medicine 2006;19(3):134-143.
Cook PA, Symms M, Boulby PA, Alexander DC. Optimal acquisition orders of diffusion-weighted MRI measurements. Journal of
Magnetic Resonance Imaging 2007;25(5):1051-1058.
Koay CG, Chang L-C, Carew JD, Pierpaoli C, Basser PJ. A unifying theoretical and algorithmic framework for least squares methods of
estimation in diffusion tensor imaging. Journal of Magnetic Resonance 2006;182(1):115-125.
Koay CG, Özarslan E, Basser PJ. A signal transformational framework for breaking the noise floor and its applications in MRI. Journal of
Magnetic Resonance 2009;197(4):108-119.
Deriche R, Calder J, Descoteaux M. Optimal real-time Q-ball imaging using regularized Kalman filtering with incremental orientation sets.
Medical Image Analysis 2009;13(4):564-579.
Nett EJ, Johnson KM, Frydrychowicz A, Del Rio AM, Schrauben E, Francois CJ, Wieben O. Four-dimensional phase contrast MRI with
accelerated dual velocity encoding. Journal of Magnetic Resonance Imaging 2012;35(6):1462-1471.
Koay CG, Nevo U, Chang LC, Pierpaoli C, Basser PJ. The elliptical cone of uncertainty and its normalized measures in diffusion tensor
imaging. IEEE Transactions on Medical Imaging 2008;27(6):834-846.

FIGURE CAPTIONS

29

Figure 1. (A) A generator (in blue) and its surrounding generators (in red) within a spherical cap
(blue closed curve) of predetermined radius. (B) The smallest convex region formed by the
surrounding generators that contains the generator. (C) Finding the smallest convex region that
contains the generator is equivalent to finding the convex hull of the surrounding generators in
the stereographic projection.

30

Figure 2. (A) Percent of relative error of initial generators obtained from the analytically exact
spiral scheme and of the final generators or centroids obtained from the proposed method
compared to the state-of-the-art method in terms of modified electrostatic energies of N
antipodally symmetric point sets with N ranges from 16 to 356. (B) The performance of the
proposed method in terms of execution time (in seconds) per iteration of the proposed method
as a function of the number of points on the upper hemisphere has a linear trend with respect to
N. (C) The frequency of invocation of Step 1, which is an O(N2) algorithm, as a function of N,
number of points on the upper hemisphere. (D) Total number of iterations as a function of N with
instances of invocation of Step1 throughout the iteration are shown in red. The highest peak is
at (N=324, 4161 iterations).

31

Figure 3. (A) 888 deterministic generators taken from the analytically exact spiral scheme and
their Voronoi regions on the upper hemisphere. 888 final (B) real and (C) virtual generators
from the proposed method. (D) 1776 antipodally symmetric points on the sphere obtained from
(B) and (C).

32

Figure 4. Pseudometrically constrained Voronoi tessellations of the sphere with 16000
generators (not shown in the figure) on the upper hemisphere.

33

Figure 5. Reconstructed images of the analytical phantom generated from different sampling
schemes at various slice locations. At the standard resolution of 128x128x128 with z = -0.181,
the images generated by the method of Wong and Roos and the generalized spiral scheme
have more noticeable ringing artifacts around the bright region than that of the proposed
method. The image generated by the tessellated icosahedral scheme has more streak artifacts
in the dark regions than that of the proposed method. At the same resolution with z = 0.228, the
images generated from these methods have more ringing artifacts than that of the proposed
method. Similar patterns of artifacts showed up more noticeably at the high resolution of
256x256x256. At this resolution, the images generated from the proposed method have less
ringing artifacts than those from other methods.

