INTERSECTION TEST AND BLOSSOMING PERTURBATION FOR DISK PARAMETRIC CURVES AND BALL PARAMETRIC SURFACES

Errors in curve and surface representation due to inaccuracies in the data are considered and accounted for by introducing disk parametric curves and ball parametric surfaces. Intersection test algorithms and interval extensions using blossoming are discussed for each of the three cases of Bézier curves, tensor product surfaces, and triangular patches. A stability analysis is also performed for each of the three cases. It is shown that under certain restrictions disk Bézier curves and triangular ball Bézier patches are stable with respect to perturbations of the control disks (balls); whereas tensor product ball Bézier surfaces are in general not.


Introduction
The curves and surfaces employed in classical geometric modeling represent the shapes of objects.This representation is an abstraction where it is assumed that the curves have no width and that the surfaces have zero thickness.This idealization does not carry over to the real world.There are restrictions on how precisely a curve can be rendered, there are restrictions on possible tool movements when a physical implementation is machined and there are restrictions on material forming to mention a few of the sources of real world problems in implementing abstract designs.
As an example consider the effects of drift and restriction of movement of a tool when the tool is used to machine an object whose boundary is defined by a mathematical curve or surface.As the tool operates, a small change in local position may result in a global deviation from the intended object.As a result the topology of the object might change, the object cannot be integrated used for its intended purpose and the object might not have the intended mechanical properties.A more extensive discussion of problems with machining accuracy is found in [18] and especially in the references therein.
A second example is the problem of how to handle the modeling of an object by means of curves and surfaces given that imprecise information about the object has been provided.If the curve is defined by control points, then how should these be adjusted to take the imprecise information into account in the best possible manner and at the same time to take into account the offset problem occurring in practical computations and operations?
In both of the above examples, curves of nonzero width, and surfaces having nonzero thickness can be employed.They can also be used to fit double boundary lines and double boundary faces.Further applications of such curves and surfaces can be found in [7,13,19,20].
In this paper previous work by the authors on fitting of parametric curves [10], fitting fat curves [9], and disk Bézier curves [8] is augmented by further results on disk curves, tensor product ball surfaces, and triangular ball patches.For each case, bisection algorithms suitable for intersection tests have been included (Sections 3,6,9).Complexity analyses and implementation details are not included for the algorithms.Interval extensions are, however, defined (Sections 4, 7, 10) as stability analyses are (Sections 5, 8, 11).

Disk and ball arithmetics
Interval analysis is a tool for inclusion and estimation of errors and uncertainties in numerical calculations [1].In one dimension, intervals are closed bounded subsets of the real line R and the set of such intervals is denoted by I(R).An arithmetic can be defined on I(R) and certain algebraic properties are observed (see, e.g., [16]) and an analysis can be developed; see Moore [11].The extension of one-dimensional intervals to higher dimensions in the space I m (R) is usually done componentwise in order to take advantage of the properties in one dimension as far as possible.An m-dimensional interval in I m (R) (also called a box) is therefore described by 2m real parameters (the componentwise lower and upper bounds).
An alternate way of extending intervals to higher dimensions is given by disks in R 2 and balls in R m (m > 2).Denoting the set of nonnegative reals by R + , a ball(disk) in R m is defined to be the set described by m + 1 real parameters.The set of balls(disks) in R m is denoted by D m (R).If B ∈ D m (R), then we also write where r = midB is the center of the ball(disk), and ρ = radB, the radius of the ball (disk). If Q. Lin and J. G. Rokne 3 Real ball(disk) arithmetic was discussed in [6] under the name "Hypernormballe."Further details of this arithmetic are found in [8].

