// Author : Luc Jaulin
    #include <list>
    #include "paving.h"

box(*F_sivia) (const box&);  // Global variable used in SIVIA;
paving X_sivia;


paving::paving()                                                      
    {           box p;
                B.clear(); left.clear();  right.clear();  val.clear();
                B.push_back(p);
                left.push_back(-1);
                right.push_back(-1);
                val.push_back(iperhaps);
    };



paving::paving(box& p,iboolean b)
    {           B.clear(); left.clear();  right.clear();  val.clear();
                B.push_back(p);
                left.push_back(-1);
                right.push_back(-1);
                val.push_back(b);
    };

//    paving::~paving ();

//----------------------------------------------------------------------
ostream& operator<<(ostream& os, const paving& X)
{   cout << "paving :"<<"\t dim="<<Size(X.B[0])<<"\n";
    /*for (int i=0;i<X.B.size();i++)
       {
	os <<i<< "\t : left= "<<X.left[i]<<",\t right= "<<X.right[i]<<",\t val= "<<X.val[i]<<", \t X=["<< X.B[i][1]<<","<< X.B[i][2]<<"] \n";
       }*/
       return (os);
}
//----------------------------------------------------------------------
paving& paving::Clear(iboolean b)        // After Clear(), the paving contains a single box.
    {  	box P;
        P=B[0];
        B.clear();
        left.clear();
        right.clear();
        val.clear();
        B.push_back(P);
        left.push_back(-1);
        right.push_back(-1);
        val.push_back(b);
        return (*this);
    }
//----------------------------------------------------------------------
paving& paving::Remove_sons(int i)
    {  	if (i<0) return (*this);
        Remove_sons(left[i]);
        Remove_sons(right[i]);
        if (left[i]!=-1)  val[left[i]]=empty;
        if (right[i]!=-1)  val[right[i]]=empty;
        left[i]=-1;         right[i]=-1;
        return (*this);
    }

//----------------------------------------------------------------------
paving& paving::Reunite(int i)
{  if ((i<0)||(left[i]==-1)) return (*this);
   Reunite(left[i]);
   Reunite(right[i]);
   if (val[right[i]]==val[left[i]])
            { val[i]=val[left[i]];
              if ((left[left[i]]==-1)&&(left[right[i]]==-1)) Remove_sons(i);
            };
   if (i==0) Clean();
   return (*this);
}
//----------------------------------------------------------------------
paving& paving::Clean()    //garbage collector
{  paving X(*this);
   Clear(iperhaps);
   list<int> LX;
   list<int> LY;
   LX.push_back(0);
   LY.push_back(0);
   while (!LX.empty())
     { int i=LX.back();
       LX.pop_back();
       int j=LY.back();
       LY.pop_back();
       val[j]=X.val[i];
       B[j]=X.B[i];
       if (X.left[i]!=-1)
           { LX.push_back(X.right[i]);
             LX.push_back(X.left[i]);
             Expand(j);
             LY.push_back(right[j]);
             LY.push_back(left[j]);
           }
     };
   return (*this);
}
//----------------------------------------------------------------------
paving& paving::Expand(int i)
    {  	box P1,P2;
        if (left[i]==-1)
        {   Bisect(B[i],P1,P2);
            int n=B.size();
            left[i]=n;              right[i]=n+1;
            left.push_back(-1);     right.push_back(-1);
            left.push_back(-1);     right.push_back(-1);
            val.push_back(val[i]);  val.push_back(val[i]);
            B.push_back(P1);        B.push_back(P2);
        }
        return (*this);
    }


//----------------------------------------------------------------------
iboolean  Inside(paving& Z, box& X,int k=0)   // returns 0, if outside, 1 if inside, empty if X is empty, and iperhaps otherwize
{  if (X.IsEmpty())
        return empty; // A cause de l'Union ci-dessous
   if (!Subset(X,Z.B[k]))
      return iperhaps; // n'arrive jamais, sauf pour k=0
   if (Z.val[k]!=iperhaps) return Z.val[k];
   int kleft=Z.left[k];
   int kright=Z.right[k];
   if ((kleft!=-1)||(kright!=-1))
   { box Xleft =X & Z.B[kleft];
     box Xright=X & Z.B[kright];
     iboolean I1= Inside(Z,Xleft,kleft);
     iboolean I2= Inside(Z,Xright,kright);

     return(Union(I1,I2));
   }
   return (iperhaps);
}
//----------------------------------------------------------------------
void Contract_distance_gt_paving(paving& A, double z, box& X,int k) //Contract X with respect to distance to a subpaving lower than z.
{  if (X.IsEmpty())  return;
   if (A.val[k]==ifalse)
       {X=EmptyBox(X);return;}                  // NE MARCHE PAS ! A REFAIRE !
   int kleft=A.left[k];
   int kright=A.right[k];
   if (((A.val[k]==iperhaps)&&(kleft==-1)&&(kright==-1))||(A.val[k]==itrue))
   {  interval X1=X[1];
      interval X2=X[2];
      interval A1=(A.B[k])[1];
      interval A2=(A.B[k])[2];
      interval Z(0,z*z);    //X1^2+X2^2=Z
      interval D1=X1-A1;
      interval D2=X2-A2;
      interval D12=Sqr(D1);
      interval D22=Sqr(D2);
      Cplus(Z,D12,D22,1);
      Cplus(Z,D12,D22,-1);
      Csqr(D22,D2,-1);
      Csqr(D12,D1,-1);
      Cmoins(D2,X2,A2,-1);
      Cmoins(D1,X1,A1,-1);
      X=box(X1,X2);
      return;
   }
     box X1(X);
     box X2(X);
     Contract_distance_gt_paving(A,z,X1,kleft);
     Contract_distance_gt_paving(A,z,X2,kright);
     X=Union(X1,X2);
     return;
}
//----------------------------------------------------------------------
// Operation between pavings
void  op_binaire(paving& Z, paving& X,paving& Y, iboolean(*op)(iboolean,iboolean),int k,int i,int j)
{  if (k==0) {Z.B[0]=X.B[0]; Z.Clear(iperhaps);};
   iboolean vali = op(X.val[i],Y.val[j]);
   Z.val[k]=vali;
   int ileft=X.left[i];
   int jleft=Y.left[j];
   int iright=X.right[i];
   int jright=Y.right[j];
   if (((ileft!=-1)||(jleft!=-1))&&(vali==iperhaps))
       { if (ileft==-1) {ileft=i; iright=i;}
         if (jleft==-1) {jleft=j; jright=j;}
         Z.Expand(k);
         op_binaire(Z,X,Y,op,Z.left[k],ileft,jleft);
         op_binaire(Z,X,Y,op,Z.right[k],iright,jright);
       }
   if (k==0) {Z.Clean();Z.Reunite();};
   return;
}
//----------------------------------------------------------------------
// Unary operation between pavings
void  op_unaire(paving& Z, paving& X, iboolean(*op)(iboolean),int k,int i)
{  if (k==0) {Z.B[0]=X.B[0]; Z.Clear(iperhaps);};
   iboolean vali = op(X.val[i]);
   Z.val[k]=vali;
   int ileft=X.left[i];
   int iright=X.right[i];
   if ((ileft!=-1)&&(vali==iperhaps))
       { Z.Expand(k);
         op_unaire(Z,X,op,Z.left[k],ileft);
         op_unaire(Z,X,op,Z.right[k],iright);
       }
   if (k==0) {Z.Clean();Z.Reunite();};
   return;
}


