333 lines
10 KiB
Rust
333 lines
10 KiB
Rust
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
|
|
//sort by {minx,maxx,miny,maxy,minz,maxz} (6 lists)
|
|
//find the sets that minimizes the sum of surface areas
|
|
//splitting is done when the minimum split sum of surface areas is larger than the node's own surface area
|
|
|
|
//start with bisection into octrees because a bad bvh is still 1000x better than no bvh
|
|
//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),
|
|
}
|
|
impl<N,L> RecursiveContent<N,L>{
|
|
pub fn empty()->Self{
|
|
Self::Branch(Vec::new())
|
|
}
|
|
}
|
|
pub struct BvhNode<L>{
|
|
content:RecursiveContent<BvhNode<L>,L>,
|
|
aabb:Aabb,
|
|
}
|
|
impl<L> BvhNode<L>{
|
|
pub fn empty()->Self{
|
|
Self{
|
|
content:RecursiveContent::empty(),
|
|
aabb:Aabb::default(),
|
|
}
|
|
}
|
|
pub fn sample_aabb<F:FnMut(&L)>(&self,aabb:&Aabb,f:&mut F){
|
|
match &self.content{
|
|
RecursiveContent::Leaf(model)=>f(model),
|
|
RecursiveContent::Branch(children)=>for child in children{
|
|
//this test could be moved outside the match statement
|
|
//but that would test the root node aabb
|
|
//you're probably not going to spend a lot of time outside the map,
|
|
//so the test is extra work for nothing
|
|
if aabb.intersects(&child.aabb){
|
|
child.sample_aabb(aabb,f);
|
|
}
|
|
},
|
|
}
|
|
}
|
|
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)
|
|
}
|
|
pub fn into_visitor<F:FnMut(L)>(self,f:&mut F){
|
|
match self.content{
|
|
RecursiveContent::Leaf(model)=>f(model),
|
|
RecursiveContent::Branch(children)=>for child in children{
|
|
child.into_visitor(f)
|
|
},
|
|
}
|
|
}
|
|
}
|
|
|
|
pub fn generate_bvh<T>(boxen:Vec<(T,Aabb)>)->BvhNode<T>{
|
|
generate_bvh_node(boxen,false)
|
|
}
|
|
|
|
fn generate_bvh_node<T>(boxen:Vec<(T,Aabb)>,force:bool)->BvhNode<T>{
|
|
let n=boxen.len();
|
|
if force||n<20{
|
|
let mut aabb=Aabb::default();
|
|
let nodes=boxen.into_iter().map(|b|{
|
|
aabb.join(&b.1);
|
|
BvhNode{
|
|
content:RecursiveContent::Leaf(b.0),
|
|
aabb:b.1,
|
|
}
|
|
}).collect();
|
|
BvhNode{
|
|
content:RecursiveContent::Branch(nodes),
|
|
aabb,
|
|
}
|
|
}else{
|
|
let mut sort_x=Vec::with_capacity(n);
|
|
let mut sort_y=Vec::with_capacity(n);
|
|
let mut sort_z=Vec::with_capacity(n);
|
|
for (i,(_,aabb)) in boxen.iter().enumerate(){
|
|
let center=aabb.center();
|
|
sort_x.push((i,center.x));
|
|
sort_y.push((i,center.y));
|
|
sort_z.push((i,center.z));
|
|
}
|
|
sort_x.sort_by_key(|&(_,c)|c);
|
|
sort_y.sort_by_key(|&(_,c)|c);
|
|
sort_z.sort_by_key(|&(_,c)|c);
|
|
let h=n/2;
|
|
let median_x=sort_x[h].1;
|
|
let median_y=sort_y[h].1;
|
|
let median_z=sort_z[h].1;
|
|
//locate a run of values equal to the median
|
|
//partition point gives the first index for which the predicate evaluates to false
|
|
let first_index_eq_median_x=sort_x.partition_point(|&(_,x)|x<median_x);
|
|
let first_index_eq_median_y=sort_y.partition_point(|&(_,y)|y<median_y);
|
|
let first_index_eq_median_z=sort_z.partition_point(|&(_,z)|z<median_z);
|
|
let first_index_gt_median_x=sort_x.partition_point(|&(_,x)|x<=median_x);
|
|
let first_index_gt_median_y=sort_y.partition_point(|&(_,y)|y<=median_y);
|
|
let first_index_gt_median_z=sort_z.partition_point(|&(_,z)|z<=median_z);
|
|
//pick which side median value copies go into such that both sides are as balanced as possible based on distance from n/2
|
|
let partition_point_x=if n.abs_diff(2*first_index_eq_median_x)<n.abs_diff(2*first_index_gt_median_x){first_index_eq_median_x}else{first_index_gt_median_x};
|
|
let partition_point_y=if n.abs_diff(2*first_index_eq_median_y)<n.abs_diff(2*first_index_gt_median_y){first_index_eq_median_y}else{first_index_gt_median_y};
|
|
let partition_point_z=if n.abs_diff(2*first_index_eq_median_z)<n.abs_diff(2*first_index_gt_median_z){first_index_eq_median_z}else{first_index_gt_median_z};
|
|
//this ids which octant the boxen is put in
|
|
let mut octant=vec![0;n];
|
|
for &(i,_) in &sort_x[partition_point_x..]{
|
|
octant[i]+=1<<0;
|
|
}
|
|
for &(i,_) in &sort_y[partition_point_y..]{
|
|
octant[i]+=1<<1;
|
|
}
|
|
for &(i,_) in &sort_z[partition_point_z..]{
|
|
octant[i]+=1<<2;
|
|
}
|
|
//generate lists for unique octant values
|
|
let mut list_list=Vec::with_capacity(8);
|
|
let mut octant_list=Vec::with_capacity(8);
|
|
for (i,(data,aabb)) in boxen.into_iter().enumerate(){
|
|
let octant_id=octant[i];
|
|
let list_id=if let Some(list_id)=octant_list.iter().position(|&id|id==octant_id){
|
|
list_id
|
|
}else{
|
|
let list_id=list_list.len();
|
|
octant_list.push(octant_id);
|
|
list_list.push(Vec::new());
|
|
list_id
|
|
};
|
|
list_list[list_id].push((data,aabb));
|
|
}
|
|
let mut aabb=Aabb::default();
|
|
if list_list.len()==1{
|
|
generate_bvh_node(list_list.remove(0),true)
|
|
}else{
|
|
BvhNode{
|
|
content:RecursiveContent::Branch(
|
|
list_list.into_iter().map(|b|{
|
|
let node=generate_bvh_node(b,false);
|
|
aabb.join(&node.aabb);
|
|
node
|
|
}).collect()
|
|
),
|
|
aabb,
|
|
}
|
|
}
|
|
}
|
|
}
|