DIRECTORY IO, Real, Rope, Matrix3d, Spline3d, Vector3d;

Spline3dImpl: CEDAR PROGRAM
IMPORTS Matrix3d, Real, Vector3d
EXPORTS Spline3d
~ BEGIN
Matrix:				TYPE ~ Matrix3d.Matrix;
MatrixRep:			TYPE ~ Matrix3d.MatrixRep;

Pair:					TYPE ~ Matrix3d.Pair;
Quad:					TYPE ~ Matrix3d.Quad;
Line:					TYPE ~ Vector3d.Line;

Triple:				TYPE ~ Spline3d.Triple;
TripleSequence:		TYPE ~ Spline3d.TripleSequence;
TripleSequenceRec:	TYPE ~ Spline3d.TripleSequenceRec;

RealSequence:		TYPE ~ Spline3d.RealSequence;
RealSequenceRec:	TYPE ~ Spline3d.RealSequenceRec;

KnotSequence:		TYPE ~ Spline3d.KnotSequence;
KnotSequenceRep:	TYPE ~ Spline3d.KnotSequenceRep;

Bezier:				TYPE ~ Spline3d.Bezier;
Bspline:				TYPE ~ Spline3d.Bspline;
Hermite:				TYPE ~ Spline3d.Hermite; 

Coeffs:				TYPE ~ Spline3d.Coeffs; 
CoeffsRep:			TYPE ~ Spline3d.CoeffsRep;
CoeffsSequence:		TYPE ~ Spline3d.CoeffsSequence;
CoeffsSequenceRec:	TYPE ~ Spline3d.CoeffsSequenceRec;

TPosSequence:		TYPE ~ Spline3d.TPosSequence;
TPosSequenceRec:	TYPE ~ Spline3d.TPosSequenceRec;
InterpolateCyclic: PUBLIC PROC [knots: KnotSequence, tension: REAL _ 1.0]
RETURNS [CoeffsSequence] ~ {
coeffs:				CoeffsSequence;
nCoeffs:				NAT _ MAX[1, knots.length-1];	-- number of curve intervals
lastK:					NAT _ nCoeffs-1;					-- max for interval index
h, a, b, c, d, r, s:	RealSequence;

IF knots.length < 1 THEN RETURN[NIL];

IF coeffs=NIL OR coeffs.maxLength<nCoeffs THEN coeffs _ NEW[CoeffsSequenceRec[nCoeffs]];

IF knots.length < 3 THEN {
[coeffs[0][3][0], coeffs[0][3][1], coeffs[0][3][2]] _ knots[0];
FOR c: NAT IN[0..2] DO FOR r: NAT IN[0..3] DO coeffs[0][c][r] _ 0.0; ENDLOOP; ENDLOOP;
IF knots.length = 2 THEN {
coeffs[0][2][0] _ knots[1].x-knots[0].x;
coeffs[0][2][1] _ knots[1].y-knots[0].y;
coeffs[0][2][2] _ knots[1].z-knots[0].z;
};
RETURN[coeffs];
};

h _ NEW[RealSequenceRec[nCoeffs]];
a _ NEW[RealSequenceRec[nCoeffs]];
b _ NEW[RealSequenceRec[nCoeffs]];
d _ NEW[RealSequenceRec[nCoeffs]];
c _ NEW[RealSequenceRec[nCoeffs]];
r _ NEW[RealSequenceRec[nCoeffs]];
s _ NEW[RealSequenceRec[nCoeffs]];

FOR i: NAT IN [0..nCoeffs) DO
[coeffs[i][3][0], coeffs[i][3][1], coeffs[i][3][2]] _ knots[i];
ENDLOOP;

h _ ComputeChords[knots, h];

a[0] _ 2.0*(h[0]+h[lastK]);
FOR i: NAT IN [1..lastK) DO a[i] _ 2.0*(h[i-1]+h[i])-h[i-1]*h[i-1]/a[i-1]; ENDLOOP;
c[0] _ h[lastK];
FOR i: NAT IN [1..lastK) DO c[i] _ -h[i-1]*c[i-1]/a[i-1]; ENDLOOP;

FOR i: NAT IN [0..3) DO
p2: REAL;
FOR k: NAT IN [0..lastK] DO
d[k] _ SELECT i FROM
0			=> (knots[k+1].x-knots[k].x)/h[k],
1			=> (knots[k+1].y-knots[k].y)/h[k],
ENDCASE	=> (knots[k+1].z-knots[k].z)/h[k];
ENDLOOP;

b[0] _ 6.0*(d[0]-d[lastK]);
FOR k: NAT IN [1..lastK) DO b[k] _ 6.0*(d[k]-d[k-1]); ENDLOOP;
FOR k: NAT IN [1..lastK) DO b[k] _ b[k]-h[k-1]*b[k-1]/a[k-1]; ENDLOOP;

r[lastK] _ 1.0;
s[lastK] _ 0.0;
FOR k: NAT DECREASING IN [0..lastK) DO
r[k] _ -(h[k]*r[k+1]+c[k])/a[k];
s[k] _ (b[k]-h[k]*s[k+1])/a[k];
ENDLOOP;

p2 _ 6.0*(d[lastK]-d[lastK-1])-h[lastK]*s[0]-h[lastK-1]*s[lastK-1];
p2 _ p2/(h[lastK]*r[0]+h[lastK-1]*r[lastK-1]+2.0*(h[lastK]+h[lastK-1]));
coeffs[lastK][1][i] _ 0.5*p2;
FOR k: NAT IN [0..lastK) DO coeffs[k][1][i] _ 0.5*(r[k]*p2+s[k]); ENDLOOP;

FOR k: NAT IN [0..lastK) DO
coeffs[k][0][i] _ (coeffs[k+1][1][i]-coeffs[k][1][i])/3.0;
coeffs[k][2][i] _ d[k]-h[k]*(2.0*coeffs[k][1][i]+coeffs[k+1][1][i])/3.0;
ENDLOOP;
coeffs[lastK][0][i] _ (coeffs[0][1][i]-coeffs[lastK][1][i])/3.0;
coeffs[lastK][2][i] _ d[lastK]-h[lastK]*(2.0*coeffs[lastK][1][i]+coeffs[0][1][i])/3.0;

FOR k: NAT IN [0..lastK] DO
coeffs[k][0][i] _ h[k]*h[k]*coeffs[k][0][i];
coeffs[k][1][i] _ h[k]*h[k]*coeffs[k][1][i];
coeffs[k][2][i] _ h[k]*coeffs[k][2][i];
ENDLOOP;

ENDLOOP;

RETURN[coeffs];
};

Interpolate: PUBLIC PROC [knots: KnotSequence, tan0, tan1: Triple _ [0.0, 0.0, 0.0], tension: REAL _ 1.0, coeffs: CoeffsSequence _ NIL] RETURNS [CoeffsSequence] ~ {
SELECT knots.length FROM
0 => RETURN[coeffs];
1, 2 => {
len: REAL _ Vector3d.Distance[knots[0], knots[1]];
IF coeffs = NIL OR coeffs.maxLength < 1 THEN coeffs _ NEW[CoeffsSequenceRec[1]];
coeffs[0] _ CoeffsFromHermite[
[knots[0],
 knots[1],
 Vector3d.Mul[tan0, len],
 Vector3d.Mul[tan1, len]],
coeffs[0]];
RETURN[coeffs];
};
ENDCASE => {
chords: RealSequence _ ComputeChords[knots, NIL];
tangents: TripleSequence _ ComputeTangents[knots, chords, NIL, tan0, tan1];
RETURN[ComputeCoeffs[knots, tangents, tension, chords, coeffs]];
};
};

ComputeChords: PROC [knots: KnotSequence, chords: RealSequence _ NIL]
RETURNS [RealSequence] ~ {
IF chords = NIL OR chords.maxLength < knots.length
THEN chords _ NEW[RealSequenceRec[knots.length]];
FOR n: NAT IN[0..knots.length-1) DO
IF (chords[n] _ Vector3d.Distance[knots[n+1], knots[n]]) = 0.0 THEN chords[n] _ 1.0;
ENDLOOP;
RETURN[chords];
};

ComputeTangents: PROC [knots: KnotSequence, chords: RealSequence, tangents: TripleSequence _ NIL, tan0, tan1: Triple] RETURNS [TripleSequence] ~ {
N: TripleSequence _ NEW[TripleSequenceRec[knots.length]];		-- nonzero elements of M
B: TripleSequence _ NEW[TripleSequenceRec[knots.length]];		-- a control matrix
IF tangents = NIL OR tangents.maxLength < knots.length
THEN tangents _ NEW[TripleSequenceRec[knots.length]];
SetEndConditions[knots, N, B, tangents, chords, tan0, tan1];
GaussianEliminate[knots, N, B, tangents, chords];
RETURN[tangents];
};

