include/realroot/fatarcs.hpp

Go to the documentation of this file.

#include <stack>
#include <fstream>
#include <string>
#include <vector>

#include <realroot/Interval.hpp>
#include <realroot/ring_bernstein_tensor.hpp>
#include <realroot/ring_monomial_tensor.hpp>
#include <realroot/tensor_convert.hpp>
#include <realroot/fatarcs_fcts.hpp>

#pragma once

#define REP    poly.rep()
#define SIZE   poly.size()
#define TODOUBLE as<double>

using namespace mmx;

//BOX structure
template<class C>
struct domain
  { 
    typedef Interval<C>         int_t;
    typedef Seq< Interval<C> >  seqint_t;
    typedef std::vector<C>      vec_t;

//CONSTRUCTORS

    seqint_t     I;
   

    domain( ){ }; 
    domain( int n ){ int_t unit( C(0),C(1) ); for(int i = 0; i < n; i++){I<<unit;}; }; 
    domain( int_t intxy ){ I<<intxy<<intxy; };
    domain( int_t intx , int_t inty ){ I<<intx<<inty; };
 
//FCTS
    C dim(){
      int d=0; 
      for(unsigned i = 0; i < I.size(); i++){d++;}
      return(d);
    }
 

    C diam(){
      C d=C(0); 
      for(unsigned i = 0; i < I.size(); i++){
          d +=I[i].size()*I[i].size();}
      return(sqrt(d));
    }
 
  vec_t llc( )
    {vec_t v(I.size());
      for(unsigned i = 0; i < I.size(); i++){
          v[i]= I[i].m;};
      return v;
    };

  vec_t urc( )
    { vec_t v(I.size());
      for(unsigned i = 0; i < I.size(); i++){
          v[i]= I[i].M;};
      return v;
    };

  vec_t center( )
    { vec_t v(I.size());
      for(unsigned i = 0; i < I.size(); i++){
          v[i]= I[i].center();};
      return v;
    };


  Seq<C> delta( )
    {Seq<C> d; 
      for(unsigned i = 0; i < I.size(); i++){
          d<<I[i].size();}
      return(d);
    };
 
  void print(int digits)
    { std::cout.precision(digits);
      std::cout<<"["<< TODOUBLE(I[0].m)<<", "<< TODOUBLE(I[0].M)<<"]x["<< TODOUBLE(I[1].m)<<", "<< TODOUBLE(I[1].M)<<"],"<<std::endl;
      // std::cout << "box midpt: (" <<(llc()[0]+urc()[0])/C(2)<<","<<(llc()[1]+urc()[1])/C(2)<<")"<<std::endl;
      std::cout << "diam: " <<diam()<<std::endl;
    };


  void split(domain<C> * ch)
    { 
     
      C cx=I[0].center() ;
      C cy=I[1].center() ;
     
      Interval<C> Ixl(I[0].m , cx), Ixr(cx , I[0].M), Iyl(I[1].m , cy), Iyr(cy, I[1].M);

     
      ch[0]=domain(Ixl,Iyl);
      ch[1]=domain(Ixr,Iyl);
      ch[2]=domain(Ixl,Iyr);
      ch[3]=domain(Ixr,Iyr);
    };/*subdivison of a box*/
    /*
Seq< vec_t > cornerpts( )
    { Seq< vec_t > cp;
      cp<<vec(I[0].m,I[1].m)<<vec(I[0].M,I[1].m)<<vec(I[0].m,I[1].M)<<vec(I[0].M,I[1].M);
     return cp;
     };*/

};//end of boxes


//arc_rep structure

template<class C>
struct arc_rep
{

  typedef C                                              coeff_t;
  typedef std::vector<coeff_t>                           vec_t;
  typedef Seq<coeff_t>                                   seq_t;
  typedef domain<coeff_t>                                box_t;
  typedef polynomial<  coeff_t ,with<Bernstein>  >       bernstein_t;
  typedef polynomial<  coeff_t ,with<MonomialTensor> >   monom_t;
 
 

   vec_t pts[3];
    
  //CONSTRUCTORS

  arc_rep(){pts[0]=vec_t(); pts[1]=vec_t(); pts[2]=vec_t();};

  arc_rep( vec_t p0,  vec_t p1, vec_t p2){
    if(p0.size()!=0 && p1.size()!=0 &&  p2.size()!=0 )
      {if(matrat(v_minus(p1,p0),v_minus(p2,p0))[1]> coeff_t(0)){pts[0]=p2; pts[1]=p1; pts[2]=p0;}
        else{pts[0]=p0; pts[1]=p1; pts[2]=p2;}; }else{pts[0]=p0; pts[1]=p1; pts[2]=p2;}
  };

