//#####################################################################
// Copyright 2002-2007, Robert Bridson, Ronald Fedkiw, Eran Guendelman, Geoffrey Irving, Sergey Koltakov, Neil Molino, Igor Neverov, Duc Nguyen, Avi Robinson-Mosher, Craig Schroeder,
//      Andrew Selle, Jerry Talton, Joseph Teran.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
#include <Core/Math_Tools/RANGE.h>
#include <Core/Log/LOG.h>
#include <Core/Matrices/DIAGONAL_MATRIX.h>
#include <Core/Random_Numbers/RANDOM_NUMBERS.h>
#include <Core/Read_Write/READ_WRITE.h>
#include <Core/Vectors/VECTOR.h>
namespace PhysBAM{
//#####################################################################
// Function Normal
//#####################################################################
template<class TV> TV RANGE<TV>::
Normal(const int aggregate) const
{
    assert(aggregate>=0 && aggregate<d*2);
    TV direction=TV::Axis_Vector(aggregate>>1);
    return (aggregate&1)?direction:-direction;
}
//#####################################################################
// Function Surface
//#####################################################################
template<class TV> TV RANGE<TV>::
Surface(const TV& X) const
{
    if(!Lazy_Inside(X)) return Clamp(X);
    TV min_side=X-min_corner,max_side=max_corner-X,result=X;
    int min_argmin=min_side.Arg_Min(),max_argmin=max_side.Arg_Min();
    if(min_side(min_argmin)<=max_side(max_argmin)) result(min_argmin)=min_corner(min_argmin);
    else result(max_argmin)=max_corner(max_argmin);
    return result;
}
//#####################################################################
// Function Signed_Distance
//#####################################################################
template<class TV> typename TV::SCALAR RANGE<TV>::
Signed_Distance(const TV& X) const
{
    TV lengths=Edge_Lengths();
    if(lengths.Contains(0)) PHYSBAM_NOT_IMPLEMENTED();
    TV phi=abs(X-Center())-(T).5*lengths;
    if(!phi.All_Less_Equal(TV())) return phi.Componentwise_Max(TV()).Magnitude();
    return phi.Max();
}
//#####################################################################
// Function Normal
//#####################################################################
template<class TV> TV RANGE<TV>::
Normal(const TV& X) const
{
    if(Lazy_Inside(X)){
        TV phis_max=X-max_corner,phis_min=min_corner-X;
        int axis=phis_min.Componentwise_Max(phis_max).Arg_Max();
        return T(phis_max[axis]>phis_min[axis]?1:-1)*TV::Axis_Vector(axis);}
    else{
        TV phis_max=X-max_corner,phis_min=min_corner-X;
        TV normal;
        for(int i=0;i<d;i++){
            T phi=max(phis_min(i),phis_max(i));
            normal(i)=phi>0?(phis_max(i)>phis_min(i)?phi:-phi):0;}
        return normal.Normalized();}
}
//#####################################################################
// Function Hessian
//#####################################################################
template<class TV> SYMMETRIC_MATRIX<typename TV::SCALAR,TV::m> RANGE<TV>::
Hessian(const TV& X) const
{
    if(Lazy_Inside(X)) return SYMMETRIC_MATRIX<T,TV::m>();
    TV phis_max=X-max_corner,phis_min=min_corner-X;
    TV normal,D_normal_diag;

    for(int i=0;i<d;i++){
        T phi=max(phis_min(i),phis_max(i));
        if(phi>0){
            normal(i)=phis_max(i)>phis_min(i)?phi:-phi;
            D_normal_diag(i)=1;}}
    T norm=normal.Normalize();
    return (DIAGONAL_MATRIX<T,TV::m>(D_normal_diag)-Outer_Product(normal))/norm;
}
//#####################################################################
// Function Fill_Uniform
//#####################################################################
template<class TV> void RANGE<TV>::
Fill_Uniform(RANDOM_NUMBERS<T>& rand,const RANGE<TV>& box)
{
    rand.Fill_Uniform(min_corner,box.min_corner,box.max_corner);
    rand.Fill_Uniform(max_corner,min_corner,box.max_corner);
}
//#####################################################################
// Function Name
//#####################################################################
template<class TV> std::string RANGE<TV>::
Name()
{
    return LOG::sprintf("RANGE<VECTOR<T,%d>",d);
}
//#####################################################################
// Function Read
//#####################################################################
template<class TV> template<class RW> void RANGE<TV>::
Read(UNTYPED_ISTREAM<RW>& input)
{
    if(is_same<RW,char>::value)
        Read_Binary_Group(input,min_corner,max_corner);
    else
        for(int i=0;i<TV::m;i++)
            Read_Binary(input,min_corner(i),max_corner(i));
}
//#####################################################################
// Function Write
//#####################################################################
template<class TV> template<class RW> void RANGE<TV>::
Write(UNTYPED_OSTREAM<RW>& output) const
{
    if(is_same<RW,char>::value)
        Write_Binary_Group(output,min_corner,max_corner);
    else
        for(int i=0;i<TV::m;i++)
            Write_Binary(output,min_corner(i),max_corner(i));
}
//#####################################################################
#define INSTANTIATION_HELPER(T,d) \
    template std::string RANGE<VECTOR<T,d> >::Name(); \
    template VECTOR<T,d> RANGE<VECTOR<T,d> >::Normal(const int) const; \
    template VECTOR<T,d> RANGE<VECTOR<T,d> >::Normal(const VECTOR<T,d>&) const; \
    template VECTOR<T,d> RANGE<VECTOR<T,d> >::Surface(const VECTOR<T,d>&) const; \
    template SYMMETRIC_MATRIX<T,d> RANGE<VECTOR<T,d> >::Hessian(const VECTOR<T,d>&) const; \
    template VECTOR<T,d>::SCALAR RANGE<VECTOR<T,d> >::Signed_Distance(const VECTOR<T,d>&) const;

INSTANTIATION_HELPER(float,1);
INSTANTIATION_HELPER(float,2);
INSTANTIATION_HELPER(float,3);
INSTANTIATION_HELPER(double,1);
INSTANTIATION_HELPER(double,2);
INSTANTIATION_HELPER(double,3);
template void RANGE<VECTOR<double,1> >::Read(UNTYPED_ISTREAM<double>&);
template void RANGE<VECTOR<double,1> >::Read(UNTYPED_ISTREAM<float>&);
template void RANGE<VECTOR<double,1> >::Read(UNTYPED_ISTREAM<char>&);
template void RANGE<VECTOR<double,1> >::Write(UNTYPED_OSTREAM<double>&) const;
template void RANGE<VECTOR<double,1> >::Write(UNTYPED_OSTREAM<float>&) const;
template void RANGE<VECTOR<double,1> >::Write(UNTYPED_OSTREAM<char>&) const;
template void RANGE<VECTOR<double,2> >::Read(UNTYPED_ISTREAM<double>&);
template void RANGE<VECTOR<double,2> >::Read(UNTYPED_ISTREAM<float>&);
template void RANGE<VECTOR<double,2> >::Read(UNTYPED_ISTREAM<char>&);
template void RANGE<VECTOR<double,2> >::Write(UNTYPED_OSTREAM<double>&) const;
template void RANGE<VECTOR<double,2> >::Write(UNTYPED_OSTREAM<float>&) const;
template void RANGE<VECTOR<double,2> >::Write(UNTYPED_OSTREAM<char>&) const;
template void RANGE<VECTOR<double,3> >::Read(UNTYPED_ISTREAM<double>&);
template void RANGE<VECTOR<double,3> >::Read(UNTYPED_ISTREAM<float>&);
template void RANGE<VECTOR<double,3> >::Read(UNTYPED_ISTREAM<char>&);
template void RANGE<VECTOR<double,3> >::Write(UNTYPED_OSTREAM<double>&) const;
template void RANGE<VECTOR<double,3> >::Write(UNTYPED_OSTREAM<float>&) const;
template void RANGE<VECTOR<double,3> >::Write(UNTYPED_OSTREAM<char>&) const;
template void RANGE<VECTOR<float,1> >::Read(UNTYPED_ISTREAM<double>&);
template void RANGE<VECTOR<float,1> >::Read(UNTYPED_ISTREAM<float>&);
template void RANGE<VECTOR<float,1> >::Read(UNTYPED_ISTREAM<char>&);
template void RANGE<VECTOR<float,1> >::Write(UNTYPED_OSTREAM<double>&) const;
template void RANGE<VECTOR<float,1> >::Write(UNTYPED_OSTREAM<float>&) const;
template void RANGE<VECTOR<float,1> >::Write(UNTYPED_OSTREAM<char>&) const;
template void RANGE<VECTOR<float,2> >::Read(UNTYPED_ISTREAM<double>&);
template void RANGE<VECTOR<float,2> >::Read(UNTYPED_ISTREAM<float>&);
template void RANGE<VECTOR<float,2> >::Read(UNTYPED_ISTREAM<char>&);
template void RANGE<VECTOR<float,2> >::Write(UNTYPED_OSTREAM<double>&) const;
template void RANGE<VECTOR<float,2> >::Write(UNTYPED_OSTREAM<float>&) const;
template void RANGE<VECTOR<float,2> >::Write(UNTYPED_OSTREAM<char>&) const;
template void RANGE<VECTOR<float,3> >::Read(UNTYPED_ISTREAM<double>&);
template void RANGE<VECTOR<float,3> >::Read(UNTYPED_ISTREAM<float>&);
template void RANGE<VECTOR<float,3> >::Read(UNTYPED_ISTREAM<char>&);
template void RANGE<VECTOR<float,3> >::Write(UNTYPED_OSTREAM<double>&) const;
template void RANGE<VECTOR<float,3> >::Write(UNTYPED_OSTREAM<float>&) const;
template void RANGE<VECTOR<float,3> >::Write(UNTYPED_OSTREAM<char>&) const;
template void RANGE<VECTOR<int,1> >::Read(UNTYPED_ISTREAM<double>&);
template void RANGE<VECTOR<int,1> >::Read(UNTYPED_ISTREAM<float>&);
template void RANGE<VECTOR<int,1> >::Read(UNTYPED_ISTREAM<char>&);
template void RANGE<VECTOR<int,1> >::Write(UNTYPED_OSTREAM<double>&) const;
template void RANGE<VECTOR<int,1> >::Write(UNTYPED_OSTREAM<float>&) const;
template void RANGE<VECTOR<int,1> >::Write(UNTYPED_OSTREAM<char>&) const;
template void RANGE<VECTOR<int,2> >::Read(UNTYPED_ISTREAM<double>&);
template void RANGE<VECTOR<int,2> >::Read(UNTYPED_ISTREAM<float>&);
template void RANGE<VECTOR<int,2> >::Read(UNTYPED_ISTREAM<char>&);
template void RANGE<VECTOR<int,2> >::Write(UNTYPED_OSTREAM<double>&) const;
template void RANGE<VECTOR<int,2> >::Write(UNTYPED_OSTREAM<float>&) const;
template void RANGE<VECTOR<int,2> >::Write(UNTYPED_OSTREAM<char>&) const;
template void RANGE<VECTOR<int,3> >::Read(UNTYPED_ISTREAM<double>&);
template void RANGE<VECTOR<int,3> >::Read(UNTYPED_ISTREAM<float>&);
template void RANGE<VECTOR<int,3> >::Read(UNTYPED_ISTREAM<char>&);
template void RANGE<VECTOR<int,3> >::Write(UNTYPED_OSTREAM<float>&) const;
template void RANGE<VECTOR<int,3> >::Write(UNTYPED_OSTREAM<char>&) const;
template void RANGE<VECTOR<int,3> >::Write(UNTYPED_OSTREAM<double>&) const;
template void RANGE<VECTOR<double,1> >::Fill_Uniform(RANDOM_NUMBERS<double>&,RANGE<VECTOR<double,1> > const&);
template void RANGE<VECTOR<double,2> >::Fill_Uniform(RANDOM_NUMBERS<double>&,RANGE<VECTOR<double,2> > const&);
template void RANGE<VECTOR<double,3> >::Fill_Uniform(RANDOM_NUMBERS<double>&,RANGE<VECTOR<double,3> > const&);
template void RANGE<VECTOR<float,1> >::Fill_Uniform(RANDOM_NUMBERS<float>&,RANGE<VECTOR<float,1> > const&);
template void RANGE<VECTOR<float,2> >::Fill_Uniform(RANDOM_NUMBERS<float>&,RANGE<VECTOR<float,2> > const&);
template void RANGE<VECTOR<float,3> >::Fill_Uniform(RANDOM_NUMBERS<float>&,RANGE<VECTOR<float,3> > const&);
}
