src/MULTIPR/MULTIPR_DecimationAccel.cxx

   1 //  Copyright (C) 2007-2008  CEA/DEN, EDF R&D
   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 // Partitioning/decimation module for the SALOME v3.2 platform
  20 //
  21 /**
  22  * \file    MULTIPR_DecimationAccel.cxx
  23  *
  24  * \brief   see MULTIPR_DecimationAccel.hxx
  25  *
  26  * \author  Olivier LE ROUX - CS, Virtual Reality Dpt
  27  *
  28  * \date    01/2007
  29  */
  30
  31 //*****************************************************************************
  32 // Includes section
  33 //*****************************************************************************
  34
  35 #include "MULTIPR_DecimationAccel.hxx"
  36 #include "MULTIPR_PointOfField.hxx"
  37 #include "MULTIPR_Exceptions.hxx"
  38
  39 #include <cmath>
  40 #include <iostream>
  41
  42 using namespace std;
  43
  44
  45 namespace multipr
  46 {
  47
  48
  49 //*****************************************************************************
  50 // Class DecimationAccel implementation
  51 //*****************************************************************************
  52
  53
  54 ostream& operator<<(ostream& pOs, DecimationAccel& pA)
  55 {
  56     pOs << "DecimationAccel:" << endl;
  57     return pOs;
  58 }
  59
  60
  61 //*****************************************************************************
  62 // Class DecimationAccelGrid
  63 //*****************************************************************************
  64
  65 DecimationAccelGrid::DecimationAccelGrid()
  66 {
  67     mGrid = NULL;
  68
  69     reset();
  70 }
  71
  72
  73 DecimationAccelGrid::~DecimationAccelGrid()
  74 {
  75     reset();
  76 }
  77
  78
  79 void DecimationAccelGrid::reset()
  80 {
  81     mNum = 0;
  82
  83     mMin[0] = numeric_limits<med_float>::quiet_NaN();
  84     mMin[1] = numeric_limits<med_float>::quiet_NaN();
  85     mMin[2] = numeric_limits<med_float>::quiet_NaN();
  86
  87     mMax[0] = numeric_limits<med_float>::quiet_NaN();
  88     mMax[1] = numeric_limits<med_float>::quiet_NaN();
  89     mMax[2] = numeric_limits<med_float>::quiet_NaN();
  90
  91     mInvLen[0] = numeric_limits<med_float>::quiet_NaN();
  92     mInvLen[1] = numeric_limits<med_float>::quiet_NaN();
  93     mInvLen[2] = numeric_limits<med_float>::quiet_NaN();
  94
  95     mSize[0] = 0;
  96     mSize[1] = 0;
  97     mSize[2] = 0;
  98
  99     if (mGrid != NULL)
 100     {
 101         delete[] mGrid;
 102         mGrid = NULL;
 103     }
 104
 105     mFlagPrintAll = false;
 106 }
 107
 108
 109 void DecimationAccelGrid::create(const std::vector<PointOfField>& pPts)
 110 {
 111     // check if grid have been initialized
 112     if (mSize[0] == 0) throw IllegalStateException("call setSize() before", __FILE__, __LINE__);
 113
 114     // compute bbox of the grid
 115     computeBBox(pPts);
 116
 117     // allocate the grid
 118     int size = mSize[0] * mSize[1] * mSize[2];
 119     mGrid = new vector<PointOfField>[size];
 120
 121     // fill the grid
 122     mNum = pPts.size();
 123     for (int i = 0 ; i < mNum ; i++)
 124     {
 125         vector<PointOfField>& cell = getCell(pPts[i]);
 126         cell.push_back(pPts[i]);
 127     }
 128 }
 129
 130
 131 void DecimationAccelGrid::configure(const char* pArgv)
 132 {
 133     // check arguments
 134     if (pArgv == NULL) throw NullArgumentException("", __FILE__, __LINE__);
 135
 136     int sizeX = 0;
 137     int sizeY = 0;
 138     int sizeZ = 0;
 139
 140     int ret = sscanf(pArgv, "%d %d %d", &sizeX, &sizeY, &sizeZ);
 141
 142     if (ret != 3) throw IllegalArgumentException("expected 3 parameters", __FILE__, __LINE__);
 143     if (sizeX <= 0) throw IllegalArgumentException("number of cells along X-axis must be > 0", __FILE__, __LINE__);
 144     if (sizeY <= 0) throw IllegalArgumentException("number of cells along Y-axis must be > 0", __FILE__, __LINE__);
 145     if (sizeZ <= 0) throw IllegalArgumentException("number of cells along Z-axis must be > 0", __FILE__, __LINE__);
 146
 147     reset();
 148
 149     mSize[0] = sizeX;
 150     mSize[1] = sizeY;
 151     mSize[2] = sizeZ;
 152 }
 153
 154
 155 void DecimationAccelGrid::getCellCoord(
 156         med_float pX, med_float pY, med_float pZ,
 157         int* pIx, int* pIy, int* pIz) const
 158 {
 159     med_float cx = (pX - mMin[0]) * mInvLen[0];
 160     med_float cy = (pY - mMin[1]) * mInvLen[1];
 161     med_float cz = (pZ - mMin[2]) * mInvLen[2];
 162
 163     // clamp all indices to avoid overflow
 164     *pIx = med_int(cx);
 165     if (*pIx >= mSize[0]) *pIx = mSize[0] - 1;
 166     else if (*pIx < 0) *pIx = 0;
 167
 168     *pIy = med_int(cy);
 169     if (*pIy >= mSize[1]) *pIy = mSize[1] - 1;
 170     else if (*pIy < 0) *pIy = 0;
 171
 172     *pIz = med_int(cz);
 173     if (*pIz >= mSize[2]) *pIz = mSize[2] - 1;
 174     else if (*pIz < 0) *pIz = 0;
 175 }
 176
 177
 178 int DecimationAccelGrid::getCellIndex(int pI, int pJ, int pK) const
 179 {
 180     int index = pK * (mSize[0] * mSize[1]) + pJ * mSize[0] + pI;
 181
 182     return index;
 183 }
 184
 185
 186 vector<PointOfField>& DecimationAccelGrid::getCell(const PointOfField& pPt)
 187 {
 188     int ix, iy, iz;
 189
 190     getCellCoord(
 191         pPt.mXYZ[0], pPt.mXYZ[1], pPt.mXYZ[2],
 192         &ix, &iy, &iz);
 193
 194     int index = getCellIndex(ix, iy, iz);
 195
 196     return mGrid[index];
 197 }
 198
 199
 200 vector<PointOfField> DecimationAccelGrid::findNeighbours(
 201     med_float pCenterX,
 202     med_float pCenterY,
 203     med_float pCenterZ,
 204     med_float pRadius) const
 205 {
 206     //---------------------------------------------------------------------
 207     // Determine the axis aligned bounding box of the sphere ((x, y, z), r)
 208     //---------------------------------------------------------------------
 209     med_float sphereBBoxMin[3];
 210     med_float sphereBBoxMax[3];
 211
 212     sphereBBoxMin[0] = pCenterX - pRadius;
 213     sphereBBoxMin[1] = pCenterY - pRadius;
 214     sphereBBoxMin[2] = pCenterZ - pRadius;
 215
 216     sphereBBoxMax[0] = pCenterX + pRadius;
 217     sphereBBoxMax[1] = pCenterY + pRadius;
 218     sphereBBoxMax[2] = pCenterZ + pRadius;
 219
 220     //---------------------------------------------------------------------
 221     // Determine the cells of the grid intersected by the sphere ((x, y, z), r)
 222     // => range of cells are [iMinCell[0], iMaxCell[0]] x [iMinCell[1], iMaxCell[1]] x [iMinCell[2], iMaxCell[2]]
 223     //---------------------------------------------------------------------
 224     int iMinCell[3];
 225     int iMaxCell[3];
 226
 227     getCellCoord(sphereBBoxMin[0], sphereBBoxMin[1], sphereBBoxMin[2],
 228             &iMinCell[0], &iMinCell[1], &iMinCell[2]);
 229
 230     getCellCoord(sphereBBoxMax[0], sphereBBoxMax[1], sphereBBoxMax[2],
 231             &iMaxCell[0], &iMaxCell[1], &iMaxCell[2]);
 232
 233     //---------------------------------------------------------------------
 234     // Collect points of the field which are in the sphere
 235     //---------------------------------------------------------------------
 236     vector<PointOfField> res;
 237
 238     if (isnan(pCenterX) || isnan(pCenterY) || isnan(pCenterZ))
 239     {
 240         return res;
 241     }
 242
 243     // for all the cells in the grid intersected by the sphere ((x, y, z), r)
 244     for (int i = iMinCell[0] ; i <= iMaxCell[0] ; i++)
 245     {
 246         for (int j = iMinCell[1] ; j <= iMaxCell[1] ; j++)
 247         {
 248             for (int k = iMinCell[2] ; k <= iMaxCell[2] ; k++)
 249             {
 250                 int idCell = getCellIndex(i, j, k);
 251
 252                 vector<PointOfField>& cell = mGrid[idCell];
 253
 254                 // for all the points in the current cell
 255                 for (vector<PointOfField>::const_iterator itPoint = cell.begin() ; itPoint != cell.end() ; itPoint++)
 256                 {
 257                     const PointOfField& currentPt = *itPoint;
 258
 259                     // test if currentPt is in the sphere ((x, y, z), r)
 260                     med_float vecX = currentPt.mXYZ[0] - pCenterX;
 261                     med_float vecY = currentPt.mXYZ[1] - pCenterY;
 262                     med_float vecZ = currentPt.mXYZ[2] - pCenterZ;
 263
 264                     med_float norm = med_float( sqrt( vecX*vecX + vecY*vecY + vecZ*vecZ ) );
 265                     if (norm < pRadius)
 266                     {
 267                         // only add the point if it is different from (x, y, z)
 268                         if ((currentPt.mXYZ[0] != pCenterX) ||
 269                             (currentPt.mXYZ[1] != pCenterY) ||
 270                             (currentPt.mXYZ[2] != pCenterZ))
 271                         {
 272                             res.push_back(currentPt);
 273                         }
 274                     }
 275                 }
 276             }
 277         }
 278     }
 279
 280     return res;
 281 }
 282
 283
 284 void DecimationAccelGrid::computeBBox(const std::vector<PointOfField>& pPts)
 285 {
 286     for (int itDim = 0 ; itDim < 3 ; itDim++)
 287     {
 288         mMin[itDim] = numeric_limits<med_float>::max();
 289         mMax[itDim] = -mMin[itDim];
 290     }
 291
 292     for (unsigned i = 0 ; i < pPts.size() ; i++)
 293     {
 294         for (int itDim = 0 ; itDim < 3 ; itDim++)
 295         {
 296             med_float coord = pPts[i].mXYZ[itDim];
 297             if (coord < mMin[itDim]) mMin[itDim] = coord;
 298             if (coord > mMax[itDim]) mMax[itDim] = coord;
 299         }
 300     }
 301
 302     mInvLen[0] = med_float(mSize[0]) / (mMax[0] - mMin[0]);
 303     mInvLen[1] = med_float(mSize[1]) / (mMax[1] - mMin[1]);
 304     mInvLen[2] = med_float(mSize[2]) / (mMax[2] - mMin[2]);
 305 }
 306
 307
 308 ostream& operator<<(ostream& pOs, DecimationAccelGrid& pG)
 309 {
 310     pOs << "DecimationAccelGrid:" << endl;
 311     pOs << "    Num=" << pG.mNum << endl;
 312     pOs << "    Size=" << pG.mSize[0] << " x " << pG.mSize[1] << " x " << pG.mSize[2] << endl;
 313     pOs << "    BBox=[" << pG.mMin[0] << " ; " << pG.mMax[0] << "] x [" << pG.mMin[1] << " ; " << pG.mMax[1] << "] x [" << pG.mMin[2] << " ; " << pG.mMax[2] << "]" << endl;
 314     pOs << "    Inv len.=" << pG.mInvLen[0] << " ; " << pG.mInvLen[1] << " ; " << pG.mInvLen[2] << endl;
 315
 316     if (pG.mFlagPrintAll)
 317     {
 318         int checkNumCells = 0;
 319         int numCells = pG.mSize[0] * pG.mSize[1] * pG.mSize[2];
 320         for (int i = 0 ; i < numCells ; i++)
 321         {
 322             vector<PointOfField>& cell = pG.mGrid[i];
 323             cout << "    Cell " << i << ": #=" << cell.size() << ": " << endl;
 324             for (unsigned j = 0 ; j < cell.size() ; j++)
 325             {
 326                 cout << "        " << cell[j] << endl;
 327             }
 328             checkNumCells += cell.size();
 329         }
 330
 331         if (pG.mNum != checkNumCells) throw IllegalStateException("", __FILE__, __LINE__);
 332     }
 333
 334     return pOs;
 335 }
 336
 337
 338 } // namespace multipr
 339
 340 // EOF