  arc_rep( vec_t p[3]){
    if(p[0].size()!=0 && p[1].size()!=0 &&  p[2].size()!=0 ){
      if(matrat( v_minus(p[1],p[0]),v_minus(p[2],p[0]) )[1]> coeff_t(0)){pts[0]=p[2]; pts[1]=p[1]; pts[2]=p[0];}
      else{pts[0]=p[0]; pts[1]=p[1]; pts[2]=p[2];}; }else{pts[0]=p[0]; pts[1]=p[1]; pts[2]=p[2];} 
  };
 
  //FCTS

 
  coeff_t weight()
  {return(dot(v_minus(pts[0],pts[1]),v_minus(pts[1],pts[2]))/
          (v_length(v_minus(pts[1],pts[2]))*v_length(v_minus(pts[1],pts[0]))));};
  
 
  vec_t midc(){   
    vec_t pp=v_minus(pts[0],pts[2]); pp=rotl(pp);
    vec_t m=v_mul(coeff_t( sign(dot(pp, v_minus(pts[1],pts[2]) ))),pp);
    coeff_t w=weight(); 
  
    vec_t d;
    if( (coeff_t(1)-w*w) > coeff_t(0) && w!=coeff_t(0) ){
      d=v_plus(v_div(v_plus(pts[0],pts[2]),coeff_t(2)),v_div(v_mul(coeff_t(sqrt( coeff_t(1)-w*w )),m),(coeff_t(2)*w))); }
    else {
      d=v_div(v_plus(pts[0],pts[2]),coeff_t(2)); };
    return(d);

  };/*middle control point in bb rep.*/



  seq_t arc_eq(){ 

  
    vec_t cvec=matrat(v_minus(pts[0],pts[1]),v_minus(pts[2],pts[1]));  
     
    coeff_t eh2=cvec[1];
    coeff_t ehx=dot(cvec, vec(pts[2][1]-pts[0][1],(coeff_t)(-1)*pts[0][0]-pts[2][0]));
    coeff_t ehy=dot(cvec, vec(pts[0][0]-pts[2][0],(coeff_t)(-1)*pts[0][1]-pts[2][1]));  
    coeff_t ehc=dot(cvec, vec(pts[0][1]*pts[2][0]-pts[0][0]*pts[2][1],
                                         pts[0][0]*pts[2][0]+pts[0][1]*pts[2][1])); 

    seq_t eh; eh<<eh2<<ehx<<ehy<<ehc;
     
  return(eh);
  }; /*coeff vector of monom rep: eh2*(x^2+y^2)+ehx*x+ehy*y+ehc=0*/

  

 vec_t critpt(int var){ 

  vec_t pt;
   Seq<coeff_t> coeffx,coeffy,wvec; 
   pt=vec(coeff_t(-1),coeff_t(-1));
   coeffx<<pts[0][0]<< midc()[0]*weight()<< pts[2][0];
   coeffy<<pts[0][1]<< midc()[1]*weight()<<pts[2][1]; 
   wvec<<coeff_t(1)<< weight()<<coeff_t(1);  
   coeff_t t, cx,cy,w;

   if (var==0 && (coeffx[2]-coeffx[0])!=0)
     {
       t=(coeffx[0]-coeffx[1])/(coeffx[2]-coeffx[0]);
     }
  else if (var==1 && (coeffy[2]-coeffy[0])!=0)
     {
       t=(coeffy[0]-coeffy[1])/(coeffy[2]-coeffy[0]);
     };

   if(is_unit1(t)){
   tensor::eval(cx,seq2b(coeffx).rep(),t);
   tensor::eval(cy,seq2b(coeffy).rep(),t);
   tensor::eval(w,seq2b(wvec).rep(),t);
   pt=vec( cx/w , cy/w );   
   }
   //  std::cout <<pt<<std::endl;
  return(pt);
  }; /*tang. pt in var*/


