ESyS-Particle  4.0.1
RandomSpherePacker.hpp
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 };