synaps/msolve/Bezier2d.h

00001 #ifndef SYNAPS_SOLVE_BEZIER2D_H
00002 #define SYNAPS_SOLVE_BEZIER2D_H
00003 
00004 
00005 #include <synaps/init.h>
00006 #include <synaps/base/Seq.h>
00007 #include <synaps/linalg.h>
00008 #include <synaps/mpoly.h>
00009 #include <synaps/solve/Bezier1d.h>
00010 
00011 
00012 __BEGIN_NAMESPACE_SYNAPS
00013 
00014 
00016 typedef  MatrDse<array::rep2d<long double> > MAT;
00017 
00018 
00019 
00021 struct box2d
00022 {
00023 
00024   long double x0, x1, y0, y1;
00025   typedef box2d self_t;
00026 
00027   box2d (): x0(0),x1(1),y0(0),y1(1) {}
00028   box2d(const box2d & b1): x0(b1.x0),x1(b1.x1),y0(b1.y0),y1(b1.y1) {}
00029 
00030   self_t operator=(const box2d & b1)
00031     {
00032       x0=b1.x0; x1=b1.x1; y0=b1.y0; y1=b1.y1;
00033       return *this;
00034     }
00035 };
00036 
00037 std::ostream& operator<<(std::ostream & os, const box2d & b)
00038 {
00039   os<<"[["<<b.x0<<","<<b.x1<<"],["<<b.y0<<","<<b.y1<<"]]";
00040   return os;
00041 }
00042 
00043 
00044 /* Functions used to change the appropriate co-ordinates of the box for 
00045 different splits of the matrix, depending on whether it is for the first or 
00046 the second newly generated matrix, that is has to change the coordinates for.
00047 PARAMETERS: @b@ :box whose coordinates have to be changed. 
00048 */
00049 void rsplit1box2d(box2d &b)
00050 {
00051 b.x1=(b.x0+b.x1)/2;
00052 }
00053 
00054 
00055 void rsplit2box2d(box2d &b)
00056 {
00057 b.x0=(b.x0+b.x1)/2;
00058 }
00059 
00060 void csplit1box2d(box2d &b)
00061 {
00062 b.y1=(b.y0+b.y1)/2;
00063 }
00064 
00065 void csplit2box2d(box2d &b)
00066 {
00067 b.y0=(b.y0+b.y1)/2;
00068 }
00069 
00070 
00076 int chkbox2d(box2d b,double e)
00077 {
00078   int ret=0;
00079   if((b.x1-b.x0)<e) ret+=1;
00080   if((b.y1-b.y0)<e) ret+=2;
00081   return(ret);
00082 }
00083 
00087 template<class T>
00088 long double isolate(const T &l, long double c0, long double c1,
00089                     const double eps)
00090 {
00091   if((c1-c0)<eps)
00092     {return (c1+c0)/2;}
00093   else if (*l.begin()==0) return c0;
00094   else 
00095     {
00096       array::rep1d<double> l0(l.size()), l1(l);
00097       subdivise(l0,l1);
00098       if(Variation(l0)==1) 
00099         return isolate(l0,c0,(c0+c1)/2,eps);
00100       else if(Variation(l1)==1)
00101         return isolate(l1,(c0+c1)/2,c1,eps);
00102       else if(*l1.begin()==0) return (c0+c1)/2;
00103       else 
00104         {cout<<"Pb in the isolation"<<endl; return (c0+c1)/2;}
00105   }
00106 }
00107 
00108 
00114 struct SPL
00115 {
00116   long double prec;
00117   SPL(): prec(0.0001)  {}
00118   SPL(long double eps): prec(eps) {}
00119   SPL(const SPL & spl):prec(spl.prec) {}
00120       
00121   template<class T> int operator()(const T & ma, const box2d & b) const;
00122 };
00123 
00124 
00144 template<class T>
00145 int SPL::operator()(const T & ma,const box2d & b) const  
00146 {
00147  
00148  int split=0, arc=0;
00149  int nr=ma.nbrow()-1;
00150  int nc=ma.nbcol()-1;
00151  int rowvar0,rowvar1,colvar1,colvar0;
00152 
00153 
00154  MAT::row_type rt;
00155  MAT::col_type ct;
00156 
00157  MAT::row_type r0=Row(ma,0);
00158  MAT::row_type r1=Row(ma,nr);  
00159  MAT::col_type c0=Col(ma,0);
00160  MAT::col_type c1=Col(ma,nc);
00161 
00162  rowvar0=Variation(r0);
00163  rowvar1=Variation(r1);
00164  colvar0=Variation(c0);
00165  colvar1=Variation(c1);
00166  // cout <<"-> "<<rowvar0<<" "<<rowvar1<<" "<<colvar0<<" "<<colvar1<<" : ";
00167 
00168  //Spliting information; Arc type.
00169  split=0; arc=0;
00170 
00171 
00172  //Shape size information.
00173  int chk=chkbox2d(b,prec); 
00174  //cout <<chk<<" ";
00175 
00176  // small boxin both direction
00177  if(chk==3){
00178    if(rowvar0) arc+=1; 
00179    if(rowvar1) arc+=2;
00180    if(colvar0) arc+=4;
00181    if(colvar1) arc+=8; 
00182    if(arc==13 || arc==14)  
00183      return 12;
00184    else if (arc==7 || arc==11) 
00185      return 3;
00186    else if(arc)
00187      return 15;
00188    else 
00189      return 0;
00190  } else if (chk==1) {
00191    if(rowvar0) arc+=1; 
00192    if(rowvar1) arc+=2;
00193    if(colvar0==1) arc+=4; else if(colvar0>1) return -2;
00194    if(colvar1==1) arc+=8; else if(colvar1>1) return -2;
00195    if(arc==7 || arc==11 || arc==13 || arc==14 || arc==15) return -2;
00196    if(arc) return arc;
00197  } else if(chk==2) {
00198    if(rowvar0==1) arc+=1; else if(rowvar0>1) return -1;
00199    if(rowvar1==1) arc+=2; else if(rowvar1>1) return -1;
00200    if(colvar0) arc+=4;
00201    if(colvar1) arc+=8;
00202    if(arc==7 || arc==11 || arc==13 || arc==14 || arc==15) return -1;
00203    if(arc) return arc;
00204  } else {
00205    for(int i=0;i<=nr;i++)
00206      {
00207        rt=Row(ma,i);
00208        if (Variation(rt)>1) return -1;
00209      }
00210    for(int i=0;i<=nc;i++)
00211      {
00212        ct=Col(ma,i);
00213        if (Variation(ct)>1) return -2;
00214      }
00215    if(rowvar0==1) arc+=1; else if(rowvar0>1) return -1;
00216    if(rowvar1==1) arc+=2; else if(rowvar1>1) return -1;
00217    if(colvar0==1) arc+=4; else if(colvar0>1) return -2;
00218    if(colvar1==1) arc+=8; else if(colvar1>1) return -2;
00219 
00220    if(arc==7)  return -1;
00221    else if(arc==11) return -1;
00222    else if(arc==13) return -2;
00223    else if(arc==14) return -2;
00224    else if(arc==15) return -2;   
00225    else if(arc==0)  return 0;
00226    return arc;
00227  }
00228  return 0;
00229 }
00230 
00231  
00232 
00233 
00241 template<class T,class S>
00242 void subdivision(vector<box2d>& sol,
00243                   const T& b, const box2d & box,
00244                   int s1, int s2, 
00245                   const S & spl)
00246 {
00247   box2d box0(box), box1(box);
00248   
00249   int d1= b.nbrow();
00250   int d2= b.nbcol();
00251  
00252   MAT::row_type r0=Row(b,0), r1=Row(b,(d1-1));  
00253   MAT::col_type c0=Col(b,0), c1=Col(b,(d2-1));  
00254 
00255   int oui=spl(b,box);
00256   //  cout<<oui<<"\n"<<endl;
00257 
00258   if(oui==-1)
00259     {
00260       MAT b0(d1,d2);
00261       MAT b1(b.rep());
00262       
00263       for(int i=0;i<d1;i++)
00264        {
00265          MAT::row_type r0=Row(b0,i), r1=Row(b1,i);
00266          subdivise(r0,r1);
00267        }
00268      
00269       rsplit1box2d(box0);
00270       rsplit2box2d(box1);
00271       subdivision(sol,b0,box0,s1,s2,spl);
00272       subdivision(sol,b1,box1,s1,s2,spl);
00273     }
00274   else if(oui==-2)
00275     {
00276       MAT b0(d1,d2);
00277       MAT b1(b.rep());
00278       for(int i=0;i<d2;i++)
00279         {
00280           MAT::col_type c0=Col(b0,i), c1=Col(b1,i);
00281           subdivise(c0,c1);
00282         }
00283      csplit1box2d(box0);
00284      csplit2box2d(box1);
00285      
00286      subdivision(sol,b0,box0,s1,s2,spl);
00287      subdivision(sol,b1,box1,s1,s2,spl); 
00288     }  
00289   else if(oui>0) {
00290     
00291     box2d sbox(box);
00292     
00293     int c=0;
00294     // Arc through S
00295     if(oui%2)
00296       {
00297         array::rep1d<double> l0(d2);
00298         for(int i=0;i<d2;++i) {l0[i]=r0[i];}
00299         sbox.x0= isolate(l0,box.x0,box.x1,spl.prec);
00300         sbox.y0= box.y0;
00301         c++;
00302       }
00303     oui/=2;
00304     // Arc through N
00305     if(oui%2)
00306       {
00307         array::rep1d<double> l0(d2);
00308         for(int i=0;i<d2;++i) {l0[i]=r1[i];}
00309         if(c){
00310           sbox.x1= isolate(l0,box.x0,box.x1,spl.prec);
00311           sbox.y1= box.y1;
00312         } else {
00313           sbox.x0= isolate(l0,box.x0,box.x1,spl.prec);
00314           sbox.y0= box.y1;
00315         }
00316         c++;
00317       }
00318     oui/=2;
00319     // Arc through E
00320     if(oui%2)
00321       {
00322         array::rep1d<double> l0(d1);
00323         for(int i=0;i<d1;++i) {l0[i]=c0[i];}
00324         if(c){
00325           sbox.x1= box.x0;
00326           sbox.y1= isolate(l0,box.y0,box.y1,spl.prec);
00327         } else {
00328           sbox.x0= box.x0;
00329           sbox.y0= isolate(l0,box.y0,box.y1,spl.prec);
00330         }
00331         c++;
00332       }
00333     oui /=2;
00334     // Arc through 0
00335     if(oui%2)
00336       {
00337         array::rep1d<double> l0(d1);
00338         for(int i=0;i<d1;++i) {l0[i]=c1[i];}
00339         if(c){
00340           sbox.x1= box.x1;
00341           sbox.y1= isolate(l0,box.y0,box.y1,spl.prec);
00342         } else {
00343           sbox.x0= box.x1;
00344           sbox.y0= isolate(l0,box.y0,box.y1,spl.prec);
00345         }
00346         c++;
00347       }
00348 
00349     if(c>1) {
00350       sbox.x0=sbox.x0/(1-sbox.x0);
00351       sbox.x1=sbox.x1/(1-sbox.x1);
00352       sbox.y0=sbox.y0/(1-sbox.y0);
00353       sbox.y1=sbox.y1/(1-sbox.y1);
00354       sbox.x0*=s1; sbox.x1*=s1;
00355       sbox.y0*=s2; sbox.y1*=s2;
00356       //      cout <<".... "<<sbox<<endl;
00357       sol.push_back(sbox);
00358     }
00359   }
00360 }
00361  
00362 
00363 //----------------------------------------------------------------------
00364 template<class MT>
00365 void shrink(MT & ma, box2d & b, long double rg1, long double rg2)
00366 {
00367   for(unsigned i=0;i<ma.nbrow();i++)
00368     for(unsigned j=0; j<ma.nbcol()-1;j++)
00369       for(int k= ma.nbcol()-1;k>(int)j;k--)
00370         {
00371           ma(i,k)= rg1*ma(i,k)+(1-rg1)*ma(i,k-1);
00372         }
00373   for(unsigned i=0;i<ma.nbcol();i++)
00374     for(unsigned j=0; j<ma.nbrow()-1;j++)
00375       for(int k= ma.nbrow()-1;k>(int)j;k--)
00376         {
00377           ma(k,i)= rg1*ma(k,i)+(1-rg1)*ma(k-1,i);
00378         }
00379   b.x1=rg1; b.y1=rg2;
00380 }
00381 
00382 //----------------------------------------------------------------------
00389 template<class MT, class SPLIT>
00390 Seq<box2d> Draw(const MT & p,
00391                 const SPLIT& split,
00392                 long double rg1, long double rg2)
00393 {
00394   
00395   //cout <<"matrix\n"<<b<<endl;
00396   //initialised with [0,1]x[0,1]
00397 
00398   int d1=p.nbrow();
00399   int d2=p.nbcol();
00400 
00401   Seq<box2d> sol;//= new vector<box2d>;
00402 
00403   box2d box1, box2, box3, box4;
00404   MT bt1(p), bt2(p), bt3(p), bt4(p);
00405 
00406   //cout <<"                                 [0,+oo[x[0,+oo["<<endl;
00407   shrink(bt1,box1,rg1,rg2);  cout <<bt1<<endl;
00408   subdivision(sol.rep(),bt1,box1,1,1,split);
00409 
00410   //changing x to -x
00411   //cout <<"                                 ]-oo, 0]x[0,+oo["<<endl;
00412   for (int i=0;i<d1;i++)
00413     for (int j=1;j<d2;j+=2) {bt3(i,j)*=-1;} 
00414   shrink(bt3,box3,rg1,rg2);  cout <<bt2<<endl;
00415   subdivision(sol.rep(),bt3,box3,-1,1,split);
00416   
00417   //changing y to -y 
00418   //cout <<"                                 [0,+oo[x]-oo,0]"<<endl;
00419   for (int i=1;i<d1;i+=2)
00420     for (int j=0;j<d2;j++) {bt2(i,j)*=-1;} 
00421   shrink(bt2,box2,rg1,rg2);  cout <<bt3<<endl;
00422   subdivision(sol.rep(),bt2,box2,1,-1,split); 
00423   
00424   //changing x to -x , y to -y
00425   //cout <<"                                 ]-oo,0]x]-oo,0]"<<endl;
00426    for (int i=0;i<d1;i++)
00427     for (int j=1;j<d2;j+=2) {bt4(i,j)*=-1;} 
00428    for (int i=1;i<d1;i+=2)
00429      for (int j=0;j<d2;j++) {bt4(i,j)*=-1;} 
00430   shrink(bt4,box4,rg1,rg2);  cout <<bt4<<endl;
00431   subdivision(sol.rep(),bt4,box4,-1,-1,split);
00432 
00433  return sol;
00434 }
00435 
00436 //----------------------------------------------------------------------
00437 template<class R,class O>
00438 const MAT build_matrix(const MPol<R,O> & p)
00439 {
00440   int d0= Degree(p,0)+1, d1=Degree(p,1)+1;
00441   MAT  A(d1,d0);
00442   MATRIX::init(A,0);
00443   for(typename MPol<R,O>::const_iterator it=p.begin();it!=p.end();it++)
00444     {
00445       A((*it)[1],(*it)[0])= it->coeff();
00446     }
00447   for (int k=0;k<d0;k++)    
00448    for (int j=0;j<d1;j++)
00449    { 
00450     A(k,j)/= (binomial(k,(d0-1))*binomial(j,(d1-1)));  
00451    }
00452 
00453   return A;
00454 }
00455 //----------------------------------------------------------------------
00461 template<class R, class O, class SPLIT>
00462 Seq<box2d> Draw(const MPol<R,O> & p,
00463                 const SPLIT& spl,
00464                 long double rg1=0.93, long double rg2=0.93)
00465 {
00466   MAT A = build_matrix(p);
00467   return Draw(A,spl,rg1,rg2);
00468 }
00469 
00470 
00471 __END_NAMESPACE_SYNAPS
00472 
00473 
00474 //----------------------------------------------------------------------
00475 #endif // SYNAPS_SOLVE_BEZIER2D_H
00476 

SYNAPS DOCUMENTATION
logo