The intersection of disk curves
In many applications it is important to know if two curves intersect or not.As an example consider a piece of cloth which is intended for a garment.The garment outlines are defined by curves on an abstraction of the piece of cloth.If the outlines do not intersect, the cutting process is deemed feasible and when the curves are placed on the cloth it is assumed that the garment pieces can be cut out of the cloth.In reality, there are several sources of error such as tool movement, workpiece movement, and numerical evaluation errors that may invalidate the cutting process; even though the abstract design was feasible.In order to guard against such problems, the intersection test for the curves is normally made with a tolerance.Two curves p 1 (t) and p 2 (t) are therefore considered, and it is assumed that both curves are defined over t ∈ [0,1] without loss of generality.In order to guarantee that the curves do not intersect, we require that the minimum distance between the curves satisfies min where δ * is a tolerance parameter.This can then be rephrased in terms of disk parametric curves.A disk parametric curve consists of a parametric axis curve r(t) and an associated variable radius ρ(t), and it is written as B(t) = r(t),ρ(t) .As t varies, the disk r(t),ρ(t) sweeps out a domain whose inner and outer envelopes are r(t) and r(t) as shown in Figure 3.1.The interval extension of B(t) to a parametric interval T is written as B(T) where it is assumed that the extension is computed as, for example, a centered form (see [11,17]).
The curve intersection problem considered above can then be stated in terms of parametric disk curves B 1 (t) = p 1 (t),δ * , B 2 (t) = p 2 (t),δ * , t ∈ [0,1], and requiring that these disk curves do not intersect.If the disk curves do not intersect, then we are guaranteed that the axis curves do not intersect up to the assumption on the possible errors of representation.
We cast the intersection problem in a slightly more general setting where the two disk curves are allowed to have width as a function of the parameter t.Let therefore , where r i (t) ∈ B i (t) is the axis curve and ρ i (t) is the variable perturbation of the axis curve due to an error, i = 1,2.We will test to see whether or not holds.If the property does not hold, then we cannot guarantee that the nonintersection of the axis curves will remain valid when the error perturbations are taken into account and hence the topology of the curves might be changed by the perturbations.
From the above, we can see that our problem can be viewed as an intersection test between two disk curves.This is shown in Figure 3.2.In the following, we will use circular arithmetic to construct a test algorithm for the intersection of disk curves.
Note that Furthermore for any subintervals S ⊆ , T ⊆ -, it is assumed that a computational procedure exists for finding inclusions for The proposed algorithm is based on an exhaustion method, that is, bisection together with interval extension evaluations are used to progressively delete subintervals of [0,1] which correspond to disk segments guaranteed not to intersect.The algorithm is included here for completeness even though the principle of the algorithm is well known.
Because of finite machine precision, we have to prescribe an admissible error bound EPS which is used to terminate the algorithm when the width of subintervals is less than EPS.
We use lists andto store the subintervals of [0,1] for B 1 (s) and B 2 (t) that must be processed further.
After each bisection, we put the uncertain subintervals that must be tested further at the end of the list.This means that the widest subinterval remaining is always the first on each list.For convenience, we let and for any B 1 ∈ Ꮾ 1 (), B 2 ∈ Ꮾ 2 (-), we define (2) Set S = first item on , T = first item on -.
(5) Bisect S so that then enter S i onto the list as the last item.(7) If w(T) < EPS then go to (11).( 8) Delete T from -. (9) Bisect T so that If the algorithm terminates with either orbeing empty, then the two disk curves B 1 (s), B 2 (t) do not intersect.Conversely if the algorithm terminates with both andbeing nonempty, then B 1 (s) ∩ B 2 (t) = ∅ holds to the precision EPS, that is, the two disk curves B 1 (s), B 2 (t) intersect within that precision.In this case it might be necessary to rerun the algorithm with smaller EPS in order to possibly get a more precise result.

