# sv-lncs - Bio-Grid - University of Connecticut

Software and s/w Development

Dec 1, 2013 (4 years and 7 months ago)

111 views

Parallel Computation of the Minimum Separation
Distance of B
é
zier Curves and Surfaces

Lauren Bissett
1
, Nicholas Woodfield
1

1

University of Connecticut

Department of Computer Science and Engineering

Storrs CT, 06269

{lauren.bissett, nicholas.woodfield}@uc
onn.edu

Abstract.

The detection of the minimum separation distance between any two
points on a geometric surface is a tedious and time
-
consuming process.
Algorithms exist to determine this distance, but they are slow and not practical
for use. However, the
se algorithms can be easily parallelized to improve
performance. By using modern parallel programming architectures, the
speedups for the detection of this distance make it fast enough for practical use.

Keywords:
B
é
zier surfaces, B
é
zier curves, temporal
aliasing, minimum
separation distance, double normals, CUDA

1 Introduction

The underlying principle behind computer animation is the continuous movement
between image frames to give the illusion of motion. When an unexpected change
occurs between frames,

it disrupts the fluidity of motion for the viewer. This problem
is known as temporal aliasing.

One method of preventing temporal aliasing is the detection of imminent or
existing self
-
intersections in geometric surfaces. This can be accomplished by
calcul
ating the minimum separation distance, or the smallest length of separation of
any two points. This distance, if sufficiently small, will indicate the presence of a self
-
intersection.

1.1 Minimum Separation Distance

Although the minimum separation dista
nce is roughly defined as the smallest segment
between any two points on a geometric surface, there are some subtleties. A pair of
points can only be considered a potential segment if they are
doubly normal
.
A
doubly normal segment is such that both its en
dpoints are normal to the geometric
object. Therefore, the minimum separation distance is defined as the the minimum of
all lengths of all double normal segments.

1.2 Calculation of the Minimum Separation Distance

The general algorithm to find the minimu
m separation distance consists of three main
parts:

Generate candidate doubly normal segments

Refine results using Newton's method

Find the segment of minimum length

The geometric object is first searched for potential double normal segments, by
exhaus
tively testing pairs of points. This results in a number of potential estimates for
Newton's Method. Once the candidates are iterated through Newton's Method, their
lengths are compared, and the smallest is considered the minimum separation
distance.

1.3

Geometry

The geometry chosen for this research were B
é
zier curves and B
é
zier surfaces. The
algorithms can also handle composite curves and surfaces.

1.4 Parallelization

This algorithm is ideal for parallelization; the exhaustive checking of large sample

sizes of points in finding candidate doubly normal segments is time
-
consuming in
procedural implementations. Parallelization of this step would decrease run times
exponentially.

2 B
é
zier Curves

2.1 Generation of Candidate Seed Pairs

The generation of s
eed pairs on a given curve is an important step, as good initial
estimates for Newton's method are ideal for finding the correct result.

A pair of endpoints will be considered doubly normal if they satisfy the following
equations:

[
c
s

c
t
]

c'
s
=
0
[
c
s

c
t
]

c'
t
=
0

(
1
)

Where s,t
Є

[0,1]

Initial Estimate Algorithm for Newton's Method.

To produce initial estimates for
Newton's method,
it is first necessary to define a sample size (
n
). This sample size
determines the number of sampl
e points to take along the curve.

The following algorithm generates the initial estimates; the doubly normal test
referenced in the code is detailed below.

Initial Estimate Algorithm For Newton's Method [2]

Input: n, and sub
-
curves ci , where i = 1, .
. . , m

1. For k = 1, . . . , n − 1

s = 1/n + (k
-

1)/n.

2. For = 1, . . . , n − 1

t = 1/n + (
-

1)/n.

3. For i = 1, . . . , m

4. For j = 1, . . . , m

5. If i = j and s = t, cull versus doubly normal test.

6.

If i != j, cull versus doubly normal test.

Output: all segments not culled.

Doubly Normal Test.
The initial estimate algorithm generates pair values, but these
points must then be tested to determine if they are doubly normal, and if they are a
double nor
mal segment.

A key component of this double normal test is the user
-
defined tolerance values, e
1

and e
2
. These tolerances are used to determine if the segment is sufficiently close to
being normal to the curve. It is therefore important to chose values
that will not result
in too many or too few candidates found.

The test, given the two initial points p
1

and p
2
, checks if the endpoints are
sufficiently normal if:

