fatarcs_fcts.hpp

Go to the documentation of this file.

#include <vector>
#include <realroot/Interval.hpp>
#include <realroot/Seq.hpp>
#define SMALL  C(10e-20)
#define REP    poly.rep()
#define SIZE   poly.size()

#pragma once

//#include <realroot/fatarcs.hpp>

template<class C>
struct arc_rep;
template<class C>
struct domain;
template<class C>
struct box_rep;

using namespace mmx;

//vec_t fcts

template<class C>
std::vector<C> vec (const C c1, const C c2) {
  std::vector<C> a(2, C(0));
   a[0]= c1; a[1]= c2;
  return a;
}

template<class C>
std::vector<C> rotl(std::vector<C> p){
  std::vector<C> v; v=vec(C(-1)*p[1],p[0]); return(v);
};       /*rotate left*/

template<class C>
C dot(std::vector<C> p ,std::vector<C> q){
  C v=C(0);
 for(unsigned i=0; i<p.size(); i++){
   v +=p[i]*q[i];
 } 
return(v);
};

template<class C>
std::vector<C> v_plus(std::vector<C> p ,std::vector<C> q){
  std::vector<C> v(p.size(),C(0));
 assert( p.size()==q.size() );
 for(unsigned i=0; i<p.size(); i++){
   v[i]=p[i]+q[i];
 } 
return(v);
};

template<class C>
std::vector<C> v_minus(std::vector<C> p ,std::vector<C> q){
  std::vector<C> v(p.size(),C(0));
  assert( p.size()==q.size() );
 for(unsigned i=0; i<p.size(); i++){
   v[i]=p[i]-q[i];
 } 
return(v);
};

template<class C>
std::vector<C> v_div(std::vector<C> p ,C q){
  if (q!=C(0)){
    std::vector<C>  v(p.size(),C(0));
    for(unsigned i=0; i<p.size(); i++){
      v[i]=p[i]/q;}; 
    return(v);}else{return std::vector<C>();}
};

template<class C>
std::vector<C> v_mul(C q, std::vector<C> p){
  std::vector<C>  v(p.size(),C(0));
 for(unsigned i=0; i<p.size(); i++){
   v[i]=p[i]*q;
 } 
return(v);
};

template<class C>
C v_length(std::vector<C> p){
  return(sqrt(dot(p,p)));
};   /*length of vector*/

template<class C>
std::vector<C> matrat(std::vector<C> v1,std::vector<C> v2){
  std::vector<C> mm; mm=vec(v1[0]*v2[0]+v1[1]*v2[1],v1[0]*v2[1]-v2[0]*v1[1]); 
  return(mm);
};

template<class C>
std::vector<C> crossrat(std::vector<C> v1,std::vector<C> v2,std::vector<C> v3,std::vector<C> v4){ 
  std::vector<C> mm=matrat(matrat(v_minus(v1,v2),v_minus(v3,v2)),matrat(v_minus(v1,v4),v_minus(v3,v4))); 
return(mm);
};


//solver and test fcts
template<class C>
bool is_zero(C p){
  if(p==C(0)){return true;}
  else{return false;}
  };

template<class C>
bool is_small(C p){
  if(p==C(0)){return true;}
  else{return false;}
  };

/*
template<class C>
bool is_small(C p){
  if(p<SMALL && p>C(-1)*SMALL){return true;}
  else{return false;}
  };
*/

template<class C>
Seq<C> solve2(C a, C b, C c){
  Seq<C> r; r<<C(-1)<<C(-1);  
  C root1,root2;  
  if(is_small(a)){if(!is_small(b)){r[0]=C(-1)*c/b;r[1]=C(-1)*c/b;}}
  else{if(sign(b*b-C(4)*a*c) >= 0){
      root1=(C(-1)*b+sqrt(b*b-C(4)*a*c))/(C(2)*a);
      root2=(C(-1)*b-sqrt(b*b-C(4)*a*c))/(C(2)*a);
      r[0]=root1;r[1]=root2;}
    }
 
  return(r);
};/*solve(a*x^2+b*x+c=0)*/