//----------------------------------------------------------------------
paving& paving::Sivia(iboolean(*test)(const box&),iboolean(*op)(iboolean,iboolean),double eps)
{   list<int> L;
    L.push_back(0);
    while (!L.empty())
     { int i=L.front();    L.pop_front();
       iboolean testBi=test(B[i]);
       iboolean vali = op(val[i],testBi);
       if (vali!=iperhaps)
         { Remove_sons(i);  }
       else
           if  ((B[i].Width()>eps) && (testBi==iperhaps))    //Luc
             {  Expand(i);
                L.push_back(left[i]);
                L.push_back(right[i]);
             }
       val[i]=vali;
     };
   Clean();
   Reunite();
   return (*this);
}

iboolean Test_sivia(const box& P)     //used only in Sivia below
{   box Y=F_sivia(P);
    return (Inside(X_sivia,Y,0));
}

//----------------------------------------------------------------------
paving& paving::Sivia(box(*F)(const box&),paving& X,iboolean(*op)(iboolean,iboolean),double eps)
{  X_sivia=X;
   F_sivia=F;
   Sivia(Test_sivia,op,eps);
}





//----------------------------------------------------------------------
paving  Sqcap(paving& X,paving& Y)
{    paving Z;
     op_binaire(Z,X,Y,Inter);
     return (Z);
}
//----------------------------------------------------------------------
void  Set_Contractor_Inter(paving& Z,paving& X,paving& Y)   //project with respect to Z = X inter Y
{    //op_binaire(Z,X,Y,Inter);  //Z=X inter Y

     //D= paving(X.B[0]);
     //op_binaire(Z,X,Y,Restrict);  //Z=X inter Y


     //paving Z1=Z;
//     op_binaire(D,0,X,0,Y,0,0);     // and
//     op_binaire(Z,0,Z1, 0,D,0,2);   // intersect
     //A FINIR
}

//----------------------------------------------------------------------
void  Set_Contractor_Overlap(paving& X,paving& Y)   //project with respect to X and Y overlap
{    paving XinterY;
     op_binaire(XinterY,X,Y,And);
     if (XinterY.val[0]==ifalse)
     {  X.Clear(empty);  Y.Clear(empty);}
}
//----------------------------------------------------------------------
iboolean Overlap(paving& X,paving& Y)   //returns if  X and Y overlap
{    paving XinterY;
     op_binaire(XinterY,X,Y,And);
     if (XinterY.val[0]==ifalse) return ifalse;
     for (int k=0;k<XinterY.B.size();k++)
        if (XinterY.val[0]==itrue) return itrue;
     return iperhaps;
}

//----------------------------------------------------------------------
void  Set_Contractor_Subset(paving& X,paving& Y)   //project with respect to X \subset Y
{    paving Z;
     // Solution 1
     //op_binaire(Z,X,Y,leq);  //X=X inter Y
     //X=Z;
     //op_binaire(Z,Y,X,geq);  //Y=X union Y
     //Y=Z ;
     //  ------
     // Solution 2
     op_binaire(Z,X,Y,And);  //X=X inter Y
     X=Sqcap(X,Z);
     op_binaire(Z,Y,X,Or);  //Y=X union Y
     Y=Sqcap(Y,Z);
}
//----------------------------------------------------------------------
void  Set_Contractor_Disjoint(paving& X,paving& Y)   //project with respect to X inter Y = empty
{    paving All(X.B[0],iperhaps);
     paving Z1;
     op_binaire(Z1,All,Y,Restrict);
     paving Z2;
     Z2=Sqcap(X,Z1);

     paving Z3;
     op_binaire(Z3,All,X,Restrict);
     paving Z4;
     Z4=Sqcap(Y,Z3);

     X=Z2; Y=Z4;
}

