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