00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include "Surface.h"
00021
00022 #include <fstream>
00023 #include <set>
00024 #include <math.h>
00025
00026 #define PI 3.14159265
00027
00028 #include <GTL/components.h>
00029
00030 #include <TarisApplication.h>
00031
00032 Surface::Surface()
00033 {
00034 }
00035
00036 Surface::Surface( const vector<Point>& lpo, const vector<Triangle>& ltri ):
00037 pair< vector<Point>, vector<Triangle> >( lpo, ltri )
00038 {
00039 }
00040
00041 Surface::Surface( const Surface& surf ):
00042 pair< vector<Point>, vector<Triangle> >( surf.first, surf.second )
00043 {
00044 }
00045
00046 Surface::~Surface()
00047 {
00048 }
00049
00050 PointVector Surface::getPoints() const
00051 {
00052 return this->first ;
00053 }
00054
00055 TriangleVector Surface::getTriangles() const
00056 {
00057 return this->second ;
00058 }
00059
00060 void Surface::setPoints( const PointVector& pv )
00061 {
00062 this->first = pv ;
00063 }
00064
00065 void Surface::setTriangles( const TriangleVector& tv )
00066 {
00067 this->second = tv ;
00068 }
00069
00070 inline double distance2( double x1, double y1, double z1, double x2, double y2, double z2 )
00071 {
00072 if( ( (x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2) ) <= 0 )
00073 cerr << "Error: Distancia euclidiana menor a cero !!!" << endl ;
00074 return (x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2) ;
00075 }
00076
00077
00078 double Surface::computeArea(){
00079
00080 double area = 0.0 ;
00081
00082 for( TriangleVector::iterator it = this->second.begin(); it != this->second.end(); it++ ){
00083
00084
00085 double a = sqrt( distance2( it->getPoint1().getX(), it->getPoint1().getY(), it->getPoint1().getZ(),
00086 it->getPoint2().getX(), it->getPoint2().getY(), it->getPoint2().getZ() ) ) ;
00087
00088 double b = sqrt( distance2( it->getPoint1().getX(), it->getPoint1().getY(), it->getPoint1().getZ(),
00089 it->getPoint3().getX(), it->getPoint3().getY(), it->getPoint3().getZ() ) ) ;
00090
00091 double c = sqrt( distance2( it->getPoint2().getX(), it->getPoint2().getY(), it->getPoint2().getZ(),
00092 it->getPoint3().getX(), it->getPoint3().getY(), it->getPoint3().getZ() ) );
00093
00094 area += 0.25*sqrt( (a + (b + c))*(c - (a - b))*(c + (a - b))*(a + (b - c)) ) ;
00095 }
00096
00097 return area ;
00098 }
00099
00105 void Surface::printOogl( ostream& os )
00106 {
00107 if( TarisApplication::Instance()->getDebugLevel() >= 1 )
00108 cout << "Saving file in oogl format ... " ;
00109
00110 os << "OFF" << endl ;
00111 os << first.size() << " " << second.size() << endl ;
00112
00113 map<int, int> ids ;
00114 int i=0 ;
00115 for( PointVector::iterator it = first.begin(); it != first.end(); it++, i++ ){
00116 os << it->getX() << " " << it->getY() << " " << it->getZ() << endl ;
00117 ids[it->getId()] = i ;
00118 }
00119
00120 for( TriangleVector::iterator it = second.begin(); it != second.end(); it++ ){
00121 os << "3 " ;
00122 os << ids[ it->getPoint1().getId() ] << " " ;
00123 os << ids[ it->getPoint2().getId() ] << " " ;
00124 os << ids[ it->getPoint3().getId() ] << endl ;
00125 }
00126
00127 if( TarisApplication::Instance()->getDebugLevel() >= 1 )
00128 cerr << "OK" << endl ;
00129 }
00130
00135 SurfacesVector Surface::computeComponents()
00136 {
00137 SurfacesVector output ;
00138
00139 if( TarisApplication::Instance()->getDebugLevel() >= 2 ){
00140 cerr << endl ;
00141 cerr << "**************************************************************" << endl ;
00142 cerr << endl ;
00143 cerr << " ALGORITMO DE BUSQUEDA DE COMPONENTES CONEXAS " << endl ;
00144 cerr << " -------------------------------------------- " << endl ;
00145 cerr << endl ;
00146 cerr << endl ;
00147 cerr << " Construyendo grafo de superficie ... " ;
00148 }
00149
00150 graph g ;
00151 map<int, node> nodeId ;
00152
00153
00154
00155
00156 int i=0 ;
00157 for( PointVector::iterator it = first.begin(); it != first.end(); it++, i++ ){
00158 node n = g.new_node() ;
00159 nodeId[i] = n ;
00160 }
00161
00162
00163
00164
00165
00166 for( TriangleVector::iterator it = second.begin(); it != second.end(); it++ ){
00167 g.new_edge( nodeId[it->getPoint1().getId()], nodeId[it->getPoint2().getId()] ) ;
00168 g.new_edge( nodeId[it->getPoint2().getId()], nodeId[it->getPoint3().getId()] ) ;
00169 g.new_edge( nodeId[it->getPoint3().getId()], nodeId[it->getPoint1().getId()] ) ;
00170 }
00171
00172 if( TarisApplication::Instance()->getDebugLevel() >= 2 ){
00173 cerr << "OK" << endl ;
00174
00175 cerr << endl ;
00176 cerr << " Busqueda de componentes: " << endl ;
00177 cerr << " -------------------------" << endl ;
00178 cerr << endl ;
00179 }
00180
00181 components algol ;
00182
00183
00184
00185
00186
00187
00188 if( TarisApplication::Instance()->getDebugLevel() >= 2 )
00189 cerr << " Verificando condiciones del algoritmo ... " ;
00190
00191 if (algol.check(g) != algorithm::GTL_OK){
00192 if( TarisApplication::Instance()->getDebugLevel() >= 2 )
00193 cerr << "OK" << endl ;
00194 }else{
00195 if( TarisApplication::Instance()->getDebugLevel() >= 2 )
00196 cerr << "Fall贸!!!!!!" << endl ;
00197 }
00198
00199 if( TarisApplication::Instance()->getDebugLevel() >= 2 )
00200 cerr << " Ejecutando el algoritmo ... " ;
00201
00202 if (algol.run(g) == dfs::GTL_OK){
00203 if( TarisApplication::Instance()->getDebugLevel() >= 2 )
00204 cerr << "OK" << endl ;
00205 }else{
00206 if( TarisApplication::Instance()->getDebugLevel() >= 2 )
00207 cerr << "Fall贸!!!!!!" << endl ;
00208 }
00209
00210 if( TarisApplication::Instance()->getDebugLevel() >= 2 )
00211 cerr << " Reconstruyendo componentes ... " ;
00212
00213
00214
00215
00216
00217 set<Triangle> setTriangleTMP ;
00218
00219
00220 for( components::component_iterator it2 = algol.components_begin(); it2 != algol.components_end(); it2++ ){
00221
00222
00223
00224
00225
00226
00227
00229 if( it2->first.size() > 10 ){
00230 Surface outputSurfaceTMP ;
00231
00232 for( list<node>::iterator it3 = it2->first.begin(); it3 != it2->first.end(); it3++ ){
00233
00234
00235
00236
00237
00238
00239 outputSurfaceTMP.first.push_back( first[it3->id()] ) ;
00240
00241 for( TriangleVector::iterator it1 = second.begin(); it1 != second.end(); it1++ ){
00242
00243 if( it1->contains( it3->id() ) ){
00244 setTriangleTMP.insert( *it1 ) ;
00245 }
00246
00247 }
00248 }
00249
00250
00251
00252
00253
00254
00255 for( set<Triangle>::iterator it4 = setTriangleTMP.begin(); it4 != setTriangleTMP.end(); it4++ )
00256 outputSurfaceTMP.second.push_back( *it4 ) ;
00257
00258
00259
00260
00261 output.push_back( outputSurfaceTMP ) ;
00262
00263 setTriangleTMP.clear() ;
00264 }
00265
00266 }
00267
00268 if( TarisApplication::Instance()->getDebugLevel() >= 2 ){
00269 cerr << "OK" << endl ;
00270
00271 cerr << " N煤mero de componentes conexas encontradas = " << output.size() << endl ;
00272 cerr << endl ;
00273 cerr << " Breve descripci贸n: " << endl ;
00274 cerr << endl ;
00275 cerr << " COMPONENTE\tPUNTOS\tTRIANGULOS " << endl ;
00276 cerr << " -----------------------------------" << endl ;
00277 for( SurfacesVector::iterator it = output.begin(); it != output.end(); it++ ){
00278 cerr << " \t" << i << "\t\t" << it->first.size() << "\t" << it->second.size() << endl ;
00279 }
00280
00281 cerr << endl ;
00282 cerr << "**************************************************************" << endl ;
00283 cerr << endl ;
00284 }
00285
00286 return output ;
00287 }
00288
00293 bool Surface::contains( Surface s )
00294 {
00295 return s.isContainedIn( *this , 0, 0) ;
00296 }
00297
00315 bool Surface::isContainedIn( Surface s , int a, int b)
00316 {
00317 int i = 0 ;
00318 int insideVector[3] = {0, 0, 0} ;
00319 bool inside = false ;
00320
00321 for( TriangleVector::iterator it = s.second.begin(); it != s.second.end(); it++ ){
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332 double xPromInternal = ( this->second[a].getPoint1().getX() +
00333 this->second[a].getPoint2().getX() +
00334 this->second[a].getPoint3().getX() )/3.0 ;
00335
00336 double yPromInternal = ( this->second[a].getPoint1().getY() +
00337 this->second[a].getPoint2().getY() +
00338 this->second[a].getPoint3().getY() )/3.0 ;
00339
00340 double zPromInternal = ( this->second[a].getPoint1().getZ() +
00341 this->second[a].getPoint2().getZ() +
00342 this->second[a].getPoint3().getZ() )/3.0 ;
00343
00344
00345 double xPromExternal = ( s.second[b].getPoint1().getX() +
00346 s.second[b].getPoint2().getX() +
00347 s.second[b].getPoint3().getX() )/3.0 ;
00348
00349 double yPromExternal = ( s.second[b].getPoint1().getY() +
00350 s.second[b].getPoint2().getY() +
00351 s.second[b].getPoint3().getY() )/3.0 ;
00352
00353 double zPromExternal = ( s.second[b].getPoint1().getZ() +
00354 s.second[b].getPoint2().getZ() +
00355 s.second[b].getPoint3().getZ() )/3.0 ;
00356
00357
00358
00359
00360
00361
00362
00363
00364 int j = a+1 ;
00365 while( fabs(xPromExternal-xPromInternal) < 1e-6 ){
00366
00367 xPromInternal = ( this->second[j].getPoint1().getX() +
00368 this->second[j].getPoint2().getX() +
00369 this->second[j].getPoint3().getX() )/3.0 ;
00370
00371 yPromInternal = ( this->second[j].getPoint1().getY() +
00372 this->second[j].getPoint2().getY() +
00373 this->second[j].getPoint3().getY() )/3.0 ;
00374
00375 zPromInternal = ( this->second[j].getPoint1().getZ() +
00376 this->second[j].getPoint2().getZ() +
00377 this->second[j].getPoint3().getZ() )/3.0 ;
00378
00379 j++ ;
00380
00381 }
00382
00383
00384
00385
00386
00387
00388 double a = ( it->getPoint2().getY() - it->getPoint1().getY() )*( it->getPoint3().getZ() - it->getPoint1().getZ() )-
00389 ( it->getPoint2().getZ() - it->getPoint1().getZ() )*( it->getPoint3().getY() - it->getPoint1().getY() ) ;
00390
00391 double b = ( it->getPoint2().getZ() - it->getPoint1().getZ() )*( it->getPoint3().getX() - it->getPoint1().getX() )-
00392 ( it->getPoint2().getX() - it->getPoint1().getX() )*( it->getPoint3().getZ() - it->getPoint1().getZ() ) ;
00393
00394 double c = ( it->getPoint2().getX() - it->getPoint1().getX() )*( it->getPoint3().getY() - it->getPoint1().getY() )-
00395 ( it->getPoint2().getY() - it->getPoint1().getY() )*( it->getPoint3().getX() - it->getPoint1().getX() ) ;
00396
00397 double d = a*it->getPoint1().getX()+b*it->getPoint1().getY()+c*it->getPoint1().getZ() ;
00398
00399
00400
00401
00402
00403
00404
00405
00406 double x = (d-
00407 b*
00408 (
00409 (xPromInternal*yPromInternal-xPromInternal*yPromExternal)
00410 /
00411 ( xPromExternal-xPromInternal )
00412 +
00413 yPromInternal
00414 )
00415 -
00416 c*
00417 (
00418 (xPromInternal*zPromInternal-xPromInternal*zPromExternal)
00419 /
00420 ( xPromExternal-xPromInternal )
00421 +
00422 zPromInternal
00423 )
00424
00425 )
00426 /
00427 (
00428 a+b*( yPromExternal - yPromInternal )
00429 /
00430 ( xPromExternal - xPromInternal )
00431 +
00432 c*( zPromExternal - zPromInternal )
00433 /
00434 ( xPromExternal - xPromInternal )
00435 ) ;
00436
00437
00438
00439 double y = ( x - xPromInternal )*( yPromExternal - yPromInternal )/
00440 ( xPromExternal - xPromInternal ) + yPromInternal;
00441
00442 double z = ( x - xPromInternal )*( zPromExternal - zPromInternal )/
00443 ( xPromExternal - xPromInternal ) + zPromInternal;
00444
00445
00446
00447
00448 double dp1 = distance2( x, y, z, it->getPoint1().getX(), it->getPoint1().getY(), it->getPoint1().getZ() ) ;
00449 double dp2 = distance2( x, y, z, it->getPoint2().getX(), it->getPoint2().getY(), it->getPoint2().getZ() ) ;
00450 double dp3 = distance2( x, y, z, it->getPoint3().getX(), it->getPoint3().getY(), it->getPoint3().getZ() ) ;
00451
00452
00453 double d12 = distance2( it->getPoint1().getX(), it->getPoint1().getY(), it->getPoint1().getZ(),
00454 it->getPoint2().getX(), it->getPoint2().getY(), it->getPoint2().getZ() ) ;
00455
00456 double d13 = distance2( it->getPoint1().getX(), it->getPoint1().getY(), it->getPoint1().getZ(),
00457 it->getPoint3().getX(), it->getPoint3().getY(), it->getPoint3().getZ() ) ;
00458
00459 double d23 = distance2( it->getPoint2().getX(), it->getPoint2().getY(), it->getPoint2().getZ(),
00460 it->getPoint3().getX(), it->getPoint3().getY(), it->getPoint3().getZ() ) ;
00461
00462
00463 double alpha_12 = acos( ( dp1 + dp2 - d12 )/( 2.0*sqrt(dp1*dp2) ) )*180.0/PI ;
00464 double alpha_13 = acos( ( dp1 + dp3 - d13 )/( 2.0*sqrt(dp1*dp3) ) )*180.0/PI ;
00465 double alpha_23 = acos( ( dp2 + dp3 - d23 )/( 2.0*sqrt(dp2*dp3) ) )*180.0/PI ;
00466
00467
00468 double totalAlpha = alpha_12 + alpha_13 + alpha_23 ;
00469
00470
00471
00472
00473
00474
00475
00476 if( ( fabs( totalAlpha - 360.0 ) < 1e-6 ) ){
00477 if( !inside ){
00478 insideVector[0] = ( (xPromExternal-xPromInternal)>=0 ? 1:-1 ) ;
00479 insideVector[1] = ( (yPromExternal-yPromInternal)>=0 ? 1:-1 ) ;
00480 insideVector[2] = ( (zPromExternal-zPromInternal)>=0 ? 1:-1 ) ;
00481 inside = true ;
00482 i++ ;
00483 }else{
00484 int insideVector2[3] ;
00485
00486 insideVector2[0] = ( (x-xPromInternal)>=0 ? 1:-1 ) ;
00487 insideVector2[1] = ( (y-yPromInternal)>=0 ? 1:-1 ) ;
00488 insideVector2[2] = ( (z-zPromInternal)>=0 ? 1:-1 ) ;
00489
00490 if( ( insideVector[0] == insideVector2[0] ) &&
00491 ( insideVector[1] == insideVector2[1] ) &&
00492 ( insideVector[2] == insideVector2[2] )){
00493 i++ ;
00494 }
00495 }
00496 }
00497 }
00498
00499 return ( ( i%2 == 0 ) ? false : true ) ;
00500 }
00501
00502