INTERFACE Triangulation; (* Triangulated tridimensional meshes for automatic topology visualization. This interface contain essentially procedures created by J. Stolfi and R. Marcone (see the copyright and authorship futher down), modified extensively by L. Lozada for the tridimensional case. *) IMPORT Octf, Random, R3, LR4, Wr; (* === ELEMENTS === *) TYPE CARD = CARDINAL; BOOL = BOOLEAN; VisitProc = Octf.VisitProc; SRBits = Octf.SRBits; Pair = Octf.Pair; Face = Octf.Face; Edge = Octf.Edge; Node = Octf.Node; (* A "Triangulation.Pair" is an "Octf.Pair" whose "facetedge" field is a "Triangulation.FacetEdge". *) TYPE FacetEdge <: PublicFacetEdge; PublicFacetEdge = Octf.FacetEdge OBJECT mark: BOOL := FALSE; (* Mark the facetedge, used by RefineT.*) vh: Vertex; (* Vertex medial associated to facetedge.*) old: CARD; METHODS init(): FacetEdge; (* initializes "self" as an isolated facetedge of triangulation *) END; TYPE Vertex <: PublicVertex; PublicVertex = Node OBJECT exists: BOOL := TRUE; (* FALSE for ghost vertices *) fixed: BOOL := FALSE; (* TRUE if position is fixed *) xmark: BOOL := FALSE; (* Used by MakePolyhedronTopology *) color: R3.T := R3.T{0.0,0.0,0.0};(* Color for painting *) transp: R3.T := R3.T{0.0,..}; (* Transp. coefficient for painting *) radius: REAL := 0.020; (* Vertex radius for drawing *) label: TEXT := "VV"; (* label vertex, useful in bar.sub *) END; Vertices = ARRAY OF Vertex; Polyhedron <: PublicPolyhedron; PublicPolyhedron = Node OBJECT exists: BOOL := TRUE; (* FALSE for ghost tetrahedrons *) ymark: BOOL := FALSE; (* Used by ????? *) (* vertex: REF ARRAY OF Node; (* Vertices defining the polyhedron *) *) color : R3.T := R3.T{1.0,..}; (* Color for painting *) transp: R3.T := R3.T{1.0,..}; (* Transp. coefficient for painting *) root: INTEGER := -1; (* save the "root" tetrahedron *) END; (* === ELEMENT CREATION === *) PROCEDURE MakeFacetEdge(): Pair; (* Creates a new unattached facetedge with distinct endpoints. *) PROCEDURE MakeVertex(): Vertex; (* Creates an unattached vertex record. *) PROCEDURE MakePolyhedron(): Polyhedron; (* Creates an unattached polyhedron record. *) (* === PAIR PROPERTIES === *) PROCEDURE Org(a: Pair): Node; (* The origin of pair "a". *) PROCEDURE SetOrg(a: Pair; n: Node); (* Set the origin of pair "a" to be "n". *) PROCEDURE SetAllOrgs(a: Pair; n: Node); (* Does "SetOrg(t,n)" for some pairs "t" with the same origin as "a". *) PROCEDURE Set(a: Pair; n: Node); (* Does "SetOrg(t,n)" for all adjacents pairs "t" reaches by chains Fnext Enexts and Clocks. *) PROCEDURE Pneg(a: Pair): Node; (* The polyhedron negative of "a". *) PROCEDURE SetPneg(a: Pair; n: Node); (* Set the polyhedron negative of "a" to "n". Same as: "SetOrg(Sdual(a),n)". *) PROCEDURE SetNextPneg(a: Pair; n: Node); (* Does "SetPneg(t,n)" for all pairs "t" with the same face component. *) PROCEDURE SetAllPneg(a: Pair; n: Node); (* Set the pairs facetdeges (12) belonging to the same negative polyhedron "Pneg" equal to "n". *) PROCEDURE Ppos(a: Pair): Node; (* The polyhedron positive of "a". *) PROCEDURE SetPpos(a: Pair; n: Node); (* Sets the polyhedron positive of "a" to "n". Same as: "SetOrg(Clock(Sdual(a),n) = SetOrg(Tors(a),n)". *) PROCEDURE SetAllPpos(a: Pair; n: Node); (* Does "SetPpos(t,n)" for all pairs "t" with same Ppos as "a". *) PROCEDURE SetNextPpos(a : Pair; n: Node); (* Set the pairs facetdeges (12) belong to same polyhedron positive "Ppos" equal to "n". *) PROCEDURE OrgV(a: Pair): Vertex; (* The origin of pair "a", narrowed to type "Vertex". Corresponding to one vertex of primal subdivision: "C". *) PROCEDURE DesV(a: Pair): Vertex; (* The destination of pair "a", narrowed to type "Vertex". Corresponding to one vertex of primal subdivision: "C". *) PROCEDURE PnegP(a : Pair): Polyhedron; (* The polyhedron negative of pair "a", narrowed to type "Polyhedron". Corresponding to one vertex of dual subdivision: "C'". *) PROCEDURE PposP(a : Pair): Polyhedron; (* The polyhedron positive of pair "a", narrowed to type "Polyhedron". Corresponding to one vertex of dual subdivision: "C'". *) PROCEDURE TetraNegVertices(a: Pair): ARRAY [0..3] OF Vertex; (* The vertices of the tetrahedron "PnegP(a)". *) PROCEDURE TetraPosVertices(a: Pair): ARRAY [0..3] OF Vertex; (* The vertices of the tetrahedron "PposP(a)". *) PROCEDURE TetraNegPosVertices(a: Pair): ARRAY [0..4] OF Node; (* The vertices of the tetrahedra "PnegP(a)" and "PposP(a)". *) PROCEDURE TetraFaces(a: Pair): ARRAY [0..3] OF Face; (* The four faces of the tetrahedron "PnegP(a)". *) PROCEDURE TetraEdges(a: Pair): ARRAY [0..5] OF Edge; (* The six edges of the tetrahedron "PnegP(a)". *) PROCEDURE FaceEdges(a: Pair): ARRAY [0..2] OF Edge; (* The three edges of the component face of the pair "a". *) PROCEDURE EdgeIsBorder(a: Pair): BOOL; (* TRUE if the edge of "a" belongs to the manifold's border (i.e. is on the boundary of a NIL polyhedron). *) PROCEDURE FaceIsBorder(a: Pair): BOOL; (* TRUE iff the face of "a" belongs to the manifold's border. Same as "Ppos(a) = NIL OR Pneg(a) = NIL". *) (* === CONSTRUCTION TOOLS === *) PROCEDURE MakeTetraTopo(nx, ny: CARD): ARRAY [0..7] OF Pair; (* Builds a topological tetrahedron subdivided radialy. The number of tetrahedra cells is "nx" by "ny". Returns the corners facetedges pairs (8). *) PROCEDURE EmphasizeTetrahedron(a, b: Pair; n: CARD); (* Emphasizes the original elements of a tetrahedron produced by the MakeTetraTopo(order,order) procedure. *) PROCEDURE Glue( a,b : Pair; n: CARD; setorg: BOOL := TRUE; ) : Pair; (* Glue two topological tetrahedra, by identification the "n" triangular facets on the topological boundary of tetrahedra. The traversal of the topological boundary will be realized by Quad- Edge's functions reducided to functions FacetEdge.The identification of triangular facets will be done by the "meld" procedure, removing the pairs facetedges on the chain "b" and updating the relations between vertices and polyhedra. If SetOrg = TRUE then the SetAllOrgs is executed. *) (* === GLOBAL PROCEDURES === *) PROCEDURE NumberVertices(a: Pair): CARD; (* Enumerates all (primal) vertices reachable from "a" by chains of "Onext", "Onext_1" and "Clock", and assigns them distinct serial numbers from 0 up. Returns the number of vertices found.Use the procedure travesal "EnumVertices" of library "libm3triang". *) TYPE Topology = RECORD NV: INTEGER; (* Number of vertex *) NE: CARD; (* Number of edges *) NF: CARD; (* Number of faces *) NP: INTEGER; (* Number of polyhedra *) NFE: INTEGER; (* Number of facetedges *) (* der: CARD; (* Edge Ring Degree *) *) (* bdr: CARD; (* 0 = without border, 1 = with border, 2 = cells are octahedra *) *) (* Elements: *) vertex: REF Vertices; edge: REF ARRAY OF Edge; face: REF ARRAY OF Face; polyhedron: REF ARRAY OF Polyhedron; facetedge: REF ARRAY OF Pair; (* one pair for each facetedge *) out: REF ARRAY OF Pair; (* one pair for each vertex where: Org(out[v]) = vertex[v] *) region: REF ARRAY OF Pair; (* one pair for each polyhedron where: Org(region[r]) = polyhedron[r] *) END; TopoDict = RECORD OldVertNum : REF ARRAY OF CARD; OldEdgeNum : REF ARRAY OF CARD; OldFaceNum : REF ARRAY OF CARD; OldPolyNum : REF ARRAY OF CARD; END; PROCEDURE MakeTopology(a: Pair): Topology; (* Builds the tabular representation of the tridimensional map reachable from "a". Assumes that "NumberVertices", "Octf.NumberEdges", "Octf.NumberFacets", "Octf.NumberFacetEdges" have been called, after the topology of the map was modified for the last time. *) TYPE PolyhedronTopology = RECORD NV: CARD; (* Number of vertices *) NE: CARD; (* Number of edges *) NF: CARD; (* Number of faces *) (* The following Pairs have Pneg = the polyhedron in question *) vRef: REF ARRAY OF Pair; (* One Pair out of each vertex *) eRef: REF ARRAY OF Pair; (* One Pair along each edge *) fRef: REF ARRAY OF Pair; (* One Pair on the boundary of each face *) END; PROCEDURE MakePolyhedronTopology(a: Pair): PolyhedronTopology; (* Returns the elements of the boundary of the polyhedron Pneg(Dual(a)) = Org(a). Note that `a.facetedge.edge' is an edge of the dual map. *) TYPE TetraCorners = ARRAY[0..11] OF Pair; PROCEDURE GetTetraCorners(a : Pair) : TetraCorners; (* This procedure builds an array of 12 facetedge pairs, such all elements on the array have the same handness and the same negative polyhedron "Pneg". The pair "a" in the argument must be a pair that represent one tetrahe- dron of $T$ (in the primal space). *) PROCEDURE CollectTetrahedra(READONLY tp: Topology) : REF ARRAY OF Pair; (* Returns a list "t" with one facetedge on each tetrahedron from tp in numerical order, with consistent orientations whenever possible. The facetedges will be such that Pneg(t[i]).num = i. *) (* === LOW-LEVEL ENUMERATION PROCEDURES === *) PROCEDURE CollectFaceEdges(a: Pair; VAR re: REF ARRAY OF Pair; VAR ne: CARDINAL); (* Enumerates the edges of face "a.face" (with successive "Enext"s) and stores them in "re[0..ne-1]". Also sets "ne" and (re)allocates "re^" as needed. *) PROCEDURE CollectPolyhedronFaces(a: Pair; VAR rf: REF ARRAY OF Pair; VAR nf: CARDINAL); (* Enumerates the faces of polyhedron "Pneg(a)" and stores them in "rf[0..nf-1]". Also sets "nf" and (re)allocates "rf^" as needed. Assumes that the "xmark" bits of all face records are FALSE. *) PROCEDURE CollectPolyhedronEdges( a: Pair; VAR re: REF ARRAY OF Pair; VAR ne: CARDINAL; VAR rf: REF ARRAY OF Pair; (* Working area *) ); (* Collects all edges of the polyhedron "Pneg(a)", returns their representative pairs in "re[0..ne-1]". Also sets "ne" and (re)allocates "re^" as needed. Assumes all "xmark" bits are FALSE. The "rf" vector is a working area, expanded as needed; may be NIL initially. *) PROCEDURE CollectPolyhedronVertices( a: Pair; VAR rv: REF ARRAY OF Pair; VAR nv: CARDINAL; VAR rf: REF ARRAY OF Pair; (* Working area *) ); (* Collects all vertices of the polyhedron "Pneg(a)", returns their "OrgV"-representatives in "re[0..ne-1]". Also sets "ne" and (re)allocates "re^" as needed. Assumes all "xmark" bits are FALSE. The "rf" vector is a working area, expanded as needed; may be NIL initially. *) (* === INPUT/OUTPUT === *) PROCEDURE DegreeOfVertex(a: Pair) : CARD ; (* Compute the degree of vertex that is Org(a) (i.e. the number of component edges incident to vertex). *) TYPE AdjacencyMatrix = ARRAY OF ARRAY OF BOOL; PROCEDURE MakeAdjacencyMatrix(READONLY top: Topology): REF AdjacencyMatrix; (* Builds the adjacency matrix for the topology "top". *) PROCEDURE TriviallyIsomorphic(READONLY ta, tb: Topology): BOOL; (* True iff "ta" and "tb" are topologically isomorphic, with the trivial isomorphism (that is, if elements with same index have the same topological relationship in both). *) PROCEDURE GetVariableVertices( READONLY top: Topology; VAR vr: ARRAY OF BOOL; ); (* Sets "vr[v] := TRUE" for every vertex "v" that is not fixed. *) (* === GEOMETRIC TOOLS === *) TYPE Coords = ARRAY OF LR4.T; PROCEDURE InitCoords(coins: Random.T; VAR c: Coords; radius: LONGREAL := 1.0d0); (* Fills c with random coordinates in the range [-radius __ +radius]. *) PROCEDURE GenCoords(READONLY t: Topology): REF Coords; (* Allocates a coordinate vector "c", and initializes it using "InitCoords" with random numbers in "[-1.0 __ +1.0]". *) PROCEDURE Displace(READONLY top: Topology; d: LR4.T; VAR c: Coords); (* Displaces all existing vertices by "d". *) PROCEDURE Scale(READONLY top: Topology; s: LONGREAL; VAR c: Coords); (* Scales the coordinates of all existing vertices by "s". *) PROCEDURE FaceCross(a: Pair; READONLY c: Coords): LR4.T; (* A vector approximately perpendicular to component face of pair "a", this vector is computing by the mean of two perpendicular vectors (with orientation topological consistent) for two tetrahedrons inci- dents in this face. Returns the unit vector if the face has atributte "exists=FALSE". *) PROCEDURE FaceNormal(a: Pair; READONLY c: Coords): LR4.T; (* Normal of component face of "a"; same as "LR4.Dir(FaceCross(a, c))". Returns an arbitrary unit vector if the face has zero area. *) PROCEDURE PolyCross(a: Pair; READONLY c: Coords): LR4.T; (* A vector approximately perpendicular to Pneg of pair "a". Returns the unit vector if the polyhedron has atributte "exists=FALSE". *) PROCEDURE PolyNormal(a: Pair; READONLY c: Coords): LR4.T; (* Normal of "Pneg(a)"; same as "LR4.Dir(PolyCross(a, c))". Returns an arbitrary unit vector if the face has zero area *) PROCEDURE EdgeCross(a: Pair; READONLY c: Coords): LR4.T; (* A vector approximately perpendicular to component edge of pair "a", this vector is computing by the mean of perpendicular vectors (with orientation topological consistent) to all tetrahedrons incidents in this edge. Returns the unit vector if the edge has atributte "exists=FALSE". *) PROCEDURE EdgeNormal(a: Pair; READONLY c: Coords): LR4.T; (* Normal of component edge of "a"; same as "LR4.Dir(EdgeCross(a, c))". Returns an arbitrary unit vector if the edge has zero lenght. *) PROCEDURE FaceBarycenter(a: Pair; READONLY c: Coords): LR4.T; (* The barycenter of face component of pair "a". *) PROCEDURE TetraBarycenter(a: Pair; READONLY c: Coords): LR4.T; (* The barycenter of the negative tetrahedron of pair "a". *) PROCEDURE VertexCross(a: Pair; READONLY c: Coords; READONLY top: Topology): LR4.T; (* A vector approximately orthogonal to the Org(a), whose length is proportional the volume of the polyhedron defined by the "existing" neighbors of "Org(a)". *) PROCEDURE VertexNormal(a: Pair; READONLY c: Coords; READONLY top: Topology): LR4.T; (* Estimated normal at Org(a), considering only neighbors that exist; same as "LR4.Dir(VertexCross(a, c))". Returns an arbitrary unit vector if "VertexCross(a, c)" is zero. *) TYPE Quadp = ARRAY [0..3] OF Pair; Quadv = ARRAY [0..3] OF Vertex; Triv = ARRAY [0..2] OF Vertex; Trip = ARRAY [0..2] OF Pair; PROCEDURE Neighbors( v: Vertex; READONLY top: Topology; ): REF Vertices; (* Neighbor vertices of "v", for any topology. *) PROCEDURE StarOfVertex(a: Pair; READONLY top: Topology): REF ARRAY OF Quadv; (* Return the vertex set that conform every tetrahedron belonging to star of OrgV(a). The vertex set is sorted such as, the first element is OrgV(a), the second is OrgV(Enext(a)), the third OrgV(Enext_1(a)) and the las element is OrgV(Enext_1(Fnext_1(a))). *) PROCEDURE NumberPolyOfStar(quadv: REF ARRAY OF Quadv): CARD; (* Return the number of polyhedrons that belong to star of OrgV(a). *) PROCEDURE ComputeAllVertexNormals(READONLY top: Topology; READONLY c: Coords ): REF ARRAY OF LR4.T; (* Returns a vector with the result of VertexNormal applied to each vertex of "top". *) PROCEDURE ComputeAllEdgeNormals(READONLY top: Topology; READONLY c: Coords ): REF ARRAY OF LR4.T; (* Returns a vector with the result of EdgeNormal applied to each edge of "top". *) PROCEDURE ComputeAllFaceNormals(READONLY top: Topology; READONLY c: Coords ): REF ARRAY OF LR4.T; (* Returns a vector with the result of Face Normal applied to each face of "top". *) PROCEDURE ComputeAllPolyhedronNormals(READONLY top: Topology; READONLY c: Coords): REF ARRAY OF LR4.T; (* Returns a vector with the result of PolyNormal applied to each polyhe- dron of "top". *) (* === GEOMETRIC AVERAGES === *) (* In the following procedures, when "all = TRUE" the quantities are computed over all elements of the appropriate dimension; otherwise they are computed only over those elements with "existing = TRUE". *) PROCEDURE Barycenter(READONLY top: Topology; READONLY c: Coords; all: BOOL): LR4.T; (* Returns the barycenter of all vertices. *) PROCEDURE PartialBarycenter(READONLY v: Vertices; READONLY c: Coords; all: BOOL): LR4.T; (* Barycenter of the vertices "v[i]". Previously called "NeighborBarycenter". *) PROCEDURE MeanVertexDistance( READONLY top: Topology; READONLY c: Coords; READONLY ctr: LR4.T; all: BOOL ): LONGREAL; (* The average distance of vertices from "ctr", in the root-mean-square sense; that is, "sqrt(sum(norm(c[v]-ctr)^2)/N)". *) PROCEDURE MaxVertexDistance( READONLY top: Topology; READONLY c: Coords; READONLY ctr: LR4.T; all: BOOL; ): LONGREAL; (* The maximum distance of a vertex from "ctr". *) PROCEDURE MeanEdgeLength(READONLY top: Topology; READONLY c: Coords; all: BOOL): LONGREAL; (* The average length of edges, in the root-mean-square sense; that is "sqrt(sum(dist(c[org(e)], c[dst(e)])^2)/N)". *) PROCEDURE MeanPolyhedronNormal(READONLY top: Topology; READONLY c: Coords; all: BOOL): LR4.T; (* The average of all tetrahedra normals, normalized to unit length. *) PROCEDURE MeanThickness( READONLY top: Topology; READONLY c: Coords; READONLY ctr: LR4.T; READONLY norm: LR4.T; all: BOOL; ): LONGREAL; (* Average vertex deviation from the hyperplane that passes through "ctr" and is orthogonal to "norm", in the root mean square sense. *) PROCEDURE NormalizeVertexDistances( READONLY top: Topology; VAR c: Coords; all: BOOL; ); (* Shifts and scales all vertices so that they have barycenter (0,0,0,0) and unit mean square distance from the origin. The "all" parameter is used to compute the barycenter and the mean square distance. *) PROCEDURE NormalizeEdgeLengths( READONLY top: Topology; VAR c: Coords; all: BOOL ); (* Shifts and scales all existing vertices so that they have barycenter (0,0,0,0), and the mean square length of existing edges is 1.0. The "all" parameter is used to compute the barycenter and the mean edge length. *) (* === INPUT/OUTPUT === *) TYPE TopCom = RECORD top: Topology; comments: TEXT; END; PROCEDURE WriteTopology( name: TEXT; READONLY top: Topology; comments: TEXT := " "; ); (* Writes "top", and "comments" to file disk in a format that can be read back. The file will have the given "name" with ".tp" appended. *) PROCEDURE WriteState( name: TEXT; READONLY top: Topology; READONLY c: Coords; comments: TEXT := " "; ); (* Writes the coordinates geometrics to file disk in a format that can be read back. The file will have the given "name" with ".st" appended. *) PROCEDURE WriteMaterials( name: TEXT; READONLY top: Topology; comments: TEXT := " "; ro_te: BOOL := FALSE; ); (* Writes the materials properties to file disk in a format that can be read back. The file will have the given "name" with ".ma" appended. *) PROCEDURE WriteDualTopology( name: TEXT; READONLY top: Topology; comments: TEXT := " "; ); (* Writes the dual of "top", and "comments" to file disk in a format that can be read back. The file will have the given "name" with ".tp" appended. *) PROCEDURE WriteStDe( wr: Wr.T; READONLY c: Coords; (* Vertex coordinates *) READONLY Dc: Coords; (* Vertex coordinates derivates *) prec: CARD := 4; (* Significant figures to use for "c" and "Dc" *) comments: TEXT := ""; ); (* Writes the coordinates "c" and derivatives "Dc" to file "wr" in the ".sd" format. *) PROCEDURE ReadToMa(name: TEXT; ro_te: BOOL := FALSE): TopCom; (* A more general procedure to read two disk files created by "WriteTopology" and "WriteMaterials". The files must have the given "name" with ".tp" and ".ma". *) PROCEDURE ReadState(name: TEXT): REF Coords; (* Reads a disk file created by "WriteState". The file must have the given "name" with ".st" appended. *) END Triangulation. (**************************************************************************) (* *) (* Copyright (C) 2001 Universidade Estadual de Campinas (UNICAMP) *) (* *) (* Authors: *) (* L. P. Lozada & J. Stolfi - UNICAMP *) (* *) (* This file can be freely used, distributed, and modified, provided that *) (* this copyright and authorship notice is included in every copy or *) (* derived version. *) (* *) (* DISCLAIMER: This software is offered ``as is'', without any guarantee *) (* as to fitness for any particular purpose. Neither the copyright *) (* holder nor the authors or their employers can be held responsible *) (* for any damages that may result from its use. *) (* *) (* Last edited on 2001-05-21 02:22:21 by stolfi *) (**************************************************************************)