p
1,
p
2

c'
s
<e
1

p
1,
p
2

c'
t
<e
2

(
2
)

If these two inequal
ities are both true, the segment is considered doubly normal and
will be used as an initial estimate for Newton's method.

2.2 Newton's Method

The initial estimates (samples s and t) are then run through Newton's method for
refinement. Newton's method is ef
ficient and produces decent results.

The following equations compute Newton's method for B
é
zier curves:

[
s
n+
1
t
n+
1
]
=
[
s
n
t
n
]

J

1
s
n
,t
n
[
f
s
n
,t
n
g
s
n
,t
n
]
n=
0,1,2,
.
.
.

(
3
)

until the inverse Jacobian matrix (J
-
1
)

is less than some user
-
defined threshold va
lue
(e), Where f and g are the left
-
hand side of the functions in equation 2.1. The
following is the calculation of the Jacobian matrix:

J
s
n
,t
n
=
[

f

s

f

t

g

s

g

t
]

(
4
)

J

1
=
1

f

s

g

t

f

t

g

s
[

g

t

f

t

g

s

f

s
]

(
5
)

And:

f

s
=
[
c
s

c
t
]

c''
s
+c'
s

c'
s

g

s
=c'
s

c'
t

f

t
=

c'
t

c

g

t
=
[
c
s

c
t
]

c''
t
+c'
t

c'
t

(
6
)

After Newton's method has been run on all the seed pairs, the segment of minimum
length is found (through a simple search) and this value is considered the minimum
sepa
ration distance.

2.3 Results

With B
é
zier curves, the CUDA implementation had an approximately 200% speedup
compared to the procedural C code.

Table 1.

Results for B
é
zier curves

n

e1

e2

C code

CUDA code

512

0.4

0.4

28s

154.76ms

512

0.2

0.2

26s

154ms

2
00

0.4

0.4

4.3s

46.6ms

10

0.4

0.4

15ms

1.81ms

3 B
é
zier Surfaces

The algorithm used for B
é
zier curves naturally extends to B
é
zier surfaces, merely
introducing two new variables: u and v.

A pair of endpoints will be considered doubly normal if they sati
sfy the following
equations:

[
S
u
1,
v
1

S
u
2,
v
2
]

S

u
1
=
0
[
S
u
1,
v
1

S
u
2,
v
2
]

S

v
1
=
0
[
S
u
1,
v
1

S
u
2,
v
2
]

S

u
2
=
0
[
S
u
1,
v
1

S
u
2,
v
2
]

S

v
2
=
0

(
7
)

3.1 Generation of Candidate Seed Pairs

Initial Estimate Algorithm for Newton's Method.
To produce initial estimates for
Newton's method with surfaces, again def
ining a sample size (
n
) is necessary,
determining the number of sample points to take along the surface.

Initial Estimate Algorithm for Newton's Method

Input:
n, and sub
-
meshes ci , where i = 1, . . . , m

1. for i = 1 to n

1

u
1

= i / n

2. for j = 1

to n

1

v
1

= j / n

3. for k = 1 to n

1

u
2

= k / n

4. for l = 1 to n

1

v
2

= l / n

5. For i = 1, . . . , m

6. For j = 1, . . . , m

7. If i = j and u
1

= u
2

&& v
1

= v
2
, cull v
ersus doubly normal test.

8. If i != j, cull versus doubly normal test.

Output: All segments not culled

Double Normal Test.

The double normal test is similar to the curve case, in which
each u, v endpoint is determined if sufficiently dou
ble normal, within a certain
threshold.

[
S
u
1,
v
1

S
u
2,
v
2
]

S

u
1
<e
1
[
S
u
1,
v
1

S
u
2,
v
2
]

S

v
1
<e
2
[
S
u
1,
v
1

S
u
2,
v
2
]

S

u
2
<e
3
[
S
u
1,
v
1

S
u
2,
v
2
]

S

v
2
<e
4

(
8
)

3.2

Newton's Method

The initial estimate results (samples u
1
, u
2
, v
1
, v
2
) are then run through the bivariate
Newton's method. The following equations co
mpute Newton's method for B
ézier
surfaces:

[
u
1
n+
1
v
1
n+
1
u
2
n+
1
v
2
n+
1
]
=
[
u
1
n
v
1
n
u
2
n
v
2
n
]

J

1
u
1
n
,v
1
n
,u
2
n
,v
2
n
F
u
1
n
,v
1
n
,u
2
n
,v
2
n

