DIRECTORY Basics, CedarProcess, FileNames, FS, G2dBasic, G3dBasic, G3dIO, G3dVector, G3dVoxel, IO, Real, RefTab, Rope, RuntimeError;

G3dVoxelImpl: CEDAR PROGRAM
IMPORTS Basics, CedarProcess, FileNames, FS, G3dIO, IO, Real, RefTab, Rope, RuntimeError, G3dVector
EXPORTS G3dVoxel

~ BEGIN

Nat3:						TYPE ~ G2dBasic.Nat3;
Triple:					TYPE ~ G2dBasic.Triple;
RealSequence:			TYPE ~ G2dBasic.RealSequence;
RealSequenceRep:		TYPE ~ G2dBasic.RealSequenceRep;
Box:						TYPE ~ G3dBasic.Box;
Grid:						TYPE ~ G3dVoxel.Grid;
GridRep:					TYPE ~ G3dVoxel.GridRep;
RealProc:					TYPE ~ G3dVoxel.RealProc;
Reals2d:					TYPE ~ G3dVoxel.Reals2d;
Reals3d:					TYPE ~ G3dVoxel.Reals3d;
Reals2dRep:				TYPE ~ G3dVoxel.Reals2dRep;
Reals3dRep:				TYPE ~ G3dVoxel.Reals3dRep;
Refs:						TYPE ~ G3dVoxel.Refs;
Refs2d:					TYPE ~ G3dVoxel.Refs2d;
Refs3d:					TYPE ~ G3dVoxel.Refs3d;
Voxels:					TYPE ~ G3dVoxel.Voxels;
VoxelsRep:				TYPE ~ G3dVoxel.VoxelsRep;
ROPE:						TYPE ~ Rope.ROPE;
Error: PUBLIC SIGNAL [reason: ATOM] = CODE;
CreateVoxels: PUBLIC PROC [sparse: BOOL ¬ TRUE, size: Nat3 ¬ [0, 0, 0], name: ROPE ¬ NIL]
RETURNS [voxels: Voxels]
~ {
voxels ¬ NEW[VoxelsRep ¬ [name: name, type: IF sparse THEN hash ELSE array, size: size]];
IF sparse
THEN voxels.hashTable ¬ RefTab.Create[equal: Equal, hash: Hash]
ELSE {
ENABLE RuntimeError.BoundsFault => Error[$sizeTooGreat];
voxels.array ¬ NEW[Refs3d[size.x]];
FOR i: NAT IN [0..size.x) DO
voxels.array[i] ¬ NEW[Refs2d[size.y]];
FOR j: NAT IN [0..size.y) DO
voxels.array[i][j] ¬ NEW[Refs[size.z]];
ENDLOOP;
ENDLOOP;
};
};

Store: PUBLIC PROC [voxels: Voxels, address: Nat3, voxel: REF] ~ {
ENABLE RuntimeError.BoundsFault => Error[$addressOutOfRange];
IF voxels = NIL THEN Error[$noVoxels];
IF voxels.type = hash
THEN [] ¬ RefTab.Store[voxels.hashTable, NEW[Nat3 ¬ address], voxel]
ELSE voxels.array[address.x][address.y][address.z] ¬ voxel;
};

RetrieveFromVoxels: PUBLIC PROC [voxels: Voxels, address: Nat3] RETURNS [voxel: REF] ~ {
ENABLE RuntimeError.BoundsFault => Error[$addressOutOfRange];
IF voxels = NIL THEN Error[$noVoxels];
IF voxels.type = hash
THEN {
refNat3: REF Nat3 ¬ ObtainRefNat3[];
refNat3­ ¬ address;
voxel ¬ RefTab.Fetch[voxels.hashTable, refNat3].val;
ReleaseRefNat3[refNat3];
}
ELSE RETURN[voxels.array[address.x][address.y][address.z]];
};

Equal: PROC [key1, key2: REF ANY] RETURNS [BOOL] ~ {
RETURN [NARROW[key1, REF Nat3]­ = NARROW[key2, REF Nat3]­];
};

