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.