<<>> <> <> <> 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.