ESyS-Particle  2.3.2
SphereFitter.h
Go to the documentation of this file.
1 
2 // //
3 // Copyright (c) 2003-2014 by The University of Queensland //
4 // Centre for Geoscience Computing //
5 // http://earth.uq.edu.au/centre-geoscience-computing //
6 // //
7 // Primary Business: Brisbane, Queensland, Australia //
8 // Licensed under the Open Software License version 3.0 //
9 // http://www.opensource.org/licenses/osl-3.0.php //
10 // //
12 
13 
14 #ifndef ESYS_LSMSPHEREFITTER_H
15 #define ESYS_LSMSPHEREFITTER_H
16 
17 #include "Geometry/Sphere3d.h"
18 #include "Geometry/Sphere2d.h"
20 
21 #include <stdexcept>
22 #include <sstream>
23 
24 namespace esys
25 {
26  namespace lsm
27  {
28  template <typename TmplFitTraits>
30  {
31  public:
32  typedef TmplFitTraits FitTraits;
33  typedef typename FitTraits::Validator Validator;
34  typedef typename FitTraits::Particle Particle;
36  typedef typename FitTraits::Plane3D Plane3D;
37 
38  SphereFitter(const std::string &name, Validator &validator)
39  : m_pValidator(&validator),
41  m_getFitCount(0),
43  m_name(name)
44  {
45  }
46 
47  virtual ~SphereFitter()
48  {
49  }
50 
51  virtual Particle getFitParticle(
52  const Particle &particle,
53  const ParticleVector &neighbours,
54  const Plane3D &plane
55  ) = 0;
56 
58  {
59  return Particle::INVALID;
60  }
61 
62  void incrGetFit()
63  {
64  m_getFitCount++;
65  }
66 
68  {
70  }
71 
73  {
75  }
76 
77  const std::string &getName() const
78  {
79  return m_name;
80  }
81 
82  void write(std::ostream &oStream) const
83  {
84  oStream
85  << (getName() + ": ")
86  << "fit count = " << m_getFitCount
87  << ", success = " << m_successfulFitCount
88  << ", fail = " << m_failedFitCount;
89  }
90 
91  std::string toString() const
92  {
93  std::stringstream sStream;
94  write(sStream);
95  return sStream.str();
96  }
97 
98  bool particleIsValid(const Particle &particle) const
99  {
100  return getValidator().particleIsValid(particle);
101  }
102 
103  protected:
105  {
106  return *m_pValidator;
107  }
108 
109  const Validator &getValidator() const
110  {
111  return *m_pValidator;
112  }
113  private:
118  std::string m_name;
119  };
120 
121  template <typename TmplFitTraits>
122  inline std::ostream &operator<<(
123  std::ostream &oStream,
124  const SphereFitter<TmplFitTraits> &fitter
125  )
126  {
127  fitter.write(oStream);
128  return oStream;
129  }
130 
131  template <typename TmplFitTraits>
132  class MoveToSurfaceFitter : public SphereFitter<TmplFitTraits>
133  {
134  public:
136  typedef typename Inherited::Validator Validator;
137  typedef typename Inherited::Particle Particle;
139  typedef typename Inherited::Plane3D Plane3D;
141  : Inherited("Mv to Surface", validator)
142  {
143  }
144 
146  const Particle &stationary,
147  const Particle &movable
148  )
149  {
150  Particle moved = movable;
151  const Vec3 centreDiff = movable.getPos() - stationary.getPos();
152  const double centreDist = centreDiff.norm();
153  if (centreDist > 0.0) {
154  const Vec3 newCentrePos =
155  stationary.getPos()
156  +
157  (centreDiff/(centreDist))*(stationary.getRad() + movable.getRad());
158  moved.moveTo(newCentrePos);
159  }
160  return moved;
161  }
162 
164  const Particle &particle,
165  const ParticleVector &neighbours,
166  const Plane3D &plane
167  )
168  {
169  this->incrGetFit();
170  Particle newParticle = particle;
171  if (neighbours.size() > 0) {
172  if (
173  (particle.getPos() - neighbours[0]->getPos()).norm()
174  <
175  (particle.getRad() + neighbours[0]->getRad())
176  ) {
177  newParticle = moveToSurface(*(neighbours[0]), particle);
178  }
179  }
180  if (newParticle.isValid() && !this->particleIsValid(newParticle)) {
181  newParticle = this->getInvalidParticle();
182  this->incrFailedFit();
183  } else if (newParticle.isValid()) {
184  this->incrSuccessfulFit();
185  }
186  return newParticle;
187  }
188  };
189 
190  template <typename TmplFitTraits>
191  class ThreeDSphereFitter : public SphereFitter<TmplFitTraits>
192  {
193  public:
195  typedef typename Inherited::Validator Validator;
196  typedef typename Inherited::Particle Particle;
198  typedef typename Inherited::Plane3D Plane3D;
199 
201  : Inherited(" 3D Fitter", validator)
202  {
203  }
204 
206  const Particle& Po,
207  const ParticleVector &particleVector
208  )
209  {
210  Particle fittedParticle = this->getInvalidParticle();
211  Vec3 M;
212  double r;
213  int id=Po.getID();
214 
215  if (particleVector.size() > 3) {
216  const bool foundFit =
218  particleVector[0]->getPos(),
219  particleVector[1]->getPos(),
220  particleVector[2]->getPos(),
221  particleVector[3]->getPos(),
222  particleVector[0]->getRad(),
223  particleVector[1]->getRad(),
224  particleVector[2]->getRad(),
225  particleVector[3]->getRad(),
226  M,
227  r
228  );
229  if (foundFit) {
230  fittedParticle = Particle(M, r, id);
231  }
232  } else {
233  throw std::runtime_error("findAFit: particleVector argument has fewer than 4 elements.");
234  }
235  return fittedParticle;
236  }
237 
239  const Particle &particle,
240  const ParticleVector &neighbours,
241  const Plane3D &plane
242  )
243  {
244  this->incrGetFit();
245  Particle newParticle = this->getInvalidParticle();
246  if(neighbours.size() > 3){ // at least 4 neighbors
247  const Particle &Pi = *(neighbours[0]);
248  const double ndist=(particle.getPos()-Pi.getPos()).norm();
249  if (ndist > 0.0) {
250  newParticle = particle;
251  if (ndist < Pi.getRad()){ // if Po inside Pi -> move Po to the surface of Pi
252  newParticle.moveTo(
253  Pi.getPos()+((particle.getPos()-Pi.getPos())*(Pi.getRad()/ndist))
254  );
255  }
256  const double dist_p = plane.sep(newParticle.getPos());
257  const double dist_3 = (newParticle.getPos()-neighbours[3]->getPos()).norm()-neighbours[3]->getRad();
258  if (dist_p > dist_3) { // 4th particle closer than plane -> fit 4 particles
259  newParticle = findAFit(newParticle, neighbours);
260  }
261  else {
262  newParticle = this->getInvalidParticle();
263  }
264  }
265  }
266  if (newParticle.isValid() && !this->particleIsValid(newParticle)) {
267  newParticle = this->getInvalidParticle();
268  this->incrFailedFit();
269  } else if (newParticle.isValid()) {
270  this->incrSuccessfulFit();
271  }
272  return newParticle;
273  }
274  };
275 
276  template <typename TmplFitTraits>
277  class TwoDSphereFitter : public SphereFitter<TmplFitTraits>
278  {
279  public:
281  typedef typename Inherited::Validator Validator;
282  typedef typename Inherited::Particle Particle;
284  typedef typename Inherited::Plane3D Plane3D;
285 
287  : Inherited(" 2D Fitter", validator)
288  {
289  }
290 
292  const Particle& Po,
293  const ParticleVector &particleVector
294  )
295  {
296  Particle fittedParticle = this->getInvalidParticle();
297  Vec3 M;
298  double r;
299  int id=Po.getID();
300 
301  if (particleVector.size() > 2) {
302  const bool foundFit =
304  particleVector[0]->getPos(),
305  particleVector[1]->getPos(),
306  particleVector[2]->getPos(),
307  particleVector[0]->getRad(),
308  particleVector[1]->getRad(),
309  particleVector[2]->getRad(),
310  M,
311  r
312  );
313  if (foundFit) {
314  fittedParticle = Particle(M, r, id);
315  }
316  } else {
317  throw std::runtime_error("findAFit: particleVector argument has fewer than 3 elements.");
318  }
319  return fittedParticle;
320  }
321 
323  const Particle &particle,
324  const ParticleVector &neighbours,
325  const Plane3D &plane
326  )
327  {
328  this->incrGetFit();
329  Particle newParticle = this->getInvalidParticle();
330  if(neighbours.size() > 2){ // at least 3 neighbors
331  const Particle &Pi = *(neighbours[0]);
332  const double ndist=(particle.getPos()-Pi.getPos()).norm();
333  if (ndist > 0.0) {
334  newParticle = particle;
335  if (ndist < Pi.getRad()){ // if Po inside Pi -> move Po to the surface of Pi
336  newParticle.moveTo(
337  Pi.getPos()+((particle.getPos()-Pi.getPos())*(Pi.getRad()/ndist))
338  );
339  }
340  const double dist_p = plane.sep(newParticle.getPos());
341  const double dist_2 = (newParticle.getPos()-neighbours[2]->getPos()).norm()-neighbours[2]->getRad();
342  if (dist_p > dist_2) { // 3rd particle closer than plane -> fit 3 particles
343  newParticle = findAFit(newParticle, neighbours);
344  }
345  else {
346  newParticle = this->getInvalidParticle();
347  }
348  }
349  }
350  if (newParticle.isValid() && !this->particleIsValid(newParticle)) {
351  newParticle = this->getInvalidParticle();
352  this->incrFailedFit();
353  } else if (newParticle.isValid()) {
354  this->incrSuccessfulFit();
355  }
356  return newParticle;
357  }
358  };
359 
360  template <typename TmplFitTraits>
361  class TwoDSphereSphereFitter : public SphereFitter<TmplFitTraits>
362  {
363  public:
365  typedef typename Inherited::Validator Validator;
366  typedef typename Inherited::Particle Particle;
368  typedef typename Inherited::Plane3D Plane3D;
369 
371  Validator &validator,
372  const BoundingSphere &sphere
373  )
374  : Inherited(" 2D Sph Fitter", validator),
375  m_sphere(sphere)
376  {
377  }
378 
380  const Particle& Po,
381  const ParticleVector &particleVector
382  )
383  {
384  Particle fittedParticle = this->getInvalidParticle();
385  Vec3 M;
386  double r;
387  int id=Po.getID();
388 
389  if (particleVector.size() >= 2) {
390  const bool foundFit =
393  particleVector[0]->getPos(),
394  particleVector[1]->getPos(),
395  -m_sphere.getRadius(),
396  particleVector[0]->getRad(),
397  particleVector[1]->getRad(),
398  M,
399  r
400  );
401  if (foundFit) {
402  fittedParticle = Particle(M, r, id);
403  }
404  } else {
405  throw std::runtime_error(
406  "TwoDSphereSphereFitter::findAFit: particleVector "
407  "argument has fewer than 2 elements."
408  );
409  }
410  return fittedParticle;
411  }
412 
414  const Particle &particle,
415  const ParticleVector &neighbours,
416  const Plane3D &plane
417  )
418  {
419  this->incrGetFit();
420  Particle newParticle = this->getInvalidParticle();
421  if (neighbours.size() >= 2) { // 2 neighbours: try 2 particles + bounding-sphere
422  const Particle &Pi = *(neighbours[0]);
423  const double ndist=(particle.getPos()-Pi.getPos()).norm();
424  if (
425  (ndist > 0.0)
426  &&
427  (
428  (neighbours.size() == 2)
429  ||
430  (
431  fabs((m_sphere.getCentre()-particle.getPos()).norm()-m_sphere.getRadius())
432  <
433  (particle.getPos()-neighbours[2]->getPos()).norm()-neighbours[2]->getRad()
434  )
435  )
436  ) {
437  newParticle = particle;
438  if (ndist < Pi.getRad()) { // if Po inside Pi -> move Po to the surface of Pi
439  newParticle.moveTo(
440  Pi.getPos() + ((particle.getPos() - Pi.getPos()) * (Pi.getRad()/ndist))
441  );
442  }
443  newParticle = findAFit(newParticle, neighbours);
444  }
445  }
446  if (newParticle.isValid() && !this->particleIsValid(newParticle)) {
447  newParticle = this->getInvalidParticle();
448  this->incrFailedFit();
449  } else if (newParticle.isValid()) {
450  this->incrSuccessfulFit();
451  }
452  return newParticle;
453  }
454 
455  private:
457  };
458 
459  template <typename TmplFitTraits>
460  class ThreeDSphereSphereFitter : public SphereFitter<TmplFitTraits>
461  {
462  public:
464  typedef typename Inherited::Validator Validator;
465  typedef typename Inherited::Particle Particle;
467  typedef typename Inherited::Plane3D Plane3D;
468 
470  Validator &validator,
471  const BoundingSphere &sphere
472  )
473  : Inherited(" 3D Sph Fitter", validator),
474  m_sphere(sphere)
475  {
476  }
477 
479  const Particle& Po,
480  const ParticleVector &particleVector
481  )
482  {
483  Particle fittedParticle = this->getInvalidParticle();
484  Vec3 M;
485  double r;
486  int id=Po.getID();
487 
488  if (particleVector.size() > 2) {
489  const bool foundFit =
492  particleVector[0]->getPos(),
493  particleVector[1]->getPos(),
494  particleVector[2]->getPos(),
495  -m_sphere.getRadius(),
496  particleVector[0]->getRad(),
497  particleVector[1]->getRad(),
498  particleVector[2]->getRad(),
499  M,
500  r
501  );
502  if (foundFit) {
503  fittedParticle = Particle(M, r, id);
504  }
505  } else {
506  throw
507  std::runtime_error(
508  "ThreeDSphereSphereFitter::findAFit: particleVector argument"
509  " has fewer than 3 elements."
510  );
511  }
512  return fittedParticle;
513  }
514 
516  const Particle &particle,
517  const ParticleVector &neighbours,
518  const Plane3D &plane
519  )
520  {
521  this->incrGetFit();
522  Particle newParticle = this->getInvalidParticle();
523  if (neighbours.size() >= 3) { // 3 neighbours: try 3 particles + bounding-sphere
524  const Particle &Pi = *(neighbours[0]);
525  const double ndist=(particle.getPos()-Pi.getPos()).norm();
526  if (
527  (ndist > 0.0)
528  &&
529  (
530  (neighbours.size() == 3)
531  ||
532  (
533  fabs((m_sphere.getCentre()-particle.getPos()).norm()-m_sphere.getRadius())
534  <
535  (particle.getPos()-neighbours[3]->getPos()).norm()-neighbours[3]->getRad()
536  )
537  )
538  ) {
539  newParticle = particle;
540  if (ndist < Pi.getRad()) { // if Po inside Pi -> move Po to the surface of Pi
541  newParticle.moveTo(
542  Pi.getPos() + ((particle.getPos() - Pi.getPos()) * (Pi.getRad()/ndist))
543  );
544  }
545  newParticle = findAFit(newParticle, neighbours);
546  }
547  }
548  if (newParticle.isValid() && !this->particleIsValid(newParticle)) {
549  newParticle = this->getInvalidParticle();
550  this->incrFailedFit();
551  } else if (newParticle.isValid()) {
552  this->incrSuccessfulFit();
553  }
554  return newParticle;
555  }
556 
557  private:
559  };
560 
561  template <typename TmplFitTraits>
562  class TwoDPlaneSphereFitter : public SphereFitter<TmplFitTraits>
563  {
564  public:
566  typedef typename Inherited::Validator Validator;
567  typedef typename Inherited::Particle Particle;
569  typedef typename Inherited::Plane3D Plane3D;
570 
572  : Inherited(" 2D Plane", validator)
573  {
574  }
575 
577  const Particle& particle,
578  const ParticleVector &particleVector,
579  const Plane3D& plane
580  )
581  {
582  Particle fittedParticle = this->getInvalidParticle();
583  Vec3 M;
584  double r;
585  const int id = particle.getID();
586 
587  if (particleVector.size() > 1) {
588  const bool foundFit =
590  particleVector[0]->getPos(),
591  particleVector[1]->getPos(),
592  plane.GetO(),
593  plane.GetW(),
594  particleVector[0]->getRad(),
595  particleVector[1]->getRad(),
596  M,
597  r
598  );
599  if (foundFit) {
600  fittedParticle = Particle(M, r, id);
601  }
602  } else {
603  throw std::runtime_error("findAFit: particleVector vector contains less than 2 particles.");
604  }
605 
606  return fittedParticle;
607  }
608 
610  const Particle &particle,
611  const ParticleVector &neighbours,
612  const Plane3D &plane
613  )
614  {
615  this->incrGetFit();
616  Particle newParticle = this->getInvalidParticle();
617  if (neighbours.size() >= 2) { // 2 neighbors -> try 2 particles + plane
618  const Particle &Pi = *(neighbours[0]);
619  const double ndist=(particle.getPos()-Pi.getPos()).norm();
620  if (
621  (ndist > 0.0)
622  &&
623  (
624  (neighbours.size() == 2)
625  ||
626  (
627  plane.sep(particle.getPos())
628  <
629  (particle.getPos()-neighbours[2]->getPos()).norm()-neighbours[2]->getRad()
630  )
631  )
632  ) {
633  newParticle = particle;
634  if (ndist < Pi.getRad()) { // if Po inside Pi -> move Po to the surface of Pi
635  newParticle.moveTo(
636  Pi.getPos() + ((particle.getPos() - Pi.getPos()) * (Pi.getRad()/ndist))
637  );
638  }
639  newParticle = findAFit(newParticle, neighbours, plane);
640  }
641  }
642  if (newParticle.isValid() && !this->particleIsValid(newParticle)) {
643  newParticle = this->getInvalidParticle();
644  this->incrFailedFit();
645  } else if (newParticle.isValid()) {
646  this->incrSuccessfulFit();
647  }
648  return newParticle;
649  }
650  };
651 
652  template <typename TmplFitTraits>
653  class ThreeDPlaneSphereFitter : public SphereFitter<TmplFitTraits>
654  {
655  public:
657  typedef typename Inherited::Validator Validator;
658  typedef typename Inherited::Particle Particle;
660  typedef typename Inherited::Plane3D Plane3D;
661 
663  : Inherited(" 3D Plane", validator)
664  {
665  }
666 
668  const Particle& particle,
669  const ParticleVector &particleVector,
670  const Plane3D& plane
671  )
672  {
673  Particle fittedParticle = this->getInvalidParticle();
674  Vec3 M;
675  double r;
676  const int id = particle.getID();
677 
678  if (particleVector.size() > 2) {
679  const bool foundFit =
681  particleVector[0]->getPos(),
682  particleVector[1]->getPos(),
683  particleVector[2]->getPos(),
684  plane.GetO(),
685  plane.GetW(),
686  particleVector[0]->getRad(),
687  particleVector[1]->getRad(),
688  particleVector[2]->getRad(),
689  M,
690  r
691  );
692  if (foundFit) {
693  fittedParticle = Particle(M, r, id);
694  }
695  } else {
696  throw std::runtime_error("findAFit: particleVector vector contains less than 3 particles.");
697  }
698 
699  return fittedParticle;
700  }
701 
703  const Particle &particle,
704  const ParticleVector &neighbours,
705  const Plane3D &plane
706  )
707  {
708  this->incrGetFit();
709  Particle newParticle = this->getInvalidParticle();
710  if (neighbours.size() >= 3) { // 3 neighbors -> try 3 particles + plane
711  const Particle &Pi = *(neighbours[0]);
712  const double ndist=(particle.getPos()-Pi.getPos()).norm();
713  if (
714  (ndist > 0.0)
715  &&
716  (
717  (neighbours.size() == 3)
718  ||
719  (
720  plane.sep(particle.getPos())
721  <
722  (particle.getPos()-neighbours[3]->getPos()).norm()-neighbours[3]->getRad()
723  )
724  )
725  ) {
726  newParticle = particle;
727  if (ndist < Pi.getRad()) { // if Po inside Pi -> move Po to the surface of Pi
728  newParticle.moveTo(
729  Pi.getPos() + ((particle.getPos() - Pi.getPos()) * (Pi.getRad()/ndist))
730  );
731  }
732  newParticle = findAFit(newParticle, neighbours, plane);
733  }
734  }
735  if (newParticle.isValid() && !this->particleIsValid(newParticle)) {
736  newParticle = this->getInvalidParticle();
737  this->incrFailedFit();
738  } else if (newParticle.isValid()) {
739  this->incrSuccessfulFit();
740  }
741  return newParticle;
742  }
743  };
744  }
745 }
746 
747 #include "Geometry/SphereFitter.hpp"
748 
749 #endif