src/INTERP_KERNEL/DirectedBoundingBox.cxx

   1 // Copyright (C) 2009-2021  OPEN CASCADE
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // File      : DirectedBoundingBox.cxx
  20 // Created   : Mon Apr 12 14:41:22 2010
  21 // Author    : Edward AGAPOV (eap)
  22
  23 #include "DirectedBoundingBox.hxx"
  24
  25 #include "InterpolationUtils.hxx"
  26
  27 #define __TENSOR(i,j) tensor[(i)*_dim+(j)]
  28 #define __AXIS(i)     (&_axes[(i)*_dim])
  29 #define __MIN(i)      _minmax[i*2]
  30 #define __MAX(i)      _minmax[i*2+1]
  31 #define __DMP(msg) \
  32   //  cout << msg << endl
  33
  34 using namespace std;
  35
  36 namespace
  37 {
  38   //================================================================================
  39   /*!
  40    * \brief Add point coordinates to inertia tensor in 3D space
  41    */
  42   //================================================================================
  43
  44   inline void addPointToInertiaTensor3D(const double*   coord,
  45                                         const double*   gc,
  46                                         vector<double>& tensor)
  47   {
  48     // we fill the upper triangle of tensor only
  49     const int _dim = 3;
  50     double x = coord[0] - gc[0], y = coord[1] - gc[1], z = coord[2] - gc[2];
  51     __TENSOR(0,0) += y*y + z*z;
  52     __TENSOR(1,1) += x*x + z*z;
  53     __TENSOR(2,2) += x*x + y*y;
  54     __TENSOR(0,1) -= x*y;
  55     __TENSOR(0,2) -= x*z;
  56     __TENSOR(1,2) -= y*z;
  57   }
  58   //================================================================================
  59   /*!
  60    * \brief Add point coordinates to inertia tensor in 2D space
  61    */
  62   //================================================================================
  63
  64   inline void addPointToInertiaTensor2D(const double*   coord,
  65                                         const double*   gc,
  66                                         vector<double>& tensor)
  67   {
  68     // we fill the upper triangle of tensor only
  69     const int _dim = 2;
  70     double x = coord[0] - gc[0], y = coord[1] - gc[1];
  71     __TENSOR(0,0) += y*y;
  72     __TENSOR(1,1) += x*x;
  73     __TENSOR(0,1) -= x*y;
  74   }
  75
  76   //================================================================================
  77   /*!
  78    * \brief Find eigenvectors of tensor using Jacobi's method
  79    */
  80   //================================================================================
  81
  82   bool JacobiEigenvectorsSearch( const int _dim, vector<double>& tensor, vector<double>& _axes)
  83   {
  84     if ( _dim == 3 )
  85       {
  86         __DMP( "Tensor : {"
  87                << "{ "<<__TENSOR(0,0) << ", "<<__TENSOR(0,1) << ", "<<__TENSOR(0,2) << "} "
  88                << "{ "<<__TENSOR(1,0) << ", "<<__TENSOR(1,1) << ", "<<__TENSOR(1,2) << "} "
  89                << "{ "<<__TENSOR(2,0) << ", "<<__TENSOR(2,1) << ", "<<__TENSOR(2,2) << "}} ");
  90       }
  91     else
  92       {
  93         __DMP( "Tensor : {"
  94           << "{ "<<__TENSOR(0,0) << ", "<<__TENSOR(0,1) << "} "
  95           << "{ "<<__TENSOR(1,0) << ", "<<__TENSOR(1,1) << "}} ");
  96       }
  97
  98     const int maxRot = 5*_dim*_dim; // limit on number of rotations
  99     const double tol = 1e-9;
 100
 101     // set _axes to identity
 102     int i,j;
 103     for ( i = 0; i < _dim; ++i )
 104       for ( j = 0; j < _dim; ++j )
 105         __AXIS(i)[j] = ( i==j ? 1. : 0 );
 106
 107     bool solved = false;
 108     for ( int iRot = 0; iRot < maxRot; ++ iRot )
 109       {
 110         // find max off-diagonal element of the tensor
 111         int k = 0, l = 0;
 112         double max = 0.;
 113         for ( i = 0; i < _dim-1; ++i )
 114           for ( j = i+1; j < _dim; ++j )
 115             if ( fabs( __TENSOR(i,j)) > max )
 116               max = fabs( __TENSOR(i,j) ), k = i, l = j;
 117         solved = ( max < tol );
 118         if ( solved )
 119           break;
 120
 121         // Rotate to make __TENSOR(k,l) == 0
 122
 123         double diff = __TENSOR(l,l) - __TENSOR(k,k);
 124         double t; // tangent of rotation angle
 125         if ( fabs(__TENSOR(k,l)) < abs(diff)*1.0e-36)
 126           {
 127             t = __TENSOR(k,l)/diff;
 128           }
 129         else
 130           {
 131             double phi = diff/(2.0*__TENSOR(k,l));
 132             t = 1.0/(abs(phi) + sqrt(phi*phi + 1.0));
 133             if ( phi < 0.0) t = -t;
 134           }
 135         double c = 1.0/sqrt(t*t + 1.0); // cosine of rotation angle
 136         double s = t*c; // sine of rotation angle
 137         double tau = s/(1.0 + c);
 138         __TENSOR(k,k) -= t*__TENSOR(k,l);
 139         __TENSOR(l,l) += t*__TENSOR(k,l);
 140         __TENSOR(k,l) = 0.0;
 141
 142 #define __ROTATE(T,r1,c1,r2,c2) \
 143 { \
 144 int i1 = r1*_dim+c1, i2 = r2*_dim+c2; \
 145 double t1 = T[i1], t2 = T[i2]; \
 146 T[i1] -= s * ( t2 + tau * t1);\
 147 T[i2] += s * ( t1 - tau * t2);\
 148 }
 149         for ( i = 0; i < k; ++i )       // Case of i < k
 150           __ROTATE(tensor, i,k,i,l);
 151
 152         for ( i = k+1; i < l; ++i )     // Case of k < i < l
 153           __ROTATE(tensor, k,i,i,l);
 154
 155         for ( i = l + 1; i < _dim; ++i )   // Case of i > l
 156           __ROTATE(tensor, k,i,l,i);
 157
 158         for ( i = 0; i < _dim; ++i )       // Update transformation matrix
 159           __ROTATE(_axes, i,k,i,l);
 160       }
 161
 162     __DMP( "Solved = " << solved );
 163     if ( _dim == 3 ) {
 164       __DMP( " Eigen " << __TENSOR(0,0)<<", "<<__TENSOR(1,1)<<", "<<__TENSOR(2,2) );
 165       for ( int ii=0; ii <3; ++ii )
 166         __DMP( ii << ": " << __AXIS(ii)[0] << ", " << __AXIS(ii)[1] << ", " << __AXIS(ii)[2] );
 167     }
 168     else {
 169       __DMP( " Eigen " << __TENSOR(0,0) << ", " << __TENSOR(1,1) );
 170       for ( int ii=0; ii <2; ++ii )
 171         __DMP( ii << ": " << __AXIS(ii)[0] << ", " << __AXIS(ii)[1] );
 172     }
 173
 174     return solved;
 175   }
 176
 177   //================================================================================
 178   /*!
 179    * \brief Return true if two minmaxes do not intersect
 180    */
 181   //================================================================================
 182
 183   inline bool isMinMaxOut(const double* minmax1,
 184                           const double* minmax2,
 185                           int           dim)
 186   {
 187     for ( int i = 0; i < dim; ++i )
 188       {
 189         if ( minmax1[i*2] > minmax2[i*2+1] ||
 190              minmax1[i*2+1] < minmax2[i*2] )
 191           return true;
 192       }
 193     return false;
 194   }
 195
 196 } // noname namespace
 197
 198 namespace INTERP_KERNEL
 199 {
 200
 201   //================================================================================
 202   /*!
 203    * \brief Creates empty box intended to further initialization via setData()
 204    */
 205   //================================================================================
 206
 207   DirectedBoundingBox::DirectedBoundingBox():_dim(0)
 208   {
 209   }
 210
 211   //================================================================================
 212   /*!
 213    * \brief Creates bounding box of a mesh
 214    *  \param pts - coordinates of points in full interlace
 215    *  \param numPts - number of points in the mesh
 216    *  \param dim - space dimension
 217    */
 218   //================================================================================
 219
 220   DirectedBoundingBox::DirectedBoundingBox(const double*  pts,
 221                                            const unsigned numPts,
 222                                            const unsigned dim)
 223     : _dim(dim), _axes(dim*dim), _minmax(2*dim)
 224   {
 225     // init box extremities
 226     for ( unsigned i = 0; i < _dim; ++i )
 227       _minmax[1+i*2] = -numeric_limits<double>::max(),
 228         _minmax[i*2] =  numeric_limits<double>::max();
 229
 230     if ( numPts < 1 ) return;
 231
 232     __DMP( "DirectedBoundingBox " << __MYID );
 233
 234     const double* coord = pts;
 235     const double* coordEnd = coord + numPts * dim;
 236
 237     // compute gravity center of points
 238     double gc[3] = {0,0,0};
 239     if ( dim > 1 )
 240       {
 241         for ( coord = pts; coord < coordEnd; )
 242           for ( int i = 0; i < (int)dim; ++i )
 243             gc[i] += *coord++;
 244         for ( int j = 0; j < (int)dim; ++j )
 245           gc[j] /= numPts;
 246
 247       }
 248
 249     // compute axes and box extremities
 250     vector<double> tensor( dim * dim, 0.);
 251     switch ( dim )
 252       {
 253       case 3:
 254         for ( coord = pts; coord < coordEnd; coord += dim )
 255           addPointToInertiaTensor3D( coord, gc, tensor );
 256
 257         //computeAxes3D(tensor);
 258         JacobiEigenvectorsSearch(_dim, tensor, _axes);
 259
 260         for ( coord = pts; coord < coordEnd; coord += dim )
 261           addPointToBox( coord );
 262
 263         break;
 264
 265       case 2:
 266         for ( coord = pts; coord < coordEnd; coord += dim )
 267           addPointToInertiaTensor2D( coord, gc, tensor );
 268
 269         //computeAxes2D(tensor);
 270         JacobiEigenvectorsSearch(_dim, tensor, _axes);
 271
 272         for ( coord = pts; coord < coordEnd; coord += dim )
 273           addPointToBox( coord );
 274
 275         break;
 276
 277       default:
 278         for ( coord = pts; coord < coordEnd; coord += dim )
 279           {
 280             if ( *coord < _minmax[0] ) _minmax[0] = *coord;
 281             if ( *coord > _minmax[1] ) _minmax[1] = *coord;
 282           }
 283       }
 284   }
 285
 286   //================================================================================
 287   /*!
 288    * \brief Creates bounding box of an element
 289    *  \param pts - coordinates of points of element
 290    *  \param numPts - number of points in the element
 291    *  \param dim - space dimension
 292    */
 293   //================================================================================
 294
 295   DirectedBoundingBox::DirectedBoundingBox(const double** pts,
 296                                            const unsigned numPts,
 297                                            const unsigned dim)
 298     : _dim(dim), _axes(dim*dim), _minmax(2*dim)
 299   {
 300     // init box extremities
 301     for ( unsigned i = 0; i < _dim; ++i )
 302       _minmax[1+i*2] = -numeric_limits<double>::max(),
 303         _minmax[i*2] =  numeric_limits<double>::max();
 304
 305     if ( numPts < 1 ) return;
 306
 307     __DMP( "DirectedBoundingBox " << __MYID );
 308
 309     // compute gravity center of points
 310     double gc[3] = {0,0,0};
 311     if ( dim > 1 )
 312       {
 313         for ( unsigned i = 0; i < numPts; ++i )
 314           for ( int j = 0; j < (int)dim; ++j )
 315             gc[j] += pts[i][j];
 316         for ( int j = 0; j < (int)dim; ++j )
 317           gc[j] /= numPts;
 318       }
 319
 320     // compute axes and box extremities
 321     vector<double> tensor( dim * dim, 0.);
 322     switch ( dim )
 323       {
 324       case 3:
 325         for ( unsigned i = 0; i < numPts; ++i )
 326           addPointToInertiaTensor3D( pts[i], gc, tensor );
 327
 328         //computeAxes3D(tensor);
 329         JacobiEigenvectorsSearch(_dim, tensor, _axes);
 330
 331         for ( unsigned i = 0; i < numPts; ++i )
 332           addPointToBox( pts[i] );
 333
 334         break;
 335       case 2:
 336         for ( unsigned i = 0; i < numPts; ++i )
 337           addPointToInertiaTensor2D( pts[i], gc, tensor );
 338
 339         //computeAxes2D(tensor);
 340         JacobiEigenvectorsSearch(_dim, tensor, _axes);
 341
 342         for ( unsigned i = 0; i < numPts; ++i )
 343           addPointToBox( pts[i] );
 344
 345         break;
 346       default:
 347         for ( unsigned i = 0; i < numPts; ++i )
 348           {
 349             if ( pts[i][0] < _minmax[0] ) _minmax[0] = pts[i][0];
 350             if ( pts[i][0] > _minmax[1] ) _minmax[1] = pts[i][0];
 351           }
 352         _axes[0] = 1.0;
 353       }
 354   }
 355
 356   //================================================================================
 357   /*!
 358    * \brief Compute eigenvectors of inertia tensor
 359    */
 360   //================================================================================
 361
 362   // void DirectedBoundingBox::computeAxes3D(const std::vector<double>& tensor)
 363 //   {
 364 //     // compute principal moments of inertia which are eigenvalues of the tensor
 365 //     double eig[3];
 366 //     {
 367 //       // coefficients of polynomial equation det(tensor-eig*I) = 0
 368 //       double a = -1;
 369 //       double b = __TENSOR(0,0)+__TENSOR(1,1)+__TENSOR(2,2);
 370 //       double c =
 371 //         __TENSOR(0,1)*__TENSOR(0,1) +
 372 //         __TENSOR(0,2)*__TENSOR(0,2) +
 373 //         __TENSOR(1,2)*__TENSOR(1,2) -
 374 //         __TENSOR(0,0)*__TENSOR(1,1) -
 375 //         __TENSOR(0,0)*__TENSOR(2,2) -
 376 //         __TENSOR(1,1)*__TENSOR(2,2);
 377 //       double d =
 378 //         __TENSOR(0,0)*__TENSOR(1,1)*__TENSOR(2,2) -
 379 //         __TENSOR(0,0)*__TENSOR(1,2)*__TENSOR(1,2) -
 380 //         __TENSOR(1,1)*__TENSOR(0,2)*__TENSOR(0,2) -
 381 //         __TENSOR(2,2)*__TENSOR(0,1)*__TENSOR(0,1) +
 382 //         __TENSOR(0,1)*__TENSOR(0,2)*__TENSOR(1,2)*2;
 383
 384 //       // find eigenvalues which are roots of characteristic polynomial
 385 //       double x = (3*c/a - b*b/(a*a))/3;
 386 //       double y = (2*b*b*b/(a*a*a) - 9*b*c/(a*a) + 27*d/a)/27;
 387 //       double z = y*y/4 + x*x*x/27;
 388
 389 //       double i = sqrt(y*y/4 - z) + 1e-300;
 390 //       double j = -pow(i,1/3.);
 391 //       double y2 = -y/(2*i);
 392 //       if ( y2 > 1.0) y2 = 1.; else if ( y2 < -1.0) y2 = -1.;
 393 //       double k = acos(y2);
 394 //       double m = cos(k/3);
 395 //       double n = sqrt(3)*sin(k/3);
 396 //       double p = -b/(3*a);
 397
 398 //       eig[0] = -2*j*m + p;
 399 //       eig[1] = j *(m + n) + p;
 400 //       eig[2] = j *(m - n) + p;
 401 //     }
 402 //     // compute eigenvector of the tensor at each eigenvalue
 403 //     // by solving system [tensor-eig*I]*[axis] = 0
 404 //     bool ok = true;
 405 //     __DMP( "Tensor : {"
 406 //          << "{ "<<__TENSOR(0,0) << ", "<<__TENSOR(0,1) << ", "<<__TENSOR(0,2) << "} "
 407 //          << "{ "<<__TENSOR(1,0) << ", "<<__TENSOR(1,1) << ", "<<__TENSOR(1,2) << "} "
 408 //          << "{ "<<__TENSOR(2,0) << ", "<<__TENSOR(2,1) << ", "<<__TENSOR(2,2) << "}} ");
 409 //     for ( int i = 0; i < 3 && ok; ++i ) // loop on 3 eigenvalues
 410 //       {
 411 //         // [tensor-eig*I]
 412 //         double T[3][3]=
 413 //           {{ __TENSOR(0,0)-eig[i],__TENSOR(0,1),       __TENSOR(0,2),      },
 414 //            { __TENSOR(0,1),       __TENSOR(1,1)-eig[i],__TENSOR(1,2),      },
 415 //            { __TENSOR(0,2),       __TENSOR(1,2),       __TENSOR(2,2)-eig[i]}};
 416 //         // The determinant of T is zero, so that the equations are not linearly independent.
 417 //         // Therefore, we assign an arbitrary value (1.) to i-th component of eigenvector
 418 //         // and use two of the equations to compute the other two components
 419 //         double M[2][3], sol[2];
 420 //         for ( int j = 0, c = 0; j < 3; ++j )
 421 //           if ( i == j )
 422 //             M[0][2] = -T[0][j], M[1][2] = -T[1][j];
 423 //           else
 424 //             M[0][c] = T[0][j], M[1][c] = T[1][j], c++;
 425
 426 //         ok = solveSystemOfEquations<2>( M, sol );
 427
 428 //         double* eigenVec = __AXIS(i);
 429 //         for ( int j = 0, c = 0; j < 3; ++j )
 430 //           eigenVec[j] = ( i == j ) ? 1. : sol[c++];
 431
 432 //         // normilize
 433 //         double size = sqrt(eigenVec[0]*eigenVec[0] +
 434 //                            eigenVec[1]*eigenVec[1] +
 435 //                            eigenVec[2]*eigenVec[2] );
 436 //         if ((ok = (size > numeric_limits<double>::min() )))
 437 //           {
 438 //             eigenVec[0] /= size;
 439 //             eigenVec[1] /= size;
 440 //             eigenVec[2] /= size;
 441 //           }
 442 //       }
 443 //     if ( !ok )
 444 //       {
 445 //         __DMP( " solve3EquationSystem() - KO " );
 446 //         _axes = vector<double>( _dim*_dim, 0);
 447 //         __AXIS(0)[0] = __AXIS(1)[1] = __AXIS(2)[2] = 1.;
 448 //       }
 449 //     __DMP( " Eigen " << eig[0] << ", " << eig[1] << ", " << eig[2] );
 450 //     for ( int i=0; i <3; ++i )
 451 //       __DMP( i << ": " << __AXIS(i)[0] << ", " << __AXIS(i)[1] << ", " << __AXIS(i)[2] );
 452
 453 //     double* a0 = __AXIS(0), *a1 = __AXIS(1);
 454 //     double cross[3] = { a0[1]*a1[2]-a1[1]*a0[2],
 455 //                         a0[2]*a1[0]-a1[2]*a0[0],
 456 //                         a0[0]*a1[1]-a1[0]*a0[1] };
 457 //     __DMP( " Cross a1^a2 " << cross[0] << ", " << cross[1] << ", " << cross[2] );
 458 //   }
 459
 460   //================================================================================
 461   /*!
 462    * \brief Compute eigenvectors of inertia tensor
 463    */
 464   //================================================================================
 465
 466   // void DirectedBoundingBox::computeAxes2D(const std::vector<double>& tensor)
 467 //   {
 468 //     // compute principal moments of inertia which are eigenvalues of the tensor
 469 //     // by solving square equation det(tensor-eig*I)
 470 //     double X = (__TENSOR(0,0)+__TENSOR(1,1))/2;
 471 //     double Y = sqrt(4*__TENSOR(0,1)*__TENSOR(0,1) +
 472 //                     (__TENSOR(0,0)-__TENSOR(1,1)) * (__TENSOR(0,0)-__TENSOR(1,1)))/2;
 473 //     double eig[2] =
 474 //       {
 475 //         X + Y,
 476 //         X - Y
 477 //       };
 478 //     // compute eigenvector of the tensor at each eigenvalue
 479 //     // by solving system [tensor-eig*I]*[axis] = 0
 480 //     bool ok = true;
 481 //     for ( int i = 0; i < 2 && ok; ++i )
 482 //       {
 483 //         // [tensor-eig*I]
 484 //         double T[2][2]=
 485 //           {{ __TENSOR(0,0)-eig[i],__TENSOR(0,1)        },
 486 //            { __TENSOR(0,1),       __TENSOR(1,1)-eig[i] }};
 487
 488 //         // The determinant of T is zero, so that the equations are not linearly independent.
 489 //         // Therefore, we assign an arbitrary value (1.) to i-th component of eigenvector
 490 //         // and use one equation to compute the other component
 491 //         double* eigenVec = __AXIS(i);
 492 //         eigenVec[i] = 1.;
 493 //         int j = 1-i;
 494 //         if ((ok = ( fabs( T[j][j] ) > numeric_limits<double>::min() )))
 495 //           eigenVec[j] = -T[j][i] / T[j][j];
 496 //       }
 497 //     if ( !ok )
 498 //       {
 499 //         _axes = vector<double>( _dim*_dim, 0);
 500 //         __AXIS(0)[0] = __AXIS(1)[1] = 1.;
 501 //       }
 502 //   }
 503
 504   //================================================================================
 505   /*!
 506    * \brief Convert point coordinates into local coordinate system of the box
 507    */
 508   //================================================================================
 509
 510   void DirectedBoundingBox::toLocalCS(const double* p, double* pLoc) const
 511   {
 512     switch ( _dim )
 513       {
 514       case 3:
 515         pLoc[0] = dotprod<3>( p, __AXIS(0));
 516         pLoc[1] = dotprod<3>( p, __AXIS(1));
 517         pLoc[2] = dotprod<3>( p, __AXIS(2));
 518         break;
 519       case 2:
 520         pLoc[0] = dotprod<2>( p, __AXIS(0));
 521         pLoc[1] = dotprod<2>( p, __AXIS(1));
 522         break;
 523       default:
 524         pLoc[0] = p[0];
 525       }
 526   }
 527
 528   //================================================================================
 529   /*!
 530    * \brief Convert point coordinates from local coordinate system of the box to global CS
 531    */
 532   //================================================================================
 533
 534   void DirectedBoundingBox::fromLocalCS(const double* p, double* pGlob) const
 535   {
 536     switch ( _dim )
 537       {
 538       case 3:
 539         pGlob[0] = p[0] * __AXIS(0)[0] + p[1] * __AXIS(1)[0] + p[2] * __AXIS(2)[0];
 540         pGlob[1] = p[0] * __AXIS(0)[1] + p[1] * __AXIS(1)[1] + p[2] * __AXIS(2)[1];
 541         pGlob[2] = p[0] * __AXIS(0)[2] + p[1] * __AXIS(1)[2] + p[2] * __AXIS(2)[2];
 542         break;
 543       case 2:
 544         pGlob[0] = p[0] * __AXIS(0)[0] + p[1] * __AXIS(1)[0];
 545         pGlob[1] = p[0] * __AXIS(0)[1] + p[1] * __AXIS(1)[1];
 546         break;
 547       default:
 548         pGlob[0] = p[0];
 549       }
 550   }
 551
 552   //================================================================================
 553   /*!
 554    * \brief Enlarge box size by given value
 555    */
 556   //================================================================================
 557
 558   void DirectedBoundingBox::enlarge(const double tol)
 559   {
 560     for ( unsigned i = 0; i < _dim; ++i )
 561       __MIN(i) -= tol, __MAX(i) += tol;
 562   }
 563
 564   //================================================================================
 565   /*!
 566    * \brief Return coordinates of corners of bounding box
 567    */
 568   //================================================================================
 569
 570   void DirectedBoundingBox::getCorners(std::vector<double>& corners,
 571                                        const double*        minmax) const
 572   {
 573     int iC, nbCorners = 1;
 574     for ( int i=0;i<(int)_dim;++i ) nbCorners *= 2;
 575     corners.resize( nbCorners * _dim );
 576     // each coordinate is filled with either min or max, nbSwap is number of corners
 577     // after which min and max swap
 578     int nbSwap = nbCorners/2;
 579     for ( unsigned i = 0; i < _dim; ++i )
 580       {
 581         iC = 0;
 582         while ( iC < nbCorners )
 583           {
 584             for (int j = 0; j < nbSwap; ++j, ++iC ) corners[iC*_dim+i] = minmax[i*2];
 585             for (int j = 0; j < nbSwap; ++j, ++iC ) corners[iC*_dim+i] = minmax[i*2+1];
 586           }
 587         nbSwap /= 2;
 588       }
 589   }
 590
 591   //================================================================================
 592   /*!
 593    * \brief Test if this box intersects with the other
 594    *  \retval bool - true if there is no intersection
 595    */
 596   //================================================================================
 597
 598   bool DirectedBoundingBox::isDisjointWith(const DirectedBoundingBox& box) const
 599   {
 600     if ( _dim < 1 || box._dim < 1 ) return false; // empty box includes all
 601     if ( _dim == 1 )
 602       return isMinMaxOut( &box._minmax[0], &this->_minmax[0], _dim );
 603
 604     // boxes are disjoined if their minmaxes in local CS of either of boxes do not intersect
 605     for ( int isThisCS = 0; isThisCS < 2; ++isThisCS )
 606       {
 607         const DirectedBoundingBox* axisBox   = isThisCS ? this : &box;
 608         const DirectedBoundingBox* cornerBox = isThisCS ? &box : this;
 609
 610         // find minmax of cornerBox in the CS of axisBox
 611
 612         DirectedBoundingBox mmBox((double*)0,0,_dim); //!< empty box with CS == axisBox->_axes
 613         mmBox._axes = axisBox->_axes;
 614
 615         vector<double> corners;
 616         getCorners( corners, &cornerBox->_minmax[0] );
 617
 618         double globCorner[3];
 619         for ( std::size_t iC = 0, nC = corners.size()/_dim; iC < nC; ++iC)
 620           {
 621             cornerBox->fromLocalCS( &corners[iC*_dim], globCorner );
 622             mmBox.addPointToBox( globCorner );
 623           }
 624         if ( isMinMaxOut( &mmBox._minmax[0], &axisBox->_minmax[0], _dim ))
 625           return true;
 626       }
 627     return false;
 628   }
 629
 630   //================================================================================
 631   /*!
 632    * \brief Test if this box intersects with an non-directed box
 633    *  \retval bool - true if there is no intersection
 634    */
 635   //================================================================================
 636
 637   bool DirectedBoundingBox::isDisjointWith(const double* box) const
 638   {
 639     if ( _dim < 1 ) return false; // empty box includes all
 640     if ( _dim == 1 )
 641       return isMinMaxOut( &_minmax[0], box, _dim );
 642
 643     // boxes are disjoined if their minmaxes in local CS of either of boxes do not intersect
 644
 645     // compare minmaxes in locals CS of this directed box
 646     {
 647       vector<double> cornersOther;
 648       getCorners( cornersOther, box );
 649       DirectedBoundingBox mmBox((double*)0,0,_dim); //!< empty box with CS == this->_axes
 650       mmBox._axes = this->_axes;
 651       for ( std::size_t iC = 0, nC = cornersOther.size()/_dim; iC < nC; ++iC)
 652         mmBox.addPointToBox( &cornersOther[iC*_dim] );
 653
 654       if ( isMinMaxOut( &mmBox._minmax[0], &this->_minmax[0], _dim ))
 655         return true;
 656     }
 657
 658     // compare minmaxes in global CS
 659     {
 660       vector<double> cornersThis;
 661       getCorners( cornersThis, &_minmax[0] );
 662       DirectedBoundingBox mmBox((double*)0,0,_dim); //!< initailized _minmax
 663       double globCorner[3];
 664       for ( std::size_t iC = 0, nC = cornersThis.size()/_dim; iC < nC; ++iC)
 665         {
 666           fromLocalCS( &cornersThis[iC*_dim], globCorner );
 667           for ( int i = 0; i < (int)_dim; ++i )
 668             {
 669               if ( globCorner[i] < mmBox._minmax[i*2] )   mmBox._minmax[i*2] = globCorner[i];
 670               if ( globCorner[i] > mmBox._minmax[i*2+1] ) mmBox._minmax[i*2+1] = globCorner[i];
 671             }
 672         }
 673       if ( isMinMaxOut( &mmBox._minmax[0], box, _dim ))
 674         return true;
 675     }
 676     return false;
 677   }
 678
 679   //================================================================================
 680   /*!
 681    * \brief Return true if given point is out of this box
 682    */
 683   //================================================================================
 684
 685   bool DirectedBoundingBox::isOut(const double* point) const
 686   {
 687     if ( _dim < 1 ) return false; // empty box includes all
 688
 689     double pLoc[3];
 690     toLocalCS( point, pLoc );
 691     bool out = isLocalOut( pLoc );
 692 #ifdef _DEBUG_
 693     switch (_dim)
 694       {
 695       case 3:
 696         __DMP(__MYID<<": "<<point[0]<<", "<<point[1]<<", "<<point[2]<<" "<<(out?"OUT":"IN"));break;
 697       case 2:
 698         __DMP(__MYID<<": "<<point[0]<<", "<<point[1]<<" "<<(out?"OUT":"IN"));break;
 699       case 1:
 700         __DMP(__MYID<<": "<<point[0]<<" "<<(out?"OUT":"IN"));break;
 701       }
 702 #endif
 703     return out;
 704   }
 705
 706   //================================================================================
 707   /*!
 708    * \brief Return array of internal data
 709    */
 710   //================================================================================
 711
 712   vector<double> DirectedBoundingBox::getData() const
 713   {
 714     vector<double> data(1, _dim);
 715     if ( _dim > 0 )
 716     {
 717       data.insert( data.end(), &_axes[0], &_axes[0] + _axes.size());
 718       data.insert( data.end(), &_minmax[0], &_minmax[0] + _minmax.size());
 719     }
 720     if ( data.size() < (unsigned)dataSize( _dim ))
 721       data.resize( dataSize( _dim ), 0 );
 722     return data;
 723   }
 724
 725   //================================================================================
 726   /*!
 727    * \brief Initializes self with data retrieved via getData()
 728    */
 729   //================================================================================
 730
 731   void DirectedBoundingBox::setData(const double* data)
 732   {
 733     _dim = unsigned( *data++ );
 734     if ( _dim > 0 )
 735       {
 736         _axes.assign( data, data+_dim*_dim ); data += _dim*_dim;
 737         _minmax.assign( data, data+2*_dim );
 738       }
 739     else
 740       {
 741         _axes.clear();
 742         _minmax.clear();
 743       }
 744   }
 745
 746   //================================================================================
 747   /*!
 748    * \brief Return size of internal data returned by getData() depending on space dim
 749    */
 750   //================================================================================
 751
 752   int DirectedBoundingBox::dataSize(int dim)
 753   {
 754     return 1 + dim*dim + 2*dim; // : _dim + _axes + _minmax
 755   }
 756 }