Back to index

salome-med  6.5.0
DirectedBoundingBox.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2009-2012  OPEN CASCADE
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 // File      : DirectedBoundingBox.cxx
00020 // Created   : Mon Apr 12 14:41:22 2010
00021 // Author    : Edward AGAPOV (eap)
00022 
00023 #include "DirectedBoundingBox.hxx"
00024 
00025 #include "InterpolationUtils.hxx"
00026 
00027 #define __TENSOR(i,j) tensor[(i)*_dim+(j)]
00028 #define __AXIS(i)     (&_axes[(i)*_dim])
00029 #define __MIN(i)      _minmax[i*2]
00030 #define __MAX(i)      _minmax[i*2+1]
00031 #define __MYID        (long(this)%10000)
00032 #define __DMP(msg) \
00033   //  cout << msg << endl
00034 
00035 using namespace std;
00036 
00037 namespace
00038 {
00039   //================================================================================
00043   //================================================================================
00044 
00045   inline void addPointToInertiaTensor3D(const double*   coord,
00046                                         const double*   gc,
00047                                         vector<double>& tensor)
00048   {
00049     // we fill the upper triangle of tensor only
00050     const int _dim = 3;
00051     double x = coord[0] - gc[0], y = coord[1] - gc[1], z = coord[2] - gc[2];
00052     __TENSOR(0,0) += y*y + z*z;
00053     __TENSOR(1,1) += x*x + z*z;
00054     __TENSOR(2,2) += x*x + y*y;
00055     __TENSOR(0,1) -= x*y;
00056     __TENSOR(0,2) -= x*z;
00057     __TENSOR(1,2) -= y*z;
00058   }
00059   //================================================================================
00063   //================================================================================
00064 
00065   inline void addPointToInertiaTensor2D(const double*   coord,
00066                                         const double*   gc,
00067                                         vector<double>& tensor)
00068   {
00069     // we fill the upper triangle of tensor only
00070     const int _dim = 2;
00071     double x = coord[0] - gc[0], y = coord[1] - gc[1];
00072     __TENSOR(0,0) += y*y;
00073     __TENSOR(1,1) += x*x;
00074     __TENSOR(0,1) -= x*y;
00075   }
00076 
00077   //================================================================================
00081   //================================================================================
00082 
00083   bool JacobiEigenvectorsSearch( const int _dim, vector<double>& tensor, vector<double>& _axes)
00084   {
00085     if ( _dim == 3 )
00086       __DMP( "Tensor : {"
00087            << "{ "<<__TENSOR(0,0) << ", "<<__TENSOR(0,1) << ", "<<__TENSOR(0,2) << "} "
00088            << "{ "<<__TENSOR(1,0) << ", "<<__TENSOR(1,1) << ", "<<__TENSOR(1,2) << "} "
00089            << "{ "<<__TENSOR(2,0) << ", "<<__TENSOR(2,1) << ", "<<__TENSOR(2,2) << "}} ");
00090     else
00091       __DMP( "Tensor : {"
00092            << "{ "<<__TENSOR(0,0) << ", "<<__TENSOR(0,1) << "} "
00093            << "{ "<<__TENSOR(1,0) << ", "<<__TENSOR(1,1) << "}} ");
00094 
00095     const int maxRot = 5*_dim*_dim; // limit on number of rotations
00096     const double tol = 1e-9;
00097 
00098     // set _axes to identity
00099     int i,j;
00100     for ( i = 0; i < _dim; ++i )
00101       for ( j = 0; j < _dim; ++j )
00102         __AXIS(i)[j] = ( i==j ? 1. : 0 );
00103 
00104     bool solved = false;
00105     for ( int iRot = 0; iRot < maxRot; ++ iRot )
00106       {
00107         // find max off-diagonal element of the tensor
00108         int k = 0, l = 0;
00109         double max = 0.;
00110         for ( i = 0; i < _dim-1; ++i )
00111           for ( j = i+1; j < _dim; ++j )
00112             if ( fabs( __TENSOR(i,j)) > max )
00113               max = fabs( __TENSOR(i,j) ), k = i, l = j;
00114         solved = ( max < tol );
00115         if ( solved )
00116           break;
00117 
00118         // Rotate to make __TENSOR(k,l) == 0
00119 
00120         double diff = __TENSOR(l,l) - __TENSOR(k,k);
00121         double t; // tangent of rotation angle
00122         if ( fabs(__TENSOR(k,l)) < abs(diff)*1.0e-36)
00123           {
00124             t = __TENSOR(k,l)/diff;
00125           }
00126         else
00127           {
00128             double phi = diff/(2.0*__TENSOR(k,l));
00129             t = 1.0/(abs(phi) + sqrt(phi*phi + 1.0));
00130             if ( phi < 0.0) t = -t;
00131           }
00132         double c = 1.0/sqrt(t*t + 1.0); // cosine of rotation angle
00133         double s = t*c; // sine of rotation angle
00134         double tau = s/(1.0 + c);
00135         __TENSOR(k,k) -= t*__TENSOR(k,l);
00136         __TENSOR(l,l) += t*__TENSOR(k,l);
00137         __TENSOR(k,l) = 0.0;
00138 
00139 #define __ROTATE(T,r1,c1,r2,c2) \
00140 { \
00141 int i1 = r1*_dim+c1, i2 = r2*_dim+c2; \
00142 double t1 = T[i1], t2 = T[i2]; \
00143 T[i1] -= s * ( t2 + tau * t1);\
00144 T[i2] += s * ( t1 - tau * t2);\
00145 }
00146         for ( i = 0; i < k; ++i )       // Case of i < k
00147           __ROTATE(tensor, i,k,i,l);
00148 
00149         for ( i = k+1; i < l; ++i )     // Case of k < i < l
00150           __ROTATE(tensor, k,i,i,l);
00151 
00152         for ( i = l + 1; i < _dim; ++i )   // Case of i > l
00153           __ROTATE(tensor, k,i,l,i);
00154 
00155         for ( i = 0; i < _dim; ++i )       // Update transformation matrix
00156           __ROTATE(_axes, i,k,i,l);
00157       }
00158 
00159     __DMP( "Solved = " << solved );
00160     if ( _dim == 3 ) {
00161       __DMP( " Eigen " << __TENSOR(0,0)<<", "<<__TENSOR(1,1)<<", "<<__TENSOR(2,2) );
00162       for ( int ii=0; ii <3; ++ii )
00163         __DMP( ii << ": " << __AXIS(ii)[0] << ", " << __AXIS(ii)[1] << ", " << __AXIS(ii)[2] );
00164     }
00165     else {
00166       __DMP( " Eigen " << __TENSOR(0,0) << ", " << __TENSOR(1,1) );
00167       for ( int ii=0; ii <2; ++ii )
00168         __DMP( ii << ": " << __AXIS(ii)[0] << ", " << __AXIS(ii)[1] );
00169     }
00170 
00171     return solved;
00172   }
00173 
00174   //================================================================================
00178   //================================================================================
00179 
00180   inline bool isMinMaxOut(const double* minmax1,
00181                           const double* minmax2,
00182                           int           dim)
00183   {
00184     for ( int i = 0; i < dim; ++i )
00185       {
00186         if ( minmax1[i*2] > minmax2[i*2+1] ||
00187              minmax1[i*2+1] < minmax2[i*2] )
00188           return true;
00189       }
00190     return false;
00191   }
00192 
00193 } // noname namespace
00194 
00195 namespace INTERP_KERNEL
00196 {
00197 
00198   //================================================================================
00202   //================================================================================
00203 
00204   DirectedBoundingBox::DirectedBoundingBox():_dim(0)
00205   {
00206   }
00207 
00208   //================================================================================
00215   //================================================================================
00216 
00217   DirectedBoundingBox::DirectedBoundingBox(const double*  pts,
00218                                            const unsigned numPts,
00219                                            const unsigned dim)
00220     : _dim(dim), _axes(dim*dim), _minmax(2*dim)
00221   {
00222     // init box extremities
00223     for ( unsigned i = 0; i < _dim; ++i )
00224       _minmax[1+i*2] = -numeric_limits<double>::max(),
00225         _minmax[i*2] =  numeric_limits<double>::max();
00226 
00227     if ( numPts < 1 ) return;
00228 
00229     __DMP( "DirectedBoundingBox " << __MYID );
00230 
00231     const double* coord = pts;
00232     const double* coordEnd = coord + numPts * dim;
00233 
00234     // compute gravity center of points
00235     double gc[3] = {0,0,0};
00236     if ( dim > 1 )
00237       {
00238         for ( coord = pts; coord < coordEnd; )
00239           for ( int i = 0; i < (int)dim; ++i )
00240             gc[i] += *coord++;
00241         for ( int j = 0; j < (int)dim; ++j )
00242           gc[j] /= numPts;
00243 
00244       }
00245 
00246     // compute axes and box extremities
00247     vector<double> tensor( dim * dim, 0.);
00248     switch ( dim )
00249       {
00250       case 3:
00251         for ( coord = pts; coord < coordEnd; coord += dim )
00252           addPointToInertiaTensor3D( coord, gc, tensor );
00253 
00254         //computeAxes3D(tensor);
00255         JacobiEigenvectorsSearch(_dim, tensor, _axes);
00256 
00257         for ( coord = pts; coord < coordEnd; coord += dim )
00258           addPointToBox( coord );
00259 
00260         break;
00261 
00262       case 2:
00263         for ( coord = pts; coord < coordEnd; coord += dim )
00264           addPointToInertiaTensor2D( coord, gc, tensor );
00265 
00266         //computeAxes2D(tensor);
00267         JacobiEigenvectorsSearch(_dim, tensor, _axes);
00268 
00269         for ( coord = pts; coord < coordEnd; coord += dim )
00270           addPointToBox( coord );
00271 
00272         break;
00273 
00274       default:
00275         for ( coord = pts; coord < coordEnd; coord += dim )
00276           {
00277             if ( *coord < _minmax[0] ) _minmax[0] = *coord;
00278             if ( *coord > _minmax[1] ) _minmax[1] = *coord;
00279           }
00280       }
00281   }
00282 
00283   //================================================================================
00290   //================================================================================
00291 
00292   DirectedBoundingBox::DirectedBoundingBox(const double** pts,
00293                                            const unsigned numPts,
00294                                            const unsigned dim)
00295     : _dim(dim), _axes(dim*dim), _minmax(2*dim)
00296   {
00297     // init box extremities
00298     for ( unsigned i = 0; i < _dim; ++i )
00299       _minmax[1+i*2] = -numeric_limits<double>::max(),
00300         _minmax[i*2] =  numeric_limits<double>::max();
00301 
00302     if ( numPts < 1 ) return;
00303 
00304     __DMP( "DirectedBoundingBox " << __MYID );
00305 
00306     // compute gravity center of points
00307     double gc[3] = {0,0,0};
00308     if ( dim > 1 )
00309       {
00310         for ( unsigned i = 0; i < numPts; ++i )
00311           for ( int j = 0; j < (int)dim; ++j )
00312             gc[j] += pts[i][j];
00313         for ( int j = 0; j < (int)dim; ++j )
00314           gc[j] /= numPts;
00315       }
00316 
00317     // compute axes and box extremities
00318     vector<double> tensor( dim * dim, 0.);
00319     switch ( dim )
00320       {
00321       case 3:
00322         for ( unsigned i = 0; i < numPts; ++i )
00323           addPointToInertiaTensor3D( pts[i], gc, tensor );
00324 
00325         //computeAxes3D(tensor);
00326         JacobiEigenvectorsSearch(_dim, tensor, _axes);
00327 
00328         for ( unsigned i = 0; i < numPts; ++i )
00329           addPointToBox( pts[i] );
00330 
00331         break;
00332       case 2:
00333         for ( unsigned i = 0; i < numPts; ++i )
00334           addPointToInertiaTensor2D( pts[i], gc, tensor );
00335 
00336         //computeAxes2D(tensor);
00337         JacobiEigenvectorsSearch(_dim, tensor, _axes);
00338 
00339         for ( unsigned i = 0; i < numPts; ++i )
00340           addPointToBox( pts[i] );
00341 
00342         break;
00343       default:
00344         for ( unsigned i = 0; i < numPts; ++i )
00345           {
00346             if ( pts[i][0] < _minmax[0] ) _minmax[0] = pts[i][0];
00347             if ( pts[i][0] > _minmax[1] ) _minmax[1] = pts[i][0];
00348           }
00349         _axes[0] = 1.0;
00350       }
00351   }
00352 
00353   //================================================================================
00357   //================================================================================
00358 
00359   // void DirectedBoundingBox::computeAxes3D(const std::vector<double>& tensor)
00360 //   {
00361 //     // compute principal moments of inertia which are eigenvalues of the tensor
00362 //     double eig[3];
00363 //     {
00364 //       // coefficients of polynomial equation det(tensor-eig*I) = 0
00365 //       double a = -1;
00366 //       double b = __TENSOR(0,0)+__TENSOR(1,1)+__TENSOR(2,2);
00367 //       double c =
00368 //         __TENSOR(0,1)*__TENSOR(0,1) +
00369 //         __TENSOR(0,2)*__TENSOR(0,2) +
00370 //         __TENSOR(1,2)*__TENSOR(1,2) -
00371 //         __TENSOR(0,0)*__TENSOR(1,1) -
00372 //         __TENSOR(0,0)*__TENSOR(2,2) -
00373 //         __TENSOR(1,1)*__TENSOR(2,2);
00374 //       double d =
00375 //         __TENSOR(0,0)*__TENSOR(1,1)*__TENSOR(2,2) -
00376 //         __TENSOR(0,0)*__TENSOR(1,2)*__TENSOR(1,2) -
00377 //         __TENSOR(1,1)*__TENSOR(0,2)*__TENSOR(0,2) -
00378 //         __TENSOR(2,2)*__TENSOR(0,1)*__TENSOR(0,1) +
00379 //         __TENSOR(0,1)*__TENSOR(0,2)*__TENSOR(1,2)*2;
00380 
00381 //       // find eigenvalues which are roots of characteristic polynomial
00382 //       double x = (3*c/a - b*b/(a*a))/3;
00383 //       double y = (2*b*b*b/(a*a*a) - 9*b*c/(a*a) + 27*d/a)/27;
00384 //       double z = y*y/4 + x*x*x/27;
00385 
00386 //       double i = sqrt(y*y/4 - z) + 1e-300;
00387 //       double j = -pow(i,1/3.);
00388 //       double y2 = -y/(2*i);
00389 //       if ( y2 > 1.0) y2 = 1.; else if ( y2 < -1.0) y2 = -1.;
00390 //       double k = acos(y2);
00391 //       double m = cos(k/3);
00392 //       double n = sqrt(3)*sin(k/3);
00393 //       double p = -b/(3*a);
00394 
00395 //       eig[0] = -2*j*m + p;
00396 //       eig[1] = j *(m + n) + p;
00397 //       eig[2] = j *(m - n) + p;
00398 //     }
00399 //     // compute eigenvector of the tensor at each eigenvalue
00400 //     // by solving system [tensor-eig*I]*[axis] = 0
00401 //     bool ok = true;
00402 //     __DMP( "Tensor : {"
00403 //          << "{ "<<__TENSOR(0,0) << ", "<<__TENSOR(0,1) << ", "<<__TENSOR(0,2) << "} "
00404 //          << "{ "<<__TENSOR(1,0) << ", "<<__TENSOR(1,1) << ", "<<__TENSOR(1,2) << "} "
00405 //          << "{ "<<__TENSOR(2,0) << ", "<<__TENSOR(2,1) << ", "<<__TENSOR(2,2) << "}} ");
00406 //     for ( int i = 0; i < 3 && ok; ++i ) // loop on 3 eigenvalues
00407 //       {
00408 //         // [tensor-eig*I]
00409 //         double T[3][3]=
00410 //           {{ __TENSOR(0,0)-eig[i],__TENSOR(0,1),       __TENSOR(0,2),      },
00411 //            { __TENSOR(0,1),       __TENSOR(1,1)-eig[i],__TENSOR(1,2),      },
00412 //            { __TENSOR(0,2),       __TENSOR(1,2),       __TENSOR(2,2)-eig[i]}};
00413 //         // The determinant of T is zero, so that the equations are not linearly independent.
00414 //         // Therefore, we assign an arbitrary value (1.) to i-th component of eigenvector
00415 //         // and use two of the equations to compute the other two components
00416 //         double M[2][3], sol[2];
00417 //         for ( int j = 0, c = 0; j < 3; ++j )
00418 //           if ( i == j )
00419 //             M[0][2] = -T[0][j], M[1][2] = -T[1][j];
00420 //           else
00421 //             M[0][c] = T[0][j], M[1][c] = T[1][j], c++;
00422 
00423 //         ok = solveSystemOfEquations<2>( M, sol );
00424 
00425 //         double* eigenVec = __AXIS(i);
00426 //         for ( int j = 0, c = 0; j < 3; ++j )
00427 //           eigenVec[j] = ( i == j ) ? 1. : sol[c++];
00428 
00429 //         // normilize
00430 //         double size = sqrt(eigenVec[0]*eigenVec[0] +
00431 //                            eigenVec[1]*eigenVec[1] +
00432 //                            eigenVec[2]*eigenVec[2] );
00433 //         if ((ok = (size > numeric_limits<double>::min() )))
00434 //           {
00435 //             eigenVec[0] /= size;
00436 //             eigenVec[1] /= size;
00437 //             eigenVec[2] /= size;
00438 //           }
00439 //       }
00440 //     if ( !ok )
00441 //       {
00442 //         __DMP( " solve3EquationSystem() - KO " );
00443 //         _axes = vector<double>( _dim*_dim, 0);
00444 //         __AXIS(0)[0] = __AXIS(1)[1] = __AXIS(2)[2] = 1.;
00445 //       }
00446 //     __DMP( " Eigen " << eig[0] << ", " << eig[1] << ", " << eig[2] );
00447 //     for ( int i=0; i <3; ++i )
00448 //       __DMP( i << ": " << __AXIS(i)[0] << ", " << __AXIS(i)[1] << ", " << __AXIS(i)[2] );
00449 
00450 //     double* a0 = __AXIS(0), *a1 = __AXIS(1);
00451 //     double cross[3] = { a0[1]*a1[2]-a1[1]*a0[2],
00452 //                         a0[2]*a1[0]-a1[2]*a0[0],
00453 //                         a0[0]*a1[1]-a1[0]*a0[1] };
00454 //     __DMP( " Cross a1^a2 " << cross[0] << ", " << cross[1] << ", " << cross[2] );
00455 //   }
00456 
00457   //================================================================================
00461   //================================================================================
00462 
00463   // void DirectedBoundingBox::computeAxes2D(const std::vector<double>& tensor)
00464 //   {
00465 //     // compute principal moments of inertia which are eigenvalues of the tensor
00466 //     // by solving square equation det(tensor-eig*I)
00467 //     double X = (__TENSOR(0,0)+__TENSOR(1,1))/2;
00468 //     double Y = sqrt(4*__TENSOR(0,1)*__TENSOR(0,1) +
00469 //                     (__TENSOR(0,0)-__TENSOR(1,1)) * (__TENSOR(0,0)-__TENSOR(1,1)))/2;
00470 //     double eig[2] =
00471 //       {
00472 //         X + Y,
00473 //         X - Y
00474 //       };
00475 //     // compute eigenvector of the tensor at each eigenvalue
00476 //     // by solving system [tensor-eig*I]*[axis] = 0
00477 //     bool ok = true;
00478 //     for ( int i = 0; i < 2 && ok; ++i )
00479 //       {
00480 //         // [tensor-eig*I]
00481 //         double T[2][2]=
00482 //           {{ __TENSOR(0,0)-eig[i],__TENSOR(0,1)        },
00483 //            { __TENSOR(0,1),       __TENSOR(1,1)-eig[i] }};
00484 
00485 //         // The determinant of T is zero, so that the equations are not linearly independent.
00486 //         // Therefore, we assign an arbitrary value (1.) to i-th component of eigenvector
00487 //         // and use one equation to compute the other component
00488 //         double* eigenVec = __AXIS(i);
00489 //         eigenVec[i] = 1.;
00490 //         int j = 1-i;
00491 //         if ((ok = ( fabs( T[j][j] ) > numeric_limits<double>::min() )))
00492 //           eigenVec[j] = -T[j][i] / T[j][j];
00493 //       }
00494 //     if ( !ok )
00495 //       {
00496 //         _axes = vector<double>( _dim*_dim, 0);
00497 //         __AXIS(0)[0] = __AXIS(1)[1] = 1.;
00498 //       }
00499 //   }
00500 
00501   //================================================================================
00505   //================================================================================
00506 
00507   void DirectedBoundingBox::toLocalCS(const double* p, double* pLoc) const
00508   {
00509     switch ( _dim )
00510       {
00511       case 3:
00512         pLoc[0] = dotprod<3>( p, __AXIS(0));
00513         pLoc[1] = dotprod<3>( p, __AXIS(1));
00514         pLoc[2] = dotprod<3>( p, __AXIS(2));
00515         break;
00516       case 2:
00517         pLoc[0] = dotprod<2>( p, __AXIS(0));
00518         pLoc[1] = dotprod<2>( p, __AXIS(1));
00519         break;
00520       default:
00521         pLoc[0] = p[0];
00522       }
00523   }
00524 
00525   //================================================================================
00529   //================================================================================
00530 
00531   void DirectedBoundingBox::fromLocalCS(const double* p, double* pGlob) const
00532   {
00533     switch ( _dim )
00534       {
00535       case 3:
00536         pGlob[0] = p[0] * __AXIS(0)[0] + p[1] * __AXIS(1)[0] + p[2] * __AXIS(2)[0];
00537         pGlob[1] = p[0] * __AXIS(0)[1] + p[1] * __AXIS(1)[1] + p[2] * __AXIS(2)[1];
00538         pGlob[2] = p[0] * __AXIS(0)[2] + p[1] * __AXIS(1)[2] + p[2] * __AXIS(2)[2];
00539         break;
00540       case 2:
00541         pGlob[0] = p[0] * __AXIS(0)[0] + p[1] * __AXIS(1)[0];
00542         pGlob[1] = p[0] * __AXIS(0)[1] + p[1] * __AXIS(1)[1];
00543         break;
00544       default:
00545         pGlob[0] = p[0];
00546       }
00547   }
00548 
00549   //================================================================================
00553   //================================================================================
00554 
00555   void DirectedBoundingBox::enlarge(const double tol)
00556   {
00557     for ( unsigned i = 0; i < _dim; ++i )
00558       __MIN(i) -= tol, __MAX(i) += tol;
00559   }
00560 
00561   //================================================================================
00565   //================================================================================
00566 
00567   void DirectedBoundingBox::getCorners(std::vector<double>& corners,
00568                                        const double*        minmax) const
00569   {
00570     int iC, nbCorners = 1;
00571     for ( int i=0;i<(int)_dim;++i ) nbCorners *= 2;
00572     corners.resize( nbCorners * _dim );
00573     // each coordinate is filled with either min or max, nbSwap is number of corners
00574     // after which min and max swap
00575     int nbSwap = nbCorners/2;
00576     for ( unsigned i = 0; i < _dim; ++i )
00577       {
00578         iC = 0;
00579         while ( iC < nbCorners )
00580           {
00581             for (int j = 0; j < nbSwap; ++j, ++iC ) corners[iC*_dim+i] = minmax[i*2];
00582             for (int j = 0; j < nbSwap; ++j, ++iC ) corners[iC*_dim+i] = minmax[i*2+1];
00583           }
00584         nbSwap /= 2;
00585       }
00586   }
00587 
00588   //================================================================================
00593   //================================================================================
00594 
00595   bool DirectedBoundingBox::isDisjointWith(const DirectedBoundingBox& box) const
00596   {
00597     if ( _dim < 1 || box._dim < 1 ) return false; // empty box includes all
00598     if ( _dim == 1 )
00599       return isMinMaxOut( &box._minmax[0], &this->_minmax[0], _dim );
00600 
00601     // boxes are disjoined if their minmaxes in local CS of either of boxes do not intersect
00602     for ( int isThisCS = 0; isThisCS < 2; ++isThisCS )
00603       {
00604         const DirectedBoundingBox* axisBox   = isThisCS ? this : &box;
00605         const DirectedBoundingBox* cornerBox = isThisCS ? &box : this;
00606 
00607         // find minmax of cornerBox in the CS of axisBox
00608 
00609         DirectedBoundingBox mmBox((double*)0,0,_dim); 
00610         mmBox._axes = axisBox->_axes;
00611 
00612         vector<double> corners;
00613         getCorners( corners, &cornerBox->_minmax[0] );
00614 
00615         double globCorner[3];
00616         for ( int iC = 0, nC = corners.size()/_dim; iC < nC; ++iC)
00617           {
00618             cornerBox->fromLocalCS( &corners[iC*_dim], globCorner );
00619             mmBox.addPointToBox( globCorner );
00620           }
00621         if ( isMinMaxOut( &mmBox._minmax[0], &axisBox->_minmax[0], _dim ))
00622           return true;
00623       }
00624     return false;
00625   }
00626 
00627   //================================================================================
00632   //================================================================================
00633 
00634   bool DirectedBoundingBox::isDisjointWith(const double* box) const
00635   {
00636     if ( _dim < 1 ) return false; // empty box includes all
00637     if ( _dim == 1 )
00638       return isMinMaxOut( &_minmax[0], box, _dim );
00639 
00640     // boxes are disjoined if their minmaxes in local CS of either of boxes do not intersect
00641 
00642     // compare minmaxes in locals CS of this directed box
00643     {
00644       vector<double> cornersOther;
00645       getCorners( cornersOther, box );
00646       DirectedBoundingBox mmBox((double*)0,0,_dim); 
00647       mmBox._axes = this->_axes;
00648       for ( int iC = 0, nC = cornersOther.size()/_dim; iC < nC; ++iC)
00649         mmBox.addPointToBox( &cornersOther[iC*_dim] );
00650 
00651       if ( isMinMaxOut( &mmBox._minmax[0], &this->_minmax[0], _dim ))
00652         return true;
00653     }
00654 
00655     // compare minmaxes in global CS
00656     {
00657       vector<double> cornersThis;
00658       getCorners( cornersThis, &_minmax[0] );
00659       DirectedBoundingBox mmBox((double*)0,0,_dim); 
00660       double globCorner[3];
00661       for ( int iC = 0, nC = cornersThis.size()/_dim; iC < nC; ++iC)
00662         {
00663           fromLocalCS( &cornersThis[iC*_dim], globCorner );
00664           for ( int i = 0; i < (int)_dim; ++i )
00665             {
00666               if ( globCorner[i] < mmBox._minmax[i*2] )   mmBox._minmax[i*2] = globCorner[i];
00667               if ( globCorner[i] > mmBox._minmax[i*2+1] ) mmBox._minmax[i*2+1] = globCorner[i];
00668             }
00669         }
00670       if ( isMinMaxOut( &mmBox._minmax[0], box, _dim ))
00671         return true;
00672     }
00673     return false;
00674   }
00675 
00676   //================================================================================
00680   //================================================================================
00681 
00682   bool DirectedBoundingBox::isOut(const double* point) const
00683   {
00684     if ( _dim < 1 ) return false; // empty box includes all
00685 
00686     double pLoc[3];
00687     toLocalCS( point, pLoc );
00688     bool out = isLocalOut( pLoc );
00689 #ifdef _DEBUG_
00690     switch (_dim)
00691       {
00692       case 3:
00693         __DMP(__MYID<<": "<<point[0]<<", "<<point[1]<<", "<<point[2]<<" "<<(out?"OUT":"IN"));break;
00694       case 2:
00695         __DMP(__MYID<<": "<<point[0]<<", "<<point[1]<<" "<<(out?"OUT":"IN"));break;
00696       case 1:
00697         __DMP(__MYID<<": "<<point[0]<<" "<<(out?"OUT":"IN"));break;
00698       }
00699 #endif
00700     return out;
00701   }
00702 
00703   //================================================================================
00707   //================================================================================
00708 
00709   vector<double> DirectedBoundingBox::getData() const
00710   {
00711     vector<double> data(1, _dim);
00712     if ( _dim > 0 )
00713     {
00714       data.insert( data.end(), &_axes[0], &_axes[0] + _axes.size());
00715       data.insert( data.end(), &_minmax[0], &_minmax[0] + _minmax.size());
00716     }
00717     if ( data.size() < (unsigned)dataSize( _dim ))
00718       data.resize( dataSize( _dim ), 0 );
00719     return data;
00720   }
00721 
00722   //================================================================================
00726   //================================================================================
00727 
00728   void DirectedBoundingBox::setData(const double* data)
00729   {
00730     _dim = unsigned( *data++ );
00731     if ( _dim > 0 )
00732       {
00733         _axes.assign( data, data+_dim*_dim ); data += _dim*_dim;
00734         _minmax.assign( data, data+2*_dim );
00735       }
00736     else
00737       {
00738         _axes.clear();
00739         _minmax.clear();
00740       }
00741   }
00742 
00743   //================================================================================
00747   //================================================================================
00748 
00749   int DirectedBoundingBox::dataSize(int dim)
00750   {
00751     return 1 + dim*dim + 2*dim; // : _dim + _axes + _minmax
00752   }
00753 }