From 6bb170c7ff8d5fccd00d30b03e63159d6abcaf1b Mon Sep 17 00:00:00 2001 From: Chris Hodapp Date: Sun, 24 May 2020 22:43:32 -0400 Subject: [PATCH] Half-made DCEL; half-broken continuous transform thingy --- README.md | 10 +- src/dcel.rs | 313 ++++++++++++++++++++++++++++++++++++++++++++++++ src/examples.rs | 30 ++++- src/lib.rs | 8 +- src/rule.rs | 298 ++++++++++++++++++++++++++------------------- 5 files changed, 532 insertions(+), 127 deletions(-) create mode 100644 src/dcel.rs diff --git a/README.md b/README.md index 45b3c62..0819acd 100644 --- a/README.md +++ b/README.md @@ -2,8 +2,12 @@ ## Highest priority: +- Work on my doubly-connected edge list so I can complete some other + things below! - Implement the continuous parametric transformations from 2020-05-07 in my notes. This will require some new abstractions. +- Document `parametric_mesh` better. It is pretty hairy and could + also benefit from some modularity. - Get identical or near-identical meshes to `ramhorn_branch` from Python. (Should just be a matter of tweaking parameters.) - Look at performance. @@ -56,6 +60,10 @@ ## Research Areas +- [Geometry and Algorithms for Computer Aided Design (Hartmann)](https://www2.mathematik.tu-darmstadt.de/~ehartmann/cdgen0104.pdf) +- https://en.wikipedia.org/wiki/Surface_triangulation +- https://www.cs.cmu.edu/~quake/triangle.html + ## Reflections & Quick Notes - Generalizing to space curves moves this away from the "discrete @@ -63,4 +71,4 @@ for discrete automata. - If you *pre* multiply a transformation: you are transforming the entire global space. If you *post* multiply: you are transforming - the current local space. \ No newline at end of file + the current local space. diff --git a/src/dcel.rs b/src/dcel.rs new file mode 100644 index 0000000..e9b951b --- /dev/null +++ b/src/dcel.rs @@ -0,0 +1,313 @@ +use std::fmt; + +use crate::mesh::{Mesh}; +use crate::xform::{Vertex}; + +/// Doubly-connected edge list mesh (or a half-edge mesh), +/// parametrized over some vertex type. +#[derive(Clone, Debug)] +pub struct DCELMesh { + pub verts: Vec>, + pub faces: Vec, + pub halfedges: Vec, + + pub num_verts: usize, + pub num_faces: usize, + pub num_halfedges: usize, +} + +/// A vertex of a mesh, combined with an arbitrary half-edge that has +/// this vertex as its origin. This is always relative to some parent +/// Mesh. +#[derive(Clone, Debug)] +pub struct DCELVertex { + /// The vertex itself. + pub v: V, + /// A half-edge (given as an index into 'halfedges'); + /// arbitrary, but `mesh.halfedges[halfedge] = v` must be true + pub halfedge: usize, +} + +/// A face, given as a half-edge that lies on its boundary (and must +/// traverse it counter-clockwise). This is always relative to some +/// parent Mesh, as in Vertex. +#[derive(Clone, Debug)] +pub struct DCELFace { + /// A boundary half-edge of this face (given as an index into + /// 'halfedges'). + pub halfedge: usize, +} + +/// A half-edge, given in terms of its origin vertex, the face that the +/// half-edge lies on the boundary of, its optional "twin" half-edge that +/// lies on an adjacent face, and previous and next half-edges (to +/// traverse the boundaries of the face). This is always relative to +/// some parent Mesh, as in Vertex and Face. +#[derive(Clone, Debug)] +pub struct DCELHalfEdge { + /// Origin vertex (given as an index into 'verts') + pub vert: usize, + /// Face this half-edge lies on the boundary of (given as an index + /// into 'faces') + pub face: usize, + /// If false, ignore twin_halfedge. (If this is true, then it must + /// also be true for the twin.) + pub has_twin: bool, + /// The twin half-edge (given as an index into 'halfedges'). + /// The twin of the twin must point back to this HalfEdge. + pub twin_halfedge: usize, + /// The next half-edge on the boundary (given as an index into + /// 'halfedges'). 'prev_halfedge' of this half-edge must point + /// back to this same HalfEdge. Repeatedly following 'next_halfedge' + /// must also lead back to this same HalfEdge. + pub next_halfedge: usize, + /// The previous half-edge on the boundary (given as an index into + /// 'halfedges'). 'next_halfedge' of this half-edge must point + /// back to this HalfEdge. Repeatedly following 'prev_halfedge' + /// must also lead back to this same HalfEdge. + pub prev_halfedge: usize, +} + +impl fmt::Display for DCELMesh { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + + let v_strs: Vec = self.verts.iter().enumerate().map(|(i,v)| { + format!("v{}={:?}", i, v.v) + }).collect(); + let v_str = v_strs.join(","); + + let f_strs: Vec = self.faces.iter().enumerate().map(|(i,f)| { + format!("f{}=e{}", i, f.halfedge) + }).collect(); + let f_str = f_strs.join(","); + + let e_strs: Vec = self.halfedges.iter().enumerate().map(|(i,h)| { + let twin = if h.has_twin { + format!(" tw{}", h.twin_halfedge) + } else { + String::from("") + }; + format!("e{}=v{} f{}{} n{} p{}", i, h.vert, h.face, twin, h.next_halfedge, h.prev_halfedge) + }).collect(); + let e_str = e_strs.join(","); + + write!(f, "DCELMesh({} verts, {}; {} faces, {}; {} halfedges, {})", + self.num_verts, v_str, + self.num_faces, f_str, + self.num_halfedges, e_str) + } +} + +impl DCELMesh { + + pub fn new() -> DCELMesh { + DCELMesh { + verts: vec![], + faces: vec![], + halfedges: vec![], + num_verts: 0, + num_faces: 0, + num_halfedges: 0, + } + } + + pub fn face_to_halfedges(&self, face_idx: usize) -> Vec { + let mut edges: Vec = vec![]; + let start_idx = self.faces[face_idx].halfedge; + edges.push(start_idx); + + let mut idx = self.halfedges[start_idx].next_halfedge; + while start_idx != idx { + edges.push(idx); + idx = self.halfedges[idx].next_halfedge; + } + return edges; + } + + pub fn face_to_verts(&self, face_idx: usize) -> Vec { + self.face_to_halfedges(face_idx).iter().map(|e| { + self.halfedges[*e].vert + }).collect() + } + + /// Add a face with no connections + pub fn add_face(&mut self, verts: [V; 3]) -> usize { + // Vertices will be at indices v_n, v_n+1, v_n+2: + let v_n = self.num_verts; + // The face will be at index f_n: + let f_n = self.num_faces; + // The half-edges will be at indices e_n, e_n+1, e_n+2: + let e_n = self.num_halfedges; + + // Half-edges and vertices can be inserted both at once: + for i in 0..3 { + let n = (i + 1) % 3; + let p = (i + 2) % 3; + self.halfedges.push(DCELHalfEdge { + vert: v_n + i, + face: f_n, + has_twin: false, + twin_halfedge: 0, + next_halfedge: e_n + n, + prev_halfedge: e_n + p, + }); + self.verts.push(DCELVertex { + v: verts[i], + halfedge: e_n + i, + }); + } + self.num_halfedges += 3; + self.num_verts += 3; + + // Finally, add the face (any halfedge is fine): + self.faces.push(DCELFace { halfedge: e_n }); + self.num_faces += 1; + + f_n + } + + /// Add a face that lies on an existing boundary - i.e. one half-edge + /// has a twin half-edge already on the mesh. As this gives two + /// vertices, only one other vertex needs specified. + pub fn add_face_twin1(&mut self, twin: usize, vert: V) -> usize { + // 'vert' will be at index v_n: + let v_n = self.num_verts; + // The face will be at index f_n: + let f_n = self.num_faces; + // The half-edges will be at indices e_n, e_n+1, e_n+2: + let e_n = self.num_halfedges; + + // Note the reversal of direction + let twin_halfedge = &self.halfedges[twin]; + let v1 = self.halfedges[twin_halfedge.next_halfedge].vert; + let v2 = twin_halfedge.vert; + + // Insert 'twin' half-edge first: + self.halfedges.push(DCELHalfEdge { + vert: v1, + face: f_n, + has_twin: true, + twin_halfedge: twin, + next_halfedge: e_n + 1, + prev_halfedge: e_n + 2, + }); + self.halfedges[twin].has_twin = true; + self.halfedges[twin].twin_halfedge = e_n; + self.halfedges.push(DCELHalfEdge { + vert: v2, + face: f_n, + has_twin: false, + twin_halfedge: 0, + next_halfedge: e_n + 2, + prev_halfedge: e_n, + }); + self.halfedges.push(DCELHalfEdge { + vert: v_n, + face: f_n, + has_twin: false, + twin_halfedge: 0, + next_halfedge: e_n, + prev_halfedge: e_n + 1, + }); + + self.num_halfedges += 3; + + // Since the 2nd halfedge we inserted (e_n + 1) has origin v_n: + self.verts.push(DCELVertex { + v: vert, + halfedge: e_n + 1, + }); + self.num_verts += 1; + + // Finally, add the face (any halfedge is fine): + self.faces.push(DCELFace { halfedge: e_n }); + self.num_faces += 1; + + f_n + } + + /// Add a face that lies on two connected boundaries - i.e. two of its + /// half-edges have twins already on the mesh. + /// + /// Twin half-edges should be given in counter-clockwise order; that + /// is, for the resultant face, one half-edge's twin will be twin1, and + /// the next half-edge's twin will be twin2. + /// Also: `self.halfedges[twin2].next_halfedge` must equal `twin1`. + pub fn add_face_twin2(&mut self, twin1: usize, twin2: usize) -> usize { + // The face will be at index f_n: + let f_n = self.num_faces; + // The half-edges will be at indices e_n, e_n+1, e_n+2: + let e_n = self.num_halfedges; + + // The origin vertex is 'between' the two edges, but because their + // order is reversed (as twins), this is twin1's origin: + let twin_halfedge = &self.halfedges[twin1]; + let v1 = twin_halfedge.vert; + let v2 = self.halfedges[twin_halfedge.next_halfedge].vert; + // Final vertex is back around to twin2's origin: + let v3 = self.halfedges[twin2].vert; + + self.halfedges.push(DCELHalfEdge { + vert: v1, + face: f_n, + has_twin: true, + twin_halfedge: twin1, + next_halfedge: e_n + 1, + prev_halfedge: e_n + 2, + }); // index e_n + self.halfedges[twin1].has_twin = true; + self.halfedges[twin1].twin_halfedge = e_n; + self.halfedges.push(DCELHalfEdge { + vert: v2, + face: f_n, + has_twin: false, + twin_halfedge: 0, + next_halfedge: e_n + 2, + prev_halfedge: e_n, + }); // index e_n + 1 + self.halfedges.push(DCELHalfEdge { + vert: v3, + face: f_n, + has_twin: true, + twin_halfedge: twin2, + next_halfedge: e_n, + prev_halfedge: e_n + 1, + }); // index e_n + 2 + self.halfedges[twin2].has_twin = true; + self.halfedges[twin2].twin_halfedge = e_n + 2; + self.num_halfedges += 3; + + // Finally, add the face (any halfedge is fine): + self.faces.push(DCELFace { halfedge: e_n }); + self.num_faces += 1; + + f_n + } +} + +pub fn convert_mesh(m: &DCELMesh) -> Mesh { + let n = m.faces.len(); + let mut faces: Vec = vec![0; 3 * n]; + + for i in 0..n { + + let e0 = m.faces[i].halfedge; + let h0 = &m.halfedges[e0]; + faces[3*i + 0] = h0.vert; + let e1 = h0.next_halfedge; + let h1 = &m.halfedges[e1]; + faces[3*i + 1] = h1.vert; + let e2 = h1.next_halfedge; + let h2 = &m.halfedges[e2]; + faces[3*i + 2] = h2.vert; + if h2.next_halfedge != e0 { + panic!(format!("Face {}: half-edges {},{},{} return to {}, not {}", + i, e0, e1, e2, h2.next_halfedge, e0)); + } + } + + Mesh { + verts: m.verts.iter().map(|e| e.v).collect(), + faces: faces, + } +} \ No newline at end of file diff --git a/src/examples.rs b/src/examples.rs index 1205ecf..94e8619 100644 --- a/src/examples.rs +++ b/src/examples.rs @@ -11,6 +11,7 @@ use crate::mesh::{Mesh, MeshFunc, VertexUnion, vert_args}; use crate::xform::{Transform, Vertex, vertex, Mat4, id}; use crate::rule::{Rule, RuleFn, RuleEval, Child}; use crate::prim; +use crate::dcel; /* pub fn cube_thing() -> Rule<()> { @@ -1266,6 +1267,8 @@ pub fn test_parametric() -> Mesh { vertex( 0.5, 0.5, 0.0), vertex( 0.5, -0.5, 0.0), ]; + //let base_verts = util::subdivide_cycle(&base_verts, 2); + //let base_verts = util::subdivide_cycle(&base_verts, 16); let t0 = 0.0; let t1 = 16.0; @@ -1275,5 +1278,30 @@ pub fn test_parametric() -> Mesh { scale((0.8).powf(t)) }; - crate::rule::parametric_mesh(base_verts, xform, t0, t1, 0.0001) + crate::rule::parametric_mesh(base_verts, xform, t0, t1, 0.005) +} + +pub fn test_dcel(fname: &str) { + let mut mesh: dcel::DCELMesh = dcel::DCELMesh::new(); + let f1 = mesh.add_face([ + vertex(-0.5, -0.5, 0.0), + vertex(-0.5, 0.5, 0.0), + vertex( 0.5, 0.5, 0.0), + ]); + let f2 = mesh.add_face_twin1(mesh.faces[f1].halfedge, vertex(0.0, 0.0, 1.0)); + + let vl1 = mesh.face_to_verts(f1); + println!("verts = {:?}", vl1); + let vl2 = mesh.face_to_verts(f2); + println!("verts = {:?}", vl2); + + //let f3 = mesh.add_face_twin2(); + + println!("DCEL mesh = {}", mesh); + + let mesh_conv = dcel::convert_mesh(&mesh); + + println!("Mesh = {:?}", mesh_conv); + + mesh_conv.write_stl_file(fname).unwrap(); } \ No newline at end of file diff --git a/src/lib.rs b/src/lib.rs index 9281a3e..7116043 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -6,6 +6,7 @@ pub mod prim; pub mod util; pub mod xform; pub mod examples; +pub mod dcel; //pub use crate::examples; //pub use crate::openmesh::test_thing; @@ -138,7 +139,12 @@ mod tests { #[test] fn test_parametric() { - examples::test_parametric().write_stl_file("test_parametric.stl"); + examples::test_parametric().write_stl_file("test_parametric.stl").unwrap(); + } + + #[test] + fn test_dcel() { + examples::test_dcel("test_dcel.stl"); } } // need this for now: diff --git a/src/rule.rs b/src/rule.rs index f5e0781..1849fc8 100644 --- a/src/rule.rs +++ b/src/rule.rs @@ -373,18 +373,45 @@ pub fn parametric_mesh(frame: Vec, f: F, t0: f32, t1: f32, max_err: f // neighboring ones. mesh_idx: usize, // Index of this vertex in the mesh // (If it's in 'frontier', the vertex is in the mesh). - neighbor1: usize, // Index of 'frontier' of one neighbor - neighbor2: usize, // Index of 'frontier' of other neighbor + neighbor: [usize; 2], // Indices of 'frontier' of each neighbor + side_faces: [(usize,bool); 2], // The two faces (given as an index + // into 'faces' below) which contain the edge from this vertex to + // its neighbors. If the bool is false, there is no face. }; + // A face that is still undergoing subdivision. This is used for a stack + // in the main loop. 'verts' gives the index of vertices of the face, + // and 'shared_faces' gives neighboring faces, i.e. those which share + // an edge with the face. This always refers to a face that is already + // in 'faces' - but this face may be replaced in the course of + // processing. + #[derive(Clone, Debug)] + struct tempFace { + // Indices into 'verts' below: + verts: [usize; 3], + // The parameter values corresponding to 'verts': + ts: [f32; 3], + // The 'frame' vertices (i.e. vertex at f(t0)) corresponding + // to 'verts': + frame_verts: [Vertex; 3], + // Index into 'faces' below for the starting vertex: + face: usize, + // If the bool is true, this gives an index into 'faces' below for + // a face that shares an edge with this face; if the bool is false, + // disregard (there is no neighbor here). This goes in a specific + // order: the face sharing edge (verts[0], verts[1]), then edge + // (verts[1], verts[2]), then edge (verts[2], verts[0]). + shared_faces: [(usize, bool); 3], + } + // Init 'frontier' with each 'frame' vertex, and start it at t=t0. let mut frontier: Vec = frame.iter().enumerate().map(|(i,v)| frontierVert { vert: *v, t: t0, frame_vert: *v, mesh_idx: i, - neighbor1: (i - 1) % n, - neighbor2: (i + 1) % n, + neighbor: [(i - 1) % n, (i + 1) % n], + side_faces: [(0, false); 2], // initial vertices have no faces }).collect(); // Every vertex in 'frontier' has a trajectory it follows - which is // simply the position as we transform the original vertex by f(t), @@ -399,19 +426,20 @@ pub fn parametric_mesh(frame: Vec, f: F, t0: f32, t1: f32, max_err: f let mut verts: Vec = frame.clone(); let mut faces: Vec = vec![]; - let mut count = 0; + // 'stack' is cleared at every iteration below, and contains faces that + // may still require subdivision. + let mut stack: Vec = vec![]; while !frontier.is_empty() { - - for (i,f) in frontier.iter().enumerate() { - println!("DEBUG: frontier[{}]: vert={},{},{} t={} mesh_idx={} neighbors={},{}", i, f.vert.x, f.vert.y, f.vert.z, f.t, f.mesh_idx, f.neighbor1, f.neighbor2); + for (i, f) in frontier.iter().enumerate() { + println!("DEBUG: frontier[{}]: vert={},{},{} t={} mesh_idx={} neighbor={:?}", i, f.vert.x, f.vert.y, f.vert.z, f.t, f.mesh_idx, f.neighbor); } // Pick a vertex to advance. // // Heuristic for now: pick the 'furthest back' (lowest t) - let (i,v) = { - let (i,v) = frontier.iter().enumerate().min_by(|(i,f), (j, g)| + let (i, v) = { + let (i, v) = frontier.iter().enumerate().min_by(|(i, f), (j, g)| f.t.partial_cmp(&g.t).unwrap_or(std::cmp::Ordering::Equal)).unwrap(); (i, v.clone()) }; @@ -420,9 +448,15 @@ pub fn parametric_mesh(frame: Vec, f: F, t0: f32, t1: f32, max_err: f if v.t >= t1 { break; } + // TODO: Fix boundary behavior here and make sure final topology + // is right. println!("DEBUG: Moving vertex {}, {:?} (t={}, frame_vert={:?})", i, v.vert, v.t, v.frame_vert); + // Move this vertex further along, i.e. t + dt. (dt is set by + // the furthest we can go while remaining within 'err', i.e. when we + // make our connections we look at how far points on the *edges* + // diverge from the trajectory of the continuous transformation). let mut dt = (t1 - t0) / 100.0; let vf = v.frame_vert; for iter in 0..100 { @@ -432,9 +466,9 @@ pub fn parametric_mesh(frame: Vec, f: F, t0: f32, t1: f32, max_err: f // // If we assume some continuity in f, then we can guess that // the worst error occurs at the midpoint of the edge: - let edge_mid = 0.5*(f(v.t).mtx + f(v.t + dt).mtx)*vf; + let edge_mid = 0.5 * (f(v.t).mtx + f(v.t + dt).mtx) * vf; // ...relative to the trajectory midpoint: - let traj_mid = f(v.t + dt/2.0).mtx * vf; + let traj_mid = f(v.t + dt / 2.0).mtx * vf; let err = (edge_mid - traj_mid).norm(); //println!("DEBUG iter {}: dt={}, edge_mid={:?}, traj_mid={:?}, err={}", iter, dt, edge_mid, traj_mid, err); @@ -460,133 +494,149 @@ pub fn parametric_mesh(frame: Vec, f: F, t0: f32, t1: f32, max_err: f let t = v.t + dt; let v_next = f(t).mtx * vf; + // DEBUG + /* + let l1 = (v.vert - v_next).norm(); + let l2 = (frontier[v.neighbor1].vert - v.vert).norm(); + let l3 = (frontier[v.neighbor2].vert - v.vert).norm(); + println!("l1={} l2={} l3={}", l1, l2, l3); + */ + // Add this vertex to our mesh: - let pos = verts.len(); + let v_next_idx = verts.len(); verts.push(v_next); + // and add two faces: + let face_idx = faces.len(); + faces.append(&mut vec![ + v_next_idx, v.mesh_idx, frontier[v.neighbor[0]].mesh_idx, // face_idx + v.mesh_idx, v_next_idx, frontier[v.neighbor[1]].mesh_idx, // face_idx + 1 + ]); // Replace this vertex in the frontier: frontier[i] = frontierVert { vert: v_next, frame_vert: vf, - mesh_idx: pos, + mesh_idx: v_next_idx, t: t, - neighbor1: v.neighbor1, - neighbor2: v.neighbor2, + neighbor: v.neighbor, + side_faces: [(face_idx, true), (face_idx + 1, true)], }; - // We might update neighbor1/neighbor2 below. + // Note above that we added two edges that include v_next_idx: + // (frontier[v.neighbor[0]].mesh_idx, v_next_idx) + // (frontier[v.neighbor[1]].mesh_idx, v_next_idx) + // hence, side_faces has the respective face for each. - // The starting vertex (v.vert) is already in the mesh, but - // the two 'neighbor' vertices need to be considered: - let neighbors = [v.neighbor1, v.neighbor2]; + // Also add these faces to the stack of triangles to check for + // subdivision. They may be replaced. + let f0 = &frontier[v.neighbor[0]]; + stack.push(tempFace { + verts: [v_next_idx, v.mesh_idx, f0.mesh_idx], + ts: [t, v.t, f0.t], + frame_verts: [vf, vf, f0.frame_vert], + face: face_idx, + shared_faces: [(face_idx + 1, true), (0, false), f0.side_faces[1]], + }); + let f1 = &frontier[v.neighbor[1]]; + stack.push(tempFace { + verts: [v.mesh_idx, v_next_idx, f1.mesh_idx], + ts: [v.t, t, f1.t], + frame_verts: [vf, vf, f1.frame_vert], + face: face_idx + 1, + shared_faces: [(face_idx, true), f1.side_faces[0], (0, false)], + }); + // Note that vf appears several times in frame_verts because + // several vertices sit in the same trajectory (thus, same 'frame' + // vertex). - // We can consider two faces - one on each side of the 'main' - // edge (v, v_next), going to vertex side1_idx or side2_idx. - // We may also need to subdivide one or both of these faces - // along the edge (v_next, side1_idx) or (v_next, side2_idx). - for (j,neighbor_idx) in neighbors.iter().enumerate() { - let neighbor = &frontier[*neighbor_idx]; - let side_idx = neighbor.mesh_idx; - // To do so, we can consider the error on the edge - // (v_next, neighbor), using mostly the same 'midpoint' method - // as each iteration when computing dt. - // First, find edge midpoint: - let edge_mid = 0.5*(v_next + neighbor.vert); - // To find the 'interpolated' trajectory we need to interpolate - // back at the original frame (we may need this anyway for other - // computations later): - let frame_mid = 0.5*(vf + neighbor.frame_vert); - // ...then trace these to find 'trajectory' midpoint; - let t_mid = v.t + dt/2.0; - let xf = f(t_mid).mtx; - let traj_mid = xf * frame_mid; - // and finally, error: - let mid_err = (edge_mid - traj_mid).norm(); - println!("DEBUG: j={} considering edge {},{}, edge_mid={:?} traj_mid={:?} mid_err={}", j, i, *neighbor_idx, edge_mid, traj_mid, mid_err); - - // DEBUG: This is disabled until fixed - if false && mid_err > max_err { - println!("Error over threshold, splitting ({},{})", i, neighbor_idx); - // Error is high enough we should split. - let pos_new = verts.len(); - verts.push(traj_mid); - let new_neighbor = frontier.len(); - frontier.push(frontierVert { - vert: traj_mid, - frame_vert: frame_mid, - mesh_idx: pos_new, - t: t_mid, - neighbor1: *neighbor_idx, - neighbor2: i, - }); - // We've put a new vertex in between (see neighbor1/2) - // so we have to update this information. - let mut patch_neighbor = |idx: usize, from: usize, to: usize| { - println!("DEBUG: patch_neighbor of idx {} from {} to {}", idx, from, to); - if frontier[idx].neighbor1 == from { frontier[idx].neighbor1 = to; } - if frontier[idx].neighbor2 == from { frontier[idx].neighbor2 = to; } - }; - patch_neighbor(i, *neighbor_idx, new_neighbor); - patch_neighbor(*neighbor_idx, i, new_neighbor); - - let mut new_faces = if j % 2 == 0 { - vec![ - v.mesh_idx, pos, pos_new, - v.mesh_idx, pos_new, side_idx, - ] - } else { - vec![ - pos, v.mesh_idx, pos_new, - pos_new, v.mesh_idx, side_idx, - ] - }; - faces.append(&mut new_faces); + while !stack.is_empty() { + let face = stack.pop().unwrap(); + println!("DEBUG: Examining face: {:?}", face); + let v0 = verts[face.verts[0]]; + let v1 = verts[face.verts[1]]; + let v2 = verts[face.verts[2]]; + let d01 = (v0 - v1).norm(); + let d02 = (v0 - v2).norm(); + let d12 = (v1 - v2).norm(); + // Law of cosines: + let cos0 = (d01*d01 + d02*d02 - d12*d12) / (2.0 * d01 * d02); + let cos1 = (d01*d01 + d12*d12 - d02*d02) / (2.0 * d01 * d12); + let cos2 = (d02*d02 + d12*d12 - d01*d01) / (2.0 * d02 * d12); + // TODO: Perhaps use https://en.wikipedia.org/wiki/Circumscribed_circle#Cartesian_coordinates + // to go by circumradius instead. Or - find it by law of sines? + println!("DEBUG: d01={} d02={} d12={} cos0={} cos1={} cos2={}", d01, d02, d12, cos0, cos1, cos2); + if (cos0 < 0.0 || cos0 > 0.7 || + cos1 < 0.0 || cos1 > 0.7 || + cos2 < 0.0 || cos2 > 0.7) { + println!("DEBUG: Angles out of range!"); } else { - // Error is within limit. Just add one face. - - let mut new_faces = if j % 2 == 0 { - vec![v.mesh_idx, pos, side_idx] - } else { - vec![pos, v.mesh_idx, side_idx] - }; - faces.append(&mut new_faces); - // TODO: This is likely the wrong 'general' way to handle - // winding order but for n=2 it probably is fine + println!("DEBUG: Angles OK"); } - } + // TODO: Figure out how to subdivide in this case - count += 1; - if count > 50 { - break; - } + // The triangle forms a plane. Get this plane's normal vector. + let a = (v0 - v1).xyz(); + let b = (v0 - v2).xyz(); + let normal = a.cross(&b).normalize(); + // Make a new point that is on the surface, but roughly a + // midpoint in parameter space. Exact location isn't crucial. + let t_mid = (face.ts[0] + face.ts[1] + face.ts[2]) / 3.0; + let v_mid = (face.frame_verts[0] + face.frame_verts[1] + face.frame_verts[2]) / 3.0; + let p = f(t_mid).mtx * v_mid; + + let d = p.xyz().dot(&normal); + + println!("DEBUG: t_mid={} v_mid={},{},{} p={},{},{}", t_mid, v_mid.x, v_mid.y, v_mid.z, p.x, p.y, p.z); + println!("DEBUG: d={}", d); + + // DEBUG + /* + let n = verts.len(); + verts.push(p); + faces.push(face.verts[0]); + faces.push(face.verts[1]); + faces.push(n); + */ + if (d <= max_err) { + // This triangle is already in the mesh, and already popped + // off of the stack. We're done. + continue; + } + println!("DEBUG: d > err, splitting this triangle"); + + // The face has 3 edges. We split each of them (making a new + // vertex in the middle), and then make 3 new edges - one + // between each pair of new vertices - to replace the face + // with four smaller faces. + + // This split is done in 'parameter' space: + let pairs = [(0,1), (1,2), (0,2)]; + let mut mids = pairs.iter().map(|(i,j)| { + let t = (face.ts[*i] + face.ts[*j]) / 2.0; + let v = (face.frame_verts[*i] + face.frame_verts[*j]) / 2.0; + f(t).mtx * v + }).collect(); + + // DEBUG + let n = verts.len(); + // Index n+0 is (0,1), n+1 is (1,2), n+2 is (0,2) + verts.append(&mut mids); + faces[face.face] = n + 0; + faces[face.face + 1] = n + 1; + faces[face.face + 2] = n + 2; + faces.extend_from_slice(&[ + face.verts[0], n + 0, n + 2, + face.verts[1], n + 1, n + 0, + face.verts[2], n + 2, n + 1, + //n + 0, n + 1, n + 2, + ]); + + // TODO: Look at shared_faces because these must be split too. + + // TODO: Do I have to correct *other* things in the stack too? + // If they refer to the same face I may be invalidating + // something here! + } } - // Move this vertex further along, i.e. t + dt. (dt is set by - // the furthest we can go while remaining within 'err', i.e. when we - // make our connections we look at how far points on the *edges* - // diverge from the trajectory of the continuous transformation). - - // Add this vertex to the mesh, and connect it to: the vertex we - // started with, and the two neighbors of that vertex. - - // Repeat at "Pick a vertex...". - - // Don't move t + dt past t1. Once a frontier vertex is placed at - // that value of t, remove it. - - - // Missing: Anything about when to subdivide an edge. - // If I assume a good criterion of "when" to subdivide an edge, the - // "how" is straightforward: find the edge's two neighbors in the - // frontier. Trace them back to their 'original' vertices at t=t0 - // (these should just be stored alongside each frontier member), - // produce an interpolated vertex. Produce an interpolated t from - // respective t of the two neighbors in the frontier; use that t - // to move the 'interpolated' vertex along its trajectory. - // - // Add new vertex to mesh (and make the necessary connections) - // and to frontier. - - // But still missing from that: When do I collapse a subdivision - // back down? return Mesh { verts, faces }; } \ No newline at end of file