TorusKnotImpl.mesa

 Dobkin August 8, 1986 2:12:03 pm PDT

This is a procedure for creating triangles which cover the surface of the mathematical shape known as the torus knot. The torus knot can be thought of as a mathematical shape created by wrapping a piece of cord about a torus (or doughnut or sphere with one handle). The wrapping is governed by three parameters, A1,A2,A3. A2 represents the number of times that the cord travels around the corcumference of the torus. A1 represnts the number of times the cord makes a complete revolution about the center of the torus and A3 represents the avergae of the inner and outer radii of the torus. We assume here that the inner radius is A3 -1 and the outer radius is A3 - 1.

The triangles are then create by viewing the knot as the curve resulting from a mapping from [0,1] into R3 (Euclidian 3 space) and uniformly taking samples of the curve -- a total of N such samples are taken. Next, at each sample, a circle is approximated by an M-gon to simulate the effect of a tube passing around the curve in 3-space. The Radius variable controls the radius of the M-gonal cross-sections. Each cross-section is taken as the plane perpendicular to the tangent to the curve at that sample point. Axes represnting x and y coordinate axes are then created in these planes. Adjacent cross-sections are lined up so that the the two center points (i.e. samples off the curve) and the two unit-vectors in the x direction (a total of 4 points) are coplanar. In pictures, this makes the spines appear to move straight. At the closing of the knot, it is not possible to do this oriaentation exactly, but an approximation is done which isn't noticable.

The user is also able to generate only fragments of knots, and to transform the knot once it has been generated. The end result can either be appended to a file or be written to a new file. The usage of particular aruments is described in detail below.

I am aware of three remaining things that should be done.
1. The KnotDescriptorRec RECORD should be enhanced to include a transformation matrix. In this case, it would be necessary to change the code in MakeRow to apply the transform to each Point3 geneerated into Retrow. The correct way to do this would be to use the Matrix3 and Vector3 interfaces. In this case, it would be possible to eliminate some of this code which is duplicated there if the Point3 TYPE was changed to the Triple TYPE.

2. The code should allow it's output to be generated to either the current format file or a file that fits into Frank Crow's ThreeDWorld. This is straightforward. In Retrow, when points are generated, they need to be numbered. Then, PutTriangle would be modified to print out triples of integers rather than triples of triples of floating point numbers. This would cause the inner loop of RowofTriangles to Change it's form to something which PutTriangles by vertex indices. It is actually the case that PR[J] is expressible as (I-1)*M + J and TR[J] is expressible as I*M + J for suitable index variable I.

3. FInally, there is a problem in closing the knots. Having taken a sample and the slope of the curve at that point, a plane perpendicular to the slope vector is given x and y axes (Fx and Fy in the code) which are used to generate the M-gonal cross-section. In all planes after the first, the x and y axes are selected to line up with those in the preceeding one. However, this cannot occur when matching together the last and first cross-sections. Thus, when these two M-gons are triangulated together, a twisted ruled cylinder can result. This needs to be fixed by introducing another subroutine to do this patching. To see this effect, view the default knot at the default parameters of hidesweep with a 135 degree rotation about the y-axis.

This code written by David Dobkin, comments to princeton!dpd@seismo.

DIRECTORY
FS USING [StreamOpen],
IO USING [Close, PutF, real, STREAM],
Real USING [Float,Fix],
RealFns USING [SinDeg, CosDeg,SqRt],
Random USING [RandomStream,Create,NextInt],
Rope USING [ROPE];

KnotPicture: CEDAR PROGRAM

IMPORTS FS, IO, Real, RealFns, Random =

BEGIN OPEN FS, IO, Real, RealFns, Random, Rope;

KnotDescriptorRec: TYPE = RECORD [

 N : INTEGER ← 600, -- Number of samples taken on the curve

 M : INTEGER ← 24,  -- Number of samples taken at each cross-section

 Radius : REAL ← .4, -- Radius of tube passed along the curve (or M-gon at cross-sections).

 A1 : INTEGER ← 4, -- Number of times the knot goes around the center of the torus

 A2 : INTEGER ← 7 , -- Number of times the knot goes around the circumference of the torus

 A3 : REAL ← 2.5 ,  -- Central radius of the tours (average of inner and outer radii)

 FileOutName: Rope.ROPE ← "trefoil.1400.txt", -- Name of file where triangles are to be written

 FileAppend: BOOLEAN ← FALSE, -- Flag tells whether to append to file or create it.

 Partial: BOOLEAN ← FALSE, -- This allows only part of the knot to be generated

 StartAt: BOOLEAN ← TRUE, -- set this flag if initial cross-sections will NOT be generated

 StartAtPoint: INTEGER ← 2, -- First cross-section to be generated

 StopAt: BOOLEAN ← TRUE, -- set this flag if final cross-sections will NOT be generated

 StopAtPoint: INTEGER ← 4 -- Last cross-section to be generated

 ];

