ESyS-Particle  4.0.1
ParticleFitter.h
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