<<>> <> <> <> DIRECTORY G2dBasic, G3dBasic, G3dOctree, G3dPolygon, G3dVector, ImplicitAdapt, ImplicitDefs, ImplicitPoints, ImplicitPolygons, IO, MessageWindow, RealFns, Rope; ImplicitAdaptImpl: CEDAR MONITOR IMPORTS G2dBasic, G3dOctree, G3dPolygon, G3dVector, ImplicitPoints, ImplicitPolygons, IO, MessageWindow, RealFns EXPORTS ImplicitAdapt ~ BEGIN Triple: TYPE ~ G3dBasic.Triple; TripleSequence: TYPE ~ G3dBasic.TripleSequence; TripleSequenceRep: TYPE ~ G3dBasic.TripleSequenceRep; ActiveFaces: TYPE ~ G3dOctree.ActiveFaces; Cross: TYPE ~ G3dOctree.Cross; CrossedEdges: TYPE ~ G3dOctree.CrossedEdges; CrossPolygonProc: TYPE ~ G3dOctree.CrossPolygonProc; Cube: TYPE ~ G3dOctree.Cube; Corner: TYPE ~ G3dOctree.Corner; Corners: TYPE ~ G3dOctree.Corners; CubeProc: TYPE ~ G3dOctree.CubeProc; CubeStack: TYPE ~ G3dOctree.CubeStack; DirectionCorners: TYPE ~ G3dOctree.DirectionCorners; Edge: TYPE ~ G3dOctree.Edge; Face: TYPE ~ G3dOctree.Face; FourCorners: TYPE ~ G3dOctree.FourCorners; FourOctants: TYPE ~ G3dOctree.FourOctants; Octant: TYPE ~ G3dOctree.Octant; Octree: TYPE ~ G3dOctree.Octree; OctreeRep: TYPE ~ G3dOctree.OctreeRep; TwoCorners: TYPE ~ G3dOctree.TwoCorners; ThreeCubes: TYPE ~ G3dOctree.ThreeCubes; FourDirections: TYPE ~ G3dOctree.FourDirections; Direction: TYPE ~ G3dOctree.Direction; Neighborhood: TYPE ~ G3dOctree.Neighborhood; TwoFaces: TYPE ~ G3dOctree.TwoFaces; TwoOctants: TYPE ~ G3dOctree.TwoOctants; CrossSequence: TYPE ~ ImplicitDefs.CrossSequence; MacroPolygons: TYPE ~ ImplicitDefs.MacroPolygons; NormalProc: TYPE ~ ImplicitDefs.NormalProc; StatusProc: TYPE ~ ImplicitDefs.StatusProc; ValueProc: TYPE ~ ImplicitDefs.ValueProc; ROPE: TYPE ~ Rope.ROPE; <> << The Algorithm>> <> <> <> <> <> <> <> <> <> <> <> <> <> <<>> <> <> <<>> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <<>> SubdivideAbort: ERROR = CODE; Subdivide: PUBLIC PROC [ root: Cube, adaptLimit: NAT, flatLimit: REAL, valueProc: ValueProc, threshold: REAL ¬ 1.0, normalProc: NormalProc ¬ NIL, statusProc: StatusProc ¬ NIL, clientData: REF ANY ¬ NIL] RETURNS [msg: ROPE] ~ { SubdivideAndMore: PROC [parent: Cube] ~ { <> <> <> <> <> MaybeInheritSurfacePoint: PROC [edge: Edge] ~ { kid: Corner ~ directionCorners[G3dOctree.DirectionFromEdge[edge]]; par: Corner ~ G3dOctree.EdgeCorners[parent, edge].c1; IF par.inside = kid.inside THEN SELECT G3dOctree.AxisFromEdge[edge] FROM x => kid.lCross ¬ par.lCross; y => kid.bCross ¬ par.bCross; z => kid.nCross ¬ par.nCross; ENDCASE; }; MaybeSubdivideEdgeNeighbors: PROC [edge: Edge] ~ { twoCorners: TwoCorners ~ G3dOctree.EdgeCorners[parent, edge]; IF twoCorners.c0.inside = twoCorners.c1.inside THEN { kidCorner: Corner ~ directionCorners[G3dOctree.DirectionFromEdge[edge]]; IF twoCorners.c0.inside # kidCorner.inside THEN { cs: ThreeCubes ¬ G3dOctree.EdgeNeighbors[parent, edge]; IF cs.c0 # NIL AND cs.c0.terminal THEN AddToBottomOfStack[cs.c0]; IF cs.c1 # NIL AND cs.c1.terminal THEN AddToBottomOfStack[cs.c1]; IF cs.c2 # NIL AND cs.c2.terminal THEN AddToBottomOfStack[cs.c2]; }; }; }; MaybeSubdivideFaceNeighbor: PROC [face: Face] ~ { d: FourDirections ~ SELECT face FROM l => [lb, lt, ln, lf], r => [rb, rt, rn, rf], b => [lb, rb, bn, bf], t => [lt, rt, tn, tf], n => [ln, rn, bn, tn], ENDCASE => [lf, rf, bf, tf]; bool: BOOL ~ directionCorners[d.d0].inside; IF bool = directionCorners[d.d1].inside AND bool = directionCorners[d.d2].inside AND bool = directionCorners[d.g3d].inside AND bool # directionCorners[G3dOctree.DirectionFromFace[face]].inside THEN { nbor: Cube ~ neighbors[G3dOctree.DirectionFromFace[face]]; IF nbor # NIL AND nbor.terminal THEN AddToBottomOfStack[nbor]; }; }; RestrictOctree: PROC [parent: Cube] ~ { FOR direction: Direction IN Direction DO IF direction # c AND direction # none AND grandNeighbors[direction] # NIL AND grandNeighbors[direction].terminal THEN AddToBottomOfStack[grandNeighbors[direction]]; ENDLOOP; }; directionCorners: DirectionCorners; neighbors, grandNeighbors: Neighborhood; neighbors ¬ G3dOctree.Neighbors[parent]; grandNeighbors ¬ G3dOctree.Neighbors[parent.parent]; G3dOctree.Subdivide[parent]; ImplicitPoints.SetKidValues[parent, valueProc, threshold, clientData]; directionCorners ¬ G3dOctree.GetDirectionCorners[parent]; FOR edge: Edge IN Edge DO MaybeInheritSurfacePoint[edge]; ENDLOOP; IF parent.parent # NIL THEN RestrictOctree[parent]; -- restricted octree FOR edge: Edge IN Edge DO -- edge anomaly MaybeSubdivideEdgeNeighbors[edge]; ENDLOOP; FOR face: Face IN Face DO -- face anomaly MaybeSubdivideFaceNeighbor[face]; ENDLOOP; FOR octant: Octant IN Octant DO IF ShouldSubdivide[parent.kids[octant]] THEN AddToBottomOfStack[parent.kids[octant]]; ENDLOOP; }; ShouldSubdivide: PROC [cube: Cube] RETURNS [should: BOOL] ~ { IF cube = NIL OR NOT cube.terminal OR cube.level >= levelLimit THEN RETURN[FALSE] ELSE { mp: MacroPolygons ¬ MacroPolygonize[cube, valueProc, threshold, clientData]; polygon: CrossSequence ¬ mp.polygon; should ¬ SELECT mp.nPolygons FROM none => FALSE, one => TooRound[polygon, cosFlatLimit, cosPerpFlatLimit, valueProc, threshold, normalProc, cube, clientData], ENDCASE => TRUE; }; }; CheckStatusProc: PROC [cube: Cube, atom: ATOM] ~ { IF statusProc = NIL THEN RETURN; cube.refAny ¬ atom; IF statusProc[cube] = $Abort THEN ERROR SubdivideAbort; }; TerminalProc: CubeProc ~ { terminalLevel ¬ cube.level; levelLimit ¬ terminalLevel+adaptLimit; RETURN[FALSE]; }; MaybeSubdivide: CubeProc ~ { nTerminal ¬ nTerminal+1; CheckStatusProc[cube, $Dot]; IF ShouldSubdivide[cube] THEN AddToBottomOfStack[cube]; }; AddToBottomOfStack: PROC [cube: Cube] ~ { IF cube = NIL OR NOT cube.terminal OR cube.level >= levelLimit THEN RETURN; G3dOctree.WriteBottomOfCubeStack[cube, cubeStack]; CheckStatusProc[cube, $Cross]; }; ProcessTopOfStack: PROC ~ { cube: Cube ¬ G3dOctree.ReadTopOfCubeStack[cubeStack]; IF cube = NIL OR NOT cube.terminal OR cube.level >= levelLimit THEN RETURN; counts[cube.level] ¬ counts[cube.level]+1; SubdivideAndMore[cube]; CheckStatusProc[cube, $X]; }; terminalLevel, levelLimit, nTerminal: NAT ¬ 0; counts: ARRAY [0..100) OF NAT ¬ ALL[0]; cosPerpFlatLimit: REAL ¬ RealFns.CosDeg[90.0-flatLimit]; cosFlatLimit: REAL ¬ RealFns.CosDeg[flatLimit]; cubeStack: CubeStack ¬ G3dOctree.NewCubeStack[20000]; ImplicitPoints.SetTerminalCubeValues[root, valueProc, threshold, clientData]; G3dOctree.ApplyToTerminal[root, TerminalProc]; G3dOctree.ApplyToLevel[root, MaybeSubdivide, terminalLevel]; WHILE NOT G3dOctree.CubeStackEmpty[cubeStack] DO ProcessTopOfStack[! SubdivideAbort => EXIT]; ENDLOOP; msg ¬ IO.PutFR["%g/%g", IO.int[counts[terminalLevel]], IO.int[nTerminal]]; FOR n: NAT IN (terminalLevel..100) WHILE counts[n] # 0 DO msg ¬ IO.PutFR["%g,%g/%g", IO.rope[msg], IO.int[counts[n]], IO.int[8*counts[n-1]]]; ENDLOOP; msg ¬ IO.PutFR1["%g adapted", IO.rope[msg]]; }; <> << The Algorithm>> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> < 1 polygons => subdivision of cube and its neighbors, why?>> <> <> <> <> <> <> <> <> <> <> <> <<>> <> <<19 new corners (1 center of parent, 12 midpoints of edges, and 6 face centers).>> <<>> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <> <<>> TrackAbort: ERROR = CODE; Track: PUBLIC PROC [ seedSize: REAL, surfacePoint: Triple, valueProc: ValueProc, threshold: REAL ¬ 1.0, normalProc: NormalProc ¬ NIL, statusProc: StatusProc ¬ NIL, clientData: REF ANY ¬ NIL, adaptLimit: NAT ¬ 0, flatLimit: REAL ¬ 30.0, tolerance: REAL ¬ 0.0001] RETURNS [octree: Octree] ~ { GetActiveFaces: PROC [cube: Cube] RETURNS [activeFaces: ActiveFaces] ~ { FOR e: Edge IN Edge DO corners: TwoCorners ¬ G3dOctree.EdgeCorners[cube, e]; IF corners.c0.inside # corners.c1.inside AND NOT corners.c0.outOfRange AND NOT corners.c1.outOfRange THEN { twoFaces: TwoFaces ~ G3dOctree.EdgeFaces[e]; activeFaces[twoFaces.f0] ¬ activeFaces[twoFaces.f1] ¬ TRUE; }; ENDLOOP; }; ShouldSubdivideTest: PROC [cube: Cube] RETURNS [should: BOOL ¬ FALSE] ~ { crossPolygonProc: CrossPolygonProc ~ { InnerTest: PROC RETURNS [shouldSubdivide: BOOL ¬ FALSE] ~ { SELECT TRUE FROM nPolygon # 0 => RETURN[TRUE]; -- more than one polygon in cube polygon.length >= 3 => { <> center, normal, p1: Triple; FOR n: NAT IN [0..polygon.length) DO points[n] ¬ polygon[n].point; ENDLOOP; points.length ¬ polygon.length; center ¬ G3dPolygon.PolygonCenter[points]; normal ¬ ImplicitPoints.GetPointNormal[ center, valueProc, valueProc[center, clientData], delta, clientData]; FOR n: NAT IN [0..polygon.length) DO IF G3dVector.Dot[polygon[n].normal, normal] < cosFlatLimit THEN RETURN[TRUE]; ENDLOOP; <> p1 ¬ polygon[polygon.length-1].point; FOR n: NAT IN [0..polygon.length) DO p0: Triple ¬ p1; v: Triple ¬ G3dVector.Sub[p1 ¬ polygon[n].point, p0]; v ¬ G3dVector.Unit[v]; IF G3dVector.Dot[normal, v] > cosPerpFlatLimit THEN RETURN[TRUE]; ENDLOOP; }; ENDCASE; }; IF (should ¬ InnerTest[]) THEN RETURN[FALSE]; }; IF TRUE THEN RETURN[FALSE]; IF 0.5*cube.size < sizeLimit THEN RETURN[FALSE]; IF cube.size # cubeSize THEN delta ¬ 0.001*(cubeSize ¬ cube.size); [] ¬ ImplicitPoints.CheckAndGetCrossedEdges[ cube, valueProc, threshold, NProc, tolerance, clientData]; [] ¬ ImplicitPolygons.CrossPolygonizeCube[cube, crossPolygonProc, crossSequence, FALSE]; }; ShouldCoalesce: PROC [parent: Cube] RETURNS [should: BOOL] ~ { <> IF parent = NIL OR parent.terminal THEN RETURN[FALSE]; IF (should ¬ NOT ShouldSubdivideTest[parent]) THEN { <> FOR edge: Edge IN Edge DO twoOctants: TwoOctants ~ G3dOctree.EdgeOctants[edge]; state: BOOL ¬ parent.corners[twoOctants.o0].inside; IF parent.corners[twoOctants.o1].inside = state THEN { o: Octant ¬ twoOctants.o1; k: Cube ¬ parent.kids[twoOctants.o0]; IF k = NIL THEN {k ¬ parent.kids[twoOctants.o1]; o ¬ twoOctants.o0}; IF k # NIL AND k.corners[o].inside # state THEN RETURN[FALSE]; }; ENDLOOP; <> FOR face: Face IN Face DO fourOctants: FourOctants ¬ G3dOctree.FaceOctants[face]; state: BOOL ¬ parent.corners[fourOctants.o0].inside; IF parent.corners[fourOctants.o1].inside = state AND parent.corners[fourOctants.o2].inside = state AND parent.corners[fourOctants.o3].inside = state THEN { o: Octant ¬ fourOctants.o2; k: Cube ¬ parent.kids[fourOctants.o0]; IF k = NIL THEN {k ¬ parent.kids[fourOctants.o1]; o ¬ fourOctants.o3}; IF k = NIL THEN {k ¬ parent.kids[fourOctants.o2]; o ¬ fourOctants.o0}; IF k = NIL THEN {k ¬ parent.kids[fourOctants.o3]; o ¬ fourOctants.o1}; IF k # NIL AND k.corners[o].inside # state THEN RETURN[FALSE]; }; ENDLOOP; }; }; DefaultNormalProc: NormalProc ~ { RETURN[ImplicitPoints.GetPointNormal[point, valueProc, value, delta, clientData]]; }; AddToBottomOfStack: PROC [cube: Cube] ~ { IF (counter ¬ counter+1) > 10 THEN { counter ¬ 0; MessageWindow.Append[IO.PutFR1["\tstack size: %g", IO.int[cubeStack.size]], TRUE]; }; IF cube = NIL OR NOT cube.terminal OR cube.size < sizeLimit THEN RETURN; ImplicitPoints.SetCornerValues[cube, valueProc, threshold, clientData]; IF NOT ImplicitPoints.IsCrossed[cube] THEN RETURN; G3dOctree.WriteBottomOfCubeStack[cube, cubeStack]; cube.refAny ¬ $AddedToStack; IF statusProc # NIL AND statusProc[cube] = $Abort THEN ERROR TrackAbort; IF NOT G3dOctree.FullyPointed[G3dOctree.Root[cube]] THEN ERROR; }; InnerSubdivide: PROC [parent: Cube] ~ { <> <> <> MaybeInheritCross: PROC [edge: Edge] ~ { kid: Corner ~ directionCorners[G3dOctree.DirectionFromEdge[edge]]; par: Corner ~ G3dOctree.EdgeCorners[parent, edge].c1; IF par.inside = kid.inside THEN SELECT G3dOctree.AxisFromEdge[edge] FROM x => kid.lCross ¬ par.lCross; y => kid.bCross ¬ par.bCross; z => kid.nCross ¬ par.nCross; ENDCASE; }; MaybeSubdivideEdgeNeighbors: PROC [edge: Edge] ~ { twoCorners: TwoCorners ~ G3dOctree.EdgeCorners[parent, edge]; IF twoCorners.c0.inside = twoCorners.c1.inside THEN { kidCorner: Corner ~ directionCorners[G3dOctree.DirectionFromEdge[edge]]; IF twoCorners.c0.inside # kidCorner.inside THEN { ds: G3dOctree.ThreeDirections ¬ G3dOctree.EdgeDirections[edge]; c0: Cube ¬ neighbors[ds.d0]; c1: Cube ¬ neighbors[ds.d1]; c2: Cube ¬ neighbors[ds.d2]; IF c0 # NIL AND c0.terminal THEN InnerSubdivide[c0]; IF c1 # NIL AND c1.terminal THEN InnerSubdivide[c1]; IF c2 # NIL AND c2.terminal THEN InnerSubdivide[c2]; }; }; }; MaybeSubdivideFaceNeighbor: PROC [face: Face] ~ { d: FourDirections ¬ G3dOctree.FaceDirections[face]; bool: BOOL ~ directionCorners[d.d0].inside; IF bool = directionCorners[d.d1].inside AND bool = directionCorners[d.d2].inside AND bool = directionCorners[d.g3d].inside AND bool # directionCorners[G3dOctree.DirectionFromFace[face]].inside THEN { neighbor: Cube ~ neighbors[G3dOctree.DirectionFromFace[face]]; IF neighbor # NIL AND neighbor.terminal THEN InnerSubdivide[neighbor]; }; }; directionCorners: DirectionCorners; neighbors: Neighborhood ¬ G3dOctree.Neighbors[parent]; G3dOctree.Subdivide[parent]; ImplicitPoints.SetKidValues[parent, valueProc, threshold, clientData]; FOR o: Octant IN Octant DO k: Cube ¬ parent.kids[o]; IF ShouldSubdivideTest[k] THEN InnerSubdivide[k] ELSE AddToBottomOfStack[k]; ENDLOOP; directionCorners ¬ G3dOctree.GetDirectionCorners[parent]; FOR edge: Edge IN Edge DO MaybeInheritCross[edge]; ENDLOOP; FOR edge: Edge IN Edge DO -- test for edge anomaly MaybeSubdivideEdgeNeighbors[edge]; ENDLOOP; FOR face: Face IN Face DO -- test for face anomaly MaybeSubdivideFaceNeighbor[face]; ENDLOOP; }; ProcessTopOfStack: PROC [cubeStack: CubeStack] ~ { cube: Cube ¬ G3dOctree.ReadTopOfCubeStack[cubeStack]; IF (counter ¬ counter+1) > 10 THEN { counter ¬ 0; MessageWindow.Append[IO.PutFR1["\tstack size: %g", IO.int[cubeStack.size]], TRUE]; }; <> IF cube.terminal AND cube.refAny # $Coalesced THEN { activeFaces: ActiveFaces ¬ GetActiveFaces[cube]; FOR face: Face IN Face DO <> IF cube.refAny = $Coalesced THEN RETURN; IF activeFaces[face] THEN { neighbor: Cube ¬ G3dOctree.FaceNeighbor[cube, face]; IF neighbor = NIL THEN { new: Cube; <> <> octree.root ¬ G3dOctree.AddCube[cube, face]; new ¬ G3dOctree.FaceNeighbor[cube, face]; IF NOT new.terminal THEN ERROR; FOR p: Cube ¬ new.parent, p.parent WHILE p # NIL DO <> <> IF p.refAny = $Coalesced THEN { Coalesce[p]; RETURN; }; ENDLOOP; ConsiderCube[new]; } <> ELSE IF neighbor.refAny = NIL THEN ConsiderCube[neighbor]; }; ENDLOOP; }; <> }; Coalesce: PROC [parent: Cube] ~ { <> Inner: CubeProc ~ {cube.refAny ¬ $Coalesced}; G3dOctree.ApplyToKids[parent, Inner]; }; ConsiderCube: PROC [new: Cube] ~ { IF new.refAny = $Coalesced THEN ERROR; new.refAny ¬ $Considered; IF ShouldSubdivideTest[new] THEN InnerSubdivide[new] ELSE IF ShouldCoalesce[new.parent] THEN { Coalesce[new.parent]; <> new.parent.terminal ¬ TRUE; AddToBottomOfStack[new.parent]; } ELSE AddToBottomOfStack[new]; }; RemoveCoalescedCubes: PROC [root: Cube] ~ { IsCoalesced: PROC [cube: Cube] RETURNS [b: BOOL ¬ TRUE] ~ { Test: CubeProc ~ {IF cube.refAny # $Coalesced THEN ERROR}; IF (b ¬ cube.refAny = $Coalesced) THEN G3dOctree.ApplyToKids[cube, Test]; }; Inner: PROC [cube: Cube] ~ { cube.terminal ¬ TRUE; FOR o: Octant IN Octant DO k: Cube ¬ cube.kids[o]; IF k = NIL THEN LOOP; IF IsCoalesced[k] THEN cube.kids[o] ¬ NIL ELSE { cube.terminal ¬ FALSE; Inner[k]; }; ENDLOOP; }; Inner[root]; }; CheckForNeededSubdivisions: PROC [root: Cube] ~ { Inner: CubeProc ~ { state: BOOL; neighbors: Neighborhood ¬ G3dOctree.Neighbors[cube]; TestNeighborIfDivide: PROC [d: Direction] RETURNS [b: BOOL] ~ { n: Cube ¬ neighbors[d]; b ¬ n # NIL AND NOT n.terminal AND n.refAny # $Coalesced; IF b THEN { c: Corner ¬ G3dOctree.GetDirectionCorner[n, G3dOctree.OppositeDirection[d]]; b ¬ c # NIL AND state # c.inside; }; }; TestEdgeNeighborIfDivide: PROC [d0, d1: Direction] RETURNS [b: BOOL] ~ { n: Cube ¬ neighbors[d0]; b ¬ n # NIL AND NOT n.terminal AND n.refAny # $Coalesced; IF b THEN { d: Direction ¬ G3dOctree.AddDirection[d1, G3dOctree.OppositeDirection[d0]]; c: Corner ¬ G3dOctree.GetDirectionCorner[n, d]; b ¬ c # NIL AND state # c.inside; }; }; FOR edge: Edge IN Edge DO twoCorners: TwoCorners ~ G3dOctree.EdgeCorners[cube, edge]; state ¬ twoCorners.c0.inside; IF state = twoCorners.c1.inside THEN { ds: G3dOctree.ThreeDirections ¬ G3dOctree.EdgeDirections[edge]; IF TestEdgeNeighborIfDivide[ds.d0, ds.d1] THEN GOTO subdiv; IF TestEdgeNeighborIfDivide[ds.d1, ds.d0] THEN GOTO subdiv; IF TestNeighborIfDivide[ds.d2] THEN GOTO subdiv; }; ENDLOOP; FOR face: Face IN Face DO fourCorners: FourCorners ¬ G3dOctree.FaceCorners[cube, face]; state ¬ fourCorners.c0.inside; IF state = fourCorners.c1.inside AND state = fourCorners.c2.inside AND state = fourCorners.c3.inside AND TestNeighborIfDivide[G3dOctree.DirectionFromFace[face]] THEN GOTO subdiv; ENDLOOP; EXITS subdiv => { InnerSubdivide[cube]; FOR o: Octant IN Octant DO [] ¬ Inner[cube.kids[o]]; ENDLOOP; }; }; G3dOctree.ApplyToTerminal[root, Inner]; }; NProc: NormalProc ¬ IF normalProc # NIL THEN normalProc ELSE DefaultNormalProc; counter, nRecursiveTrackFaces, nRecurseFailures: NAT ¬ 0; cubeStack: CubeStack ¬ G3dOctree.NewCubeStack[20000]; cosFlatLimit: REAL ¬ RealFns.CosDeg[flatLimit]; cosPerpFlatLimit: REAL ¬ RealFns.CosDeg[90.0-flatLimit]; delta, cubeSize: REAL ¬ 0.0; crossSequence: CrossSequence ¬ G3dOctree.ObtainCrossSequence[]; points: TripleSequence ¬ NEW[TripleSequenceRep[12]]; sizeLimit: REAL ¬ seedSize; FOR n: NAT IN [0..adaptLimit) DO sizeLimit ¬ 0.5*sizeLimit; ENDLOOP; octree ¬ NEW[OctreeRep]; octree.root ¬ G3dOctree.NewCube[seedSize, surfacePoint]; ConsiderCube[octree.root]; WHILE NOT G3dOctree.CubeStackEmpty[cubeStack] DO ProcessTopOfStack[cubeStack ! TrackAbort => EXIT]; ENDLOOP; CheckForNeededSubdivisions[octree.root]; RemoveCoalescedCubes[octree.root]; IF NOT G3dOctree.FullyPointed[octree.root] THEN ERROR; IF NOT G3dOctree.CubeOk[octree.root] THEN ERROR; G3dOctree.ReleaseCrossSequence[crossSequence]; G3dOctree.SetOctreeFields[octree]; }; <> <> <> <> <> <> <> <> <> <> <<};>> <> <<};>> <> Roundness: PUBLIC PROC [ cube: Cube, valueProc: ValueProc, threshold: REAL ¬ 1.0, normalProc: NormalProc, clientData: REF ANY ¬ NIL] RETURNS [max: REAL ¬ 0.0] ~ { v, normal: Triple ¬ []; polygon: CrossSequence ¬ MacroPolygonize[cube, valueProc, threshold, clientData].polygon; IF polygon = NIL OR polygon.length = 0 THEN RETURN; normal ¬ G3dVector.Unit[NormalFromCrossSequence[polygon]]; <> IF polygon.length > 3 THEN { p1: Triple ¬ polygon[polygon.length-1].point; FOR i: NAT IN [0..polygon.length) DO p0: Triple ¬ p1; v ¬ G3dVector.Unit[G3dVector.Sub[p1 ¬ polygon[i].point, p0]]; max ¬ MAX[max, 90.0-G2dBasic.ArcCos[G3dVector.Dot[normal, v]]]; ENDLOOP; }; <> FOR n: NAT IN [0..polygon.length) DO cross: Cross ¬ polygon[n]; IF cross.normal = [] THEN { cross.normal ¬ IF normalProc # NIL THEN normalProc[cross.point, cross.value, clientData] ELSE ImplicitPoints.GetPointNormal[ cross.point, valueProc, cross.value, 0.01*cube.size, clientData]; }; max ¬ MAX[max, G2dBasic.ArcCos[G3dVector.Dot[normal, cross.normal]]]; ENDLOOP; G3dOctree.ReleaseCrossSequence[polygon]; }; NormalFromCrossSequence: PUBLIC PROC [polygon: CrossSequence] RETURNS [n: Triple ¬ []] ~ { p1: Triple ¬ polygon[polygon.length-1].point; FOR i: NAT IN[0..polygon.length) DO -- Newells' method p0: Triple ¬ p1; p1 ¬ polygon[i].point; n.x ¬ n.x+(p0.y-p1.y)*(p0.z+p1.z); n.y ¬ n.y+(p0.z-p1.z)*(p0.x+p1.x); n.z ¬ n.z+(p0.x-p1.x)*(p0.y+p1.y); ENDLOOP; }; LinearDeviation: PUBLIC PROC [cube: Cube] RETURNS [error: REAL ¬ 0.0] ~ { corners: Corners ~ cube.corners; F: PROC [i, j, k: NAT] RETURNS [REAL] ~ { o: Octant; d: Direction ¬ c; d ¬ G3dOctree.AddDirection[d, IF i = 0 THEN l ELSE r]; d ¬ G3dOctree.AddDirection[d, IF j = 0 THEN b ELSE t]; d ¬ G3dOctree.AddDirection[d, IF k = 0 THEN n ELSE f]; o ¬ G3dOctree.OctantFromDirection[d]; IF NOT corners[o].valueSet THEN ERROR; RETURN[corners[o].value]; }; f000: REAL ~ corners[lbn].value; f001: REAL ~ corners[lbf].value; f010: REAL ~ corners[ltn].value; f011: REAL ~ corners[ltf].value; f100: REAL ~ corners[rbn].value; f101: REAL ~ corners[rbf].value; f110: REAL ~ corners[rtn].value; f111: REAL ~ corners[rtf].value; a: REAL ¬ 0.25*((f100+f101+f110+f111)-(f000+f001+f010+f011)); b: REAL ~ 0.25*((f010+f011+f110+f111)-(f000+f001+f100+f101)); c: REAL ~ 0.25*((f001+f011+f101+f111)-(f000+f100+f010+f110)); d: REAL ~ 0.5*f000+0.25*(f100+f010+f001-f111); FOR i: NAT IN [0..1] DO FOR j: NAT IN [0..1] DO FOR k: NAT IN [0..1] DO dif: REAL ~ a*i+b*j+c*k+d-F[i, j, k]; error ¬ error+dif*dif; ENDLOOP; ENDLOOP; ENDLOOP; error ¬ error/cube.size; }; TooRound: PUBLIC PROC [ polygon: CrossSequence, cosFlatLimit: REAL, cosPerpFlatLimit: REAL, valueProc: ValueProc, threshold: REAL ¬ 1.0, normalProc: NormalProc, cube: Cube, clientData: REF ANY] RETURNS [tooRound: BOOL ¬ FALSE] ~ { p0, p1, v, normal: Triple ¬ []; p1 ¬ polygon[polygon.length-1].point; FOR i: NAT IN[0..polygon.length) DO -- Newells' method p0 ¬ p1; p1 ¬ polygon[i].point; normal ¬ [ normal.x+(p0.y-p1.y)*(p0.z+p1.z), normal.y+(p0.z-p1.z)*(p0.x+p1.x), normal.z+(p0.x-p1.x)*(p0.y+p1.y)]; ENDLOOP; normal ¬ G3dVector.Unit[normal]; <> IF polygon.length > 3 THEN { p1 ¬ polygon[polygon.length-1].point; FOR i: NAT IN [0..polygon.length) DO p0 ¬ p1; v ¬ G3dVector.Unit[G3dVector.Sub[p1 ¬ polygon[i].point, p0]]; IF G3dVector.Dot[normal, v] > cosPerpFlatLimit THEN {tooRound ¬ TRUE; EXIT}; ENDLOOP; }; <> IF tooRound = FALSE THEN FOR n: NAT IN [0..polygon.length) DO cross: Cross ¬ polygon[n]; IF cross.normal = [] THEN { cross.normal ¬ IF normalProc # NIL THEN normalProc[cross.point, cross.value, clientData] ELSE ImplicitPoints.GetPointNormal[ cross.point, valueProc, cross.value+threshold, 0.01*cube.size, clientData]; }; IF G3dVector.Dot[normal, cross.normal] < cosFlatLimit THEN {tooRound ¬ TRUE; EXIT}; ENDLOOP; }; <> <> <<};>> <<>> MacroPolygonize: PUBLIC PROC [ cube: Cube, valueProc: ValueProc, threshold: REAL ¬ 1.0, clientData: REF ANY] RETURNS [macro: MacroPolygons] ~ { crossedEdges: CrossedEdges ¬ ImplicitPoints.CheckAndGetCrossedEdges[cube, valueProc, threshold,,, clientData]; FOR e: Edge IN Edge DO IF crossedEdges[e].cIn # NIL THEN GOTO found; REPEAT found => { edge, start: Edge ¬ e; done: ARRAY Edge OF BOOL ¬ ALL[FALSE]; crosses: CrossSequence ¬ G3dOctree.ObtainCrossSequence[]; DO edge ¬ ImplicitPoints.NextCWEdgeCrossed[edge, crossedEdges]; crosses[crosses.length] ¬ ImplicitPoints.CrossedPoint[crossedEdges[edge], valueProc, threshold,,,, clientData]; crosses.length ¬ crosses.length+1; done[edge] ¬ TRUE; IF edge = start THEN { macro.polygon ¬ crosses; macro.nPolygons ¬ one; EXIT; }; ENDLOOP; FOR e: Edge IN Edge DO IF crossedEdges[e].cIn # NIL AND NOT done[e] THEN { G3dOctree.ReleaseCrossSequence[crosses]; macro.nPolygons ¬ many; RETURN; }; ENDLOOP; }; FINISHED => { macro.nPolygons ¬ none; RETURN; }; ENDLOOP; }; END.