perm filename CURVE.SAI[S,AIL] blob
sn#044816 filedate 1973-05-29 generic text, type T, neo UTF8
00100
00200
00300 ⊃ Here we declare some arrays that are used by several procedures
00400 and some parameters which are set in INITIAL;
00500
00600
00700
00800 SHORT INTEGER ARRAY HIST[-1000:1000];
00900 SHORT REAL ARRAY PEAKS[0:30];
01000 DEFINE LINEPTS(I) = { EDGES[EDGE_LIMIT-(I)] };
01100 SHORT REAL THFACTOR,CFACTOR,NDC2,NDTH2,NDS3,NDS32,NDOK;
01200 SHORT INTEGER NMIN,NPEAKS,NELLIPS,NVERT,MINELLS,RMAX;
01300 SHORT INTEGER MAXRES,MINRES,NHIGH_1,FLAG;
01400 SHORT INTEGER MAXT,MAXC,MAXR,MAXDR,RHIGH,BIT_FACTOR;
01500 SHORT INTEGER INDEX, EI, L_INDEX, F_INDEX, LINE_INDEX;
01600 SHORT INTEGER MINVAL;
01700 SHORT REAL NDRAD,NDRSQ;
01800
01900
02000
02100 DEFINE PTH =EDGE_TH -EBS+1, PC=EDGE_C -EBS+1,
02200 POPX=EDGE_X-EBS+1, POPY=EDGE_Y-EBS+1;
02300
02400 DEFINE PI32={(1.5*PI1)};
02500
02600 DEFINE FPOINT="20",THET="1",CEE="2",SIGNX="11",SIGNY="12",
02700 THET1="1",THET2="2",XELLIPSE="9",YELLIPSE="10",RELLIPSE="11",
02800 THRANGE="12",CHIVALUE="7",VERT_END1="3",VERT_END2="4",
02900 THETXY=FPOINT-3,
03000 VERT_END=FPOINT-2, DIRSIN=FPOINT-1, DIRCOS=FPOINT;
03100
03200
03300 SIMPLE INTEGER PROCEDURE SIGN(REAL X);
03400 RETURN(X/ABS(X));
03500
03600
03700
03800
03900 SIMPLE REAL PROCEDURE GET_ANGLE(SHORT REAL CX,CY);
04000
04100 ⊃ Here we return the angle defined by delta-x and delta-y;
04200 ⊃ Phi's range from -π/2 to 3π/2;
04300
04400 BEGIN
04500 SHORT REAL PHI;
04600 IF ABS(CX) ≥.10@-8 THEN
04700 BEGIN
04800 PHI ← ATAN(CY/CX);
04900 IF CX<0 THEN PHI←PHI+PI1;
05000 END ELSE
05100 IF CY≥0 THEN PHI ← PIO2 ELSE PHI ← -PIO2;
05200 RETURN(PHI);
05300 END;
05400
05500
00100
00200
00300 SIMPLE PROCEDURE ELLIP_DISP(INTEGER N);
00400
00500 ⊃ Here we display ellipse "N" which is stored in array ELLIPS.
00600 Also "A(N)" or "E(N)" is displayed;
00700
00800 BEGIN "EDISP"
00900 SHORT REAL X0,Y0,R0,THETA1,THETA2,TXX,TYY;
01000 LABEL NUM_ARC;
01100 X0←ELLIPS[N,XELLIPSE];
01200 Y0←ELLIPS[N,YELLIPSE];
01300 R0←ELLIPS[N,RELLIPSE];
01400 THETA1←ELLIPS[N,THET1];
01500 THETA2←ELLIPS[N,THET2];
01600 TXX←TX(X0); TYY←TY(Y0);
01700 IF THETA1=THETA2 THEN
01800 MKCIRCLE(R0*SCAL,TXX,TYY,25) ELSE
01900 MKARC(THETA1+PI1,THETA2+PI1,R0*SCAL,TXX,TYY,25);
02000 NUM_ARC: AIVECT(TXX,TYY);
02100 IF CAL_COMP THEN DPYBIG(6);
02200 IF THETA1=THETA2 THEN DPYSST("E"&CVS(N)) ELSE
02300 DPYSST("A"&CVS(N));
02400 END "EDISP";
02500
00100
00200
00300
00400 BOOLEAN PROCEDURE CURVETEST(SHORT INTEGER N; SHORT REAL X0,Y0,NDTH,RADIUS;
00500 SHORT REAL ARRAY LTEMP,ELANG);
00600
00700 ⊃ Here we determine if there are a set of edge points in array
00800 LTEMP that form an arc with center near (X0,Y0) and THETA
00900 value near the angle of the tangent;
01000
01100 BEGIN "CURVETEST"
01200 SHORT INTEGER I,J,CYC,MEAN,ISAVE,K,MIN,SUMN, SUM_PLUS,SUM_MINUS;
01300 SHORT INTEGER TEST3,IMAX,IMIN,MAX,AI,AK;
01400 SHORT INTEGER MAX2,ISAVE2,LOWEST,DEL,LOWT;
01500 BOOLEAN TEST, GAP1_FOUND;
01600 STRING PTEST;
01700 SHORT REAL TEMP,SUMXCUBE,SUMYCUBE,SUMXYSQ,SUMYXSQ;
01800 SHORT REAL RESID,THETA,MAXANG,MAXSAVE,ANGAP;
01900 SHORT REAL MG_LENGTH, GAP, DEN, RESMIN, RESMAX;
02000 SHORT REAL CHI,RES,RADSQ,NRES,DX,DY,NDTHEL;
02100 SHORT REAL SUMX,SUMY,SUMXY,SUMXSQ,SUMYSQ,SUM_RESID;
02200 SHORT REAL XFIT,YFIT,RFIT,DELX,DELY,DELXSQ,DELYSQ;
02300 SHORT REAL XOLD,YOLD,FF,GG,DFDX,DFDY,DGDX,DGDY,DELTA_X,DELTA_Y;
02400 LABEL AGAIN,AFTER,AGAIN2,AGAIN3,AGAIN4,AGAIN5,LAST1,LAST2,POOR_ELLIPSE;
02500 LABEL AVER,TLINE,PART1,PART2,PART3,FIRSTK,LASTK, BEFORE;
02600 LABEL LAS1,LAS2,BEF1,NEWT,FILL,CONT_FILL,RID;
02700
02800
02900
03000
03100 SIMPLE PROCEDURE ZERO_SUM;
03200 BEGIN
03300 SUMXCUBE←0; SUMYCUBE←0; SUMX←SUMY←SUMXYSQ←0;
03400 SUMXSQ←SUMYSQ←SUMYXSQ←SUMXY←0;
03500 END;
03600
03700
03800
03900 SIMPLE PROCEDURE SUM_STORE(SHORT INTEGER I);
04000 BEGIN
04100 DELX←LTEMP[I+POPX]-XFIT;
04200 DELY←LTEMP[I+POPY]-YFIT;
04300 DELXSQ←DELX*DELX;
04400 DELYSQ←DELY*DELY;
04500 SUMX←SUMX+DELX;
04600 SUMY←SUMY+DELY;
04700 SUMXY←SUMXY+DELX*DELY;
04800 SUMYXSQ←SUMYXSQ+DELY*DELXSQ;
04900 SUMXYSQ←SUMXYSQ+DELX*DELYSQ;
05000 SUMXSQ←SUMXSQ+DELXSQ;
05100 SUMYSQ←SUMYSQ+DELYSQ;
05200 SUMXCUBE←SUMXCUBE+DELX*DELXSQ;
05300 SUMYCUBE←SUMYCUBE+DELY*DELYSQ;
05400 END;
05500
05600
05700
05800 SIMPLE PROCEDURE RES_ORDER;
05900 BEGIN
06000 SHORT INTEGER I,J,LASTPT;
06100 LASTPT←INDEX;
06200 SUMN←0; ZERO_SUM;
06300 FOR I←EBS STEP EBS UNTIL LASTPT DO
06400 BEGIN "INDE"
06500 DX←LTEMP[I+POPX]-XFIT;
06600 DY←LTEMP[I+POPY]-YFIT;
06700 NRES←(DX↑2+DY↑2-RADSQ)/NDRESC;
06800 IF NRES<RESMIN∨NRES>RESMAX THEN
06900 BEGIN
07000 FOR J←-3 STEP 1 UNTIL 0 DO
07100 LTEMP[I+J]↔LTEMP[LASTPT+J];
07200 I←I-EBS;
07300 LASTPT←LASTPT-EBS;
07400 END ELSE
07500 BEGIN
07600 SUMN←SUMN+1;
07700 SUM_STORE(I);
07800 END;
07900 END "INDE";
08000 F_INDEX←LASTPT;
08100 END;
08200
08300
08400
08500 SIMPLE PROCEDURE THET_ORDER;
08600 BEGIN
08700 SHORT INTEGER I,J,ISAVE,LOWT,LOWEST;
08800 SHORT REAL MIN;
08900 LABEL AGAIN3;
09000
09100 LOWEST←1;
09200 AGAIN3: MIN←1000;
09300 FOR I←LOWEST STEP 1 UNTIL L_INDEX DO
09400 IF ELANG[I]<MIN THEN
09500 BEGIN ISAVE←I; MIN←ELANG[I]; END;
09600 IF ISAVE≠LOWEST THEN
09700 BEGIN
09800 ELANG[ISAVE]↔ELANG[LOWEST];
09900 ISAVE←ISAVE*EBS;
10000 LOWT←LOWEST*EBS;
10100 FOR J←POPX STEP 1 UNTIL POPX+EBS-1 DO
10200 LTEMP[ISAVE+J]↔LTEMP[LOWT+J];
10300 END;
10400 LOWEST←LOWEST+1;
10500 IF LOWEST≠L_INDEX THEN GO TO AGAIN3;
10600 END;
10700
10800
10900
11000
11100 SIMPLE PROCEDURE GAP_CHECK;
11200 BEGIN
11300 SHORT INTEGER I;
11400 FOR I←2 STEP 1 UNTIL L_INDEX DO
11500 IF ELANG[I]-ELANG[I-1]>MAXANG THEN
11600 BEGIN ISAVE←I; MAXANG←ELANG[I]-ELANG[I-1]; END;
11700 END;
11800
11900
12000
12100 SIMPLE INTEGER PROCEDURE TEST_END_GAP;
12200 BEGIN
12300 SHORT INTEGER I,J,K,N,MAX;
12400 LABEL FIRST,BACK_CHECK,LAST;
12500 ⊃ OUTSTR(CRLF&" IN TEST_END_GAP ");
12600 ANGAP←NDS/RFIT;
12700 N←0;
12800 F_INDEX←INDEX;
12900 FOR J←2 STEP 1 UNTIL L_INDEX DO
13000 IF ANGAP>ELANG[J]-ELANG[J-1] THEN
13100 GO TO FIRST ELSE N←N+1;
13200 FIRST: IF N=0 THEN GO TO BACK_CHECK;
13300 MAX←(J-1)*EBS;
13400 FOR I←EBS STEP EBS UNTIL MAX DO
13500 BEGIN
13600 ⊃ NOW WE PUT THE BAD-FIT EDGE POINTS BACK IN ARRAY EDGES;
13700 EDGE_INDEX ← EDGE_INDEX +1;
13800 EI ← EDGE_INDEX*EDGE_BLSIZE;
13900 EDGES[EI+PTH] ← LTEMP[I+PTH];
14000 EDGES[EI+PC] ← LTEMP[I+PC];
14100 EDGES[EI+POPX] ← LTEMP[I+POPX];
14200 EDGES[EI+POPY] ← LTEMP[I+POPY];
14300 L_INDEX ← L_INDEX -1;
14400 INDEX ← INDEX -EDGE_BLSIZE;
14500 END;
14600 K←0;
14700 FOR I←MAX+EBS STEP EBS UNTIL F_INDEX DO
14800 BEGIN
14900 K←K+EBS;
15000 AI←I DIV EBS; AK←K DIV EBS;
15100 LTEMP[K+PTH]←LTEMP[I+PTH];
15200 LTEMP[K+PC]←LTEMP[I+PC];
15300 LTEMP[K+POPX]←LTEMP[I+POPX];
15400 LTEMP[K+POPY]←LTEMP[I+POPY];
15500 ELANG[AK]←ELANG[AI];
15600 END;
15700 BACK_CHECK: IF L_INDEX≤NHIGH THEN RETURN(-1);
15800 FOR I←L_INDEX STEP -1 UNTIL 2 DO
15900 IF ANGAP>ELANG[I]- ELANG[I-1] THEN GO TO LAST ELSE
16000 BEGIN
16100 ⊃ NOW WE PUT THE BAD-FIT EDGE POINTS BACK IN ARRAY EDGES;
16200 EDGE_INDEX ← EDGE_INDEX +1;
16300 EI ← EDGE_INDEX*EDGE_BLSIZE;
16400 EDGES[EI+PTH] ← LTEMP[I*EBS+PTH];
16500 EDGES[EI+PC] ← LTEMP[I*EBS+PC];
16600 EDGES[EI+POPX] ← LTEMP[I*EBS+POPX];
16700 EDGES[EI+POPY] ← LTEMP[I*EBS+POPY];
16800 L_INDEX ← L_INDEX -1;
16900 INDEX ← INDEX -EDGE_BLSIZE;
17000 N←N+1;
17100 END;
17200 PRINT(N);
17300 LAST: IF L_INDEX≤NHIGH THEN RETURN(-1);
17400 IF N>2 THEN RETURN(+1) ELSE RETURN(0);
17500 END;
17600
17700
17800
17900
18000
18100
18200 RESID←0; NRES←100;
18300 XFIT←X0;
18400 YFIT←Y0;
18500 RADSQ←RADIUS↑2;
18600 INDEX ← INDEX -EDGE_BLSIZE;
18700 L_INDEX ← INDEX DIV EDGE_BLSIZE;
18800
18900 BEF1: CYC←-1; TEST3←0; NDTHEL←NDTH; DEL←1;
19000 BEFORE:
19100 SUM_PLUS←SUM_MINUS←0;
19200 FOR J←-995 STEP 1 UNTIL 995 DO HIST[J]←0;
19300 ⊃ HERE WE START 2 CYCLES TO SEPARATE EDGE_POINTS BY THEIR:
19400 THETA-VALUES AND RESIDUALS;
19500 AGAIN: CYC←CYC+1; ZERO_SUM; SUMN←0;
19600 MINRES←1000; MAXRES←-1000;
19700 AGAIN2: FOR I←EDGE_BLSIZE STEP EDGE_BLSIZE UNTIL L_INDEX*EDGE_BLSIZE DO
19800 BEGIN "KGET"
19900 DX←LTEMP[I+POPX]-XFIT;
20000 DY←LTEMP[I+POPY]-YFIT;
20100 THETA←GET_ANGLE(DX,DY);
20200 TEMP←ABS(LTEMP[I+PTH]-THETA);
20300 WHILE TEMP>PIO2 DO TEMP←TEMP-PI1;
20400 IF ABS(TEMP)<NDTHEL THEN
20500 BEGIN
20600 RES←DX↑2+DY↑2-RADSQ;
20700 NRES←RES/NDRESC;
20800 ⊃ OUTSTR(" RES="&CVG(RES)&CRLF);
20900 IF CYC=0 THEN
21000 BEGIN
21100 IF NRES>-990.∧NRES<990. THEN
21200 BEGIN
21300 J←NRES;
21400 IF J<MINRES THEN MINRES←J
21500 ELSE IF J>MAXRES THEN MAXRES←J;
21600 HIST[J]←HIST[J]+1;
21700 END
21800 ELSE
21900 BEGIN
22000 IF NRES<-990. THEN SUM_MINUS←SUM_MINUS+1
22100 ELSE SUM_PLUS←SUM_PLUS+1;
22200 END;
22300 GO TO LASTK;
22400 END;
22500 IF NRES<RESMIN∨NRES>RESMAX THEN GO TO AFTER;
22600 ⊃ OUTSTR(CRLF&"RESIDUAL="&CVG(RESID)&" SUMN="&CVS(SUMN));
22700 FIRSTK: SUMN←SUMN+1;
22800 SUM_STORE(I);
22900 ⊃ OUTSTR(CRLF&"RESIDUAL="&CVG(RESID)&
23000 " X="&CVG(LTEMP[I+POPX])&
23100 " Y="&CVG(LTEMP[I+POPY])&
23200 " TH="&CVG(LTEMP[I+PTH])&
23300 " C="&CVG(LTEMP[I+PC]));
23400
23500 END ELSE
23600 AFTER: IF CYC>0 THEN
23700 BEGIN
23800 ⊃ NOW WE PUT THE BAD-FIT EDGE POINTS BACK IN ARRAY EDGES;
23900 EDGE_INDEX ← EDGE_INDEX +1;
24000 EI ← EDGE_INDEX*EDGE_BLSIZE;
24100 EDGES[EI+PTH] ← LTEMP[I+PTH];
24200 EDGES[EI+PC] ← LTEMP[I+PC];
24300 EDGES[EI+POPX] ← LTEMP[I+POPX];
24400 EDGES[EI+POPY] ← LTEMP[I+POPY];
24500 LTEMP[I+PTH]←LTEMP[INDEX+PTH];
24600 LTEMP[I+PC]←LTEMP[INDEX+PC];
24700 LTEMP[I+POPX]←LTEMP[INDEX+POPX];
24800 LTEMP[I+POPY]←LTEMP[INDEX+POPY];
24900 I ← I - EDGE_BLSIZE;
25000 L_INDEX ← L_INDEX -1;
25100 INDEX ← INDEX -EDGE_BLSIZE;
25200 END;
25300 LASTK: END "KGET";
25400 F_INDEX←INDEX;
25500 IF CYC=0 THEN
25600 BEGIN
25700 MAX←0;
25800 PRINT(MINRES); PRINT(MAXRES);
25900 PRINT(SUM_MINUS); PRINT(SUM_PLUS);
26000 FOR J←MINRES STEP DEL UNTIL MAXRES DO
26100 IF HIST[J]≥MAX THEN
26200 BEGIN
26300 IF HIST[J]=MAX∧ABS(J)>ABS(ISAVE) THEN CONTINUE;
26400 MAX←HIST[J]; ISAVE←J;
26500 END;
26600 H_DISPLAY(MINRES-2,DEL,MAXRES+2,MAX,HIST);
26700 IF MAX<RHIGH THEN GO TO POOR_ELLIPSE;
26800 IF ¬TEST3∧(ISAVE<-MINVAL∨ISAVE>MINVAL) THEN
26900 TEST3←-1 ELSE TEST3←0;
27000
27100 PART1: FOR J←ISAVE STEP DEL UNTIL MAXRES DO
27200 IF HIST[J]>0 THEN IMAX←J
27300 ELSE IF HIST[J+DEL]=0∧HIST[J+2*DEL]=0∧HIST[J+3*DEL]=0 THEN
27400 GO TO PART2;
27500
27600 PART2: FOR J←ISAVE STEP -DEL UNTIL MINRES DO
27700 IF HIST[J]>0 THEN IMIN←J
27800 ELSE IF HIST[J-DEL]=0∧HIST[J-2*DEL]=0∧HIST[J-3*DEL]=0 THEN
27900 GO TO PART3;
28000
28100 PART3: RESMIN←IMIN-1;
28200 RESMAX←IMAX+1;
28300 PRINT(RESMIN); PRINT(RESMAX);
28400 H_DISPLAY(RESMIN-2,DEL,RESMAX+2,MAX,HIST);
28500 DPYCLR;
28600 DPYSET(BUF);
28700 INIT_DOMAIN(X1,Y2,X2,Y1);
28800 BOUNDARY(X1,Y2,X2,Y1);
28900 AIVECT(-300,420);
29000 DPYSST("CYCLE=0 BEFORE OMITTING BAD POINTS");
29100 EDGE_DISP(0,L_INDEX,X1,Y1,X2,Y2,LTEMP);
29200 DPYOUT(1);
29300 INCHWL;
29400 INCHWL;
29500 DPYCLR;
29600 IF CYC≠0 THEN RETURN(0);
29700 SUMN←0;
29800 FOR J←IMIN STEP DEL UNTIL IMAX DO
29900 SUMN← SUMN+HIST[J];
30000 OUTSTR(CRLF); PRINT(SUMN);
30100 IF SUMN≤NHIGH THEN GO TO POOR_ELLIPSE;
30200 IF ¬TEST3 THEN GO TO AGAIN ELSE RES_ORDER;
30300
30400 END;
30500
30600
30700
30800 ⊃ THIS IS USED WHEN THERE ARE TOO FEW GOOD EDGE POINTS;
30900 IF RESMAX<-9000. THEN RETURN(0);
31000 AVER: IF SUMN≤NHIGH∧¬TEST3 THEN GO TO POOR_ELLIPSE;
31100
31200 NEWT:
31300 RADSQ←(SUMXSQ+SUMYSQ)/SUMN;
31400 DFDX←2.*(SUMX*SUMX/SUMN-SUMXSQ);
31500 DGDY←2.*(SUMY*SUMY/SUMN-SUMYSQ);
31600 FF←SUMXCUBE+SUMXYSQ-SUMX*RADSQ;
31700 GG←SUMYCUBE+SUMYXSQ-SUMY*RADSQ;
31800 DFDY←DGDX←2.*(SUMX*SUMY/SUMN-SUMXY);
31900 DEN←DGDY*DFDX-DGDX*DFDY;
32000 DELTA_Y←(FF*DGDX-GG*DFDX)/DEN;
32100 DELTA_X←(GG*DFDY-FF*DGDY)/DEN;
32200 XOLD←XFIT; YOLD←YFIT;
32300 XFIT←XOLD+DELTA_X;
32400 YFIT←YOLD+DELTA_Y;
32500 RFIT←SQRT(RADSQ);
32600 OUTSTR(CRLF); PRINT(RADSQ); PRINT(XOLD); PRINT(YOLD);
32700 OUTSTR(CRLF); PRINT(RFIT); PRINT(XFIT); PRINT(YFIT);
32800 OUTSTR(CRLF); PRINT(FF); PRINT(GG);
32900 IF ABS(XFIT-XOLD)+ABS(YFIT-YOLD)>NDOK THEN
33000 BEGIN
33100 ZERO_SUM;
33200 FOR I←EBS STEP EBS UNTIL F_INDEX DO
33300 SUM_STORE(I);
33400 GO TO NEWT;
33500 END;
33600
33700 IF TEST3 THEN BEGIN CYC←-1; GO TO BEFORE; END;
33800
33900 RFIT←SQRT(RADSQ);
34000 OUTSTR(CRLF); PRINT(RFIT);
34100 DPYSET(BUF);
34200 AIVECT(-300,420);
34300 DPYSST("SHOWING FIT");
34400 INIT_DOMAIN(X1,Y2,X2,Y1);
34500 BOUNDARY(X1,Y2,X2,Y1);
34600 EDGE_DISP(0,L_INDEX,X1,Y1,X2,Y2,LTEMP);
34700 DPYOUT(1);
34800 INCHWL;
34900 MKCIRCLE(RFIT*SCAL,TX(XFIT),TY(YFIT),25);
35000 DPYOUT(1);
35100 INCHWL; INCHWL;
35200 DPYCLR;
35300 LAST1: SUM_RESID←0;
35400 FOR I←EBS STEP EBS UNTIL INDEX DO
35500 BEGIN
35600 K←I DIV EBS;
35700 DELX←LTEMP[I+POPX]-XFIT;
35800 DELY←LTEMP[I+POPY]-YFIT;
35900 ELANG[K]←GET_ANGLE(DELX,DELY);
36000 DELXSQ←DELX*DELX;
36100 DELYSQ←DELY*DELY;
36200 RESID←(DELXSQ+DELYSQ-RADSQ)↑2;
36300 SUM_RESID←SUM_RESID+RESID;
36400 END;
36500 CHI←SQRT(SUM_RESID/(SUMN-3));
36600 OUTSTR(CRLF&"X0="&CVG(X0)
36700 &" XFIT="&CVG(XFIT)
36800 &" Y0="&CVG(Y0)
36900 &" YFIT="&CVG(YFIT)&"
37000 R0="&CVG(RADIUS)
37100 &" RFIT="&CVG(RFIT)
37200 &" CHI="&CVG(CHI));
37300 IF CHI>NDCHI THEN GO TO POOR_ELLIPSE ELSE GO TO TLINE;
37400 POOR_ELLIPSE: GO TO RID;
37500 ⊃ SINCE THINGS ARE NOT WORKING OUT, WE TRY SPLITTING THE EDGE_POINTS
37600 IN ORDER TO TAKE CARE OF CASE WHEN 2 ELLIPSES GIVE US A POOR CENTER;
37700 IF CYC=0∧(SUM_MINUS>NHIGH∨SUM_PLUS>NHIGH) THEN
37800 BEGIN
37900 IF DISP_POINTS THEN
38000 BEGIN
38100 DPYSET(BUF);
38200 AIVECT(-300,420);
38300 DPYSST("BEFORE SPLIT");
38400 INIT_DOMAIN(X1,Y2,X2,Y1);
38500 BOUNDARY(X1,Y2,X2,Y1);
38600 MKCIRCLE(RFIT*SCAL,TX(XFIT),TY(YFIT),25);
38700 EDGE_DISP(0,L_INDEX,X1,Y1,X2,Y2,LTEMP);
38800 DPYOUT(1);
38900 INCHWL;
39000 INCHWL;
39100 DPYCLR;
39200 END;
39300 CYC←-1;
39400 ⊃ SPLIT(TH_AVE,C_AVE,SINTH,COSTH,LTEMP,SUMTH,SUMC,SUMN);
39500 OUTSTR(CRLF&" IN SPLIT ");
39600 GO TO BEFORE;
39700 END
39800 ELSE
39900 BEGIN
40000 IF DISP_POINTS THEN BEGIN "QUESTION" OUTSTR(CRLF&"THIS IS A POOR FIT --
40100 ONLY "&CVS(SUMN)&" GOOD POINTS. CHI = "&CVG(CHI)&"
40200 DO YOU WANT TO SEE THE EDGES -- Y OR N?");
40300 PTEST ← INCHWL;
40400 IF EQU(PTEST,"Y") THEN
40500 BEGIN
40600 DPYSET(BUF);
40700 INIT_DOMAIN(X1,Y2,X2,Y1);
40800 BOUNDARY(X1,Y2,X2,Y1);
40900 EDGE_DISP(0,L_INDEX,X1,Y1,X2,Y2,LTEMP);
41000 DPYOUT(1);
41100 INCHWL;
41200 INCHWL;
41300 DPYCLR;
41400 END; END "QUESTION";
41500 ⊃ WE ARE GOING TO GET RID OF THIS ELLIPSE;
41600 RID: RESMAX←-10000; CYC←2;
41700 ELLIPS[N,2]←ELLIPS[NELLIPS-1,2];
41800 ELLIPS[N,1]←ELLIPS[NELLIPS-1,1];
41900 NELLIPS←NELLIPS-1;
42000 GO TO AGAIN2;
42100 END;
42200 TLINE: ⊃ IF ¬DISP_POINTS THEN GO TO LAST2;
42300 DPYCLR;
42400 DPYSET(BUF);
42500 AIVECT(-300,420);
42600 DPYSST("CHI="&CVG(CHI)&" NO. OF POINTS="&CVG(SUMN));
42700 BOUNDARY(X1,Y2,X2,Y1);
42800 EDGE_DISP(0,L_INDEX,X1,Y1,X2,Y2,LTEMP);
42900 DPYOUT(1);
43000 INCHWL;
43100 INCHWL;
43200 LAST2: GAP1_FOUND←0;
43300 THET_ORDER;
43400 IF ¬DISP_POINTS THEN GO TO AGAIN4;
43500 DPYSET(BUF);
43600 BOUNDARY(X1,Y2,X2,Y1);
43700 FOR I←EBS STEP EBS UNTIL INDEX DO
43800 BEGIN
43900 OPX←LTEMP[I+POPX];
44000 OPY←LTEMP[I+POPY];
44100 TEMP←LTEMP[I+PTH];
44200 CX←COS(TEMP);
44300 CY←SIN(TEMP);
44400 DISP_EDGE;
44500 DPYOUT(1);
44600 END;
44700 DPYCLR;
44800 ⊃ FOR I←1 STEP 1 UNTIL L_INDEX DO
44900 OUTSTR(CRLF&" I="&CVS(I)&" X="&CVG(LTEMP[I*EBS+POPX])
45000 &" Y="&CVG(LTEMP[I*EBS+POPY])
45100 &" THETA="&CVG(ELANG[I]));
45200 AGAIN4: MEAN←L_INDEX DIV 2;
45300 MG_LENGTH←NDGAP*(ELANG[L_INDEX]-ELANG[1]);
45400 TEMP←NDS3/RFIT;
45500 IF MG_LENGTH>TEMP THEN MG_LENGTH←TEMP;
45600 OUTSTR(CRLF&"INDEX="&CVS(INDEX)
45700 &" MEAN="&CVS(MEAN)
45800 &" MG_LENGTH="&CVG(MG_LENGTH));
45900
46000 MAXSAVE←MAXANG←PIT2+ELANG[1]-ELANG[L_INDEX];
46100 ISAVE←1;
46200 GAP_CHECK;
46300 IF MAXANG<MG_LENGTH THEN
46400 BEGIN
46500 ELLIPS[N,THET1]←ELLIPS[N,THET2]←ELANG[1];
46600 ELLIPS[N,THRANGE]←PIT2;
46700 GO TO CONT_FILL; ⊃ We have found a circle;
46800 END;
46900 AGAIN3: IF MAXANG≠MAXSAVE THEN
47000 BEGIN
47100 ⊃ Order points starting at end of largest gap;
47200 FOR I←1 STEP 1 UNTIL ISAVE-1 DO
47300 ELANG[I]←ELANG[I]+PIT2;
47400 THET_ORDER;
47500 END;
47600 BEGIN
47700 TEMP←TEST_END_GAP;
47800 IF TEMP=1 THEN BEGIN CYC←-1; GO TO BEFORE; END;
47900 IF TEMP=-1 THEN GO TO POOR_ELLIPSE;
48000 MAXANG←0;
48100 GAP_CHECK;
48200 IF MAXANG<MG_LENGTH THEN GO TO FILL;
48300 END;
48400 AGAIN5: F_INDEX←INDEX;
48500 FOR J←2 STEP 1 UNTIL L_INDEX DO
48600 BEGIN
48700 GAP←ELANG[J]-ELANG[J-1];
48800 IF GAP>MG_LENGTH THEN
48900 BEGIN "ELIMINATE"
49000 GAP1_FOUND←-1;
49100 MAX←(J-1)*EBS;
49200 IF DISP_POINTS THEN OUTSTR(CRLF&"GAP FOUND N="&CVS(N)
49300 &" GAP="&CVG(GAP)&" MG_LENGTH="&CVG(MG_LENGTH));
49400 IF J<MEAN THEN
49500 BEGIN "FIRST"
49600 FOR I←EBS STEP EBS UNTIL MAX DO
49700 BEGIN
49800 ⊃ NOW WE PUT THE BAD-FIT EDGE POINTS BACK IN ARRAY EDGES;
49900 EDGE_INDEX ← EDGE_INDEX +1;
50000 EI ← EDGE_INDEX*EDGE_BLSIZE;
50100 EDGES[EI+PTH] ← LTEMP[I+PTH];
50200 EDGES[EI+PC] ← LTEMP[I+PC];
50300 EDGES[EI+POPX] ← LTEMP[I+POPX];
50400 EDGES[EI+POPY] ← LTEMP[I+POPY];
50500 L_INDEX ← L_INDEX -1;
50600 INDEX ← INDEX -EDGE_BLSIZE;
50700 END;
50800 K←0;
50900 FOR I←MAX+EBS STEP EBS UNTIL F_INDEX DO
51000 BEGIN
51100 K←K+EBS;
51200 AI←I DIV EBS; AK←K DIV EBS;
51300 LTEMP[K+PTH]←LTEMP[I+PTH];
51400 LTEMP[K+PC]←LTEMP[I+PC];
51500 LTEMP[K+POPX]←LTEMP[I+POPX];
51600 LTEMP[K+POPY]←LTEMP[I+POPY];
51700 ELANG[AK]←ELANG[AI];
51800 END;
51900 END "FIRST"
52000 ELSE
52100 FOR I←MAX+EBS STEP EBS UNTIL F_INDEX DO
52200 BEGIN
52300 ⊃ NOW WE PUT THE BAD-FIT EDGE POINTS BACK IN ARRAY EDGES;
52400 EDGE_INDEX ← EDGE_INDEX +1;
52500 EI ← EDGE_INDEX*EDGE_BLSIZE;
52600 EDGES[EI+PTH] ← LTEMP[I+PTH];
52700 EDGES[EI+PC] ← LTEMP[I+PC];
52800 EDGES[EI+POPX] ← LTEMP[I+POPX];
52900 EDGES[EI+POPY] ← LTEMP[I+POPY];
53000 L_INDEX ← L_INDEX -1;
53100 INDEX ← INDEX -EDGE_BLSIZE;
53200 END;
53300 IF L_INDEX ≤ NHIGH THEN GO TO POOR_ELLIPSE;
53400 MG_LENGTH←NDGAP*(ELANG[L_INDEX]-ELANG[1]);
53500 TEMP←NDS32/RFIT;
53600 IF MG_LENGTH<TEMP THEN MG_LENGTH←TEMP;
53700 GO TO AGAIN5;
53800 END "ELIMINATE";
53900 END;
54000 IF GAP1_FOUND THEN BEGIN CYC←-1; GO TO BEFORE; END;
54100 OUTSTR(CRLF&" NO GAPS FOUND THIS CYCLE");
54200 ⊃ Filling array ellips;
54300 FILL: ELLIPS[N,THET1]←ELANG[1];
54400 ELLIPS[N,THET2]←ELANG[L_INDEX];
54500 ELLIPS[N,THRANGE]←ELANG[L_INDEX]-ELANG[1];
54600 CONT_FILL: ELLIPS[N,XELLIPSE]←XFIT;
54700 ELLIPS[N,YELLIPSE]←YFIT;
54800 ELLIPS[N,RELLIPSE]←RFIT;
54900 ELLIPS[N,5]←SUMN;
55000 ELLIPS[N,CHIVALUE]←CHI;
55100 ELLIPS[N,VERT_END1]←1;
55200 ELLIPS[N,VERT_END2]←L_INDEX;
55300 DPYCLR;
55400 DPYSET(BUF);
55500 AIVECT(-300,420);
55600 DPYSST("CHI="&CVG(CHI)&" NO. OF POINTS="&CVG(SUMN));
55700 BOUNDARY(X1,Y2,X2,Y1);
55800 EDGE_DISP(0,L_INDEX,X1,Y1,X2,Y2,LTEMP);
55900 DPYOUT(1);
56000 INCHWL;
56100 INCHWL;
56200 RETURN(-1);
56300 END "CURVETEST";
56400
00100
00200
00300
00400
00500 PROCEDURE CONNELL(SHORT INTEGER N;SHORT REAL ARRAY ELANG);
00600 BEGIN "CONNELL"
00700 SHORT INTEGER I,J,INDEX1,INDEX2, ISAVE,JHUMP;
00800 SHORT INTEGER TEMP,LOWEST;
00900 SHORT REAL THETA1,THETA2,X0,Y0,R0,ANGAP32,MIN,DEN,GAP;
01000 BOOLEAN SOLID;
01100 LABEL OMIT;
01200 OUTSTR(CRLF&" IN CONNELL ");
01300 X0←ELLIPS[N,XELLIPSE];
01400 Y0←ELLIPS[N,YELLIPSE];
01500 R0←ELLIPS[N,RELLIPSE];
01600 THETA1←ELLIPS[N,THET1];
01700 THETA2←ELLIPS[N,THET2];
01800 INDEX1←ELLIPS[N,VERT_END1];
01900 ⊃ STARTING LOCATION OF EDGES IN ARRAY ELANG;
02000 INDEX2←ELLIPS[N,VERT_END2];
02100 ⊃ ENDING LOCATION OF EDGES IN ARRAY ELANG;
02200 ELLIPS[N,VERT_END1]←-1;
02300 ELLIPS[N,VERT_END2]←-1;
02400 ELLIPS[N,13]←-1;
02500 ELLIPS[N,14]←-1;
02600 ELLIPS[N,15]←-1;
02700 ELLIPS[N,16]←-1;
02800 ANGAP32←NDS32/R0;
02900 SOLID←-1; JHUMP←FPOINT+3;
03000
03100 ELLIPS[N,FPOINT+1]←X0+R0*COS(THETA1);
03200 ELLIPS[N,FPOINT+2]←Y0+R0*SIN(THETA1);
03300 ISAVE←1;
03400 FOR I←2 STEP 1 UNTIL INDEX2 DO
03500 BEGIN
03600 GAP←ELANG[I]-ELANG[I-1];
03700 IF (SOLID∧GAP<ANGAP32)∨(¬SOLID∧GAP>ANGAP32) THEN CONTINUE;
03800 ELLIPS[N,JHUMP]← SOLID;
03900 ELLIPS[N,JHUMP+1]← X0+ R0*COS(ELANG[I]);
04000 ELLIPS[N,JHUMP+2]← Y0+ R0*SIN(ELANG[I]);
04100 JHUMP←JHUMP+3;
04200 SOLID←¬SOLID;
04300 ISAVE←ISAVE+1;
04400 END;
04500 ⊃ TO END THE LINE WE DO THE FOLLOWING;
04600 ELLIPS[N,JHUMP]← SOLID;
04700 ELLIPS[N,JHUMP+1]← X0+R0*COS(ELANG[INDEX2]);
04800 ELLIPS[N,JHUMP+2]← Y0+R0*SIN(ELANG[INDEX2]);
04900 ISAVE←ISAVE+1;
05000 ⊃ FILLING ARRAY HUMPS;
05100 ELLIPS[N,5]←ISAVE; ⊃ THIS IS THE NUMBER OF POINTS
05200 ALONG THE LINE;
05300 ELLIPS[N,6]←FPOINT+3*(ISAVE-1); ⊃ THIS IS LOCATION JUST
05400 BEFORE X-LOCATION OF LAST POINT;
05500 ⊃ FOR I←INDEX1 STEP 1 UNTIL INDEX2 DO
05600 OUTSTR(CRLF&"THETA="&CVG(ELANG[I])
05700 &" X="&CVG(SIN(ELANG[I]))
05800 &" Y="&CVG(COS(ELANG[I])));
05900 JHUMP←JHUMP+3;
06000 ⊃ FOR I←1 STEP 1 UNTIL JHUMP DO
06100 OUTSTR(CRLF&" ELLIPS -- I="&CVS(I)&" "&CVG(ELLIPS[N,I]));
06200 IF ¬DISP_POINTS THEN GO TO OMIT;
06300 DPYSET(BUF);
06400 BOUNDARY(X1,Y2,X2,Y1);
06500 ELLIP_DISP(N);
06600 DPYOUT(1);
06700 INCHWL;
06800 INCHWL;
06900 OMIT:
07000 END "CONNELL";
07100
07200
07300
00100
00200
00300 BOOLEAN PROCEDURE LINER(SHORT REAL X0,Y0,RADIUS);
00400
00500 ⊃ HERE THE EDGE-POINTS ARE SEPARATED INTO A GROUP FOR EACH LINE;
00600
00700 BEGIN "LINER"
00800 SHORT INTEGER N,I;
00900 SHORT REAL TH_AVE,C_AVE,TEMP;
01000 LABEL BEF1,LAST;
01100 SHORT REAL ARRAY LTEMP[1:LT_LIMIT],ELANG[1:LT_LIMIT DIV 4];
01200
01300 INDEX ← EDGE_BLSIZE;
01400 ⊃ HERE WE TAKE ALL EDGE_POINTS AND STORE THEM IN ARRAY LTEMP;
01500 FOR I←EDGE_BLSIZE STEP EDGE_BLSIZE UNTIL EDGE_INDEX*EDGE_BLSIZE DO
01600 BEGIN "LGET"
01700 BEGIN
01800 LTEMP[INDEX+PTH] ← EDGES[I+PTH];
01900 LTEMP[INDEX+PC] ← EDGES[I+PC];
02000 LTEMP[INDEX+POPX] ← EDGES[I+POPX];
02100 LTEMP[INDEX+POPY] ← EDGES[I+POPY];
02200 EI ← EDGE_INDEX*EDGE_BLSIZE;
02300 EDGES[I+POPX] ← EDGES[EI+POPX];
02400 EDGES[I+POPY] ← EDGES[EI+POPY];
02500 EDGES[I+PC] ← EDGES[EI+PC];
02600 EDGES[I+PTH] ← EDGES[EI+PTH];
02700 EDGE_INDEX ← EDGE_INDEX -1;
02800 INDEX ← INDEX + EDGE_BLSIZE;
02900 I ← I - EDGE_BLSIZE;
03000 END;
03100 IF INDEX≥LT_LIMIT THEN GO TO BEF1;
03200 END "LGET";
03300 BEF1:
03400 L_INDEX←INDEX DIV EBS;
03500 MINVAL←80;
03600 NDS3←3.*NDS;
03700 NDS32←1.5*NDS;
03800 RHIGH←2;
03900 NHIGH←5;
04000 NDOK←.05;
04100 NELLIPS←1;
04200 IF CURVETEST(1,X0,Y0,NDTH,RADIUS,LTEMP,ELANG) THEN
04300 BEGIN CONNELL(1,ELANG); NELLS←NELLS+1; END;
04400 IF ¬DISP_POINTS THEN GO TO LAST;
04500 DPYSET(BUF);
04600 LINE_FIND(EDGE_INDEX,EDGES);
04700 POINTER(TH_AVE,NDTH,C_AVE,NDC);
04800 LINE_FIND(L_INDEX,LTEMP);
04900 AIVECT(-300,420);
05000 DPYSST("TH_AVE="&CVG(TH_AVE)&" C_AVE="&CVG(C_AVE));
05100 DPYOUT(1);
05200 INCHWL;
05300 INCHWL;
05400 DPYCLR;
05500 DPYSET(BUF);
05600 AIVECT(-300,420);
05700 DPYSST("CHI="&CVG(ELLIPS[N,7])&" NO. OF POINTS="&CVS(ELLIPS[N,8]));
05800 INIT_DOMAIN(X1,Y2,X2,Y1);
05900 BOUNDARY(X1,Y2,X2,Y1);
06000 EDGE_DISP(0,L_INDEX,X1,Y1,X2,Y2,LTEMP);
06100 DPYOUT(1);
06200 INCHWL;
06300 LEDGE_INDEX ← LINE_INDEX DIV EDGE_BLSIZE;
06400 INIT_DOMAIN(X1,Y2,X2,Y1);
06500 BOUNDARY(X1,Y2,X2,Y1);
06600 EDGE_DISP(-1,LEDGE_INDEX,X1,Y1,X2,Y2,EDGES);
06700 ⊃ THIS WILL RESULT IN A DISPLAY OF THE LINE EDGE_POINTS;
06800 DPYOUT(1);
06900 INCHWL;
07000 INCHWL;
07100 DPYCLR;
07200 LAST:
07300 RETURN(-1);
07400 END "LINER";
07500
07600