ESyS-Particle
4.0.1
|
00001 00002 // // 00003 // Copyright (c) 2003-2011 by The University of Queensland // 00004 // Earth Systems Science Computational Centre (ESSCC) // 00005 // http://www.uq.edu.au/esscc // 00006 // // 00007 // Primary Business: Brisbane, Queensland, Australia // 00008 // Licensed under the Open Software License version 3.0 // 00009 // http://www.opensource.org/licenses/osl-3.0.php // 00010 // // 00012 00013 00014 #include "Foundation/console.h" 00015 #include "Foundation/StringUtil.h" 00016 #include "Geometry/RandomSpherePacker.h" 00017 #include "Geometry/SphereFitter.h" 00018 #include "Geometry/ParticleComparer.h" 00019 00020 #include <algorithm> 00021 #include <stdexcept> 00022 #include <float.h> 00023 00024 namespace esys 00025 { 00026 namespace lsm 00027 { 00028 template<typename TmplTraits> 00029 SphereFittedPIterator<TmplTraits>::SphereFittedPIterator( 00030 Packer &packer, 00031 int maxInsertionFailures, 00032 const BoundingSphere &bSphere 00033 ) : 00034 m_pPacker(&packer), 00035 m_fitterPtrVector(), 00036 m_maxInsertionFailures(maxInsertionFailures), 00037 m_lastFailCount(0), 00038 m_successCount(0), 00039 m_next(Fitter::getInvalidParticle()), 00040 m_bSphere(bSphere) 00041 { 00042 if (m_next.isValid()) 00043 { 00044 throw 00045 std::runtime_error( 00046 std::string("isValid method returns true for INVALID particle:") 00047 + 00048 StringUtil::toString(m_next) 00049 ); 00050 } 00051 initialiseFitterPtrVector(); 00052 } 00053 00054 template<typename TmplTraits> 00055 int SphereFittedPIterator<TmplTraits>::getMaxInsertionFailures() const 00056 { 00057 return m_maxInsertionFailures; 00058 } 00059 00060 template<typename TmplTraits> 00061 const typename SphereFittedPIterator<TmplTraits>::Packer & 00062 SphereFittedPIterator<TmplTraits>::getPacker() const 00063 { 00064 return *m_pPacker; 00065 } 00066 00067 template<typename TmplTraits> 00068 typename SphereFittedPIterator<TmplTraits>::Packer & 00069 SphereFittedPIterator<TmplTraits>::getPacker() 00070 { 00071 return *m_pPacker; 00072 } 00073 00074 template<typename TmplTraits> 00075 const BoundingSphere & 00076 SphereFittedPIterator<TmplTraits>::getBSphere() const 00077 { 00078 return m_bSphere; 00079 } 00080 00081 template<typename TmplTraits> 00082 void 00083 SphereFittedPIterator<TmplTraits>::initialiseFitterPtrVector() 00084 { 00085 FitterPtrVector fitters; 00086 fitters.push_back(FitterPtr(new Move2SurfaceFitter(getPacker()))); 00087 if (getPacker().is2d()) 00088 { 00089 fitters.push_back(FitterPtr(new TwoDSFitter(getPacker()))); 00090 fitters.push_back(FitterPtr(new TwoDSSphereFitter(getPacker(), getBSphere()))); 00091 } 00092 else 00093 { 00094 fitters.push_back(FitterPtr(new ThreeDFitter(getPacker()))); 00095 fitters.push_back(FitterPtr(new ThreeDSSphereFitter(getPacker(), getBSphere()))); 00096 } 00097 m_fitterPtrVector = fitters; 00098 } 00099 00100 template<typename TmplTraits> 00101 typename SphereFittedPIterator<TmplTraits>::FitterPtrVector & 00102 SphereFittedPIterator<TmplTraits>::getFitterPtrVector() 00103 { 00104 return m_fitterPtrVector; 00105 } 00106 00107 template<typename TmplTraits> 00108 const typename SphereFittedPIterator<TmplTraits>::FitterPtrVector & 00109 SphereFittedPIterator<TmplTraits>::getFitterPtrVector() const 00110 { 00111 return m_fitterPtrVector; 00112 } 00113 00114 template<typename TmplTraits> 00115 Vec3 SphereFittedPIterator<TmplTraits>::getRandomPoint() const 00116 { 00117 return getPacker().getRandomPoint(); 00118 } 00119 00120 template<typename TmplTraits> 00121 typename SphereFittedPIterator<TmplTraits>::Particle 00122 SphereFittedPIterator<TmplTraits>::getCandidateParticle( 00123 const Vec3 &point 00124 ) 00125 { 00126 return getPacker().getCandidateParticle(point); 00127 } 00128 00129 template<typename TmplTraits> 00130 typename SphereFittedPIterator<TmplTraits>::ParticleVector 00131 SphereFittedPIterator<TmplTraits>::getClosestNeighbours( 00132 const Particle& particle, 00133 int numClosest 00134 ) 00135 { 00136 return getPacker().getClosestNeighbours(particle, numClosest); 00137 } 00138 00139 template<typename TmplTraits> 00140 typename SphereFittedPIterator<TmplTraits>::Particle & 00141 SphereFittedPIterator<TmplTraits>::generateNext() 00142 { 00143 m_next = Fitter::getInvalidParticle(); 00144 if (m_lastFailCount < getMaxInsertionFailures()) 00145 { 00146 int numFails = 0; 00147 //int insertCount = 0; 00148 FitterPtrVector fitters = getFitterPtrVector(); 00149 Particle particle = Fitter::getInvalidParticle(); 00150 Particle fitParticle = particle; 00151 Plane plane=Plane(Vec3(-1.0, 0, 0), Vec3(DBL_MAX/2, DBL_MAX/2, DBL_MAX/2)); 00152 ParticleVector neighbours; 00153 while ((!fitParticle.isValid()) && (numFails < getMaxInsertionFailures())) 00154 { 00155 particle = getCandidateParticle(getRandomPoint()); 00156 neighbours = getClosestNeighbours(particle, 4); 00157 00158 for ( 00159 typename FitterPtrVector::iterator it = getFitterPtrVector().begin(); 00160 ( 00161 (it != getFitterPtrVector().end()) 00162 && 00163 (!fitParticle.isValid()) 00164 ); 00165 it++ 00166 ) 00167 { 00168 fitParticle = (*it)->getFitParticle(particle, neighbours, plane); 00169 } 00170 numFails++; 00171 } 00172 m_lastFailCount = numFails; 00173 m_successCount += ((fitParticle.isValid()) ? 1 : 0); 00174 m_next = fitParticle; 00175 } 00176 return m_next; 00177 } 00178 00179 template<typename TmplTraits> 00180 bool SphereFittedPIterator<TmplTraits>::hasNext() 00181 { 00182 return (m_next.isValid() || generateNext().isValid()); 00183 } 00184 00185 template<typename TmplTraits> 00186 typename SphereFittedPIterator<TmplTraits>::Particle 00187 SphereFittedPIterator<TmplTraits>::next() 00188 { 00189 if (!(m_next.isValid())) 00190 { 00191 generateNext(); 00192 } 00193 const Particle next = m_next; 00194 m_next = Fitter::getInvalidParticle(); 00195 return next; 00196 } 00197 00198 template<typename TmplTraits> 00199 void SphereFittedPIterator<TmplTraits>::logInfo() 00200 { 00201 console.Info() 00202 << "BoundingSphere: minPt = " << getPacker().getBBox().getMinPt() 00203 << ", maxPt = " << getPacker().getBBox().getMaxPt() << "\n" 00204 << "BoundingSphere: sizes = " << getPacker().getBBox().getSizes() << "\n"; 00205 00206 for ( 00207 typename FitterPtrVector::iterator it = getFitterPtrVector().begin(); 00208 (it != getFitterPtrVector().end()); 00209 it++ 00210 ) { 00211 console.Info() << (*(*it)).toString() << "\n"; 00212 } 00213 console.Info() << "Total successful fits = " << m_successCount << "\n"; 00214 } 00215 00216 00217 00218 00219 00220 00221 00222 00223 00224 00225 00226 00227 00228 00229 00230 00231 00232 00233 00234 00235 00236 00237 00238 00239 00240 00241 00242 template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap> 00243 RandomSpherePacker<TPartGen,TCubicPackWrap>::RandomSpherePacker( 00244 ParticleGeneratorPtr particleGeneratorPtr, 00245 ParticlePoolPtr particlePoolPtr, 00246 NTablePtr nTablePtr, 00247 const BoundingSphere &bSphere, 00248 double tolerance, 00249 double cubicPackRadius, 00250 int maxInsertionFailures, 00251 bool do2d 00252 ) 00253 : Inherited( 00254 particleGeneratorPtr, 00255 particlePoolPtr, 00256 nTablePtr, 00257 do2d ? bSphere.get2dBBox() : bSphere.getBBox(), 00258 BoolVector(3,false), 00259 tolerance, 00260 cubicPackRadius 00261 ), 00262 m_bSphere(bSphere), 00263 m_maxInsertionFailures(maxInsertionFailures) 00264 { 00265 } 00266 00267 template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap> 00268 RandomSpherePacker<TPartGen,TCubicPackWrap>::~RandomSpherePacker() 00269 { 00270 } 00271 00272 template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap> 00273 const BoundingSphere & 00274 RandomSpherePacker<TPartGen,TCubicPackWrap>::getBSphere( 00275 ) const 00276 { 00277 return m_bSphere; 00278 } 00279 00280 template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap> 00281 double 00282 RandomSpherePacker<TPartGen,TCubicPackWrap>::getRandom( 00283 double minVal, 00284 double maxVal 00285 ) const 00286 { 00287 return minVal + (maxVal-minVal)*rng::s_zeroOneUniform(); 00288 } 00289 00290 template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap> 00291 Vec3 RandomSpherePacker<TPartGen,TCubicPackWrap>::getRandomPoint() const 00292 { 00293 return 00294 Vec3( 00295 getRandom(this->getBBox().getMinPt().X(), this->getBBox().getMaxPt().X()), 00296 getRandom(this->getBBox().getMinPt().Y(), this->getBBox().getMaxPt().Y()), 00297 getRandom(this->getBBox().getMinPt().Z(), this->getBBox().getMaxPt().Z()) 00298 ); 00299 } 00300 00301 template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap> 00302 typename RandomSpherePacker<TPartGen,TCubicPackWrap>::ParticleVector 00303 RandomSpherePacker<TPartGen,TCubicPackWrap>::getClosestNeighbours( 00304 const Particle& particle, 00305 int numClosest 00306 ) 00307 { 00308 ParticleVector neighbourVector = 00309 this->getNTable().getUniqueNeighbourVector( 00310 particle.getPos(), 00311 this->getParticleGenerator().getMaxFitRadius() + this->getTolerance() 00312 ); 00313 std::sort( 00314 neighbourVector.begin(), 00315 neighbourVector.end(), 00316 ParticleComparer<Particle>(particle) 00317 ); 00318 00319 if (static_cast<int>(neighbourVector.size()) > numClosest) { 00320 neighbourVector.erase( 00321 neighbourVector.begin() + numClosest, 00322 neighbourVector.end() 00323 ); 00324 } 00325 00326 return neighbourVector; 00327 } 00328 00329 template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap> 00330 bool RandomSpherePacker<TPartGen,TCubicPackWrap >::particleIsValid( 00331 const Particle &particle 00332 ) const 00333 { 00334 return 00335 ( 00336 this->getParticleGenerator().isValidFitRadius(particle.getRad()) 00337 && 00338 particleFitsInBSphereWithNeighbours(particle) 00339 ); 00340 } 00341 00342 template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap> 00343 int RandomSpherePacker<TPartGen,TCubicPackWrap>::getMaxInsertionFailures() const 00344 { 00345 return m_maxInsertionFailures; 00346 } 00347 00348 template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap> 00349 void RandomSpherePacker<TPartGen,TCubicPackWrap>::generateRandomFill() 00350 { 00351 StuffedParticleIterator it = 00352 StuffedParticleIterator( 00353 *this, 00354 getMaxInsertionFailures(), 00355 getBSphere() 00356 ); 00357 while (it.hasNext()) 00358 { 00359 const Particle p = it.next(); 00360 this->createAndInsertParticle(p); 00361 } 00362 it.logInfo(); 00363 } 00364 00365 template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap> 00366 bool 00367 RandomSpherePacker<TPartGen,TCubicPackWrap>::particleFitsInBSphere( 00368 const Particle &particle 00369 ) const 00370 { 00371 return 00372 getBSphere().contains( 00373 BoundingSphere(particle.getPos(), particle.getRadius()), 00374 this->getTolerance() 00375 ); 00376 } 00377 00378 template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap> 00379 bool 00380 RandomSpherePacker<TPartGen,TCubicPackWrap>::particleFitsInBSphereWithNeighbours( 00381 const Particle &particle 00382 ) const 00383 { 00384 return 00385 ( 00386 particleFitsInBSphere(particle) 00387 && 00388 this->particleFitsWithNeighbours(particle) 00389 ); 00390 } 00391 00392 template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap> 00393 void RandomSpherePacker<TPartGen,TCubicPackWrap>::generateCubicPackingInSphere() 00394 { 00395 GridIterator pointIt = GridIterator(this->getBBox(), this->getCubicPackingRadius()); 00396 while (pointIt.hasNext()) { 00397 const Particle candidate = 00398 this->getCandidateParticle(pointIt.next(), this->getCubicPackingRadius()); 00399 if (particleFitsInBSphereWithNeighbours(candidate)) { 00400 this->createAndInsertParticle(candidate); 00401 } 00402 } 00403 } 00404 00405 template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap> 00406 void RandomSpherePacker<TPartGen,TCubicPackWrap>::generate() 00407 { 00408 generateCubicPackingInSphere(); 00409 generateRandomFill(); 00410 } 00411 }; 00412 };