G3dVoxelImpl.mesa
Copyright Ó 1988, 1992 by Xerox Corporation. All rights reserved.
Bloomenthal, July 15, 1992 11:18 pm PDT
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;
Errors
Error: PUBLIC SIGNAL [reason: ATOM] = CODE;
Storage and Retrieval of Voxel Data
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];
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);
c ¬ Basics.BITSHIFT[Basics.BITSHIFT[Basics.BITAND[a.x, 15], 4]+Basics.BITAND[a.y, 15], 4]+Basics.BITAND[a.z, 15];
};
Storage and Retrieval of Function Values Within A Grid
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]];
};
Support
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.