Hash: PROC [key: REF ANY] RETURNS [c: CARDINAL] ~ {
a: REF Nat3 ~ NARROW[key];
c ¬ Basics.BITSHIFT[Basics.BITSHIFT[Basics.BITAND[a.x, 15], 4]+Basics.BITAND[a.y, 15], 4]+Basics.BITAND[a.z, 15];
};
CreateGrid: PUBLIC PROC [
values: Reals3d,					-- values at the grid points
box: Box ¬ [],					-- domain of the grid, presuming equal spacing
x, y, z: RealSequence ¬ NIL,	-- domain of the grid, presuming unequal spacing
clientData: REF ¬ NIL]
RETURNS [grid: Grid]
~ {
grid ¬ NEW[GridRep ¬ [x: x, y: y, z: z, values: values, clientData: clientData]];
IF box = []
THEN {
BadSeq: PROC [r: RealSequence] RETURNS [BOOL ¬ FALSE] ~ {
IF r = NIL OR r.length = 0 THEN RETURN[TRUE];
FOR n: NAT IN [1..r.length) DO IF r[n] <= r[n-1] THEN RETURN[TRUE]; ENDLOOP;
};
IF BadSeq[x] OR BadSeq[y] OR BadSeq[z] THEN Error[$badIndices];
grid.evenlySpaced ¬ FALSE;
IF values.length # x.length OR
values[0].length # y.length OR
values[0][0].length # z.length THEN Error[$indexMismatch];
grid.res ¬ [x.length, y.length, z.length];
grid.box ¬ MinMaxOfReals[x, y, z];
}
ELSE {
IF x # NIL OR y # NIL OR z # NIL THEN Error[$badIndices];
grid.evenlySpaced ¬ TRUE;
grid.res ¬ [values.length, values[0].length, values[0][0].length];
SetGridBox[grid, box];
};
};

MinMaxOfReals: PROC [x, y, z: RealSequence] RETURNS [Box] ~ {
RETURN[[[x[0], y[0], z[0]], [x[x.length-1], y[y.length-1], z[z.length-1]]]];
};

SetGridBox: PUBLIC PROC [grid: Grid, box: Box] ~ {
IF grid.evenlySpaced
THEN {
grid.scale ¬ [
REAL[grid.res.x]/(box.max.x-box.min.x),
REAL[grid.res.y]/(box.max.y-box.min.y),
REAL[grid.res.z]/(box.max.z-box.min.z)]
}
ELSE {
SetReals: PROC [reals: RealSequence, newMin: REAL] ~ {
oldMin: REAL ¬ reals[0];
FOR n: NAT IN [0..reals.length) DO
reals[n] ¬ (reals[n]-oldMin)*scale+newMin;
ENDLOOP;
};
nowDelta: Triple ¬ G3dVector.Sub[grid.box.max, grid.box.min];
newDelta: Triple ¬ G3dVector.Sub[box.max, box.min];
scales: Triple ¬ G3dVector.DivVectors[newDelta, nowDelta];
scale: REAL ¬ MIN[scales.x, scales.y, scales.z];
center: Triple ¬ G3dVector.Midpoint[box.min, box.max];
SetReals[grid.x, center.x-0.5*scale*nowDelta.x];
SetReals[grid.y, center.y-0.5*scale*nowDelta.y];
SetReals[grid.z, center.z-0.5*scale*nowDelta.z];
grid.box ¬ MinMaxOfReals[grid.x, grid.y, grid.z];
};
};

AllocateGridValues: PUBLIC PROC [res: Nat3] RETURNS [reals3d: Reals3d] ~ {
reals3d ¬ NEW[Reals3dRep[res.x]];
reals3d.length ¬ res.x;
FOR i: NAT IN [0..res.x) DO
reals2d: Reals2d ¬ reals3d[i] ¬ NEW[Reals2dRep[res.y]];
reals2d.length ¬ res.y;
FOR j: NAT IN [0..res.y) DO
reals: RealSequence ¬ reals2d[j] ¬ NEW[RealSequenceRep[res.z]];
reals.length ¬ res.z;
ENDLOOP;
ENDLOOP;
};

