//#####################################################################
// Copyright 2002-2006, Robert Bridson, Ronald Fedkiw, Eran Guendelman, Geoffrey Irving, Eftychios Sifakis.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class PLANE
//#####################################################################
#ifndef __PLANE__
#define __PLANE__

#include <Core/Basic_Geometry/BASIC_GEOMETRY_FORWARD.h>
#include <Core/Matrices/SYMMETRIC_MATRIX.h>
#include <Core/Meshes/MESHES_FORWARD.h>
#include <Core/Vectors/VECTOR.h>
namespace PhysBAM{

template<class T>
class PLANE
{
    typedef VECTOR<T,3> TV;
public:
    typedef TV VECTOR_T;

    TV normal;
    TV x0; // point on the plane

    PLANE()
        :normal(0,1,0),x0(0,0,0)
    {}

    PLANE(const TV& normal_input,const TV& x1_input)
        :normal(normal_input),x0(x1_input)
    {}

    PLANE(const TV& x1_input,const TV& x2_input,const TV& x3_input)
    {
        Specify_Three_Points(x1_input,x2_input,x3_input);
    }

    void Specify_Three_Points(const TV& x1_input,const TV& x2_input,const TV& x3_input)
    {normal=Normal(x1_input,x2_input,x3_input);x0=x1_input;}

    static TV Normal(const TV& x0,const TV& x1,const TV& x2)
    {return (x1-x0).Cross(x2-x0).Normalized();}

    template<class T_ARRAY>
    static TV Normal(const T_ARRAY& X)
    {
        static_assert(T_ARRAY::fixed_m<0 || T_ARRAY::fixed_m==3);
        if constexpr(T_ARRAY::fixed_m<0){assert(X.Size()==3);}
        return Normal(X(0),X(1),X(2));
    }

    TV Normal() const
    {return normal;}

    TV Normal(const TV& X) const
    {return normal;}

    static TV Normal_Direction(const TV& x0,const TV& x1,const TV& x2)
    {return (x1-x0).Cross(x2-x0);} // can have any magnitude

    template<class T_ARRAY>
    static TV Normal_Direction(const T_ARRAY& X)
    {
        static_assert(T_ARRAY::fixed_m<0 || T_ARRAY::fixed_m==3);
        if constexpr(T_ARRAY::fixed_m<0){assert(X.Size()==3);}
        return Normal_Direction(X(0),X(1),X(2));
    }

    T Signed_Distance(const TV& location) const
    {return normal.Dot(location-x0);}

    // inside is the half space behind the normal
    bool Inside(const TV& location,const T thickness_over_two) const
    {return Signed_Distance(location)<=-thickness_over_two;}

    bool Lazy_Inside(const TV& location) const
    {return Signed_Distance(location)<=0;}

    bool Outside(const TV& location,const T thickness_over_two) const
    {return !Inside(location,-thickness_over_two);}

    bool Lazy_Outside(const TV& location) const
    {return !Lazy_Inside(location);}

    bool Boundary(const TV& location,const T thickness_over_two) const
    {return abs(Signed_Distance(location))<thickness_over_two;}

    // closest point on the surface
    TV Surface(const TV& location) const
    {return location-Signed_Distance(location)*normal;}

    bool Segment_Intersection(const TV& endpoint1,const TV& endpoint2,T& interpolation_fraction) const
    {return Segment_Plane_Intersection(endpoint1,endpoint2,interpolation_fraction);}

    RANGE<TV> Bounding_Box() const
    {return RANGE<TV>::Full_Box();}

    SYMMETRIC_MATRIX<T,TV::m> Hessian(const TV& X) const
    {return SYMMETRIC_MATRIX<T,TV::m>();}

    static std::string Name()
    {return "PLANE<T>";}

//#####################################################################
    bool Segment_Plane_Intersection(const TV& endpoint1,const TV& endpoint2,T& interpolation_fraction) const;
    template<class RW> void Read(UNTYPED_ISTREAM<RW>& input);
    template<class RW> void Write(UNTYPED_OSTREAM<RW>& output) const;
    void Fill_Uniform(RANDOM_NUMBERS<T>& rand,const RANGE<TV>& box);
//#####################################################################
};

template<class T> class RAY_HIT;
template<class TV> class RAY;
template<class T,class TV> PLANE<T> Scaled(const PLANE<T>& b,const TV& center,T scale);
template<class T,class TV> PLANE<T> Transformed(const PLANE<T>& b,const FRAME<TV>& f);

template<class T> TRIANGULATED_SURFACE<T>* Generate_Triangles(const PLANE<T>& plane,
    T tesselated_size=1000,int elements_per_side=1);
template<class T> TRIANGULATED_SURFACE<T>* Tessellate_Boundary(const PLANE<T>& plane,
    T tesselated_size=1000,int elements_per_side=1)
{return Generate_Triangles(plane,tesselated_size,elements_per_side);}

template<class T> RAY_HIT<T> Intersects(const RAY<VECTOR<T,3> >& ray,
    const PLANE<T>& plane,T thickness=0);
template<class T> VECTOR<RAY_HIT<T>,2> Intersection_Interval(
    const RAY<VECTOR<T,3> >& ray,const PLANE<T>& plane,T thickness=0);
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);
template<class T> RAY_HIT<T> Rectangle_Intersects(const RAY<VECTOR<T,3> >& ray,
    const PLANE<T>& plane1,const PLANE<T>& plane2,const PLANE<T>& plane3,
    const PLANE<T>& plane4,const PLANE<T>& plane5,const T thickness_over_two=0);
template<class T> void
All_Intersections(ARRAY<RAY_HIT<T> >& hits,const RAY<VECTOR<T,3> >& ray,const PLANE<T>& plane, T isocontour=0);
}
#endif
