shape: include/shape/ssi/ssiqtsl.hpp Source File

00001 # ifndef SHAPE_SSIQTSL_HPP
00002 # define SHAPE_SSIQTSL_HPP
00003 # include <cmath>
00004 # include <vector>
00005 # include <fstream>
00006 # include <shape/ssi/ssi_def.hpp>
00007 # include <shape/ssi/ssiqts.hpp>
00008 # include <shape/ssi/ssiqts_tri_tri.hpp>
00009 
00010 # define SSIQTS ssiqts<C,V>
00011 # define SSIQTSL ssiqtsl<C,V>
00012 # define TMPL template<class C, class V>
00013 
00014 //# define ParametricSurface shape::parametric_surface<double>
00015 
00016 namespace mmx {
00017 
00018 template<class C, class V= shape::default_env>
00019 struct ssiqtsl : SSIQTS
00020 {
00021     typedef typename SSIQTS::vector3 vector3;
00022     typedef typename shape::use<shape::ssi_def,C,V>::ParametricSurface  ParametricSurface;
00023     typedef fxv<double,2> vector2;
00024 
00025     std::vector< vector2 > lpts;
00026     std::vector< vector2 > rpts;
00027     std::vector< vector3 > spcs;
00028     std::vector< double  > errs;
00029     std::vector< vector3 > spcsorig;
00030     double errmax;
00031     void gmvdump();
00032     void gpdump();
00033 
00034     ssiqtsl (ParametricSurface * srfa, ParametricSurface * srfb, int n );
00035 
00036     static void space2prm(vector2 & pa, vector2 & pb,  const SSIQTSL::vector3& sa,
00037                           const SSIQTSL::vector3& sb,
00038                           const SSIQTSL::vector3& base,
00039                           const SSIQTSL::vector3& pu,
00040                           const SSIQTSL::vector3& pv );
00041     static void merge( typename SSIQTS::aabb3 & box, const typename SSIQTS::aabb3& a, const typename SSIQTS::aabb3 & b );
00042     static void solve_cf( vector3 ** qpa, vector3 ** qpb, int a, int b, SSIQTSL & ssi  );
00043     static void sample_curve3d( ParametricSurface * srfa, ParametricSurface * srfb, SSIQTSL& ssi );
00044     static void solve( SSIQTSL & ssi );
00045 };
00046 
00047 
00048 void POLYLINE( std::ostream& o, double * src, int nv, float * colors, int nc )
00049 {
00050     o << "{ = VECT\n";
00051     o << 1 << " " << nv << " " << nc << std::endl;
00052     o << nv << std::endl;
00053     o << nc << std::endl;
00054     for ( int i = 0; i < nv; i++, src += 3 )
00055         o << src[0] << " " << src[1] << " " << src[2] << std::endl;
00056     for ( int i = 0; i < nc; i++ )
00057         o << colors[i*4] << " " << colors[i*4+1] << " " << colors[i*4+2] << " " << colors[i*4+3] << std::endl;
00058     o << "}\n";
00059 };
00060 
00061 void UVSMP( std::ostream& o, double * src, unsigned m, unsigned n )
00062 {
00063 #define _ij_(i,j) (i)*n+j
00064     o << " { = MESH\n";
00065     o << m << " " << n << std::endl;
00066     for ( unsigned i = 0; i < m*n; i ++ )
00067         o << (src+i*3)[0] << " " <<  (src+i*3)[1] << " " << (src+i*3)[2] << std::endl;
00068     o << " }\n";
00069 };
00070 
00071 #define __up__   0
00072 #define __down__ 1
00073 /* return the up triangle of quad q */
00074 #define __up__triangle__(q)         q[0],q[3],q[2]
00075 /* return the down triangle of quad q */
00076 #define __down__triangle__(q)       q[0],q[1],q[2]
00077 /* return the parametric plane of up triangle of q */
00078 #define __up__param_triangle__(q)   q[3],q[2],q[0]
00079 /* return the parametric plane of down triangle of q */
00080 #define __down__param_triangle__(q) q[1],q[0],q[2]
00081 #define __convert_order__(i,point) { point[(i+1)%2] = 1.0-point[(i+1)%2]; }
00082 
00083 /* check for the intersection beetween {UP,DOWN}/{UP,DOWN} triangles */
00084 #define __ssiqtsl__triangle_triangle_case__(trig0,trig1)                            \
00085 {                                                                          \
00086     intersection =                                                          \
00087     intersectp_triangles3_isegment                                     \
00088     ( coplanar, seg[0], seg[1],                                         \
00089     trig0##triangle__(qa),                                            \
00090     trig1##triangle__(qb), double(1e-9) );                            \
00091     if ( intersection )                                                \
00092 {                                                                      \
00093     if ( !coplanar )                                                     \
00094 {                                                                  \
00095     SSIQTSL::space2prm(/* points in parameter space */                        \
00096     seg0[0], seg0[1],                                       \
00097     /* corresponding points in R3 */                        \
00098     seg[0],  seg[1],                                        \
00099     /* parametric plane = base + 2 vectors  */              \
00100     trig0##param_triangle__(qa)                             \
00101     );                                                      \
00102     SSIQTSL::space2prm( seg1[0], seg1[1],                                     \
00103     seg[0], seg[1],                                       \
00104     trig1##param_triangle__(qb) );                         \
00105     /* retrieve the correct offset coordinates in quad */            \
00106     __convert_order__(trig0,seg0[0]);                                \
00107     __convert_order__(trig0,seg0[1]);                                \
00108     __convert_order__(trig1,seg1[0]);                                \
00109     __convert_order__(trig1,seg1[1]);                                \
00110     \
00111     double * left_uvals = ssi.smpa->m_uvals;                         \
00112     double * left_vvals = ssi.smpa->m_vvals;                         \
00113     double * righ_uvals = ssi.smpb->m_uvals;                         \
00114     double * righ_vvals = ssi.smpb->m_vvals;                         \
00115     for ( int i = 0; i < 2; i ++ )                                   \
00116 {                                                              \
00117     seg0[i][0] = left_uvals[a/ssi.m] + seg0[i][0] * (left_uvals[a/ssi.m+1]-left_uvals[a/ssi.m]);\
00118     seg0[i][1] = left_vvals[a%ssi.m] + seg0[i][1] * (left_vvals[a%ssi.m+1]-left_vvals[a%ssi.m]);\
00119 };                                                             \
00120     for ( int i = 0; i < 2; i ++ )                                   \
00121 {                                                              \
00122     seg1[i][0] = righ_uvals[b/ssi.m] + seg1[i][0] * (righ_uvals[b/ssi.m+1]-righ_uvals[b/ssi.m]);\
00123     seg1[i][1] = righ_vvals[b%ssi.m] + seg1[i][1] * (righ_vvals[b%ssi.m+1]-righ_vvals[b%ssi.m]);\
00124 };                                                             \
00125     ssi.lpts.push_back(seg0[0]);                                     \
00126     ssi.lpts.push_back(seg0[1]);                                     \
00127     ssi.rpts.push_back(seg1[0]);                                     \
00128     ssi.rpts.push_back(seg1[1]);                                     \
00129 }                                                                  \
00130 }                                                                      \
00131 }
00132 
00133 TMPL
00134 //static
00135 void SSIQTSL::space2prm( SSIQTSL::vector2 & pa,
00136                          SSIQTSL::vector2 & pb,
00137                          const SSIQTSL::vector3& sa,
00138                          const SSIQTSL::vector3& sb,
00139                          const SSIQTSL::vector3& base,
00140                          const SSIQTSL::vector3& pu,
00141                          const SSIQTSL::vector3& pv )
00142 {
00143     /* T(u,v) = base + u*bu +v*bv
00144      =>
00145      spc[0] - base[0] = delta[0]  = / bu[0] bv[0]\  / u \
00146      spc[1] - base[1] = delta[1]  = | bu[1] bv[1]|  |   |
00147      spc[2] - base[2] = delta[2]  = \ bu[2] bv[2]/  \ v /
00148   */
00149     typename SSIQTSL::vector3 bu;
00150     for ( int i = 0; i < 3; i ++ ) bu[i] = pu[i]-base[i];
00151     typename SSIQTSL::vector3 bv;
00152     for ( int i = 0; i < 3; i ++ ) bv[i] = pv[i]-base[i];
00153     double muu, mvv, muv;
00154     muu = 0;
00155     for ( int i = 0; i < 3; i ++ ) muu += bu[i]*bu[i];
00156     mvv = 0;
00157     for ( int i = 0; i < 3; i ++ ) mvv += bv[i]*bv[i];
00158     muv = 0;
00159     for ( int i = 0; i < 3; i ++ ) muv += bu[i]*bv[i];
00160     double detm = muu*mvv - muv*muv;
00161     typename SSIQTSL::vector3 delta;
00162     double x, y;
00163     for ( int k = 0; k < 3; k ++ ) delta[k] = sa[k]-base[k];
00164     x = 0;
00165     for ( int k = 0; k < 3; k ++ ) x += bu[k]*delta[k];
00166     y = 0;
00167     for ( int k = 0; k < 3; k ++ ) y += bv[k]*delta[k];
00168     pa[0] = (mvv * x - muv * y)/detm;
00169     pa[1] = (muu * y - muv * x)/detm;
00170     for ( int k = 0; k < 3; k ++ ) delta[k] = sb[k]-base[k];
00171     x = 0;
00172     for ( int k = 0; k < 3; k ++ ) x += bu[k]*delta[k];
00173     y = 0;
00174     for ( int k = 0; k < 3; k ++ ) y += bv[k]*delta[k];
00175     pb[0] = (mvv * x - muv * y)/detm;
00176     pb[1] = (muu * y - muv * x)/detm;
00177 };
00178 
00179 TMPL
00180 //static
00181 void SSIQTSL::merge( typename SSIQTS::aabb3 & box, const typename SSIQTS::aabb3& a, const typename SSIQTS::aabb3 & b )
00182 {
00183     for ( unsigned i = 0; i < 3; i ++ )
00184     {
00185         box[i][0] = std::min(a[i][0],b[i][0]);
00186         box[i][1] = std::max(a[i][1],b[i][1]);
00187     };
00188 };
00189 
00190 TMPL
00191 //static
00192 void SSIQTSL::solve_cf( typename SSIQTSL::vector3 ** qpa, typename SSIQTSL::vector3 ** qpb, int a, int b, SSIQTSL & ssi  )
00193 {
00194     // scaling
00195     typename SSIQTSL::aabb3 zoom;
00196     double mx = 0;
00197     merge(zoom,ssi.boxa[a],ssi.boxb[b]);
00198     /*
00199   std::cout << "boxa\n";
00200   for ( int i = 0; i < 3; i ++ )
00201     std::cout << boxa[a][i][0] << " " << boxa[a][i][1] << std::endl;
00202   std::cout << "boxb\n";
00203   for ( int i = 0; i < 3; i ++ )
00204     std::cout << boxb[b][i][0] << " " << boxb[b][i][1] << std::endl;
00205   */
00206     for ( int k = 0; k < 3; k ++ ) {
00207         zoom[k][1] -= zoom[k][0];
00208         if ( zoom[k][1] > mx ) mx = zoom[k][1];
00209     };
00210 
00211     double sc = double(1.0)/mx;
00212     typename SSIQTSL::vector3 qa[4];
00213     typename SSIQTSL::vector3 qb[4];
00214 
00215     for ( int i = 0; i < 4; i ++ )
00216         for ( int k = 0; k < 3; k ++ )
00217         {
00218             qa[i][k] = ((*(qpa[i]))[k] - zoom[k][0]) * sc;
00219             qb[i][k] = ((*(qpb[i]))[k] - zoom[k][0]) * sc;
00220         };
00221 
00222     typename SSIQTSL::vector3 seg [2];
00223     typename SSIQTSL::vector2 seg0[2];
00224     typename SSIQTSL::vector2 seg1[2];
00225     bool coplanar = false;
00226     bool intersection;
00227     __ssiqtsl__triangle_triangle_case__(__up__,__up__);
00228     __ssiqtsl__triangle_triangle_case__(__up__,__down__);
00229     __ssiqtsl__triangle_triangle_case__(__down__,__up__);
00230     __ssiqtsl__triangle_triangle_case__(__down__,__down__);
00231 };
00232 
00233 
00234 TMPL
00235 void SSIQTSL::solve( SSIQTSL & ssi )
00236 {
00237     typename SSIQTSL::vector3 * sleft  = (typename SSIQTSL::vector3*)(&(ssi.smpa->m_svals)[0]);
00238     typename SSIQTSL::vector3 * sright = (typename SSIQTSL::vector3*)(&(ssi.smpb->m_svals)[0]);
00239 
00240 
00241     for( int * cf = ssi.m_bcf; cf != ssi.m_ecf; cf += *cf )
00242     {
00243         if ( cf[cfbg_] != -1 )
00244         {
00245             int lu = cf[cflq_]/ssi.m;
00246             int lv = cf[cflq_]%ssi.m;
00247             typename SSIQTSL::vector3 * qpl[4];
00248             qpl[0] = sleft  + lu*(ssi.m+1)+lv;
00249             qpl[1] = sleft  + (lu+1)*(ssi.m+1)+lv;
00250             qpl[2] = qpl[1] + 1;
00251             qpl[3] = qpl[0] + 1;
00252             int * rpq = cf+cfbg_;
00253             while ( *rpq != -1 && rpq != cf+*cf )
00254             {
00255                 int ru = (*rpq/ssi.m);
00256                 int rv = (*rpq%ssi.m);
00257                 typename SSIQTSL::vector3 * qpr[4];
00258                 qpr[0] = sright + ru*(ssi.m+1)+rv;
00259                 qpr[1] = sright + (ru+1)*(ssi.m+1)+rv;
00260                 qpr[2] = qpr[1] + 1;
00261                 qpr[3] = qpr[0] + 1;
00262                 solve_cf(qpl,qpr,cf[cflq_],*rpq,ssi);
00263                 rpq++;
00264             };
00265         }
00266     };
00267 };
00268 
00269 
00270 TMPL
00271 void SSIQTSL::gpdump()
00272 {
00273     std::ofstream gpl("SSIQTSL-L.dat");
00274     std::ofstream gpr("SSIQTSL-R.dat");
00275 
00276     for ( unsigned i = 0; i < lpts.size(); i ++ )
00277         gpl << lpts[i][0] << " " << lpts[i][1] << std::endl;
00278     for ( unsigned i = 0; i < rpts.size(); i ++ )
00279         gpr << rpts[i][0] << " " << rpts[i][1] << std::endl;
00280 };
00281 
00282 TMPL
00283 void SSIQTSL::gmvdump()
00284 {
00285     std::ofstream o("SSIQTSL.gmv");
00286     o << "LIST\n";
00287     float colora[] = { 0.0, 0.0, 1.0 };
00288     float colorb[] = { 0.0, 1.0, 0.0 };
00289     o << "appearance {\n material{\n diffuse " <<
00290          colora[0] << " " <<  colora[1] << " " << colora[2] << "\n";
00291     o << "}\n}\n";
00292     UVSMP(o,&(this->smpa->m_svals[0][0]),this->m+1,this->m+1);
00293     o << "appearance {\n material{\n diffuse " <<
00294          colorb[0] << " " <<  colorb[1] << " " << colorb[2] << "\n";
00295     o << "}\n}\n";
00296     UVSMP(o,&(this->smpb->m_svals[0][0]),this->m+1,this->m+1);
00297     float color[] = { 1.0, 0.0, 0.0, 1.0 };
00298     for ( unsigned i = 0; i < spcs.size(); i += 2 )
00299     {
00300         color[1] = double(i % 2);
00301         color[2] = double((i+1)%2);
00302         POLYLINE(o,&(spcs[i][0]),2,color,1);
00303     };
00304 };
00305 
00306 TMPL
00307 void SSIQTSL::sample_curve3d( ParametricSurface * srfa,
00308                               ParametricSurface * srfb, SSIQTSL& ssi )
00309 {
00310     //  std::cout << " left = " << ssi.lpts.size() << std::endl;
00311     int nl = ssi.lpts.size();
00312     //  std::cout << " right = " << ssi.rpts.size() << std::endl;
00313     ssi.spcs.resize( nl );
00314     std::vector<typename SSIQTSL::vector3> srfb_samples(nl);
00315     shape::use<shape::ssi_def,C,V>::eval(srfa, (double*)(&(ssi.spcs[0])), (const double*)(&ssi.lpts[0]), ssi.lpts.size() );
00316     shape::use<shape::ssi_def,C,V>::eval(srfb, (double*)(&(srfb_samples[0])), (const double*)(&ssi.rpts[0]), ssi.rpts.size() );
00317     //  std::vector<double> errs;
00318     //  ssi.errs.resize( ssi.lpts.size() );
00319     ssi.errmax = 0;
00320     for ( unsigned i = 0; i < ssi.lpts.size(); i ++ )
00321     {
00322 
00323         //      int lk0, lk1, stat;
00324         //      SSIQTSL::vector3 pa,pb;
00325         typename SSIQTSL::vector3 delta;
00326         double s = 0;
00327         //      std::cout << ssi.spcs[i] << std::endl;
00328         //      std::cout << srfb_samples[i] << std::endl;
00329         for ( int k =0; k < 3; k ++ )
00330             delta[k] = ssi.spcs[i][k]-srfb_samples[i][k];
00331         s = max_abs( delta );
00332         for ( int k = 0; k < 3; k ++ )
00333             ssi.spcs[i][k] = (ssi.spcs[i][k]+srfb_samples[i][k])/2;
00334         if ( s > ssi.errmax ) ssi.errmax = s;
00335     };
00336 };
00337 
00338 TMPL
00339 SSIQTSL::ssiqtsl( ParametricSurface *  srfa,
00340                   ParametricSurface *  srfb,
00341                   int n ) : SSIQTS(srfa,srfb,n,1,1)
00342 {
00343     solve(*this);
00344     sample_curve3d( srfa, srfb, *this );
00345     //  std::cout << errmax << std::endl;
00346     // gmvdump();
00347 };
00348 
00349 
00350 // gmv::stream & operator<<( gmv::stream & o,
00351 //                        const ssi_linear & ssi )
00352 // {
00353 //   float colora[] = { 0.0, 0.0, 1.0 };
00354 //   float colorb[] = { 0.0, 1.0, 0.0 };
00355 //   o.matcolor(colora);
00356 //   o << (*ssi.smpa);
00357 //   o << (*ssi.smpb);
00358 //   float color[] = { 1.0, 0.0, 0.0, 1.0 };
00359 //   for ( int i = 0; i < ssi.spcs.size(); i += 2 )
00360 //     {
00361 //       //    color[1] = double(i % 2);
00362 //       //      color[2] = double((i+1)%2);
00363 //       o.polyline( &(ssi.spcs[i][0]), 2, color, 1);
00364 //     };
00365 //   return o;
00366 // };
00367 
00368 } //namespace mmx
00369 
00370 # undef TMPL
00371 # undef SSIQTS
00372 # undef SSIQTSL
00373 # endif
00374 
00375