The interval extension of a disk Bézier curve
In this section the effect of errors in curve representation is considered for a Bézier curve.
Let {b i }, i = 1,2,...,n be a series of control points.A planar Bézier curve defined on these control points can be written as where see, for example, [3].
There are two sources of errors in the computation of a Bézier curve.These are perturbations of the control points b i and perturbations of the parameter t.
In order to consider the perturbations of control points b i , we use control disks B i = b i ,ρ i instead of the original b i 's obtaining a disk Bézier curve: This means that a disk Bézier curve can be computed from two scalar Bézier curves forming the center and radius of the disk curve, respectively.In order to consider perturbations of the parameter t, we have to consider the interval extension of a disk curve, and we use a one-dimensional disk parameter T = [t − δ,t + δ] instead of the original t, where δ > 0 is an error radius.
The polar form of the curve [2] is an efficient approach for computing the value of a Bézier curve.This can be extended to the computation of a disk curve via the polar form of the disk curve.If F is the polar form of b(t) and if t is written as where b i = F(0,...,0 j ,1,...,1 i ).
Q. Lin and J. G. Rokne 7 Consider now the extension F(T,...,T) of the polar form F to an interval T ⊆ [0,1].This is a type of nested form which is an improvement over the extension of the Bézier form which is defined by Due to the subdistributivity of interval or disk arithmetic, we have and it is also of lower computational complexity when evaluated in a nested fashion.The de Casteljau recursive process described in, for example, [3, page 45] can be used to compute F(T,...,T): In the above process, T is an interval and B r i is a disk; hence the result B n 0 is a disk, that is, The polar form extension F(T,...,T) is an inclusion for the disk Bézier curve which in turn includes the Bézier curve perturbed with errors.Letting we have (4.10)

Intersection test and blossoming perturbation
We can also view T as a disk so that it can be expressed as t,δ ; thus using (2.4) and where 0 (4.12) When the perturbations ρ i is the control points and the perturbation δ of the evaluation point tends to zero, we would expect that the value of the Bézier curve be obtained.This is indeed the case as can be seen from (4.12) and we have the following theorem for the extension of the polar form. ) This means that F(T,...,T) converges to b(t).

Disk analysis for error generated from blossoming
Blossoming is a technique popularized in [2,14,15] for studying polynomials in computer-aided geometric design by replacing them by simple multilinear functions.A blossoming algorithm is now used to estimate the error propagation in computation of a Bézier curve in order to analyze its computational stability.The stability of polynomials expressed in the Bernstein form, which forms the basis for Bézier curves, was shown to be computationally stable in [4,5].
Q. Lin and J. G. Rokne 9 A general disk curve can be generated by means of the following formula: where be the disk perturbations of α r i , β r i .We get (5.3) Considering only the error radii, we have (5.4) In order to estimate the error bound further, we assume that (5.5) We then have (5.6) This shows that the perturbation error vanishes when ρ 0 and τ → 0.
C 1 continuous at join of B 1 (t) and B 2 (t) Figure 5.1.The error propagation.
We now consider error propagation for a disk Bézier spline.Let the disk Bézier spline be given, where and where the axis lines of B k−1 (t), B k (t) (k = 2,...,m) are C 1 continuous; see Figure 5.1.We now consider how the error perturbation of B 1 (t) will affect the other B k (t) (k = 1).When the error of B 1 (t) is evaluated at the join with B 2 (t), the value ρ (2)  0 will change from the original value to ρ (2)  0 + Δρ 0 which gives rise to an additional error radius (1 − t + δ) n Δρ 0 for B 2 (t).It can be seen that this extra error increases further to (1 + δ) n Δρ 0 at the initial point of B 2 (t) (t = 0) due to floating point computational errors.This extra error Δ 0 is bounded by δ n Δρ 0 (0 < δ < 1) at the terminal of B 2 (t) (t = 1).This means that although the perturbation increases at the initial point, it is damped during the computational process.In a similar manner this error is transmitted to B 3 (t) (t = 1) and so on.At the terminal segment B m (t), the extra error quantity Δρ 0 generated from B 1 (t) results in δ mn Δρ 0 .
As we can see from above, for the disk Bézier curves, the extra perturbation generated from a segment of a curve will be propagated to each of the following segments of the curve along with the control points.This perturbation is attenuated as long as the computation is precise enough, that is, if 0 < δ < 1.The construction of a disk Bézier curve is therefore stable.

