//#####################################################################
// Copyright 2002-2007, Doug Enright, Ronald Fedkiw, Frederic Gibou, Geoffrey Irving, Sergey Koltakov, Neil Molino, Igor Neverov, Andrew Selle, 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.
//#####################################################################
#include <Core/Level_Sets/LEVELSET.h>
#include <Core/Arrays/FACE_ARRAYS.h>
#include <Core/Basic_Geometry/RAY.h>
#include <Core/Basic_Geometry/RAY_HIT.h>
#include <Core/Boundaries/BOUNDARY.h>
#include <Core/Grids/CELL_ITERATOR.h>
#include <Core/Grids/FACE_ITERATOR.h>
#include <Core/Grids/NODE_ITERATOR.h>
#include <Core/Implicit_Objects/IMPLICIT_OBJECT.h>
#include <Core/Interpolation/MULTILINEAR_INTERPOLATION.h>
#include <Core/Level_Sets/FAST_MARCHING_METHOD_UNIFORM.h>
#include <Core/Level_Sets/LEVELSET_UTILITIES.h>
#include <Core/Level_Sets/Visit_Ray_Grid.h>
#include <Core/Math_Tools/INTERVAL.h>
#include <Core/Nonlinear_Equations/ITERATIVE_SOLVER.h>
#include <Core/Polynomials/CUBIC.h>
#include <Core/Polynomials/QUADRATIC.h>
#include <Core/Read_Write/READ_WRITE.h>
#include <Core/Vectors/ZERO_VECTOR.h>
namespace PhysBAM{
template<class TV> LINEAR_INTERPOLATION_NODE<TV,typename TV::SCALAR> LEVELSET<TV>::interpolation_default;
template<class TV> LINEAR_INTERPOLATION_NODE<TV,TV> LEVELSET<TV>::normal_interpolation_default;
//#####################################################################
// Constructor
//#####################################################################
template<class TV> LEVELSET<TV>::
LEVELSET(const int number_of_ghost_cells_input)
    :LEVELSET(*new ARRAY<T,IV>,number_of_ghost_cells_input)
{
    owns_phi=true;
}
//#####################################################################
// Constructor
//#####################################################################
template<class TV> LEVELSET<TV>::
LEVELSET(ARRAY<T,IV>& phi_input,const int number_of_ghost_cells_input)
    :levelset_callbacks(0),face_velocities_valid_mask_current(0),boundary_default(*new BOUNDARY<TV,T>),
    phi(phi_input),normals(0),curvature(0),cell_range(0),number_of_ghost_cells(number_of_ghost_cells_input)
{
    Set_Small_Number();
    Set_Max_Time_Step();
    boundary=&boundary_default;
    interpolation=&interpolation_default;
    curvature_interpolation=&interpolation_default;
    normal_interpolation=&normal_interpolation_default;
    secondary_interpolation=0;
    Set_Band_Width();
}
//#####################################################################
// Destructor
//#####################################################################
template<class TV> LEVELSET<TV>::
~LEVELSET()
{
    if(owns_phi) delete &phi;
    delete &boundary_default;
    delete normals;
    delete curvature;
    delete cell_range;
}
//#####################################################################
// Function CFL
//#####################################################################
template<class TV> typename TV::SCALAR LEVELSET<TV>::
CFL(const ARRAY<T,FACE_INDEX<TV::m> >& face_velocities) const
{
    T dt_convection=0;
    for(CELL_ITERATOR<TV> iterator(grid);iterator.Valid();iterator.Next()){IV cell=iterator.Cell_Index();
        T local_V_norm=0;
        for(int axis=0;axis<TV::m;axis++)
            local_V_norm+=grid.one_over_dX[axis]*maxabs(face_velocities(axis,grid.First_Face_Index_In_Cell(axis,cell)),face_velocities(axis,grid.Second_Face_Index_In_Cell(axis,cell)));
        dt_convection=max(dt_convection,local_V_norm);}
    return 1/max(dt_convection,1/max_time_step);
}
//#####################################################################
// Function CFL
//#####################################################################
template<class TV> typename TV::SCALAR LEVELSET<TV>::
CFL(const ARRAY<TV,IV>& velocity) const
{
    T dt_convection=0;
    for(CELL_ITERATOR<TV> iterator(grid);iterator.Valid();iterator.Next()){IV cell=iterator.Cell_Index();
        dt_convection=max(dt_convection,velocity(cell).Dot(grid.one_over_dX));}
    return 1/max(dt_convection,1/max_time_step);
}
//#####################################################################
// Function Iterative_Find_Interface
//#####################################################################
template<class TV> TV LEVELSET<TV>::
Iterative_Find_Interface(TV left,TV right,const int iterations) const
{
    T phi_left=Phi(left),phi_right=Phi(right);
    TV interface=Linear(left,right,LEVELSET_UTILITIES<T>::Theta(phi_left,phi_right));
    int phi_left_sign=(phi_left<=0?-1:1),phi_right_sign=(phi_right<=0?-1:1);
    for(int i=0;i<iterations;i++){
        T phi=Phi(interface);
        int phi_sign=(phi<=0?-1:1);
        if(phi_left_sign*phi_sign<0){
            right=interface;phi_right=phi;phi_right_sign=phi_sign;
            interface=Linear(left,interface,LEVELSET_UTILITIES<T>::Theta(phi_left,phi));}
        else if(phi_right_sign*phi_sign<0){
            left=interface;phi_left=phi;phi_left_sign=phi_sign;
            interface=Linear(interface,right,LEVELSET_UTILITIES<T>::Theta(phi,phi_right));}
        else break;}
    return interface;
}
//#####################################################################
// Function Compute_Gradient
//#####################################################################
template<class TV> void LEVELSET<TV>::
Compute_Gradient(ARRAY<TV,IV>& gradient,const T time) const
{
    TV one_over_two_dx=(T).5*grid.one_over_dX;
    int ghost_cells=3;
    ARRAY<T,IV> phi_ghost(grid.Node_Indices(ghost_cells));
    boundary->Fill_Ghost_Nodes(grid,phi,phi_ghost,ghost_cells);
    gradient.Resize(grid.Node_Indices(ghost_cells-1));
    for(CELL_ITERATOR<TV> iterator(grid,ghost_cells-1);iterator.Valid();iterator.Next()){
        IV cell_index=iterator.Cell_Index();
        for(int axis=0;axis<TV::m;axis++){
            IV axis_vector=IV::Axis_Vector(axis);
            gradient(cell_index)(axis)=(phi_ghost(cell_index+axis_vector)-phi_ghost(cell_index-axis_vector))*one_over_two_dx(axis);}}
}
//#####################################################################
// Function Gradient
//#####################################################################
template<class TV> TV LEVELSET<TV>::
Gradient(const TV& location) const
{
    TV G;
    for(int d=0;d<TV::m;d++){
        TV a(location),b(location);
        a(d)-=grid.dX(d);
        b(d)+=grid.dX(d);
        G(d)=(T).5*(Phi(b)-Phi(a))*grid.one_over_dX(d);}
    return G;
}
//#####################################################################
// Function Normal
//#####################################################################
template<class TV> TV LEVELSET<TV>::
Normal(const TV& location) const
{
    if(normals) return normal_interpolation->Interpolate_Node(grid,*normals,location).Normalized();
    return Gradient(location).Normalized();
}
//#####################################################################
// Function Extended_Normal
//#####################################################################
template<class TV> TV LEVELSET<TV>::
Extended_Normal(const TV& location) const
{
    if(normals) return normal_interpolation->Interpolate_Node(grid,*normals,grid.domain.Clamp(location)).Normalized();
    TV N;
    for(int d=0;d<TV::m;d++){
        TV a(location),b(location);
        a(d)-=grid.dX(d);
        b(d)+=grid.dX(d);
        N(d)=(T).5*(Extended_Phi(b)-Extended_Phi(a))*grid.one_over_dX(d);}
    return N.Normalized();
}
//#####################################################################
// Function Compute_Normals
//#####################################################################
// note that sqrt(phix^2+phiy^2+phiz^2)=1 if it's a distance function
template<class TV> void LEVELSET<TV>::
Compute_Normals(const T time)
{
    int ghost_cells=3;
    ARRAY<T,IV> phi_ghost(grid.Node_Indices(ghost_cells));
    boundary->Fill_Ghost_Nodes(grid,phi,phi_ghost,ghost_cells);

    if(!normals) normals=new ARRAY<TV,IV>(grid.Node_Indices(ghost_cells-1));
    for(CELL_ITERATOR<TV> iterator(grid,ghost_cells-1);iterator.Valid();iterator.Next()){
        const IV& cell = iterator.Cell_Index();
        int index=phi_ghost.Standard_Index(iterator.Cell_Index());
        TV& N((*normals)(cell));
        for(int d=0;d<TV::m;d++)
            N(d)=(T).5*(phi_ghost.array(index+phi_ghost.stride(d))-phi_ghost.array(index-phi_ghost.stride(d)))*grid.one_over_dX(d);
        N.Normalize();}
}
//#####################################################################
// Function Hessian
//#####################################################################
template<class TV> SYMMETRIC_MATRIX<typename TV::SCALAR,TV::m> LEVELSET<TV>::
Hessian(const TV& X) const
{
    SYMMETRIC_MATRIX<T,TV::m> H;
    for(int i=0;i<TV::m;i++){
        TV A(X),B(X);
        A(i)-=grid.dX(i);
        B(i)+=grid.dX(i);
        H(i,i)=(Phi(B)-2*Phi(X)+Phi(A))*sqr(grid.one_over_dX(i));}

    for(int i=0;i<TV::m;i++){
        TV A(X),B(X);
        A(i)-=grid.dX(i);
        B(i)+=grid.dX(i);
        for(int j=i+1;j<TV::m;j++){
            TV C(A),D(A),E(B),F(B);
            C(j)-=grid.dX(j);
            D(j)+=grid.dX(j);
            E(j)-=grid.dX(j);
            F(j)+=grid.dX(j);
            H(i,j)=(Phi(F)-Phi(E)-Phi(D)+Phi(A))*(T).25*grid.one_over_dX(i)*grid.one_over_dX(j);}}

    return H;
}
//#####################################################################
// Function Compute_Cell_Minimum_And_Maximum
//#####################################################################
template<class TV> void LEVELSET<TV>::
Compute_Cell_Minimum_And_Maximum(const bool recompute_if_exists)
{
    if(!recompute_if_exists && cell_range) return;
    RANGE<IV> pruned_range(phi.domain);
    pruned_range.max_corner-=1;
    if(!cell_range) cell_range=new ARRAY<INTERVAL<T>,IV>(pruned_range);
    for(RANGE_ITERATOR<TV::m> it(pruned_range);it.Valid();it.Next()){
        int index=phi.Standard_Index(it.index);
        INTERVAL<T> interval(phi.array(index));
        for(int i=1;i<(1<<TV::m);i++){
            int ind=index;
            for(int d=0;d<TV::m;d++)
                ind+=((i>>d)&1)*phi.stride(d);
            interval.Enlarge_To_Include_Point(phi.array(ind));}
        (*cell_range)(it.index)=interval;}
}
//#####################################################################
// Function Hessian
//#####################################################################
template<class TV> SYMMETRIC_MATRIX<typename TV::SCALAR,TV::m> LEVELSET<TV>::
Hessian(const ARRAY<T,IV>& phi_input,const IV& index) const
{
    int s_index=phi_input.Standard_Index(index);
    SYMMETRIC_MATRIX<T,TV::m> ddphi;
    for(int a=0;a<TV::m;a++)
        for(int b=0;b<a;b++)
            ddphi(a,b)=(T).25*(phi_input.array(s_index+phi_input.stride(a)+phi_input.stride(b))-phi_input.array(s_index+phi_input.stride(a)-phi_input.stride(b))
                -phi_input.array(s_index-phi_input.stride(a)+phi_input.stride(b))+phi_input.array(s_index-phi_input.stride(a)-phi_input.stride(b)))*grid.one_over_dX(a)*grid.one_over_dX(b);

    for(int a=0;a<TV::m;a++)
        ddphi(a,a)=(phi_input.array(s_index+phi_input.stride(a))-2*phi_input.array(s_index)+phi_input.array(s_index-phi_input.stride(a)))*sqr(grid.one_over_dX(a));
    return ddphi;
}
//#####################################################################
// Function Gradient
//#####################################################################
template<class TV> TV LEVELSET<TV>::
Gradient(const ARRAY<T,IV>& phi_input,const IV& index) const
{
    int s_index=phi_input.Standard_Index(index);
    TV dphi;
    for(int a=0;a<TV::m;a++)
        dphi(a)=(T).5*(phi_input.array(s_index+phi_input.stride(a))-phi_input.array(s_index-phi_input.stride(a)))*grid.one_over_dX(a);
    return dphi;
}
//#####################################################################
// Function Compute_Curvature
//#####################################################################
template<class TV> typename TV::SCALAR LEVELSET<TV>::
Compute_Curvature(const ARRAY<T,IV>& phi_input,const IV& index) const
{
    TV dphi=Gradient(phi_input,index);
    SYMMETRIC_MATRIX<T,TV::m> ddphi=Hessian(phi_input,index);

    T norm_squared=dphi.Magnitude_Squared(),norm=sqrt(norm_squared);
    if(norm<small_number) return LEVELSET_UTILITIES<T>::Sign(phi_input(index))/grid.dX.Min();
    T curvature=(dphi.Dot(ddphi*dphi)-norm_squared*ddphi.Trace())/(norm*norm_squared);
    return minmag(curvature,LEVELSET_UTILITIES<T>::Sign(curvature)/grid.dX.Min());
}
//#####################################################################
// Function Compute_Curvature
//#####################################################################
// kappa = - DIV(normal), negative for negative phi inside, positive for positive phi inside, sqrt(phix^2+phiy^2+phiy^2)=1 for distance functions
template<class TV> void LEVELSET<TV>::
Compute_Curvature(const T time)
{
    int ghost_cells=3;
    ARRAY<T,IV> phi_ghost(grid.Node_Indices(ghost_cells));
    boundary->Fill_Ghost_Nodes(grid,phi,phi_ghost,ghost_cells);

    if(!curvature) curvature=new ARRAY<T,IV>(grid.Node_Indices(ghost_cells-1));
    for(CELL_ITERATOR<TV> iterator(grid,ghost_cells-1);iterator.Valid();iterator.Next())
        (*curvature)(iterator.Cell_Index())=Compute_Curvature(phi_ghost,iterator.Cell_Index());
}
//#####################################################################
// Function Compute_Curvature
//#####################################################################
template<class TV> typename TV::SCALAR LEVELSET<TV>::
Compute_Curvature(const TV& location) const
{
    TV l2=(location-(T).5*grid.dX-grid.domain.min_corner)*grid.one_over_dX;
    IV cell(floor(l2));
    TV w(l2-TV(cell));
    if(!phi.Valid_Index(cell-1)) return 0;
    if(!phi.Valid_Index(cell+2)) return 0;

    T k[1<<TV::m];
    for(int i=0;i<(1<<TV::m);i++){
        IV index(cell);
        for(int d=0;d<TV::m;d++)
            if(i&(1<<d))
                index(d)++;
        k[i]=Compute_Curvature(phi,index);}

    return Linear(k,w);
}
//#####################################################################
// Function Fast_Marching_Method
//#####################################################################
template<class TV> void LEVELSET<TV>::
Fast_Marching_Method(const T time,const T stopping_distance,const ARRAY<IV>* seed_indices,const bool add_seed_indices_for_ghost_cells,int process_sign)
{
    Get_Signed_Distance_Using_FMM(phi,time,stopping_distance,seed_indices,add_seed_indices_for_ghost_cells,process_sign);
}
//#####################################################################
// Function Get_Signed_Distance_Using_FMM
//#####################################################################
template<class TV> void LEVELSET<TV>::
Get_Signed_Distance_Using_FMM(ARRAY<T,IV>& signed_distance,const T time,const T stopping_distance,const ARRAY<IV>* seed_indices,const bool add_seed_indices_for_ghost_cells,int process_sign)
{
    const int ghost_cells=2*number_of_ghost_cells+1;
    ARRAY<T,IV> phi_ghost(grid.Node_Indices(ghost_cells),no_init);
    boundary->Fill_Ghost_Nodes(grid,phi,phi_ghost,ghost_cells);
    FAST_MARCHING_METHOD_UNIFORM<TV> fmm;
    fmm.seed_indices=seed_indices;
    if(seed_indices) fmm.correct_interface_phi=false;
    fmm.add_seed_indices_for_ghost_cells=add_seed_indices_for_ghost_cells;
    fmm.process_sign=process_sign;
    fmm.Fast_Marching_Method(grid,at_nodes,ghost_cells,phi_ghost,stopping_distance);

    ARRAY<T,IV>::Get(signed_distance,phi_ghost);
    boundary->Apply_Boundary_Condition_Node(grid,signed_distance);
}
//#####################################################################
// Function Approximate_Surface_Size
//#####################################################################
// calculates the approximate perimeter using delta functions
template<class TV> typename TV::SCALAR LEVELSET<TV>::
Approximate_Surface_Size(const T interface_thickness,const T time) const
{
    ARRAY<T,IV> phi_ghost(grid.Node_Indices(number_of_ghost_cells));
    boundary->Fill_Ghost_Nodes(grid,phi,phi_ghost,number_of_ghost_cells);
    T interface_half_width=interface_thickness*grid.dX.Max()/2,surface_size=0;
    for(RANGE_ITERATOR<TV::m> it(grid.Node_Indices());it.Valid();it.Next())
        surface_size+=LEVELSET_UTILITIES<T>::Delta(phi_ghost(it.index),interface_half_width)*Gradient(phi,it.index).Magnitude();
    return surface_size*grid.dX.Size();
}
//#####################################################################
// Function Extended_Phi
//#####################################################################
template<class TV> auto LEVELSET<TV>::
Extended_Phi(const TV& location) const -> T
{
    TV clamped_location(grid.domain.Clamp(location));
    T phi_value=interpolation->Interpolate_Node(grid,phi,clamped_location);
    T magnitude_squared=(clamped_location-location).Magnitude_Squared();
    if(magnitude_squared) phi_value=sqrt(magnitude_squared+sqr(max((T)0,phi_value)));
    return phi_value;
}
//#####################################################################
// Function Lazy_Inside
//#####################################################################
template<class TV> bool LEVELSET<TV>::
Lazy_Inside(const TV& clamped_X,const T contour_value) const
{
    assert(cell_range);
    IV index=grid.Node_LC(clamped_X,number_of_ghost_cells,true);
    
    
    const INTERVAL<T>& range=(*cell_range)(index);
    if(range.min_corner>contour_value) return false;
    else if(range.max_corner<=contour_value) return true;
    return interpolation->Interpolate_Node(grid,phi,clamped_X)<=contour_value;
}
//#####################################################################
// Function Lazy_Inside_And_Value
//#####################################################################
template<class TV> bool LEVELSET<TV>::
Lazy_Inside_And_Value(const TV& clamped_X,T& phi_value,const T contour_value) const
{
    assert(cell_range);
    IV index=grid.Node_LC(clamped_X,number_of_ghost_cells,true);
    if((*cell_range)(index).min_corner>contour_value) return false;
    phi_value=interpolation->Interpolate_Node(grid,phi,clamped_X);
    return phi_value<=contour_value;
}
//#####################################################################
// Function Lazy_Inside_Extended_Levelset
//#####################################################################
template<class TV> bool LEVELSET<TV>::
Lazy_Inside_Extended_Levelset(const TV& unclamped_X,const T contour_value) const
{
    if(grid.domain.Lazy_Inside_Half_Open(unclamped_X))
        return Lazy_Inside(unclamped_X,contour_value);
    if(contour_value<=0) return false;
    T sqr_contour_value=sqr(contour_value);
    TV clamped_X(grid.domain.Clamp(unclamped_X));
    T magnitude_squared=(unclamped_X-clamped_X).Magnitude_Squared();
    if(magnitude_squared>sqr_contour_value) return false;
    IV index=grid.Node_LC(clamped_X,number_of_ghost_cells,true);
    T phi_value=interpolation->Interpolate_Node(grid,phi,clamped_X);
    return magnitude_squared+sqr(max((T)0,phi_value))<=sqr_contour_value;
}
//#####################################################################
// Function Lazy_Inside_Extended_Levelset_And_Value
//#####################################################################
template<class TV> bool LEVELSET<TV>::
Lazy_Inside_Extended_Levelset_And_Value(const TV& unclamped_X,T& phi_value,
    const T contour_value) const
{
    assert(cell_range);
    TV clamped_X(grid.domain.Clamp(unclamped_X));
    T magnitude_squared=(unclamped_X-clamped_X).Magnitude_Squared();
    if(contour_value <= 0){if(magnitude_squared > 0) return false;}
    else if(magnitude_squared > sqr(contour_value)) return false;
    IV index=grid.Node_LC(clamped_X,number_of_ghost_cells,true);
    if((*cell_range)(index).min_corner>contour_value) return false;
    phi_value=interpolation->Interpolate_Node(grid,phi,clamped_X);
    if(magnitude_squared>0) phi_value=sqrt(magnitude_squared+sqr(max((T)0,phi_value)));
    return phi_value<=contour_value;
}
//#####################################################################
// Function Lazy_Outside_Extended_Levelset_And_Value
//#####################################################################
template<class TV> bool LEVELSET<TV>::
Lazy_Outside_Extended_Levelset_And_Value(const TV& unclamped_X,T& phi_value,
    const T contour_value) const
{
    assert(cell_range);
    TV clamped_X(grid.domain.Clamp(unclamped_X));
    T magnitude_squared=(unclamped_X-clamped_X).Magnitude_Squared();
    IV index=grid.Node_LC(clamped_X,number_of_ghost_cells,true);
    if(magnitude_squared==0 && (*cell_range)(index).max_corner<=contour_value) return false;
    phi_value=interpolation->Interpolate_Node(grid,phi,clamped_X);
    if(magnitude_squared>0) phi_value=sqrt(magnitude_squared+sqr(max((T)0,phi_value)));
    return phi_value>contour_value;
}
//#####################################################################
// Function Read
//#####################################################################
template<class TV> template<class RW> void LEVELSET<TV>::
Read(UNTYPED_ISTREAM<RW>& input)
{
    Read_Binary_Group(input,grid,phi);
}
//#####################################################################
// Function Write
//#####################################################################
template<class TV> template<class RW> void LEVELSET<TV>::
Write(UNTYPED_OSTREAM<RW>& output) const
{
    Write_Binary_Group(output,grid,phi);
}
//#####################################################################
// Function All_Intersections
//#####################################################################
template<class TV,class T> int
All_Intersections(const TV& E,const TV& D,T tmax,const VECTOR<T,(1<<TV::m)>& phi,T* ts)
{
    if constexpr(TV::m==1)
    {
        T rex=E.x,rdx=D.x;
        T p0=phi(0),p1=phi(1)-p0,x0=p1*rex+p0,a=p1*rdx;
        if(x0==0)
        {
            *ts=0;
            return 1;
        }
        if((a*tmax+x0<0)!=(x0<0))
        {
            *ts=-x0/a;
            return 1;
        }
        return 0;
    }
    else if constexpr(TV::m==2)
    {
        T rex=E.x,rdx=D.x;
        T rey=E.y,rdy=D.y;
        T p0=phi(0),p1=phi(1)-p0,p2=phi(2)-p0,p3=phi(3)-phi(2)-p1;
        T y0=p2*rey+p0,y1=p3*rey+p1,x0=y1*rex+y0,x2=p3*rex+p2;
        QUADRATIC<T> quad(p3*rdy*rdx,y1*rdx+x2*rdy,x0);
        quad.Compute_Roots_In_Interval(0,tmax);
        for(int i=0;i<quad.roots;i++)
            ts[i]=quad.root[i];
        return quad.roots;
    }
    else
    {
        T rex=E.x,rdx=D.x;
        T rey=E.y,rdy=D.y;
        T rez=E.z,rdz=D.z;
        T p0=phi(0),p1=phi(1)-p0,p2=phi(2)-p0,p4=phi(4)-p0,rdyz=rdy*rdz;
        T p54=phi(5)-phi(4),p5=p54-p1,p3=phi(3)-phi(2)-p1,p6=phi(6)-phi(4)-p2;
        T p7=phi(7)-phi(6)-p54-p3;
        T z0=p4*rez+p0,z1=p5*rez+p1,z2=p6*rez+p2,z3=p7*rez+p3;
        T y0=z2*rey+z0,y1=z3*rey+z1,y4=p6*rey+p4,y5=p7*rey+p5;
        T x0=y1*rex+y0,x2=z3*rex+z2,x4=y5*rex+y4,x6=p7*rex+p6;
        CUBIC<T> cubic(p7*rdx*rdyz,(z3*rdy+y5*rdz)*rdx+x6*rdyz,
            y1*rdx+x2*rdy+x4*rdz,x0);
        cubic.Compute_Roots_In_Interval(0,tmax);
        for(int i=0;i<cubic.roots;i++)
            ts[i]=cubic.root[i];
        return cubic.roots;
    }
}
//#####################################################################
// Function All_Intersections
//#####################################################################
template<class TV,class T> int
All_Intersections(const RAY<TV>& ray,const GRID<TV>& grid,
    const VECTOR<int,TV::m>& cell,const VECTOR<T,(1<<TV::m)>& phi,T* ts)
{
    T t0,t1;
    auto domain=grid.Cell_Domain(cell);
    if(!Get_Intersection_Range(ray,domain,t0,t1))
        return {};
    TV dir=ray.direction*grid.one_over_dX;
    TV A=(ray.Point(t0)-domain.min_corner)*grid.one_over_dX;
    int num_int=All_Intersections(A,dir,t1-t0,phi,ts);
    for(int i=0;i<num_int;i++)
        ts[i]=t0+ts[i];
    return num_int;
}
//#####################################################################
// Function Intersects
//#####################################################################
template<class TV> RAY_HIT<typename TV::SCALAR>
Intersects(const RAY<TV>& ray,const LEVELSET<TV>& levelset,typename TV::SCALAR isocontour)
{
    assert(levelset.cell_range);
    typedef typename TV::SCALAR T;
    typedef VECTOR<int,TV::m> IV;
    RAY_HIT<T> hit;
    Visit_Ray_Grid(ray,levelset.Native_Grid(),levelset.number_of_ghost_cells,
        [&,isocontour](const IV& cell,T t0,T t1)
        {
            if(!(*levelset.cell_range)(cell).Contains(isocontour)) return true;
            VECTOR<T,(1<<TV::m)> phi;
            for(int i=0;i<phi.m;i++)
                phi(i)=levelset.phi(cell+Binary_Counts<TV::m>(i));
            T ts[TV::m]={};
            int num_int=All_Intersections(ray,levelset.Native_Grid(),cell,phi-isocontour,ts);
            if(num_int==0) return true;
            hit.Set(ts[0]);
            return false;
        });
    return hit;
}
//#####################################################################
// Function All_Intersections
//#####################################################################
template<class TV,class T> void
All_Intersections(ARRAY<RAY_HIT<T> >& hits,const RAY<TV>& ray,
    const LEVELSET<TV>& levelset,const T isocontour)
{
    assert(levelset.cell_range);
    typedef VECTOR<int,TV::m> IV;
    Visit_Ray_Grid(ray,levelset.Native_Grid(),levelset.number_of_ghost_cells,
        [&,isocontour](const IV& cell,T t0,T t1)
        {
            if(!(*levelset.cell_range)(cell).Contains(isocontour)) return true;
            VECTOR<T,(1<<TV::m)> phi;
            for(int i=0;i<phi.m;i++)
                phi(i)=levelset.phi(cell+Binary_Counts<TV::m>(i));
            T ts[TV::m]={};
            int num_int=All_Intersections(ray,levelset.Native_Grid(),cell,phi-isocontour,ts);
            for(int i=0;i<num_int;i++)
                hits.Append(RAY_HIT<T>{ts[i]});
            return false;
        });
}
template class LEVELSET<VECTOR<float,1> >;
template class LEVELSET<VECTOR<float,2> >;
template class LEVELSET<VECTOR<float,3> >;
template class LEVELSET<VECTOR<double,1> >;
template class LEVELSET<VECTOR<double,2> >;
template class LEVELSET<VECTOR<double,3> >;
template void LEVELSET<VECTOR<double,1> >::Read(UNTYPED_ISTREAM<double>&);
template void LEVELSET<VECTOR<double,1> >::Read(UNTYPED_ISTREAM<float>&);
template void LEVELSET<VECTOR<double,1> >::Read(UNTYPED_ISTREAM<char>&);
template void LEVELSET<VECTOR<double,1> >::Write(UNTYPED_OSTREAM<double>&) const;
template void LEVELSET<VECTOR<double,1> >::Write(UNTYPED_OSTREAM<float>&) const;
template void LEVELSET<VECTOR<double,1> >::Write(UNTYPED_OSTREAM<char>&) const;
template void LEVELSET<VECTOR<double,2> >::Read(UNTYPED_ISTREAM<double>&);
template void LEVELSET<VECTOR<double,2> >::Read(UNTYPED_ISTREAM<float>&);
template void LEVELSET<VECTOR<double,2> >::Read(UNTYPED_ISTREAM<char>&);
template void LEVELSET<VECTOR<double,2> >::Write(UNTYPED_OSTREAM<double>&) const;
template void LEVELSET<VECTOR<double,2> >::Write(UNTYPED_OSTREAM<float>&) const;
template void LEVELSET<VECTOR<double,2> >::Write(UNTYPED_OSTREAM<char>&) const;
template void LEVELSET<VECTOR<double,3> >::Read(UNTYPED_ISTREAM<double>&);
template void LEVELSET<VECTOR<double,3> >::Read(UNTYPED_ISTREAM<float>&);
template void LEVELSET<VECTOR<double,3> >::Read(UNTYPED_ISTREAM<char>&);
template void LEVELSET<VECTOR<double,3> >::Write(UNTYPED_OSTREAM<double>&) const;
template void LEVELSET<VECTOR<double,3> >::Write(UNTYPED_OSTREAM<float>&) const;
template void LEVELSET<VECTOR<double,3> >::Write(UNTYPED_OSTREAM<char>&) const;
template void LEVELSET<VECTOR<float,1> >::Read(UNTYPED_ISTREAM<double>&);
template void LEVELSET<VECTOR<float,1> >::Read(UNTYPED_ISTREAM<float>&);
template void LEVELSET<VECTOR<float,1> >::Read(UNTYPED_ISTREAM<char>&);
template void LEVELSET<VECTOR<float,1> >::Write(UNTYPED_OSTREAM<double>&) const;
template void LEVELSET<VECTOR<float,1> >::Write(UNTYPED_OSTREAM<float>&) const;
template void LEVELSET<VECTOR<float,1> >::Write(UNTYPED_OSTREAM<char>&) const;
template void LEVELSET<VECTOR<float,2> >::Read(UNTYPED_ISTREAM<double>&);
template void LEVELSET<VECTOR<float,2> >::Read(UNTYPED_ISTREAM<float>&);
template void LEVELSET<VECTOR<float,2> >::Read(UNTYPED_ISTREAM<char>&);
template void LEVELSET<VECTOR<float,2> >::Write(UNTYPED_OSTREAM<double>&) const;
template void LEVELSET<VECTOR<float,2> >::Write(UNTYPED_OSTREAM<float>&) const;
template void LEVELSET<VECTOR<float,2> >::Write(UNTYPED_OSTREAM<char>&) const;
template void LEVELSET<VECTOR<float,3> >::Read(UNTYPED_ISTREAM<double>&);
template void LEVELSET<VECTOR<float,3> >::Read(UNTYPED_ISTREAM<float>&);
template void LEVELSET<VECTOR<float,3> >::Read(UNTYPED_ISTREAM<char>&);
template void LEVELSET<VECTOR<float,3> >::Write(UNTYPED_OSTREAM<double>&) const;
template void LEVELSET<VECTOR<float,3> >::Write(UNTYPED_OSTREAM<float>&) const;
template void LEVELSET<VECTOR<float,3> >::Write(UNTYPED_OSTREAM<char>&) const;
template RAY_HIT<double> Intersects<VECTOR<double,1> >(
    RAY<VECTOR<double,1> > const&,LEVELSET<VECTOR<double,1> > const&,double);
template RAY_HIT<double> Intersects<VECTOR<double,2> >(
    RAY<VECTOR<double,2> > const&,LEVELSET<VECTOR<double,2> > const&,double);
template RAY_HIT<double> Intersects<VECTOR<double,3> >(
    RAY<VECTOR<double,3> > const&,LEVELSET<VECTOR<double,3> > const&,double);
template RAY_HIT<float> Intersects<VECTOR<float,1> >(
    RAY<VECTOR<float,1> > const&,LEVELSET<VECTOR<float,1> > const&,float);
template RAY_HIT<float> Intersects<VECTOR<float,2> >(
    RAY<VECTOR<float,2> > const&,LEVELSET<VECTOR<float,2> > const&,float);
template RAY_HIT<float> Intersects<VECTOR<float,3> >(
    RAY<VECTOR<float,3> > const&,LEVELSET<VECTOR<float,3> > const&,float);
template void All_Intersections<VECTOR<double,1>,double>(ARRAY<RAY_HIT<double>,int>&,RAY<VECTOR<double,1> > const&,LEVELSET<VECTOR<double,1> > const&,double);
template void All_Intersections<VECTOR<double,2>,double>(ARRAY<RAY_HIT<double>,int>&,RAY<VECTOR<double,2> > const&,LEVELSET<VECTOR<double,2> > const&,double);
template void All_Intersections<VECTOR<double,3>,double>(ARRAY<RAY_HIT<double>,int>&,RAY<VECTOR<double,3> > const&,LEVELSET<VECTOR<double,3> > const&,double);
template void All_Intersections<VECTOR<float,1>,float>(ARRAY<RAY_HIT<float>,int>&,RAY<VECTOR<float,1> > const&,LEVELSET<VECTOR<float,1> > const&,float);
template void All_Intersections<VECTOR<float,2>,float>(ARRAY<RAY_HIT<float>,int>&,RAY<VECTOR<float,2> > const&,LEVELSET<VECTOR<float,2> > const&,float);
template void All_Intersections<VECTOR<float,3>,float>(ARRAY<RAY_HIT<float>,int>&,RAY<VECTOR<float,3> > const&,LEVELSET<VECTOR<float,3> > const&,float);
}
