//#####################################################################
// Copyright 2002, 2003, Robert Bridson, Ronald Fedkiw, Eran Guendelman.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class PLANE
//##################################################################### 
#include <Core/Basic_Geometry/PLANE.h>
#include <Core/Basic_Geometry/RAY.h>
#include <Core/Basic_Geometry/RAY_HIT.h>
#include <Core/Implicit_Objects/INTERSECTION_HELPERS.h>
#include <Core/Math_Tools/RANGE.h>
#include <Core/Matrices/FRAME.h>
#include <Core/Meshes/SEGMENTED_CURVE.h>
#include <Core/Meshes/TRIANGULATED_SURFACE.h>
#include <Core/Random_Numbers/RANDOM_NUMBERS.h>
#include <Core/Read_Write/READ_WRITE.h>
namespace PhysBAM{
//#####################################################################
// Segment_Plane_Intersection
//#####################################################################
template<class T> bool PLANE<T>::
Segment_Plane_Intersection(const TV& endpoint1,const TV& endpoint2,T& interpolation_fraction) const
{
    T denominator=normal.Dot(endpoint2-endpoint1);
    if(!denominator){interpolation_fraction=FLT_MAX;return false;} // parallel
    interpolation_fraction=normal.Dot(x0-endpoint1)/denominator;
    return (interpolation_fraction>=0 && interpolation_fraction<=1);
}
//#####################################################################
// Function Fill_Uniform
//#####################################################################
template<class T> void PLANE<T>::
Fill_Uniform(RANDOM_NUMBERS<T>& rand,const RANGE<TV>& box)
{
    rand.Fill_Uniform(x0,box.min_corner,box.max_corner);
    normal=rand.template Get_Direction<TV>();
}
//#####################################################################
// Function Read
//#####################################################################
template<class T> template<class RW> void PLANE<T>::
Read(UNTYPED_ISTREAM<RW>& input)
{
    Read_Binary_Group(input,normal,x0);
}
//#####################################################################
// Function Write
//#####################################################################
template<class T> template<class RW> void PLANE<T>::
Write(UNTYPED_OSTREAM<RW>& output) const
{
    Write_Binary_Group(output,normal,x0);
}
//#####################################################################
// Function Scaled
//#####################################################################
template<class T,class TV> PLANE<T>
Scaled(const PLANE<T>& b,const TV& center,T scale)
{
    return PLANE<T>(b.normal,(b.x0-center)*scale+center);
}
//#####################################################################
// Function Transformed
//#####################################################################
template<class T,class TV> PLANE<T>
Transformed(const PLANE<T>& b,const FRAME<TV>& f)
{
    return PLANE<T>(f.r.Rotate(b.normal),f*b.x0);
}
//#####################################################################
// Function Generate_Triangles
//#####################################################################
template<class T> TRIANGULATED_SURFACE<T>*
Generate_Triangles(const PLANE<T>& plane,T tesselated_size,int elements_per_side)
{
    typedef VECTOR<T,3> TV;
    TRIANGULATED_SURFACE<T>* surface=TRIANGULATED_SURFACE<T>::Create();
    int m=elements_per_side+1,n=elements_per_side;
    surface->particles.Add_Elements(m*m);
    surface->mesh.Initialize_Herring_Bone_Mesh(m,m);
    TV u=tesselated_size*plane.normal.Unit_Orthogonal_Vector(),v=plane.normal.Cross(u);
    for(int j=0,p=0;j<=n;j++)
        for(int i=0;i<=n;i++)
            surface->particles.X(p++)=plane.x0+u*((T)2*j/n-1)-v*((T)2*i/n-1);
    surface->Update_Number_Nodes();
    return surface;
}
//#####################################################################
// Function Intersects
//#####################################################################
template<class T> RAY_HIT<T> Intersects(const RAY<VECTOR<T,3> >& ray,const PLANE<T>& plane,T isocontour)
{
    T distance=plane.normal.Dot(ray.endpoint-plane.x0)-isocontour;
    T rate_of_approach=-plane.normal.Dot(ray.direction);
    if(rate_of_approach*distance<=0) return {}; // no intersection
    return RAY_HIT<T>{distance/rate_of_approach};
}
//#####################################################################
// Function Intersection_Interval
//#####################################################################
template<class T> VECTOR<RAY_HIT<T>,2>
Intersection_Interval(const RAY<VECTOR<T,3> >& ray,const PLANE<T>& plane,T isocontour)
{
    T distance=plane.normal.Dot(ray.endpoint-plane.x0)-isocontour;
    T rate_of_approach=-plane.normal.Dot(ray.direction);
    VECTOR<RAY_HIT<T>,2> v;
    if(distance<0)
    {
        v(0).Set(0);
        if(rate_of_approach<0) v(1).Set(distance/rate_of_approach);
    }
    else if(rate_of_approach>0) v(0).Set(distance/rate_of_approach);
    return v;
}
//#####################################################################
// Function Intersects
//#####################################################################
template<class T> RAY_HIT<T> Intersects(const RAY<VECTOR<T,3> >& ray,
    const PLANE<T>& plane, const T thickness_over_two,const T distance,
    const T rate_of_approach)
{
    RAY_HIT<T> hit;
    if(distance>-thickness_over_two && distance<thickness_over_two)
        return RAY_HIT<T>{0}; // within the boundary
    if(rate_of_approach*distance<=0) return {}; // no intersection
    if(rate_of_approach>0) hit.Set((distance-thickness_over_two)/rate_of_approach);
    else hit.Set((distance+thickness_over_two)/rate_of_approach);
    if(hit.dist>FLT_MAX || hit.dist<-FLT_MAX) return {}; // guard against overflow
    return hit;
}
//#####################################################################
// Function Intersects
//#####################################################################
template<class T> RAY_HIT<T> Rectangle_Intersects(const RAY<VECTOR<T,3> >& ray,
    const PLANE<T>& plane,const PLANE<T>& bounding_plane1,
    const PLANE<T>& bounding_plane2,const PLANE<T>& bounding_plane3,
    const PLANE<T>& bounding_plane4,const T thickness_over_two)
{
    auto hit=Intersects(ray,plane,thickness_over_two);
    if(!hit.Hit()) return hit;
    VECTOR<T,3> point=ray.Point(hit.dist);
    if(bounding_plane1.Outside(point,thickness_over_two)) return {};
    if(bounding_plane2.Outside(point,thickness_over_two)) return {};
    if(bounding_plane3.Outside(point,thickness_over_two)) return {};
    if(bounding_plane4.Outside(point,thickness_over_two)) return {};
    return hit;
}
//#####################################################################
// Function All_Intersections
//#####################################################################
template<class T> void
All_Intersections(ARRAY<RAY_HIT<T> >& hits,const RAY<VECTOR<T,3> >& ray,const PLANE<T>& plane, T isocontour)
{
    auto iv=Intersection_Interval(ray,plane,isocontour);
    hits=iv.Array_View(0,Actual_Length(iv));
}
//#####################################################################
template class PLANE<float>;
template class PLANE<double>;
template PLANE<double> Scaled<double,VECTOR<double,3> >(PLANE<double> const&,VECTOR<double,3> const&,double);
template PLANE<double> Transformed<double,VECTOR<double,3> >(PLANE<double> const&,FRAME<VECTOR<double,3> > const&);
template PLANE<float> Scaled<float,VECTOR<float,3> >(PLANE<float> const&,VECTOR<float,3> const&,float);
template PLANE<float> Transformed<float,VECTOR<float,3> >(PLANE<float> const&,FRAME<VECTOR<float,3> > const&);
template TRIANGULATED_SURFACE<double>* Generate_Triangles<double>(PLANE<double> const&,double,int);
template TRIANGULATED_SURFACE<float>* Generate_Triangles<float>(PLANE<float> const&,float,int);
template RAY_HIT<double> Intersects<double>(RAY<VECTOR<double,3> > const&,PLANE<double> const&,double);
template RAY_HIT<float> Intersects<float>(RAY<VECTOR<float,3> > const&,PLANE<float> const&,float);
template VECTOR<RAY_HIT<double>,2> Intersection_Interval<double>(RAY<VECTOR<double,3> > const&,PLANE<double> const&,double);
template VECTOR<RAY_HIT<float>,2> Intersection_Interval<float>(RAY<VECTOR<float,3> > const&,PLANE<float> const&,float);
template RAY_HIT<float> Intersects(const RAY<VECTOR<float,3> >&,const PLANE<float>&,const float,const float,const float);
template RAY_HIT<float> Rectangle_Intersects(const RAY<VECTOR<float,3> >&,const PLANE<float>&,const PLANE<float>&,const PLANE<float>&,const PLANE<float>&,const PLANE<float>&,const float);
template RAY_HIT<double> Intersects(const RAY<VECTOR<double,3> >&,const PLANE<double>&,const double,const double,const double);
template RAY_HIT<double> Rectangle_Intersects(const RAY<VECTOR<double,3> >&,const PLANE<double>&,const PLANE<double>&,const PLANE<double>&,const PLANE<double>&,const PLANE<double>&,const double);
template void All_Intersections<double>(ARRAY<RAY_HIT<double>,int>&,RAY<VECTOR<double,3> > const&,PLANE<double> const&,double);
template void All_Intersections<float>(ARRAY<RAY_HIT<float>,int>&,RAY<VECTOR<float,3> > const&,PLANE<float> const&,float);
template void PLANE<double>::Read(UNTYPED_ISTREAM<char>&);
template void PLANE<double>::Read(UNTYPED_ISTREAM<double>&);
template void PLANE<double>::Read(UNTYPED_ISTREAM<float>&);
template void PLANE<double>::Write(UNTYPED_OSTREAM<char>&) const;
template void PLANE<double>::Write(UNTYPED_OSTREAM<double>&) const;
template void PLANE<double>::Write(UNTYPED_OSTREAM<float>&) const;
template void PLANE<float>::Read(UNTYPED_ISTREAM<char>&);
template void PLANE<float>::Read(UNTYPED_ISTREAM<double>&);
template void PLANE<float>::Read(UNTYPED_ISTREAM<float>&);
template void PLANE<float>::Write(UNTYPED_OSTREAM<char>&) const;
template void PLANE<float>::Write(UNTYPED_OSTREAM<double>&) const;
template void PLANE<float>::Write(UNTYPED_OSTREAM<float>&) const;
}
