src/MULTIPR/MULTIPR_DecimationAccel.cxx

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