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