//#####################################################################
// Copyright 2002-2008, Robert Bridson, Ronald Fedkiw, Eran Guendelman, Jon Gretarsson, Andrew Selle.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
#include <Core/Basic_Geometry/RAY.h>
#include <Core/Basic_Geometry/RAY_ACCEL.h>
#include <Core/Basic_Geometry/RAY_HIT.h>
namespace PhysBAM{
//#####################################################################
// Function Intersects
//#####################################################################
template<class TV> RAY_HIT<typename TV::SCALAR>
Intersects(const RAY<TV>& ray,const RANGE<TV>& box,typename TV::SCALAR isocontour)
{
    typedef typename TV::SCALAR T;
    T start_t,end_t;
    if(!Get_Intersection_Range(ray,box.Thickened(isocontour),start_t,end_t))
        return {};
    if(start_t>0) return RAY_HIT<T>{start_t};
    return RAY_HIT<T>{end_t};
}
//#####################################################################
// Function Intersection_Interval
//#####################################################################
template<class TV,class T> VECTOR<RAY_HIT<T>,2>
Intersection_Interval(const RAY<TV>& ray,const RANGE<TV>& box,T isocontour)
{
    typedef RAY_HIT<T> H;
    T start_t,end_t;
    if(!Get_Intersection_Range(ray,box.Thickened(isocontour),start_t,end_t))
        return {};
    return {H{start_t},H{end_t}};
}
//#####################################################################
// Function All_Intersections
//#####################################################################
template<class TV,class T> void
All_Intersections(ARRAY<RAY_HIT<T> >& hits,const RAY<TV>& ray,
    const RANGE<TV>& box,T isocontour)
{
    typedef RAY_HIT<T> H;
    T start_t,end_t;
    if(Get_Intersection_Range(ray,box.Thickened(isocontour),start_t,end_t))
    {
        hits.Append(H{start_t});
        hits.Append(H{end_t});
    }
}
//#####################################################################
// Function Get_Intersection_Range
//#####################################################################
template<class TV> bool Get_Intersection_Range(const RAY<TV>& ray,const RANGE<TV>& box,typename TV::SCALAR& start_t,typename TV::SCALAR& end_t)
{
    typedef typename TV::SCALAR T;
    // This comes from a paper "An efficient and Robust Ray-Box Intersection Algorithm" by Williams, Barrus, Morley, and Shirley
    // http://www.cs.utah.edu/~rmorley/pubs/box.pdf
    RAY_ACCEL<TV> ra(ray);
    TV extremes[2];
    extremes[0]=box.min_corner; // low extreme
    extremes[1]=box.max_corner; // high extreme
    T t_min=(extremes[ra.direction_is_negative.x].x-ray.endpoint.x)*ra.inverse_direction.x;
    T t_max=(extremes[1-ra.direction_is_negative.x].x-ray.endpoint.x)*ra.inverse_direction.x;
    if constexpr(TV::m>=2)
    {
        T t_y_min=(extremes[ra.direction_is_negative.y].y-ray.endpoint.y)*ra.inverse_direction.y;
        T t_y_max=(extremes[1-ra.direction_is_negative.y].y-ray.endpoint.y)*ra.inverse_direction.y;
        if(t_min>t_y_max || t_y_min>t_max) return false;
        if(t_y_min>t_min)t_min=t_y_min;
        if(t_y_max<t_max)t_max=t_y_max;
        if constexpr(TV::m>=3)
        {
            T t_z_min=(extremes[ra.direction_is_negative.z].z-ray.endpoint.z)*ra.inverse_direction.z;
            T t_z_max=(extremes[1-ra.direction_is_negative.z].z-ray.endpoint.z)*ra.inverse_direction.z;
            if(t_min>t_z_max || t_z_min>t_max) return false;
            if(t_z_min>t_min)t_min=t_z_min;
            if(t_z_max<t_max)t_max=t_z_max;
        }
    }
    if(t_max<=0) return false;
    start_t=max((T)0,t_min);
    end_t=t_max;
    return true;
}
template RAY_HIT<float> Intersects(const RAY<VECTOR<float,1> >&,const RANGE<VECTOR<float,1> >&,float);
template RAY_HIT<float> Intersects(const RAY<VECTOR<float,2> >&,const RANGE<VECTOR<float,2> >&,float);
template RAY_HIT<float> Intersects(const RAY<VECTOR<float,3> >&,const RANGE<VECTOR<float,3> >&,float);
template bool Get_Intersection_Range(const RAY<VECTOR<float,1> >&,const RANGE<VECTOR<float,1> >&,float&,float&);
template bool Get_Intersection_Range(const RAY<VECTOR<float,2> >&,const RANGE<VECTOR<float,2> >&,float&,float&);
template bool Get_Intersection_Range(const RAY<VECTOR<float,3> >&,const RANGE<VECTOR<float,3> >&,float&,float&);
template RAY_HIT<double> Intersects(const RAY<VECTOR<double,1> >&,const RANGE<VECTOR<double,1> >&,double);
template RAY_HIT<double> Intersects(const RAY<VECTOR<double,2> >&,const RANGE<VECTOR<double,2> >&,double);
template RAY_HIT<double> Intersects(const RAY<VECTOR<double,3> >&,const RANGE<VECTOR<double,3> >&,double);
template bool Get_Intersection_Range(const RAY<VECTOR<double,1> >&,const RANGE<VECTOR<double,1> >&,double&,double&);
template bool Get_Intersection_Range(const RAY<VECTOR<double,2> >&,const RANGE<VECTOR<double,2> >&,double&,double&);
template bool Get_Intersection_Range(const RAY<VECTOR<double,3> >&,const RANGE<VECTOR<double,3> >&,double&,double&);
template void All_Intersections(ARRAY<RAY_HIT<double> >&,const RAY<VECTOR<double,1> >&,const RANGE<VECTOR<double,1> >&,double);
template void All_Intersections(ARRAY<RAY_HIT<double> >&,const RAY<VECTOR<double,2> >&,const RANGE<VECTOR<double,2> >&,double);
template void All_Intersections(ARRAY<RAY_HIT<double> >&,const RAY<VECTOR<double,3> >&,const RANGE<VECTOR<double,3> >&,double);
template void All_Intersections(ARRAY<RAY_HIT<float> >&,const RAY<VECTOR<float,1> >&,const RANGE<VECTOR<float,1> >&,float);
template void All_Intersections(ARRAY<RAY_HIT<float> >&,const RAY<VECTOR<float,2> >&,const RANGE<VECTOR<float,2> >&,float);
template void All_Intersections(ARRAY<RAY_HIT<float> >&,const RAY<VECTOR<float,3> >&,const RANGE<VECTOR<float,3> >&,float);
template VECTOR<RAY_HIT<double>,2> Intersection_Interval<VECTOR<double,2>,double>(RAY<VECTOR<double,2> > const&,RANGE<VECTOR<double,2> > const&,double);
template VECTOR<RAY_HIT<double>,2> Intersection_Interval<VECTOR<double,3>,double>(RAY<VECTOR<double,3> > const&,RANGE<VECTOR<double,3> > const&,double);
template VECTOR<RAY_HIT<float>,2> Intersection_Interval<VECTOR<float,2>,float>(RAY<VECTOR<float,2> > const&,RANGE<VECTOR<float,2> > const&,float);
template VECTOR<RAY_HIT<float>,2> Intersection_Interval<VECTOR<float,3>,float>(RAY<VECTOR<float,3> > const&,RANGE<VECTOR<float,3> > const&,float);
}