SetEndConditions: PROC [knots: KnotSequence, N, B, tangents: TripleSequence, chords: RealSequence, tan0, tan1: Triple] ~ {
z: NAT _ knots.length-1;
IF Vector3d.Null[tan0] THEN {
N[0] _ [0.0, 1.0, 0.5];
B[0] _ Vector3d.Mul[Vector3d.Sub[knots[1], knots[0]], 1.5/chords[0]];
}
ELSE {
tangents[0] _ B[0] _ tan0;
N[0] _ [0.0, 1.0, 0.0];
};
IF Vector3d.Null[tan1] THEN {
N[z] _ [2.0, 4.0, 0.0];
B[z] _ Vector3d.Mul[Vector3d.Sub[knots[z], knots[z-1]], 6.0/chords[z-1]];
}
ELSE {
tangents[z] _ B[z] _ tan1;
N[z] _ [0.0, 1.0, 0.0];
};
};

GaussianEliminate: PROC [knots: KnotSequence, N, B, tangents: TripleSequence, chords: RealSequence] ~ {
nPts: NAT _ knots.length;

FOR n: NAT IN[1..nPts-1) DO
l0: REAL _ chords[n-1];
l1: REAL _ chords[n];
N[n] _ [l1, l0+l0+l1+l1, l0];
B[n] _ Vector3d.Mul[
Vector3d.Add[
Vector3d.Mul[Vector3d.Sub[knots[n+1], knots[n]], l0*l0],
Vector3d.Mul[Vector3d.Sub[knots[n], knots[n-1]], l1*l1]],
3.0/(l0*l1)];
ENDLOOP;

FOR n: NAT IN[1..nPts) DO											-- gaussian elimination
d, q: REAL;
IF N[n].x = 0.0 THEN LOOP;
d _ N[n-1].y/N[n].x;
N[n] _ Vector3d.Sub[Vector3d.Mul[N[n], d], N[n-1]];
B[n] _ Vector3d.Sub[Vector3d.Mul[B[n], d], B[n-1]];
q _ 1.0/N[n].y;
N[n] _ Vector3d.Mul[N[n], q];
B[n] _ Vector3d.Mul[B[n], q];
ENDLOOP;

tangents[nPts-1] _ Vector3d.Div[B[nPts-1], N[nPts-1].y];		-- back substitution
FOR n: NAT IN[2..nPts] DO
tangents[nPts-n] _
Vector3d.Div[
Vector3d.Sub[
B[nPts-n],
Vector3d.Mul[tangents[nPts-n+1], N[nPts-n].z]],
N[nPts-n].y];
ENDLOOP;
};

ComputeCoeffs: PROC [knots: KnotSequence, tangents: TripleSequence, tension: REAL, chords: RealSequence, coeffs: CoeffsSequence _ NIL] RETURNS [CoeffsSequence] ~ {
nCoeffs: NAT _ knots.length-1;

IF coeffs = NIL OR coeffs.maxLength < nCoeffs
THEN coeffs _ NEW[CoeffsSequenceRec[nCoeffs]];

IF tension # 1.0 THEN FOR n: NAT IN[0..nCoeffs] DO
tangents[n] _ Vector3d.Mul[tangents[n], tension];
ENDLOOP;
FOR m: NAT IN[0..nCoeffs) DO
l: REAL _ chords[m];
dknots: Triple _ Vector3d.Sub[knots[m+1], knots[m]];
a: Triple _ Vector3d.Mul[tangents[m], l];
b: Triple _ Vector3d.Add[Vector3d.Mul[tangents[m+1], l], a];
c: Triple _ Vector3d.Add[Vector3d.Mul[dknots, -2.0], b];
d: Triple _ Vector3d.Sub[Vector3d.Sub[Vector3d.Mul[dknots, 3.0], b], a];
IF coeffs[m] = NIL THEN coeffs[m] _ NEW[CoeffsRep];
coeffs[m][0] _ [c.x, c.y, c.z, 0.0];
coeffs[m][1] _ [d.x, d.y, d.z, 0.0];
coeffs[m][2] _ [a.x, a.y, a.z, 0.0];
coeffs[m][3] _ [knots[m].x, knots[m].y, knots[m].z, 1.0];
ENDLOOP;

RETURN[coeffs];
};
hermiteToCoeffs: Matrix _ NEW[MatrixRep _ [
			[ 2.0,	-2.0,	 1.0,	 1.0],
			[-3.0,	 3.0,	-2.0,	-1.0],
			[ 0.0,	 0.0,	 1.0,	 0.0],
			[ 1.0,	 0.0,	 0.0,	 0.0]]];

coeffsToHermite: Matrix _ NEW[MatrixRep _ [
			[0.0,	 0.0,	 0.0,	 1.0],
			[1.0,	 1.0,	 1.0,	 1.0],
			[0.0,	 0.0,	 1.0,	 0.0],
			[3.0,	 2.0,	 1.0,	 0.0]]];

bezierToCoeffs: Matrix _ NEW[MatrixRep _ [
			[-1.0,	 3.0,	-3.0,	 1.0],
			[ 3.0,	-6.0,	 3.0,	 0.0],
			[-3.0,	 3.0,	 0.0,	 0.0],
			[ 1.0,	 0.0,	 0.0,	 0.0]]];

coeffsToBezier: Matrix _ NEW[MatrixRep _ [
			[0.0,	 0.0,		 0.0,		 1.0],
			[0.0,	 0.0,	 	1.0/3.0,	 1.0],
			[0.0,	 1.0/3.0,	 2.0/3.0,	 1.0],
			[1.0,	 1.0,	 	1.0,		 1.0]]];

bsplineToCoeffs: Matrix _ NEW[MatrixRep _ [
			[-1.0/6.0,	 3.0/6.0,	-3.0/6.0,	 1.0/6.0],
			[ 3.0/6.0,	-6.0/6.0,	 3.0/6.0,	 0.0/6.0],
			[-3.0/6.0,	 0.0/6.0,	 3.0/6.0,	 0.0/6.0],
			[ 1.0/6.0,	 4.0/6.0,	 1.0/6.0,	 0.0/6.0]]];

coeffsToBspline: Matrix _ NEW[MatrixRep _ [
			[0.0,	  2.0/3.0,	-1.0,	 1.0],
			[0.0,	-1.0/3.0,	 0.0,	 1.0],
			[0.0,	  2.0/3.0,	 1.0,	 1.0],
			[1.0,	 11.0/3.0,	 2.0,	 1.0]]];

ConvertToCoeffs: PROC [p0, p1, p2, p3: Triple, convert, out: Matrix, hermite: BOOL _ FALSE]
RETURNS [Coeffs] ~ {
w: REAL _ IF hermite THEN 0.0 ELSE 1.0;
m: Matrix _ Matrix3d.ObtainMatrix[];
m^ _ [
[p0.x,	p0.y,	p0.z,	1.0],
[p1.x,	p1.y,	p1.z,	1.0],
[p2.x,	p2.y,	p2.z,	w],
[p3.x,	p3.y,	p3.z,	w]];
IF out = NIL THEN out _ NEW[CoeffsRep];
[] _ Matrix3d.Mul[convert, m, out];
Matrix3d.ReleaseMatrix[m];
RETURN[out];
};

ConvertFromCoeffs: PROC [c: Coeffs, convert: Matrix] RETURNS [p0, p1, p2, p3: Triple] ~ {
m: Matrix _ Matrix3d.ObtainMatrix[];
[] _ Matrix3d.Mul[convert, c, m];
p0 _ [m[0][0], m[0][1], m[0][2]];
p1 _ [m[1][0], m[1][1], m[1][2]];
p2 _ [m[2][0], m[2][1], m[2][2]];
p3 _ [m[3][0], m[3][1], m[3][2]];
Matrix3d.ReleaseMatrix[m];
RETURN[p0, p1, p2, p3];
};

CoeffsFromBezier: PUBLIC PROC [b: Bezier, out: Coeffs _ NIL] RETURNS [Coeffs] ~ {
RETURN[ConvertToCoeffs[b.b0, b.b1, b.b2, b.b3, bezierToCoeffs, out]];
};
CoeffsFromBspline: PUBLIC PROC [b: Bspline, out: Coeffs _ NIL] RETURNS [Coeffs] ~ {
RETURN[ConvertToCoeffs[b.b0, b.b1, b.b2, b.b3, bsplineToCoeffs, out]];
};
CoeffsFromHermite: PUBLIC PROC [h: Hermite, out: Coeffs _ NIL] RETURNS [Coeffs] ~ {
RETURN[ConvertToCoeffs[h.p0, h.p1, h.tan0, h.tan1, hermiteToCoeffs, out, TRUE]];
};

