//#####################################################################
// Copyright 2002-2007, Ronald Fedkiw, Eran Guendelman, Geoffrey Irving, Nipun Kwatra, Avi Robinson-Mosher, Craig Schroeder, Tamar Shinar, Jerry Talton.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class ORIENTED_BOX
//#####################################################################
#ifndef __ORIENTED_BOX__
#define __ORIENTED_BOX__

#include <Core/Basic_Geometry/BASIC_GEOMETRY_FORWARD.h>
#include <Core/Log/DEBUG_UTILITIES.h>
#include <Core/Math_Tools/RANGE.h>
#include <Core/Matrices/DIAGONAL_MATRIX.h>
#include <Core/Matrices/MATRIX.h>
#include <Core/Matrices/MATRIX_POLICY.h>
#include <Core/Matrices/SYMMETRIC_MATRIX.h>
#include <Core/Meshes/MESHES_FORWARD.h>
namespace PhysBAM{

template<class TV>
class ORIENTED_BOX
{
    typedef typename TV::SCALAR T;
    enum WORKAROUND{d=TV::m};
public:
    typedef TV VECTOR_T;

    TV corner; // root corner of the box
    MATRIX<T,d> edges; // principal edges of the box, emanating from the corner

    ORIENTED_BOX()
        :corner(TV()),edges(MATRIX<T,d>::Identity_Matrix())
    {}

    ORIENTED_BOX(const TV& corner_input,const MATRIX<T,d>& edges_input)
        :corner(corner_input),edges(edges_input)
    {}

    ORIENTED_BOX(const RANGE<TV>& box,const ROTATION<TV>& rotation)
    {
        corner=rotation.Rotate(box.min_corner);
        edges=rotation.Rotation_Matrix()*DIAGONAL_MATRIX<T,TV::m>(box.Edge_Lengths());
    }

    ORIENTED_BOX(const RANGE<TV>& box,const ROTATION<TV>& rotation,const TV& corner_input)
        :corner(corner_input)
    {
        edges=rotation.Rotation_Matrix()*DIAGONAL_MATRIX<T,TV::m>(box.Edge_Lengths());
    }

    ORIENTED_BOX(const RANGE<TV>& box,const FRAME<TV>& frame)
    {
        corner=frame*box.min_corner;
        edges=frame.r.Rotation_Matrix()*DIAGONAL_MATRIX<T,TV::m>(box.Edge_Lengths());
    }

    TV Center() const
    {
        TV u;
        for(int i=0;i<TV::m;i++) u+=edges.Column(i);
        return corner+(T).5*u;
    }

    void Scale_About_Center(const T factor)
    {
        TV center=Center();
        edges*=factor;
        corner=(corner-center)*factor+center;
    }

    void Scale_About_Center(const TV factor)
    {
        TV center=Center();
        edges=edges*DIAGONAL_MATRIX<T,TV::m>(factor);
        corner=(corner-center)*factor+center;
    }

    ORIENTED_BOX Scaled_About_Center(const T factor) const
    {ORIENTED_BOX box(*this);box.Scale_About_Center(factor);return box;}

    ORIENTED_BOX Scaled_About_Center(const TV edge_factor)
    {ORIENTED_BOX box(*this);box.Scale_About_Center(edge_factor);return box;}

    RANGE<TV> Axis_Aligned_Bounding_Box() const
    {RANGE<TV> box(corner);
    for(int i=0;i<d;i++) box.Enlarge_By_Sign(edges.Column(i));
    return box;}

    RANGE<TV> Bounding_Box() const // for templatization purposes
    {return Axis_Aligned_Bounding_Box();}

//#####################################################################
    bool Lazy_Inside(const TV &location) const;
    bool Intersection(const ORIENTED_BOX& box) const;
    bool Intersection(const RANGE<TV>& box) const;
    T Signed_Distance(const TV& X) const;
    TV Normal(const TV& X) const;
    SYMMETRIC_MATRIX<T,TV::m> Hessian(const TV& X) const;
    bool Separating_Test(const ORIENTED_BOX& box,const TV& plane_normal_direction) const;
    bool Separating_Test(const RANGE<TV>& box,const TV& plane_normal_direction) const;
    static std::string Name();
    // Handles a zero-length edge
    ORIENTED_BOX Thickened(const T thickness) 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);
private:
    void Project_Points_Onto_Line(const TV& direction,T& line_min,T& line_max) const;
    
public:
    T Volume() const; 
//#####################################################################
};
template<class TV>
inline std::ostream& operator<<(std::ostream& output_stream,const ORIENTED_BOX<TV>& box)
{
    output_stream<<"("<<box.corner<<")   (";
    for(int i=0;i<TV::m-1;i++) output_stream<<box.edges.Column(i)<<" : ";
    output_stream<<box.edges.Column(TV::m-1)<<")";
    return output_stream;
}

template<class T> class RAY_HIT;
template<class TV> class RAY;

template<class T,class TV> ORIENTED_BOX<TV>
Scaled(const ORIENTED_BOX<TV>& b,const TV& center,T scale);
template<class TV> ORIENTED_BOX<TV>
Transformed(const ORIENTED_BOX<TV>& b,const FRAME<TV>& f);
template<class TV> ORIENTED_BOX<TV>
Transformed(const RANGE<TV>& b,const FRAME<TV>& f)
{return ORIENTED_BOX<TV>(b,f);}

template<class T> TRIANGULATED_SURFACE<T>*
Generate_Triangles(const ORIENTED_BOX<VECTOR<T,3> >& box,VECTOR<int,3> res=(VECTOR<int,3>()+1));
// Note: isocontour is only approximate
template<class TV> RAY_HIT<typename TV::SCALAR>
Intersects(const RAY<TV>& ray,const ORIENTED_BOX<TV>& box, typename TV::SCALAR isocontour=0);
template<class T> RAY_HIT<T>
Fuzzy_Intersects(const RAY<VECTOR<T,2> >& ray,
    const ORIENTED_BOX<VECTOR<T,2> >& box,const T segment_intersect_epsilon=0);
template<class TV,class T> void
All_Intersections(ARRAY<RAY_HIT<T> >& hits,const RAY<TV>& ray,const ORIENTED_BOX<TV>& box, T isocontour=0);
template<class TV> VECTOR<RAY_HIT<typename TV::SCALAR>,2>
Intersection_Interval(const RAY<TV>& ray,const ORIENTED_BOX<TV>& box,
    typename TV::SCALAR isocontour=0);
}
#endif