The intersection ball surfaces
The results in Sections 3-5 are now first generalized to ball surfaces and then later to triangular ball patches.
Q. Lin and J. G. Rokne 11 The topological relation between two given surfaces may be changed by a perturbation.In practise we therefore often require that the minimum distance between two surfaces p 1 (s,t), p 2 (s,t), s,t ∈ [0,1], satisfies min and hence we need to consider the intersection of two ball surfaces with thickness δ * .In this section we consider a more general version of the problem.Two ball surfaces where each ρ i (s,t) ∈ B i (s,t) has to be understood as the error perturbation applied to the surfaces r i (s,t).We want to check if holds or not.If the relationship does not hold, then r 1 (s,t) and r 2 (s,t) may intersect.
From above we can see that our problem is equivalent to the problem of testing for intersection of two ball surfaces.In the following we will use the ball arithmetic to construct a test algorithm for this intersection computation.
Note that where the notation X (S,T) is used to define the variable on the left side of the equation.Let B i (X) = B i (S,T), S,T ⊆ I, i = 1,2, be interval extensions of the ball surfaces B i (s,t), s,t ∈ I, i = 1,2.The proposed algorithm, along the lines of Algorithm 3.1, is based on the exhaustion method which consists of bisection together with interval extension evaluations which are used to progressively delete rectangular subregions of I 2 which correspond to tensor product ball surface pieces guaranteed not to intersect.
Because of the finite machine precision, we also have to prescribe an error bound EPS used to terminate the process.Given a rectangular subregion X = (S,T) of I 2 as in Figure 6.1, the bisection method consists of subdividing X into two equal rectangles X 1 ∪ X 2 = X at the midpoint max{s − s,t − t} of the longer side of the rectangle.
We use lists ᐄ and ᐅ to store subregions of I 2 that must be processed further.After each bisection, we put the uncertain subregions that must be tested further at the end of their list.This means that the largest subregion remaining is always first on a list.
For convenience, we let and for any The test algorithm is as follows.
(5) Bisect X so that X then enter X i onto the list ᐄ as the last item.(7) If w(Y ) < EPS then go to (11) If the algorithm terminates with ᐄ or ᐅ empty, then the two ball surfaces B 1 (s,t), B 2 (u,v) do not intersect (to the precision EPS).Conversely, if both ᐄ and ᐅ are not empty, then B 1 (s,t) ∩ B 2 (u,v) = ∅ holds to the precision EPS, that is, the two ball surfaces intersect.In this case it is necessary to return to the algorithm with a smaller EPS if a more precise result is required.

The interval extension of a tensor product ball Bézier surface
In this section we consider the particular case of tensor product ball surfaces.
Let {b i j }, i, j = 0,1,...,n be a series of 3D-control points.A tensor product Bézier surface defined on these control points is expressed as As in the case of Bézier curves, we consider errors both due to perturbations of the control points b i j and due to the perturbations of the parameters s and t.
In order to reflect the perturbations of the control points b i j , we use a control ball B i j = b i j ,ρ i j instead of the original b i j .This results in a ball Bézier surface: This means that a ball Bézier surface can be computed from two scalar Bézier surfaces forming the center and radius of the ball surface, respectively.In order to include the perturbation of the parameters s,t, we have to consider the interval extension of a ball surface.We use S = [s − ,s + ] and T = [t − δ,t + δ] instead of the original s,t, where ,δ > 0 are the error radii.
Let F be the polar form of b(s,t), and let s be written as ).
Consider now the interval extension F(S,...,S;T,...,T) of the polar form F over (S,T), S,T ⊆ I.This is a type of nested form which is an improvement of the Bézier form which is defined as In practical computations, we use the following de Casteljau recursive process (see [3]) to obtain F(S,...,S;T,...,T): Output B n,n 0,0 .Return.(7.6)In the above process, S,T are symmetric intervals and B r,r i, j is a ball; hence the result B n,n 0,0 is a ball, that is, B n,n 0,0 = F(S,...,S;T,...,T).Let Since S,T are disks, they can be written as s, , t,δ , and we have When the perturbations of the control points and the error radii , δ of the evaluation point tend to zero, we would expect that the value of the tensor product Bézier surface to be obtained.This is indeed the case as can be seen from ( 7.10), and we have the following theorem for the extension of the polar form.
For the extension to polar form the following theorem follows.

