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         *pIx = med_int(cx);
 146         if (*pIx >= mSize[0]) *pIx = mSize[0] - 1;
 147         else if (*pIx < 0) *pIx = 0;
 148
 149         *pIy = med_int(cy);
 150         if (*pIy >= mSize[1]) *pIy = mSize[1] - 1;
 151         else if (*pIy < 0) *pIy = 0;
 152
 153         *pIz = med_int(cz);
 154         if (*pIz >= mSize[2]) *pIz = mSize[2] - 1;
 155         else if (*pIz < 0) *pIz = 0;
 156 }
 157
 158
 159 int DecimationAccelGrid::getCellIndex(int pI, int pJ, int pK) const
 160 {
 161         int index = pK * (mSize[0] * mSize[1]) + pJ * mSize[0] + pI;
 162
 163         return index;
 164 }
 165
 166
 167 vector<PointOfField>& DecimationAccelGrid::getCell(const PointOfField& pPt)
 168 {
 169         int ix, iy, iz;
 170
 171         getCellCoord(
 172                 pPt.mXYZ[0], pPt.mXYZ[1], pPt.mXYZ[2],
 173                 &ix, &iy, &iz);
 174
 175         int index = getCellIndex(ix, iy, iz);
 176
 177         return mGrid[index];
 178 }
 179
 180
 181 vector<PointOfField> DecimationAccelGrid::findNeighbours(
 182         med_float pCenterX,
 183         med_float pCenterY,
 184         med_float pCenterZ,
 185         med_float pRadius) const
 186 {
 187         //---------------------------------------------------------------------
 188         // Determine the axis aligned bounding box of the sphere ((x, y, z), r)
 189         //---------------------------------------------------------------------
 190         med_float sphereBBoxMin[3];
 191         med_float sphereBBoxMax[3];
 192
 193         sphereBBoxMin[0] = pCenterX - pRadius;
 194         sphereBBoxMin[1] = pCenterY - pRadius;
 195         sphereBBoxMin[2] = pCenterZ - pRadius;
 196
 197         sphereBBoxMax[0] = pCenterX + pRadius;
 198         sphereBBoxMax[1] = pCenterY + pRadius;
 199         sphereBBoxMax[2] = pCenterZ + pRadius;
 200
 201         //---------------------------------------------------------------------
 202         // Determine the cells of the grid intersected by the sphere ((x, y, z), r)
 203         // => range of cells are [iMinCell[0], iMaxCell[0]] x [iMinCell[1], iMaxCell[1]] x [iMinCell[2], iMaxCell[2]]
 204         //---------------------------------------------------------------------
 205         int iMinCell[3];
 206         int iMaxCell[3];
 207
 208         getCellCoord(sphereBBoxMin[0], sphereBBoxMin[1], sphereBBoxMin[2],
 209                 &iMinCell[0], &iMinCell[1], &iMinCell[2]);
 210
 211         getCellCoord(sphereBBoxMax[0], sphereBBoxMax[1], sphereBBoxMax[2],
 212                 &iMaxCell[0], &iMaxCell[1], &iMaxCell[2]);
 213
 214         //---------------------------------------------------------------------
 215         // Collect points of the field which are in the sphere
 216         //---------------------------------------------------------------------
 217         vector<PointOfField> res;
 218
 219 /*
 220 // debug
 221 cout << "Center = " << pCenterX << " " << pCenterY << " " << pCenterZ << endl;
 222 cout << "Radius = " << pRadius << endl;
 223 cout << "Visited cells : [" << iMinCell[0] << " ; " << iMaxCell[0] << "] x ["<< iMinCell[1] << " ; " << iMaxCell[1] << "] x [" << iMinCell[2] << " ; " << iMaxCell[2] << "]" << endl;
 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 //printf("DEBUG: visited cell(%d %d %d) -> %d\n", i, j, k, idCell);
 236
 237                                 vector<PointOfField>& cell = mGrid[idCell];
 238
 239                                 // for all the points in the current cell
 240                                 for (vector<PointOfField>::const_iterator itPoint = cell.begin() ; itPoint != cell.end() ; itPoint++)
 241                                 {
 242                                         const PointOfField& currentPt = *itPoint;
 243
 244                                         // test if currentPt is in the sphere ((x, y, z), r)
 245                                         med_float vecX = currentPt.mXYZ[0] - pCenterX;
 246                                         med_float vecY = currentPt.mXYZ[1] - pCenterY;
 247                                         med_float vecZ = currentPt.mXYZ[2] - pCenterZ;
 248
 249                                         med_float norm = med_float( sqrt( vecX*vecX + vecY*vecY + vecZ*vecZ ) );
 250                                         if (norm < pRadius)
 251                                         {
 252                                                 // only add the point if it is different from (x, y, z)
 253                                                 if ((currentPt.mXYZ[0] != pCenterX) ||
 254                                                     (currentPt.mXYZ[1] != pCenterY) ||
 255                                                     (currentPt.mXYZ[2] != pCenterZ))
 256                                                 {
 257                                                         res.push_back(currentPt);
 258                                                 }
 259                                         }
 260                                 }
 261                         }
 262                 }
 263         }
 264
 265         return res;
 266 }
 267
 268
 269 void DecimationAccelGrid::computeBBox(const std::vector<PointOfField>& pPts)
 270 {
 271         for (int itDim = 0 ; itDim < 3 ; itDim++)
 272         {
 273                 mMin[itDim] = numeric_limits<med_float>::max();
 274                 mMax[itDim] = -mMin[itDim];
 275         }
 276
 277         for (unsigned i = 0 ; i < pPts.size() ; i++)
 278         {
 279                 for (int itDim = 0 ; itDim < 3 ; itDim++)
 280                 {
 281                         med_float coord = pPts[i].mXYZ[itDim];
 282                         if (coord < mMin[itDim]) mMin[itDim] = coord;
 283                         if (coord > mMax[itDim]) mMax[itDim] = coord;
 284                 }
 285         }
 286
 287         mInvLen[0] = med_float(mSize[0]) / (mMax[0] - mMin[0]);
 288         mInvLen[1] = med_float(mSize[1]) / (mMax[1] - mMin[1]);
 289         mInvLen[2] = med_float(mSize[2]) / (mMax[2] - mMin[2]);
 290 }
 291
 292
 293 ostream& operator<<(ostream& pOs, DecimationAccelGrid& pG)
 294 {
 295         pOs << "DecimationAccelGrid:" << endl;
 296         pOs << "    Num=" << pG.mNum << endl;
 297         pOs << "    Size=" << pG.mSize[0] << " x " << pG.mSize[1] << " x " << pG.mSize[2] << endl;
 298         pOs << "    BBox=[" << pG.mMin[0] << " ; " << pG.mMax[0] << "] x [" << pG.mMin[1] << " ; " << pG.mMax[1] << "] x [" << pG.mMin[2] << " ; " << pG.mMax[2] << "]" << endl;
 299         pOs << "    Inv len.=" << pG.mInvLen[0] << " ; " << pG.mInvLen[1] << " ; " << pG.mInvLen[2] << endl;
 300
 301         if (pG.mFlagPrintAll)
 302         {
 303                 int checkNumCells = 0;
 304                 int numCells = pG.mSize[0] * pG.mSize[1] * pG.mSize[2];
 305                 for (int i = 0 ; i < numCells ; i++)
 306                 {
 307                         vector<PointOfField>& cell = pG.mGrid[i];
 308                         cout << "    Cell " << i << ": #=" << cell.size() << ": " << endl;
 309                         for (unsigned j = 0 ; j < cell.size() ; j++)
 310                         {
 311                                 cout << "        " << cell[j] << endl;
 312                         }
 313                         checkNumCells += cell.size();
 314                 }
 315
 316                 if (pG.mNum != checkNumCells) throw IllegalStateException("", __FILE__, __LINE__);
 317         }
 318
 319         return pOs;
 320 }
 321
 322
 323 } // namespace multipr
 324
 325 // EOF