(
9
)

Where:

F=
[
F
1
F
2
G
1
G
2
]

(
10
)

And:

F
1
=
[
S
u
1,
v
1

S
u
2,
v
2
]

S

u
1
=
0
F
2
=
[
S
u
1,
v
1

S
u
2,
v
2
]

S

v
1
=
0
G
1
=
[
S
u
1,
v
1

S
u
2,
v
2
]

S

u
2
=
0
G
2
=
[
S
u
1,
v
1

S
u
2,
v
2
]

S

v
2
=
0

(
11
)

And where:

J
u
1,
v
1,
u
2,
v
2
=
[

F
1

u
1

F
1

v
1

F
1

u
2

F
1

v
2

F
2

u
1

F
2

v
1

F
2

u
2

F
2

v
2

G
1

u
1

G
1

v
1

G
1

u
2

G
1

v
2

G
2

u
1

G
2

v
1

G
2

u
2

G
2

v
2
]

(
12
)

The sixteen partial derivatives are listed in the appendix.

4 Parallelization with CUDA

One of the main goals of this research was to parallelize the algorithms used to find
the minimum separation distance. Luckily, the algorithms are easily adapted to make
use of parallel programming. The parallel architecture chosen was NVIDIA's CUDA
programming language, which allows parallel programming on the graphical
processing unit (GP
U).

4.1 CUDA Background

NVIDIA developed CUDA as a parallel programming language that allows easy
development on the GPU for the programmer, allowing “supercomputing for the
masses” [7]. Normally GPU programming is very difficult, but CUDA extends
elegantl
y as part of C. CUDA
-
enabled GPU's support vastly multi
-
applications, allowing the launch of thousands of simultaneous threads of execution
[3]. It also is a cost
-
effective method of fast, parallel computing without investing
enormous sums of mone
y in supercomputers. The only downside to CUDA is that is
hardware
-
specific, only running on NVIDIA, CUDA
-
enabled GPU's.

The CUDA architecture has three levels of organization: threads, blocks and the
grid. Threads are the basic component of execution,
and have their own register and
local memory. Threads are then organized into blocks;
block is dependent on the architecture
constant and texture memory. The grid contains all the bloc
ks, and can be 1 or 2
-
dimensional. Register and shared memory are fastest, whereas local, constant, global
and texture memory are slower.

4.2 CUDA Specs

The machine used for the results presented in this paper used a QuadroFX 5800 GPU
on a machine running
Windows Vista.

4.3 Parallelization of B
é
zier Curves

The algorithm for the generation of candidate seed pairs (described in section 2.1)
uses four loops, or O(n
4
) complexity, making run times extremely slow for large
sample sizes.

To reduce the run time fo
r this step, each thread handles one sample point on the
curve, and compares it to every other point on all other curves; if it passes the double
normal test, Newton's method is performed.

The Main Kernel.

The main kernel containing the CUDA code perform
s two steps:
generation of candidate pairs, and Newton's method. Each thread in the kernel

The entire input composite curve, represented as an array of control points

A two
-
dimensional output array to hold the double normal segmen
ts

The sample size

The number of curves

Tolerance values (e, e
1
, e
2
)

As stated, each thread handles one sample point on one curve, and checks for double
normals along its curve as well as all the others in the composite. Therefore a 1D grid
used, arranged into blocks. If the sample size is smaller than the allowed
number of threads per block, each curve is represented as one block of threads.
However, if the sample size if larger, the curve must be split up over multiple blocks.

If a threa
d finds double normal segment(s), it then runs Newton's method on that
pair. If more than one double normal was found in that thread, it compares the newest
one found to the previously found segment. If smaller, it is saved and the old one
thrown out. At
the end of execution, the thread saves the shortest segment in the two
dimensional array, whose size is the number of blocks multiplied by the number of
threads. Once all threads have finished executing, the output array will contain the
shortest double no

Secondary Steps.

Once the main kernel has finished finding all double normals, a
second kernel is launched. The purpose of this kernel is to sort through the returned
double normal segments and find the one of minimum length.

It works by first finding
the minimum length from each block. It returns these values, and then they are
iterated over to find the global minimum.

4.4 Parallelization of B
é
zier Surfaces

Modifying the existing CUDA code from B
é
zier curves to surfaces invol
ves not just
modifying the algorithm and equations used, but also changing the layout of the code.
Instead of a 1D array of threads, now a 2D one is used.