Ball analysis of error generated from blossoming for tensor product surfaces
A general ball surface generated by means of the following recursive formula: where is considered, where the ball blossoming is B n,n 0,0 .The blossoming algorithm is used to account for the error propagation in the computation of a ball Bézier surface.To analyze the errors in the algorithm the following balls are defined: When the above perturbations are entered into (8.1),we get b r,r i, j ,ρ r,r i, .
If we only consider the error radius, we have In order to further estimate the error bound, the quantities in the error radius have to be bounded as follows: The error radius can then be estimated as The error propagation among tensor product ball Bézier surfaces is considered next.Let the ball Bézier surface be given, where Moreover, we require that B k−1,l (s,t) and B k,l (s,t) (∀l) are C 1 continuous on their boundary curve, and that B k,l−1 (s,t) and B k,l (s,t) (∀k) are C 1 continuous on their boundary curve; see also Figure 8.1.
Q. Lin and J. G. Rokne 19 We wish to assess the effect of an error perturbation of B 11 (s,t) on the other B k,l (s,t) (k,l not equal to 1 at the same time).
When the error of B 11 (s,t) is transmitted to the curve linked to B 12 (s,t), the value of ρ (1,2)  0, j changes from its original value to ρ (1,2) 0, j + Δρ 0 , which then gives rise to a change of the global error radius of B 12 (s,t).In this case the change of the error of B 11 (s,t) will add an extra quantity that is no larger than j+q=n n j q to the error radius of B 12 (s,t).As this extra error is transmitted to the two boundary curves (s = 1,t = 1) of B 12 (s,t), it will give rise to an increase in the error radii of those curves bounded by When (s = 1, t = 1), that is, at the control ball center b (1,2)  n,n , the extra error quantity is as long as the computational precision < (1 + 2δ) −1 satisfies < 1/3.It should be noted that when s = 0, t = 1, the extra error quantity at the control ball b (1,2)  0,n is As can be seen from the above, the extra perturbation generated from a ball Bézier surface patch will be propagated to each of the following surface patches along with the control balls.In the general case, if the computation is not precise enough, that is, ,δ are not sufficiently small, we can find a control ball b (1,2)  n,n , where the extra error tends to enlarge at this ball.Therefore the construction of tensor product ball Bézier surfaces may be computationally unstable, that is, the construction for curves is absolutely stable but that for surfaces is not.

The intersection of triangular ball patches
We sometimes have to consider the intersection of triangular ball patches in the applications where triangular patches are defined in, for example, [3].
Letting the parameters (s,t) be the barycentric coordinates with regard to some domain triangle, we require that the minimum distance between two triangular patches p 1 (s,t), p 2 (s,t), s,t ∈ [0,1], satisfies min ∀s,t,u,v∈[0,1] p 1 (s,t) − p 2 (s,t)P ≥ 2δ * . (9.1) We must therefore study the intersection of two triangular ball patches B 1 (s,t), B 2 (s,t), s,t ∈ [0,1], both having thickness δ * .We are given B i (s,t) = r i (s,t), ρ i (s,t) , i = 1,2, s,t ∈ [0,1], where each ρ i (s,t) ∈ B i (s,t) will be understood as a deformation of r i (s,t) due to an error, and we are going to consider whether holds.If not, it is possible that r 1 (s,t) and r 2 (s,t) will intersect because of the error perturbation.
From the above we can see that the problem can be viewed as one of intersecting two triangular ball patches.In the following we will therefore construct a test algorithm for the intersection of two triangular ball patches.
Note that  Let B i (S,T), S,T ⊆ I, i = 1,2, be an interval extension of the triangular ball patches B i (s,t), s,t ∈ I, i = 1,2.The proposed algorithm, which is another generalization of Algorithm 3.1, is based on bisection together with interval extension evaluations which are used to delete subregions of (I,I) which correspond to the two pieces of ball patches guaranteed not to intersect.For a subregion X = (S,T) of (I,I), the bisection method is the same as for tensor product surfaces, that is, X is partitioned into two subregions We use a list ᐄ to store the subregions of (I,I) for B 1 (X) that must be processed further, and use another list ᐅ to store the subregions of (I,I) for B 2 (Y ) that must be processed further.
After each bisection we put the uncertain subregions that must be tested further at the end of their list.Q. Lin and J. G. Rokne 21 For convenience, we let and for any B 1 ∈ Ꮾ 1 (ᐄ), B 2 ∈ Ꮾ 2 (ᐅ), we define The test algorithm is as follows.(2) Set X = first item on ᐄ, Y = first item on ᐅ.
(5) Bisect X so that X