template<class C>
bool is_unit1(C p){
  if(p<C(0) || p>C(1) ){return false;}
  else{return true;}
};/*1dim unit test*/


template<class C>
bool is_unit2( std::vector<C> p){
  if( is_unit1( p[0] ) && is_unit1( p[1] ) ){return true;}
  else{return false;}
};/*2dim unit test*/


//bezier poly comp. fcts

int binomial(int n, int p) {
   int n1=1;
   if(p!=0) {

     if((n-p)<p){ p=n-p;};
     int n2=n;
     for(int i=1;i<=p;i++){ n1=n1*n2/i; n2--;};
    }
 
    return(n1);
}/*binomials*/


template<class C>
polynomial< C ,with<Bernstein> > bbasis(C co, int n, int deg)
  {
    polynomial< C ,with<Bernstein> > ti(co*C(binomial(n,deg)),deg,0), t("1-x0"), tj(C(1),0,0); 

    for(int i=0; i<n-deg; i++){ mul(tj,t); }; 

    mul(ti, tj);
    return(ti);
  }/*BB basis polynomials*/


template<class C>
void pow( polynomial< C ,with<Bernstein> > & poly, int n)
  {
    polynomial< C ,with<Bernstein> > t=poly;

    for(int i=1; i<n; i++){ mul(poly, t); }; 
  }/*nth power of a bbpoly*/



template<class C>
polynomial< C ,with<Bernstein> >  seq2b( Seq<C>  coeffseq )
  { 
    int n=coeffseq.size()-1;
    polynomial<  C ,with<Bernstein> > pol, pi;

    for(int i=0; i<=n; i++ ){
      pi=bbasis(coeffseq[i],n,i); 
      add( pol , pi);}

    return(pol);
  };/*coeff vector -> bbpoly*/


template<class C>
C pint( polynomial<  C ,with<MonomialTensor> >  poly)
  {
    int n=SIZE;
    C res=C(0);     

    for(int i=0; i<n; i++ ){
      res += REP[i]/C(i+1) ;}

    return(res);
  };/*poly integration on (0,1)*/

    
template<class C>
C min_coeff( polynomial< C ,with<Bernstein> > poly){

  C  min =REP[0];

  for ( unsigned i = 1; i < SIZE; i++ ){
    if( min > REP[i] ){ 
      min = REP[i];
    }

  }
  return min;
};/*minimal coeff*/


template<class C>
C max_coeff( polynomial< C ,with<Bernstein> > poly){

  C max =REP[0];

  for ( unsigned i = 1; i < SIZE; i++ ){
    if( max < REP[i] ){ 
      max = REP[i];
    }

  }
  return max;
};/*max coeff*/



template<class C>
C abs_max_coeff( polynomial< C ,with<Bernstein> > poly){
  C min =min_coeff(poly);
  C max =max_coeff(poly);

  if( max*max < min*min){ max =C(-1)*min; };
 
  return max;
};/* max abs coeff*/



template<class C>
Seq<C> corner_values(polynomial< C ,with<Bernstein> > p){
  C c00,c10,c01,c11;

  eval( c00, p.rep(),  vec(C(0),C(0)) );
  eval( c10, p.rep(),  vec(C(1),C(0)) );
  eval( c01, p.rep(),  vec(C(0),C(1)) );
  eval( c11, p.rep(),  vec(C(1),C(1)) );
  Seq<C> cval; cval<<c00<<c10<<c01<<c11;
  return(cval);

};

//arc fcts

