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
00045
00046
00047
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
00167
00168
00169 split=0; arc=0;
00170
00171
00172
00173 int chk=chkbox2d(b,prec);
00174
00175
00176
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
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
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
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
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
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
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
00396
00397
00398 int d1=p.nbrow();
00399 int d2=p.nbcol();
00400
00401 Seq<box2d> sol;
00402
00403 box2d box1, box2, box3, box4;
00404 MT bt1(p), bt2(p), bt3(p), bt4(p);
00405
00406
00407 shrink(bt1,box1,rg1,rg2); cout <<bt1<<endl;
00408 subdivision(sol.rep(),bt1,box1,1,1,split);
00409
00410
00411
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
00418
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
00425
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