Q     Why is the tracer written in Oberon? Why does it run on an obscure OS rather than Linux or Windows?

A     Who are you and why are you bothering me?

http://www.youtube.com/watch?v=sZDCRgvjgbs

This shows the normal-shaded spheres vasterized into mspace by the code in the modules below. Extremes of scale and detail are possible with this tracer. I was deeply influenced by JJ Gibson's ecological approach to visual perception and wanted a realtime renderer that could handle any amount of geometric detail. Socalled photorealism was not the object. 

At first I thought of explaining the tracer in the manner of a teacher, starting with theoretical motivations and simple examples, illustrated helpfully. Realizing in time that as I anticipated no readers, or perhaps only a single ideal reader, and in either case the effort at pedagogical soundness would be wasted, I begin with two complete modules of code: neither the first written nor the first imported at runtime; yet if I were to explain this code adequately (which I shall not) the whole of the tracer would be made lucid.

MODULE srM2Space;
IMPORT srBase, srE, srFRep, Out := DebugLog;

CONST
M=2;
MMO=1;
CUBERADIUS= 0.7;  (*change to sqrt(3)/2 ???*)

TYPE PT = srBase.PT;
TYPE COLOR = srBase.COLOR;
TYPE Ray = srBase.Ray;
TYPE Voxel = srBase.Voxel;
TYPE FR=srBase.FR;;
TYPE NCUBE=RECORD
filled,death,passable: BOOLEAN;
normal: PT;
color:COLOR;
mirror:REAL;
END;

TYPE ARM = ARRAY M,M,M OF Voxel;
TYPE NRM = ARRAY M,M,M OF NCUBE;

TYPE cell* = OBJECT(srFRep.MSV);

VAR
blox*: ARM;
nblox:  NRM;
red,green,blue,black: REAL;
swallow*: BOOLEAN;

PROCEDURE & init;
BEGIN
iscell:=TRUE;
complex:=TRUE;
passable:=TRUE;
swallow:=FALSE;
allowdeathray:=TRUE;
END init;

PROCEDURE camnotify;
BEGIN
IF swallow THEN Out.Char('!'); srBase.world:=SELF END
END camnotify;

PROCEDURE setcolor* (r,g,b,bl: REAL);
BEGIN
red:= r;
green := g;
blue := b;
black:= bl;
END setcolor;

PROCEDURE bounds* (i, j, k: LONGINT; VAR out: BOOLEAN);
BEGIN
IF (i < 0) OR (i > MMO) OR (j < 0) OR (j > MMO) OR (k < 0) OR (k > MMO) THEN
out := TRUE
ELSE
out := FALSE
END
END bounds;

PROCEDURE fill*(v: Voxel);
VAR
i,j,k: INTEGER;
BEGIN
FOR i := 0 TO MMO DO FOR j := 0 TO MMO DO FOR k:= 0 TO MMO DO
blox[i,j,k] := v
END END END
END fill;

PROCEDURE serp*(v,w: Voxel);
BEGIN
blox[0,0,0]:=v;
blox[0,1,0]:=v;
blox[1,0,0]:=v;
blox[1,1,0]:=v;
blox[0,0,1]:=w;
blox[0,1,1]:=w;
blox[1,0,1]:=w;
blox[1,1,1]:=w;
END serp;

PROCEDURE slant*(u,v,w: Voxel);
BEGIN
fill(v);
blox[0,0,0]:=v;
blox[1,1,1]:=w;
END slant;

PROCEDURE fillwithprobability*(v: Voxel; p: REAL);
VAR
i,j,k: INTEGER;
BEGIN
FOR i := 0 TO MMO DO FOR j := 0 TO MMO DO FOR k:= 0 TO MMO DO
IF srBase.rand.Uniform()<p THEN blox[i,j,k] := v END
END END END
END fillwithprobability;

PROCEDURE fillchequer*(v,w: Voxel);
VAR
i,j,k: INTEGER;
BEGIN
FOR i := 0 TO MMO DO FOR j := 0 TO MMO DO FOR k:= 0 TO MMO DO
IF ODD(i+j+k) THEN blox[i,j,k] := v ELSE blox[i,j,k] := w END
END END END
END fillchequer;

PROCEDURE fillcqr2*(v,w: Voxel);
VAR
i,j,k: INTEGER;
BEGIN
FOR i := 0 TO MMO DO FOR j := 0 TO MMO DO FOR k:= 0 TO MMO DO
fillchequer(v,w)
END END END
END fillcqr2;

PROCEDURE fillcqr3*(v,w: Voxel);
VAR
i,j,k: INTEGER;
BEGIN
FOR i := 0 TO MMO DO FOR j := 0 TO MMO DO FOR k:= 0 TO MMO DO
fillcqr2(v,w)
END END END
END fillcqr3;

PROCEDURE ncolor(VAR ray: Ray; cube:NCUBE);
VAR
dot,r,g,b: REAL;
PROCEDURE reflect(VAR m,n:PT);
BEGIN
dot := m.x*n.x+m.y*n.y+m.z*n.z;
n.x:= 2*n.x*dot; n.y := 2*n.y*dot; n.z := 2*n.z*dot;
m.x := m.x-n.x; m.y := m.y-n.y; m.z := m.z-n.z;
END reflect;

BEGIN
IF (cube.mirror>0.1)THEN
reflect(ray.dxyz,cube.normal);
dot:=ABS(dot);
r:= (cube.color.red * ray.ra*dot*cube.mirror);
g:= (cube.color.green* ray.ga*dot*cube.mirror);
b:= (cube.color.blue * ray.ba*dot*cube.mirror);
ray.r := ray.r + r;
ray.g := ray.g + g;
ray.b := ray.b + b;
ray.ra:= ray.ra - r - 0.1;
ray.ga:= ray.ga - g - 0.1;
ray.ba:= ray.ba - b - 0.1;
IF ray.dxyz.x < 0 THEN ray.di := FALSE  ELSE ray.di := TRUE END;
IF ray.dxyz.y < 0 THEN ray.dj := FALSE  ELSE ray.dj := TRUE END;
IF ray.dxyz.z < 0 THEN ray.dk := FALSE  ELSE ray.dk := TRUE END;
ELSE
dot := ABS(cube.normal.x*ray.dxyz.x + cube.normal.y*ray.dxyz.y+ cube.normal.z*ray.dxyz.z);
IF dot<0.4 THEN dot:=0.4 END;
ray.r := ray.r + cube.color.red * ray.ra*dot;
ray.g := ray.g + cube.color.green * ray.ga*dot;
ray.b := ray.b + cube.color.blue * ray.ba*dot;
ray.terminate:=TRUE
END
END ncolor;