template<class C>
std::vector<C> approx(polynomial< C ,with<Bernstein> > poly , std::vector<C> ep1, std::vector<C> ep2, int MTH){
  polynomial< C ,with<MonomialTensor> > xt, yt, res, monp, mp("x0-x0^2");
  polynomial< C ,with<Bernstein> > bxt, byt; 
    C t, I; 
    Seq<C> sols;
    std::vector<C> movec;
   
    if (ep2.size()==0 ){movec=ep1;}
    else{ 
 tensor::assign(monp.rep(),poly.rep());
   
 Seq<C> endptsx,endptsy;

    endptsx<<ep2[0]<<ep1[0];
    endptsy<<ep2[1]<<ep1[1];
    bxt=seq2b(endptsx); tensor::assign(xt.rep(),bxt.rep());
    byt=seq2b(endptsy); tensor::assign(yt.rep(),byt.rep());
       
    Seq<polynomial< C ,with<MonomialTensor> > > param; param<<xt<<yt;

    tensor::eval(res, monp.rep(), param, 1);
  

   
  C c0,c1;
  tensor::eval(c0,res.rep(),C(0));
  tensor::eval(c1,res.rep(),C(1));
  mul(res,mp);

   I=pint(res);
   C u=C(15)*I-C(0.75)*(c0+c1);
   sols=solve2(c0+c1-C(2)*u,C(2)*(u-c0),c0);
 

   if(is_unit1(sols[0]))
     {t=sols[0]; movec=vec(ep1[0]*t+ep2[0]*(C(1)-t),ep1[1]*t+ep2[1]*(C(1)-t));}
   else if(is_unit1(sols[1]))
     {t=sols[1]; movec=vec(ep1[0]*t+ep2[0]*(C(1)-t),ep1[1]*t+ep2[1]*(C(1)-t));}
   else{movec=std::vector<C>();}

    }
 
   return movec;
  }; /*approx bivar bbpoly along line*/


template<class C>
arc_rep< C >  median( box_rep< C > box, Seq<std::vector<C> > list, int MTH ){
  std::vector<C> endp[2], mpts[3];
  int s=0;
  arc_rep<C> medc;
  for(int i=0; i<3; i++){mpts[i]=std::vector<C>();};

  
  mpts[0]=approx(box.poly, list[0], list[1], MTH);
  mpts[1]=approx(box.poly, list[2], list[3], MTH);
 
    if( mpts[0].size()!=0 && mpts[1].size()!=0 ){
     
      std::vector<C> mid, rvec, v; 
      mid=v_div(v_plus(mpts[0],mpts[1]),C(2));
      rvec=rotl(v_minus(mpts[0],mpts[1])); 
       C bval[4]; 
    
       if(!is_small(rvec[0])){  
     bval[0]=rvec[1]*C(-1)*mid[0]/rvec[0]+mid[1];
     bval[1]=rvec[1]*(C(1)-mid[0])/rvec[0]+mid[1];}else{ bval[0]=C(-1); bval[1]=C(-1);}

       if(!is_small(rvec[1])){ 
     bval[2]=rvec[0]*C(-1)*mid[1]/rvec[1]+mid[0];
     bval[3]=rvec[0]*(C(1)-mid[1])/rvec[1]+mid[0];}else{ bval[2]=C(-1); bval[3]=C(-1);}
   
       if(is_unit1(bval[0]) && s<2){endp[s]= vec(C(0),bval[0]); s++;};
       if(is_unit1(bval[1]) && s<2){endp[s]= vec(C(1),bval[1]); s++;};
       if(is_unit1(bval[2]) && s<2){endp[s]= vec(bval[2],C(0)); s++;};
       if(is_unit1(bval[3]) && s<2){endp[s]= vec(bval[3],C(1)); s++;};
       mpts[2]=mpts[1]; 
  
       if(s==2){ mpts[1]=approx(box.poly, endp[0], endp[1], MTH);}

       if(mpts[1].size()!=0){ medc=arc_rep<C>(mpts);}else{medc=arc_rep<C>();}
       //else{ for(int i=0; i<3; i++){mpts[i]= vec(C(-1),C(-1));};}
      
    }else{ medc=arc_rep<C>(); };


    //std::cout <<"med: "<<mpts[0]<<" "<<mpts[1]<<" "<<mpts[2]<<std::endl;
   
     return(medc);
     }; /*median arc computation*/


  template<class C>
  void circ_is( std::vector<C> * ispts, arc_rep<C> c1, arc_rep<C> c2){

    Seq<C> eh1= c1.arc_eq(), eh2= c2.arc_eq(), sx, a;    
    C Aeh,Ceh; 
   if(!is_small(eh1[1]) && !is_small(eh1[2]) && !is_small(eh2[1]) && !is_small(eh2[2])){

    if(is_small(eh1[0]) && is_small(eh2[0])){

      Aeh=(C(-1)*eh2[1])/(eh2[2]); 
      Ceh=(C(-1)*eh2[3])/(eh2[2]);
      sx[0]=(C(-1)*eh1[2]*Ceh-eh1[3])/(eh1[2]*Aeh+eh1[1]);
        
          ispts[0]=vec(sx[0],Aeh*sx[0]+Ceh);
          ispts[1]=vec(sx[0],Aeh*sx[0]+Ceh);
      }
    else{
      if(!is_small(eh1[0]) && !is_small(eh2[0])){
    eh1=eh1/eh1[0];
    eh2=eh2/eh2[0];

    Aeh=(eh2[1]-eh1[1])/(eh1[2]-eh2[2]);
    Ceh=(eh2[3]-eh1[3])/(eh1[2]-eh2[2]);}   
      else{
    if(is_small(eh1[0]) && !is_small(eh2[0])){a=eh1; eh1=eh2; eh2=a;};
      
    eh1=eh1/eh1[0];

    Aeh=(C(-1)*eh2[1])/(eh2[2]); 
        Ceh=(C(-1)*eh2[3])/(eh2[2]);};
 


      sx=solve2(C(1)+Aeh*Aeh, C(2)*Aeh*Ceh+eh1[1]+eh1[2]*Aeh, Ceh*Ceh+eh1[2]*Ceh+eh1[3]);
     
      ispts[0]=vec(sx[0],Aeh*sx[0]+Ceh);
      ispts[1]=vec(sx[1],Aeh*sx[1]+Ceh); 
    
    }
   }
   else{  ispts[0]=vec(C(-1),C(-1));
     ispts[1]=vec(C(-1),C(-1)); 
   };
};/*cicrcle intersections (numerical behaviour!)*/


  template<class C>  
 bool is_arc(arc_rep<C> c){
    if( c.pts[0].size()!=0 && c.pts[1].size()!=0 && c.pts[2].size()!=0 ){return true;}
  else{return false;}
};