The interval extension of a triangular ball Bézier patch
Let {b i jk }, i + j + k = n, be a set of 3D-control points.A triangular Bézier patch defined by these control points is expressed as where In order to consider the error perturbation of the control points b i jk , we use a control ball B i jk = b i jk ,ρ i jk instead of the original b i jk .This results in a ball Bézier patch In order to further consider the error perturbation of the parameter s,t, we consider the interval extension of a ball patch.We use S = [s − ,s + ] and T = [t − δ,t + δ] instead of the original s and t.Since the parameters (s,t) are the barycentric coordinates with regard to some domain triangle Δ = Δ(λ,μ,θ), a point z on Δ can be represented as Consider the interval extension F Δ (S,...,S,T,...,T) of the polar form F Δ on the (S,T), S,T ⊆ I.This is a type of nested form which is an improvement of the following Bézier extension: In practical computations we use the following de Casteljau recursive process (see [3]) to Q. Lin and J. G. Rokne 23 obtain F Δ (S,...,S,T,...,T): That is, B n 0,0,0 = F Δ (S,...,S,T,...,T).Letting we have (10.10) By considering S,T as disks they can be expressed as s, , t,δ , and we have where it is assumed that 0 ≤ s + , t + δ,1 − s − t + + δ ≤ 1. Therefore F Δ (S,...,S,T ...,T) (10.12) For the extension of polar form the following theorem follows.

Ball analysis of error generated from blossoming for triangular patches
We first use the blossoming algorithm to take account of the error propagation in the computation of a triangular ball Bézier patch so as to analyze its computational stability.
A general triangular ball patch can be generated by means of the following recursive formula: where 2) The ball blossoming B n 0,0,0 is obtained from this computation.Letting  where be given.Moreover, two adjacent patches B v1 (s,t) and B v2 (s,t)(v 1 = v 2 ) are assumed to be C 1 continuous on their common boundary curve.We now wish to discuss how the error perturbation of B 1 (s,t) disrupts the other B l (s,t)(l = 1).
When the error of B 1 (s,t) is transmitted to the curve B 2 (0,t) ⊆ j+k=n n i j k t j (1 − t) k b (2)  0 jk , j+k=n n i j k t j (1 − t) k ρ (2)   0 jk (11.10) which is linked up, B 2 (s,t) the value of ρ (2) 0, j,k will change from the original value to ρ (2) 0, j,k + Δρ 0 which then gives rise to a global change in the error radius of B 2 (s,t).In this case the change in the error radius of B 1 (s,t) will result in a quantity that is no more than j+k=n n i j k (t + δ) j (1 − s − t + + δ) k • Δρ 0 = (1 − s + + 2δ) n • Δρ 0 (11.11)added to the error radius of B 2 (s,t).As this extra error is transmitted to the control point b (2)  n,0,0 (s = 1,t = 0) of B 2 (s,t), this will give rise to changes in the error radius at that point, which are no more than ( + δ) n • Δρ 0 .It is easy to see that we can guarantee ( + δ) n • Δρ 0 < Δρ 0 (11.12) as long as + δ < 1 is valid.
Q. Lin and J. G. Rokne 27 As we can see from above, in the triangular ball Bézier patches, the perturbation generated from a piece of ball patch will be transmitted to each of the following pieces of patches along with the control points.This perturbation is attenuated, however.Therefore the construction of triangular ball Bézier patches is stable.

Conclusion
In this paper disk parametric curves, tensor ball surfaces, and triangular ball patches were considered from the point of including errors in representation.Algorithms for intersection testing were given for each case, and analyses for errors were considered.It was found that disk parametric curves and triangular ball patches were stable under bounded error perturbations, whereas tensor ball surfaces were not.
then enter T j onto the listas the last item.(11)If w(S) or w(T) > EPS then go to (2).(12) End.
then enter Y j onto the list ᐅ as the last item.(11)If w(X) or w(Y ) > EPS then go to (2).(12) End.
at the midpoint max{s − s,t − t} of the interval component with the longer length.For instance, if s − s < t − t, then