PointFromIJK: PUBLIC PROC [grid: Grid, i, j, k: NAT] RETURNS [p: Triple] ~ {
p ¬ IF grid.evenlySpaced
THEN [
grid.box.min.x+REAL[i]*(grid.box.max.x-grid.box.min.x)/REAL[grid.res.x-1],
grid.box.min.y+REAL[j]*(grid.box.max.y-grid.box.min.y)/REAL[grid.res.y-1],
grid.box.min.z+REAL[k]*(grid.box.max.z-grid.box.min.z)/REAL[grid.res.z-1]]
ELSE [grid.x[i], grid.y[j], grid.z[k]];
};
RetrieveFromGrid: PUBLIC PROC [grid: Grid, q: Triple] RETURNS [v: REAL] ~ {
InBox: PROC [p: Triple] RETURNS [BOOL] ~ {
RETURN[p.x IN [box.min.x..box.max.x] AND
       p.y IN [box.min.y..box.max.y] AND
       p.z IN [box.min.z..box.max.z]];
};
InlineCompute: PROC [ix0, iy0, iz0, ix1, iy1, iz1: NAT, xF, yF, zF: REAL] RETURNS [REAL]
~ INLINE {
x0: Reals2d ¬ grid.values[ix0];
x1: Reals2d ¬ grid.values[ix1];
v000: REAL ¬ x0[iy0][iz0];
v001: REAL ¬ x0[iy0][iz1];
v010: REAL ¬ x0[iy1][iz0];
v011: REAL ¬ x0[iy1][iz1];
v100: REAL ¬ x1[iy0][iz0];
v101: REAL ¬ x1[iy0][iz1];
v110: REAL ¬ x1[iy1][iz0];
v111: REAL ¬ x1[iy1][iz1];
v00: REAL ¬ v000+zF*(v001-v000);			-- x0y0
v01: REAL ¬ v010+zF*(v011-v010);			-- x0y1
v10: REAL ¬ v100+zF*(v101-v100);			-- x1y0
v11: REAL ¬ v110+zF*(v111-v110);			-- x1y1
v0: REAL ¬ v00+yF*(v01-v00);				-- x0
v1: REAL ¬ v10+yF*(v11-v10);				-- x1
RETURN[v0+xF*(v1-v0)];						-- tri-linearly interpolated value
};
box: Box ¬ grid.box;
IF NOT InBox[q] THEN Error[$outOfRange];
IF grid.evenlySpaced
THEN {
x: REAL ¬ (q.x-grid.box.min.x)*grid.scale.x;
y: REAL ¬ (q.y-grid.box.min.y)*grid.scale.y;
z: REAL ¬ (q.z-grid.box.min.z)*grid.scale.z;
ix0: NAT ¬ Real.Floor[x];
iy0: NAT ¬ Real.Floor[y];
iz0: NAT ¬ Real.Floor[z];
v ¬ InlineCompute[ix0, iy0, iz0, ix0+1, iy0+1, iz0+1, x-ix0, y-iy0, z-iz0];
}
ELSE {
GetIndex0: PROC [reals: RealSequence, r: REAL] RETURNS [NAT] ~ {
FOR n: NAT IN [1..reals.length) DO
IF reals[n] > r THEN RETURN[n-1];
REPEAT
FINISHED => RETURN[reals.length-2];	-- r at max, so i1 = length-1, i0 = length-2
ENDLOOP;
};
ix0: NAT ¬ GetIndex0[grid.x, q.x];
iy0: NAT ¬ GetIndex0[grid.y, q.y];
iz0: NAT ¬ GetIndex0[grid.z, q.z];
ix1: NAT ¬ ix0+1;
iy1: NAT ¬ iy0+1;
iz1: NAT ¬ iz0+1;
x0Grid: REAL ¬ grid.x[ix0];
y0Grid: REAL ¬ grid.y[iy0];
z0Grid: REAL ¬ grid.z[iz0];
xF: REAL ¬ (q.x-x0Grid)/(grid.x[ix1]-x0Grid);
yF: REAL ¬ (q.y-y0Grid)/(grid.y[iy1]-y0Grid);
zF: REAL ¬ (q.z-z0Grid)/(grid.z[iz1]-z0Grid);
v ¬ InlineCompute[ix0, iy0, iz0, ix1, iy1, iz1, xF, yF, zF];
};
};
ReadGridFromFile: PUBLIC PROC [fileName: ROPE] RETURNS [g: Grid] ~ {
s: IO.STREAM ¬ FS.StreamOpen[FileNames.ResolveRelativePath[fileName]
! FS.Error => GOTO NoOpen];
IF s # NIL THEN {
ReadLocations: PROC [keyWord: ROPE] RETURNS [reals: RealSequence ¬ NIL] ~ {
line ¬ G3dIO.FindKeyWord[s, keyWord, FALSE, 20, line];
line ¬ G3dIO.GetLine[s, line];
reals ¬ NEW[RealSequenceRep[G3dIO.NWordsInRope[line.rope]]];
FOR n: NAT IN [0..reals.length ¬ reals.maxLength) DO
reals[n] ¬ G3dIO.GetReal[line];
ENDLOOP;
};
line: G3dIO.Line ¬ G3dIO.FindKeyWord[s, "NumberOfIndices",, 20, G3dIO.ObtainLine[]];
resi: NAT ¬ IO.GetCard[s];
resj: NAT ¬ IO.GetCard[s];
resk: NAT ¬ IO.GetCard[s];
x: RealSequence ¬ ReadLocations["XLocations" ! G3dIO.Error => GOTO BadLocation];
y: RealSequence ¬ ReadLocations["YLocations" ! G3dIO.Error => GOTO BadLocation];
z: RealSequence ¬ ReadLocations["ZLocations" ! G3dIO.Error => GOTO BadLocation];
values: Reals3d ¬ AllocateGridValues[[resi, resj, resk]];
line ¬ G3dIO.FindKeyWord[s, "Values", FALSE, 20, line];
G3dIO.ReleaseLine[line];
FOR i: NAT IN [0..resi) DO
reals2d: Reals2d ¬ values[i];
FOR j: NAT IN [0..resj) DO
reals: RealSequence ¬ reals2d[j];
CedarProcess.CheckAbort[];
FOR k: NAT IN [0..resk) DO
reals[k] ¬ IO.GetReal[s];
ENDLOOP;
ENDLOOP;
ENDLOOP;
IO.Close[s];
g ¬ CreateGrid[values, [], x, y, z];
};
EXITS
NoOpen => Error[$NoSuchFile];
BadLocation => Error[$UnspecifiedDomain];
};