template<class C> 
 bool is_infatarc(std::vector<C> p, arc_rep<C> c1, arc_rep<C> c2){
  if(is_arc(c1)&&is_arc(c2)){
    if( crossrat(c1.pts[2],c1.pts[1],c1.pts[0],p)[1]*crossrat(c2.pts[2],c2.pts[1],c2.pts[0],p)[1]<C(0))
      {return true;}else{return false;}
  }
  else{return false;}
};


template<class C> 
Seq< std::vector<C> > extpts(arc_rep<C> mc1, C r1, arc_rep<C> mc2, C r2){

  arc_rep<C>
  c1p=mc1.offset(r1),c1m=mc1.offset(C(-1)*r1),
  c2p=mc2.offset(r2),c2m=mc2.offset(C(-1)*r2);
  std::vector<C> p[2];
  Seq< std::vector<C> > ispts;

  if(is_arc(c1p) && is_arc(c2p)){ circ_is(p,c1p,c2p); if(is_unit2(p[0])){ispts<<p[0];}; if(is_unit2(p[1])){ispts<<p[1];};};
  if(is_arc(c1p) && is_arc(c2m)){ circ_is(p,c1p,c2m); if(is_unit2(p[0])){ispts<<p[0];}; if(is_unit2(p[1])){ispts<<p[1];};};
  if(is_arc(c1m) && is_arc(c2m)){ circ_is(p,c1m,c2m); if(is_unit2(p[0])){ispts<<p[0];}; if(is_unit2(p[1])){ispts<<p[1];};};
  if(is_arc(c1m) && is_arc(c2p)){ circ_is(p,c1m,c2p); if(is_unit2(p[0])){ispts<<p[0];}; if(is_unit2(p[1])){ispts<<p[1];};};
 
  if(is_arc(c1p) && is_infatarc(c1p.pts[0],c2m,c2p)){ispts<<c1p.pts[0];};
  if(is_arc(c1p) && is_infatarc(c1p.pts[2],c2m,c2p)){ispts<<c1p.pts[2];};
  if(is_arc(c1m) && is_infatarc(c1m.pts[0],c2m,c2p)){ispts<<c1m.pts[0];};
  if(is_arc(c1m) && is_infatarc(c1m.pts[2],c2m,c2p)){ispts<<c1m.pts[2];};

  if(is_arc(c2p) && is_infatarc(c2p.pts[0],c1m,c1p)){ispts<<c2p.pts[0];};
  if(is_arc(c2p) && is_infatarc(c2p.pts[2],c1m,c1p)){ispts<<c2p.pts[2];};
  if(is_arc(c2m) && is_infatarc(c2m.pts[0],c1m,c1p)){ispts<<c2m.pts[0];};
  if(is_arc(c2m) && is_infatarc(c2m.pts[2],c1m,c1p)){ispts<<c2m.pts[2];};
 
  if(is_arc(c1p) && is_infatarc(c1p.critpt(0),c2m,c2p)){ispts<<c1p.critpt(0);}; 
  if(is_arc(c1p) && is_infatarc(c1p.critpt(1),c2m,c2p)){ispts<<c1p.critpt(1);}; 
  if(is_arc(c1m) && is_infatarc(c1m.critpt(0),c2m,c2p)){ispts<<c1m.critpt(0);}; 
  if(is_arc(c1m) && is_infatarc(c1m.critpt(1),c2m,c2p)){ispts<<c1m.critpt(1);}; 

  if(is_arc(c2p) && is_infatarc(c2p.critpt(0),c1m,c1p)){ispts<<c2p.critpt(0);}; 
  if(is_arc(c2p) && is_infatarc(c2p.critpt(1),c1m,c1p)){ispts<<c2p.critpt(1);}; 
  if(is_arc(c2m) && is_infatarc(c2m.critpt(0),c1m,c1p)){ispts<<c2m.critpt(0);}; 
  if(is_arc(c2m) && is_infatarc(c2m.critpt(1),c1m,c1p)){ispts<<c2m.critpt(1);};

  
  Seq< std::vector<C> > cpts; cpts<<vec(C(0),C(0));cpts<<vec(C(0),C(1));cpts<<vec(C(1),C(0));cpts<<vec(C(1),C(1));

  for(unsigned i=0; i<cpts.size(); i++){ if( is_infatarc(cpts[i],c1m,c1p) && 
                          is_infatarc(cpts[i],c2m,c2p) ){ ispts<<cpts[i]; };
    }; 
 
  return ispts;
};