PROCEDURE Shade (VAR ray: Ray);
VAR
oldxyz, newxyz, xyz: srBase.PT;
ijk: srBase.IPT;
drx, dry, drz, dr,rr,gr,br,blr, X: REAL;
out,shadenil,A,B,C: BOOLEAN;
v: Voxel;
vdepth: REAL;
last:BOOLEAN;
BEGIN
ray.scale := ray.scale*M;
oldxyz := ray.xyz;
X:= ray.length*ray.scale;
IF (imposter#NIL)&(X>srBase.DTL) THEN
imposter.Shade(ray)
ELSE
xyz.x := ray.lxyz.x * M  - ray.ddxyz.x;
xyz.y := ray.lxyz.y * M  - ray.ddxyz.y;
xyz.z := ray.lxyz.z * M  - ray.ddxyz.z;
srE.E(xyz,ijk);
bounds(ijk.i,ijk.j,ijk.k,out);
IF ~out  THEN
v := blox[ijk.i,ijk.j,ijk.k];
IF (v#NIL) THEN
ray.lxyz.x := ABS(xyz.x - ijk.i);
ray.lxyz.y := ABS(xyz.y - ijk.j);
ray.lxyz.z := ABS(xyz.z - ijk.k);
v.Shade(ray);
ELSIF nblox[ijk.i,ijk.j,ijk.k].filled THEN
ncolor(ray,nblox[ijk.i,ijk.j,ijk.k])
END
END;
IF (ray.ra<0.1)&(ray.ga<0.1)&(ray.ba<0.1) THEN ray.terminate:=TRUE END;
IF ~ray.terminate THEN
REPEAT
IF ray.di  THEN
drx := ( (ijk.i + 1) - xyz.x) / ray.dxyz.x
ELSE
drx :=  (ijk.i -  xyz.x) / ray.dxyz.x
END;
IF ray.dj THEN
dry := ( (ijk.j + 1) - xyz.y) / ray.dxyz.y
ELSE
dry :=  (ijk.j - xyz.y) / ray.dxyz.y
END;
IF ray.dk  THEN
drz := ( (ijk.k + 1) - xyz.z) / ray.dxyz.z
ELSE
drz :=  (ijk.k - xyz.z) / ray.dxyz.z
END;
A:=drx<dry; B:=drx<drz; C:=dry<drz;
IF A&B THEN
dr := drx;
IF ray.di THEN
INC(ijk.i, 1);
ray.face := 1; ray.normal:= srBase.Face[0]
ELSE
INC(ijk.i, -1);
ray.face := 4; ray.normal:= srBase.Face[3]
END;
newxyz.x := xyz.x + drx * ray.dxyz.x; newxyz.y := xyz.y + drx * ray.dxyz.y; newxyz.z  := xyz.z + drx * ray.dxyz.z
ELSIF A&~B THEN
dr := drz;
IF ray.dk THEN
INC(ijk.k, 1);
ray.face := 3; ray.normal:= srBase.Face[2]
ELSE
INC(ijk.k, -1);
ray.face := 6; ray.normal:= srBase.Face[5]
END;
newxyz.x := xyz.x + drz * ray.dxyz.x; newxyz.y := xyz.y + drz * ray.dxyz.y; newxyz.z  := xyz.z + drz * ray.dxyz.z
ELSIF C THEN
dr := dry;
IF ray.dj THEN
INC(ijk.j, 1);
ray.face := 2; ray.normal:= srBase.Face[1]
ELSE
INC(ijk.j, -1);
ray.face := 5; ray.normal:= srBase.Face[4]
END;
newxyz.x := xyz.x + dry * ray.dxyz.x; newxyz.y := xyz.y + dry * ray.dxyz.y; newxyz.z  := xyz.z+ dry * ray.dxyz.z
ELSE
dr := drz;
IF ray.dk  THEN
INC(ijk.k, 1);
ray.face := 3; ray.normal:= srBase.Face[2]
ELSE
INC(ijk.k, -1);
ray.face := 6; ray.normal:= srBase.Face[5]
END;
newxyz.x := xyz.x + drz * ray.dxyz.x; newxyz.y := xyz.y + drz * ray.dxyz.y; newxyz.z  := xyz.z + drz * ray.dxyz.z
END;
vdepth:=srBase.distance(newxyz,xyz);
xyz:=newxyz;
ray.length:=ray.length+vdepth/ray.scale;
ray.xyz.x:=ray.xyz.x+xyz.x/ray.scale; ray.xyz.y:=ray.xyz.y+xyz.y/ray.scale; ray.xyz.z:=ray.xyz.z+xyz.z/ray.scale;
IF TRUE THEN
rr := red*vdepth; gr := green*vdepth; br := blue*vdepth; blr:=black*vdepth;
ray.r := ray.r + rr;
ray.g:= ray.g + gr;
ray.b := ray.b + br;
ray.ra := (ray.ra -rr) -blr;
ray.ga := (ray.ga -gr) -blr;
ray.ba := (ray.ba -br )-blr;
END;
bounds(ijk.i,ijk.j,ijk.k, out);
IF ~out THEN
v := blox[ijk.i,ijk.j,ijk.k];
IF (v#NIL)  THEN
ray.lxyz.x := ABS(xyz.x - ijk.i);
ray.lxyz.y := ABS(xyz.y - ijk.j);
ray.lxyz.z := ABS(xyz.z - ijk.k);
v.Shade(ray);
ELSIF nblox[ijk.i,ijk.j,ijk.k].filled THEN
ncolor(ray,nblox[ijk.i,ijk.j,ijk.k])
END;
END;
IF (ray.ra<0.1)&(ray.ga<0.1)&(ray.ba<0.1) THEN ray.terminate:=TRUE END;
UNTIL  ray.terminate OR out;
END END;
ray.scale := ray.scale/M;
ray.xyz := oldxyz;
END Shade;

PROCEDURE probe(x,y,z: REAL):Voxel;
VAR
X,Y,Z: REAL;
i,j,k: LONGINT;
BEGIN
srBase.clamp3(x,y,z);
X := x*M; Y := y*M; Z := z*M;
i := ENTIER(X);
j := ENTIER(Y);
k := ENTIER(Z);
IF nblox[i,j,k].filled THEN RETURN(SELF) END;
IF blox[i,j,k]#NIL THEN RETURN(blox[i,j,k].probe(X-i, Y-j, Z-k)) END;
RETURN(SELF);
END probe;

PROCEDURE passprobe(x,y,z: REAL):BOOLEAN;
VAR
X,Y,Z: REAL;
i,j,k: LONGINT;
BEGIN
srBase.clamp3(x,y,z);
X := x*M; Y := y*M; Z := z*M;
i := ENTIER(X);
j := ENTIER(Y);
k := ENTIER(Z);
IF blox[i,j,k]#NIL THEN RETURN(blox[i,j,k].passprobe(X-i, Y-j, Z-k)) END;
IF nblox[i,j,k].filled THEN RETURN(nblox[i,j,k].passable) END;
RETURN(TRUE);
END passprobe;

PROCEDURE proberay(VAR ray: Ray):Voxel;
VAR
oldxyz: srBase.PT;
ijk: srBase.IPT;
drx, dry, drz, X: REAL;
di, dj, dk: INTEGER;
out: BOOLEAN;
v, probedv: Voxel;
done: BOOLEAN;
BEGIN
ray.scale := ray.scale*M;
oldxyz := ray.xyz;
X:= ray.length*ray.scale;
IF  (X>srBase.DTL) THEN
RETURN(imposter)  (* may be NIL *)
ELSE
ray.xyz.x := ray.lxyz.x * M  - ray.dxyz.x / 1000000 ;
ray.xyz.y := ray.lxyz.y * M  - ray.dxyz.y / 1000000 ;
ray.xyz.z := ray.lxyz.z * M  - ray.dxyz.z / 1000000 ;
srE.E(ray.xyz,ijk);
bounds(ijk.i,ijk.j,ijk.k, out);
NEW(probedv);
IF ~out THEN
v := blox[ijk.i,ijk.j,ijk.k];
IF  v # NIL THEN
IF v.complex THEN
ray.lxyz.x := ABS(ray.xyz.x - ijk.i);
ray.lxyz.y := ABS(ray.xyz.y - ijk.j);
ray.lxyz.z := ABS(ray.xyz.z - ijk.k);
probedv:=v.proberay(ray);
ELSE
probedv:=v;
done:=TRUE
END
END;
END;
IF ~done THEN
IF ray.dxyz.x < 0 THEN di := - 1  ELSE di := 1 END;
IF ray.dxyz.y < 0 THEN dj := - 1  ELSE dj := 1 END;
IF ray.dxyz.z< 0 THEN dk := - 1  ELSE dk := 1 END;
REPEAT
IF di > 0 THEN
drx := ( (ijk.i + 1) - ray.xyz.x) / ray.dxyz.x
ELSE
drx :=  (ijk.i -  ray.xyz.x) / ray.dxyz.x
END;
IF dj > 0 THEN
dry := ( (ijk.j + 1) - ray.xyz.y) / ray.dxyz.y
ELSE
dry :=  (ijk.j - ray.xyz.y) / ray.dxyz.y
END;
IF dk > 0 THEN
drz := ( (ijk.k + 1) - ray.xyz.z) / ray.dxyz.z
ELSE
drz :=  (ijk.k - ray.xyz.z) / ray.dxyz.z
END;
IF (drx < dry) THEN
IF (drx < drz ) THEN
INC(ijk.i, di);
IF di > 0 THEN
ray.face := 1; ray.normal:= srBase.Face[0]
ELSE
ray.face := 4; ray.normal:= srBase.Face[3]
END;
ray.xyz.x := ray.xyz.x + drx * ray.dxyz.x; ray.xyz.y := ray.xyz.y + drx * ray.dxyz.y; ray.xyz.z  := ray.xyz.z + drx * ray.dxyz.z
ELSE
INC(ijk.k, dk);
IF dk > 0 THEN
ray.face := 3; ray.normal:= srBase.Face[2]
ELSE
ray.face := 6; ray.normal:= srBase.Face[5]
END;
ray.xyz.x := ray.xyz.x + drz * ray.dxyz.x; ray.xyz.y := ray.xyz.y + drz * ray.dxyz.y; ray.xyz.z  := ray.xyz.z + drz * ray.dxyz.z
END
ELSIF (dry < drz) THEN
INC(ijk.j, dj);
IF dj > 0 THEN
ray.face := 2; ray.normal:= srBase.Face[1]
ELSE
ray.face := 5; ray.normal:= srBase.Face[4]
END;
ray.xyz.x := ray.xyz.x + dry * ray.dxyz.x; ray.xyz.y := ray.xyz.y + dry * ray.dxyz.y; ray.xyz.z  := ray.xyz.z+ dry * ray.dxyz.z
ELSE
INC(ijk.k, dk);
IF dk > 0 THEN
ray.face := 3; ray.normal:= srBase.Face[2]
ELSE
ray.face := 6; ray.normal:= srBase.Face[5]
END;
ray.xyz.x := ray.xyz.x + drz * ray.dxyz.x; ray.xyz.y := ray.xyz.y + drz * ray.dxyz.y; ray.xyz.z  := ray.xyz.z + drz * ray.dxyz.z
END;
bounds(ijk.i,ijk.j,ijk.k, out);
IF ~out THEN
v := blox[ijk.i,ijk.j,ijk.k];
IF v #NIL THEN
IF v.complex THEN
ray.lxyz.x := ABS(ray.xyz.x - ijk.i);
ray.lxyz.y := ABS(ray.xyz.y - ijk.j);
ray.lxyz.z := ABS(ray.xyz.z - ijk.k);
probedv:= v.proberay(ray);
ELSE
probedv:=v;
done:=TRUE
END;
END;
END;
UNTIL  out;
END END;
ray.scale := ray.scale*M;
ray.xyz := oldxyz;
IF probedv=NIL THEN
NEW(probedv)
END;
RETURN v;
END proberay;


PROCEDURE deathray(VAR ray: Ray);
VAR
oldxyz: srBase.PT;
ijk: srBase.IPT;
drx, dry, drz: REAL;
di, dj, dk: INTEGER;
out: BOOLEAN;
v: Voxel;
killed: BOOLEAN;
BEGIN
oldxyz := ray.xyz;
ray.scale := ray.scale/M;
ray.xyz.x := ray.lxyz.x * M  - ray.dxyz.x / 1000000 ;
ray.xyz.y := ray.lxyz.y * M  - ray.dxyz.y / 1000000 ;
ray.xyz.z := ray.lxyz.z * M  - ray.dxyz.z / 1000000 ;
srE.E(ray.xyz,ijk);
bounds(ijk.i,ijk.j,ijk.k, out);
IF ~out THEN
v := blox[ijk.i,ijk.j,ijk.k];
IF  v # NIL THEN
IF v.complex THEN
ray.lxyz.x := ABS(ray.xyz.x - ijk.i);
ray.lxyz.y := ABS(ray.xyz.y - ijk.j);
ray.lxyz.z := ABS(ray.xyz.z - ijk.k);
v.deathray(ray);
IF ray.changed THEN killed := TRUE END;
END
END
END;
IF ~killed THEN
IF ray.dxyz.x < 0 THEN di := - 1  ELSE di := 1 END;
IF ray.dxyz.y < 0 THEN dj := - 1  ELSE dj := 1 END;
IF ray.dxyz.z< 0 THEN dk := - 1  ELSE dk := 1 END;
REPEAT
IF di > 0 THEN
drx := ( (ijk.i + 1) - ray.xyz.x) / ray.dxyz.x
ELSE
drx :=  (ijk.i -  ray.xyz.x) / ray.dxyz.x
END;
IF dj > 0 THEN
dry := ( (ijk.j + 1) - ray.xyz.y) / ray.dxyz.y
ELSE
dry :=  (ijk.j - ray.xyz.y) / ray.dxyz.y
END;
IF dk > 0 THEN
drz := ( (ijk.k + 1) - ray.xyz.z) / ray.dxyz.z
ELSE
drz :=  (ijk.k - ray.xyz.z) / ray.dxyz.z
END;
IF (drx < dry) THEN
IF (drx < drz ) THEN
INC(ijk.i, di);
IF di > 0 THEN
ray.face := 1; ray.normal:= srBase.Face[0]
ELSE
ray.face := 4; ray.normal:= srBase.Face[3]
END;
ray.xyz.x := ray.xyz.x + drx * ray.dxyz.x; ray.xyz.y := ray.xyz.y + drx * ray.dxyz.y; ray.xyz.z  := ray.xyz.z + drx * ray.dxyz.z
ELSE
INC(ijk.k, dk);
IF dk > 0 THEN
ray.face := 3; ray.normal:= srBase.Face[2]
ELSE
ray.face := 6; ray.normal:= srBase.Face[5]
END;
ray.xyz.x := ray.xyz.x + drz * ray.dxyz.x; ray.xyz.y := ray.xyz.y + drz * ray.dxyz.y; ray.xyz.z  := ray.xyz.z + drz * ray.dxyz.z
END
ELSIF (dry < drz) THEN
INC(ijk.j, dj);
IF dj > 0 THEN
ray.face := 2; ray.normal:= srBase.Face[1]
ELSE
ray.face := 5; ray.normal:= srBase.Face[4]
END;
ray.xyz.x := ray.xyz.x + dry * ray.dxyz.x; ray.xyz.y := ray.xyz.y + dry * ray.dxyz.y; ray.xyz.z  := ray.xyz.z+ dry * ray.dxyz.z
ELSE
INC(ijk.k, dk);
IF dk > 0 THEN
ray.face := 3; ray.normal:= srBase.Face[2]
ELSE
ray.face := 6; ray.normal:= srBase.Face[5]
END;
ray.xyz.x := ray.xyz.x + drz * ray.dxyz.x; ray.xyz.y := ray.xyz.y + drz * ray.dxyz.y; ray.xyz.z  := ray.xyz.z + drz * ray.dxyz.z
END;
bounds(ijk.i,ijk.j,ijk.k, out);
IF ~out THEN
v := blox[ijk.i,ijk.j,ijk.k];
IF v #NIL THEN
IF v.complex & v.allowdeathray THEN
ray.lxyz.x := ABS(ray.xyz.x - ijk.i);
ray.lxyz.y := ABS(ray.xyz.y - ijk.j);
ray.lxyz.z := ABS(ray.xyz.z - ijk.k);
v.deathray(ray);
IF ray.changed THEN killed := TRUE END;
ELSE
blox[ijk.i,ijk.j,ijk.k] := NIL;
killed:=TRUE
END;
END;
IF nblox[ijk.i,ijk.j,ijk.k].death THEN
nblox[ijk.i,ijk.j,ijk.k].filled:=FALSE;
killed:=TRUE
END
END;
UNTIL  killed OR out;
END;
ray.scale := ray.scale*M;
ray.xyz := oldxyz;
END deathray;

PROCEDURE stroke*(p:PT; level: LONGINT; normal:PT; color: COLOR);
BEGIN
IF  (level>=1) & srBase.inzerodotdotonePT(p) THEN
strokerec(p, level,normal,color);
END
END stroke;

PROCEDURE strokergb*(p:PT; level: LONGINT; normal:PT; r,g,b: REAL);
VAR
color: COLOR;
BEGIN
color.red:=r; color.green:=g; color.blue:=b;
IF  (level>=1) & srBase.inzerodotdotonePT(p) THEN
strokerec(p, level,normal,color);
END
END strokergb;

PROCEDURE strokerec(p:PT; level: LONGINT; normal:PT; color: COLOR);
VAR
i,j,k: LONGINT;
c: cell;
BEGIN
pdiv(p,M);
i := ENTIER(p.x); j := ENTIER(p.y); k := ENTIER(p.z);
IF level=1 THEN
(* we're here. *)
nblox[i,j,k].normal:=normal;
nblox[i,j,k].color:=color;
nblox[i,j,k].filled:=TRUE;
ELSE
IF blox[i,j,k]=NIL THEN
NEW(c);
c.setcolor(red,green,blue,black);
blox[i,j,k]:=c;
END;
p.x:=p.x-i; p.y:=p.y-j; p.z:=p.z-k;
blox[i,j,k].strokerec(p, level-1,normal,color);
END
END strokerec;

PROCEDURE strokevoxel*(p:PT; resolution:LONGINT; voxel:Voxel);
VAR
i,j,k: LONGINT;
c: cell;
BEGIN
srBase.clamPT(p);
strokevoxelrec(p,resolution,1,voxel);
END strokevoxel;

PROCEDURE strokevoxelrec*(p:PT; resolution,scale:LONGINT; voxel:Voxel);
VAR
i,j,k,nextscale: LONGINT;
c: cell;
BEGIN
nextscale:=scale*M;
pdiv(p,M);
i := ENTIER(p.x); j := ENTIER(p.y); k := ENTIER(p.z);
IF nextscale>resolution THEN
IF blox[i,j,k]=NIL THEN blox[i,j,k]:=voxel END
ELSE
IF blox[i,j,k] = NIL THEN
NEW(c);
blox[i,j,k]:=c;
END;
p.x:=p.x-i; p.y:=p.y-j; p.z:=p.z-k;
blox[i,j,k].strokevoxelrec(p, resolution,nextscale,voxel);
END
END strokevoxelrec;

PROCEDURE erase*(p:PT; resolution:LONGINT);
VAR
i,j,k: LONGINT;
c: cell;
BEGIN
srBase.clamPT(p);
eraserec(p,resolution,1);
END erase;

PROCEDURE eraserec*(p:PT; resolution,scale:LONGINT);
VAR
i,j,k,nextscale: LONGINT;
c: cell;
BEGIN
nextscale:=scale*M;
pdiv(p,M);
i := ENTIER(p.x); j := ENTIER(p.y); k := ENTIER(p.z);
IF nextscale>resolution THEN
blox[i,j,k]:=NIL
ELSE
IF blox[i,j,k]# NIL THEN
p.x:=p.x-i; p.y:=p.y-j; p.z:=p.z-k;
blox[i,j,k].eraserec(p, resolution,nextscale)
END
END
END eraserec;

PROCEDURE clear*(p:PT; level: LONGINT);
BEGIN
IF  (level>=1) & srBase.inzerodotdotonePT(p) THEN
clearrec(p, level);
END
END clear;

PROCEDURE clearrec(p:PT; level: LONGINT);
VAR
i,j,k: LONGINT;
BEGIN
srBase.clamPT(p);
pdiv(p,M);
i := ENTIER(p.x); j := ENTIER(p.y); k := ENTIER(p.z);
IF level=1 THEN
(* we're here. *)
nblox[i,j,k].filled:=FALSE;
blox[i,j,k]:=NIL
ELSE
IF blox[i,j,k]#NIL THEN
blox[i,j,k].clearrec(p,level-1)
END
END
END clearrec;

PROCEDURE line*(a,b: PT; level: LONGINT; v: Voxel);
VAR
tx,ty,tz, dxdt, dydt, dzdt: REAL;
t: LONGINT;
delta: REAL;
n: LONGINT;
p: PT;
BEGIN
CASE level OF
1: delta := 1/M;
|2: delta := 1/M*M;
| 3: delta := 1/M*M*M;
|4: delta := 1/M*M*M*M;
ELSE
delta := 0;
END;
IF delta > 0 THEN
n := ENTIER(srBase.distance(a,b)/delta);
tx := b.x; ty := b.y; tz := b.z;
dxdt := (a.x-b.x)/n; dydt := (a.y-b.y)/n; dzdt := (a.z-b.z)/n;
FOR t := 0 TO n DO
srBase.setPT(p,tx, ty, tz);
strokevoxel(p, level,v);
tx := tx + dxdt; ty := ty + dydt; tz := tz+dzdt;
END
END
END line;

PROCEDURE linevoxel*(a,b: PT; level: LONGINT; v: Voxel);
VAR
tx,ty,tz, dxdt, dydt, dzdt: REAL;
t: LONGINT;
delta: REAL;
n: LONGINT;
p: PT;

BEGIN
CASE level OF
1: delta := 1/M;
|2: delta := 1/M*M;
| 3: delta := 1/M*M*M;
|4: delta := 1/M*M*M*M;
ELSE
delta := 0;
END;
IF delta > 0 THEN
n := ENTIER(srBase.distance(a,b)/delta);
tx := b.x; ty := b.y; tz := b.z;
dxdt := (a.x-b.x)/n; dydt := (a.y-b.y)/n; dzdt := (a.z-b.z)/n;
FOR t := 0 TO n DO
srBase.setPT(p,tx, ty, tz);
strokevoxel(p, level,v);
tx := tx + dxdt; ty := ty + dydt; tz := tz+dzdt;
END
END
END linevoxel;

PROCEDURE FRaster*( f: FR; resolution: LONGINT);
VAR
origin: PT;
BEGIN
origin.x:=0; origin.y:=0; origin.z:=0;
FRasterrec(f,resolution,origin,1);   (* origin is (0,0,0) *)
END FRaster;

PROCEDURE FRasterrec(f: FR; resolution: LONGINT; origin: PT; scale: LONGINT); (*origin is always in world space*)
VAR
i,j,k: INTEGER;
o,p:PT;
d2s,MS,TWOMS,CRDS,CRDNS:REAL;
nextscale: LONGINT;
v: Voxel;
cell: cell;
BEGIN
MS:=M*scale;
TWOMS:=2*MS;
nextscale:=scale*M;
CRDS:=CUBERADIUS/scale;
CRDNS:=CUBERADIUS/nextscale;
IF nextscale<resolution THEN
FOR i := 0 TO MMO DO FOR j := 0 TO MMO DO FOR k:= 0 TO MMO DO
p.x:=origin.x+(i+1/2)/MS; p.y:=origin.y+(j+1/2)/MS; p.z:=origin.z+(k+1/2)/MS; (*world coordinates*)
d2s:=f.d2s(p);
IF ABS(d2s) < CRDS THEN
o.x:=p.x-1/TWOMS; o.y:=p.y-1/TWOMS; o.z:=p.z-1/TWOMS; (* p is center, o is corner *)
IF blox[i,j,k]=NIL THEN
NEW(cell);
blox[i,j,k]:=cell;
cell.FRasterrec(f,resolution,o,nextscale);
ELSE
v:=blox[i,j,k];  (* compiler disallows type tests and guards on array elements *)
IF v IS srFRep.MSV THEN
WITH v:srFRep.MSV DO
v.FRasterrec(f,resolution,o,nextscale);
END
END
END
END
END END END
END;
FOR i := 0 TO MMO DO FOR j := 0 TO MMO DO FOR k:= 0 TO MMO DO
p.x:=origin.x+(i+1/2)/MS; p.y:=origin.y+(j+1/2)/MS; p.z:=origin.z+(k+1/2)/MS;
d2s:=f.d2s(p);
IF ABS(d2s)<CRDNS THEN
v:=f.voxel(p);
IF v#NIL THEN
blox[i,j,k]:=v
ELSE
nblox[i,j,k].normal:=f.normal(p);
nblox[i,j,k].color:=f.color(p);
nblox[i,j,k].mirror:=f.mirror(p);
nblox[i,j,k].death:=f.death(p);
nblox[i,j,k].passable:=f.pass(p);
nblox[i,j,k].filled:=TRUE;
END;
END;
END END END;
END FRasterrec;
END cell;

PROCEDURE pdiv(VAR p:PT; d:REAL);
BEGIN
p.x:=p.x*d;
p.y:=p.y*d;
p.z:=p.z*d;
END pdiv;

END srM2Space.

 MODULE srFRep;
IMPORT Math, srBase,srVoxel,Raster, Graphics:=WMGraphics;

TYPE
PT=srBase.PT;
COLOR=srBase.COLOR;
Voxel=srBase.Voxel;
Name = srBase.Name;
FR=srBase.FR;

(* Cubes first, then spheres *)
(*Rectangular prisms?*)

TYPE MSV*= OBJECT(Voxel); (* Base type for all mspace cell types. *)

PROCEDURE FRasterrec*(f: FR; resolution: LONGINT; origin: PT; scale: LONGINT);
END FRasterrec;

END MSV;

TYPE rprism= OBJECT(FR)

PROCEDURE&init(c1,c2:PT);
BEGIN
IF c1.x<c2.x THEN bbc1.x:=c1.x; bbc2.x:=c2.x ELSE bbc1.x:=c2.x; bbc2.x:=c1.x END;
IF c1.y<c2.y THEN bbc1.y:=c1.y; bbc2.y:=c2.y ELSE bbc1.y:=c2.y; bbc2.y:=c1.y END;
IF c1.z<c2.z THEN bbc1.z:=c1.z; bbc2.z:=c2.z ELSE bbc1.z:=c2.z; bbc2.z:=c1.z END;
END init;

PROCEDURE in*(p:PT):BOOLEAN;
VAR
a,b,c: BOOLEAN;
BEGIN
a:= ((bbc1.x<p.x)&(p.x<bbc2.x))OR((bbc1.x>p.x)&(p.x>bbc2.x)) ;
b:= ((bbc1.y<p.y)&(p.y<bbc2.y))OR((bbc1.y>p.y)&(p.y>bbc2.y)) ;
c:= ((bbc1.x<p.z)&(p.z<bbc2.z))OR((bbc1.z>p.z)&(p.z>bbc2.z)) ;
RETURN(a&b&c);
END in;

END rprism;

TYPE sphere= OBJECT(FR)
VAR
radius*, rsquared*: REAL;
center*:PT;

PROCEDURE &init(c:PT; r: REAL);
BEGIN
radius:=ABS(r);
rsquared:=r*r;
r:=r*1.00001; (* don't want the sphere to quite touch the box *);
bbc1.x:=c.x-r;
bbc1.y:=c.y-r;
bbc1.z:=c.z-r;
bbc2.x:=c.x+r;
bbc2.y:=c.y+r;
bbc2.z:=c.z+r;
center:=c;
END init;

PROCEDURE reset*(c:PT; r: REAL);
BEGIN
radius:=ABS(r);
rsquared:=r*r;
r:=r*1.00001; (* don't want the sphere to quite touch the box *);
bbc1.x:=c.x-r;
bbc1.y:=c.y-r;
bbc1.z:=c.z-r;
bbc2.x:=c.x+r;
bbc2.y:=c.y+r;
bbc2.z:=c.z+r;
center:=c;
END reset;

PROCEDURE setcenter*(x,y,z: REAL);
BEGIN
center.x:=x; center.y:=y; center.z:=z
END setcenter;

PROCEDURE d2s*(p:PT):REAL;  (* distance to surface *)
BEGIN
norml.x:= center.x-p.x;
norml.y:= center.y-p.y;
norml.z:= center.z-p.z;
RETURN((norml.x*norml.x+norml.y*norml.y+norml.z*norml.z)-rsquared)
END d2s;

END sphere;

TYPE ellipsoid*= OBJECT(FR)
VAR
radius: REAL;
center0,center1,norml0,norml1:PT;

PROCEDURE &init(a,b:PT; m: REAL);

BEGIN
center0:=a;
center1:=b;
radius:=m;
END init;

PROCEDURE setcenter*(x,y,z,a,b,c: REAL);
BEGIN
center0.x:=x; center0.y:=y; center0.z:=z;
center1.x:=x; center1.y:=y; center1.z:=z
END setcenter;

PROCEDURE d2s*(p:PT):REAL;  (* distance to surface *)
VAR
d,e:REAL;
BEGIN
norml0.x:= center0.x-p.x;
norml0.y:= center0.y-p.y;
norml0.z:= center0.z-p.z;
norml1.x:= center1.x-p.x;
norml1.y:= center1.y-p.y;
norml1.z:= center1.z-p.z;
d:= Math.sqrt(norml0.x*norml0.x+norml0.y*norml0.y+norml0.z*norml0.z);
e:= Math.sqrt(norml1.x*norml1.x+norml1.y*norml1.y+norml1.z*norml1.z);
norml.x:=(norml0.x+norml1.x)/2;
norml.y:=(norml0.y+norml1.y)/2;
norml.y:=(norml0.z+norml1.z)/2;
RETURN((d+e)-radius)
END d2s;
END ellipsoid;

TYPE cyl*= OBJECT(FR)
VAR
rsquared,lensq: REAL;
p1,p2,d:PT;

PROCEDURE &init(a,b:PT; r: REAL);

BEGIN
rsquared:=r*r;
p1:=a; p2:=b;
d.x:=p2.x-p1.x;
d.y:=p2.y-p1.y;
d.z:=p2.z-p1.z;
lensq:=d.x*d.x+d.y*d.y+d.z*d.z;
srBase.setPT(norml,1,0,0)
END init;

PROCEDURE d2s*(p:PT):REAL;  (* distance to surface *)
VAR
dot,dsq,a,b,c: REAL;
pd:PT;
BEGIN
pd.x:=p.x-p1.x;
pd.y:=p.y-p1.y;
pd.z:=p.z-p1.z;
dot:=pd.x*d.x+pd.y*d.y+pd.z*d.z;
dsq:=(pd.x*pd.x+pd.y*pd.y+pd.z*pd.z)-dot*dot/lensq;
RETURN(dsq-rsquared);
END d2s;

END cyl;

TYPE texturedcone*= OBJECT(FR)
VAR
rsquared,r2squared,lensq: REAL;
p1,p2,d:PT;
img: Raster.Image;
fmt: Raster.Format;
copy : Raster.Mode;
W,H, bpr,adr: LONGINT;

PROCEDURE &init(a,b:PT; r,r2: REAL; n: Name);
BEGIN
rsquared:=r*r;
r2squared:=r2*r2;
p1:=a; p2:=b;
d.x:=p2.x-p1.x;
d.y:=p2.y-p1.y;
d.z:=p2.z-p1.z;
lensq:=d.x*d.x+d.y*d.y+d.z*d.z;
srBase.setPT(norml,1,0,0);
Raster.InitMode(copy, Raster.srcCopy);
img :=Graphics.LoadImage(n, TRUE);
IF img#NIL THEN W := img.width-1; H:= img.height-1; END;
END init;

PROCEDURE color*(p:PT):COLOR;
VAR
a:COLOR;
r,b,g: REAL;
tx,ty: LONGINT;
pixel:Raster.Pixel;
BEGIN
IF img#NIL THEN
tx:=ENTIER(p.z*W);
ty:=ENTIER(p.y*H);
Raster.Get(img,tx,ty,pixel,copy);
r := ORD(pixel[2])/255; g := ORD(pixel[1])/255; b := ORD(pixel[0])/255;
a.red:=r;
a.green:=g;
a.blue:=b;
ELSE
a.red:=c.red;
a.green:=c.green;
a.blue:=c.blue;
END;
RETURN a
END color;

PROCEDURE d2s*(p:PT):REAL;  (* distance to surface *)
VAR
dot,dsq,a,b,c: REAL;
pd:PT;
BEGIN
pd.x:=p.x-p1.x;
pd.y:=p.y-p1.y;
pd.z:=p.z-p1.z;
dot:=pd.x*d.x+pd.y*d.y+pd.z*d.z;
dsq:=(pd.x*pd.x+pd.y*pd.y+pd.z*pd.z)-dot*dot/lensq;
RETURN(dsq-rsquared);
END d2s;

END texturedcone;

TYPE plane*= OBJECT(FR)
VAR
d: REAL;
p:PT;

PROCEDURE &init(a,b:PT);
BEGIN
p:=a;   norml:=b;
norml.x:= norml.x + (srBase.rand.Uniform()-1/2)/20;
norml.y:= norml.y + (srBase.rand.Uniform()-1/2)/20;
norml.z:= norml.z + (srBase.rand.Uniform()-1/2)/20;
srBase.normalizePT(norml);
d:= -norml.x*p.x - norml.y*p.y - norml.z*p.z;
END init;

PROCEDURE color*(p:PT):COLOR;
VAR
a:COLOR;
BEGIN
a.red:=p.x;
a.green:=p.y;
a.blue:=p.z;
RETURN a
END color;

PROCEDURE mirror*(p:PT):REAL;
BEGIN
IF srBase.rand.Uniform()>0.9 THEN
RETURN(1/2)
ELSE
RETURN(0)
END
END mirror;

PROCEDURE d2s*(m:PT):REAL;  (* distance to surface *)
VAR
BEGIN
RETURN ABS(m.x*norml.x+m.y*norml.y+m.z*norml.z+d);
END d2s;

PROCEDURE normal*(p:PT):PT;
BEGIN
RETURN norml
END normal;

END plane;

TYPE texturedplane*= OBJECT(FR)
VAR
d: REAL;
p*:PT;
img: Raster.Image;
fmt: Raster.Format;
copy : Raster.Mode;
W,H, bpr,adr: LONGINT;

PROCEDURE &init(a,b:PT; n: Name);
BEGIN
p:=a;   norml:=b;
srBase.normalizePT(norml);
d:= -norml.x*p.x - norml.y*p.y - norml.z*p.z;
Raster.InitMode(copy, Raster.srcCopy);
img :=Graphics.LoadImage(n, TRUE);
IF img#NIL THEN W := img.width-1; H:= img.height-1; END;
END init;

PROCEDURE color*(p:PT):COLOR;
VAR
a:COLOR;
r,b,g: REAL;
tx,ty: LONGINT;
pixel:Raster.Pixel;
BEGIN
IF img#NIL THEN
tx:=ENTIER(p.z*W);
ty:=ENTIER(p.y*H);
Raster.Get(img,tx,ty,pixel,copy);
r := ORD(pixel[2])/255; g := ORD(pixel[1])/255; b := ORD(pixel[0])/255;
a.red:=r;
a.green:=g;
a.blue:=b;
ELSE
a.red:=c.red;
a.green:=c.green;
a.blue:=c.blue;
END;
RETURN a
END color;

PROCEDURE d2s*(m:PT):REAL;  (* distance to surface *)
VAR
BEGIN
RETURN ABS(m.x*norml.x+m.y*norml.y+m.z*norml.z+d);
END d2s;

PROCEDURE normal*(p:PT):PT;
BEGIN
RETURN norml
END normal;

END texturedplane;

TYPE cube*= OBJECT(FR)
VAR
white:srVoxel.DiffuseVox;

PROCEDURE color*(p:PT):COLOR;
VAR
c:COLOR;
BEGIN
c.red:=p.x; c.blue:=p.y; c.red:=p.z;
RETURN c
END color;

BEGIN
NEW(white);
white.setcolor(1,1,1);
vox:=white;
END cube;

TYPE cursor* = OBJECT(FR)


END cursor;

TYPE hexgrid = OBJECT(FR)  (*accidentally created from cyl *)
VAR
radius, rsquared,lensq: REAL;
p1,p2,d:PT;

PROCEDURE &init(a,b:PT; r: REAL);

BEGIN
rsquared:=r*r;
p1:=a; p2:=b;
d.x:=p2.x-p1.x;
d.y:=p2.y-p1.y;
d.z:=p2.z-p1.z;
lensq:=d.x*d.x+d.y*d.y+d.z*d.z;
END init;

PROCEDURE d2s*(p:PT):REAL;  (* distance to surface *)
VAR
dot,dsq,a,b,c: REAL;
pd:PT;
BEGIN
pd.x:=p.x-p1.x;
pd.y:=p.y-p1.y;
pd.z:=p.z-p1.z;
dot:=pd.x*p.x+pd.y*p.y+pd.z*p.z;
srBase.normalizePT(pd);
norml:=pd;
(* IF dot<0 THEN
norml.x:= p1.x-p.x;
norml.y:= p1.y-p.y;
norml.z:= p1.z-p.z;
RETURN((norml.x*norml.x+norml.y*norml.y+norml.z*norml.z)-rsquared)
ELSIF dot>lensq THEN
norml.x:= p2.x-p.x;
norml.y:= p2.y-p.y;
norml.z:= p2.z-p.z;
RETURN((norml.x*norml.x+norml.y*norml.y+norml.z*norml.z)-rsquared)
ELSE
dsq:=(pd.x*pd.x+pd.y*pd.y+pd.z*pd.z)-dot*dot/lensq;
RETURN(dsq-rsquared);
END *)
dsq:=(pd.x*pd.x+pd.y*pd.y+pd.z*pd.z)-dot*dot/lensq;
RETURN(dsq-rsquared);
END d2s;

END hexgrid;

TYPE coloredball*=OBJECT(sphere)

PROCEDURE color*(p:PT):COLOR;
BEGIN
RETURN c
END color;

END coloredball;

TYPE blueball*=OBJECT(sphere)
VAR
c,d:COLOR;

PROCEDURE color*(p:PT):COLOR;
BEGIN
RETURN c
END color;

PROCEDURE mirror*(p:PT):REAL;
BEGIN
RETURN 0
END mirror;

BEGIN
c.red:=0.49;
c.green:=0.5;
c.blue:=0.99;
d.red:=0.59;
d.green:=0.5;
d.blue:=0.09;
END blueball;

TYPE yellowball*=OBJECT(sphere)
VAR
c,d:COLOR;

PROCEDURE color*(p:PT):COLOR;
BEGIN
IF ODD(ENTIER(p.x*734)) THEN RETURN c ELSE RETURN d END
END color;

BEGIN
c.red:=0.99;
c.green:=0.89;
c.blue:=0.5;
d.red:=0.79;
d.green:=0.99;
d.blue:=0.5;
END yellowball;

TYPE skyball*=OBJECT(sphere)
VAR
c,d:COLOR;

PROCEDURE color*(p:PT):COLOR;
VAR
b:REAL;
BEGIN
b:=center.z-p.z;
IF b<0 THEN RETURN c ELSE RETURN d END
END color;

BEGIN
c.red:=0.01;
c.green:=0.19;
c.blue:=0.99;
d.red:=0.01;
d.green:=0.09;
d.blue:=0.5;
END skyball;

TYPE yellowcyl*=OBJECT(cyl)
VAR
c,d:COLOR;

PROCEDURE color*(p:PT):COLOR;
BEGIN
RETURN c
END color;

BEGIN
c.red:=0.99;
c.green:=0.89;
c.blue:=0.09;
d.red:=0.89;
d.green:=0.99;
d.blue:=0.5;
END yellowcyl;

TYPE whiteball*=OBJECT(sphere)
VAR
white: srVoxel.DiffuseVox;
c:COLOR;

PROCEDURE color*(p:PT):COLOR;
BEGIN
RETURN c
END color;

BEGIN
c.red:=0.99;
c.green:=0.89;
c.blue:=0.99;
END whiteball;

END srFRep.