common: BvhNode::sample_ray

lib/common/src/bvh.rs

@@ -1,4 +1,10 @@
+use std::cmp::Ordering;
+use std::collections::BTreeMap;
+
 use crate::aabb::Aabb;
+use crate::ray::Ray;
+use crate::integer::{Ratio,Planar64};
+use crate::instruction::{InstructionCollector,TimedInstruction};
 
 //da algaritum
 //lista boxens
@@ -10,6 +16,96 @@ use crate::aabb::Aabb;
 //sort the centerpoints on each axis (3 lists)
 //bv is put into octant based on whether it is upper or lower in each list
 
+
+pub fn intersect_aabb(ray:&Ray,aabb:&Aabb)->Option<Ratio<Planar64,Planar64>>{
+	// n.(o+d*t)==n.p
+	// n.o + n.d * t == n.p
+	// t == (n.p - n.o)/n.d
+	let mut hit=None;
+	match ray.direction.x.cmp(&Planar64::ZERO){
+		Ordering::Less=>{
+			let rel_min=aabb.min()-ray.origin;
+			let rel_max=aabb.max()-ray.origin;
+			let dy=rel_max.x*ray.direction.y;
+			let dz=rel_max.x*ray.direction.z;
+			// x is negative, so inequalities are flipped
+			if rel_min.y*ray.direction.x>dy&&dy>rel_max.y*ray.direction.x
+			 &&rel_min.z*ray.direction.x>dz&&dz>rel_max.z*ray.direction.x{
+				let t=rel_max.x/ray.direction.x;
+				hit=Some(hit.map_or(t,|best_t|t.min(best_t)));
+			}
+		},
+		Ordering::Equal=>(),
+		Ordering::Greater=>{
+			let rel_min=aabb.min()-ray.origin;
+			let rel_max=aabb.max()-ray.origin;
+			let dy=rel_min.x*ray.direction.y;
+			let dz=rel_min.x*ray.direction.z;
+			// x is positive, so inequalities are normal
+			if rel_min.y*ray.direction.x<dy&&dy<rel_max.y*ray.direction.x
+			 &&rel_min.z*ray.direction.x<dz&&dz<rel_max.z*ray.direction.x{
+				let t=rel_min.x/ray.direction.x;
+				hit=Some(hit.map_or(t,|best_t|t.min(best_t)));
+			}
+		},
+	}
+	match ray.direction.z.cmp(&Planar64::ZERO){
+		Ordering::Less=>{
+			let rel_min=aabb.min()-ray.origin;
+			let rel_max=aabb.max()-ray.origin;
+			let dx=rel_max.z*ray.direction.x;
+			let dy=rel_max.z*ray.direction.y;
+			// z is negative, so inequalities are flipped
+			if rel_min.x*ray.direction.z>dx&&dx>rel_max.x*ray.direction.z
+			 &&rel_min.y*ray.direction.z>dy&&dy>rel_max.y*ray.direction.z{
+				let t=rel_max.z/ray.direction.z;
+				hit=Some(hit.map_or(t,|best_t|t.min(best_t)));
+			}
+		},
+		Ordering::Equal=>(),
+		Ordering::Greater=>{
+			let rel_min=aabb.min()-ray.origin;
+			let rel_max=aabb.max()-ray.origin;
+			let dx=rel_min.z*ray.direction.x;
+			let dy=rel_min.z*ray.direction.y;
+			// z is positive, so inequalities are normal
+			if rel_min.x*ray.direction.z<dx&&dx<rel_max.x*ray.direction.z
+			 &&rel_min.y*ray.direction.z<dy&&dy<rel_max.y*ray.direction.z{
+				let t=rel_min.z/ray.direction.z;
+				hit=Some(hit.map_or(t,|best_t|t.min(best_t)));
+			}
+		},
+	}
+	match ray.direction.y.cmp(&Planar64::ZERO){
+		Ordering::Less=>{
+			let rel_min=aabb.min()-ray.origin;
+			let rel_max=aabb.max()-ray.origin;
+			let dz=rel_max.y*ray.direction.z;
+			let dx=rel_max.y*ray.direction.x;
+			// y is negative, so inequalities are flipped
+			if rel_min.z*ray.direction.y>dz&&dz>rel_max.z*ray.direction.y
+			 &&rel_min.x*ray.direction.y>dx&&dx>rel_max.x*ray.direction.y{
+				let t=rel_max.y/ray.direction.y;
+				hit=Some(hit.map_or(t,|best_t|t.min(best_t)));
+			}
+		},
+		Ordering::Equal=>(),
+		Ordering::Greater=>{
+			let rel_min=aabb.min()-ray.origin;
+			let rel_max=aabb.max()-ray.origin;
+			let dz=rel_min.y*ray.direction.z;
+			let dx=rel_min.y*ray.direction.x;
+			// y is positive, so inequalities are normal
+			if rel_min.z*ray.direction.y<dz&&dz<rel_max.z*ray.direction.y
+			 &&rel_min.x*ray.direction.y<dx&&dx<rel_max.x*ray.direction.y{
+				let t=rel_min.y/ray.direction.y;
+				hit=Some(hit.map_or(t,|best_t|t.min(best_t)));
+			}
+		},
+	}
+	hit
+}
+
 pub enum RecursiveContent<N,L>{
 	Branch(Vec<N>),
 	Leaf(L),
@@ -44,6 +140,92 @@ impl<L> BvhNode<L>{
 			},
 		}
 	}
