Half-working pyramid example & roll back DCEL crap

README.md

@@ -2,10 +2,13 @@
 
 ## Highest priority:
 
+- Fix the commented-out tests in `examples.rs`.
 - Just scrap `parametric_mesh` as much as possible and use existing
   tools (e.g. OpenSubdiv) because this DCEL method is just painful for
   what it is and I have some questions on how it can even work
   theoretically.
+- Connect up the `parametric_mesh` stuff that remains, and worry about
+  perfect meshes later.
 - Get identical or near-identical meshes to `ramhorn_branch` from
   Python.  (Should just be a matter of tweaking parameters.)
 - Look at performance.

src/examples.rs

@@ -1,5 +1,5 @@
 use std::rc::Rc;
-use std::f32::consts::{FRAC_PI_2, PI};
+use std::f32::consts::{FRAC_PI_2, FRAC_PI_3, FRAC_PI_6, PI};
 use std::f32;
 
 use nalgebra::*;
@@ -14,7 +14,6 @@ use crate::prim;
 use crate::dcel;
 use crate::dcel::{DCELMesh, VertSpec};
 
-/*
 pub fn cube_thing() -> Rule<()> {
 
     // Quarter-turn in radians:
@@ -38,8 +37,8 @@ pub fn cube_thing() -> Rule<()> {
 
         let xforms = turns.iter().map(|xf| xf.scale(0.5).translate(6.0, 0.0, 0.0));
         RuleEval {
-            geom: Rc::new(prim::cube()),
+            geom: Rc::new(prim::cube().to_meshfunc()),
-            final_geom: Rc::new(prim::empty_mesh()),
+            final_geom: Rc::new(prim::empty_mesh().to_meshfunc()),
             children: xforms.map(move |xf| Child {
                 rule: self_.clone(),
                 xf: xf,
@@ -50,7 +49,6 @@ pub fn cube_thing() -> Rule<()> {
 
     Rule { eval: Rc::new(rec), ctxt: () }
 }
-*/
 
 pub fn barbs(random: bool) -> Rule<()> {
 
@@ -166,6 +164,68 @@ pub fn barbs(random: bool) -> Rule<()> {
     Rule { eval: base, ctxt: () }
 }
 
+pub fn pyramid() -> Rule<()> {
+
+    // base_verts are for a triangle in the XY plane, centered at
+    // (0,0,0), with unit sidelength:
+    let (b0, bn);
+    let rt3 = (3.0).sqrt();
+    let rt6 = (6.0).sqrt();
+    let base_verts: Vec<VertexUnion> = vec_indexed![
+    @b0 VertexUnion::Vertex(vertex( rt3/3.0,      0.0,     0.0)),
+    VertexUnion::Vertex(vertex(    -rt3/6.0, -1.0/2.0,     0.0)),
+    VertexUnion::Vertex(vertex(    -rt3/6.0,  1.0/2.0,     0.0)),
+    VertexUnion::Vertex(vertex(         0.0,      0.0, rt6/3.0)),
+    @bn,
+    ];
+
+    let test = rule_fn!(() => |_s, base_verts| {
+        RuleEval {
+            geom: Rc::new(MeshFunc {
+                verts: base_verts,
+                faces: vec![ 0, 1, 2 ], //,  0, 3, 1,  2, 3, 0,  1, 3, 2],
+            }),
+            final_geom: Rc::new(MeshFunc {
+                verts: vec![],
+                faces: vec![],
+            }),
+            children: vec![],
+        }
+    });
+
+    let base_to_side = |i| {
+        let rt3 = (3.0).sqrt();
+        let rt6 = (6.0).sqrt();
+        let v01 = Unit::new_normalize(Vector3::new(rt3/2.0, 1.0/2.0, 0.0));
+        let axis = 2.0 * (rt3 / 3.0).asin();
+        let angle = 2.0 * FRAC_PI_3 * (i as f32);
+        id().
+            rotate(&Vector3::z_axis(), angle).
+            translate(rt3/18.0, -1.0/6.0, rt6/9.0).
+            rotate(&v01, axis)
+    };
+
+    let base = rule_fn!(() => |_s, base_verts| {
+        RuleEval {
+            geom: Rc::new(MeshFunc {
+                verts: base_verts,
+                faces: vec![ 0, 1, 2,  0, 3, 1,  2, 3, 0,  1, 3, 2],
+            }),
+            final_geom: Rc::new(MeshFunc {
+                verts: vec![],
+                faces: vec![],
+            }),
+            children: vec![
+                child!(rule!(test, ()), base_to_side(0),),
+                child!(rule!(test, ()), base_to_side(1),),
+                child!(rule!(test, ()), base_to_side(2),),
+            ],
+        }
+    });
+
+    Rule { eval: base, ctxt: () }
+}
+
 /*
 // Meant to be a copy of twist_from_gen from Python &
 // automata_scratch, but has since acquired a sort of life of its own
@@ -202,14 +262,14 @@ pub fn twist(f: f32, subdiv: usize) -> Rule<()> {
 
         let seed_next = incr.transform(&seed2);
 
-        let geom = Rc::new(util::zigzag_to_parent(seed_next.clone(), n));
+        //let vc = util::centroid(&seed_next);
-        // TODO: Cleanliness fix - why not just make these return meshes?
+        //let faces = util::connect_convex(0..n, n, true);
-        let (vc, faces) = util::connect_convex(&seed_next, true);
+        let geom = util::parallel_zigzag(seed_next, 0..n, 0..n);
-        let final_geom = Rc::new(OpenMesh {
+        let final_geom = MeshFunc {
-            verts: vec![vc],
+            verts: vec![],
-            alias_verts: vec![],
+            faces: vec![],
-            faces: faces,
+            // TODO: get actual verts here
-        });
+        };
         
         let c = move |self_: Rc<Rule<()>>| -> RuleEval<()> {
             RuleEval {
@@ -231,7 +291,7 @@ pub fn twist(f: f32, subdiv: usize) -> Rule<()> {
     
     let start = move |_| -> RuleEval<()> {
         
-        let child = |incr, dx, i, ang0, div| -> (OpenMesh, Child<()>) {
+        let child = |incr, dx, i, ang0, div| -> (MeshFunc, Child<()>) {
             let xform = Transform::new().
                 rotate(&y, ang0 + (qtr / div * (i as f32))).
                 translate(dx, 0.0, 0.0);
@@ -247,7 +307,7 @@ pub fn twist(f: f32, subdiv: usize) -> Rule<()> {
             // and in the process, generate faces for these seeds:
             let (centroid, f) = util::connect_convex(&vs, false);
             vs.push(centroid);
-            (OpenMesh { verts: vs, faces: f, alias_verts: vec![] }, c)
+            (MeshFunc { verts: vs, faces: f }, c)
         };
         
         // Generate 'count' children, shifted/rotated differently:
@@ -260,7 +320,9 @@ pub fn twist(f: f32, subdiv: usize) -> Rule<()> {
     
     Rule { eval: Rc::new(start), ctxt: () }
 }
+*/
 
+/*
 #[derive(Copy, Clone)]
 pub struct NestSpiral2Ctxt {
     init: bool,
@@ -380,7 +442,7 @@ pub fn nest_spiral_2() -> Rule<NestSpiral2Ctxt> {
         },
     }
 }
-
+*/
 
 #[derive(Copy, Clone)]
 pub struct TorusCtxt {
@@ -389,6 +451,7 @@ pub struct TorusCtxt {
     stack: [Transform; 3],
 }
 
+/*
 pub fn twisty_torus() -> Rule<TorusCtxt> {
     let subdiv = 8;
     let seed = vec![
@@ -399,8 +462,10 @@ pub fn twisty_torus() -> Rule<TorusCtxt> {
     ];
     let xf = Transform::new().rotate(&Vector3::x_axis(), -0.9);
     let seed = util::subdivide_cycle(&xf.transform(&seed), subdiv);
-    
+
     let n = seed.len();
+    let geom = util::parallel_zigzag(seed, 0..n, n..(2*n));
+    // TODO: where are parent Args?
     let geom = Rc::new(util::zigzag_to_parent(seed.clone(), n));
     let (vc, faces) = util::connect_convex(&seed, true);
     let final_geom = Rc::new(OpenMesh {
@@ -481,7 +546,9 @@ pub fn twisty_torus() -> Rule<TorusCtxt> {
         },
     }
 }
+ */
 
+/*
 pub fn twisty_torus_hardcode() -> Rule<()> {
     let subdiv = 8;
     let seed = vec![
@@ -1272,7 +1339,7 @@ pub fn test_parametric() -> Mesh {
     //let base_verts = util::subdivide_cycle(&base_verts, 16);
 
     let t0 = 0.0;
-    let t1 = 15;
+    let t1 = 15.0;
     let xform = |t: f32| -> Transform {
         id().
             translate(0.0, 0.0, t/5.0).

src/lib.rs

@@ -66,20 +66,21 @@ mod tests {
         geom.transform(&xf1).write_stl_file("xform_trans_rot.stl").unwrap();
     }
 
-    /*
     // TODO: These tests don't test any conditions, so this is useful
     // short-hand to run, but not very meaningful as a test.
     #[test]
     fn cube_thing() {
         run_test(examples::cube_thing(), 3, "cube_thing3", false);
     }
-    */
 
     #[test]
     fn barbs() { run_test(examples::barbs(false), 80, "barbs", false); }
 
     #[test]
     fn barbs_random() { run_test(examples::barbs(true), 80, "barbs_random", false); }
+
+    #[test]
+    fn pyramid() { run_test(examples::pyramid(), 3, "pyramid", false); }
     /*
     #[test]
     fn twist() {

src/mesh.rs

@@ -62,7 +62,7 @@ impl Mesh {
         stl_io::write_stl(writer, triangles.iter())
     }
 
-    fn to_meshfunc(&self) -> MeshFunc {
+    pub fn to_meshfunc(&self) -> MeshFunc {
         MeshFunc {
             faces: self.faces.clone(),
             verts: self.verts.iter().map(|v| VertexUnion::Vertex(*v)).collect(),

src/rule.rs

@@ -1,12 +1,9 @@
 use std::borrow::Borrow;
 use std::rc::Rc;
 use std::f32;
-use std::collections::HashMap;
 
 use crate::mesh::{Mesh, MeshFunc, VertexUnion};
 use crate::xform::{Transform, Vertex};
-use crate::dcel;
-use crate::dcel::{DCELMesh, VertSpec};
 
 pub type RuleFn<S> = Rc<dyn Fn(Rc<Rule<S>>) -> RuleEval<S>>;
 
@@ -347,6 +344,8 @@ impl<S> RuleEval<S> {
 
 }
 
+// fa001f47d40de989da6963e442f31c278c88abc8
+
 /// Produce a mesh from a starting frame, and a function `f` which produces
 /// transformations that change continuously over its argument (the range
 /// of which is given by `t0` and `t1`).  By convention, `f(t0)` should
@@ -367,43 +366,23 @@ pub fn parametric_mesh<F>(frame: Vec<Vertex>, f: F, t0: f32, t1: f32, max_err: f
         panic!("frame must have at least 3 vertices");
     }
 
-    #[derive(Copy, Clone, Debug)]
-    struct VertexTrajectory {
-        // Vertex position
-        vert: Vertex,
-        // Parameter value; f(t) should equal vert
-        t: f32,
-        // "Starting" vertex position, i.e. at f(t0). Always either a frame
-        // vertex, or a linear combination of two neighboring ones.
-        frame_vert: Vertex,
-    };
-
-    let mut mesh: DCELMesh<VertexTrajectory> = DCELMesh::new();
-
-    #[derive(Clone, Debug)]
     struct frontierVert {
-        traj: VertexTrajectory,
+        vert: Vertex, // Vertex position
-        // If the boundaries on either side of this vertex lie on a face
+        t: f32, // Parameter value; f(t) should equal vert
-        // (which is the case for all vertices *except* the initial ones),
+        frame_idx: usize, // Index of 'frame' this sits in the trajectory of
-        // then this gives the halfedges of those boundaries. halfedges[0]
+        mesh_idx: usize, // Index of this vertex in the mesh
-        // connects the 'prior' vertex on the frontier to this, and
+        neighbor1: usize, // Index of 'frontier' of one neighbor
-        // halfedges[1] connect this to the 'next' vertex on the fronter.
+        neighbor2: usize, // Index of 'frontier' of other neighbor
-        // (Direction matters. If halfedges[0] is given, it must *end* at
-        // 'traj.vert'. If halfedges[1] is given, it must *begin* at 'traj.vert'.)
-        halfedges: [Option<usize>; 2],
-        // If this vertex is already in 'mesh', its vertex index there:
-        vert_idx: Option<usize>,
     };
 
     // Init 'frontier' with each 'frame' vertex, and start it at t=t0.
     let mut frontier: Vec<frontierVert> = frame.iter().enumerate().map(|(i,v)| frontierVert {
-        traj: VertexTrajectory {
+        vert: *v,
-            vert: *v,
+        t: t0,
-            t: t0,
+        frame_idx: i,
-            frame_vert: *v,
+        mesh_idx: i,
-        },
+        neighbor1: (i - 1) % n,
-        halfedges: [None; 2],
+        neighbor2: (i + 1) % n,
-        vert_idx: None,
     }).collect();
     // Every vertex in 'frontier' has a trajectory it follows - which is
     // simply the position as we transform the original vertex by f(t),
@@ -413,73 +392,28 @@ pub fn parametric_mesh<F>(frame: Vec<Vertex>, f: F, t0: f32, t1: f32, max_err: f
     // new triangles every time we advance one, until each vertex
     // reaches t=t1 - in a way that forms the mesh we want.
 
+    // That mesh will be built up here, starting with frame vertices:
+    // (note initial value of mesh_idx)
+    let mut verts: Vec<Vertex> = frame.clone();
+    let mut faces: Vec<usize> = vec![];
+
     while !frontier.is_empty() {
-        /*
+
-        for (i, f) in frontier.iter().enumerate() {
+        // Pick a vertex to advance.
-            println!("DEBUG: frontier[{}]: vert={},{},{} vert_idx={:?} t={} halfedges={:?}", i, f.vert.x, f.vert.y, f.vert.z, f.vert_idx, f.t, f.halfedges);
+        //
-            match f.vert_idx {
+        // Heuristic for now: pick the 'furthest back' (lowest t)
-                Some(vert) => {
+        let (i,v) = frontier.iter().enumerate().min_by(|(i,f), (j, g)|
-                    match f.halfedges[1] {
+            f.t.partial_cmp(&g.t).unwrap_or(std::cmp::Ordering::Equal)).unwrap();
-                        Some(idx) => {
+        // TODO: Make this less ugly?
-                            if mesh.halfedges[idx].vert != vert {
+
-                                println!("  Error: halfedge[1]={} starts at {}, should start at {}", idx, mesh.halfedges[idx].vert, vert);
+        if v.t >= t1 {
-                            }
+            break;
-                        },
-                        None => {},
-                    };
-                    match f.halfedges[0] {
-                        Some(idx) => {
-                            let v2 = mesh.halfedges[mesh.halfedges[idx].next_halfedge].vert;
-                            if v2 != vert {
-                                println!("  Error: halfedge[0]={} ends at {}, should start at {}", idx, v2, vert);
-                            }
-                        },
-                        None => {},
-                    };
-                },
-                None => {},
-            };
         }
-         */
 
-        // Pick a vertex to advance.  For now, this is an easy
+        println!("DEBUG: Moving vertex {}, {:?} (t={}, frame_idx={})", i, v.vert, v.t, v.frame_idx);
-        // heuristic: pick the 'furthest back' one (lowest t value).
-        let (i, v) = {
 
-            let mut t_min = t1;
+        let mut dt = (t1 - t0) / 100.0;
-            let mut i_min = 0;
+        let vf = frame[v.frame_idx];
-            let mut all_done = true;
-
-            for (i, v) in frontier.iter().enumerate() {
-                if v.traj.t < t_min {
-                    i_min = i;
-                    t_min = v.traj.t;
-                }
-                all_done &= v.traj.t >= t1;
-            }
-
-            // If no vertex can be advanced, we're done:
-            if all_done {
-                println!("All past t1");
-                break;
-            }
-            (i_min, frontier[i_min].clone())
-        };
-
-        // Indices (into 'frontier') of previous and next:
-        let i_prev = (i + n - 1) % n;
-        let i_next = (i + 1) % n;
-
-        println!("DEBUG: Moving frontier vertex {}, {:?} (t={}, frame_vert={:?})", i, v.traj.vert, v.traj.t, v.traj.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 vf = v.traj.frame_vert;
-        let vt = v.traj.t;
-        let dt_max = t1 - vt;
-        let mut dt = ((t1 - t0) / 100.0).min(dt_max);
         for iter in 0..100 {
             // Consider an edge from f(v.t)*vf to f(v.t + dt)*vf.
             // These two endpoints have zero error from the trajectory
@@ -487,292 +421,78 @@ pub fn parametric_mesh<F>(frame: Vec<Vertex>, 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(vt).mtx + f(vt + 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(vt + 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);
+            println!("DEBUG iter {}: dt={}, edge_mid={:?}, traj_mid={:?}, err={}", iter, dt, edge_mid, traj_mid, err);
 
             let r = (err - max_err).abs() / max_err;
             if r < 0.10 {
-                //println!("err close enough");
+                println!("err close enough");
-                println!("Error under threshold in {} iters, dt={}", iter, dt);
                 break;
             } else if err > max_err {
                 dt = dt / 2.0;
-                //println!("err > max_err, reducing dt to {}", dt);
+                println!("err > max_err, reducing dt to {}", dt);
             } else {
-                dt = (dt * 1.2).min(dt_max);
+                dt = dt * 1.2;
-                //println!("err < max_err, increasing dt to {}", dt);
+                println!("err < max_err, increasing dt to {}", dt);
             }
         }
-        // (note that this is just a crappy numerical approximation of
-        // dt given a desired error and it would probably be possible
-        // to do this directly given an analytical form of the
-        // curvature of f at some starting point)
 
-        let t = (vt + dt).min(t1);
+        let t = v.t + dt;
-        let v_next = VertexTrajectory {
+        let v_next = f(t).mtx * vf;
-            vert: f(t).mtx * vf,
-            t: t,
-            frame_vert: vf,
-        };
 
-        // DEBUG
+        // Add this vertex to our mesh:
-        /*
+        let pos = verts.len();
-        let l1 = (v.vert - v_next).norm();
+        verts.push(v_next);
-        let l2 = (frontier[v.neighbor1].vert - v.vert).norm();
+        // There are 3 other vertices of interest: the one we started
-        let l3 = (frontier[v.neighbor2].vert - v.vert).norm();
+        // from (v) and its two neighbors. We make two edges - one on
-        println!("l1={} l2={} l3={}", l1, l2, l3);
+        // each side of the edge (v, v_next).
-         */
-
-        // Add two faces to our mesh. They share two vertices, and thus
-        // the boundary in between those vertices.
-
-        //println!("DEBUG: Adding 1st face");
-
-        // First face: connect prior frontier vertex to 'v' & 'v_next'.
-        // f1 is that face,
-        // edge1 is: prior frontier vertex -> v_next.
-        // edge_v_next is: v_next -> v
-        let (f1, edge1, edge_v_next) = match v.halfedges[0] {
-            // However, the way we add the face depends on whether we are
-            // adding to an existing boundary or not:
-            None => {
-                let neighbor = &frontier[i_prev];
-                //println!("DEBUG: add_face()");
-                let (f1, edges) = mesh.add_face([
-                    VertSpec::New(v_next), // edges[0]: v_next -> v
-                    match v.vert_idx {
-                        None => VertSpec::New(v.traj),
-                        Some(idx) => VertSpec::Idx(idx),
-                    }, // edges[1]: v -> neighbor
-                    match neighbor.vert_idx {
-                        None => VertSpec::New(neighbor.traj),
-                        Some(idx) => VertSpec::Idx(idx),
-                    }, // edges[2]: neighbor -> v_next
-                ]);
-
-                if neighbor.vert_idx.is_none() {
-                    // If neighbor.vert_idx is None, then we had to
-                    // add its vertex to the mesh for the face we just
-                    // made - so mark it in the frontier:
-                    frontier[i_prev].vert_idx = Some(mesh.halfedges[edges[2]].vert);
-                }
-
-                (f1, edges[2], edges[0])
-            },
-            Some(edge_idx) => {
-                // edge_idx is half-edge from prior frontier vertex to 'v'
-                //println!("DEBUG: add_face_twin1({}, ({},{},{}))", edge_idx, v_next.x, v_next.y, v_next.z);
-                let (f1, edges) = mesh.add_face_twin1(edge_idx, v_next);
-                // edges[0] is twin of edge_idx, thus v -> neighbor
-                // edges[1] is neighbor -> v_next
-                // edges[2] is v_next -> v
-                (f1, edges[1], edges[2])
-            },
-        };
-        //println!("DEBUG: edge1={} edge_v_next={}", edge1, edge_v_next);
-
-        // DEBUG
-        //mesh.check();
-        //mesh.print();
-
-        //println!("DEBUG: Adding 2nd face");
-        // Second face: connect next frontier vertex to 'v' & 'v_next'.
-        // f2 is that face.
-        // edge2 is: v_next -> next frontier vertex
-        let (f2, edge2) = match v.halfedges[1] {
-            // Likewise, the way we add the second face depends on
-            // the same (but for the other side).  Regardless,
-            // they share the boundary between v_next and v.vert - which
-            // is edges1[0].
-            None => {
-                let neighbor = &frontier[i_next];
-
-                let (f2, edges) = match neighbor.vert_idx {
-                    None => {
-                        //let v = neighbor.traj.vert;
-                        //println!("DEBUG: add_face_twin1({}, ({},{},{}))", edge_v_next, v.x, v.y, v.z);
-                        let (f2, edges) = mesh.add_face_twin1(edge_v_next, neighbor.traj);
-                        // Reasoning here is identical to "If neighbor.vert_idx
-                        // is None..." above:
-                        frontier[i_next].vert_idx = Some(mesh.halfedges[edges[2]].vert);
-
-                        (f2, edges)
-                    },
-                    Some(vert_idx) => {
-                        //println!("DEBUG: add_face_twin1_ref({}, {})", edge_v_next, vert_idx);
-                        mesh.add_face_twin1_ref(edge_v_next, vert_idx)
-                    }
-                };
-                // For either branch:
-                // edges[0] is twin of edge_v_next, thus: v -> v_next
-                // edges[1] is: v_next -> neighbor
-                // egdes[2] is: neighbor -> v
-
-                (f2, edges[1])
-            },
-            Some(edge_idx) => {
-                //println!("DEBUG: add_face_twin2({}, {})", edge_idx, edge_v_next);
-                let (f2, edges) = mesh.add_face_twin2(edge_idx, edge_v_next);
-                // edges[0] is twin of edge_idx, thus: neighbor -> v
-                // edges[1] is twin of edge_v_next, thus: v -> v_next
-                // edges[2] is: v_next -> neighbor
-                //println!("DEBUG: add_face_twin2 returned {:?}", edges);
-                (f2, edges[2])
-            },
-        };
-        //println!("DEBUG: edge2={}", edge2);
-
-        // DEBUG
-        //println!("DEBUG: Added 2nd face");
-        //mesh.check();
-        //mesh.print();
-
-        // The 'shared' half-edge should start at v.vert - hence edges[1].
-
-        // and add two faces:
-        /*
         faces.append(&mut vec![
-            v_next_idx, v.mesh_idx, frontier[v.neighbor[0]].mesh_idx, // face_idx
+            v.mesh_idx, pos, frontier[v.neighbor1].mesh_idx,
-            v.mesh_idx, v_next_idx, frontier[v.neighbor[1]].mesh_idx, // face_idx + 1
+            pos, v.mesh_idx, frontier[v.neighbor2].mesh_idx,
         ]);
-         */
 
         // Replace this vertex in the frontier:
         frontier[i] = frontierVert {
-            traj: VertexTrajectory {
+            vert: v_next,
-                vert: v_next.vert,
+            frame_idx: v.frame_idx,
-                frame_vert: vf,
+            mesh_idx: pos,
-                t: t,
+            t: t,
-            },
+            neighbor1: v.neighbor1,
-            halfedges: [Some(edge1), Some(edge2)],
+            neighbor2: v.neighbor2,
-            vert_idx: Some(mesh.halfedges[edge_v_next].vert),
+        }
-        };
-        // Since 'halfedges' forms more or less a doubly-linked list, we
-        // must go to previous & next vertices and update them:
-        frontier[i_prev].halfedges[1] = Some(edge1);
-        frontier[i_next].halfedges[0] = Some(edge2);
     }
 
-    // A stack of face indices for faces in 'mesh' that are still
+    // Move this vertex further along, i.e. t + dt.  (dt is set by
-    // undergoing subdivision
+    // 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 stack: HashMap<usize, usize> = (0..mesh.num_faces).map(|i| (i, 0)).collect();
+    // Add this vertex to the mesh, and connect it to: the vertex we
+    // started with, and the two neighbors of that vertex.
 
-    let max_subdiv = 1;
+    // Repeat at "Pick a vertex...".
 
-    // This is masked off because it's just horrible:
+    // Don't move t + dt past t1.  Once a frontier vertex is placed at
-    while false && !stack.is_empty() {
+    // that value of t, remove it.
 
-        let (face, count) = match stack.iter().next() {
-            None => break,
-            Some((f, c)) => (*f, *c),
-        };
-        stack.remove(&face);
 
-        println!("DEBUG: Examining face: {:?} ({} subdivs)", face, count);
+    // 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.
 
-        if count >= max_subdiv {
+    // But still missing from that: When do I collapse a subdivision
-            continue;
+    // back down?
-        }
+    return Mesh { verts, faces };
-
-        let v_idx = mesh.face_to_verts(face);
-        if v_idx.len() != 3 {
-            panic!("Face {} has {} vertices, not 3?", face, v_idx.len());
-        }
-        let tr = [mesh.verts[v_idx[0]].v,
-                         mesh.verts[v_idx[1]].v,
-                         mesh.verts[v_idx[2]].v];
-        let (v0, v1, v2) = (tr[0].vert, tr[1].vert, tr[2].vert);
-        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 {
-            println!("DEBUG: Angles OK");
-        }
-        // TODO: Figure out how to subdivide in this case
-
-        // 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 = (tr[0].t + tr[1].t + tr[2].t) / 3.0;
-        let v_mid = (tr[0].frame_vert + tr[1].frame_vert + tr[2].frame_vert) / 3.0;
-        let p = f(t_mid).mtx * v_mid;
-
-        let d = p.xyz().dot(&normal);
-
-        if d.is_nan() {
-            println!("DEBUG: p={:?} normal={:?}", p, normal);
-            println!("DEBUG: a={:?} b={:?}", a, b);
-            println!("DEBUG: v0={:?} v1={:?} v2={:?}", v0, v1, v2);
-            //panic!("Distance is NaN?");
-            println!("DEBUG: distance is NaN?");
-            continue;
-        }
-
-        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);
-
-        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: Vec<VertexTrajectory> = pairs.iter().map(|(i,j)| {
-            let t = (tr[*i].t + tr[*j].t) / 2.0;
-            let v = (tr[*i].frame_vert + tr[*j].frame_vert) / 2.0;
-            VertexTrajectory {
-                vert: f(t).mtx * v,
-                frame_vert: v,
-                t: t,
-            }
-        }).collect();
-
-        // TODO: Get indices of added faces and put these on the stack too
-        match mesh.full_subdiv_face(face, mids) {
-            None => {},
-            Some((new, upd)) => {
-
-                // The 'updated' faces can remain on the stack if they
-                // already were there.  If they were taken off of the stack
-                // because they were too small, they are now even smaller,
-                // so it is fine to leave them off.
-                stack.extend(new.iter().map(|i| (*i, count + 1)));
-                stack.extend(upd.iter().map(|i| (*i, count + 1)));
-                // TODO: Why does the resultant mesh have holes if I do this
-                // and then increase max_subdiv?
-            },
-        }
-
-   }
-
-    //mesh.print();
-
-    return mesh.convert_mesh(|i| i.vert );
 }

src/util.rs

@@ -1,7 +1,6 @@
 use std::ops::Range;
-use crate::mesh::{Mesh, MeshFunc, VertexUnion};
+use crate::mesh::{MeshFunc, VertexUnion};
 use crate::xform::{Vertex};
-//use crate::rule::{Rule, Child};
 
 /// This is like `vec!`, but it can handle elements that are given
 /// with `@var element` rather than `element`, e.g. like