C ALGORITHM 475 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN COMM. ACM, VOL. 17, NO. 03, C P. 152. PROGRAM ACMTEST DEM00001 C DEMONSTRATION PROGRAM DEM00002 DIMENSION EYE(3),S(4),ST1(80,80,2),IS2(3,160) DEM00003 DIMENSION IOBJ(80,80) DEM00004 C USE WHOLE FRAME DEM00005 S(1)=0. DEM00006 S(2)=1. DEM00007 S(3)=0. DEM00008 S(4)=1. DEM00009 C SET EYE POSITION DEM00010 EYE(1)=250. DEM00011 EYE(2)=150. DEM00012 EYE(3)=100. DEM00013 C INITIALIZE PACKAGE DEM00014 CALL INIT3D(EYE,80,80,80,ST1,3,160,IS2,9,S) DEM00015 C CREATE AND PLOT TEST OBJECT DEM00016 DO 4 I=1,80 DEM00017 A=(I-50)**2 DEM00018 DO 3 J=1,80 DEM00019 C=(J-25)**2 DEM00020 D=IABS(J-63)+IABS(I-25) DEM00021 DO 3 K=1,80 DEM00022 C FLOOR DEM00023 IF(K.EQ.1) GO TO 1 DEM00024 C BALL DEM00025 IF(SQRT(A+C+(FLOAT(K)-25.)**2).LE.25.) GO TO 1 DEM00026 C POINT DEM00027 IF(D.GT.FLOAT(80-K)*.1875) GO TO 2 DEM00028 1 IOBJ(J,K)=1 DEM00029 GO TO 3 DEM00030 2 IOBJ(J,K)=0 DEM00031 3 CONTINUE DEM00032 4 CALL DANDR(80,80,ST1,3,160,160,IS2,9,S,IOBJ,80) DEM00033 C ADVANCE TO THE NEXT FRAME. DEM00034 CALL FRAME DEM00035 C A SECOND PICTURE WILL NOW BE CALLED USING THE SAME SIZE DEM00036 C ARRAYS AND EYE POSITION. THIS MEANS THE CALL TO INIT3D, DEM00037 C THE BIGGEST TIME CONSUMER, CAN BE SKIPPED IF THE FOLLOWING DEM00038 C FOUR LINES ARE INCLUDED. DEM00039 REWIND 9 DEM00040 DO 5 I=1,3 DEM00041 DO 5 J=1,160 DEM00042 5 IS2(I,J)=0 DEM00043 C THIS PICTURE WILL BE THE T=4 CONTOUR SURFACE OF DEM00044 C T=1/SQRT(U*U+V*V+W*W)+(.5-V)**2/SQRT(U*U+V*V). DEM00045 DO 9 I=1,80 DEM00046 U=(40.5-FLOAT(I))/79. DEM00047 UU=U*U DEM00048 DO 8 J=1,80 DEM00049 V=(FLOAT(J)-40.5)/79. DEM00050 VV=V*V DEM00051 A=1./SQRT(UU+VV) DEM00052 DO 8 K=1,80 DEM00053 C THE FOLLOWING CARD ADDS AXES. DEM00054 IF(I*J.EQ.1.OR.I*K.EQ.1.OR.J*K.EQ.1) GO TO 6 DEM00055 W=(FLOAT(K)-40.5)/79. DEM00056 IF(1./SQRT(UU+VV+W*W)+(.5-V)**2*A.LE.4.) GO TO 7 DEM00057 6 IOBJ(J,K)=1 DEM00058 GO TO 8 DEM00059 7 IOBJ(J,K)=0 DEM00060 8 CONTINUE DEM00061 9 CALL DANDR(80,80,ST1,3,160,160,IS2,9,S,IOBJ,80) DEM00062 C FLUSH PLOT BUFFER DEM00063 CALL FRAME DEM00064 STOP DEM00065 END DEM00066 SUBROUTINE INIT3D(EYE,NU,NV,NW,ST1,LX,NY,IS2,IU,S) 3D 00067 DIMENSION EYE(3),ST1(NV,NW,2),IS2(LX,NY),S(4) C C BY THOMAS WRIGHT C COMPUTING FACILITY C THE NATIONAL CENTER FOR ATMOSPHERIC RESEARCH C BOULDER, COLORADO 80302 C NCAR IS SPONSORED BY THE NATIONAL SCIENCE FOUNDATION. C C THE METHOD IS DESCRIBED IN DETAIL IN - A ONE-PASS HIDDEN- C LINE REMOVER FOR COMPUTER DRAWN THREE-SPACE OBJECTS. PROC C 1972 SUMMER COMPUTER SIMULATION CONFERENCE, 261-267, 1972. C C THIS VERSION IS FOR USE ON CDC 6000 OR 7000 COMPUTERS. C C THIS PACKAGE OF ROUTINES PLOTS 3-DIMENSIONAL OBJECTS WITH C HIDDEN PARTS NOT SHOWN. OBJECTS ARE STORED IN AN ARRAY, C WITH THE POSITION IN THE ARRAY CORRESPONDING TO A LOCATION C IN 3-SPACE AND THE VALUE OF THE ARRAY ELEMENT TELLING IF C ANY OBJECT IS PRESENT AT THE LOCATION. C C INIT3D IS AN INITIALIZATION ROUTINE FOR THIS PACKAGE. IT C IS CALLED, THEN A SEQUENCE OF CALLS ARE MADE TO DANDR TO C PRODUCE A PICTURE. C EYE AN ARRAY 3 LONG CONTAINING THE U, V, AND W COORDI- C NATES OF THE EYE POSITION. OBJECTS ARE CONSIDERED C TO BE IN A BOX WITH 2 EXTREME CORNERS AT (1,1,1) AND C (NU,NV,NW). THE EYE POSITION MUST HAVE POSITIVE C COORDINATES AWAY FROM THE COORDINATE PLANES U=0, C V=0, AND W=0. WHILE GAINING EXPERIENCE WITH THE C PACKAGE, USE EYE(1)=5*NU, EYE(2)=4*NV, EYE(3)=3*NW. C NU U DIRECTION LENGTH OF THE BOX CONTAINING THE OBJECTS C NV V DIRECTION LENGTH OF THE BOX CONTAINING THE OBJECTS C NW W DIRECTION LENGTH OF THE BOX CONTAINING THE OBJECTS C ST1 A SCRATCH ARRAY AT LEAST NV*NW*2 WORDS LONG. C LX FIRST DIMENSION OF A SCRATCH ARRAY, IS2, USED BY THE C PACKAGE FOR REMEMBERING WHERE IT SHOULD NOT DRAW. C LX=1+NX/NBPW. SEE DANDR COMMENTS FOR NX AND NBPW. C NY SECOND DIMENSION OF IS2. SEE DANDR COMMENTS. C IS2 A SCRATCH ARRAY AT LEAST LX*NY WORDS LONG. C IU UNIT NUMBER OF SCRATCH FILE FOR THE PACKAGE. ST1 C WILL BE WRITTEN NU TIMES ON THIS FILE. C S AN ARRAY 4 LONG WHICH CONTAINS THE COORDINATES OF C THE AREA WHERE THE PICTURE IS TO BE DRAWN. THAT IS, C ALL PLOTTING COORDINATES GENERATED WILL BE BOUNDED C AS FOLLOWS-- X COORDINATES WILL BE BETWEEN S(1) AND C S(2), Y COORDINATES WILL BE BETWEEN S(3) AND S(4). C TO PREVENT DISTORTION, HAVE S(2)-S(1)=S(4)-S(3). C C IF SEVERAL PICTURES ARE TO BE DRAWN WITH THE SAME SIZE C ARRAYS AND EYE POSITION AND THE USER REWINDS IU AND FILLS C IS2 WITH ZEROES, INIT3D NEED NOT BE CALLED FOR OTHER THAN C THE FIRST PICTURE. C C SET UP TRANSFORMATION ROUTINE FOR THIS LINE OF SIGHT. U=NU V=NV W=NW CALL SETORG(U*.5,V*.5,W*.5,EYE(1),EYE(2),EYE(3)) C FIND EXTREMES IN TRANSFORMED SPACE. CALL PERSPC(1.,1.,W,D,YT,D) CALL PERSPC(U,V,1.,D,YB,D) CALL PERSPC(U,1.,1.,XL,D,D) CALL PERSPC(1.,V,1.,XR,D,D) C ADJUST EXTREMES TO PREVENT DISTORTION WHEN GOING FROM C TRANSFORMED SPACE TO PLOTTER SPACE. DIF=(XR-XL-YT+YB)*.5 IF(DIF) 1,3,2 1 XL=XL+DIF XR=XR-DIF GO TO 3 2 YB=YB-DIF YT=YT+DIF 3 REWIND IU C FIND THE PLOTTER COORDINATES OF THE 3-SPACE LATTICE POINTS C1=.9*(S(2)-S(1))/(XR-XL) C2=.05*(S(2)-S(1))+S(1) C3=.9*(S(4)-S(3))/(YT-YB) C4=.05*(S(4)-S(3))+S(3) DO 5 I=1,NU U=NU+1-I DO 4 J=1,NV V=J DO 4 K=1,NW CALL PERSPC(U,V,FLOAT(K),X,Y,D) ST1(J,K,1)=C1*(X-XL)+C2 4 ST1(J,K,2)=C3*(Y-YB)+C4 C WRITE THEM ON UNIT IU. 5 WRITE(IU) ST1 REWIND IU C ZERO OUT ARRAY WHERE VISIBILITY IS REMEMBERED. DO 6 J=1,NY DO 6 I=1,LX 6 IS2(I,J)=0 RETURN END SUBROUTINE SETORG(X,Y,Z,XT,YT,ZT) 3D 00163 C C THIS ROUTINE IMPLEMENTS THE 3-SPACE TO 2-SPACE TRANSFOR- C MATION BY KUBER, SZABO AND GIULIERI, THE PERSPECTIVE C REPRESENTATION OF FUNCTIONS OF TWO VARIABLES. J. ACM 15, C 2, 193-204,1968. C SETORG ARGUMENTS C X,Y,Z ARE THE 3-SPACE COORDINATES OF THE INTERSECTION C OF THE LINE OF SIGHT AND THE IMAGE PLANE. THIS C POINT CAN BE THOUGHT OF AS THE POINT LOOKED AT. C XT,YT,ZT ARE THE 3-SPACE COORDINATES OF THE EYE POSITION. C C PERSPC ARGUMENTS C X,Y,Z ARE THE 3-SPACE COORDINATES OF A POINT TO BE C TRANSFORMED. C XT,YT THE RESULTS OF THE 3-SPACE TO 2-SPACE TRANSFOR- C MATION. C ZT NOT USED. C C STORE THE PARAMETERS OF THE SETORG CALL FOR USE WHEN C PERSPC IS CALLED. AX=X AY=Y AZ=Z EX=XT EY=YT EZ=ZT C AS MUCH COMPUTATION AS POSSIBLE IS DONE DURING EXECUTION C OF SETORG SINCE PERSPC IS CALLED THOUSANDS OF TIMES FOR C EACH CALL TO SETORG. DX=AX-EX DY=AY-EY DZ=AZ-EZ D=SQRT(DX*DX+DY*DY+DZ*DZ) COSAL=DX/D COSBE=DY/D COSGA=DZ/D AL=ACOS(COSAL) BE=ACOS(COSBE) GA=ACOS(COSGA) SINGA=SIN(GA) C THE 3-SPACE POINT LOOKED AT IS TRANSFORMED INTO (0,0) OF C THE 2-SPACE. THE 3-SPACE Z AXIS IS TRANSFORMED INTO THE C 2-SPACE Y AXIS. IF THE LINE OF SIGHT IS CLOSE TO PARALLEL C TO THE 3-SPACE Z AXIS, THE 3-SPACE Y AXIS IS CHOSEN (IN- C STEAD OF THE 3-SPACE Z AXIS) TO BE TRANSFORMED INTO THE C 2-SPACE Y AXIS. IF(SINGA.LT.0.0001) GO TO 1 R=1./SINGA ASSIGN 2 TO JUMP RETURN 1 SINBE=SIN(BE) R=1./SINBE ASSIGN 3 TO JUMP RETURN C******************** ENTRY PERSPC *********************** ENTRY PERSPC Q=D/((X-EX)*COSAL+(Y-EY)*COSBE+(Z-EZ)*COSGA) GO TO JUMP,(2,3) 2 XT=((EX+Q*(X-EX)-AX)*COSBE-(EY+Q*(Y-EY)-AY)*COSAL)*R YT=(EZ+Q*(Z-EZ)-AZ)*R RETURN 3 XT=((EZ+Q*(Z-EZ)-AZ)*COSAL-(EX+Q*(X-EX)-AX)*COSGA)*R YT=(EY+Q*(Y-EY)-AY)*R RETURN END SUBROUTINE DANDR(NV,NW,ST1,LX,NX,NY,IS2,IU,S,IOBJS,MV) 3D 00229 DIMENSION ST1(NV,NW,2),IS2(LX,NY),S(4),IOBJS(MV,NW) C C THIS ROUTINE IS CALLED NU TIMES, EACH CALL PROCESSING THE C PART OF THE PICTURE AT U=NU+1-I WHERE I IS THE NUMBER OF C THE CALL TO DANDR. THAT IS, THE PART OF THE PICTURE AT C U=NU IS PROCESSED DURING THE FIRST CALL, THE PART OF THE C PICTURE AT U=NU-1 IS PROCESSED DURING THE SECOND CALL, AND C SO ON UNTIL THE PART OF THE PICTURE AT U=1 IS PROCESSED C DURING THE LAST CALL. C NV SEE INIT3D COMMENTS. C NW SEE INIT3D COMMENTS. C ST1 SEE INIT3D COMMENTS. C LX THE NUMBER OF WORDS NEEDED TO HOLD NX BITS. ALSO, C THE FIRST DIMENSION OF IS2. C NX NUMBER OF CELLS IN THE X DIRECTION OF A MODEL OF THE C IMAGE PLANE. A SILHOUETTE OF THE PARTS OF THE PIC- C TURE PROCESSED SO FAR IS STORED IN THIS MODEL. LINES C TO BE DRAWN ARE TESTED FOR VISIBILITY BY EXAMINING C THE SILHOUETTE. LINES IN THE SILHOUETTE ARE HIDDEN. C LINES OUT OF THE SILHOUETTE ARE VISIBLE. THE SOLU- C TION IS APPROXIMATE BECAUSE THE SILHOUETTE IS NOT C FORMED EXACTLY. SEE IS2 COMMENT BELOW. C NY NUMBER OF CELLS IN THE Y DIRECTION OF THE MODEL OF C THE IMAGE PLANE. ALSO THE SECOND DIMENSION OF IS2. C IS2 AN ARRAY TO HOLD THE IMAGE PLANE MODEL. IT IS C DIMENSIONED LX BY NY. THE MODEL IS NX BY NY AND C PACKED DENSELY. IF HIDDEN LINES ARE DRAWN, DECREASE C NX AND NY (AND LX IF POSSIBLE). IF VISIBLE LINES C ARE LEFT OUT OF THE PICTURE, INCREASE NX AND NY (AND C LX IF NEED BE). AS A GUIDE, SOME EXAMPLES SHOWING C SUCCESSFUL CHOICES ARE LISTED C GIVEN NU NV NW RESULTING NX NY FROM TESTING C 100 100 60 200 200 C 60 60 60 110 110 C 40 40 40 75 75 C IU SEE INIT3D COMMENTS. C IOBJS A NV BY NW ARRAY (WITH ACTUAL FIRST DIMENSION MV IN C THE CALLING PROGRAM) DESCRIBING THE OBJECT. IF THIS C IS CALL NUMBER I TO DANDR, THE PART OF THE PICTURE C AT U=NU+1-I IS TO BE PROCESSED. IOBJS DEFINES THE C OBJECTS TO BE DRAWN IN THE FOLLOWING MANNER -- C IOBJS(J,K)=1 IF ANY OBJECT CONTAINS THE POINT C (NU+1-I,J,K) AND IOBJS(J,K)=0 OTHERWISE. C MV ACTUAL FIRST DIMENSION OF IOBJS IN THE CALLING C PROGRAM. C C************** MACHINE DEPENDANT CONSTANTS **************** C NBPW NUMBER OF BITS PER WORD C MASK AN ARRAY NBPW LONG. MASK(I)=2**(I-1), I=1,2,...,NBPW C CDC 6000 OR 7000 VERSION DIMENSION MASK(60) DATA NBPW/60/ DATA MASK/1B,2B,4B,10B,20B,40B,100B,200B,400B,1000B, A 2000B,4000B,10000B,20000B,40000B,100000B,200000B, B 400000B,1000000B,2000000B,4000000B,10000000B, C 20000000B,40000000B,100000000B,200000000B,400000000B, D 1000000000B,2000000000B,4000000000B,10000000000B, E 20000000000B,40000000000B,100000000000B, F 200000000000B,400000000000B,1000000000000B, G 2000000000000B,4000000000000B,10000000000000B, H 20000000000000B,40000000000000B,100000000000000B, I 200000000000000B,400000000000000B,1000000000000000B, J 2000000000000000B,4000000000000000B, K 10000000000000000B,20000000000000000B, L 40000000000000000B,100000000000000000B, M 200000000000000000B,400000000000000000B, N 1000000000000000000B,2000000000000000000B, O 4000000000000000000B,10000000000000000000B, P 20000000000000000000B,40000000000000000000B/ C ASSIGN 12 TO IRET C RX AND RY ARE USED TO MAP PLOTTER COORDINATES INTO THE C IMAGE PLANE MODEL. RX=(FLOAT(NX)-1.)/(S(2)-S(1)) RY=(FLOAT(NY)-1.)/(S(4)-S(3)) C READ THE RELATIVE PLOTTER COORDINATES OF THE LATTICE C POINTS FROM UNIT IU. READ(IU) ST1 C DX, DY AND DZ ARE USED TO FIND REQUIRED COORDINATES OF C NON-LATTICE POINTS. NVD2=NV/2 NWD2=NW/2 DX=(ST1(NV,NWD2,1)-ST1(1,NWD2,1))*.5/(FLOAT(NV)-1.) DY=(ST1(1,NWD2,2)-ST1(NV,NWD2,2))*.5/(FLOAT(NV)-1.) DZ=(ST1(NVD2,NW,2)-ST1(NVD2,1,2))*.5/(FLOAT(NW)-1.) C SLOPE IS USED TO DEFORM THE IMAGE PLANE MODEL SO THAT C LINES OF CONSTANT Y OF THE IMAGE MODEL HAVE THE SAME C SLOPE AS LINES OF CONSTANT U AND W IN THE PICTURE. THIS C IMPROVES THE PICTURE. SLOPE=DY/DX C THE FOLLOWING LOOPS THROUGH STATEMENT 12 GENERATE THE .5 C CONTOUR LINES IN 2-SPACE FOR THE ARRAY IOBJS (WHICH CON- C TAINS ONLY ZEROES AND ONES), TESTS THE LINES FOR VISIBIL- C ITY, AND CALLS A ROUTINE TO PLOT THE VISIBLE LINES. DO 12 I=2,NV JUMP=IOBJS(I-1,1)*8+IOBJS(I,1)*4+1 DO 12 J=2,NW X=ST1(I,J,1) Y=ST1(I,J,2) C DECIDE WHICH OF THE 16 POSSIBILITIES THIS IS. JUMP=(JUMP-1)/4+IOBJS(I-1,J)*8+IOBJS(I,J)*4+1 GO TO (12,2,4,5,7,8,3,10,10,1,8,7,5,4,2,12),JUMP C GOING TO 1 MEANS JUMP=10 WHICH MEANS ONLY THE LOWER-RIGHT C AND UPPER-LEFT ELEMENTS OF THIS CELL ARE SET TO 1. C TWO LINES SHOULD BE DRAWN, A DIAGONAL CONNECTING THE C MIDDLE OF THE BOTTOM TO THE MIDDLE OF THE RIGHT SIDE OF C THE CELL (LOWER-RIGHT LINE), AND A DIAGONAL CONNECTING THE C MIDDLE OF THE LEFT SIDE TO THE MIDDLE OF THE TOP (UPPER- C LEFT LINE) OF THE CELL. 1 ASSIGN 9 TO IRET C LOWER-RIGHT LINE 2 X1=X Y1=Y-DZ X2=X+DX Y2=Y-DY GO TO 11 C LOWER-LEFT AND UPPER-RIGHT 3 ASSIGN 6 TO IRET C LOWER-LEFT 4 X1=X Y1=Y-DZ X2=X-DX Y2=Y+DY GO TO 11 C HORIZONTAL 5 X1=X+DX Y1=Y-DY X2=X-DX Y2=Y+DY GO TO 11 C UPPER-LEFT 6 ASSIGN 12 TO IRET 7 X1=X+DX Y1=Y-DY X2=X Y2=Y+DZ GO TO 11 C VERTICAL 8 X1=X Y1=Y-DZ X2=X Y2=Y+DZ GO TO 11 9 ASSIGN 12 TO IRET C UPPER-LEFT 10 X1=X-DX Y1=Y+DY X2=X Y2=Y+DZ C TEST VISIBILITY OF THIS LINE SEGMENT. 11 IX=(X1-S(1))*RX IY=MOD(IFIX((Y1-S(3))*RY-SLOPE*FLOAT(IX))+NY,NY)+1 IBIT=MOD(IX,NBPW)+1 IX=IX/NBPW+1 C*********** .AND. USED AS A MASKING OPERATOR ************** IV=IS2(IX,IY).AND.MASK(IBIT) C IF EITHER END OF THE LINE IS AT A MARKED SPOT ON THE IMAGE C PLANE MODEL, THE LINE IS HIDDEN IF(IV.NE.0) GO TO IRET,(6,9,12) IX=(X2-S(1))*RX IY=MOD(IFIX((Y2-S(3))*RY-SLOPE*FLOAT(IX))+NY,NY)+1 IBIT=MOD(IX,NBPW)+1 IX=IX/NBPW+1 C*********** .AND. USED AS A MASKING OPERATOR ************** IV=IS2(IX,IY).AND.MASK(IBIT) IF(IV.NE.0) GO TO IRET,(6,9,12) C*************** UNDEFINED EXTERNAL REFERENCE ************** C SUBROUTINE LINE(X1,Y1,X2,Y2) IS ASSUMED TO DRAW A LINE C FROM (X1,Y1) TO (X2,Y2) CALL LINE(X1,Y1,X2,Y2) GO TO IRET,(6,9,12) 12 CONTINUE C CODE THROUGH STATEMENT 13 CREATES AN APPROXIMATION OF C THE SILHOUETTE OF THE PART OF THE PICTURE JUST DRAWN BY C MARKING THE IMAGE PLANE MODEL WHERE THE OBJECT OCCURS. DO 13 I=1,NV DO 13 J=1,NW IF(IOBJS(I,J).EQ.0) GO TO 13 IX=(ST1(I,J,1)-S(1))*RX+0.5 TWK=SLOPE*FLOAT(IX)-0.5 IY=MOD(IFIX((ST1(I,J,2)-S(3))*RY-TWK)+NY,NY)+1 IBIT=MOD(IX,NBPW)+1 IX=IX/NBPW+1 C************ .OR. USED AS A MASKING OPERATOR ************** IS2(IX,IY)=IS2(IX,IY).OR.MASK(IBIT) 13 CONTINUE RETURN END