BezierFromCoeffs: PUBLIC PROC [c: Coeffs] RETURNS [Bezier] ~ {
b0, b1, b2, b3: Triple;
[b0, b1, b2, b3] _ ConvertFromCoeffs[c, coeffsToBezier];
RETURN[[b0, b1, b2, b3]];
};
BsplineFromCoeffs: PUBLIC PROC [c: Coeffs] RETURNS [Bspline] ~ {
b0, b1, b2, b3: Triple;
[b0, b1, b2, b3] _ ConvertFromCoeffs[c, coeffsToBspline];
RETURN[[b0, b1, b2, b3]];
};
HermiteFromCoeffs: PUBLIC PROC [c: Coeffs] RETURNS [Hermite] ~ {
p0, p1, tan0, tan1: Triple;
[p0, p1, tan0, tan1] _ ConvertFromCoeffs[c, coeffsToHermite];
RETURN[[p0, p1, tan0, tan1]];
};
PerPointProc: TYPE = Spline3d.PerPointProc;

WalkBezier: PUBLIC PROC [b: Bezier, proc: PerPointProc, epsilon: REAL _ 0.05, doLast: BOOL _ TRUE] ~ {
Inner: PROC [b: Bezier] ~ {
IF FlatBezier[b, epsilon]
THEN proc[b.b0]
ELSE {
left, rite: Bezier;
[left, rite] _ SplitBezier[b];
Inner[left];
Inner[rite];
};
};
Inner[b];
IF doLast THEN proc[b.b3];
};
FwdDif: PUBLIC PROC [in: Coeffs, nSegments: INTEGER, out: Coeffs_NIL] RETURNS [Coeffs] ~ {
fwdDif: Matrix3d.Matrix _ NEW[Matrix3d.MatrixRep];
d1, d2, d3, d6: REAL;
IF nSegments = 0 THEN RETURN[NIL];
d1 _ 1.0/nSegments;
d2 _ d1*d1;
d3 _ d1*d2;
d6 _ 6.0*d3;
fwdDif^ _ [
		[0.0,	0.0,		0.0,	1.0],
		[d3,	d2,		d1,	0.0],
		[d6,	d2+d2,	0.0,	0.0],
		[d6,	0.0,		0.0,	0.0]];
RETURN[Matrix3d.Mul[fwdDif, in, out]];
};

Samples: PUBLIC PROC [c: Coeffs, nPoints: INTEGER, points: TripleSequence _ NIL] RETURNS [TripleSequence] ~ {
p: ARRAY [0..2] OF REAL;
fwdDif: Coeffs _ FwdDif[c, nPoints-1];
IF points = NIL OR points.maxLength < nPoints THEN
points _ NEW[TripleSequenceRec[nPoints]];
points[0] _ [fwdDif[0][0], fwdDif[0][1], fwdDif[0][2]];
FOR i: INTEGER IN[0..2] DO p[i] _ fwdDif[0][i]; ENDLOOP;
FOR i: INTEGER IN[1..nPoints) DO
FOR j: INTEGER IN[0..2] DO
p[j] _ p[j]+fwdDif[1][j];
fwdDif[1][j] _ fwdDif[1][j]+fwdDif[2][j];
fwdDif[2][j] _ fwdDif[2][j]+fwdDif[3][j];
ENDLOOP;
points[i] _ [p[0], p[1], p[2]];
ENDLOOP;
RETURN[points];
};

Position: PUBLIC PROC [c: Coeffs, t: REAL] RETURNS [Triple] ~ {
ret: ARRAY [0..2] OF REAL;
IF t = 0.0 THEN FOR i: NAT IN[0..2] DO ret[i] _ c[3][i]; ENDLOOP
ELSE IF t = 1.0 THEN FOR i: NAT IN[0..2] DO ret[i] _ c[0][i]+c[1][i]+c[2][i]+c[3][i]; ENDLOOP
ELSE {
t2: REAL _ t*t;
t3: REAL _ t*t2;
FOR i: NAT IN[0..2] DO ret[i] _ t3*c[0][i]+t2*c[1][i]+t*c[2][i]+c[3][i]; ENDLOOP
};
RETURN[[ret[0], ret[1], ret[2]]];
};

Velocity: PUBLIC PROC [c: Coeffs, t: REAL] RETURNS [Triple] ~ {
tan: ARRAY [0..2] OF REAL;
a: REAL _ 3.0;
b: REAL _ 2.0;
IF t = 0.0 THEN RETURN[[c[2][0], c[2][1], c[2][2]]];
IF t # 1.0 THEN {a _ a*t*t; b _ b*t};
FOR i: NAT IN[0..2] DO tan[i] _ a*c[0][i]+b*c[1][i]+c[2][i]; ENDLOOP;
RETURN[[tan[0], tan[1], tan[2]]];
};
Acceleration: PUBLIC PROC [c: Coeffs, t: REAL] RETURNS [Triple] ~ {
a: REAL = 6.0*t;
RETURN[[a*c[0][0]+c[1][0]+c[1][0], a*c[0][1]+c[1][1]+c[1][1], a*c[0][2]+c[1][2]+c[1][2]]];
};
Tangent:  PUBLIC PROC [c: Coeffs, t: REAL] RETURNS [Triple] ~  {
RETURN[Velocity[c, t]];
};

MinAcceleration: PUBLIC PROC [c: Spline3d.Coeffs] RETURNS [t: REAL] ~ {
a0: Triple _ Acceleration[c, 0.0];
axis: Triple _ Vector3d.Sub[Acceleration[c, 1.0], a0];
t _ IF Vector3d.Null[axis] THEN 0.5 ELSE -Vector3d.Dot[a0, axis]/Vector3d.Square[axis];
IF t < 0.0 THEN t _ 0.0 ELSE IF t > 1.0 THEN t _ 1.0;
};

CurvatureMag: PUBLIC PROC [c: Coeffs, t: REAL] RETURNS [REAL] ~ {
tan: Triple _ Tangent[c, t];
acc: Triple _ Acceleration[c, t];
m: REAL _ Vector3d.Mag[tan];
RETURN[IF m = 0.0 THEN 0.0 ELSE Vector3d.Mag[Vector3d.Cross[acc, tan]]/(m*m*m)];
};

Curvature: PUBLIC PROC [c: Coeffs, t: REAL] RETURNS [Triple] ~ {
tan: Triple _ Tangent[c, t];
IF tan = [0.0, 0.0, 0.0] THEN RETURN[tan] ELSE {
acc: Triple _ Acceleration[c, t];
mag: REAL _ Vector3d.Mag[tan];
mag2: REAL _ mag*mag;
RETURN Vector3d.Div[Vector3d.Cross[Vector3d.Cross[tan, acc], tan], mag2*mag2];
};
};

RefVec: PUBLIC PROC [c: Coeffs, t: REAL] RETURNS [Triple] ~ {
ref: Triple _ Curvature[c, t];
IF Vector3d.Square[ref] < 0.9 THEN ref _ Vector3d.Ortho[Tangent[c, t]];
RETURN[ref];
};
Length: PUBLIC PROC [c: Coeffs] RETURNS [REAL] ~ {
sum: REAL _ 0.0;
fd: Coeffs _ FwdDif[c, 100];
FOR i: NAT IN[0..100) DO
sum _ sum+Vector3d.Mag[[fd[1][0], fd[1][1], fd[1][2]]];
FOR j: NAT IN[0..2] DO
fd[1][j] _ fd[1][j]+fd[2][j];
fd[2][j] _ fd[2][j]+fd[3][j];
ENDLOOP;
ENDLOOP;
RETURN[sum];
};
ConvexHullArea: PUBLIC PROC [b: Bezier] RETURNS [REAL] ~ {
TwiceTriArea: PROC [p0, p1, p2: Triple] RETURNS [REAL] ~ {
RETURN[Vector3d.Mag[Vector3d.Cross[Vector3d.Sub[p1, p0], Vector3d.Sub[p2, p1]]]];
};
a1: REAL _ TwiceTriArea[b.b0, b.b1, b.b2];
a2: REAL _ TwiceTriArea[b.b2, b.b3, b.b0];
a3: REAL _ TwiceTriArea[b.b0, b.b1, b.b3];
a4: REAL _ TwiceTriArea[b.b1, b.b2, b.b3];
RETURN[0.25*(a1+a2+a3+a4)];
};

