From 6b8a7b8bc69dda2dab39226dc3829ce3cd1d1df5 Mon Sep 17 00:00:00 2001 From: Chris Hodapp Date: Thu, 1 Oct 2020 12:46:43 -0400 Subject: [PATCH] Half-working pyramid example & roll back DCEL crap --- README.md | 3 + src/examples.rs | 103 ++++++++++-- src/lib.rs | 5 +- src/mesh.rs | 2 +- src/rule.rs | 438 +++++++++--------------------------------------- src/util.rs | 3 +- 6 files changed, 172 insertions(+), 382 deletions(-) diff --git a/README.md b/README.md index b3589e4..cb091e0 100644 --- a/README.md +++ b/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. diff --git a/src/examples.rs b/src/examples.rs index 3382e65..e8b6694 100644 --- a/src/examples.rs +++ b/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()), - final_geom: Rc::new(prim::empty_mesh()), + geom: Rc::new(prim::cube().to_meshfunc()), + 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 = 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)); - // TODO: Cleanliness fix - why not just make these return meshes? - let (vc, faces) = util::connect_convex(&seed_next, true); - let final_geom = Rc::new(OpenMesh { - verts: vec![vc], - alias_verts: vec![], - faces: faces, - }); + //let vc = util::centroid(&seed_next); + //let faces = util::connect_convex(0..n, n, true); + let geom = util::parallel_zigzag(seed_next, 0..n, 0..n); + let final_geom = MeshFunc { + verts: vec![], + faces: vec![], + // TODO: get actual verts here + }; let c = move |self_: Rc>| -> 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 { }, } } - +*/ #[derive(Copy, Clone)] pub struct TorusCtxt { @@ -389,6 +451,7 @@ pub struct TorusCtxt { stack: [Transform; 3], } +/* pub fn twisty_torus() -> Rule { let subdiv = 8; let seed = vec![ @@ -399,8 +462,10 @@ pub fn twisty_torus() -> Rule { ]; 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 { }, } } + */ +/* 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). diff --git a/src/lib.rs b/src/lib.rs index 7116043..87fc8f9 100644 --- a/src/lib.rs +++ b/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() { diff --git a/src/mesh.rs b/src/mesh.rs index 79d4336..561b28b 100644 --- a/src/mesh.rs +++ b/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(), diff --git a/src/rule.rs b/src/rule.rs index 504bd9a..dc89498 100644 --- a/src/rule.rs +++ b/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 = Rc>) -> RuleEval>; @@ -347,6 +344,8 @@ impl RuleEval { } +// 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(frame: Vec, 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 = DCELMesh::new(); - - #[derive(Clone, Debug)] struct frontierVert { - traj: VertexTrajectory, - // If the boundaries on either side of this vertex lie on a face - // (which is the case for all vertices *except* the initial ones), - // then this gives the halfedges of those boundaries. halfedges[0] - // connects the 'prior' vertex on the frontier to this, and - // halfedges[1] connect this to the 'next' vertex on the fronter. - // (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; 2], - // If this vertex is already in 'mesh', its vertex index there: - vert_idx: Option, + vert: Vertex, // Vertex position + t: f32, // Parameter value; f(t) should equal vert + frame_idx: usize, // Index of 'frame' this sits in the trajectory of + mesh_idx: usize, // Index of this vertex in the mesh + neighbor1: usize, // Index of 'frontier' of one neighbor + neighbor2: usize, // Index of 'frontier' of other neighbor }; // Init 'frontier' with each 'frame' vertex, and start it at t=t0. let mut frontier: Vec = frame.iter().enumerate().map(|(i,v)| frontierVert { - traj: VertexTrajectory { - vert: *v, - t: t0, - frame_vert: *v, - }, - halfedges: [None; 2], - vert_idx: None, + vert: *v, + t: t0, + frame_idx: i, + mesh_idx: i, + neighbor1: (i - 1) % n, + neighbor2: (i + 1) % n, }).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(frame: Vec, 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 = frame.clone(); + let mut faces: Vec = vec![]; + while !frontier.is_empty() { - /* - for (i, f) in frontier.iter().enumerate() { - 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 { - Some(vert) => { - match f.halfedges[1] { - Some(idx) => { - if mesh.halfedges[idx].vert != vert { - println!(" Error: halfedge[1]={} starts at {}, should start at {}", idx, mesh.halfedges[idx].vert, vert); - } - }, - 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. + // + // Heuristic for now: pick the 'furthest back' (lowest t) + 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(); + // TODO: Make this less ugly? + + if v.t >= t1 { + break; } - */ - // Pick a vertex to advance. For now, this is an easy - // heuristic: pick the 'furthest back' one (lowest t value). - let (i, v) = { + println!("DEBUG: Moving vertex {}, {:?} (t={}, frame_idx={})", i, v.vert, v.t, v.frame_idx); - let mut t_min = t1; - let mut i_min = 0; - 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); + let mut dt = (t1 - t0) / 100.0; + let vf = frame[v.frame_idx]; 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(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(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!("Error under threshold in {} iters, dt={}", iter, dt); + println!("err close enough"); 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); - //println!("err < max_err, increasing dt to {}", dt); + dt = dt * 1.2; + 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 v_next = VertexTrajectory { - vert: f(t).mtx * vf, - t: t, - frame_vert: vf, - }; + 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 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: - /* + // Add this vertex to our mesh: + let pos = verts.len(); + verts.push(v_next); + // There are 3 other vertices of interest: the one we started + // from (v) and its two neighbors. We make two edges - one on + // each side of the edge (v, v_next). 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 + v.mesh_idx, pos, frontier[v.neighbor1].mesh_idx, + pos, v.mesh_idx, frontier[v.neighbor2].mesh_idx, ]); - */ // Replace this vertex in the frontier: frontier[i] = frontierVert { - traj: VertexTrajectory { - vert: v_next.vert, - frame_vert: vf, - t: t, - }, - halfedges: [Some(edge1), Some(edge2)], - 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); + vert: v_next, + frame_idx: v.frame_idx, + mesh_idx: pos, + t: t, + neighbor1: v.neighbor1, + neighbor2: v.neighbor2, + } } - // A stack of face indices for faces in 'mesh' that are still - // undergoing subdivision + // 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 stack: HashMap = (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: - while false && !stack.is_empty() { + // Don't move t + dt past t1. Once a frontier vertex is placed at + // 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 { - continue; - } - - 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 = 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 ); + // But still missing from that: When do I collapse a subdivision + // back down? + return Mesh { verts, faces }; } diff --git a/src/util.rs b/src/util.rs index 2d2ee2c..9e16adb 100644 --- a/src/util.rs +++ b/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