KnotDescriptor: TYPE = REF KnotDescriptorRec;

KnotCreate: PROC [ KD : KnotDescriptor ] ~ {

Debug: BOOLEAN ← FALSE;

FileDebugOut: Rope.ROPE ← "debug.txt";

EPSILON : REAL ← .0001;

TWOPI : REAL ← 6.283185 ;

rs : RandomStream ← Create[range:1000000];

Point3: TYPE = ARRAY [1..3] OF REAL;

RowofPoint3s: TYPE = ARRAY [0..100] OF Point3;

FirstRow, PrevRow, ThisRow: RowofPoint3s;

MySinCos: PROC [angle, period : INTEGER] RETURNS [s,c: REAL] ~ {

 a: REAL ← Real.Float[angle]/Real.Float[period];
b: INTEGER ← Real.Fix[a];

 ang:REAL ← 360 * (a - Real.Float[b]);

 s ← RealFns.SinDeg[ang];

 c ← RealFns.CosDeg[ang];

};

GetVals : PROC [i : INTEGER] RETURNS [v,d : Point3] ~ {
-- get values of a point on the curve and the derivative there
 SinTheta1, CosTheta1, SinTheta2, CosTheta2, Deriv1, Deriv2 : REAL;
 [SinTheta1 , CosTheta1 ] ← MySinCos[KD.A1*i,KD.N];
 [SinTheta2 , CosTheta2 ] ← MySinCos[KD.A2*i,KD.N];
 Deriv1 ← Real.Float[KD.A1]*TWOPI/Real.Float[KD.N];
 Deriv2 ← Real.Float[KD.A2]*TWOPI/Real.Float[KD.N]; 
 v[1] ← SinTheta1 * (KD.A3 + CosTheta2);
 v[2] ← CosTheta1 * (KD.A3 + CosTheta2);
 v[3] ← SinTheta2;

 d[1] ← Deriv1*CosTheta1*(KD.A3 + CosTheta2) - SinTheta1*Deriv2*SinTheta2;
 d[2] ← -(Deriv1*SinTheta1*(KD.A3 + CosTheta2) + CosTheta1*Deriv2*SinTheta2);
 d[3] ← Deriv2*CosTheta2;
};

GetFirstFrame : PROC [Dir : Point3] RETURNS [Xd, Yd : Point3] ~ {

-- find axes in first frames, other axes derive from this one.
 Dir ← Normalize[Dir];

  IF (ABS[Dir[3]] > (1.0 - EPSILON)) THEN Xd ← ZCross[Dir]  ELSE Xd ← YCross[Dir];
 Xd ← Normalize[Xd];
 Yd ← Cross[Dir,Xd];
 Yd ← Normalize[Yd];
};

GetFrame : PROC [Dir, Pos, OldPos, OldXd : Point3] RETURNS [Xd, Yd : Point3] ~ {
-- in all succeeding frames, axes are aligned to those in the first frame.

 DelFc : Point3 ← Subtract[Pos , OldPos] ;
 Xd ← Subtract[ OldXd , Scalar[Dot[Dir,OldFx]/Dot[Dir,DelFc] , DelFc] ];
 Xd ← Normalize[Xd];
 Yd ← Cross[Dir,Xd];
 Yd ← Normalize[Yd];
};

MakeRow: PROC [Center, Xd, Yd : Point3] RETURNS [ Retrow : RowofPoint3s] ~ {

-- creates a row of triangles from two M-gonal cross-sections.
CosJ, SinJ : REAL;

FOR J: INTEGER IN [0..KD.M-1] DO

 [SinJ, CosJ] ← MySinCos[J, KD.M];

 Retrow[J] ← Add[Center, Add[ Scalar[ KD.Radius*CosJ, Xd],Scalar[KD.Radius*SinJ, Yd] ] ];

An additional statement has to be added here to apply a transformation matrix.

--At this point, vertices should be output if the ThreeDWorldformat is being used for output.
 ENDLOOP;

};

CopyRow: PROC [OldRow : RowofPoint3s] RETURNS [NewRow : RowofPoint3s] ~ {
FOR J: INTEGER IN [0..KD.M-1] DO NewRow[J] ← OldRow[J]; ENDLOOP;

};

PutTriangle: PROC [p,q,r: Point3,lastone:BOOLEAN] ~ {

RowofTriangles: PROC [ PR, TR : RowofPoint3s, lastone: BOOLEAN] ~ {
FOR JJ:INTEGER IN [0..KD.M-2] DO
PutTriangle[PR[JJ],PR[JJ+1],TR[JJ+1],FALSE];
PutTriangle[PR[JJ],TR[JJ],TR[JJ+1],FALSE];
ENDLOOP;
PutTriangle[PR[KD.M-1],PR[0],TR[0],FALSE];
PutTriangle[PR[KD.M-1],TR[KD.M-1],TR[0],lastone];

};

ZCross: PROC [p: Point3] RETURNS [r: Point3] ~ {r[1] ← p[2]; r[2] ← -p[1];};

YCross: PROC [p: Point3] RETURNS [r: Point3] ~ {r[3] ← p[1]; r[1] ← -p[3];};
 

Scalar: PROC [s:REAL, p: Point3] RETURNS [q: Point3] ~ {

Dot: PROC [p, q: Point3] RETURNS [s: REAL ← 0.0] ~ {

Cross: PROC [p, q: Point3] RETURNS [r: Point3] ~ {

Normalize: PROC [p : Point3] RETURNS [q: Point3] ~ {q ← Scalar[(1.0/RealFns.SqRt[Dot[p,p]]), p]; };

Add: PROC [p, q: Point3] RETURNS [r: Point3] ~ {

 FOR i: CARDINAL IN [1..3] DO r[i] ← p[i] + q[i]; ENDLOOP; };

Subtract: PROC [p, q: Point3] RETURNS [r: Point3] ~ {

 FOR i: CARDINAL IN [1..3] DO r[i] ← p[i] - q[i]; ENDLOOP; };

Fp,Fc,Fx,Fy, OldFx, OldFc, FirstFc : Point3;

ulimit,llimit: INTEGER;

File: IO.STREAM ;

FileDebug:IO.STREAM ;

IF (KD.FileAppend) THEN File ← FS.StreamOpen[KD.FileOutName, $append] ELSE File ← FS.StreamOpen[KD.FileOutName, $create];

IF (Debug) THEN FileDebug ← FS.StreamOpen[FileDebugOut, $create];

IF (KD.Partial AND KD.StartAt) THEN llimit ← KD.StartAtPoint ELSE llimit ← 0;

[Fc , Fp ] ← GetVals[llimit];

[Fx , Fy ] ← GetFirstFrame[Fp];

FirstRow ← MakeRow[Fc,Fx,Fy];

IF (Debug) THEN PutThisRow[Fc,FirstRow];

PrevRow ← CopyRow[FirstRow];

OldFx ← Fx;

FirstFc ← OldFc ← Fc;

IF (KD.Partial AND KD.StopAt) THEN ulimit ← KD.StopAtPoint ELSE ulimit ← KD.N-1;

FOR I : INTEGER IN [llimit+1..ulimit-1] DO
 [Fc , Fp ] ← GetVals[I];
 [Fx , Fy ] ← GetFrame[Fp, Fc, OldFc, OldFx];
 ThisRow ← MakeRow[Fc,Fx,Fy];
 RowofTriangles[PrevRow, ThisRow,FALSE];
PrevRow ← CopyRow[ThisRow];

IF (Debug) THEN PutThisRow[Fc,ThisRow];
 OldFx ← Fx;
 OldFc ← Fc;
 ENDLOOP;

 [Fc, Fp] ← GetVals[ulimit];

 [Fx, Fy] ← GetFrame[Fp, Fc, OldFc, OldFx];

 ThisRow ← MakeRow[Fc,Fx,Fy];

 RowofTriangles[PrevRow, ThisRow,(KD.Partial AND (KD.StartAt OR KD.StopAt))];

 IF (Debug) THEN PutThisRow[Fc,ThisRow];

IF (NOT (KD.Partial AND (KD.StopAt OR KD.StartAt))) THEN RowofTriangles[ThisRow,FirstRow,TRUE];

File.Close[];

IF (Debug) THEN FileDebug.Close[];

};

Kd : KnotDescriptor ← NEW[KnotDescriptorRec];
Kd.FileOutName ← "big.txt";
KnotCreate[ Kd ] ;

END.