ConvexHullLength: PUBLIC PROC [b: Bezier] RETURNS [REAL] ~ {
d0: Triple _ Vector3d.Sub[b.b1, b.b0];
d1: Triple _ Vector3d.Sub[b.b2, b.b1];
d2: Triple _ Vector3d.Sub[b.b3, b.b2];
d3: Triple _ Vector3d.Sub[b.b3, b.b0];
RETURN[Vector3d.Mag[d0]+Vector3d.Mag[d1]+Vector3d.Mag[d2]+Vector3d.Mag[d3]];
};

FlatBezier: PUBLIC PROC [b: Bezier, epsilon: REAL _ 0.05] RETURNS [BOOL] ~ {
d3mag: REAL;
d0: Triple _ Vector3d.Sub[b.b1, b.b0];
d1: Triple _ Vector3d.Sub[b.b2, b.b1];
d2: Triple _ Vector3d.Sub[b.b3, b.b2];
d3: Triple _ Vector3d.Sub[b.b3, b.b0];
IF Vector3d.Dot[d3, d0] < 0.0 OR Vector3d.Dot[d3, d2] < 0.0 THEN RETURN[FALSE]; -- bulge
IF (d3mag _ Vector3d.Mag[d3]) < 0.000001 THEN RETURN[TRUE];
RETURN[
(Vector3d.Mag[d0]+Vector3d.Mag[d1]+Vector3d.Mag[d2])/Vector3d.Mag[d3]<1.0+epsilon];
};
Tiny: PUBLIC PROC [c: Coeffs, epsilon: REAL _ 0.05] RETURNS [BOOL] ~ {
RETURN[Vector3d.Distance[Position[c, 0.0], Position[c, 1.0]] < epsilon];
};
Resolution: PUBLIC PROC [c: Coeffs, epsilon: REAL] RETURNS [INTEGER] ~ {
MaxAccel: PROC [curve: Coeffs] RETURNS [REAL] ~ {
Bigger: PROC [r0, r1: REAL] RETURNS [REAL] ~ {RETURN[IF r0 > r1 THEN r0 ELSE r1]};
max: REAL _ 0.0;
a0: Triple _ Acceleration[curve, 0.0];
a1: Triple _ Acceleration[curve, 1.0];
max _ max+Bigger[ABS[a0.x], ABS[a1.x]];
max _ max+Bigger[ABS[a0.y], ABS[a1.y]];
RETURN[max];
};
RETURN[MAX[1, Real.RoundI[Real.SqRt[MaxAccel[c]/(8.0*epsilon)]]]];
};
NearestPoint: PUBLIC PROC [p: Triple, c: Coeffs, t0: REAL _ 0.0, t1: REAL _ 1.0, epsilon: REAL _ 0.01] RETURNS [cPt: Triple, t, dist: REAL] ~ {
InnerNear: PROC [t0, t1: REAL, p0, p1: Triple, d0, d1: REAL]
RETURNS [cPtI: Triple, tI, distI: REAL]~ {
tt: REAL _ 0.5*(t0+t1);
cp: Triple _ Position[c, tt];
d: REAL _ Vector3d.Distance[p, cp];
dc0: REAL _ Vector3d.Distance[p0, cp];
dc1: REAL _ Vector3d.Distance[p1, cp];
IF d < min THEN min _ d;
IF Vector3d.Distance[p0, p1] < MAX[0.0001, epsilon*d] THEN RETURN[cp, tt, d];
IF dc0 < d0 AND dc1 < d1
THEN SELECT TRUE FROM
MIN[d0, d, d1] > min => {cPtI _ cp; tI _ tt; distI _ d;};
d0 < d1 => [cPtI, tI, distI] _ InnerNear[t0, tt, p0, cp, d0, d];
ENDCASE => [cPtI, tI, distI] _ InnerNear[tt, t1, cp, p1, d, d1]
ELSE {
tt0, tt1, dist0, dist1: REAL;
pt0, pt1: Triple;
[pt0, tt0, dist0] _ InnerNear[t0, tt, p0, cp, d0, d];
[pt1, tt1, dist1] _ InnerNear[tt, t1, cp, p1, d, d1];
IF dist0 < dist1
THEN RETURN[pt0, tt0, dist0]
ELSE RETURN[pt1, tt1, dist1];
};
};
min: REAL _ 10000.0;
p0: Triple _ Position[c, t0];
p1: Triple _ Position[c, t1];
[cPt, t, dist] _ InnerNear[t0, t1, p0, p1, Vector3d.Distance[p, p0], Vector3d.Distance[p, p1]];
};

NearestLine: PUBLIC PROC [line: Line, c: Coeffs, t0: REAL _ 0.0, t1: REAL _ 1.0, epsilon: REAL _ 0.01] RETURNS [cPt, lPt: Triple, t, dist: REAL] ~ {
InnerNear: PROC [t0, t1: REAL, p0, p1: Triple, d0, d1: REAL]
RETURNS [cPtI, lPtI: Triple, tI, distI: REAL] ~ {
tt: REAL _ 0.5*(t0+t1);
cp: Triple _ Position[c, tt];
lp: Triple _ Vector3d.LinePoint[cp, line];
d: REAL _ Vector3d.Distance[lp, cp];
dc0: REAL _ Vector3d.Distance[p0, cp];
dc1: REAL _ Vector3d.Distance[p1, cp];
IF d < min THEN min _ d;
IF Vector3d.Distance[p0, p1] < MAX[0.0001, epsilon*d] THEN RETURN[cp, lp, tt, d];
IF dc0 < d0 AND dc1 < d1
THEN SELECT TRUE FROM
MIN[d0, d, d1] > min => {cPtI _ cp; lPtI _ lp; tI _ tt; distI _ d;};
d0 < d1 => [cPtI, lPtI, tI, distI] _ InnerNear[t0, tt, p0, cp, d0, d];
ENDCASE => [cPtI, lPtI, tI, distI] _ InnerNear[tt, t1, cp, p1, d, d1]
ELSE {
tt0, tt1, dist0, dist1: REAL;
cPt0, lPt0, cPt1, lPt1: Triple;
[cPt0, lPt0, tt0, dist0] _ InnerNear[t0, tt, p0, cp, d0, d];
[cPt1, lPt1, tt1, dist1] _ InnerNear[tt, t1, cp, p1, d, d1];
IF dist0 < dist1
THEN RETURN[cPt0, lPt0, tt0, dist0]
ELSE RETURN[cPt1, lPt1, tt1, dist1];
};
};
min: REAL _ 10000.0;
p0: Triple _ Position[c, t0];
p1: Triple _ Position[c, t1];
[cPt, lPt, t, dist] _
InnerNear[t0, t1, p0, p1, DistanceToLine[p0, line], DistanceToLine[p1, line]];
};

DistanceToLine: PROC [p: Triple, l: Line] RETURNS [REAL] ~ {
lp: Triple _ Vector3d.LinePoint[p, l];
RETURN[Vector3d.Distance[p, lp]];
};

NearestSpline: PUBLIC PROC [c1, c2: Coeffs, epsilon: REAL _ 0.01] RETURNS [t1, t2: REAL] ~ {
InnerNear: PROC [t0, t1: REAL, p0, p1: Triple, d0, d1: REAL] RETURNS [tI, dI: REAL] ~ {
tt: REAL _ 0.5*(t0+t1);
pp: Triple _ Position[c1, tt];
dc0: REAL _ Vector3d.Distance[p0, pp];
dc1: REAL _ Vector3d.Distance[p1, pp];
d: REAL _ NearestPoint[pp, c2, 0.0, 1.0].dist;
IF d < min THEN min _ d;
IF Vector3d.Distance[p0, p1] < MAX[0.0001, epsilon*d] THEN RETURN[tt, d];
IF dc0 < d0 AND dc1 < d1
THEN {
t, d: REAL;
SELECT TRUE FROM
MIN[d0, d, d1] > min => RETURN[tt, d];
d0 < d1 => {[t, d] _ InnerNear[t0, tt, p0, pp, d0, d]; RETURN[t, d]};
ENDCASE => {[t, d] _ InnerNear[tt, t1, pp, p1, d, d1]; RETURN[t, d]};
}
ELSE {
tt0, tt1, dist0, dist1: REAL;
[tt0, dist0] _ InnerNear[t0, tt, p0, pp, d0, d];
[tt1, dist1] _ InnerNear[tt, t1, pp, p1, d, d1];
IF dist0 < dist1
THEN RETURN[tt0, dist0]
ELSE RETURN[tt1, dist1];
};
};
min: REAL _ 10000.0;
p0: Triple _ Position[c1, 0.0];
p1: Triple _ Position[c1, 1.0];
d0: REAL _ NearestPoint[p0, c2, 0.0, 1.0, epsilon].dist;
d1: REAL _ NearestPoint[p1, c2, 0.0, 1.0, epsilon].dist;
t1 _ InnerNear[0.0, 1.0, p0, p1, d0, d1].tI;
t2 _ NearestPoint[Position[c1, t1], c2, 0.0, 1.0, epsilon].t;
};

