use strafesnet_common::integer::vec3; use strafesnet_common::integer::vec3::Vector3; use strafesnet_common::integer::{Fixed,Planar64Vec3}; use crate::model::{DirectedEdge,FEV,MeshQuery,MinkowskiMesh,MinkowskiVert}; // This algorithm is based on Lua code // written by Trey Reynolds in 2021 type Simplex=[MinkowskiVert;N]; #[derive(Clone,Copy)] enum Simplex1_3{ Simplex1(Simplex<1>), Simplex2(Simplex<2>), Simplex3(Simplex<3>), } impl Simplex1_3{ fn push_front(self,v:MinkowskiVert)->Simplex2_4{ match self{ Simplex1_3::Simplex1([v0])=>Simplex2_4::Simplex2([v,v0]), Simplex1_3::Simplex2([v0,v1])=>Simplex2_4::Simplex3([v,v0,v1]), Simplex1_3::Simplex3([v0,v1,v2])=>Simplex2_4::Simplex4([v,v0,v1,v2]), } } } enum Simplex2_4{ Simplex2(Simplex<2>), Simplex3(Simplex<3>), Simplex4(Simplex<4>), } /* local function absDet(r, u, v, w) if w then return math.abs((u - r):Cross(v - r):Dot(w - r)) elseif v then return (u - r):Cross(v - r).magnitude elseif u then return (u - r).magnitude else return 1 end end */ impl Simplex2_4{ fn det_is_zero(&self,mesh:&MinkowskiMesh)->bool{ match self{ &Self::Simplex4([p0,p1,p2,p3])=>{ let p0=mesh.vert(p0); let p1=mesh.vert(p1); let p2=mesh.vert(p2); let p3=mesh.vert(p3); (p1-p0).cross(p2-p0).dot(p3-p0)==Fixed::ZERO }, &Self::Simplex3([p0,p1,p2])=>{ let p0=mesh.vert(p0); let p1=mesh.vert(p1); let p2=mesh.vert(p2); (p1-p0).cross(p2-p0)==vec3::zero() }, &Self::Simplex2([p0,p1])=>{ let p0=mesh.vert(p0); let p1=mesh.vert(p1); p1-p0==vec3::zero() } } } } /* local function choosePerpendicularDirection(d) local x, y, z = d.x, d.y, d.z local best = math.min(x*x, y*y, z*z) if x*x == best then return Vector3.new(y*y + z*z, -x*y, -x*z) elseif y*y == best then return Vector3.new(-x*y, x*x + z*z, -y*z) else return Vector3.new(-x*z, -y*z, x*x + y*y) end end */ fn choose_perpendicular_direction(d:Planar64Vec3)->Planar64Vec3{ let x=d.x.abs(); let y=d.y.abs(); let z=d.z.abs(); if xPlanar64Vec3{ vec3::X } fn reduce1( [v0]:Simplex<1>, mesh:&MinkowskiMesh, point:Planar64Vec3, )->Reduced{ // --debug.profilebegin("reduceSimplex0") // local a = a1 - a0 let p0=mesh.vert(v0); // local p = -a let p=-(p0+point); // local direction = p let mut dir=p; // if direction.magnitude == 0 then // direction = chooseAnyDirection() if dir==vec3::zero(){ dir=choose_any_direction(); } // return direction, a0, a1 Reduced{ dir, simplex:Simplex1_3::Simplex1([v0]), } } // local function reduceSimplex1(a0, a1, b0, b1) fn reduce2( [v0,v1]:Simplex<2>, mesh:&MinkowskiMesh, point:Planar64Vec3, )->Reduced{ // --debug.profilebegin("reduceSimplex1") // local a = a1 - a0 // local b = b1 - b0 let p0=mesh.vert(v0); let p1=mesh.vert(v1); // local p = -a // local u = b - a let p=-(p0+point); let u=p1-p0; // -- modify to take into account the radiuses // local p_u = p:Dot(u) let p_u=p.dot(u); // if p_u >= 0 then if !p_u.is_negative(){ // local direction = u:Cross(p):Cross(u) let direction=u.cross(p).cross(u); // if direction.magnitude == 0 then if direction==vec3::zero(){ return Reduced{ dir:choose_perpendicular_direction(u), simplex:Simplex1_3::Simplex2([v0,v1]), }; } // -- modify the direction to take into account a0R and b0R // return direction, a0, a1, b0, b1 return Reduced{ dir:direction.narrow_1().unwrap(), simplex:Simplex1_3::Simplex2([v0,v1]), }; } // local direction = p let mut dir=p; // if direction.magnitude == 0 then if dir==vec3::zero(){ dir=choose_perpendicular_direction(u); } // return direction, a0, a1 Reduced{ dir, simplex:Simplex1_3::Simplex1([v0]), } } // local function reduceSimplex2(a0, a1, b0, b1, c0, c1) fn reduce3( [v0,mut v1,v2]:Simplex<3>, mesh:&MinkowskiMesh, point:Planar64Vec3, )->Reduced{ // --debug.profilebegin("reduceSimplex2") // local a = a1 - a0 // local b = b1 - b0 // local c = c1 - c0 let p0=mesh.vert(v0); let p1=mesh.vert(v1); let p2=mesh.vert(v2); // local p = -a // local u = b - a // local v = c - a let p=-(p0+point); let mut u=p1-p0; let v=p2-p0; // local uv = u:Cross(v) // local up = u:Cross(p) // local pv = p:Cross(v) // local uv_up = uv:Dot(up) // local uv_pv = uv:Dot(pv) let mut uv=u.cross(v); let mut up=u.cross(p); let pv=p.cross(v); let uv_up=uv.dot(up); let uv_pv=uv.dot(pv); // if uv_up >= 0 and uv_pv >= 0 then if !uv_up.is_negative()&&!uv_pv.is_negative(){ // local uvp = uv:Dot(p) let uvp=uv.dot(p); // local direction = uvp < 0 and -uv or uv let direction=if uvp.is_negative(){ -uv }else{ uv }; // return direction, a0, a1, b0, b1, c0, c1 return Reduced{ dir:direction.narrow_1().unwrap(), simplex:Simplex1_3::Simplex3([v0,v1,v2]), }; } // local u_u = u:Dot(u) // local v_v = v:Dot(v) // local uDist = uv_up/(u_u*v.magnitude) // local vDist = uv_pv/(v_v*u.magnitude) // local minDist2 = math.min(uDist, vDist) let u_dist=uv_up*v.length(); let v_dist=uv_pv*u.length(); // if vDist == minDist2 then if v_dist= 0 then if !p_u.is_negative(){ // local direction = up:Cross(u) let direction=up.cross(u); // if direction.magnitude == 0 then if direction==vec3::zero(){ // direction = uv return Reduced{ dir:uv.narrow_1().unwrap(), simplex:Simplex1_3::Simplex2([v0,v1]), }; } // return direction, a0, a1, b0, b1 return Reduced{ dir:direction.narrow_1().unwrap(), simplex:Simplex1_3::Simplex2([v0,v1]), }; } // local direction = p let dir=p; // if direction.magnitude == 0 then if dir==vec3::zero(){ // direction = uv return Reduced{ dir:uv.narrow_1().unwrap(), simplex:Simplex1_3::Simplex1([v0]), }; } // return direction, a0, a0 Reduced{ dir, simplex:Simplex1_3::Simplex1([v0]), } } // local function reduceSimplex3(a0, a1, b0, b1, c0, c1, d0, d1) fn reduce4( [v0,mut v1,mut v2,v3]:Simplex<4>, mesh:&MinkowskiMesh, point:Planar64Vec3, )->Reduce{ // --debug.profilebegin("reduceSimplex3") // local a = a1 - a0 // local b = b1 - b0 // local c = c1 - c0 // local d = d1 - d0 let p0=mesh.vert(v0); let p1=mesh.vert(v1); let p2=mesh.vert(v2); let p3=mesh.vert(v3); // local p = -a // local u = b - a // local v = c - a // local w = d - a let p=-(p0+point); let mut u=p1-p0; let mut v=p2-p0; let w=p3-p0; // local uv = u:Cross(v) // local vw = v:Cross(w) // local wu = w:Cross(u) // local uvw = uv:Dot(w) // local pvw = vw:Dot(p) // local upw = wu:Dot(p) // local uvp = uv:Dot(p) let mut uv=u.cross(v); let vw=v.cross(w); let wu=w.cross(u); let uv_w=uv.dot(w); let pv_w=vw.dot(p); let up_w=wu.dot(p); let uv_p=uv.dot(p); // if pvw/uvw >= 0 and upw/uvw >= 0 and uvp/uvw >= 0 then if !pv_w.div_sign(uv_w).is_negative() ||!up_w.div_sign(uv_w).is_negative() ||!uv_p.div_sign(uv_w).is_negative(){ // origin is contained, this is a positive detection // local direction = Vector3.new(0, 0, 0) // return direction, a0, a1, b0, b1, c0, c1, d0, d1 return Reduce::Escape([v0,v1,v2,v3]); } // local uvwSign = uvw < 0 and -1 or uvw > 0 and 1 or 0 // local uvDist = uvp*uvwSign/uv.magnitude // local vwDist = pvw*uvwSign/vw.magnitude // local wuDist = upw*uvwSign/wu.magnitude // local minDist3 = math.min(uvDist, vwDist, wuDist) let uv_dist=uv_p.mul_sign(uv_w); let vw_dist=pv_w.mul_sign(uv_w); let wu_dist=up_w.mul_sign(uv_w); let wu_len=wu.length(); let uv_len=uv.length(); let vw_len=vw.length(); if vw_dist*wu_len= 0 and uv_pv >= 0 then if !uv_up.is_negative()&&!uv_pv.is_negative(){ // local direction = uvw < 0 and uv or -uv // return direction, a0, a1, b0, b1, c0, c1 if uv_w.is_negative(){ return Reduce::Reduced(Reduced{ dir:uv.narrow_1().unwrap(), simplex:Simplex1_3::Simplex3([v0,v1,v2]), }); }else{ return Reduce::Reduced(Reduced{ dir:-uv.narrow_1().unwrap(), simplex:Simplex1_3::Simplex3([v0,v1,v2]), }); } } // local u_u = u:Dot(u) // local v_v = v:Dot(v) // local uDist = uv_up/(u_u*v.magnitude) // local vDist = uv_pv/(v_v*u.magnitude) // local minDist2 = math.min(uDist, vDist) let u_dist=uv_up*v.length(); let v_dist=uv_pv*u.length(); // if vDist == minDist2 then if v_dist= 0 then if !p_u.is_negative(){ // local direction = up:Cross(u) let direction=up.cross(u); // if direction.magnitude == 0 then if direction==vec3::zero(){ // direction = uvw < 0 and uv or -uv // return direction, a0, a1, b0, b1 if uv_w.is_negative(){ return Reduce::Reduced(Reduced{ dir:uv.narrow_1().unwrap(), simplex:Simplex1_3::Simplex2([v0,v1]), }); }else{ return Reduce::Reduced(Reduced{ dir:-uv.narrow_1().unwrap(), simplex:Simplex1_3::Simplex2([v0,v1]), }); } } // return direction, a0, a1, b0, b1 return Reduce::Reduced(Reduced{ dir:direction.narrow_1().unwrap(), simplex:Simplex1_3::Simplex2([v0,v1]), }); } // local direction = p let dir=p; // if direction.magnitude == 0 then if dir==vec3::zero(){ // direction = uvw < 0 and uv or -uv if uv_w.is_negative(){ return Reduce::Reduced(Reduced{ dir:uv.narrow_1().unwrap(), simplex:Simplex1_3::Simplex1([v0]), }); }else{ return Reduce::Reduced(Reduced{ dir:-uv.narrow_1().unwrap(), simplex:Simplex1_3::Simplex1([v0]), }); } } // return direction, a0, a1 Reduce::Reduced(Reduced{ dir, simplex:Simplex1_3::Simplex1([v0]), }) } struct Reduced{ dir:Planar64Vec3, simplex:Simplex1_3, } enum Reduce{ Escape(Simplex<4>), Reduced(Reduced), } impl Simplex2_4{ fn reduce(self,mesh:&MinkowskiMesh,point:Planar64Vec3)->Reduce{ match self{ Self::Simplex2(simplex)=>Reduce::Reduced(reduce2(simplex,mesh,point)), Self::Simplex3(simplex)=>Reduce::Reduced(reduce3(simplex,mesh,point)), Self::Simplex4(simplex)=>reduce4(simplex,mesh,point), } } } pub fn contains_point(mesh:&MinkowskiMesh,point:Planar64Vec3)->bool{ const ENABLE_FAST_FAIL:bool=true; // TODO: remove mesh negation minimum_difference::(&-mesh,point, // on_exact |is_intersecting,_simplex|{ is_intersecting }, // on_escape |_simplex|{ // intersection is guaranteed at this point true }, // fast_fail value ||false ) } //infinity fev algorithm state transition #[derive(Debug)] enum Transition{ Done,//found closest vert, no edges are better Vert(MinkowskiVert),//transition to vert } enum EV{ Vert(MinkowskiVert), Edge(crate::model::MinkowskiEdge), } impl MinkowskiMesh<'_>{ fn next_transition_vert(&self,vert_id:MinkowskiVert,best_distance_squared:&mut Fixed<2,64>,point:Planar64Vec3)->Transition{ let mut best_transition=Transition::Done; for &directed_edge_id in self.vert_edges(vert_id).as_ref(){ //is boundary uncrossable by a crawl from infinity let edge_verts=self.edge_verts(directed_edge_id.as_undirected()); //select opposite vertex let test_vert_id=edge_verts.as_ref()[directed_edge_id.parity() as usize]; //test if it's closer let diff=point-self.vert(test_vert_id); let distance_squared=diff.dot(diff); if distance_squared<*best_distance_squared{ best_transition=Transition::Vert(test_vert_id); *best_distance_squared=distance_squared; } } best_transition } fn final_ev(&self,vert_id:MinkowskiVert,best_distance_squared:&mut Fixed<2,64>,point:Planar64Vec3)->EV{ let mut best_transition=EV::Vert(vert_id); let diff=point-self.vert(vert_id); for &directed_edge_id in self.vert_edges(vert_id).as_ref(){ let edge_n=self.directed_edge_n(directed_edge_id); //is boundary uncrossable by a crawl from infinity //check if time of collision is outside Time::MIN..Time::MAX let d=edge_n.dot(diff); //test the edge let edge_nn=edge_n.dot(edge_n); if !d.is_negative()&&d<=edge_nn{ let distance_squared={ let c=diff.cross(edge_n); //wrap for speed (c.dot(c)/edge_nn).divide().wrap_2() }; if distance_squared<=*best_distance_squared{ best_transition=EV::Edge(directed_edge_id.as_undirected()); *best_distance_squared=distance_squared; } } } best_transition } fn crawl_boundaries(&self,mut vert_id:MinkowskiVert,mut best_distance_squared:Fixed<2,64>,point:Planar64Vec3)->EV{ loop{ match self.next_transition_vert(vert_id,&mut best_distance_squared,point){ Transition::Done=>return self.final_ev(vert_id,&mut best_distance_squared,point), Transition::Vert(new_vert_id)=>vert_id=new_vert_id, } } } /// This function drops a vertex down to an edge or a face if the path from infinity did not cross any vertex-edge boundaries but the point is supposed to have already crossed a boundary down from a vertex fn crawl_to_closest_fev(&self,simplex:Simplex<3>,point:Planar64Vec3)->FEV::{ let (vert_id,best_distance_squared)=simplex.into_iter().map(|vert_id|{ let diff=point-self.vert(vert_id); (vert_id,diff.dot(diff)) }).min_by_key(|&(_,d)|d).unwrap(); //start on any vertex //cross uncrossable vertex-edge boundaries until you find the closest vertex or edge //cross edge-face boundary if it's uncrossable match self.crawl_boundaries(vert_id,best_distance_squared,point){ //if a vert is returned, it is the closest point to the infinity point EV::Vert(vert_id)=>FEV::Vert(vert_id), EV::Edge(edge_id)=>{ //cross to face if the boundary is not crossable and we are on the wrong side let edge_n=self.edge_n(edge_id); // point is multiplied by two because vert_sum sums two vertices. let delta_pos=point*2-{ let &[v0,v1]=self.edge_verts(edge_id).as_ref(); self.vert(v0)+self.vert(v1) }; for (i,&face_id) in self.edge_faces(edge_id).as_ref().iter().enumerate(){ let face_n=self.face_nd(face_id).0; //edge-face boundary nd, n facing out of the face towards the edge let boundary_n=face_n.cross(edge_n)*(i as i64*2-1); let boundary_d=boundary_n.dot(delta_pos); //check if time of collision is outside Time::MIN..Time::MAX //infinity_dir can always be treated as a velocity if !boundary_d.is_positive(){ //both faces cannot pass this condition, return early if one does. return FEV::Face(face_id); } } FEV::Edge(edge_id) }, } } } #[derive(Debug)] pub struct OhNoes; pub fn closest_fev_not_inside<'a>(mesh:&MinkowskiMesh<'a>,point:Planar64Vec3)->Option>,OhNoes>>{ const ENABLE_FAST_FAIL:bool=false; // TODO: remove mesh negation minimum_difference::(&-mesh,point, // on_exact |is_intersecting,simplex|{ if is_intersecting{ return None; } // Convert simplex to FEV // Vertices must be inverted since the mesh is inverted Some(match simplex{ Simplex1_3::Simplex1([v0])=>Ok(FEV::Vert(-v0)), Simplex1_3::Simplex2([v0,v1])=>{ // invert let (v0,v1)=(-v0,-v1); // TODO: handle non-canonnical multi-edge spanning edges // dumbest stupidest brute force search let v0e=mesh.vert_edges(v0); for &v0e in v0e.as_ref(){ // check opposite vertex to see if it is v1 if mesh.edge_verts(v0e.as_undirected()).as_ref()[v0e.parity() as usize]==v1{ return Some(Ok(FEV::Edge(v0e.as_undirected()))); } } Err(OhNoes) }, Simplex1_3::Simplex3([v0,v1,v2])=>{ // invert let (v0,v1,v2)=(-v0,-v1,-v2); let p0=mesh.vert(v0); let p1=mesh.vert(v1); let p2=mesh.vert(v2); // spanning simplex normal let n=(p2-p0).cross(p1-p0); // Scan opposite vertices for coplanar ones and find a coplanar face for &v0e in mesh.vert_edges(v0).as_ref(){ // check if opposite vertex is coplanar let o0=mesh.edge_verts(v0e.as_undirected()).as_ref()[v0e.parity() as usize]; let o0p=mesh.vert(o0); if n.dot(o0p-p0).is_zero(){ // always take left face (or right idk just use the parity) let f0=mesh.edge_faces(v0e.as_undirected()).as_ref()[v0e.parity() as usize]; { } } } // Shimmy to the side until you find a face that contains the closest point Err(OhNoes) }, }) }, // on_escape |_simplex|{ // intersection is guaranteed at this point // local norm, dist, u0, u1, v0, v1, w0, w1 = expand(queryP, queryQ, a0, a1, b0, b1, c0, c1, d0, d1, 1e-5) // let simplex=refine_to_exact(mesh,simplex); None }, // fast_fail value is irrelevant and will never be returned! ||unreachable!() ) } // local function minimumDifference( // queryP, radiusP, // queryQ, radiusQ, // exitRadius, testIntersection // ) fn minimum_difference( mesh:&MinkowskiMesh, point:Planar64Vec3, on_exact:impl FnOnce(bool,Simplex1_3)->T, on_escape:impl FnOnce(Simplex<4>)->T, on_fast_fail:impl FnOnce()->T, )->T{ // local initialAxis = queryQ() - queryP() // local new_point_p = queryP(initialAxis) // local new_point_q = queryQ(-initialAxis) // local direction, a0, a1, b0, b1, c0, c1, d0, d1 let mut initial_axis=mesh.hint_point()+point; // degenerate case if initial_axis==vec3::zero(){ initial_axis=choose_any_direction(); } let last_point=mesh.farthest_vert(-initial_axis); // this represents the 'a' value in the commented code let mut last_pos=mesh.vert(last_point); let Reduced{dir:mut direction,simplex:mut simplex_small}=reduce1([last_point],mesh,point); // exitRadius = testIntersection and 0 or exitRadius or 1/0 // for _ = 1, 100 do loop{ // new_point_p = queryP(-direction) // new_point_q = queryQ(direction) // local next_point = new_point_q - new_point_p let next_point=mesh.farthest_vert(direction); let next_pos=mesh.vert(next_point); // if -direction:Dot(next_point) > (exitRadius + radiusP + radiusQ)*direction.magnitude then if ENABLE_FAST_FAIL&&direction.dot(next_pos+point).is_negative(){ return on_fast_fail(); } let simplex_big=simplex_small.push_front(next_point); // if // direction:Dot(next_point - a) <= 0 or // absDet(next_point, a, b, c) < 1e-6 if !direction.dot(next_pos-last_pos).is_positive() ||simplex_big.det_is_zero(mesh){ // Found enough information to compute the exact closest point. // local norm = direction.unit // local dist = a:Dot(norm) // local hits = -dist < radiusP + radiusQ let is_intersecting=(last_pos+point).dot(direction).is_positive(); return on_exact(is_intersecting,simplex_small); } // direction, a0, a1, b0, b1, c0, c1, d0, d1 = reduceSimplex(new_point_p, new_point_q, a0, a1, b0, b1, c0, c1) match simplex_big.reduce(mesh,point){ // if a and b and c and d then Reduce::Escape(simplex)=>{ // Enough information to conclude that the meshes are intersecting. // Topology information is computed if needed. return on_escape(simplex); }, Reduce::Reduced(reduced)=>{ direction=reduced.dir; simplex_small=reduced.simplex; }, } // next loop this will be a last_pos=next_pos; } }