+	fn populate_nodes<'a,T,F>(
+		&'a self,
+		collector:&mut InstructionCollector<&'a L,Ratio<Planar64,Planar64>>,
+		nodes:&mut BTreeMap<Ratio<Planar64,Planar64>,&'a BvhNode<L>>,
+		ray:&Ray,
+		start_time:Ratio<Planar64,Planar64>,
+		f:&F,
+	)
+		where
+			T:Ord+Copy,
+			Ratio<Planar64,Planar64>:From<T>,
+			F:Fn(&L,&Ray)->Option<T>,
+	{
+		match &self.content{
+			RecursiveContent::Leaf(leaf)=>if let Some(time)=f(leaf,ray){
+				let ins=TimedInstruction{time:time.into(),instruction:leaf};
+				if start_time.lt_ratio(ins.time)&&ins.time.lt_ratio(collector.time()){
+					collector.collect(Some(ins));
+				}
+			},
+			RecursiveContent::Branch(children)=>for child in children{
+				if child.aabb.contains(ray.origin){
+					child.populate_nodes(collector,nodes,ray,start_time,f);
+				}else{
+					// Am I an upcoming superstar?
+					if let Some(t)=intersect_aabb(ray,&child.aabb){
+						if start_time.lt_ratio(t)&&t.lt_ratio(collector.time()){
+							nodes.insert(t,child);
+						}
+					}
+				}
+			},
+		}
+	}
+	pub fn sample_ray<T,F>(
+		&self,
+		ray:&Ray,
+		start_time:T,
+		time_limit:T,
+		f:F,
+	)->Option<(T,&L)>
+		where
+			T:Ord+Copy,
+			T:From<Ratio<Planar64,Planar64>>,
+			Ratio<Planar64,Planar64>:From<T>,
+			F:Fn(&L,&Ray)->Option<T>,
+	{
+		// source of nondeterminism when Aabb boundaries are coplanar
+		let mut nodes=BTreeMap::new();
+
+		let start_time=start_time.into();
+		let time_limit=time_limit.into();
+		let mut collector=InstructionCollector::new(time_limit);
+		// break open all nodes that contain ray.origin and populate nodes with future intersection times
+		self.populate_nodes(&mut collector,&mut nodes,ray,start_time,&f);
+
+		// swim through nodes one at a time
+		while let Some((t,node))=nodes.pop_first(){
+			if collector.time()<t{
+				break;
+			}
+			match &node.content{
+				RecursiveContent::Leaf(leaf)=>if let Some(time)=f(leaf,ray){
+					let ins=TimedInstruction{time:time.into(),instruction:leaf};
+					// this lower bound can also be omitted
+					// but it causes type inference errors lol
+					if start_time.lt_ratio(ins.time)&&ins.time.lt_ratio(collector.time()){
+						collector.collect(Some(ins));
+					}
+				},
+				// break open the node and predict collisions with the child nodes
+				RecursiveContent::Branch(children)=>for child in children{
+					// Am I an upcoming superstar?
+					if let Some(t)=intersect_aabb(ray,&child.aabb){
+						// we don't need to check the lower bound
+						// because child aabbs are guaranteed to be within the parent bounds.
+						if t<collector.time(){
+							nodes.insert(t,child);
+						}
+					}
+				},
+			}
+		}
+
+		collector.take().map(|TimedInstruction{time,instruction:leaf}|(time.into(),leaf))
+	}
 	pub fn into_inner(self)->(RecursiveContent<BvhNode<L>,L>,Aabb){
 		(self.content,self.aabb)
 	}