Main Kernel
. The main kernel behaves similarly to the one used for Bezier curves.
g one parametric s, t value, each thread in this kernel handles a u, v

Composite mesh input as an array of control points

Output array that will contain found double normal segments.

Sample size.

Tolerance values.

This setup eliminates the need of the first two loops in the pseudocode listed in
section 3.1, already reducing the time complexity, in addition to CUDA's fast
runtimes.

So, similar to the curve kernel, each thread handles one u,v pair, and checks again
st
every other sub
-
mesh for double normal segments. If one is found, it is run through
Newton's method. If more than one double normal was found by the thread, it
compares the current segment to the previously saved segment, and throws away the
shorter.

If the sample size is less than the number of allowed threads per block, each sub
-
mesh is represented by one block. Otherwise, the meshes are split over multiple
blocks. Once all threads have finished execution, the kernel returns an array
containing the
shortest segment found by each thread.

Secondary Steps
. Once the kernel returns the results, a second kernel is launched to
search through the results, and return a smaller array, containing the shortest segment
from each block. From there, the array is s
earched to find the shortest segment. This
is the minimum separation distance.

Future Work

Although the method of parallelization for Bezier surfaces was developed, it was not
implemented. Therefore future work would include fully implementing the CUDA
cod
e for Bezier surfaces, and performing timing tests to determine if sufficiently fast
for practical use.

Acknowledgments

The authors were partially supported by grant NSF Grant #0755373, and supported by
University of Connecticut Biogrid REU program. Both a
uthors thank Dr. Bing Wang
and Dr. Tom Peters for their support and guidance.

References

1. Antonov, E., Bisceglio, J., Borovikov, I., Noble, N., Peters, T.J.: Convergence of Geometric
Algorithms for Graphics & Animation 8

14

2. Moore, L.F.: Computational
Topology of Spline Curves for Geometric and Molecular
Approximations 66

82

3. NVIDIA CUDA Zone,
http://
www.NVIDIA.com/cuda

4. Mathews, John H.
Numerical Methods For Computer Science, Engineering, And
Mathematics
. 1987 Prentice
-
Hall Inc: 1987. 380
-
382

5.
Miller L., Moore E.L.F., Peters T.J., Russell A.C.: Topological Neighborhoods for Spline
Curves: Practice and Theory

6. AK Peters ,

Real
-
Time Rendering 2
nd

Edition
, LTD July 2002 Thomas Akenine
-
Moller et al.

7. Dr. Dobbs CUDA, Supercomputing for the Masses

Part 1

7, http://www.ddj.com/hpc
-
high
-
performance
-
computing/207200659

8. K.E. Jordan et al.: Modeling time and topology for animation and visualization.

9. E. L. F. Moore, T. J. Peters, and J. A. Roulier. Preserving computatational topology by
subdivision

10. Patrikalakis, Nicholas M. et al.: Shape interrogation for computer aided design and
manufacturing.

Appendix

Partial Derivatives for Bivariate Newton's Method:

F
1

u
1
=

S

u
1

S

u
1
[
S
u
1,
v
1

S
u
2,
v
2
]

2
S

u
1
2

(
13
)

F
1

v
1
=

S

v
1

S

u
1

(
14
)

F
1

u
2
=

S

u
2

S

u
1

(
15
)

F
1

v
2
=

S

v
2

S

u
1

(
16
)

F
2

u
1
=

S

u
1

S

v
1

(
17
)

F
2

v
1
=

S

v
1

S

v
1
[
S
u
1,
v
1

S
u
2,
v
2
]

2
S

v
1
2

(
18
)

F
2

u
2
=

S

u
2

S

v
1

(
19
)

F
2

v
2
=

S

v
2

S

v
1

(
20
)

G
1

u
1
=

S

u
1

S

u
2

(
21
)

G
1

v
1
=

S

v
1

S

u
2

(
22
)

G
1

u
2
=

S

u
2

S

u
2
[
S
u
1,
v
1

S
u
2,
v
2
]

2
S

u
2
2

(
23
)

G
1

v
2
=

S

v
2

S

u
2

(
24
)

G
2

u
1
=

S

u
1

S

v
2

(
25
)

G
2

v
1
=

S

v
1

S

v
2

(
26
)

G
2

u
2
=

S

u
2

S

v
2

(
27
)

G
2

v
2
=

S

v
2

S

v
2
[
S
u
1,
v
1

S
u
2,
v
2
]

2
S

v
2
2

(
28
)