WriteGridToFile: PUBLIC PROC [grid: Grid, fileName: ROPE, comment: ROPE] ~ {
WriteRealSequence: PROC [keyword: ROPE, reals: RealSequence] ~ {
IO.PutF[s, "%l%g~%l\n", IO.rope["b"], IO.rope[keyword], IO.rope["B"]];
FOR n: NAT IN [0..reals.length) DO
IO.PutF1[s, "%g\t", IO.real[reals[n]]];
ENDLOOP;
IO.PutRope[s, "\n"];
};
s: IO.STREAM ¬ FS.StreamOpen[fileName, $create];
IO.PutRope[s, Rope.Cat["-- ", comment, "\n"]];
IO.PutFL[s, "%lNumberOfIndices~%l\n%g\t%g\t%g\n",
LIST[IO.rope["b"], IO.int[grid.res.x], IO.int[grid.res.y], IO.int[grid.res.z], IO.rope["B"]]];
WriteRealSequence["XLocations", grid.x];
WriteRealSequence["YLocations", grid.y];
WriteRealSequence["ZLocations", grid.z];
IO.PutF[s, "%lValues~%l\n", IO.rope["b"], IO.rope["B"]];
FOR i: NAT IN [0..grid.res.x) DO
reals2d: Reals2d ¬ grid.values[i];
FOR j: NAT IN [0..grid.res.y) DO
reals: RealSequence ¬ reals2d[j];
FOR k: NAT IN [0..grid.res.z) DO
IO.PutF1[s, "%g\t", IO.real[reals[k]]];
ENDLOOP;
IO.PutRope[s, "\n"];
ENDLOOP;
ENDLOOP;
};

GridFromVoxels: PUBLIC PROC [voxels: Voxels, realProc: RealProc] RETURNS [g: Grid] ~ {
g ¬ NEW[GridRep ¬ [evenlySpaced: TRUE]];
};
nScratchRefNat3s:	NAT ~ 10;
scratchRefNat3s:	ARRAY [0..nScratchRefNat3s) OF REF Nat3 ¬ ALL[NIL];

ObtainRefNat3: PROC RETURNS [REF Nat3] ~ {
FOR i: NAT IN [0..nScratchRefNat3s) DO
refNat3: REF Nat3 ~ scratchRefNat3s[i];
IF refNat3 # NIL THEN {
scratchRefNat3s[i] ¬ NIL;
RETURN[refNat3];
};
ENDLOOP;
RETURN[NEW[Nat3]];
};

