perm filename BUG.SAI[S,AIL]2 blob
sn#029656 filedate 1973-03-12 generic text, type C, neo UTF8
COMMENT ⊗ VALID 00013 PAGES
C REC PAGE DESCRIPTION
C00001 00001 VALID 00019 PAGES
C00003 00002 BEGIN "CURVE"
C00011 00003 ⊃ Following are the macros to unpack,packed info in arrays.
C00013 00004 ⊃ OUTER BLOCK STARTS HERE
C00015 00005 SIMPLE INTEGER PROCEDURE INT3DBD(SHORT INTEGER BD,LOW,HIGHREAL ARRAY PLANE
C00021 00006 SIMPLE PROCEDURE BESTPLANE(INTEGER IFRST,ILASTREAL ARRAY X,Y,Z,U,V,GPLANE,COORD)
C00029 00007 SIMPLE PROCEDURE LINEDPY(REAL ARRAY X,Y,Z,ANSINTEGER IFIRST,ILAST)
C00034 00008 PARABOLA
C00037 00009 INTEGER PROCEDURE AXIS (SHORT INTEGER GRPREFERENCE INTEGER BD1,BD2,LOW1F,LOW2F,
C00051 00010 ⊃ PART OF AXIS
C00062 00011 PROCEDURE MK_PIECES
C00067 00012 ⊃ MAIN PART OF OUTER
C00068 00013 MAIN PROGRAM
C00070 ENDMK
C⊗;
BEGIN "CURVE"
COMMENT MODIFICATIONS OF CURV17.SAI,TERMINATION CONDITIONS PUT IN THE AXIS
FINDER ;
REQUIRE "{}{}" DELIMITERS ;
DEFINE CRLF ={'15&'12},LF={'12},LF = {'12},
⊃={COMMENT},
ITYPE(FOO)={BEGIN OUTSTR("FOO"&" "&CVS(FOO)&" ")END},
RTYPE(FOO)={BEGIN OUTSTR("FOO"&" "&CVG(FOO)&" ")END},
OTYPE(FOO)={BEGIN OUTSTR("FOO"&" "&CVOS(FOO)&" ")END},
ATYPE(FOO)={BEGIN OUT(OUTFILE,"FOO"&" "&CVS(FOO)&" ")END},
BTYPE(FOO)={BEGIN OUT(OUTFILE,"FOO"&" "&CVG(FOO)&" ")END},
IOUT(CH,FOO)={BEGIN OUT(CH,CVS(FOO))END},
ROUT(CH,FOO)={BEGIN OUT(CH,CVG(FOO))END},
SIARRAY = { INTEGER ARRAY},
SRARRAY = { REAL ARRAY} ;
DEFINE MAXGROUPS = {25}, MAXBDS = {50} ;
INTEGER DISFINAL,DISPOINTS,DISCONT,FMODE,LMODE,
FIRSTONLY,RUNBEFORE;
REAL PDEN;
INTEGER PPOG,LPOG,FINALPOG,
COUNT,BRCHAR,EOF,FLAG,NPOINTS,NLINES,NNEW,NSEG,INDEX,WORDS,
UPSIGN,RSIGN,SCANSIGN,DEBUG,DISP ;
DEFINE TTY={1},DSKIN={2},DSKOUT={3},PTFILE={4},PLOTFILE={5},OUTFILE= {6},BDFILE={7};
DEFINE P={'17};
DEFINE LOOP={←1 STEP 1 UNTIL};
REQUIRE -1 NEW_ITEMS ;
REQUIRE "DPYSUB.HDR[1,PDQ]" SOURCE_FILE;
COMMENT FOLLOWING DECLARATIONS HAVE COME FROM DPY3D.SAI[1,GJA];
DEFINE CALFILE = {1} , DSKIN = {4} ;
DEFINE PI = {3.1415926},BIG = {1.0@30} ;
REQUIRE "SAITRG[SYS,BGB]" SOURCE_FILE;
REQUIRE "MATRIX.HDR[L,RKN]" SOURCE_FILE ;
STRING FILNAM;
INTERNAL INTEGER MAXPOINTS,MAXLINES,MAXSEGS;
INTERNAL INTEGER ARRAY DPYBUF[1:1000];
INTERNAL INTEGER POG,WHICHFRAME;
INTERNAL REAL DPYMULT,XOFFSET,YOFFSET ;
INTERNAL REAL ARRAY SHIFT,UNSHIFIT,MA,MB[1:4,1:4];
INTERNAL REAL ARRAY ZEROPLANE[1:2,1:4] ;
INTERNAL REAL ARRAY LASLOC,LASAXIS[1:4] ;
INTERNAL INTEGER ARRAY ORIENT,SSIGN[1:2];
INTERNAL REAL ARRAY CAMTRANS[1:4,1:4],CAMVECT[1:8];
INTERNAL REAL STEPANGLE ;
INTERNAL REAL SEGTHRESH;
DEFINE NOTFOUND = {-10000.0},ORTHOGONAL = {-20000.0} ,
FLAT = {-30000.0}, SPHERE = {-40000.0} ,OUTRANGE= {-50000.0} ;
⊃ INTEGER ARRAY PTS[1:MAXPOINTS]; ⊃ NOT USED ;
⊃ INTEGER ARRAY MID,BASSIZ[1:MAXSEGS]; ⊃ PTS AND BASSIZ NOT USED ANY MORE ;
⊃ PTS has coords,packed two to a word , OFFSET + 200.
BASSIZ - Packed array.Left half contains BASE(accessed by a macro).
BASE points to first location in PTS for the segment.
Right half contains SIZSEG,the no of words of PTS taken up
by this segment.
MID - Midpoint of segment, X in left half,Y in the right half.
Stored values are just the sum of XBEG and YBEG.To get
the right value,divide by two and subtract the offset 200.
LINNO - line no of the segment.
BEGEND- Coordinates of the first and last point of the segment.
Info pac
⊂!eft half contains beginning.First nine bits
contain x coord,next nine y coord.OFFSET 200.
FLSEG- Left half contains serial no of first segment in this line.
Right half contains last segment.
GRPNO- Group that this seg belongs to.-1 indicates does not belong to any.
0 indicates seg not present.(larger than SEGS2).
BASGRP- Tells where the first segment in this group is located in GRPLST.
(Gives index in GRPLST).
GRPLST- Linear list of segments.A group is contiguos no of members.No is
availaible in GRPSIZ.Starting point in BASGRP.Note some of the
entries may not belong to any group.
LINE1,SEGS1- Total no of the lines and segments in first scan.
LINE2,SEGS2- Total no of the lines and segments in both scans. ;
⊃ REQUIRE "SOBMAT[SYS,HE]" LOAD_MODULE ;
⊃ EXTERNAL PROCEDURE QSOLVE (REAL ARRAY A,X,B);
REQUIRE "FILIN[1,RKN]" LOAD_MODULE ;
EXTERNAL SIMPLE PROCEDURE FILIN(REFERENCE STRING FILE);
EXTERNAL SIMPLE PROCEDURE GET_ARBDS; ⊃ PROCEDURES TO DEFINE MAX BOUNDS FOR PTFILE;
⊃ Following are the macros to unpack,packed info in arrays.;
DEFINE
BASE(ISEG) = {(BASIZ[ISEG] LSH -18) } ,
SIZSEG(ISEG) = { (BASIZ[ISEG] LAND '777777 ) },
FIRSTSEG(LINE) = { (FLSEG[LINE] LSH - 18) } ,
LASTSEG(LINE) = { (FLSEG[LINE] LAND '777777 ) } ,
XBEG(SEG)= {( (BEGEND[SEG] LSH -27) - 200) },
YBEG(SEG)= {(((BEGEND[SEG] LSH -18) LAND '777) - 200 )},
XEND(SEG)= {(((BEGEND[SEG] LAND '777777) LSH -9)- 200) },
YEND(SEG)= {( (BEGEND[SEG] LAND '777) - 200 )} ,
MIDX(SEG)= {( MID[SEG] LSH -18 )} ,
MIDY(SEG)= {( MID[SEG] LAND '777777 ) } ,
BDNOBEG(SEG) = { (BDNO[SEG] LSH -18 )} ,
BDNOEND(SEG) = { (BDNO[SEG] LAND '777777 )} ,
BDLOCBEG(SEG) = { (BDLOC[SEG] LSH -18 )} ,
BDLOCEND(SEG) = { (BDLOC[SEG] LAND '777777 )} ,
GET_BDNO(SEG) = { (IF SEG>0 THEN BDNOBEG(SEG) ELSE BDNOEND(-SEG)) },
GET_BDLOC(SEG)= { (IF SEG>0 THEN BDLOCBEG(SEG) ELSE BDLOCEND(-SEG)) },
BDSEG(LOC)= { ( (BDARRY[LOC] LAND '777777) - '100000 ) } ,
BDMARK(LOC) = { (BDARRY[LOC] LSH -18) } ;
⊃ OUTER BLOCK STARTS HERE ;
PROCEDURE OUTER ;
BEGIN "OUTER"
INTERNAL INTEGER ARRAY LINNO,BEGEND,GRPNO,GRPLST[1:MAXSEGS];
INTERNAL INTEGER ARRAY FLSEG,MIRANG[1:MAXLINES] ;
INTERNAL REAL ARRAY COLLIN[1:MAXLINES,1:4,1:3];
INTERNAL INTEGER ARRAY GRPSIZ,BASGRP[1:MAXGROUPS];
INTERNAL REAL ARRAY WIDTH[1:MAXSEGS];
INTERNAL INTEGER LINE1,LINE2,SEGS1,SEGS2;
INTERNAL INTEGER ARRAY BASBD,SIZBD[1:MAXBDS],BDARRY[1:2*MAXSEGS];
INTERNAL INTEGER ARRAY BDNO,BDLOC[1:MAXSEGS];
INTERNAL INTEGER NUMBDS; ⊃ Actual no of boundary segments found ;
REQUIRE "CUR19A[1,RKN]" LOAD_MODULE ;
EXTERNAL PROCEDURE READPTS(STRING FILE);
EXTERNAL SIMPLE PROCEDURE READCAL ;
EXTERNAL PROCEDURE BDIN ;
EXTERNAL SIMPLE PROCEDURE DPYSEG(INTEGER ISEG);
EXTERNAL PROCEDURE DISPLAY ;
EXTERNAL SIMPLE PROCEDURE GROUP ;
EXTERNAL SIMPLE PROCEDURE GETXYBD(INTEGER BD,INDX;REFERENCE REAL X,Y);
EXTERNAL SIMPLE PROCEDURE THREED(INTEGER LINE;REAL U,V;REFERENCE REAL X,Y,Z);
EXTERNAL SIMPLE PROCEDURE TWOD(REAL X,Y,Z;REFERENCE REAL U,V);
EXTERNAL PROCEDURE GET_BREAKS ;
SIMPLE INTEGER PROCEDURE INT3DBD(SHORT INTEGER BD,LOW,HIGH;REAL ARRAY PLANE ;
REFERENCE REAL XINT,YINT,ZINT;REFERENCE INTEGER BDINDX);
BEGIN "INT3DBD"
⊃ Finds intersection of given plane with the specified boundary segment.;
⊃ This one for 3d intersections ;
SHORT REAL MD,FD,LD,FX,FY,FZ,LX,LY,LZ,MX,MY,MZ,A,B,C,D,DIFF,UINT,VINT;
SHORT INTEGER FU,FV,LU,LV,MU,MV,MID,DEBUG,BASB,FIRST,LAST,LINE,SEG1,SEG2,SIZB;
BOOLEAN INTDEB ;
LABEL L1,L2,L3;
DEFINE NEAR_LOW = {0}, NEAR_HIGH = {1}, FOUND = {-1} ;
⊃ Finds the coords of point in the current segment.You set BASSEG before
calling it.;
SIMPLE PROCEDURE COORD3(SHORT INTEGER INDX ; REFERENCE REAL X,Y,Z);
BEGIN "COORD3"
INTEGER SEG;
SEG← BDSEG(BASB+((INDX-1+SIZB) MOD SIZB)) ;
⊃ The term INDX-1+SIZB is to account for
a. negative INDX by adding SIZB to it.
b.-1 is so that INDX= SIZB does not give zero. ;
IF SEG> 0 THEN
THREED(LINE← LINNO[SEG],XBEG(SEG),YBEG(SEG),X,Y,Z)
ELSE
THREED(LINE← LINNO[-SEG],XEND(-SEG),YEND(-SEG),X,Y,Z)
END "COORD3";
IF HIGH= LOW THEN RETURN(0);
A← PLANE[1]; B←PLANE[2]; C← PLANE[3]; D← PLANE[4];
BASB← BASBD[BD];
SIZB← SIZBD[BD];
FIRST← LOW;LAST← HIGH;
COORD3(FIRST,FX,FY,FZ);
COORD3(LAST,LX,LY,LZ);
FD← A* FX + B*FY +C*FZ + D ;
LD← A* LX + B*LY +C*LZ + D ;
IF INTDEB THEN BEGIN
RTYPE(FX);RTYPE(FY);RTYPE(FZ);RTYPE(LX);RTYPE(LY);RTYPE(LZ);
RTYPE(FD);RTYPE(LD);OUTSTR(CRLF);END;
IF FD * LD ≤0 THEN BEGIN
⊃ BINARY SEARCH WILL WORK HERE;
L1: WHILE FIRST< LAST -1 DO BEGIN
DIFF← LAST - FIRST ; ⊃ So that DIFF may be FSCed by compiler
rather than FLOATed. ;
MID← FIRST + ABS(DIFF*FD/(FD - LD)) ;⊃ Note the minus sign
in (FD - LD).Eliminates taking absolute values since
FD and LD are both known to be of opposite sign.;
IF MID = FIRST THEN MID←MID +1;
IF MID = LAST THEN MID←LAST-1;
COORD3(MID,MX,MY,MZ);
MD← A* MX + B*MY + C*MZ + D ;
IF FD * MD > 0 THEN BEGIN
FIRST←MID;FD←MD;END
ELSE BEGIN
LAST←MID;LD←MD;END;
IF INTDEB THEN BEGIN
RTYPE(FD);RTYPE(LD);RTYPE(MD);OUTSTR(CRLF);
RTYPE(FIRST);RTYPE(LAST);RTYPE(MID);OUTSTR(CRLF);
END;
END; ⊃ WHILE FIRST<LAST ;
END
ELSE BEGIN
IF (HIGH-LOW)>6 THEN RETURN(IF ABS(FD)< ABS(LD) THEN NEAR_LOW ELSE NEAR_HIGH);
⊃ Both ends on the same side of the line.No intersection.;
⊃ SEARCH STEP BY STEP,IF LESS THEN FIVE POINTS ;
L2: FOR LAST← LOW +1 STEP 1 UNTIL HIGH DO BEGIN
COORD3(LAST,LX,LY,LZ);
LD← A* LX + B* LY + C*LZ +D;
IF FD*LD≤0 THEN DONE ELSE BEGIN
IF LAST=HIGH THEN RETURN(0);
⊃ ********* NOTE IT DOES NOT RETURN WHICH END IT'S NEAR ;
FIRST← LAST; FD← LD ;
END;
IF INTDEB THEN BEGIN
RTYPE(FD);RTYPE(LD);OUTSTR(CRLF);
RTYPE(FIRST);RTYPE(LAST);OUTSTR(CRLF);
END;
END; ⊃ FOR LAST← ;
END; ⊃ ELSE BEGIN ;
L3: COORD3(FIRST,FX,FY,FZ);
COORD3(LAST,LX,LY,LZ);
FD←ABS(FD);LD←ABS(LD);
⊃ Find intersection by linear interpolation between the two points.;
DIFF← LX - FX ;
XINT←( FX+ (DIFF* FD)/(FD + LD) );
DIFF← LY - FY ;
YINT←( FY+ (DIFF* FD)/(FD + LD) );
DIFF← LZ - FZ ;
ZINT←( FZ+ (DIFF* FD)/(FD + LD) );
BDINDX← IF ABS(FD)<ABS(LD) THEN FIRST ELSE LAST ;
IF INTDEB THEN BEGIN
RTYPE(XINT); RTYPE(YINT);RTYPE(ZINT);
OUTSTR(CRLF);END;
IF DISP THEN BEGIN
TWOD(XINT,YINT,ZINT,UINT,VINT);
APOINT(UINT*DPYMULT +XOFFSET,VINT*DPYMULT+YOFFSET);
DPYOUT(1) ;
END;
RETURN(FOUND);
END "INT3DBD";
SIMPLE PROCEDURE BESTPLANE(INTEGER IFRST,ILAST;REAL ARRAY X,Y,Z,U,V,GPLANE,COORD);
BEGIN "BESTPLANE"
COMMENT FIND THE BEST PLANE FITTING THE POINTS IN X,Y,Z[IFRST:ILAST]
AND RELOCATE THOSE POINTS TO THE PLANE, PLACING THE RELOCATED POINTS IN
U,V[IFRST:ILAST].
SEE NOTEBOOK, VOLUME 11, PAGES 90-92;
INTEGER I,J,N,I1,I2,I3; REAL X1,Y1,Z1,D ;
OWN REAL ARRAY PSUMS[1:4,1:4];
MZERO(PSUMS);
FOR I←IFRST STEP 1 UNTIL ILAST DO BEGIN "BUILD PSUMS MATRIX"
PSUMS[1,4]←PSUMS[1,4]+(X1←X[I]);
PSUMS[2,4]←PSUMS[2,4]+(Y1←Y[I]);
PSUMS[3,4]←PSUMS[3,4]+(Z1←Z[I]);
PSUMS[1,1]←PSUMS[1,1]+X1*X1;
PSUMS[1,2]←PSUMS[1,2]+X1*Y1;
PSUMS[1,3]←PSUMS[1,3]+X1*Z1;
PSUMS[2,2]←PSUMS[2,2]+Y1*Y1;
PSUMS[2,3]←PSUMS[2,3]+Y1*Z1;
PSUMS[3,3]←PSUMS[3,3]+Z1*Z1;
END;
PSUMS[4,4]←N←ILAST-IFRST+1;
COMMENT TRANSLATE CENTER OF GRAVITY TO ORIGIN OF COORDINATES;
FOR I LOOP 3 DO BEGIN
PSUMS[I,I]←PSUMS[I,I]-PSUMS[I,4]↑2/N;
FOR J←I+1 STEP 1 UNTIL 3 DO
PSUMS[J,I]←PSUMS[I,J]←PSUMS[I,J]-PSUMS[I,4]*PSUMS[J,4]/N;
END;
COMMENT GUESS AT ORIENTATION FOR HOMOGENEOUS SOLUTION;
I1←IF PSUMS[1,1]<PSUMS[2,2] THEN 1 ELSE 2;
I1←IF PSUMS[I1,I1]<PSUMS[3,3] THEN I1 ELSE 3;
I2←I1 MOD 3 + 1;
I3←I2 MOD 3 + 1;
GPLANE[I1]←1.0;
COMMENT HOMOGENEOUS SOLUTION;
D←PSUMS[I2,I2]*PSUMS[I3,I3] - PSUMS[I2,I3]↑2;
GPLANE[I2]←-(PSUMS[I1,I2]*PSUMS[I3,I3] - PSUMS[I1,I3]*PSUMS[I2,I3])/D;
GPLANE[I3]←-(PSUMS[I1,I3]*PSUMS[I2,I2] - PSUMS[I1,I2]*PSUMS[I2,I3])/D;
COMMENT NORMALIZE;
D←SQRT(GPLANE[1]↑2 + GPLANE[2]↑2 + GPLANE[3]↑2);
FOR I LOOP 3 DO GPLANE[I]←GPLANE[I]/D;
GPLANE[4]←-(GPLANE[1]*PSUMS[1,4] + GPLANE[2]*PSUMS[2,4]
+GPLANE[3]*PSUMS[3,4])/N;
COMMENT FORM COORDINATION MATRIX;
MCOORD(COORD,GPLANE);
COMMENT GET PROJECTION ON GPLANE ACCORDING TO INVERSE
COORDINATION MATRIX;
FOR I←IFRST STEP 1 UNTIL ILAST DO BEGIN "PROJECT"
U[I]←COORD[1,1]*X[I] COMMENT COORD[2,1] IS ZERO;
+ COORD[3,1]*Z[I];
V[I]←COORD[1,2]*X[I] + COORD[2,2]*Y[I]
+ COORD[3,2]*Z[I];
COMMENT THE SEGNO MATRIX WILL BE GARBAGE;
END;
END "BESTPLANE";
SIMPLE REAL PROCEDURE LINEFIT(INTEGER IFIRST,ILAST;REAL ARRAY U,V,ANS);
BEGIN "LINEFIT"
COMMENT COMPUTE LOW ORDER SUMS
PERFORM LEAST SQUARES ;
REAL X,Y,Z,ERROR,SUMX,SUMXX,SUMY,SUMXY,SUMYY;
INTEGER I,J,DX,DY,BADPOINT,CONSEC,MAXOUT,N ;
REAL A,B,C,LAMBDA,TEMP,BADVAL,LAST;
N← ILAST-IFIRST + 1 ;
COMMENT COMPUTE LOW ORDER SUMS;
SUMX←SUMXX←SUMY←SUMXY←SUMYY←0;
FOR I←IFIRST STEP 1 UNTIL ILAST DO BEGIN
SUMX←SUMX+(X←U[I]);
SUMY←SUMY+(Y←V[I]);
SUMXX←SUMXX+X*X;
SUMYY←SUMYY+Y*Y;
SUMXY←SUMXY+X*Y END;
COMMENT LEAST SQUARES COMPUTATION;
A←N*SUMXX-SUMX*SUMX;
B←N*SUMXY-SUMX*SUMY;
C←N*SUMYY-SUMY*SUMY;
IF B=0 THEN BEGIN
IF A>C THEN X←0 ELSE X←1.0;
Y←1.0-X END
ELSE BEGIN
LAMBDA←(A+C-((A-C)↑2+4*B*B)↑.5)/2;
X←-B;
Y←A-LAMBDA;
A←(X↑2+Y↑2)↑.5;
X←X/A;Y←Y/A END;
ANS[1]←X;
ANS[2]←Y;
ANS[3]←Z←-(SUMX*X+SUMY*Y)/N;
COMMENT OVERALL FIT;
A←SUMXX*X+SUMXY*Y+SUMX*Z;
B←SUMXY*X+SUMYY*Y+SUMY*Z;
C←SUMX*X+SUMY*Y+N*Z;
ERROR←((A*X+B*Y+C*Z)/N)↑.5;
RETURN(ERROR);
END "LINEFIT" ;
PROCEDURE LINE3D(INTEGER IFIRST,ILAST;REAL ARRAY X,Y,Z,ANS);
BEGIN "LINE3D"
INTEGER I,J,K ;
REAL DEN ;
REAL ARRAY U,V[IFIRST:ILAST];
OWN REAL ARRAY GPLANE[1:4],COORD[1:4,1:4],XL,YL,ZL,UL,VL[1:2],ANS1[1:3];
BESTPLANE(IFIRST,ILAST,X,Y,Z,U,V,GPLANE,COORD);
LINEFIT(IFIRST,ILAST,U,V,ANS1);
⊃ NOW WE HAVE THE TWOD LINE,CONVERT TO THREE-D PARAMETERS ;
IF ABS(ANS1[1])>0.5 THEN BEGIN
VL[1]← 1.0 ; VL[2]← 10.0 ;
FOR I LOOP 2 DO
UL[I]← -(ANS1[2]*VL[I] + ANS1[3])/ANS1[1] ;
END
ELSE BEGIN
UL[1]← 1.0 ; UL[2]← 10.0 ;
FOR I LOOP 2 DO
VL[I]← -(ANS1[1]* UL[I] + ANS1[3])/ANS1[2] ;
END;
FOR I LOOP 2 DO BEGIN
XL[I]← COORD[1,1]* UL[I] + COORD[1,2]* VL[I] - COORD[1,3]*GPLANE[4];
⊃ THE W COORDINATE FOR ALL POINTS IS THE DIST.
OF THE PLANE FROM THE ORIGIN (- GPLANE[4]);
YL[I]← COORD[2,1]* UL[I] + COORD[2,2]* VL[I] - COORD[2,3]*GPLANE[4];
ZL[I]← COORD[3,1]* UL[I] + COORD[3,2]* VL[I] - COORD[3,3]*GPLANE[4];
END;
ANS[1]← XL[2]- XL[1] ; ⊃ THE THREE DIRECTION COSINES ;
ANS[2]← YL[2]- YL[1] ;
ANS[3]← ZL[2]- ZL[1] ;
⊃ NORMALIZE THE COEFFICIENTS ;
DEN← SQRT(ANS[1]↑2 + ANS[2]↑2 + ANS[3]↑2 );
FOR I LOOP 3 DO ANS[I]← ANS[I]/DEN ;
ANS[4]← 1 ;
ANS[5]← XL[1] ; ⊃ THE CO-ORDINATES OF ONE POINT ON THE LINE ;
ANS[6]← YL[1] ;
ANS[7]← ZL[1] ;
ANS[8]← IFIRST ;ANS[9]← ILAST ; ⊃ STORE THE INDICES ;
END "LINE3D" ;
SIMPLE PROCEDURE LINEDPY(REAL ARRAY X,Y,Z,ANS;INTEGER IFIRST,ILAST);
BEGIN "LINEDPY"
INTEGER I,J,K,I1,I2,I3 ;
OWN REAL ARRAY XYZF,XYZL,F,L [1:3] ;
REAL DX,DY ;
DPYSET(DPYBUF);
⊃ PICK THE LARGEST DIRECTION COSINE OF THE LINE ;
I1← IF ABS(ANS[1])> ABS(ANS[2]) THEN 1 ELSE 2 ;
I1← IF ABS(ANS[I1])> ABS(ANS[3]) THEN I1 ELSE 3 ;
⊃ PUT THE LAST AND FIRST POINTS IN F AND L ;
F[1]← X[IFIRST];F[2]← Y[IFIRST]; F[3]← Z[IFIRST];
L[1]← X[ILAST]; L[2]← Y[ILAST]; L[3]← Z[ILAST];
XYZF[I1]← F[I1];
XYZL[I1]← L[I1];
I2← I1 ;
FOR I LOOP 2 DO BEGIN
I2← I2 MOD 3 + 1 ;
XYZF[I2]← ANS[I2+4] + ANS[I2]*(XYZF[I1] - ANS[I1+4])/ANS[I1] ;
XYZL[I2]← ANS[I2+4] + ANS[I2]*(XYZL[I1] - ANS[I1+4])/ANS[I1] ;
END;
TWOD(XYZF[1],XYZF[2],XYZF[3],DX,DY);
APOINT(DX*DPYMULT + XOFFSET,DY*DPYMULT + YOFFSET) ;
TWOD(XYZL[1],XYZL[2],XYZL[3],DX,DY);
AVECT(DX*DPYMULT + XOFFSET,DY*DPYMULT + YOFFSET) ;
DPYOUT(2);
END "LINEDPY" ;
SIMPLE PROCEDURE NORMAL(REAL ARRAY LINE,PLANE;REAL XP,YP,ZP);
BEGIN "NORMAL"
INTEGER I,J,K,I1,I2,I3 ;
OWN REAL ARRAY XYZ,F [1:3] ;
⊃ PICK THE LARGEST DIRECTION COSINE OF THE LINE ;
I1← IF ABS(LINE[1])> ABS(LINE[2]) THEN 1 ELSE 2 ;
I1← IF ABS(LINE[I1])> ABS(LINE[3]) THEN I1 ELSE 3 ;
F[1]← XP;F[2]← YP; F[3]← ZP ;
XYZ[I1]← F[I1];
I2← I1 ;
FOR I LOOP 2 DO BEGIN
I2← I2 MOD 3 + 1 ;
XYZ[I2]← LINE[I2+4] + LINE[I2]*(XYZ[I1] - LINE[I1+4])/LINE[I1] ;
END;
FOR I LOOP 3 DO PLANE[I]← LINE[I]; ⊃ THE DIRECTION IS CORRECT ;
⊃ NOW COMPUTE THE DISTANCE FROM THE ORIGIN ;
PLANE[4]← -(XYZ[1]*PLANE[1] + XYZ[2]* PLANE[2] + XYZ[3]* PLANE[3] ) ;
END "NORMAL" ;
SIMPLE PROCEDURE LINEEXP(REAL ARRAY LINE,X,Y,Z,RAD;
INTEGER N ,STEPSIZ ; REFERENCE REAL XNP,YNP,ZNP);
⊃ This version extrapolates to a distance proportional to the radius function;
BEGIN "LINEEXP"
INTEGER I,J,K,I1,I2,I3,N1,N2,PREVN ;
OWN REAL ARRAY XYZ,XYZF,XYZL,F,L,LP [1:3] ;
SHORT REAL R1 ,SIGN,INCR,MINSTEP;
DEFINE STEPMULT = {0.05} ;
MINSTEP← 0.2 ;
⊃ PICK THE LARGEST DIRECTION COSINE OF THE LINE ;
I1← IF ABS(LINE[1])> ABS(LINE[2]) THEN 1 ELSE 2 ;
I1← IF ABS(LINE[I1])> ABS(LINE[3]) THEN I1 ELSE 3 ;
N1← LINE[8] ; N2← LINE[9]; ⊃ INDICES OF END POINTS ;
F[1]← X[N1];F[2]← Y[N1]; F[3]← Z[N1];
L[1]← X[N2];L[2]← Y[N2]; L[3]← Z[N2];
SIGN← IF (L[I1]-F[I1])> 0 THEN 1.0 ELSE -1.0 ;
PREVN← N - STEPSIZ ;
LP[1]← X[PREVN]; LP[2]← Y[PREVN]; LP[3]← Z[PREVN] ; ⊃ COORDINATES OF THE LAST POINT ;
INCR← ABS(STEPMULT * RAD[PREVN] *STEPSIZ) MAX MINSTEP ;
IF STEPSIZ<0 THEN INCR← -INCR ;
XYZ[I1]← LP[I1] + INCR * SIGN ;
I2← I1 ;
FOR I LOOP 2 DO BEGIN
I2← I2 MOD 3 + 1 ;
XYZ[I2]← LINE[I2+4] + LINE[I2]*(XYZ[I1] - LINE[I1+4])/LINE[I1] ;
END;
XNP← XYZ[1] ;YNP← XYZ[2]; ZNP← XYZ[3] ;
END "LINEEXP" ;
COMMENT PARABOLA;
REAL SIMPLE PROCEDURE APARABOLA(SHORT INTEGER INDXL,INDXH;REAL ARRAY RAD;
REFERENCE REAL A,B,C;REFERENCE INTEGER OFFSET);
BEGIN "APARABOLA"
COMMENT FINDS THE LEAST SQUARES FIT OF THE WIDTHS OF SEGMENTS IN
GROUP TO THE EQUATION WIDTH = A*(X-OFFSET)↑2 + B*(X-OFFSET) + C
SEE NOTEBOOK, VOL 11, PAGES 89-90;
REAL N,N2,N3,N4,N5,N6,N7,N8,N9,D,Y,XY,XXY,YY,E,F,G;
INTEGER I,J;
OFFSET← INDXL -1 ;
COMMENT GET N MATRIX;
N←INDXH - INDXL +1 ;
IF N≤2 THEN BEGIN
⊃ Not enough points to fit a parabola,so fit a straight line ;
A← 0 ; C← RAD[INDXL];
B← (RAD[INDXH] - C);
RETURN(0); ⊃ No error ;
END;
N2←N*N; N3←N2*N; N4←N3*N; N5←N4*N; N6←N5*N; N7←N6*N; N8←N7*N; N9←N8*N;
COMMENT GET THE Y MATRIX;
Y←XY←XXY←YY←0; J←0;
FOR I ← INDXL STEP 1 UNTIL INDXH DO BEGIN
D←RAD[I];
J←J+1;
Y←Y+D;
XY←XY+D*J;
XXY←XXY+D*J*J;
YY←YY+D*D;
END;
COMMENT SOLVE THE EQUATIONS;
D ← N9 - 6*N7 + 9*N5 - 4*N3;
A ← (XXY * 180*(N4 - N2)
+ XY * (E ← 180*(-N5 - N4 + N3 + N2))
+ Y * (F ← 30*(N6 + 3*N5 + N4 - 3*N3 - 2*N2)))/D;
B ← (XXY * E
+ XY * 12*(16*N6 + 30*N5 - 5*N4 - 30*N3 - 11*N2)
+ Y * (G ← 18*(-2*N7 - 7*N6 - 5*N5 + 5*N4 + 7*N3 + 2*N2)))/D;
C ← (XXY * F
+ XY * G
+ Y * 3*(3*N8 + 12*N7 + 14*N6 - 13*N4 - 12*N3 - 4*N2) )/D;
COMMENT COMPUTE THE ERROR;
RETURN (( A↑2 * (6*N5 + 15*N4 + 10*N3 - N) / 30
+ 2 * A * B * (N4 + 2*N3 + N2) / 4
+ (B↑2 + 2*A*C) * (2*N3 + 3*N2 + N) / 6
+ 2 * B * C * (N2 + N) / 2
+ C↑2 * N
- 2*A*XXY - 2*B*XY - 2*C*Y + YY) / N);
END "APARABOLA";
INTEGER PROCEDURE AXIS (SHORT INTEGER GRP;REFERENCE INTEGER BD1,BD2,LOW1F,LOW2F,
HIGH1F,HIGH2F,BSIGN,PARB,PARE,TERMB,TERME,CIRC;REAL ARRAY AXX,AXY,AXZ,RAD,ANGLE);
BEGIN "AXIS"
SHORT REAL X1,Y1,Z1,X2,Y2,Z2,XN,YN,ZN,XN1,YN1,ZN1,XN2,YN2,ZN2,COR,XCOR,YCOR,ZCOR,
DX,DY,U1,V1,U2,V2,CORTHRESH ,ERR,XNP,YNP,ZNP,A,B,C,RADERR,RADTHRESH;
OWN REAL ARRAY PLANE[1:4],F,L[1:3],AXLINE,INITLIN[1:9] ;
SHORT INTEGER I,J,K,SIZG,BASG,MIDSEG,N1,N2,N11,N22,NM,SEG1,SEG2,SEGM,ISEG,LINE,
SIZE1,SIZE2,FITSIZ ,M,I1,PAR,PAR1,PAR2,IPAR,PARBEG,PAREND,
STEPSIZ,INITNO ,IFIRST,ILAST,I2,I3,BSIGN1,BSIGN2,LEN1,LEN2,OFFSET,
LOW1,HIGH1,LOW2,HIGH2,MID1,MID2 ,
LOW11,HIGH11,LOW22,HIGH22,LOW10,LOW20,HIGH10,HIGH20,INDX1,INDX2,
HSIZE1,HSIZE2,TEMP,SIGN,FINDX,LINDX,LEN ;
BOOLEAN AXDEB,TTEST ;
LABEL L1;
DEFINE TOP= "1" , BOTTOM = "2" ;
DEFINE NEAR_LOW = {0},NEAR_HIGH = {1}, FOUND = {-1} ;
DEFINE EXP_RADIUS(PAR) = {(A*(PAR-OFFSET)↑2 + B*(PAR-OFFSET) + C) };
DEFINE BDDIST(F,L,SIZE)= {(IF (TEMP←ABS(F-L) MOD SIZE)<(SIZE/2)
THEN (TEMP-1) ELSE (SIZE-TEMP-1))} ;
REAL SIMPLE PROCEDURE DIST_LINE(INTEGER BD,INDX;REAL ARRAY LINE);
BEGIN "DIST_LINE"
SHORT REAL X,Y;
GETXYBD(BD,INDX,X,Y);
RETURN(LINE[1]*X + LINE[2]*Y + LINE[3]);
END "DIST_LINE";
SIMPLE BOOLEAN PROCEDURE TERMTEST(INTEGER BD,INDX1,INDX2,SIZE,SIGN;REAL THRESH;
REFERENCE INTEGER LEN;REAL ARRAY LINE);
BEGIN "TERMTEST"
⊃ TEST WHETHER BOUNDARY EXTNDS FAR ENOUGH,BETWEEN INDX1 AND INDX2(MODULO
SIZE) IN THE DIRECTION SIGN;
SHORT INTEGER I,SSIZE,NSTEPS,INDX ;
SHORT REAL DIST ;
LEN← ABS(INDX1 - INDX2) MOD SIZE ;
IF LEN> (SIZE/2) THEN LEN← SIZE - LEN ;
IF LEN≤1 THEN RETURN(TRUE); ⊃ OBVIOUSLY NEAR BOUNDARY;
SSIZE← SIGN * (5 MIN (LEN1/5)); ⊃ CHOOSE AN APPROPRIATE SIZE ;
NSTEPS← LEN1/SSIZE ;
INDX←INDX1;
FOR I LOOP NSTEPS DO BEGIN
DIST ← DIST_LINE(BD,INDX,LINE);
IF ABS(DIST)>(5.0 MAX THRESH) THEN RETURN(FALSE);
INDX←INDX + NSTEPS ;
END;
RETURN(TRUE);
END "TERMTEST";
BOOLEAN PROCEDURE RADCHECK(INTEGER PAR,PAR1,PAR2,STEPSIZ;REAL ARRAY RAD);
BEGIN "RADCHECK"
⊃ CHECK TO SEE IF THE PRESENT RADIUS FUNCTION IS ACCEPTABLE;
REAL A,B,C,PRAD,LRAD,RADERR,RADTHRESH,DIFF,ERAD ;
SHORT INTEGER F,L,OFFSET;
DIFF← ABS( (PRAD←RAD[PAR])-(LRAD←RAD[PAR-STEPSIZ]));
IF DIFF< (0.15*LRAD MAX 2.0) THEN RETURN(TRUE) ; ⊃ ITS OK,KEEP IT;
IF DIFF> (LRAD MAX 4.0) THEN RETURN(FALSE);⊃ CHANGE IS TOO LARGE,NEED NOT CHECK FURTHER;
⊃ NOW CHECK FOR MORE GLOBAL EVALUATIONS;
⊃ FIRST FIT A PARABOLA TO THE PREVIOUS POINT ;
IF STEPSIZ>0 THEN BEGIN F←PAR1; L←PAR2-STEPSIZ ;END
ELSE BEGIN F←PAR1-STEPSIZ ;L←PAR2; END;
APARABOLA(F,L,RAD,A,B,C,OFFSET);
ERAD← EXP_RADIUS(PAR);
RETURN(IF ABS(RADERR←(ERAD- PRAD))> (RADTHRESH← (0.25*ERAD MAX 5.0)) THEN
FALSE ELSE TRUE);
END "RADCHECK";
SIMPLE BOOLEAN PROCEDURE ISTERM(SHORT INTEGER BD,INDX1,INDX2,F,L,SIGN);
BEGIN "ISTERM"
⊃ DETERMINE WHETHER THE SEGMENT FROM F TO L OF THE BD CAN BE CONSIDERED AS
TERMINATION,INDX1,INDX2 THE INDICES OF END-POINTS;
SHORT INTEGER INDX,SIZB ;
SHORT REAL X1,Y1,X2,Y2,X,Y,XN,YN,DIST,T1,T2,T3,T4,T5,T6,T7,THETA,XT,YT;
DEFINE ATHRESH = {10.0}, XTHRESH = {7.0} ; ⊃ THRESHHOLDS ALONG AND NORMAL TO AXIS;
DEFINE EDDIST(INDX) =
{ BEGIN
GETXYBD(BD,INDX,X,Y);
XN← T1* (XT← X- T5) + T2 * (YT← Y - T6);
YN← ABS(T3* XT + T4 * YT) ;
XN← IF XN< 0 THEN ABS(XN) ELSE IF XN≤ T7 THEN 0.0 ELSE XN - T7 ;
⊃ SEE WHETHER BETWEEN THE SEGMENT OR DIST FROM THE NEAR END;
END; } ;
SIZB← SIZBD[BD];
⊃ FIRST FORM A TRANSFORMATION TO ALIGN COORDINATES WITH EDGE JOINING INDX1
AND INDX2 ;
GETXYBD(BD,INDX1,X1,Y1);
GETXYBD(BD,INDX2,X2,Y2);
THETA← ATAN2(Y2-Y1,X2-X1);
T1← T4← COS(THETA);
T3← -(T2← SIN(THETA)); ⊃ COEFFICIENTS OF THE TRANSFORMATION;
T5← X1; T6← Y1; ⊃ COORDINATES OF THE NEW ORIGIN;
T7← SQRT((X1- X2)↑2 + (Y1-Y2)↑2); ⊃ LENGTH OF THE EDGE ;
INDX← F;
WHILE TRUE DO BEGIN
EDDIST(INDX,DIST); ⊃ THE DISTANCE FROM EDGE ;
IF XN > XTHRESH OR YN> ATHRESH THEN RETURN(FALSE);
INDX← INDX + SIGN ;
IF (ABS(INDX- L) MOD SIZB ) = 0 THEN DONE ;⊃ L=INDX MODULO SIZB ;
END;
RETURN(TRUE); ⊃ NONE OF THE DIST WAS LARGER THEN THRESH ;
END "ISTERM";
DPYSET(DPYBUF);
SIZG←GRPSIZ[GRP];
BASG←BASGRP[GRP];
N11←1; N22← SIZG;
⊃ Find parameters to relate to the boundaries ;
SEG1← GRPLST[BASG+N11-1];
SEG2← GRPLST[BASG+N22-1];
⊃ First find the boundary corresponding to SEG1 ;
WHILE TRUE DO BEGIN
IF N22-N11 <3 THEN BEGIN
OUTSTR("THIS GROUP NOT IN BOUNDARY LIST,I QUIT"&CRLF);
RETURN(0); END;
BD1← BDNOBEG(SEG1);BD2← BDNOEND(SEG1);⊃ Pick the opposite boundaries;
IF BD1=0 OR BD2=0 THEN BEGIN
N11←N11+1;
SEG1← GRPLST[BASG+N11-1];END
ELSE DONE;
END;
SIZE1← SIZBD[BD1];SIZE2← SIZBD[BD2];
HSIZE1← SIZE1/2; HSIZE2← SIZE2/2; ⊃ Half the size ;
LOW1← BDLOCBEG(SEG1); LOW2← BDLOCEND(SEG1);
⊃ FIND THE BOUNDARY SEGMENT CORRESPONDING TO SEG2.THIS MAY REQUIRE A SEARCH;
WHILE TRUE DO BEGIN
IF N22-N11 <3 THEN BEGIN
OUTSTR("THIS GROUP NOT IN BOUNDARY LIST,I QUIT"&CRLF);
RETURN(0); END;
HIGH1← BDLOCBEG(SEG2); HIGH2← BDLOCEND(SEG2);
IF HIGH1= 0 OR HIGH2=0 THEN BEGIN
N22← N22 -1;
SEG2← GRPLST[BASG+N22-1];END
ELSE DONE;
END;
⊃ NOW FIND A SEGMENT IN THE MIDDLE ;
NM← N11+1;
SEGM←GRPLST[BASG+NM-1];
WHILE TRUE DO BEGIN
IF N22-NM <1 THEN BEGIN
OUTSTR("THIS GROUP NOT IN BOUNDARY LIST,I QUIT"&CRLF);
RETURN(0); END;
MID1← BDLOCBEG(SEGM); MID2← BDLOCEND(SEGM);
IF MID1= 0 OR MID2=0 THEN BEGIN
NM← NM +1;
SEGM← GRPLST[BASG+NM-1];END
ELSE DONE;
END;
IF (HIGH1-MID1)*(MID1-LOW1)<0 THEN
IF LOW1≤HIGH1 THEN LOW1← LOW1 +SIZE1
ELSE HIGH1←HIGH1 + SIZE1 ; ⊃ We have just gone around;
IF (HIGH2-MID2)*(MID2-LOW2)<0 THEN
IF LOW2≤HIGH2 THEN LOW2← LOW2 +SIZE2
ELSE HIGH2←HIGH2 + SIZE2 ; ⊃ We have just gone around;
IF LOW1> HIGH1 THEN BEGIN LOW1↔HIGH1;BSIGN1← -1; END ELSE BSIGN1← 1 ;
IF LOW2> HIGH2 THEN BEGIN LOW2↔HIGH2;BSIGN2← -1; END ELSE BSIGN2← 1 ;
BSIGN← IF BSIGN1= BSIGN2 THEN 1 ELSE -1;
⊃ PICK THE MIDDLE SECTION TO WORK ON ;
IF SIZG> 16 THEN BEGIN
N1← SIZG /4 ;
N2← 3*SIZG/4 ;END
ELSE BEGIN N1←1; N2← SIZG;END;
FITSIZ← N2 - N1 +1 ;
PAR1← 1; PAR2← FITSIZ ;
FOR I← 1 STEP 1 UNTIL FITSIZ DO BEGIN
ISEG← IF BSIGN1=1 THEN GRPLST[BASG+I+N1-2] ELSE GRPLST[BASG+N2-I];
⊃ the axis always goes from LOW1 to HIGH1 ;
U1←XBEG(ISEG);V1←YBEG(ISEG);
THREED(LINNO[ISEG],U1,V1,X1,Y1,Z1);
U2←XEND(ISEG);V2←YEND(ISEG);
THREED(LINNO[ISEG],U2,V2,X2,Y2,Z2);
AXX[I]← (X1+X2)/2;AXY[I]← (Y1+Y2)/2;AXZ[I]← (Z1+Z2)/2;
DX← ((U1+U2)/2)*DPYMULT + XOFFSET ;
DY← ((V1+V2)/2)*DPYMULT + XOFFSET ;
APOINT(DX,DY);
END;
DPYOUT(1);
⊃ NOW FIT A LINE TO THIS ;
LINE3D(PAR1,PAR2,AXX,AXY,AXZ,AXLINE);
LINEDPY(AXX,AXY,AXZ,AXLINE,PAR1,FITSIZ);
⊃ NOW ITERATE ON THIS PART ;
⊃ LET US START WORKING WITH THE BOUNDARY RIGHT FROM HERE ;
CORTHRESH ← 0.1 ;
FOR J LOOP 5 DO BEGIN "INIT_ITER"
INTEGER I2,I3,INTS1,INTS2;
DEFINE ADD1 = {
⊃ NOW ADD FIVE MORE SEGMENTS ;
IF INTS1= NEAR_HIGH THEN HIGH1← HIGH1 + 5 ELSE LOW1← LOW1-5;
IF ( (BD1= BD2) AND ( (ABS((HIGH1-LOW2) MOD SIZE1)<5) OR
(ABS((HIGH2-LOW1) MOD SIZE1)<5) ))
OR (ABS((HIGH1-LOW1) MOD SIZE1)< 5) THEN DONE; };
DEFINE ADD2 = {
⊃ NOW ADD FIVE MORE SEGMENTS ;
IF INTS2= NEAR_HIGH THEN HIGH2← HIGH2+5 ELSE LOW2← LOW2-5;
IF ( (BD1= BD2) AND ( (ABS((HIGH2-LOW1) MOD SIZE2)<5) OR
(ABS((HIGH1-LOW2) MOD SIZE2)<5) ))
OR (ABS((HIGH2-LOW2) MOD SIZE2)< 5) THEN DONE; };
COR← 0 ;M← 0 ;
DPYSET(DPYBUF);
FOR I← PAR1 STEP 1 UNTIL PAR2 DO BEGIN
NORMAL(AXLINE,PLANE,AXX[I],AXY[I],AXZ[I]);
I2← 0;
IF (INT3DBD(BD1,LOW1,HIGH1,PLANE,XN1,YN1,ZN1,INDX1)= FOUND ) AND
(INT3DBD(BD2,LOW2,HIGH2,PLANE,XN2,YN2,ZN2,INDX2)= FOUND ) THEN BEGIN
⊃ We have found intersections with both boundaries ;
M← M +1 ;
IF M= 1 THEN BEGIN PAR1← I ; LOW1F←HIGH1F←INDX1;
LOW2F←HIGH2F← INDX2; END
ELSE BEGIN
HIGH1F←INDX1;
IF BSIGN>0 THEN HIGH2F←INDX2 ELSE LOW2F←INDX2;END;
PAR← PAR1 + M-1 ;
XN← (XN1+XN2)/2; YN ←(YN1+YN2)/2 ; ZN← (ZN1+ ZN2)/2;
XCOR← XN - AXX[I];YCOR← YN - AXY[I];ZCOR← ZN -AXZ[I];
IF AXDEB THEN BEGIN
ITYPE(I);ITYPE(M);RTYPE(XCOR);RTYPE(YCOR);RTYPE(ZCOR);
OUTSTR(CRLF); END;
COR ← COR + SQRT(XCOR↑2 + YCOR↑2 + ZCOR↑2 ) ;
AXX[PAR]← XN ;AXY[PAR]← YN ; AXZ[PAR]← ZN ;
RAD[PAR]←SQRT((XN1-XN2)↑2 + (YN1-YN2)↑2 + (ZN1-ZN2)↑2) ;
⊃ ************* ANGLE[PAR]← AXLINE[6];
TWOD(XN,YN,ZN,DX,DY);
APOINT(DX*DPYMULT +XOFFSET,DY*DPYMULT + YOFFSET);
END;
END;
DPYOUT(1);
PAR2← PAR ;
IF M ≤ 1 THEN BEGIN OUTSTR("NOT ENOUGH POINTS LEFT"&CRLF);
RETURN(0); END;
LINE3D(PAR1,PAR2,AXX,AXY,AXZ,AXLINE);
LINEDPY(AXX,AXY,AXZ,AXLINE,PAR1,PAR2);
IF (COR/M) < CORTHRESH THEN DONE "INIT_ITER" ;
END "INIT_ITER" ;
⊃ PART OF AXIS ;
ARRTRAN(INITLIN,AXLINE); ⊃ SAVE THE INITIAL LINE FOUND FOR LATER USE ;
INITNO ← M ;
⊃ NOW EXTEND THIS AXIS ;
⊃ VARIABLES PAR1 AND PAR2 STAND FOR THE LOWER AND UPPER INDICES OF SEGMENTS
ALREADY INCLUDED.PAR IS THE INDEX OF THE CURRENT SEGMENT.;
APARABOLA(PAR1,PAR2,RAD,A,B,C,OFFSET);
⊃ CONVENTION FOR BOUNDARY VARIBLES .
LOW10,LOW20 etc. for storing results of initial iterations.
LOW1,HIGH1 etc. for keeping track of how far the boundaries have been extended,
to avoid crossing over .
LOW11,LOW22 etc. for picking a segment of boundary to pass along for finding
intersections with(i.e. just the current extensions).
LOW1F,LOW2F etc. for storing the `finished' results ;
STEPSIZ ← 1 ; PAR← PAR2; ⊃ Initialize for loop EXTND ;
LOW1←LOW1F;HIGH1←HIGH1F;LOW10←LOW1F+1;HIGH10←HIGH1F-1;⊃ START JUST INSIDE
OF INITIAL GROUP ;
LOW2←LOW2F;HIGH2←HIGH2F;LOW20←LOW2F+1;HIGH20←HIGH2F-1;
LOW11←HIGH11←HIGH10;
IF BSIGN>0 THEN LOW22←HIGH22←HIGH20 ELSE LOW22←HIGH22←LOW20 ;
TERMB←TERME← 0 ; ⊃ More Initialization ;
WHILE TRUE DO BEGIN "EXTND"
INTEGER SIGN,SEG,SEGF,SEGL,NSEG,NDX1,NDX2;
BOOLEAN NEWBD1,NEWBD2;
DEFINE ADDSEG1 = {
⊃ NOW ADD FIVE MORE SEGMENTS ;
IF STEPSIZ> 0 THEN BEGIN
IF BD1= BD2 THEN
BEGIN IF (LEN1← BDDIST(HIGH1,LOW2,SIZE1))≤0 THEN
BEGIN ⊃ TERME←-1; DONE "ONEDIR" ;END; END
ELSE IF (LEN1← BDDIST(HIGH1,LOW1,SIZE1))≤ 0 THEN
BEGIN ⊃ CIRC←-1 ;DONE "ONEDIR";END;
LOW11← HIGH11 ;
HIGH11← HIGH1← HIGH1 + ( 5 MIN LEN1 );
END
ELSE BEGIN
IF BD1= BD2 THEN
BEGIN IF (LEN1← BDDIST(HIGH2,LOW1,SIZE1))≤0 THEN
BEGIN ⊃ TERMB←-1;DONE "ONEDIR" ;END; END
ELSE IF (LEN1← BDDIST(HIGH1,LOW1,SIZE1))≤ 0 THEN
BEGIN ⊃ CIRC←-1 ;DONE "ONEDIR";END;
HIGH11← LOW11 ;
LOW11← LOW1← LOW1-(5 MIN LEN1);
END; };
DEFINE ADDSEG2 = {
⊃ NOW ADD FIVE MORE SEGMENTS ;
IF STEPSIZ* BSIGN >0 THEN BEGIN
IF BD1= BD2 THEN
BEGIN IF (LEN2← BDDIST(HIGH2,LOW1,SIZE2))≤0 THEN BEGIN
⊃ IF STEPSIZ>0 THEN TERME←-1 ELSE TERMB←-1 ;DONE "ONEDIR";END; END
ELSE IF (LEN2← BDDIST(HIGH2,LOW2,SIZE2))≤ 0 THEN BEGIN
⊃ CIRC ← ¬1 ;DONE "ONEDIR";END;
LOW22← HIGH22 ;
HIGH22← HIGH2← HIGH2 + ( 5 MIN LEN2);
END
ELSE BEGIN
IF BD1= BD2 THEN
BEGIN IF (LEN2← BDDIST(HIGH1,LOW2,SIZE2))≤0 THEN BEGIN
⊃ IF STEPSIZ>0 THEN TERME←-1 ELSE TERMB←-1 ;DONE "ONEDIR";END; END
ELSE IF (LEN2← BDDIST(HIGH2,LOW2,SIZE2))≤ 0 THEN BEGIN
⊃ CIRC← -1;DONE "ONEDIR";END;
HIGH22← LOW22 ;
LOW22← LOW2← LOW2 -( 5 MIN LEN2);
END; };
WHILE TRUE DO BEGIN "ONEDIR"
PAR← PAR + STEPSIZ ;
IF STEPSIZ> 0 THEN BEGIN PAR2← ILAST← PAR;
IFIRST←(PAR1 MAX (PAR-4));⊃ FIVE POINTS IF AVAILAIBLE;
END
ELSE BEGIN PAR1←IFIRST← PAR;
ILAST← (PAR2 MIN (PAR+4));END;
M← M+1 ;
FOR I LOOP 1 DO BEGIN "ONESTEP"
⊃ The loop is set up to iterate only once ;
LINEEXP(AXLINE,AXX,AXY,AXZ,RAD,PAR,STEPSIZ,XNP,YNP,ZNP);
NORMAL(AXLINE,PLANE,XNP,YNP,ZNP);
I2← 0;
WHILE INT3DBD(BD1,LOW11,HIGH11,PLANE,XN1,YN1,ZN1,INDX1)≠FOUND DO BEGIN
⊃ PUT MORE SEGMENTS IN THIS GROUP ;
IF I2=3 THEN DONE "ONEDIR"; ⊃ You have tried to extend three times;
ADDSEG1;I2← I2+1 ; END;
⊃ RESET THE INDICES ;
IF STEPSIZ>0 THEN BEGIN LOW11←HIGH11←INDX1-1;HIGH1←INDX1;END
ELSE BEGIN LOW11←HIGH11←INDX1+1;LOW1←INDX1;END;
I3←0;
WHILE INT3DBD(BD2,LOW22,HIGH22,PLANE,XN2,YN2,ZN2,INDX2)≠FOUND DO BEGIN
IF I3=3 THEN DONE "ONEDIR"; ⊃ You have tried to extend three times;
ADDSEG2;I3← I3+1 ;END;
XN←(XN1+XN2)/2;YN←(YN1+YN2)/2; ZN← (ZN1+ZN2)/2 ;
AXX[PAR]← XN ;AXY[PAR]← YN ;AXZ[PAR]← ZN ;
RAD[PAR]← SQRT((XN1-XN2)↑2 + (YN1-YN2)↑2 +(ZN1-ZN2)↑2 );
TWOD(XN,YN,ZN,DX,DY);
APOINT(DX*DPYMULT+XOFFSET,DY*DPYMULT + YOFFSET);
DPYOUT(1);
⊃ IF ABS(RADERR← (EXP_RADIUS(PAR)- RAD[PAR]))> (RADTHRESH ←0.5*C)
THEN DONE "ONEDIR" ;
⊃ *********** IF NOT RADCHECK(PAR,PAR1,PAR2,STEPSIZ,RAD) THEN DONE "ONEDIR";
⊃ *********** ANGLE[PAR]← AXLINE[6];
ERR← SQRT( (XN-XNP)↑2 + (YN-YNP)↑2 +(ZN-ZNP)↑2 );
IF STEPSIZ>0 THEN BEGIN
HIGH1F←INDX1;
IF BSIGN>0 THEN BEGIN LOW22←HIGH22←INDX2-1;HIGH2←HIGH2F←INDX2;END
ELSE BEGIN LOW22←HIGH22←INDX2+1;LOW2←LOW2F←INDX2; END;
END
ELSE BEGIN
LOW1F←INDX1;
IF BSIGN>0 THEN BEGIN LOW22←HIGH22←INDX2+1;LOW2←LOW2F←INDX2;END
ELSE BEGIN LOW22←HIGH22←INDX2-1;HIGH2←HIGH2F←INDX2; END;
END;
⊃ CHECK HERE WHETHER ANY FARTHER TO GO IN THIS DIRECTION ;
IF STEPSIZ>0 THEN BEGIN
FINDX← HIGH1F;
IF BD1=BD2 THEN BEGIN LINDX← LOW1F;SIGN←1 ;END
ELSE BEGIN LINDX←IF BSIGN>0 THEN HIGH2F ELSE LOW2F;
SIGN← BSIGN ; END; END
ELSE BEGIN
FINDX← LOW1F;
IF BD1=BD2 THEN BEGIN LINDX← HIGH1F;SIGN←-1 ;END
ELSE BEGIN LINDX←IF BSIGN>0 THEN LOW2F ELSE HIGH2F;
SIGN← - BSIGN ; END; END ;
IF TTEST THEN
IF TERMTEST(BD1,FINDX,LINDX,SIZE1,SIGN,RAD[PAR],LEN,PLANE)
THEN DONE "ONEDIR";
IF ERR < 0.2 THEN DONE "ONESTEP" ELSE BEGIN
⊃ FIT A NEW CURVE ;
LINE3D(IFIRST,ILAST,AXX,AXY,AXZ,AXLINE);
LINEDPY(AXX,AXY,AXZ,AXLINE,IFIRST,ILAST);
⊃ APARABOLA(PAR1,PAR2,RAD,A,B,C,OFFSET);
END ;
END "ONESTEP" ;
END "ONEDIR";
L1:IF STEPSIZ> 0 THEN BEGIN
⊃ RESTORE THE INDICES TO HOW FAR WE HAD ACTUALLY GONE;
HIGH1← HIGH1F;
IF BSIGN>0 THEN HIGH2← HIGH2F ELSE LOW2← LOW2F;
⊃ CHECK TO SEE IF YOU WERE NEAR THE TERMINATION;
IF BD1=BD2 THEN
BEGIN IF ISTERM(BD1,HIGH1F,LOW2F,HIGH1F,LOW2F,STEPSIZ)
THEN TERME←-1;END;
⊃ CHECK FOR CIRCULARITY WILL BE MADE ONLY WHEN GOING BACKWARDS;
STEPSIZ← -STEPSIZ;
⊃ NOW RESTORE;
ARRTRAN(AXLINE,INITLIN);PAR2← PAR2 -1;M← M-1;
LOW11←HIGH11←LOW10;
IF BSIGN>0 THEN LOW22←HIGH22←LOW20 ELSE LOW22←HIGH22←HIGH20;
PAR← PAR1; END
ELSE BEGIN
IF BD1=BD2 THEN
BEGIN IF ISTERM(BD1,LOW1F,HIGH2F,LOW1F,HIGH2F,STEPSIZ)
THEN TERMB←-1;END
ELSE BEGIN
IF BSIGN>0 THEN BEGIN NDX1←LOW2F; NDX2← HIGH2F; END
ELSE BEGIN NDX1← HIGH2F;NDX2← LOW2F; END;
IF ISTERM(BD1,LOW1F,NDX1,LOW1F,HIGH1F,STEPSIZ) AND
ISTERM(BD2,LOW1F,NDX1,NDX1,NDX2,BSIGN*STEPSIZ) THEN
CIRC← -1; END;
PAR1← PAR1 +1;M← M-1 ;
DONE "EXTND" ;END;
END "EXTND" ;
⊃ NOW DISPLAY ALL THE AXIS POINTS ;
DPYSET(DPYBUF);
FOR IPAR← PAR1 STEP 1 UNTIL PAR2 DO BEGIN
TWOD(AXX[IPAR],AXY[IPAR],AXZ[IPAR],DX,DY);
APOINT(DX*DPYMULT + XOFFSET,DY*DPYMULT + YOFFSET); END;
DPYOUT(1);
PARB←PAR1; PARE← PAR2; ⊃ BEGIN AND AND INDICES ;
RETURN(PAR2-PAR1+1); ⊃ NO OF AXIS POINTS ;
END "AXIS" ;
PROCEDURE MK_PIECES;
BEGIN "MK_PIECES"
SHORT INTEGER BD1,BD2,LOW1F,LOW2F,HIGH1F,HIGH2F,BSIGN,PARB,PARE,I,J,
BASPC,SIZPC,IGRP,NPIECES,TERMB,TERME,CIRC,SIZ1,SIZ2;
REAL ARRAY AXX,AXY,AXZ,RAD,ANGLE[-100:100];
INTEGER ARRAY FIXUP[1:128];
LABEL L1,L2;
STRING FILNAM,STR ;
SIMPLE PROCEDURE FILEOUT;
BEGIN "FILEOUT"
INTEGER DUMMY ;
LABEL AGAIN;
OPEN(DSKOUT,"DSK",8,2,2,COUNT,BRCHAR,EOF);
AGAIN:OUTSTR("OUTPUT PIECE FILE NAME");
FILNAM← INCHWL;
ENTER(DSKOUT,FILNAM,FLAG);
IF FLAG≠0 THEN BEGIN OUTSTR("ENTER UNSUCCESSFUL !!!!");GO AGAIN ;END;
WORDOUT(DSKOUT,CVSIX("PIECE1"));
WORDOUT(DSKOUT,DUMMY);
WORDOUT(DSKOUT,DUMMY);
END"FILEOUT";
SIMPLE PROCEDURE FRAMEOUT;
BEGIN "FRAMEOUT"
WORDOUT(DSKOUT,BASPC);
WORDOUT(DSKOUT,SIZPC);
ARRYOUT(DSKOUT,AXX[PARB],SIZPC);
ARRYOUT(DSKOUT,AXY[PARB],SIZPC);
ARRYOUT(DSKOUT,RAD[PARB],SIZPC);
ARRYOUT(DSKOUT,ANGLE[PARB],SIZPC);
WORDOUT(DSKOUT,BD1);
WORDOUT(DSKOUT,BD2);
WORDOUT(DSKOUT,LOW1F);
WORDOUT(DSKOUT,HIGH1F);
WORDOUT(DSKOUT,LOW2F);
WORDOUT(DSKOUT,HIGH2F);
WORDOUT(DSKOUT,1);
WORDOUT(DSKOUT,BSIGN);
WORDOUT(DSKOUT,TERMB);
WORDOUT(DSKOUT,TERME);
WORDOUT(DSKOUT,CIRC);
END "FRAMEOUT";
⊃ START OF MK_PIECES ;
FILEOUT;
BASPC← 1 ; NPIECES← 0 ;
WHILE TRUE DO BEGIN
OUTSTR("TYPE GROUP TO BE TESTED");
IGRP←INTSCAN(STR←INCHWL,J);
IF IGRP<1 OR IGRP>20 THEN DONE ;
IF GRPSIZ[IGRP]≠0 THEN
IF (SIZPC← AXIS(IGRP,BD1,BD2,LOW1F,LOW2F,HIGH1F,HIGH2F,BSIGN,PARB,PARE,
TERMB,TERME,CIRC,AXX,AXY,AXZ,RAD,ANGLE))>0 THEN BEGIN
OUTSTR("PUT THIS PIECE ON THE OUTPUT FILE?"&CRLF);
IF INCHWL = "Y" THEN BEGIN
IF BSIGN<0 THEN LOW2F↔HIGH2F; ⊃ To confirm to TOB Covention;
⊃ NOW DO SOME CLEANING UP,NORMALISE INDICES TO LIE BET 1 AND SIZBD;
SIZ1← SIZBD[BD1] ; SIZ2← SIZBD[BD2] ;
IF LOW1F≤0 THEN LOW1F← LOW1F + SIZ1
ELSE IF LOW1F>SIZ1 THEN LOW1F← LOW1F-SIZ1;
IF HIGH1F<0 THEN HIGH1F← HIGH1F + SIZ1
ELSE IF HIGH1F>SIZ1 THEN HIGH1F← HIGH1F-SIZ1;
IF LOW2F≤0 THEN LOW2F← LOW2F + SIZ2
ELSE IF LOW2F>SIZ2 THEN LOW2F← LOW2F-SIZ2;
IF HIGH2F<0 THEN HIGH2F← HIGH2F + SIZ2
ELSE IF HIGH2F>SIZ2 THEN HIGH2F← HIGH2F-SIZ2;
FRAMEOUT ;
BASPC← BASPC + SIZPC;
NPIECES← NPIECES +1 ;
END;
END;
END;
CLOSE(DSKOUT);
⊃ NOW REWRITE THE HEADER ;
L1:LOOKUP(DSKOUT,FILNAM,FLAG);
IF FLAG≠0 THEN BEGIN OUTSTR("LOOKUP OF "&FILNAM&" UNSUCCESSFUL !!!"); INCHWL; END;
L2:ENTER(DSKOUT,FILNAM,FLAG);
IF FLAG≠0 THEN BEGIN OUTSTR("ENTER OF "&FILNAM&" UNSUCCESSFUL !!!"); INCHWL; END;
USETI(DSKOUT,1);
ARRYIN(DSKOUT,FIXUP[1],'200);
FIXUP[2]← NPIECES;
FIXUP[3]← BASPC-1 ; ⊃ NOTE THIS IS THE NO OF AXPTS UPTO THE POINT ;
USETO(DSKOUT,1);
ARRYOUT(DSKOUT,FIXUP[1],'200);
CLOSE(DSKOUT);
END "MK_PIECES" ;
⊃ MAIN PART OF OUTER ;
READPTS(FILNAM);
READCAL;
BDIN;
READPTS;
GET_BREAKS;
GROUP;
AGAIN: DISPLAY;
⊃ MOMENTS ;
MK_PIECES ; GO AGAIN;
END "OUTER" ;
COMMENT MAIN PROGRAM ;
SHORT INTEGER LINE,I,J,NUM ;
STRING STR;
LABEL AGAIN ;
IF ¬RUNBEFORE THEN BEGIN
POG ← 1 ;
DISFINAL←DPYTST=0; COMMENT DISPLAY IF THIS IS A III;
DISPOINTS←0;
DISCONT←0;
FMODE←1;
LMODE←1;
DPYMULT←3;
WHICHFRAME←1;
FIRSTONLY←FALSE;
PDEN←0.1;
SEGTHRESH← 12.0;
RUNBEFORE←-1;END
ELSE DPYCLR;
STEPANGLE← PI/432 ;
FINALPOG←GETPOG;
IF DISPOINTS THEN PPOG←GETPOG;
⊃ OPEN(OUTFILE,"DSK",0,0,3,0,0,EOF);
⊃ IF EOF ≠ 0 THEN OUTSTR("OPEN FAILS ON OUTFILE"&CRLF);
⊃ ENTER(OUTFILE,"DEBUG1",FLAG);
⊃ IF FLAG THEN OUTSTR("ENTER FAILS ON OUTFILE"&CRLF);
OUTSTR("DEBUG?");
IF INCHWL = "Y" THEN DEBUG←TRUE ELSE DEBUG ← FALSE ;
OUTSTR("DISPLAY?");
IF INCHWL = "Y" THEN DISP←TRUE ELSE DISP ← FALSE;
FILIN(FILNAM);
OUTER;
END "CURVE" ;