NearestPair: PUBLIC PROC [p: Pair, c: Coeffs, persp: BOOL _ FALSE, t0: REAL _ 0.0, t1: REAL _ 1.0] RETURNS [pRet: Pair, t, dist: REAL] ~ {
InnerNear: PROC [t0, t1: REAL, p0, p1: Pair, d0, d1: REAL] ~ {
tt: REAL _ 0.5*(t0+t1);
pp: Pair _ PairPosition[c, tt, persp];
dd: REAL _ PairToPairDistance[[p.x, p.y], pp];
IF PairToPairDistance[p0, p1] < 0.01*dd THEN {
pRet _ pp;
t _ tt;
dist _ dd;
RETURN;
};
IF d0 < d1
THEN InnerNear[t0, tt, p0, pp, d0, dd]
ELSE InnerNear[tt, t1, pp, p1, dd, d1];
};
p0: Pair _ PairPosition[c, t0, persp];
p1: Pair _ PairPosition[c, t1, persp];
InnerNear[t0, t1, p0, p1, PairToPairDistance[[p.x, p.y], p0], PairToPairDistance[[p.x, p.y], p1]];
};

PairToPairDistance: PROC [p0, p1: Pair] RETURNS [REAL] ~ {
a: REAL _ p0.x-p1.x;
b: REAL _ p0.y-p1.y;
RETURN[Real.SqRt[a*a+b*b]];
};

PairPosition: PROC [c: Coeffs, t: REAL, persp: BOOL _ FALSE] RETURNS [Pair] ~ {
ret: ARRAY [0..3] OF REAL;
t2, t3: REAL;
nCoords: NAT _ IF persp THEN 4 ELSE 2;
IF t = 0.0 THEN
FOR i: NAT IN[0..nCoords) DO ret[i] _ c[3][i]; ENDLOOP
ELSE IF t = 1.0 THEN
FOR i: NAT IN[0..nCoords) DO ret[i] _ c[0][i]+c[1][i]+c[2][i]+c[3][i]; ENDLOOP
ELSE {
t2 _ t*t;
t3 _ t*t2;
FOR i: NAT IN[0..nCoords) DO ret[i] _ t3*c[0][i]+t2*c[1][i]+t*c[2][i]+c[3][i]; ENDLOOP;
};
IF persp THEN {
w: REAL _ ret[3];
IF ret[0]+w < 0.0 OR ret[0]+w < 0.0 OR ret[1]+w < 0.0 OR ret[1]-w < 0.0 OR ret[2]+w < 0.0
THEN RETURN[[2000.0, 2000.0]]
ELSE RETURN[[ret[0]/w, ret[1]/w]];
};
RETURN[[ret[0], ret[1]]];
};

FurthestPoint: PUBLIC PROC [c: Coeffs] RETURNS [p: Triple, t, dist: REAL] ~ {
max: REAL _ 0.0;
base: Triple _ Position[c, 0.0];
axis: Triple _ Vector3d.Sub[Position[c, 1.0], base];
FOR tt: REAL _ 0.0, tt+0.01 WHILE tt <= 1.0 DO
pc: Triple _ Position[c, tt];
pl: Triple _ Vector3d.LinePoint[pc, [base, axis]];
d: REAL _ Vector3d.Distance[pc, pl];
IF d > max THEN {max _ d; p _ pc; t _ tt; dist _ d};
ENDLOOP;
};
csplit0: Matrix _ NEW[MatrixRep _ [
			[1.0/8.0,	0.0,		0.0,		0.0],
			[0.0,		1.0/4.0,	0.0,		0.0],
			[0.0,		0.0,		1.0/2.0,	0.0],
			[0.0,		0.0,		0.0,		1.0]]];

csplit1: Matrix _ NEW[MatrixRep _ [
			[1.0/8.0,  0.0,			0.0,     	0.0],
			[3.0/8.0,  1.0/4.0,		0.0,     	0.0],
			[3.0/8.0,  1.0/2.0,		1.0/2.0,	0.0],
			[1.0/8.0,  1.0/4.0,		1.0/2.0,	1.0]]];

SplitCurve: PUBLIC PROC [c: Coeffs, out1, out2: Coeffs _ NIL] RETURNS [c1, c2: Coeffs] ~ {
aa, bb, cc, dd, ee, ff: REAL;
c1 _ IF out1 # NIL THEN out1 ELSE NEW[CoeffsRep];
c2 _ IF out2 # NIL THEN out2 ELSE NEW[CoeffsRep];
FOR i: NAT IN[0..2] DO
c1[0][i] _	aa _ (1.0/8.0)*c[0][i];
c1[1][i] _	bb _ (1.0/4.0)*c[1][i];
c1[2][i] _	cc _ (1.0/2.0)*c[2][i];
c1[3][i] _	dd _ c[3][i];
			ee _ 3.0*aa;
c2[0][i] _	aa;
c2[1][i] _	ff _ ee+bb;
c2[2][i] _	ff+bb+cc;
c2[3][i] _	aa+bb+cc+dd;
ENDLOOP;
};
SplitBezier: PUBLIC PROC [b: Bezier] RETURNS [b1, b2: Bezier] ~ {
Ave: PROC [p, q: Triple] RETURNS [Triple] ~ INLINE {
RETURN[[0.5*(p.x+q.x), 0.5*(p.y+q.y), 0.5*(p.z+q.z)]];
};
b01: Triple _ Ave[b.b0, b.b1];
b12: Triple _ Ave[b.b1, b.b2];
b23: Triple _ Ave[b.b2, b.b3];
b0112: Triple _ Ave[b01, b12];
b1223: Triple _ Ave[b12, b23];
b01121223: Triple _ Ave[b0112, b1223];
RETURN[[b.b0, b01, b0112, b01121223], [b01121223, b1223, b23, b.b3]];
};
Subdivide: PUBLIC PROC [c: Coeffs, t: REAL] RETURNS [c1, c2: Coeffs] ~ {
RETURN[Reparameterize[c, 0.0, t], Reparameterize[c, t, 1.0]];
};
Reparameterize: PUBLIC PROC [in: Coeffs, t0, t1: REAL, out: Coeffs_NIL] RETURNS [Coeffs] ~ {
dt1: REAL _ t1-t0;
dt2: REAL _ dt1*dt1;
t02: REAL _ t0*t0;
IF out = NIL THEN out _ NEW[CoeffsRep];
FOR i: NAT IN[0..2] DO
x: REAL _ in[0][i];
y: REAL _ in[1][i];
z: REAL _ in[2][i];
x3: REAL _ 3.0*x;
out[0][i] _ x*dt1*dt2;
out[1][i] _ x3*dt2*t0+y*dt2;
out[2][i] _ x3*dt1*t02+2.0*y*dt1*t0+z*dt1;
out[3][i] _ x*t0*t02+y*t02+z*t0+in[3][i];
ENDLOOP;
FOR i: NAT IN[0..3] DO out[i][3] _ in[i][3]; ENDLOOP;
RETURN[out];
};
Tame: PUBLIC PROC [in: Coeffs, out: Coeffs _ NIL] RETURNS [Coeffs] ~ {
v0, v2, m: REAL;
d0, d2, d3, s0, s2, r0, r2: Triple;
b: Bezier _ BezierFromCoeffs[in];									-- convert to bezier
d0 _ Vector3d.Sub[b.b1, b.b0];
d2 _ Vector3d.Sub[b.b3, b.b2];
d3 _ Vector3d.Sub[b.b3, b.b0];
m _  Vector3d.Mag[d3];
v0 _ m*Vector3d.Dot[d0, d3];
v2 _ m*Vector3d.Dot[d2, d3];
s0 _ Vector3d.Mul[d3, v0];										-- projection of d0 onto d3
s2 _ Vector3d.Mul[d3, v2];										-- projection of d2 onto d3
r0 _ Vector3d.Sub[d0, s0];
r2 _ Vector3d.Sub[d2, s2];
IF Vector3d.Dot[d3, d0] < 0.0 THEN {								-- bulge at start
d0 _ r0;
s0 _ [0.0, 0.0, 0.0];
};
IF Vector3d.Dot[d3, d2] < 0.0 THEN {								-- bulge at end
d2 _ r2;
s2 _ [0.0, 0.0, 0.0];
};
IF Vector3d.Dot[r0, r2] > 0.0 THEN {								-- sides point oppositely
IF Vector3d.Square[r0] > Vector3d.Square[r2]
THEN d2 _ s2
ELSE d0 _ s0;
};
IF v0+v2 > m THEN {												-- sides too long
scale: REAL _ m/(v0+v2);
d0 _ Vector3d.Mul[d0, scale];
d2 _ Vector3d.Mul[d2, scale];
};
b.b1 _ Vector3d.Add[b.b0, d0];
b.b2 _ Vector3d.Sub[b.b3, d2];
RETURN [CoeffsFromBezier[b, out]];
};
Same: PUBLIC PROC [c1, c2: Coeffs] RETURNS [BOOL] ~ {
FOR i: NAT IN[0..3] DO
FOR j: NAT IN[0..3] DO IF c1[i][j] # c2[i][j] THEN RETURN[FALSE]; ENDLOOP;
ENDLOOP;
RETURN[TRUE];
};

