1 // Project MULTIPR, IOLS WP1.2.1 - EDF/CS
2 // Partitioning/decimation module for the SALOME v3.2 platform
5 * \file MULTIPR_DecimationAccel.cxx
7 * \brief see MULTIPR_DecimationAccel.hxx
9 * \author Olivier LE ROUX - CS, Virtual Reality Dpt
14 //*****************************************************************************
16 //*****************************************************************************
18 #include "MULTIPR_DecimationAccel.hxx"
19 #include "MULTIPR_PointOfField.hxx"
20 #include "MULTIPR_Exceptions.hxx"
31 //*****************************************************************************
32 // Class DecimationAccel implementation
33 //*****************************************************************************
36 ostream& operator<<(ostream& pOs, DecimationAccel& pA)
38 pOs << "DecimationAccel:" << endl;
43 //*****************************************************************************
44 // Class DecimationAccelGrid
45 //*****************************************************************************
47 DecimationAccelGrid::DecimationAccelGrid()
55 DecimationAccelGrid::~DecimationAccelGrid()
61 void DecimationAccelGrid::reset()
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();
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();
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();
87 mFlagPrintAll = false;
91 void DecimationAccelGrid::create(const std::vector<PointOfField>& pPts)
93 // check if grid have been initialized
94 if (mSize[0] == 0) throw IllegalStateException("call setSize() before", __FILE__, __LINE__);
96 // compute bbox of the grid
100 int size = mSize[0] * mSize[1] * mSize[2];
101 mGrid = new vector<PointOfField>[size];
105 for (int i = 0 ; i < mNum ; i++)
107 vector<PointOfField>& cell = getCell(pPts[i]);
108 cell.push_back(pPts[i]);
113 void DecimationAccelGrid::configure(const char* pArgv)
116 if (pArgv == NULL) throw NullArgumentException("", __FILE__, __LINE__);
122 int ret = sscanf(pArgv, "%d %d %d", &sizeX, &sizeY, &sizeZ);
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__);
137 void DecimationAccelGrid::getCellCoord(
138 med_float pX, med_float pY, med_float pZ,
139 int* pIx, int* pIy, int* pIz) const
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];
146 if (*pIx >= mSize[0]) *pIx = mSize[0] - 1;
147 else if (*pIx < 0) *pIx = 0;
150 if (*pIy >= mSize[1]) *pIy = mSize[1] - 1;
151 else if (*pIy < 0) *pIy = 0;
154 if (*pIz >= mSize[2]) *pIz = mSize[2] - 1;
155 else if (*pIz < 0) *pIz = 0;
159 int DecimationAccelGrid::getCellIndex(int pI, int pJ, int pK) const
161 int index = pK * (mSize[0] * mSize[1]) + pJ * mSize[0] + pI;
167 vector<PointOfField>& DecimationAccelGrid::getCell(const PointOfField& pPt)
172 pPt.mXYZ[0], pPt.mXYZ[1], pPt.mXYZ[2],
175 int index = getCellIndex(ix, iy, iz);
181 vector<PointOfField> DecimationAccelGrid::findNeighbours(
185 med_float pRadius) const
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];
193 sphereBBoxMin[0] = pCenterX - pRadius;
194 sphereBBoxMin[1] = pCenterY - pRadius;
195 sphereBBoxMin[2] = pCenterZ - pRadius;
197 sphereBBoxMax[0] = pCenterX + pRadius;
198 sphereBBoxMax[1] = pCenterY + pRadius;
199 sphereBBoxMax[2] = pCenterZ + pRadius;
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 //---------------------------------------------------------------------
208 getCellCoord(sphereBBoxMin[0], sphereBBoxMin[1], sphereBBoxMin[2],
209 &iMinCell[0], &iMinCell[1], &iMinCell[2]);
211 getCellCoord(sphereBBoxMax[0], sphereBBoxMax[1], sphereBBoxMax[2],
212 &iMaxCell[0], &iMaxCell[1], &iMaxCell[2]);
214 //---------------------------------------------------------------------
215 // Collect points of the field which are in the sphere
216 //---------------------------------------------------------------------
217 vector<PointOfField> res;
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;
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++)
229 for (int j = iMinCell[1] ; j <= iMaxCell[1] ; j++)
231 for (int k = iMinCell[2] ; k <= iMaxCell[2] ; k++)
233 int idCell = getCellIndex(i, j, k);
235 //printf("DEBUG: visited cell(%d %d %d) -> %d\n", i, j, k, idCell);
237 vector<PointOfField>& cell = mGrid[idCell];
239 // for all the points in the current cell
240 for (vector<PointOfField>::const_iterator itPoint = cell.begin() ; itPoint != cell.end() ; itPoint++)
242 const PointOfField& currentPt = *itPoint;
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;
249 med_float norm = med_float( sqrt( vecX*vecX + vecY*vecY + vecZ*vecZ ) );
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))
257 res.push_back(currentPt);
269 void DecimationAccelGrid::computeBBox(const std::vector<PointOfField>& pPts)
271 for (int itDim = 0 ; itDim < 3 ; itDim++)
273 mMin[itDim] = numeric_limits<med_float>::max();
274 mMax[itDim] = -mMin[itDim];
277 for (unsigned i = 0 ; i < pPts.size() ; i++)
279 for (int itDim = 0 ; itDim < 3 ; itDim++)
281 med_float coord = pPts[i].mXYZ[itDim];
282 if (coord < mMin[itDim]) mMin[itDim] = coord;
283 if (coord > mMax[itDim]) mMax[itDim] = coord;
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]);
293 ostream& operator<<(ostream& pOs, DecimationAccelGrid& pG)
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;
301 if (pG.mFlagPrintAll)
303 int checkNumCells = 0;
304 int numCells = pG.mSize[0] * pG.mSize[1] * pG.mSize[2];
305 for (int i = 0 ; i < numCells ; i++)
307 vector<PointOfField>& cell = pG.mGrid[i];
308 cout << " Cell " << i << ": #=" << cell.size() << ": " << endl;
309 for (unsigned j = 0 ; j < cell.size() ; j++)
311 cout << " " << cell[j] << endl;
313 checkNumCells += cell.size();
316 if (pG.mNum != checkNumCells) throw IllegalStateException("", __FILE__, __LINE__);
323 } // namespace multipr