 arc_rep<coeff_t> offset(coeff_t r){
   vec_t newp[3];
   Seq<coeff_t> mo, ehrat;
   int i;unsigned j;
   vec_t pp=rotl(v_minus(pts[0],pts[2])); 

    coeff_t l=v_length(pp);
    coeff_t w=weight();

    for( i=0; i<3; i++){newp[i]=vec_t(); };
    arc_rep<coeff_t> ncirc(newp);
    if( l!=coeff_t(0) ){  
      if( w==coeff_t(1)){ for( i=0; i<3; i++){ newp[i]=v_plus(pts[i],v_mul(r/l,pp)); } }//case of a quarter of a circle
      else if (  
               //( coeff_t(1)-w*w )>coeff_t(0) && 
               ( r + l/(coeff_t(2)*sqrt( coeff_t(1)-w*w ) )) > coeff_t(0)  
            &&  v_length( v_minus(pts[1],v_plus(pts[0],pts[2])) )!=coeff_t(0) )
        { 

            newp[0]=v_plus(pts[0],v_mul(coeff_t(2)*r/l,
                                rotl( v_minus(v_minus(v_mul((coeff_t(1)+w),pts[1]),v_mul((w+ coeff_t(0.5)),pts[0])),
                                                v_div(pts[2],coeff_t(2))
                                                )
                                      )
                                  )
                     );

            newp[1]=v_plus(pts[1],
                           v_mul(r/v_length(v_minus(pts[1],v_div(v_plus(pts[0],pts[2]),coeff_t(2)))),
                           v_minus(v_minus(pts[1],v_div(pts[0],coeff_t(2))),v_div(pts[2],coeff_t(2)))
                           )
                     );

            newp[2]=v_minus(pts[2],
                      v_mul(coeff_t(2)*r/l,
                            rotl(v_minus(v_minus(v_mul((coeff_t(1)+w),pts[1]),v_mul((w+(coeff_t)(0.5)),pts[2])),
                                         v_div(pts[0],coeff_t(2))
                                         )
                                 )
                            )
                      );
        };
   
      //std::cout <<"newp0: "<<newp[0]<<" "<<newp[1]<<" "<<newp[2]<<std::endl;
     ncirc=arc_rep<coeff_t>(newp);
     if(is_arc(ncirc)){
     Seq<vec_t> isp;
    ehrat=ncirc.arc_eq();  
    mo=solve2(ehrat[0],ehrat[2],ehrat[3]); 
    for( j=0; j<2; j++){if(is_unit1(mo[j])){isp<<vec(coeff_t(0),mo[j]);};};
    mo=solve2(ehrat[0],ehrat[2],ehrat[0]+ehrat[1]+ehrat[3]);
    for( j=0; j<2; j++){if(is_unit1(mo[j])){isp<<vec(coeff_t(1),mo[j]);};};
    mo=solve2(ehrat[0],ehrat[1],ehrat[3]);
    for( j=0; j<2; j++){if(is_unit1(mo[j])){isp<<vec(mo[j],coeff_t(0));};};
    mo=solve2(ehrat[0],ehrat[1],ehrat[0]+ehrat[2]+ehrat[3]);
    for( j=0; j<2; j++){if(is_unit1(mo[j])){isp<<vec(mo[j],coeff_t(1));};};

    if (isp.size()==0){newp[0]=vec_t(); }else{
    newp[0]=isp[0];
    for( j=1; j<isp.size(); j++)
      {if(v_length(v_minus(isp[j],pts[0]))<v_length(v_minus(newp[0],pts[0]))){newp[0]=isp[j];};};
    };

    if (isp.size()==0){newp[2]=vec_t(); }else{
    newp[2]=isp[0];
    for( j=1; j<isp.size(); j++)
      {if(v_length(v_minus(isp[j],pts[2]))<v_length(v_minus(newp[2],pts[2]))){newp[2]=isp[j];};};
    };
   ncirc=arc_rep<coeff_t>(newp);
   // std::cout <<"newp1: "<<newp[0]<<" "<<newp[1]<<" "<<newp[2]<<std::endl;
   /*if (is_arc(ncirc)){
      if( matrat(v_minus(newp[1],newp[0]),v_minus(newp[2],newp[0]))[1]>coeff_t(0) )
      { vec_t u=newp[0]; newp[0]=newp[2]; newp[2]=u;}};*/
     };
  };
 return( ncirc );
 } /*offset arc computation(with precise endpoints) - OR OUTSIDE!!!*/

};//end of struct arc_rep


//box_rep structure


template<class C>
struct box_rep
{

  typedef C                                              coeff_t;
  typedef std::vector<coeff_t>                           vec_t;
  typedef domain<coeff_t>                                box_t;
  typedef polynomial< coeff_t ,with<Bernstein> >         bernstein_t; 
  typedef polynomial< coeff_t ,with<MonomialTensor> >    monom_t;  
  typedef arc_rep<coeff_t>               arc_rep_t;
 
  
  box_t       box;
  bernstein_t poly;

  //CONSTRUCTORS

  box_rep(bernstein_t  p){box_t bb; poly=p; box=bb; };
  box_rep(bernstein_t  p, box_t bb)
  {
    bernstein_t pa=p;

    tensor::restrict(pa.rep(), 0 , bb.I[0].m , bb.I[0].M );
    tensor::restrict(pa.rep(), 1 , bb.I[1].m , bb.I[1].M );
   
    poly=pa;
    box=bb;
  };


