Surface.cpp

Go to the documentation of this file.
00001 /***************************************************************************
00002  *   Copyright (C) 2007 by Universidad Nacional de Colombia                *
00003  *   http://www.unal.edu.co                                                *
00004  *                                                                         *
00005  *   This program is free software; you can redistribute it and/or modify  *
00006  *   it under the terms of the GNU General Public License as published by  *
00007  *   the Free Software Foundation; either version 2 of the License, or     *
00008  *   (at your option) any later version.                                   *
00009  *                                                                         *
00010  *   This program is distributed in the hope that it will be useful,       *
00011  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
00012  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
00013  *   GNU General Public License for more details.                          *
00014  *                                                                         *
00015  *   You should have received a copy of the GNU General Public License     *
00016  *   along with this program; if not, write to the                         *
00017  *   Free Software Foundation, Inc.,                                       *
00018  *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
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                 // Cuadrado de las distancias entre los puntos del tri谩ngulo
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         * Construir el grafo asociado a la superficie
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         * Se construyen las aristas del grafo y al mismo tiempo
00164         * la lista de triangulos
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         * En este paso hay un problema, el algoritmo de check
00185         * siempre retorna false, as铆 que hay que ver las 
00186         * propiedades del grafo construido
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         * Este mapa se usa para no repetir triangulos 
00214         * durante la creaci贸n de cada subsuperficie
00215         * el value es solamente un valor dummy
00216         */
00217         set<Triangle> setTriangleTMP ;
00218                 
00219                 
00220         for( components::component_iterator it2 = algol.components_begin(); it2 != algol.components_end(); it2++ ){
00221                 
00222                 /**************************************************
00223                 * Solo se seleccionan los componentes con m谩s
00224                 * de 10 puntos, sino el arbol queda con mucho ruido
00225                 * HAY QUE BUSCAR LA MANERA DE QUE ESTE PARAMETRO SEA
00226                 * MODIFICABLE POR EL USUARIO!!!!!!!
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                                 * Hay que ter cuidado en esta linea ya que se esta
00236                                 * suponiendo que la numeraci贸n en el vector coincide cona la
00237                                 * numeraci贸 en el grafo
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                         * Esta parte es para llenar la lista de tri谩gulos de la
00252                         * subsuperficie creada donde ya han sido filtrados los tri谩gulos
00253                         * repetidos gracias al mapa mapTriangleTMP
00254                         */
00255                         for( set<Triangle>::iterator it4 = setTriangleTMP.begin(); it4 != setTriangleTMP.end(); it4++ )
00256                                 outputSurfaceTMP.second.push_back( *it4 ) ;
00257                         
00258                         /******************************************************************
00259                         * Se le asignan lista de superficies la que ha sido encontrada
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                 * PromInternal son las coordenas del punto de la superficie interna
00325                 * que servira para trazar la recta entre superficies. PromExternal 
00326                 * corresponde a al otro punto de la recta definido dobre la superfice 
00327                 * de afuera
00328                 */
00329 
00330                 
00331                 // Aqui entra el parametro a
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                 // Aqui entra el parametro b
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                  * Este ciclo es para evitar que xPromExternal-xPromInternal 
00360                  * sea igual a cero  y evitar que mas adelante quede una
00361                  * division por cero
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                 * A continuaci贸n se definen los parametros a, b, c y d que determinan
00385                 * la ecuaci贸n del plano definido por los puntos point1, point2 y point3
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                 * A continuaci贸n se calculan las coordenadas x, y, z que definen
00402                 * al punto "p" de intersecci贸n del plano con la recta trazada entre las dos
00403                 * superficies
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                 // Cuadrado de la distancia del punto p a los puntos que forman el tri谩ngulo
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                 // Cuadrado de las distancias entre los puntos del tri谩ngulo
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                 // Angulos formados por dos puntos del triangulo y el punto "p"
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                 * Este if determina si el punto de intereccion "p" se encuentra por 
00472                 * dentro o por fuera del triangulo, lo cual se sabe con la suma de los 
00473                 * angulos internos. Si totalAlpha es igual a 360 ...
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 

Generated on Mon May 26 20:29:46 2008 for TARIS by  doxygen 1.5.4