src/INTERP_KERNEL/DirectedBoundingBox.cxx

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