  //FCTS

bool not_empty(){

  if( sign(min_coeff(poly))*sign(max_coeff(poly)) > 0 ){
    return false;
  }
  else{ 
    return true; 
  }
}/*detect empty boxes*/


  bool edge_event(int var , int n){

  bernstein_t face_p;
  int event = 0;

  tensor::face(face_p, poly, var, n);
  
  int stc = 0;
 
  for(unsigned i = 0 ; i < face_p.size() ; i++){
    if( stc*sign(face_p.rep()[i]) == -1 ){
      stc = sign(face_p.rep()[i]); event++;
    }
    else if ( stc == 0 && sign(face_p.rep()[i])!=0 ){
      stc = sign(face_p.rep()[i]);
    };
  }
  if(event==1){return true;}else{return false;};
}/*edge event along ONE edge*/


int first_nonzero(bernstein_t list, int e){

  int elem =0;
  int n = list.size();
  int j;

  for(int i = 0; i < n ; i++){
    if( e==0 ){ j=i; }else{j=n-i-1;}
    if( coeff_t(0)!=(list[j]) ){
      elem=sign(list[j]);
      break;
    }
  } 
  return elem;

}


  bool corner_event(int c1, int c2){
  vec_t evc=vec(coeff_t(c1),coeff_t(c2));
  coeff_t v;
  eval(v, poly.rep(), evc );
  bernstein_t face_p, face_q;
 
  if ( is_small(v) ){      
    tensor::face(face_p, poly, 0, c2);
    tensor::face(face_q, poly, 1, c1);

    if (first_nonzero(face_p, c1)*first_nonzero(face_q, c2) < 0){return true; }else{ return false; };

  }else{ return false;};
}


  Seq<vec_t>event_list(){ 
    Seq<vec_t> elist;
    if(edge_event( 0, 0)){elist<<vec(coeff_t(0),coeff_t(0))<<vec(coeff_t(1),coeff_t(0));}
    if(edge_event( 0, 1)){elist<<vec(coeff_t(0),coeff_t(1))<<vec(coeff_t(1),coeff_t(1));}
    if(edge_event( 1, 0)){elist<<vec(coeff_t(0),coeff_t(0))<<vec(coeff_t(0),coeff_t(1));}
    if(edge_event( 1, 1)){elist<<vec(coeff_t(1),coeff_t(0))<<vec(coeff_t(1),coeff_t(1));}

    if(corner_event( 0, 0)){elist<<vec(coeff_t(0),coeff_t(0))<<vec_t();}
    if(corner_event( 0, 1)){elist<<vec(coeff_t(0),coeff_t(1))<<vec_t();}
    if(corner_event( 1, 0)){elist<<vec(coeff_t(1),coeff_t(0))<<vec_t();}
    if(corner_event( 1, 1)){elist<<vec(coeff_t(1),coeff_t(1))<<vec_t();}
    return elist;

} //boxes with single implicitly defined curve segment*/


coeff_t min_grad(){

  bernstein_t p0 = diff(poly,0);
  bernstein_t p1 = diff(poly,1);
  bernstein_t grad;
  bernstein_t gradx0 = p0; 
  bernstein_t gradx1 = p1;

 
  mul(gradx0.rep(),p0.rep());
  mul(gradx1.rep(),p1.rep()); 
  add(grad.rep(),gradx0.rep(),gradx1.rep());

  return(min_coeff(grad));
}


coeff_t max_eval(arc_rep_t circ){
  bernstein_t neww, bxt, byt, bwt, num, polyco;
  monom_t monomp, newnom, xt, yt, wt;    
  Seq<coeff_t> v;

  v<<coeff_t(1)<<circ.weight()<<coeff_t(1);
  neww=seq2b( v );
  v.clear();
 
    pow(neww,degree(poly,0)+degree(poly,1));

    v<<circ.pts[0][0]<<circ.midc()[0]*circ.weight()<<circ.pts[2][0];
    bxt=seq2b(v); v.clear(); 
    tensor::assign(xt.rep(),bxt.rep());
    v<<circ.pts[0][1]<<circ.midc()[1]*circ.weight()<<circ.pts[2][1];
    byt=seq2b(v); v.clear(); 
    tensor::assign(yt.rep(),byt.rep());
    v<<coeff_t(1)<<circ.weight()<<coeff_t(1);
    bwt=seq2b(v); v.clear(); 
    tensor::assign(wt.rep(),bwt.rep());
  
    Seq<monom_t> paramv; paramv<<wt<<xt<<yt;
 
    polyco=poly;
    tensor::assign(monomp.rep(),polyco.rep());
    tensor::heval(newnom, monomp.rep(), paramv ); 
  
    let::assign(num.rep(),newnom.rep());
      
   return(abs_max_coeff(num)/min_coeff(neww));

};


};//end of structure box_rep

---