ReleaseRefNat3: PROC [refNat3: REF Nat3] ~ {
FOR i: NAT IN [0..nScratchRefNat3s) DO
IF scratchRefNat3s[i] = NIL THEN {
scratchRefNat3s[i] ¬ refNat3;
RETURN;
};
ENDLOOP;
};

END.

j
G3dVoxelImpl.mesa
Copyright Σ 1988, 1992 by Xerox Corporation.  All rights reserved.
Bloomenthal, July 15, 1992 11:18 pm PDT
Errors
Storage and Retrieval of Voxel Data
Wyvill's formula: (rem(i,l)*m+rem(j,m))*n+rem(k,n) where l=m=n=16
RETURN[((a.i MOD 16)*16+(a.j MOD 16))*16+(a.k MOD 16);
Storage and Retrieval of Function Values Within A Grid


Support
Κ
Ί•NewlineDelimiter

–
"cedarcode" style™™Jšœ	Οeœ6™BJ™'J˜JšΟk	œ{˜„
J˜—šΠblœžœ
ž˜Jšžœ\˜cJšžœ	˜J˜—šœž˜J˜Jšœžœ˜ Jšœ
žœ˜#Jšœžœ˜-Jšœžœ˜2Jšœžœ˜Jšœ	žœ˜ Jšœžœ˜%Jšœžœ˜'Jšœ
žœ˜%Jšœ
žœ˜%Jšœžœ˜*Jšœžœ˜*Jšœž	œ˜ Jšœžœ˜#Jšœžœ˜#Jšœžœ˜#Jšœžœ˜(Jšžœžœžœ
˜—
headš
Οl™Jš
Οnœ
žœ
žœ
žœžœ
˜+—š
 #™#š‘œžœ
žœ
žœžœ žœžœ
˜YJšžœ˜J˜Jš	œ	žœ žœžœžœ˜Yšžœ˜	Jšžœ;˜?šžœ˜Jšžœ2˜8Jšœžœ˜#šžœžœ
žœ
ž˜Jšœžœ˜&šžœžœ
žœ
ž˜Jšœžœ˜'Jšžœ
˜—Jšžœ
˜——J˜—J˜J˜—š‘œžœ
žœ(žœ˜BJšžœ7˜=Jšžœ
žœ
žœ˜&šžœ˜Jšžœ%žœ˜DJšžœ7˜;—J˜J˜—š
‘œžœ
žœ!žœ	žœ˜XJšžœ7˜=Jšžœ
žœ
žœ˜&šžœ˜šžœ˜Jšœ	žœ˜$J˜J˜4J˜J˜
—Jšžœ
žœ0˜;—J˜J˜—š‘œžœžœ
žœžœžœ˜4Jš
žœžœžœ
žœžœ	˜;J˜J˜—š‘œžœžœ
žœžœžœ˜3Jšœžœžœ˜J™AJšžœžœ
žœžœ™6Jšœžœžœžœžœžœ
˜qJ˜——š
 6™6š‘
œžœ˜JšœΟc˜1Jšœ’.˜AJšœžœ’0˜MJšœžœžœ
˜Jšžœ
˜J˜JšœžœG˜Qšžœ	˜šžœ˜š
‘œžœžœžœžœ˜9Jšžœžœ
žœžœ
žœ
žœ˜-Jšžœžœ
žœžœžœ
žœ
žœžœ
˜LJ˜—Jšžœžœžœžœ˜?Jšœžœ
˜šžœž˜Jšœž˜Jšœžœ˜:—J˜*J˜"J˜
—šžœ˜Jšžœžœ
žœžœ
žœžœ
žœ˜9Jšœžœ
˜J˜BJ˜J˜——J˜J˜—š‘
œžœžœ
˜=JšžœF˜LJ˜J˜—š‘
œžœ
žœ˜2šžœ˜šžœ˜˜Jšžœ#˜'Jšžœ#˜'Jšžœ#˜'—J˜
—šžœ˜š‘œžœžœ˜6Jšœžœ˜šžœžœ
žœž˜"J˜*Jšžœ
˜—J˜—J˜=J˜3J˜:Jšœžœžœ˜0J˜6J˜0J˜0J˜0J˜1J˜——J˜J˜—š‘œžœ
žœ
žœ˜JJšœ
žœ˜!J˜šžœžœ
žœž˜Jšœžœ˜7J˜šžœžœ
žœž˜Jšœ#žœ˜?J˜Jšžœ
˜—Jšžœ
˜—J˜J˜—š
‘œžœ
žœžœžœ˜Lšœžœ˜šžœ˜Jšœžœ$žœ˜JJšœžœ$žœ˜JJšœžœ$žœ˜J—Jšžœ#˜'—J˜J™—š
‘œžœ
žœžœžœ˜Kš‘œžœ
žœžœ˜*šžœžœž˜(Jšœžœž˜(Jšœžœ˜&—J˜—š‘
œžœ žœžœ
žœžœ
˜XJšœžœ˜
J˜J˜Jšœžœ˜Jšœžœ˜Jšœžœ˜Jšœžœ˜Jšœžœ˜Jšœžœ˜Jšœžœ˜Jšœžœ˜Jšœžœ’˜*Jšœžœ’˜*Jšœžœ’˜*Jšœžœ’˜*Jšœžœ’˜%Jšœžœ’˜%Jšžœ’"˜>J˜—J˜Jšžœ
žœ	žœ˜(šžœ˜šžœ˜Jšœžœ%˜,Jšœžœ%˜,Jšœžœ%˜,Jšœžœ˜Jšœžœ˜Jšœžœ˜J˜KJ˜
—šžœ˜š
‘	œžœžœžœžœ˜@šžœžœ
žœž˜"Jšžœžœ
žœ˜!Jš
ž˜Jšžœžœ’,˜PJšžœ
˜—J˜—Jšœžœ˜"Jšœžœ˜"Jšœžœ˜"Jšœžœ	˜Jšœžœ	˜Jšœžœ	˜Jšœžœ˜Jšœžœ˜Jšœžœ˜Jšœžœ%˜-Jšœžœ%˜-Jšœžœ%˜-J˜<J˜——J˜J™—š
‘œžœ
žœžœžœ˜Dšœžœ
žœžœ3˜DJšœžœ
žœ	˜—šžœžœ
žœ˜š
‘
œžœžœžœžœ˜KJšœ%žœ˜6J˜Jšœžœ1˜<šžœžœ
žœ%ž˜4J˜Jšžœ
˜—J˜—J˜TJšœžœžœ˜Jšœžœžœ˜Jšœžœžœ˜Jšœ>žœ˜PJšœ>žœ˜PJšœ>žœ˜PJ˜9Jšœ&žœ˜7J˜šžœžœ
žœž˜J˜šžœžœ
žœž˜J˜!J˜šžœžœ
žœž˜Jšœžœ˜Jšžœ
˜—Jšžœ
˜—Jšžœ
˜—Jšžœ
˜J˜$J˜—š
ž˜Jš
œ˜J˜)—J˜J˜—š
‘œžœ
žœžœžœ˜Lš‘œžœžœ˜@Jšžœžœžœžœ˜Fšžœžœ
žœž˜"Jšžœžœ˜'Jšžœ
˜—Jšžœ˜J˜—Jšœžœ
žœžœ˜0Jšžœ,˜.šžœ/˜1Jšžœ
žœžœžœžœžœ
˜^—J˜(J˜(J˜(Jšžœžœžœ˜8šžœžœ
žœž˜ J˜"šžœžœ
žœž˜ J˜!šžœžœ
žœž˜ Jšžœžœ˜'Jšžœ
˜—Jšžœ˜Jšžœ
˜—Jšžœ
˜—J˜J˜—š‘œžœ
žœ&žœ˜VJšœžœžœ˜(J˜——š
 ™Jšœžœ˜šœžœžœ
žœžœ
žœ˜DJ˜—š‘
œžœ
žœžœ
˜*šžœžœ
žœž˜&Jšœ	žœ˜'šžœžœ
žœ˜Jšœžœ
˜Jšžœ
˜J˜—Jšžœ
˜—Jšžœ
žœ˜J˜J˜—š‘œžœžœ
˜,šžœžœ
žœž˜&šžœžœ
žœ˜"J˜Jšžœ
˜J˜—Jšžœ
˜—J˜—J˜—Jšžœ
˜—…—(θ8