//#####################################################################
// Copyright 2003-2008, Eran Guendelman, Geoffrey Irving, Sergey Koltakov.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class RING
//##################################################################### 
#include <Core/Basic_Geometry/RING.h>
#include <Core/Auto_Diff/AUTO_HESS.h>
#include <Core/Basic_Geometry/CYLINDER.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/TRIANGULATED_SURFACE.h>
#include <Core/Random_Numbers/RANDOM_NUMBERS.h>
#include <Core/Read_Write/READ_WRITE.h>
#include <Core/Vectors/VECTOR.h>
namespace PhysBAM{
//#####################################################################
// Function Bounding_Box
//#####################################################################
template<class T> RANGE<VECTOR<T,3> > RING<T>::
Bounding_Box() const
{
    return CYLINDER<T>(plane1.x0,plane2.x0,outer_radius).Bounding_Box();
}
//#####################################################################
// Function Signed_Distance
//#####################################################################
template<class T> T RING<T>::
Signed_Distance(const TV& X) const
{
    TV v=X-plane1.x0;
    T plane1_phi=v.Dot(plane1.normal),plane_phi=max(plane1_phi,-height-plane1_phi);
    T radius=(v-plane1_phi*plane1.normal).Magnitude();
    T radial_phi=max(radius-outer_radius,inner_radius-radius);
    return (radial_phi>0 && plane_phi>0)?VECTOR<T,2>(radial_phi,plane_phi).Magnitude():max(radial_phi,plane_phi);
}
//#####################################################################
// Function Suface
//#####################################################################
template<class T> VECTOR<T,3> RING<T>::
Surface(const TV& X) const 
{
    TV v=X-plane1.x0;
    T plane1_phi=v.Dot(plane1.normal),plane2_phi=-height-plane1_phi,plane_phi=max(plane1_phi,plane2_phi);
    TV radial_direction=v-plane1_phi*plane1.normal;
    T radius=radial_direction.Normalize();
    T outer_phi=radius-outer_radius,inner_phi=inner_radius-radius;
    T radial_phi=max(inner_phi,outer_phi);
    if(radial_phi>0 || radial_phi>plane_phi){ // closest point is on one of the cylinders
        T nearest_radius=outer_phi>inner_phi?outer_radius:inner_radius;
        return plane1.x0+clamp(plane1_phi,-height,(T)0)*plane1.normal+nearest_radius*radial_direction;}
    return plane1_phi > plane2_phi?X-plane1_phi*plane1.normal:X+plane2_phi*plane1.normal;
}
//#####################################################################
// Function Normal
//#####################################################################
template<class T> VECTOR<T,3> RING<T>::
Normal(const TV& X) const 
{
    // compute plane normal
    TV v=X-plane1.x0;
    T plane1_phi=v.Dot(plane1.normal),plane2_phi=-height-plane1_phi;
    TV plane_normal=plane1_phi>plane2_phi?plane1.normal:plane2.normal;
    T plane_phi=max(plane1_phi,plane2_phi);
    // compute radial normal
    TV radial_normal=v-plane1_phi*plane1.normal;
    T radius=radial_normal.Normalize();
    T outer_phi=radius-outer_radius,inner_phi=inner_radius-radius;
    if(inner_phi>outer_phi) radial_normal=-radial_normal;
    T radial_phi=max(inner_phi,outer_phi);
    // combine normals
    if(radial_phi>0 && plane_phi>0){
        VECTOR<T,2> weights=VECTOR<T,2>(radial_phi,plane_phi).Normalized();
        return weights.x*radial_normal+weights.y*plane_normal;}
    return radial_phi>plane_phi?radial_normal:plane_normal;
}
//#####################################################################
// Function Normal
//#####################################################################
template<class T> VECTOR<T,3> RING<T>::
Normal(const TV& X,const int aggregate) const 
{
    assert(aggregate>=0 && aggregate<4);
    if(aggregate==0){ // outer cylinder
        TV normal=X-plane1.x0;
        normal-=normal.Dot(plane1.normal)*plane1.normal;
        return normal.Normalized();}
    else if(aggregate==1){ // inner cylinder
        TV normal=X-plane1.x0;
        normal-=normal.Dot(plane1.normal)*plane1.normal;
        return -normal.Normalized();}
    else if(aggregate==2)
        return plane1.Normal();
    else 
        return plane2.Normal();
}
//#####################################################################
// Function Hessian
//#####################################################################
template<class T> SYMMETRIC_MATRIX<T,3> RING<T>::
Hessian(const TV& X) const
{
    AUTO_HESS<TV,TV> Y=AUTO_HESS<TV,TV>::From_Var(X),v=Y-plane1.x0;
    AUTO_HESS<T,TV> plane1_phi=v.Dot(plane1.normal),plane_phi=max(plane1_phi,-height-plane1_phi);
    AUTO_HESS<T,TV> radius=(v-plane1_phi*plane1.normal).Magnitude();
    AUTO_HESS<T,TV> radial_phi=max(radius-outer_radius,inner_radius-radius);
    return radial_phi.x>0 && plane_phi.x>0?sqrt(sqr(radial_phi)+sqr(plane_phi)).ddx:max(radial_phi,plane_phi).ddx;
}
//#####################################################################
// Function Fill_Uniform
//#####################################################################
template<class T> void RING<T>::
Fill_Uniform(RANDOM_NUMBERS<T>& rand,const RANGE<TV>& box)
{
    CYLINDER<T> c;
    rand.Fill_Uniform(c,box);
    T ir;
    rand.Fill_Uniform(ir,0,c.radius);
    *this=RING<T>(c.plane1.x0,c.plane2.x0,c.radius,ir);
}
//#####################################################################
// Function Lazy_Inside
//#####################################################################
template<class T> bool RING<T>::
Lazy_Inside(const TV& X) const 
{
    if(!plane1.Lazy_Inside(X) || !plane2.Lazy_Inside(X)) return false;
    TV v=X-plane1.x0;
    v-=v.Dot(plane1.normal)*plane1.normal;
    T radius_sqr=v.Magnitude_Squared();
    return radius_sqr<=sqr(outer_radius) && radius_sqr>=sqr(inner_radius);
}
//#####################################################################
// Function Lazy_Outside
//#####################################################################
template<class T> bool RING<T>::
Lazy_Outside(const TV& X) const  
{
    return !Lazy_Inside(X);
}
//#####################################################################
// Function Inside
//#####################################################################
template<class T> bool RING<T>::
Inside(const TV& X,const T thickness_over_two) const 
{
    if(plane1.Inside(X,thickness_over_two) && plane2.Inside(X,thickness_over_two)){ // inside both planes, check the cylinder
        TV distance=X-plane1.x0;
        distance-=distance.Dot(plane1.normal)*plane1.normal;
        T distance_squared = distance.Dot(distance);
        if(distance_squared <= sqr(outer_radius-thickness_over_two) && distance_squared >= sqr(inner_radius+thickness_over_two)) return true;}
    return false;
}
//#####################################################################
// Function Outside
//#####################################################################
template<class T> bool RING<T>::
Outside(const TV& X,const T thickness_over_two) const  
{
    if(plane1.Outside(X,thickness_over_two) || plane2.Outside(X,thickness_over_two)) return true;
    else{ // check the cylinder
        TV distance=X-plane1.x0;
        distance-=distance.Dot(plane1.normal)*plane1.normal;
        T distance_squared = distance.Dot(distance);
        if(distance_squared >= sqr(outer_radius+thickness_over_two) || distance_squared <= sqr(inner_radius-thickness_over_two)) return true;
        else return false;}
}
//#####################################################################
// Function Boundary
//#####################################################################
template<class T> bool RING<T>::
Boundary(const TV& X,const T thickness_over_two) const
{
    return !Inside(X,thickness_over_two) && !Outside(X,thickness_over_two);
}
//#####################################################################
// Function Name
//#####################################################################
template<class T> std::string RING<T>::
Name()
{
    return "RING<T>";
}
//#####################################################################
// Function Volume
//#####################################################################
template<class T> T RING<T>::
Volume() const
{
    return pi_v<T>*height*(sqr(outer_radius)-sqr(inner_radius));
}
//#####################################################################
// Function Scaled
//#####################################################################
template<class T,class TV> RING<T>
Scaled(const RING<T>& b,const TV& center,T scale)
{
    return RING<T>((b.plane1.x0-center)*scale+center,
        (b.plane2.x0-center)*scale+center,
        b.outer_radius*scale,b.inner_radius*scale);
}
//#####################################################################
// Function Transformed
//#####################################################################
template<class T,class TV> RING<T>
Transformed(const RING<T>& b,const FRAME<TV>& f)
{
    return RING<T>(f*b.plane1.x0,f*b.plane2.x0,b.outer_radius,b.inner_radius);
}
//#####################################################################
// Function Generate_Triangles
//#####################################################################
template<class T> TRIANGULATED_SURFACE<T>* Generate_Triangles(const RING<T>& ring,const int m,const int n,const int p)
{
    typedef VECTOR<T,3> TV;
    TRIANGULATED_SURFACE<T>* surface=TRIANGULATED_SURFACE<T>::Create();
    surface->mesh.Initialize_Torus_Mesh(2*(n+p),m);
    GEOMETRY_PARTICLES<TV>& particles=surface->particles;
    particles.Add_Elements(2*(n+p)*m);
    TV u=ring.plane1.normal.Unit_Orthogonal_Vector();
    TV v=ring.plane1.normal.Cross(u);
    TV dw=(ring.plane2.x0-ring.plane1.x0)/n;
    T dr=(ring.outer_radius-ring.inner_radius)/p;
    for(int i=0,q=0;i<m;++i){
        T t=2*pi_v<T>*i/m;
        TV radial=cos(t)*u+sin(t)*v;
        for(int j=0;j<p;j++) particles.X(q++)=ring.plane1.x0+(ring.inner_radius+j*dr)*radial;
        for(int j=0;j<n;j++) particles.X(q++)=ring.plane1.x0+j*dw+ring.outer_radius*radial;
        for(int j=0;j<p;j++) particles.X(q++)=ring.plane2.x0+(ring.outer_radius-j*dr)*radial;
        for(int j=0;j<n;j++) particles.X(q++)=ring.plane2.x0-j*dw+ring.inner_radius*radial;}
    return surface;
}
//#####################################################################
// Function Intersects
//#####################################################################
template<class T> RAY_HIT<T> Intersects(const RAY<VECTOR<T,3> >& ray,
    const RING<T>& ring,const T isocontour)
{
    return Ray_Hit(Intersection_Interval(ray,ring,isocontour));
}
//#####################################################################
// Function Intersection_Interval
//#####################################################################
template<class T> VECTOR<RAY_HIT<T>,8>
Intersection_Interval(const RAY<VECTOR<T,3> >& ray,const RING<T>& ring,
    T isocontour)
{
    CYLINDER<T> ci(ring.plane1.x0,ring.plane2.x0,ring.inner_radius),co(ci);
    co.radius=ring.outer_radius;

    auto i0=Intersection_Interval_Infinite(ray,co,isocontour);
    auto i1=Intersection_Interval_Infinite(ray,ci,isocontour);
    auto i2=Intersection_Interval(ray,ring.plane1,isocontour);
    auto i3=Intersection_Interval(ray,ring.plane2,isocontour);
    Set_Id(i0,0);
    Set_Id(i1,1);
    Set_Id(i2,2);
    Set_Id(i3,3);
    return Difference(Intersection(Intersection(i0,i2),i3),i1);
}
//#####################################################################
// Function All_Intersections
//#####################################################################
template<class TV,class T> void
All_Intersections(ARRAY<RAY_HIT<T> >& hits,const RAY<TV>& ray,
    const RING<T>& ring,const T isocontour)
{
    auto iv=Intersection_Interval(ray,ring,isocontour);
    hits=iv.Array_View(0,Actual_Length(iv));
}
//#####################################################################
// Function Read
//#####################################################################
template<class T> template<class RW> void RING<T>::
Read(UNTYPED_ISTREAM<RW>& input)
{
    Read_Binary_Group(input,plane1.x0,plane2.x0,outer_radius,inner_radius);
    Set_Endpoints(plane1.x0,plane2.x0);
}
//#####################################################################
// Function Write
//#####################################################################
template<class T> template<class RW> void RING<T>::
Write(UNTYPED_OSTREAM<RW>& output) const
{
    Write_Binary_Group(output,plane1.x0,plane2.x0,outer_radius,inner_radius);
}
//#####################################################################
template class RING<float>;
template class RING<double>;
template RING<double> Scaled<double,VECTOR<double,3> >(RING<double> const&,VECTOR<double,3> const&,double);
template RING<double> Transformed<double,VECTOR<double,3> >(RING<double> const&,FRAME<VECTOR<double,3> > const&);
template RING<float> Scaled<float,VECTOR<float,3> >(RING<float> const&,VECTOR<float,3> const&,float);
template RING<float> Transformed<float,VECTOR<float,3> >(RING<float> const&,FRAME<VECTOR<float,3> > const&);
template TRIANGULATED_SURFACE<float>* Generate_Triangles(const RING<float>&,const int,const int,const int);
template TRIANGULATED_SURFACE<double>* Generate_Triangles(const RING<double>&,const int,const int,const int);
template RAY_HIT<double> Intersects<double>(RAY<VECTOR<double,3> > const&,RING<double> const&,double);
template RAY_HIT<float> Intersects<float>(RAY<VECTOR<float,3> > const&,RING<float> const&,float);
template void All_Intersections<VECTOR<double,3>,double>(ARRAY<RAY_HIT<double>,int>&,RAY<VECTOR<double,3> > const&,RING<double> const&,double);
template void All_Intersections<VECTOR<float,3>,float>(ARRAY<RAY_HIT<float>,int>&,RAY<VECTOR<float,3> > const&,RING<float> const&,float);
template void RING<double>::Read(UNTYPED_ISTREAM<char>&);
template void RING<double>::Read(UNTYPED_ISTREAM<double>&);
template void RING<double>::Read(UNTYPED_ISTREAM<float>&);
template void RING<double>::Write(UNTYPED_OSTREAM<char>&) const;
template void RING<double>::Write(UNTYPED_OSTREAM<double>&) const;
template void RING<double>::Write(UNTYPED_OSTREAM<float>&) const;
template void RING<float>::Read(UNTYPED_ISTREAM<char>&);
template void RING<float>::Read(UNTYPED_ISTREAM<double>&);
template void RING<float>::Read(UNTYPED_ISTREAM<float>&);
template void RING<float>::Write(UNTYPED_OSTREAM<char>&) const;
template void RING<float>::Write(UNTYPED_OSTREAM<double>&) const;
template void RING<float>::Write(UNTYPED_OSTREAM<float>&) const;
}