template<class C> 
domain<C> minmax_box(Seq<std::vector<C> > ispts, domain<C> box){
  domain<C> newbox=box;
  //  C minx, maxx, miny, maxy; 
  if(ispts.size()==0) {
    newbox=domain<C>(int(0));
  } else {
  int k=0;
  C minx=ispts[0][0], maxx=ispts[0][0], miny=ispts[0][1], maxy=ispts[0][1]; 
  for(unsigned i=0; i<ispts.size(); i++){
      if(k==0){ minx=ispts[i][0]; maxx=ispts[i][0]; miny=ispts[i][1]; maxy=ispts[i][1]; k++;};
      if(ispts[i][0]>maxx){maxx=ispts[i][0];};
      if(ispts[i][0]<minx){minx=ispts[i][0];};
      if(ispts[i][1]>maxy){maxy=ispts[i][1];};
      if(ispts[i][1]<miny){miny=ispts[i][1];};
  }; 

  Interval<C>  Jx( minx, maxx), Jy( miny, maxy) ;
   if((maxx-minx)>C(0) && (maxy-miny)>C(0))
     { newbox= domain<C>( box.I[0].m + box.delta()[0]*Jx ,  box.I[1].m + box.delta()[1]*Jy ); };
  }
  return(newbox);
};/*minmax box for the arc  intersection points and arc end points*/


//end of arc_rep fcts