Almost everybody must have seen the rotating molecule demonstration produced by Acorn. If not, it was included in two parts on the RISC User discs for Volume 1 Issues 5 and 6. While very impressive, this demo only plotted a set of frames previously created using another program. This used a technique called Ray Tracing to work out the appearance of the image at any instant, and it is this basic program that is given here. In fact, the program is fairly general, and the version here produces a series of frames which can subsequently be animated to produce a rotating letter 'A' made up of coloured spheres. One frame from the sequence is, actually, one of the sample sprites supplied with RISC OS.
The program from listing 1 should be typed in and saved before running it. The sprites generated as the program runs are particularly memory hungry, and you might need to perform some reconfiguring to increase application workspace before the program will run. The program will not run on a 305.
When the program is run, you will be prompted for a start frame, which should be a number between 0 and 129 (usually 0). The program will then proceed to calculate each frame, drawing it as it goes. Before each frame is drawn, a rough view of its shape is plotted in grey. This is the image projected forward onto the two dimensional screen. At the top of the screen you will see the frame number, the time taken on this frame, and the estimated time to complete the frame. The time to completion will initially rise, and then steady off. This is because it is calculated from the amount so far done, and the time taken. However, each pixel in the frame does not take exactly the same time to be traced. Once frame 129 has been calculated, the entire sprite file (about 550K long) is saved with the name 'RaySprites'.
The entire 130 frames take several hours to draw. If you press Escape, the program saves the current position in a file called 'StopData', and also the current sprites in the file 'RaySprites'. If when the program is re-run, you specify a negative start frame, the stop position is reloaded and the process continues.
Once the program has run, you will have a 130 frame sequence on disc. This is in the form of a single sprite file containing sprites named 000, 001, 002 ... 129. To animate the sequence, the program from listing 2 can be used. Next month we will feature a program that displays the animation sequence in a multi-tasking window, automatically rescaling it to the size of the window.
The explanation which follows concentrates on the principles involved, rather than the detailed operation of this long and complex program. If you want to find out more, an excellent book on this, and all aspects of computer graphics, is: Procedural Elements for Computer Graphics by David F. Rogers published by McGraw Hill in the International student series.
Ray tracing is very much a brute-force method of rendering an image (rendering is the process of making a conceptual object appear solid). Ray tracing works by considering the plight of every ray emitted from a light source, and this is done by taking each ray, and reflecting it off the object by following the laws of physics and optics. When a ray hits the screen, the appearance of the particular pixel it hits can be determined. In practice, this proves to be very time-consuming, because the majority of rays from the light source will not pass through the screen. Therefore, in practice ray tracing is performed in reverse, by starting with each pixel on the visible screen and tracing it backwards. Each ray is reflected off surfaces of the object until it reaches the light source. This automatically caters for shadowing and reflections, but not refraction.
The way in which a ray is reflected off an object is determined by a number of factors, and the program uses a complex expression to model this. Basically, the ray is reflected about the normal to the surface at the point of intersection, with an intensity determined by the reflection model and the cosine of the solid angle between the reflected ray and the vector from the surface to the light source. Again, the recommended book explains this in great detail. To add more realism, the specular reflection model, as this function is called, allows for ambient light, that is, overall light that comes from no particular direction. However, ray tracing cannot in any way cope with the reflection of ambient light - an effect which occurs in real life.
Despite its complexity, it is relatively easy to modify the program to handle any image made up of a series of spheres. The procedures which must be altered are PROCinit and PROCSetSpheres.
PROCinit sets up a number of parameters relating to the overall sequence. The three elements of the array ls1 contain the x, y and z components of the position of the light source. The variable maxframe should be set to the total number of frames minus one. The values of width and height determine the size of area traced (in actual numbers of pixels), while increasing the value of margin limits the viewing area, with an obvious speed increase. VD is the distance of the viewer from the object. NoSpheres is the number of spheres that make up the object, and finally, res is the resolution to which the object is traced. A value greater than one gives a faster trace, but with a loss of definition. You can experiment to see the exact effect of these values.
PROCSetSpheres is the routine that actually draws the object in 3-space. This routine is called each time a frame is about to be rendered, and is passed the current frame number. The routine must then create the object by setting up the contents of twelve arrays, described below. The change of object from frame to frame must also be catered for in this routine. If you look at the program, you will see that it starts by drawing the base image into the arrays, and then applying transformations to it according to the frame number.
Each sphere is represented by a set of corresponding elements in the twelve arrays, the array subscript being the sphere number. The arrays are:
x( ) | The position of the centre of the | |
y( ) | sphere in 3-space | |
z( ) | ||
sz( ) | The square of the radius of the sphere | |
lr( ) | The red, green and blue | |
lg( ) | components of the colour of the | |
lb( ) | sphere. These are real values between 0 and 1. | |
ka( ) | Ambient reflectance | |
kn( ) | Normal illumination | |
ksp( ) | Coefficients of specular | |
kr( ) | reflection | |
ks( ) |
These last five arrays are the parameters used for the illumination and specular reflection models. Their meanings are fully explained in the recommended book, but I suggest that unless you understand them fully, you stick to the values used in the program. This will give the spheres an appearance not dissimilar to that of vinyl silk paint. If you follow the SetSpheres procedure in the program, the above explanation should become clearer.
There are many other ways in which you might like to try and modify the program. For example, you could improve the illumination model, or allow for more than one light source, or even allow for the effects of refraction which would make transparent objects possible. If you really want a challenge then you could make the program cater for non-spherical objects, but be warned, this is not a trivial task!
This month's magazine disc contains the two programs given here, and also the version of the program used to produce the famous molecule.
10 REM > Trace
20 REM Program Ray Tracer
30 REM Version A 1.00
40 REM Author Roger Wilson
50 REM RISC User April 1989
60 REM Program Subject to Copyright
70 :
80 MODE 13:OFF:ORIGIN 512,320
90 DIM ls1(2),ls2(2),l(2),ol(2),p(2),np(2),q(2),n(2),s(2),colour(2)
100 PROCinit
110 DIM SA% SASIZE%:SA%!0=SASIZE%
120 SA%!4=0:SA%!8=16:SA%!12=16
130 resX%=4*res: resY%=4*res
140 DIM sz(NoSpheres-1)
150 DIM lr(NoSpheres-1)
160 DIM lg(NoSpheres-1)
170 DIM lb(NoSpheres-1)
180 DIM ka(NoSpheres-1)
190 DIM kn(NoSpheres-1)
200 DIM ks(NoSpheres-1)
210 DIM ksp(NoSpheres-1)
220 DIM kr(NoSpheres-1)
230 DIM x(NoSpheres-1)
240 DIM y(NoSpheres-1)
250 DIM z(NoSpheres-1)
260 :
270 DIM qX(NoSpheres-1),qY(NoSpheres-1),qZ(NoSpheres-1)
280 DIM t(NoSpheres-1), b(NoSpheres-1), c(NoSpheres-1)
290 SO%=&2E:SO_New=&109:SO_Load=&10A
300 SO_Save=&10C:SO_Get=&110
310 INPUT"Start Frame ";SF
320 IFSF<0 THEN
330 A%=OPENIN"StopData"
340 INPUT#A%,SF
350 CLOSE#A%
360 PROCTMDLOAD("RaySprites")
370 ENDIF
380 PROCASS
390 ToSprite%=FALSE
400 ON ERROR PROCCAPTURE
410 linetime%=0
420 FOR frame=SF TO maxframe
430 PROCHDR
440 Start=TIME
450 PROCSetSpheres
460 PROCProjectScene
470 PROCDrawScene:PRINTTIME-Start
480 SpriteName$=RIGHT$("000"+STR$frame,3)
490 SYS "OS_SpriteOp",SO_Get,SA%,SpriteName$,0,0,0,252,252
500 NEXT
510 PROCTMDSAVE("RaySprites")
520 END
530 :
540 DEF PROCCAPTURE
550 ON ERROR OFF
560 IFERR<>17 REPORT:PRINT" at line ";ERL:END
570 MODE0
580 A%=OPENOUT"StopData"
590 PRINT#A%,frame
600 CLOSE#A%
610 PROCTMDSAVE("RaySprites")
620 END
630 :
640 DEF PROCDrawScene
650 FOR Y%=BM TO TM STEP resY%
660 PRINT TAB(cx%,cy%);Y%
670 IF linetime% THEN PRINT"Estimated time to completion=";((TM-Y%)DIVresY%+1)*linetime%DIV100" sec. "
680 PRINT"Doing this frame for ";(TIME-Start)/100" sec. "
690 T1%=TIME
700 FOR X%=LM TO RM STEP resX%
710 IF POINT(X%,Y%) THEN
720 PROCgetlp(X%,Y%)
730 PROCcheck:IFH%THEN
740 colour(0)=2/16:colour(1)=8/16
750 colour(2)=8/16
760 ELSE
770 k=0:colour()=0:ik=1
780 REPEAT
790 PROCreflect:k+=1
800 PROCcheck
810 UNTIL H% OR k>40
820 ENDIF
830 PROCplot
840 ELSE
850 PROCplot2
860 ENDIF
870 NEXT
880 linetime%=(linetime%+(TIME-T1%))/2
890 NEXT
900 ENDPROC
910 :
920 DEF PROCgetlp(X,Y)
930 REM move pixel on view screen to image space
940 l(0)=X-XC:l(1)=Y-YC:l(2)=VD
950 p(0)=l(0):p(1)=l(1):p(2)=0
960 l()=l()/MODl()
970 ENDPROC
980 :
990 DEF PROCSetSpheres
1000 REM initialise image space
1010 RESTORE
1020 IFframe<40 an=frame/40*PI:tilt1=0 ELSE an=PI
1030 IFframe>39 IF frame<90 tilt1=(frame-40)/50*PI
1040 IFframe>89 tilt1=PI:tilt2=(frame-90)/40*PI ELSE tilt2=0
1050 factor=34 : REM approx factor to calculate atomic radii
1060 N=0:C=1:RESTORE
1070 FORY=6TO0STEP-1
1080 X=0:READA$:FORJ%=1TOLENA$
1090 IFMID$(A$,J%,1)="1" THEN
1100 x(N)=X*40-120:y(N)=Y*40-120:z(N)=0
1110 CASE C OF
1120 WHEN 1:lr(N)=1:lg(N)=.25:lb(N)=.25
1130 WHEN 2:lr(N)=.25:lg(N)=1:lb(N)=.25
1140 WHEN 3:lr(N)=1:lg(N)=1:lb(N)=.25
1150 WHEN 4:lr(N)=.25:lg(N)=.25:lb(N)=1
1160 WHEN 5:lr(N)=1:lg(N)=.25:lb(N)=1
1170 WHEN 6:lr(N)=.25:lg(N)=1:lb(N)=1
1180 WHEN 7:lr(N)=1:lg(N)=1:lb(N)=1
1190 ENDCASE
1200 ka(N)=.2:kn(N)=.6:ks(N)=.55:ksp(N)=30:kr(N)=1/3
1210 sz(N)=factor^2
1220 N+=1:C+=1:IFC>7 C=1
1230 ENDIF
1240 X+=1:NEXT
1250 NEXT
1260 cosan=COSan:sinan=SINan
1270 cost2=COStilt1:sint2=SINtilt1
1280 IFframe>39IFframe<90 cost2a=COS(2*PI-tilt1):sint2a=SIN(2*PI-tilt1) ELSE cost2a=cost2:sint2a=sint2
1290 cost1=COStilt2:sint1=SINtilt2
1300 FORN=0TONoSpheres-1
1310 t1=x(N)*cosan-z(N)*sinan:REM spin
1320 t2=x(N)*sinan+z(N)*cosan
1330 x(N)=t1:z(N)=t2
1340 t1=y(N)*cost1-z(N)*sint1:REM tilt
1350 t2=y(N)*sint1+z(N)*cost1
1360 y(N)=t1:z(N)=t2
1370 IFN AND1 cos2=cost2:sin2=sint2 ELSE cos2=cost2a:sin2=sint2a
1380 t1=y(N)*cos2-x(N)*sin2:REM rotate
1390 t2=y(N)*sin2+x(N)*cos2
1400 y(N)=t1:x(N)=t2
1410 NEXT
1420 z()+=1000
1430 ENDPROC
1440 DATA 0000001
1450 DATA 0000011
1460 DATA 0000101
1470 DATA 0001111
1480 DATA 0010001
1490 DATA 0100001
1500 DATA 1000001
1510 :
1520 DEF PROCProjectScene
1530 REM Show gray shadow of scene
1540 GCOL 0,32+8+2 TINT 0:k=VD
1550 FORN=0 TO NoSpheres-1
1560 z=z(N)+VD:CIRCLE FILL k*x(N)/z+XC,k*y(N)/z+YC,k*SQR(sz(N))/z+4
1570 NEXT
1580 ENDPROC
1590 :
1600 DEF PROCplot
1610 REM plot final coloured pixel
1620 IFABS(colour(0)-colour(1))<1/32ANDABS(colour(1)-colour(2))<1/32THEN
1630 r%=FNB(colour(0)*16):g%=r%:b%=r%:t%=r%
1640 ELSE
1650 r%=FNB(colour(0)*16)
1660 g%=FNB(colour(1)*16)
1670 b%=FNB(colour(2)*16)
1680 IFABS(colour(0)-colour(1))<1/32 r%=g%
1690 IFABS(colour(0)-colour(2))<1/32 r%=b%
1700 IFABS(colour(1)-colour(2))<1/32 b%=g%
1710 t%=g%:IFt%+1
1740 GCOL((r%AND12)>>2)+(g%AND12)+((b%AND12)<<2) TINT t%<<6
1750 IFresX%=1POINTX%,Y% ELSERECTANGLEFILLX%,Y%,resX%-4,resY%-4
1760 ENDPROC
1770 :
1780 DEF FNB(X) LOCALX%:X%=INTX:IFRND(1)
1800 =X%
1810 :
1820 DEF PROCplot2
1830 REM plot background colour pixel
1840 GCOL40 TINT 0
1850 IFresX%=1POINTX%,Y% ELSERECTANGLEFILLX%,Y%,resX%-4,resY%-4
1860 ENDPROC
1870 :
1880 DEF PROCcheck
1890 REM compute intersection of normalised vector l() from position p()
1900 REM Returns H%=TRUE if no hit, else H%=FALSE, L=distance, M%=sphere number
1910 H%=TRUE:qX()=p(0)-x():qY()=p(1)-y():qZ()=p(2)-z()
1920 b()=l(0)*qX():t()=l(1)*qY():b()=b()+t():t()=l(2)*qZ():b()=b()+t()
1930 c()=qX()*qX():t()=qY()*qY():c()=c()+t():t()=qZ()*qZ():c()=c()+t()
1940 CALL CODE%,c(),sz(),b(),H%,L,M%
1950 ENDPROC
1960 :
1970 REM Code performs:
1980 FORN%=0TONoSpheres-1:IFt(N%)<0.0NEXT:ENDPROC
1990 t=SQRt(N%):K=-b(N%)-t:t-=b(N%):IFK<t t=K
2000 IFt<1E-3NEXT:ENDPROC
2010 IFH%L=t:M%=N%:H%=FALSE ELSEIFt<L L=t:M%=N%
2020 NEXT:ENDPROC
2030 :
2040 DEF PROCreflect
2050 REM compute new reflected ray. Also compute colour effects.
2060 s(0)=x(M%):s(1)=y(M%):s(2)=z(M%)
2070 q()=l()*L:np()=p()+q():n()=np()-s()
2080 n()=n()/MODn()
2090 q()=n()*l():q()=(2*SUMq())*n()
2100 ol()=l():l()=ls1():d=L:M1%=M%:s()=1E-4*n():p()=np()+n():PROCcheck
2110 IFH%=TRUE THEN
2120 s()=n()*l():s()=(2*SUMs())*n():l()=l()-s():l()=l()/MODl()
2130 s()=n()*ls1():B1=kn(M1%)*SUMs():IFB1<0 B1=0
2140 s()=l()*ol():B2=ks(M1%)*SUMs()^ksp(M1%):IFB2<0 B2=0
2150 s(0)=lr(M1%):s(1)=lg(M1%):s(2)=lb(M1%)
2160 s()=s()*(ka(M1%)+B1):s()=s()+B2
2170 ELSE
2180 s(0)=lr(M1%):s(1)=lg(M1%):s(2)=lb(M1%):s()=s()*ka(M1%)
2190 ENDIF
2200 s()=s()*ik
2210 colour()=colour()+s():ik=ik*kr(M1%)
2220 l()=ol()-q():l()=l()/MODl():H%=FALSE:p()=np()
2230 ENDPROC
2240 :
2250 DEF PROCASS
2260 DIM CODE% 2500
2270 P%=CODE%
2280 [ OPT 0
2290 STR R14,[R0]:MOV PC,R14
2300 ]:A%=CODE%+20:CALLCODE%
2310 F1STA=CODE%!20+22*4:F1LDA=F1STA+4
2320 F1ADD=F1LDA+4:F1XSUB=F1ADD+4
2330 F1MUL=F1XSUB+4:F1XDIV=F1MUL+4
2340 FLOATB=F1XDIV+4:FINTEGB=FLOATB+4:FSQRT=FINTEGB+4
2350 PROCB(F1LDA):PROCB(F1MUL):PROCB(F1ADD)
2360 PROCB(F1XSUB):PROCB(FSQRT):PROCB(F1STA)
2370 FACC=0:FACCX=1:FGRD=2:FSIGN=3:SP=13
2380 REM Call with CALL CODE%,c(),sz(),b(),H%,L,M%
2390 REM c array, Sz array, b array, H%
hit/nohit, L length, n which sphere
2400 ADDRn=0:ADDRL=ADDRn+8:ADDRH=ADDRL+8:ADDRb=ADDRH+8:ADDRSz=ADDRb+8:ADDRc=ADDRSz+8
2410 FORZ=0TO2STEP2
2420 P%=CODE%
2430 [OPT Z
2440 STMFD SP !,{R0,R14}
2450 MOV R10,R9 ;arg list
2460 MOV R11,#0 ;N%
2470 LDR R12,[R10,#ADDRb] ;b()
2480 LDR R12,[R12]
2490 ADD R12,R12,#12 ;base of b()
2500 .LOOP ADD R9,R11,R11,LSL #2
2510 ADD R9,R9,R12 ;address% of b(N%)
2520 BL F1LDA
2530 BL F1MUL
2540 LDR R7,[R10,#ADDRSz]
2550 LDR R7,[R7]
2560 ADD R7,R7,#12
2570 ADD R6,R11,R11,LSL #2
2580 ADD R9,R6,R7
2590 BL F1ADD
2600 LDR R7,[R10,#ADDRc]
2610 LDR R7,[R7]
2620 ADD R7,R7,#12
2630 ADD R6,R11,R11,LSL #2
2640 ADD R9,R6,R7
2650 BL F1XSUB
2660 TEQ FSIGN,#0
2670 BPL LOOPEND
2680 MOV FSIGN,#0
2690 BL FSQRT
2700 ADD R6,R11,R11,LSL #2
2710 ADD R9,R6,R12 ;address% of b(N%)
2720 BL F1ADD
2730 TEQ FACC,#0
2740 EORNE FSIGN,FSIGN,#&80000000
2750 TEQ FSIGN,#0
2760 BMI LOOPEND
2770 CMP FACCX,#&70
2780 BCC LOOPEND ;if t<1e-5 break
2790 LDR R7,[R10,#ADDRH]
2800 LDR R6,[R7]
2810 CMP R6,#0
2820 BEQ NOHIT
2830 MOV R6,#0 ;if H%
2840 STR R6,[R7] ;H%=0
2850 LDR R9,[R10,#ADDRL]
2860 BL F1STA ;L=t
2870 STR R11,[SP] ;M%=N%
2880 B LOOPEND
2890 .NOHIT STMFD SP !,{R0,R1,R2,R3} ;else
2900 LDR R9,[R10,#ADDRL]
2910 BL F1LDA
2920 LDMFD SP !,{R4,R5,R6,R7}
2930 BL FCMP ;test L>t
2940 BCC LOOPEND ;ift<L
2950 MOV R0,R4
2960 MOV R1,R5
2970 MOV R2,R6
2980 MOV R3,R7
2990 BL F1STA ;L=t
3000 STR R11,[SP] ;M%=N%
3010 .LOOPEND ADD R11,R11,#1
3020 LDR R7,[R12,#-12]
3030 CMP R7,R11
3040 BCS LOOP
3050 LDMFD SP !,{R0}
3060 LDR R9,[R10,#ADDRn]
3070 STR R0,[R9]
3080 CMP R0,R0
3090 LDMFD SP !,{PC}
3100 .FCMP CMP (FSIGN+4),FSIGN ;check r0-r3 against r4-r7
3110 BNE FCMPDONE
3120 TEQ FSIGN,#0
3130 BMI FCMPNE
3140 CMP FACCX,(FACCX+4)
3150 CMPEQ FACC,(FACC+4)
3160 MOV PC,R14
3170 .FCMPNE CMP (FACCX+4),FACCX
3180 CMPEQ (FACC+4),FACC
3190 .FCMPDONE MOV PC,R14
3200 ]
3210 NEXT
3220 ENDPROC
3230 :
3240 DEF PROCB(RETURN A)
3250 A=A+8+(!A AND&FFFFFF)*4
3260 ENDPROC
3270 :
3280 DEF PROCTMDSAVE(A$)
3290 SYS "OS_SpriteOp",SO_Save,SA%,A$
3300 ENDPROC
3310 :
3320 DEF PROCTMDLOAD(A$)
3330 SYS "OS_SpriteOp",SO_Load,SA%,A$
3340 ENDPROC
3350 :
3360 DEF PROCHDR
3370 CLS
3380 PRINT"Computing frame ";frame
3390 PRINT"Y coordinate=";:cx%=POS:cy%=VPOS
3400 ENDPROC
3410 :
3420 DEF PROCinit
3430 ls1(0)=-30:ls1(1)=20:ls1(2)=-30
3440 ls1()=ls1()/MODls1()
3450 margin=0:maxframe=129
3460 width=64:height=64
3470 LM=0+margin:RM=width*4-1-margin
3480 BM=0+margin:TM=height*4-1-margin
3490 XC=width*2:YC=height*2:VD=3800
3500 NoSpheres=15
3510 SASIZE%=130*(width*height+256)
3520 res=1
3530 ENDPROC
10 REM >Display
20 REM Display ray traced sequence
30 X%=OPENIN"RaySprites":T%=EXT#X%+16
40 CLOSE#X%:DIM sp T%:!sp=T%:sp!4=0
50 sp!8=16:sp!12=16
60 SYS "OS_SpriteOp",&10A,sp,"RaySprites"
70 MODE 13:OFF
80 REPEAT
90 FOR frame=0 TO 129
100 SYS "OS_SpriteOp",&122,sp,RIGHT$("000"+STR$frame,3),512,384
110 WAIT:WAIT
120 NEXT
130 UNTIL FALSE