debugView: Coeffs _ NIL;
DebugView: PUBLIC PROC [view: Coeffs] ~ {debugView _ view};

debugStream: IO.STREAM _ NIL;
DebugStream: PUBLIC PROC [stream: IO.STREAM] ~ {debugStream _ stream};

END.
..
References:
/Cedar/CedarChest6.0/CubicSplinePackage/*.mesa
Rogers and Adams: Mathematical Elements for Computer Graphics
Notes:
ςSpline3dImpl.mesa
Copyright c 1985 by Xerox Corporation.  All rights reserved.
Bloomenthal, May 23, 1986 1:15:03 am PDT
Basic Types
Interpolation Procedures

From /Cedar/CedarChest6.0/CubicSplinePackage/RegularALSplineImpl.mesa, in which
coeffs.t3, .t2, .t1, and .t0 referred to the 3rd, 2nd and 1st derivatives, and position of
the curve.  In Spline3d, these correspond to coeffs[0], [1], [2], and [3].

set cubic constants:
compute parametrization (chord length approximation):
cyclic end conditions: precompute coefficients a and c:
compute each dimension separately:
compute constant coefficient:
cyclic end conditions: compute coefficients b, r, & s:
compute coefficient t2 (second derivative/2):

note: coeffs[lastK+1].t2.[i] = coeffs[0].t2[i]
compute coefficients t3 (third derivative/6) & and t1 (first derivative):

note again: coeffs[lastK+1].t2[i] = coeffs[0].t2[i] (0 for natural end conditions)
reparametrize coefficients to interval [0..1]:
Mostly from Rogers and Adams, with tension added:

coeffs[m] _ CoeffsFromHermite[
[knots[m],
 knots[m+1],
 Vector3d.Mul[tangents[m], chords[m]],
 Vector3d.Mul[tangents[m+1], chords[m]]],
coeffs[m]];

Conversion Procedures




Evaluation Procedures



1 sqrt, 3 divides, 14 multiplies, 6 adds
Size Procedures



This produces inconsistent results:
Search Procedures
This doesn't always work right.
Return point on c1 closest to c2 within t0,t1
Subdivision Procedures



Miscellaneous procedures

Curvature: PUBLIC PROC [c: Coeffs, t: REAL] RETURNS [Triple] ~ {
1 sqrt, 3 divides, 8 multiplies, 7 adds but this gives wrong length
tan: Triple _ Tangent[c, t];
IF tan = [0.0, 0.0, 0.0] THEN RETURN[tan] ELSE {
acc: Triple _ Acceleration[c, t];
u: REAL _ tan.x*tan.x+tan.y*tan.y+tan.z*tan.z;
RETURN Vector3d.Div[
Vector3d.Sub[
Vector3d.Mul[acc, u],
Vector3d.Mul[tan, Vector3d.Dot[acc, tan]]],
Real.SqRt[u*u*u]];
};
};

TPosFromCoeffs: PUBLIC PROC [c: Coeffs, num: NAT, t0: REAL _ 0.0, t1: REAL _ 1.0, tpos: TPosSequence _ NIL] RETURNS [TPosSequence] ~ {
Tri: TYPE ~ RECORD [t0, t1, t, area: REAL, p0, p1, p: Triple];
TriSequenceRec: TYPE ~ RECORD [element: SEQUENCE length: NAT OF Tri];
n, nTris: NAT _ 0;
tris: REF TriSequenceRec _ NEW[TriSequenceRec[num]];
LargestTri: PROC RETURNS [n: NAT] ~ {
max: REAL _ 0.0;
FOR i: NAT IN[0..nTris) DO IF tris[i].area > max THEN {max _ tris[i].area; n _ i} ENDLOOP;
};
AddTri: PROC [tri: Tri] ~ {tris[nTris] _ tri; nTris _ nTris+1};
RelArea: PROC [p0, p1, p2: Triple] RETURNS [REAL] ~ {
RETURN[Vector3d.Square[Vector3d.Cross[Vector3d.Sub[p1, p0], Vector3d.Sub[p2, p0]]]];
};
NewTri: PROC [t0, t1: REAL, p0, p1: Triple] RETURNS [Tri] ~ {
t: REAL _ 0.5*(t0+t1);
p: Triple _ Position[c, t];
RETURN[[t0, t1, t, RelArea[p0, p, p1], p0, p1, p]];
};
p0: Triple _ Position[c, t0];
p1: Triple _ Position[c, t1];
IF tpos = NIL THEN tpos _ NEW[TPosSequenceRec[num]];
AddTri[NewTri[0.0, 1.0, p0, p1]];
THROUGH [0..num-4] DO
tri: Tri _ tris[n _ LargestTri[]];
tris[n] _ NewTri[tri.t0, tri.t, tri.p0, tri.p];
AddTri[NewTri[tri.t, tri.t1, tri.p, tri.p1]];
ENDLOOP;
tpos[0] _ [t0, p0];
tpos[num-1] _ [t1, p1];
FOR n: NAT IN[1..nTris] DO
isave: NAT;
min: REAL _ 1.0;
FOR i: NAT IN[0..nTris) DO IF tris[i].t < min THEN {min _ tris[i].t; isave _ i} ENDLOOP;
tpos[n] _ [min, tris[isave].p];
tris[isave].t _ 1.0;
ENDLOOP;
RETURN[tpos];
};

Κ!˜
codeš
œ™Jšœ
Οm
œ1™<J™(—J˜JšΟk	œ
žœ+˜7J˜šΠlnœž
˜Kšžœ˜ Kšžœ	˜—Kšœž˜
headš
œ™Jšœžœ˜"Jšœžœ˜'J˜Jšœ
žœ˜Jšœ
žœ˜Jšœ
žœ˜Iproc˜Jšœžœ˜"Jšœžœ˜0Jšœžœ˜5J˜Jšœžœ˜,Jšœžœ˜1J˜Jšœžœ˜,Jšœžœ˜1J˜Mšœžœ˜"Mšœžœ˜$Mšœžœ˜%J˜Jšœžœ˜#Jšœ
žœ˜'Jšœžœ˜0Jšœžœ˜5J˜Jšœžœ˜,Jšœžœ˜1—š
œ™šΟnœžœ
žœ žœ˜IJšžœ˜J™™OJ™ZJ™JJ™—J˜JšœžœžœΟc˜FJšœžœ‘˜9J˜"J˜Kšžœžœ
žœ
žœ˜%K˜Kš
žœžœ
žœžœ
žœ˜XK˜šžœžœ˜J˜?Jšžœžœ
žœžœ
žœžœ
žœžœžœžœ
˜Všžœžœ˜Jš
œ(˜(Jš
œ(˜(Jš
œ(˜(J˜—Jšžœ	˜Jš
œ˜—J˜Jšœžœ˜"Jšœžœ˜"Jšœžœ˜"Jšœžœ˜"Jšœžœ˜"Jšœžœ˜"Jšœžœ˜"J˜Jš
œ™šžœžœ
žœž˜Jš
œ?˜?Jšžœ
˜—J˜Jš
œ5™5J˜J˜Jš
œ7™7Jš
œ˜Jš
žœžœ
žœžœ/žœ
˜SJš
œ˜Jš
žœžœ
žœžœžœ
˜BJ˜Jš
œ"™"šžœžœ
žœž˜Jšœžœ
˜	Jš
œ™šžœžœ
žœž˜šœžœž˜Jš
œ&˜&Jš
œ&˜&Jšžœ#˜*—Jšžœ
˜——˜Jš
œ6™6J˜Jš
žœžœ
žœžœžœ
˜>Jš
žœžœ
žœžœ#žœ
˜FJ˜J˜J˜š	žœžœ
ž
œ
žœž˜&J˜ J˜Jšžœ
˜——˜Jš
œ-™-J˜CJ˜HJ™Jš
œ.™.J˜Jš
žœžœ
žœžœ&žœ
˜J—˜Jš
œI™Išžœžœ
žœž˜J˜:J˜HJšžœ
˜—J™Jš
œR™RJ˜@J˜V—˜Jš
œ.™.šžœžœ
žœž˜J˜,J˜,J˜'Jšžœ
˜——˜Jšžœ
˜—J˜Jšžœ	˜˜J˜——š
 œžœFžœ!žœžœ˜€
J™1šžœž˜Jšœžœ	˜š
œ	˜	Jšœžœ)˜2Kš
žœ
žœ
žœžœ
žœ˜Pš
œ˜Jš
œ
˜
Jš
œ
˜
Jš
œ˜Jš
œ˜Jš
œ˜—Jšžœ	˜J˜—šžœ˜Jšœ,žœ˜1Jšœ:žœ˜KJšžœ:˜@J˜——J˜J˜—š 
œžœ.žœ
˜EJšžœ˜šžœ	žœ
žœ ˜2Kšžœ
žœ ˜1—šžœžœ
žœž˜#Jšžœ=žœ˜TJšžœ
˜—Jšžœ	˜Jš
œ˜J˜—š œžœHžœžœ˜’
Jšœžœ$‘˜SJšœžœ$‘˜Nšžœžœ
žœ"˜6Kšžœžœ"˜5—J˜<J˜1Jšžœ˜J˜J˜—š œžœd˜zJšœžœ˜šžœžœ˜Jš
œ˜J˜EJ˜
—šžœ˜Jš
œ˜J˜J˜—šžœžœ˜Jš
œ˜J˜IJ˜
—šžœ˜Jš
œ˜J˜J˜—Jš
œ˜J˜—š œžœP˜gJšœžœ˜J˜šžœžœ
žœž˜Jšœžœ˜Jšœžœ
˜J˜š
œ˜š
œ
˜
Jš
œ8˜8Jš
œ9˜9—J˜
—Jšžœ
˜—J˜š	žœžœ
žœ
žœ‘˜;Jšœžœ
˜Jšžœžœ
žœ
˜Jš
œ˜J˜3J˜3Jš
œ˜J˜J˜Jšžœ
˜—J˜Jšœ:‘˜Nšžœžœ
žœ
ž˜š
œ˜š
œ
˜
š
œ
˜
Jš
œ
˜
Jš
œ/˜/—Jš
œ
˜
——Jšžœ
˜—Jš
œ˜J˜—š
 
œžœ:žœ1žœžœ˜£
Kšœ	žœ˜K˜šžœ
žœ
žœ˜-Kšžœ
žœ˜.—K˜šžœžœ
žœžœ
žœ
ž˜2Kš
œ1˜1Kšžœ
˜K™—šžœžœ
žœ
ž˜š
œ™Jš
œ
™
Jš
œ™Jš
œ&™&Jš
œ)™)Jš
œ™J™—Jšœžœ
˜Jš
œ4˜4Jš
œ)˜)Jš
œ<˜<Jš
œ8˜8Jš
œH˜HJšžœ
žœ
žœ
žœ˜3Jš
œ$˜$Jš
œ$˜$Jš
œ$˜$Jš
œ9˜9Jšžœ
˜—J˜Jšžœ	˜Jš
œ˜——š
œ™šœžœ˜+Jš
œ˜Jš
œ˜Jš
œ˜Jš
œ˜J˜—šœžœ˜+Jš
œ˜Jš
œ˜Jš
œ˜Jš
œ˜J˜—šœžœ˜*Jš
œ˜Jš
œ˜Jš
œ˜Jš
œ˜J˜—šœžœ˜*Jš
œ˜Jš
œ ˜ Jš
œ#˜#Jš
œ˜J˜—šœžœ˜+Jš
œ,˜,Jš
œ,˜,Jš
œ,˜,Jš
œ.˜.J˜—šœžœ˜+Jš
œ ˜ Jš
œ˜Jš
œ ˜ Jš
œ"˜"—J˜š œžœ9žœžœ
˜[Kšžœ
˜Kš	œžœžœ	žœžœ˜'Kš
œ$˜$˜J˜J˜J˜J˜—Kšžœžœ
žœžœ˜'Kš
œ#˜#K˜Kšžœ˜K˜K˜—š œžœžœ˜YKš
œ$˜$Kš
œ!˜!K˜!K˜!K˜!K˜!K˜Kšžœ˜K˜K˜—š œžœžœžœ
˜QKšžœ?˜EK˜K™—š œžœžœžœ
˜SKšžœ@˜FK˜K™—š
 œžœ
žœžœžœ
˜SKšžœCžœ˜PK˜K˜—š œžœ
žœ
˜>Kš
œ˜Kš
œ8˜8Kšžœ˜K˜K™—š œžœ
žœ˜@Kš
œ˜Kš
œ9˜9Kšžœ˜K˜K™—š œžœ
žœ˜@Kš
œ˜Kš
œ=˜=Kšžœ˜K˜——š
œ™Kš œžœ˜+K˜š 
œžœ
žœ*žœžœžœ˜fš œžœ˜šžœ˜Kšžœ˜šžœ˜K˜K˜K˜K˜K˜——K˜—K˜	Kšžœžœ˜K˜J™—š œžœ
žœžœžœžœ
˜ZKšœžœ˜2Kšœžœ
˜Jšžœžœ
žœ
žœ˜"K˜K˜K˜K˜˜Jš
œ˜Jš
œ˜Jš
œ˜Jš
œ˜—Jšžœ ˜&K˜K˜—šΠbnΟbžœ
ž£
œžœžœ
ž£
œ˜mJšœžœžœ
žœ
˜Jš
œ&˜&šžœ
žœ
žœž˜2Jšœ	žœ˜)—J˜7Jš
žœžœ
žœžœžœ
˜8šžœžœ
žœ
ž˜ šžœžœ
žœž˜J˜J˜)J˜)Jšžœ
˜—J˜Jšžœ
˜—Jšžœ	˜J˜J˜—š œžœžœžœ
˜?Jšœžœžœ
žœ
˜Jš
žœ	žœ
žœžœ
žœžœž˜@Jšžœ
žœ	žœ
žœžœ
žœžœ*ž˜]šžœ˜Jšœžœ˜Jšœžœ˜Jš	žœžœ
žœžœ3ž˜PJš
œ˜—Jšžœ˜!J˜J˜—š œžœžœžœ
˜?Kšœžœžœ
žœ
˜Kšœžœ˜Kšœžœ˜Kšžœ	žœ
žœ˜4Kšžœ	žœ˜%Kš
žœžœ
žœžœ'žœ
˜EKšžœ˜!K˜K™—š œžœžœžœ
˜CKšœžœ	˜KšžœT˜ZK˜K™—š
 œžœ
žœžœžœ˜@Kšžœ˜Kš
œ˜—K˜š œžœžœžœ˜GJ˜"J˜6Jšœžœžœžœ/˜WJš
žœ	žœ	žœ
žœ	žœ	˜5J˜J˜—š œžœ
žœžœžœžœ˜AJš
œ˜Jš
œ!˜!Jšœžœ˜Jšžœ
žœ	žœžœ1˜PJ˜J˜—š 	œžœžœžœ
˜@J™(Jš
œ˜šžœžœ
žœžœ˜0Jš
œ!˜!Jšœžœ˜Jšœžœ˜JšžœH˜NJ˜—J˜J˜—š œžœžœžœ
˜=Jš
œ˜Jšžœžœ%˜GJšžœ˜Jš
œ˜——š
œ™š œžœ
žœžœ˜2Kšœžœ˜Kš
œ˜šžœžœ
žœ	ž˜Kš
œ7˜7šžœžœ
žœž˜K˜K˜Kšžœ
˜—Kšžœ
˜—Kšžœ˜K˜K™—š œžœ
žœžœ˜:š œžœžœžœ˜:JšžœK˜QJ˜—Jšœžœ"˜*Jšœžœ"˜*Jšœžœ"˜*Jšœžœ"˜*Jšžœ˜K˜K˜—š
 œžœ
žœ
žœžœ˜<K˜&K˜&K˜&K˜&KšžœF˜LK˜K˜—š
 
œžœžœ	žœžœ˜LKšœžœ
˜Kš
œ&˜&K˜&K˜&K˜&Kšžœžœžœ
žœ
žœ‘˜XKšžœ'žœ
žœ
žœ˜;šžœ
˜Kš
œS˜S—K˜K™—š œžœ
žœžœ	žœžœ˜FKšžœB˜HK˜K™—š
 
œžœžœžœžœ˜HJ™#š œžœžœ
žœ˜1Jš œžœ
žœžœžœžœ
žœ	žœžœ˜RJšœžœ˜Jš
œ&˜&Jš
œ&˜&Jšœžœžœ˜'Jšœžœžœ˜'Jšžœ˜Jš
œ˜—Jšžœ
žœ8˜BJ˜——š
œ™š £œžœ
žœžœžœžœžœžœ˜
š 	œžœ
žœžœ
˜<Jšžœžœ˜*Jšœžœ˜J˜Jšœžœ˜#Jšœžœ˜&Jšœžœ˜&Jšžœ	žœ	˜Jšžœžœžœ
žœ˜Mšžœ
žœ	˜šžœ
žœ
žœ
ž˜Jšžœ6˜9Jš
œ@˜@Jšžœ8˜?—šžœ˜Jšœžœ
˜Jš
œ˜Jš
œ5˜5Jš
œ5˜5šžœ˜Jšžœ
žœ˜Jšžœ
žœ˜—J˜——J˜—Jšœžœ˜J˜J˜Jš
œ_˜_˜K˜——š œžœ
žœžœžœžœ	žœžœ˜”
š 	œžœ
žœžœ
˜<Jšžœžœ˜1Jšœžœ˜J˜Jš
œ*˜*Jšœžœ˜$Jšœžœ˜&Jšœžœ˜&Jšžœ	žœ	˜Jšžœžœžœ
žœ˜Qšžœ
žœ	˜šžœ
žœ
žœ
ž˜JšžœA˜DJš
œF˜FJšžœ>˜E—šžœ˜Jšœžœ
˜Jš
œ˜J˜<J˜<šžœ˜Jšžœ
žœ˜#Jšžœ
žœ˜$—J˜——J˜—Jšœžœ˜Jš
œ˜J˜š
œ˜Jš
œN˜N—˜K˜——š œžœžœžœ˜<J˜&Jšžœ˜!J˜J˜—š
 £œžœ
žœžœžœ
žœ˜\J™š 	œžœ
žœžœžœ
žœ˜WJ™-Jšœžœ˜J˜Jšœžœ˜&Jšœžœ˜&Jšœžœ'˜.Jšžœ	žœ	˜Jšžœžœžœ
žœ˜Išžœ
žœ	˜šžœ˜Jšœžœ
˜šžœ
žœ
ž˜Jšžœžœ˜&Jšœ7žœ˜EJšžœ0žœ˜E—J˜
—šžœ˜Jšœžœ
˜Jš
œ0˜0Jš
œ0˜0šžœ˜Jšžœ
žœ˜Jšžœ
žœ
˜—J˜——J˜—Jšœžœ˜J˜J˜Jšœžœ0˜8Jšœžœ0˜8Jš
œ,˜,Jš
œ=˜=˜K˜——š ’œžœ
žœžœžœžœžœžœžœ˜Š
š 	œžœ
žœžœ˜>Jšœžœ˜Jš
œ&˜&Jšœžœ&˜.šžœ&žœ˜.Jš
œ
˜
Jš
œ˜Jš
œ
˜
Jšžœ
˜J˜—šžœ˜
Jšžœ"˜&Jšžœ#˜'—J˜—Jš
œ&˜&Jš
œ&˜&Jš
œb˜b˜K˜——š’œžœžœžœ˜:Jšœžœ
˜Jšœžœ
˜Jšžœ˜J˜J˜—š œžœžœ	žœžœžœ˜OJšœžœžœ
žœ
˜Jšœžœ
˜
Jš	œ	žœžœžœžœ˜&šžœ	ž˜Jš	žœžœ
žœ
žœž˜6—šžœ
žœ	ž˜Jš	žœžœ
žœ
žœ+ž˜N—šžœ˜Jš
œ	˜	Jš
œ
˜
Jš	žœžœ
žœ
žœ3ž˜WJš
œ˜—šžœžœ˜Jšœžœ
˜š
žœžœžœžœžœ˜YJšžœ
žœ˜Jšžœ
žœ˜"—J˜—Jšžœ˜J˜J˜—š 
œžœ
žœžœ˜MKšœžœ˜Kš
œ ˜ K˜4šžœžœžœ
ž˜.K˜Kš
œ2˜2Kšœžœ˜$Kšžœ	žœ%˜4Kšžœ
˜—J˜——š
œ™šœžœ˜#Jš
œ˜Jš
œ˜Jš
œ˜Jš
œ˜J˜—šœžœ˜#Jš
œ$˜$Jš
œ'˜'Jš
œ&˜&Jš
œ(˜(—K˜š
 
œžœ
žœ"žœžœ˜ZJšœžœ
˜Jšœžœžœ
žœžœ
žœ˜1Jšœžœžœ
žœžœ
žœ˜1šžœžœ
žœž˜J˜"J˜"J˜"J˜J˜J˜J˜J˜J˜Jšžœ
˜—J˜J™—š œžœ
žœ˜Aš œžœžœžœ˜4Kšžœ0˜6K˜—Kš
œ˜Kš
œ˜Kš
œ˜Kš
œ˜Kš
œ˜Kš
œ&˜&Kšžœ?˜EK˜K™—š 	œžœžœžœ˜HKšžœ7˜=K˜K™—š œžœ
žœžœžœžœ
˜\Kšœžœ	˜Kšœžœ˜Kšœžœ	˜Kš
žœžœ
žœž
œ
žœ˜'šžœžœ
žœž˜Kšœžœ˜Kšœžœ˜Kšœžœ˜Kšœžœ	˜K˜K˜K˜*K˜)Kšžœ
˜—Kš
žœžœ
žœžœžœ
˜5Kšžœ˜K˜——š
œ™š
 œžœ
žœžœžœ
˜FJšœžœ
˜Jš
œ#˜#Jšœ*‘˜>J˜J˜J˜J˜J˜J˜Jšœ$‘˜?Jšœ$‘˜?Jš
œ˜Jš
œ˜šžœžœ
‘˜=J˜J˜J˜—šžœžœ
‘˜;J˜J˜J˜—šžœžœ
‘˜Ešžœ*˜,Jšžœ˜Jšžœ	˜
—J˜—šžœžœ‘˜0Jšœžœ
˜J˜J˜J˜—J˜J˜Jšžœ˜"J˜J™—š œžœžœžœ˜5šžœžœ
žœž˜Mšžœžœ
žœžœ
žœžœ
žœ
žœž	˜JMšžœ
˜—Mšžœ
žœ˜
M˜—M˜Mšœžœ
˜Mš 	œžœ
žœ%˜;M˜Mšœ
žœ
žœžœ
˜Mš œžœ
žœ
žœ˜FM˜—Kšžœ
˜K˜˜Jš
œ.˜.JšœΟz+˜=—L˜š 	œžœžœžœ
™@J™CJš
œ™šžœžœ
žœžœ™0Jš
œ!™!Jšœžœ'™.šžœ™š
œ
™
Jš
œ™Jš
œ+™+—Jš
œ™—J™—J™J™—š œžœ
žœžœžœžœžœžœ™†
Jšœžœžœžœ™>Jšœžœžœžœ	žœ
žœ™EJšœ
žœ™Jšœžœžœ™4š 
œžœ
žœžœ™%Jšœžœ™Jšžœžœ
žœžœ
žœžœžœ
™ZJ™—Jš œžœ3™?š œžœžœžœ™5JšžœN™TJ™—š œžœ
žœžœ
™=Jšœžœ™J™Jšžœ-™3J™—J™J™Jšžœžœ
žœžœ™4J™!šžœž™Jš
œ"™"Jš
œ/™/Jš
œ-™-Jšžœ
™—J™J™šžœžœ
žœž™Jšœžœ
™Jšœžœ™Jšžœžœ
žœžœžœžœ
™XJ™J™Jšžœ
™—Jšžœ™
J™J™——…—cF‘: