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 #ifndef ESYS_LSMPARTICLEFITTER_H 00014 #define ESYS_LSMPARTICLEFITTER_H 00015 00016 #include "Geometry/RandomBlockGenerator.h" 00017 #include "Geometry/Sphere3d.h" 00018 #include "Geometry/Sphere2d.h" 00019 00020 namespace esys 00021 { 00022 namespace lsm 00023 { 00024 class ParticleFitter 00025 { 00026 public: 00027 typedef RandomBlockGenerator::ParticleVector ParticleVector; 00028 00029 ParticleFitter(RandomBlockGenerator &blockGenerator) 00030 : m_pGenerator(&blockGenerator), 00031 m_successfulFitCount(0), 00032 m_getFitCount(0), 00033 m_failedFitCount(0) 00034 { 00035 } 00036 00037 virtual ~ParticleFitter() {} 00038 00039 virtual SimpleParticle getFitParticle( 00040 const SimpleParticle &particle, 00041 const ParticleVector &neighbours, 00042 const Plane &plane 00043 ) = 0; 00044 00045 static const SimpleParticle INVALID; 00046 00047 void incrGetFit() 00048 { 00049 m_getFitCount++; 00050 } 00051 00052 void incrFailedFit() 00053 { 00054 m_failedFitCount++; 00055 } 00056 00057 void incrSuccessfulFit() 00058 { 00059 m_successfulFitCount++; 00060 } 00061 00062 virtual std::string getName() const = 0; 00063 00064 void write(std::ostream &oStream) const 00065 { 00066 oStream 00067 << (getName() + ": ") 00068 << "fit count = " << m_getFitCount 00069 << ", success = " << m_successfulFitCount 00070 << ", fail = " << m_failedFitCount; 00071 } 00072 00073 std::string toString() const 00074 { 00075 std::stringstream sStream; 00076 write(sStream); 00077 return sStream.str(); 00078 } 00079 00080 virtual bool particleFits(const SimpleParticle &particle) const 00081 { 00082 return getGenerator().particleFits(particle); 00083 } 00084 00085 protected: 00086 RandomBlockGenerator &getGenerator() 00087 { 00088 return *m_pGenerator; 00089 } 00090 00091 const RandomBlockGenerator &getGenerator() const 00092 { 00093 return *m_pGenerator; 00094 } 00095 private: 00096 RandomBlockGenerator *m_pGenerator; 00097 int m_successfulFitCount; 00098 int m_getFitCount; 00099 int m_failedFitCount; 00100 }; 00101 00102 inline std::ostream &operator<<(std::ostream &oStream, const ParticleFitter &fitter) 00103 { 00104 fitter.write(oStream); 00105 return oStream; 00106 } 00107 00108 class MoveToSurfaceFitter : public ParticleFitter 00109 { 00110 public: 00111 MoveToSurfaceFitter(RandomBlockGenerator &blockGenerator) 00112 : ParticleFitter(blockGenerator) 00113 { 00114 } 00115 00116 virtual std::string getName() const 00117 { 00118 return "Mv to Surface"; 00119 } 00120 00121 SimpleParticle moveToSurface( 00122 const SimpleParticle &stationary, 00123 const SimpleParticle &movable 00124 ) 00125 { 00126 SimpleParticle moved = movable; 00127 const Vec3 centreDiff = movable.getPos() - stationary.getPos(); 00128 const double centreDist = centreDiff.norm(); 00129 if (centreDist > 0.0) { 00130 const Vec3 newCentrePos = 00131 stationary.getPos() 00132 + 00133 (centreDiff/(centreDist))*(stationary.getRad() + movable.getRad()); 00134 moved.moveTo(newCentrePos); 00135 } 00136 return moved; 00137 } 00138 00139 virtual SimpleParticle getFitParticle( 00140 const SimpleParticle &particle, 00141 const ParticleVector &neighbours, 00142 const Plane &plane 00143 ) 00144 { 00145 incrGetFit(); 00146 SimpleParticle newParticle = particle; 00147 if (neighbours.size() > 0) { 00148 if ( 00149 (particle.getPos() - neighbours[0]->getPos()).norm() 00150 < 00151 (particle.getRad() + neighbours[0]->getRad()) 00152 ) { 00153 newParticle = moveToSurface(*(neighbours[0]), particle); 00154 } 00155 } 00156 if (newParticle.isValid() && !particleFits(newParticle)) { 00157 newParticle = INVALID; 00158 incrFailedFit(); 00159 } else if (newParticle.isValid()) { 00160 incrSuccessfulFit(); 00161 } 00162 return newParticle; 00163 } 00164 }; 00165 00166 class ThreeDParticleFitter : public ParticleFitter 00167 { 00168 public: 00169 ThreeDParticleFitter(RandomBlockGenerator &blockGenerator) 00170 : ParticleFitter(blockGenerator) 00171 { 00172 } 00173 00174 virtual std::string getName() const 00175 { 00176 return " 3D Fitter"; 00177 } 00178 00179 SimpleParticle findAFit( 00180 const SimpleParticle& Po, 00181 const ParticleVector &particleVector 00182 ) 00183 { 00184 SimpleParticle fittedParticle = SimpleParticle::INVALID; 00185 Vec3 M; 00186 double r; 00187 int id=Po.getID(); 00188 00189 if (particleVector.size() > 3) { 00190 const bool foundFit = 00191 Sphere3D::FillIn( 00192 particleVector[0]->getPos(), 00193 particleVector[1]->getPos(), 00194 particleVector[2]->getPos(), 00195 particleVector[3]->getPos(), 00196 particleVector[0]->getRad(), 00197 particleVector[1]->getRad(), 00198 particleVector[2]->getRad(), 00199 particleVector[3]->getRad(), 00200 M, 00201 r 00202 ); 00203 if (foundFit) { 00204 fittedParticle = SimpleParticle(M, r, id); 00205 } 00206 } else { 00207 throw std::runtime_error("findAFit: particleVector argument has fewer than 4 elements."); 00208 } 00209 return fittedParticle; 00210 } 00211 00212 virtual SimpleParticle getFitParticle( 00213 const SimpleParticle &particle, 00214 const ParticleVector &neighbours, 00215 const Plane &plane 00216 ) 00217 { 00218 incrGetFit(); 00219 SimpleParticle newParticle = INVALID; 00220 if(neighbours.size() > 3){ // at least 4 neighbors 00221 const SimpleParticle &Pi = *(neighbours[0]); 00222 const double ndist=(particle.getPos()-Pi.getPos()).norm(); 00223 if (ndist > 0.0) { 00224 newParticle = particle; 00225 if (ndist < Pi.getRad()){ // if Po inside Pi -> move Po to the surface of Pi 00226 newParticle.moveTo( 00227 Pi.getPos()+((particle.getPos()-Pi.getPos())*(Pi.getRad()/ndist)) 00228 ); 00229 } 00230 const double dist_p = plane.sep(newParticle.getPos()); 00231 const double dist_3 = (newParticle.getPos()-neighbours[3]->getPos()).norm()-neighbours[3]->getRad(); 00232 if (dist_p > dist_3) { // 4th particle closer than plane -> fit 4 particles 00233 newParticle = findAFit(newParticle, neighbours); 00234 } 00235 else { 00236 newParticle = INVALID; 00237 } 00238 } 00239 } 00240 if (newParticle.isValid() && !particleFits(newParticle)) { 00241 newParticle = INVALID; 00242 incrFailedFit(); 00243 } else if (newParticle.isValid()) { 00244 incrSuccessfulFit(); 00245 } 00246 return newParticle; 00247 } 00248 }; 00249 00250 class TwoDParticleFitter : public ParticleFitter 00251 { 00252 public: 00253 TwoDParticleFitter(RandomBlockGenerator &blockGenerator) 00254 : ParticleFitter(blockGenerator) 00255 { 00256 } 00257 00258 virtual std::string getName() const 00259 { 00260 return " 2D Fitter"; 00261 } 00262 00263 SimpleParticle findAFit( 00264 const SimpleParticle& Po, 00265 const ParticleVector &particleVector 00266 ) 00267 { 00268 SimpleParticle fittedParticle = SimpleParticle::INVALID; 00269 Vec3 M; 00270 double r; 00271 int id=Po.getID(); 00272 00273 if (particleVector.size() > 2) { 00274 const bool foundFit = 00275 Sphere2D::FillIn( 00276 particleVector[0]->getPos(), 00277 particleVector[1]->getPos(), 00278 particleVector[2]->getPos(), 00279 particleVector[0]->getRad(), 00280 particleVector[1]->getRad(), 00281 particleVector[2]->getRad(), 00282 M, 00283 r 00284 ); 00285 if (foundFit) { 00286 fittedParticle = SimpleParticle(M, r, id); 00287 } 00288 } else { 00289 throw std::runtime_error("findAFit: particleVector argument has fewer than 3 elements."); 00290 } 00291 return fittedParticle; 00292 } 00293 00294 virtual SimpleParticle getFitParticle( 00295 const SimpleParticle &particle, 00296 const ParticleVector &neighbours, 00297 const Plane &plane 00298 ) 00299 { 00300 incrGetFit(); 00301 SimpleParticle newParticle = INVALID; 00302 if(neighbours.size() > 2){ // at least 3 neighbors 00303 const SimpleParticle &Pi = *(neighbours[0]); 00304 const double ndist=(particle.getPos()-Pi.getPos()).norm(); 00305 if (ndist > 0.0) { 00306 newParticle = particle; 00307 if (ndist < Pi.getRad()){ // if Po inside Pi -> move Po to the surface of Pi 00308 newParticle.moveTo( 00309 Pi.getPos()+((particle.getPos()-Pi.getPos())*(Pi.getRad()/ndist)) 00310 ); 00311 } 00312 const double dist_p = plane.sep(newParticle.getPos()); 00313 const double dist_2 = (newParticle.getPos()-neighbours[2]->getPos()).norm()-neighbours[2]->getRad(); 00314 if (dist_p > dist_2) { // 4th particle closer than plane -> fit 4 particles 00315 newParticle = findAFit(newParticle, neighbours); 00316 } 00317 else { 00318 newParticle = INVALID; 00319 } 00320 } 00321 } 00322 if (newParticle.isValid() && !particleFits(newParticle)) { 00323 newParticle = INVALID; 00324 incrFailedFit(); 00325 } else if (newParticle.isValid()) { 00326 incrSuccessfulFit(); 00327 } 00328 return newParticle; 00329 } 00330 }; 00331 00332 class TwoDPlaneParticleFitter : public ParticleFitter 00333 { 00334 public: 00335 TwoDPlaneParticleFitter(RandomBlockGenerator &blockGenerator) 00336 : ParticleFitter(blockGenerator) 00337 { 00338 } 00339 00340 virtual std::string getName() const 00341 { 00342 return " 2D Plane"; 00343 } 00344 00345 SimpleParticle findAFit( 00346 const SimpleParticle& particle, 00347 const ParticleVector &particleVector, 00348 const Plane& plane 00349 ) 00350 { 00351 SimpleParticle fittedParticle = SimpleParticle::INVALID; 00352 Vec3 M; 00353 double r; 00354 const int id = particle.getID(); 00355 00356 if (particleVector.size() > 1) { 00357 const bool foundFit = 00358 Sphere2D::FillInWP( 00359 particleVector[0]->getPos(), 00360 particleVector[1]->getPos(), 00361 plane.GetO(), 00362 plane.GetW(), 00363 particleVector[0]->getRad(), 00364 particleVector[1]->getRad(), 00365 M, 00366 r 00367 ); 00368 if (foundFit) { 00369 fittedParticle = SimpleParticle(M, r, id); 00370 } 00371 } else { 00372 throw std::runtime_error("findAFit: particleVector vector contains less than 2 particles."); 00373 } 00374 00375 return fittedParticle; 00376 } 00377 00378 virtual SimpleParticle getFitParticle( 00379 const SimpleParticle &particle, 00380 const ParticleVector &neighbours, 00381 const Plane &plane 00382 ) 00383 { 00384 incrGetFit(); 00385 SimpleParticle newParticle = INVALID; 00386 if (neighbours.size() >= 2) { // 2 neighbors -> try 2 particles + plane 00387 const SimpleParticle &Pi = *(neighbours[0]); 00388 const double ndist=(particle.getPos()-Pi.getPos()).norm(); 00389 if ( 00390 (ndist > 0.0) 00391 && 00392 ( 00393 (neighbours.size() == 2) 00394 || 00395 ( 00396 plane.sep(particle.getPos()) 00397 < 00398 (particle.getPos()-neighbours[2]->getPos()).norm()-neighbours[2]->getRad() 00399 ) 00400 ) 00401 ) { 00402 newParticle = particle; 00403 if (ndist < Pi.getRad()) { // if Po inside Pi -> move Po to the surface of Pi 00404 newParticle.moveTo( 00405 Pi.getPos() + ((particle.getPos() - Pi.getPos()) * (Pi.getRad()/ndist)) 00406 ); 00407 } 00408 newParticle = findAFit(newParticle, neighbours, plane); 00409 } 00410 } 00411 if (newParticle.isValid() && !particleFits(newParticle)) { 00412 newParticle = INVALID; 00413 incrFailedFit(); 00414 } else if (newParticle.isValid()) { 00415 incrSuccessfulFit(); 00416 } 00417 return newParticle; 00418 } 00419 }; 00420 00421 class ThreeDPlaneParticleFitter : public ParticleFitter 00422 { 00423 public: 00424 ThreeDPlaneParticleFitter(RandomBlockGenerator &blockGenerator) 00425 : ParticleFitter(blockGenerator) 00426 { 00427 } 00428 00429 virtual std::string getName() const 00430 { 00431 return " 3D Plane"; 00432 } 00433 00434 SimpleParticle findAFit( 00435 const SimpleParticle& particle, 00436 const ParticleVector &particleVector, 00437 const Plane& plane 00438 ) 00439 { 00440 SimpleParticle fittedParticle = SimpleParticle::INVALID; 00441 Vec3 M; 00442 double r; 00443 const int id = particle.getID(); 00444 00445 if (particleVector.size() > 2) { 00446 const bool foundFit = 00447 Sphere3D::FillInWP( 00448 particleVector[0]->getPos(), 00449 particleVector[1]->getPos(), 00450 particleVector[2]->getPos(), 00451 plane.GetO(), 00452 plane.GetW(), 00453 particleVector[0]->getRad(), 00454 particleVector[1]->getRad(), 00455 particleVector[2]->getRad(), 00456 M, 00457 r 00458 ); 00459 if (foundFit) { 00460 fittedParticle = SimpleParticle(M, r, id); 00461 } 00462 } else { 00463 throw std::runtime_error("findAFit: particleVector vector contains less than 3 particles."); 00464 } 00465 00466 return fittedParticle; 00467 } 00468 00469 virtual SimpleParticle getFitParticle( 00470 const SimpleParticle &particle, 00471 const ParticleVector &neighbours, 00472 const Plane &plane 00473 ) 00474 { 00475 incrGetFit(); 00476 SimpleParticle newParticle = INVALID; 00477 if (neighbours.size() >= 3) { // 3 neighbors -> try 3 particles + plane 00478 const SimpleParticle &Pi = *(neighbours[0]); 00479 const double ndist=(particle.getPos()-Pi.getPos()).norm(); 00480 if ( 00481 (ndist > 0.0) 00482 && 00483 ( 00484 (neighbours.size() == 3) 00485 || 00486 ( 00487 plane.sep(particle.getPos()) 00488 < 00489 (particle.getPos()-neighbours[3]->getPos()).norm()-neighbours[3]->getRad() 00490 ) 00491 ) 00492 ) { 00493 newParticle = particle; 00494 if (ndist < Pi.getRad()) { // if Po inside Pi -> move Po to the surface of Pi 00495 newParticle.moveTo( 00496 Pi.getPos() + ((particle.getPos() - Pi.getPos()) * (Pi.getRad()/ndist)) 00497 ); 00498 } 00499 newParticle = findAFit(newParticle, neighbours, plane); 00500 } 00501 } 00502 if (newParticle.isValid() && !particleFits(newParticle)) { 00503 newParticle = INVALID; 00504 incrFailedFit(); 00505 } else if (newParticle.isValid()) { 00506 incrSuccessfulFit(); 00507 } 00508 return newParticle; 00509 } 00510 }; 00511 }; 00512 }; 00513 00514 #endif