13 #ifndef ESYS_LSMPARTICLEFITTER_H
14 #define ESYS_LSMPARTICLEFITTER_H
16 #include "Geometry/RandomBlockGenerator.h"
17 #include "Geometry/Sphere3d.h"
18 #include "Geometry/Sphere2d.h"
27 typedef RandomBlockGenerator::ParticleVector ParticleVector;
30 : m_pGenerator(&blockGenerator),
31 m_successfulFitCount(0),
41 const ParticleVector &neighbours,
57 void incrSuccessfulFit()
59 m_successfulFitCount++;
62 virtual std::string getName()
const = 0;
64 void write(std::ostream &oStream)
const
68 <<
"fit count = " << m_getFitCount
69 <<
", success = " << m_successfulFitCount
70 <<
", fail = " << m_failedFitCount;
73 std::string toString()
const
75 std::stringstream sStream;
82 return getGenerator().particleFits(particle);
97 int m_successfulFitCount;
102 inline std::ostream &operator<<(std::ostream &oStream,
const ParticleFitter &fitter)
104 fitter.write(oStream);
116 virtual std::string getName()
const
118 return "Mv to Surface";
127 const Vec3 centreDiff = movable.getPos() - stationary.getPos();
128 const double centreDist = centreDiff.norm();
129 if (centreDist > 0.0) {
130 const Vec3 newCentrePos =
133 (centreDiff/(centreDist))*(stationary.getRad() + movable.getRad());
134 moved.moveTo(newCentrePos);
141 const ParticleVector &neighbours,
147 if (neighbours.size() > 0) {
149 (particle.getPos() - neighbours[0]->getPos()).norm()
151 (particle.getRad() + neighbours[0]->getRad())
153 newParticle = moveToSurface(*(neighbours[0]), particle);
156 if (newParticle.isValid() && !particleFits(newParticle)) {
157 newParticle = INVALID;
159 }
else if (newParticle.isValid()) {
174 virtual std::string getName()
const
181 const ParticleVector &particleVector
189 if (particleVector.size() > 3) {
190 const bool foundFit =
192 particleVector[0]->getPos(),
193 particleVector[1]->getPos(),
194 particleVector[2]->getPos(),
195 particleVector[3]->getPos(),
196 particleVector[0]->getRad(),
197 particleVector[1]->getRad(),
198 particleVector[2]->getRad(),
199 particleVector[3]->getRad(),
207 throw std::runtime_error(
"findAFit: particleVector argument has fewer than 4 elements.");
209 return fittedParticle;
214 const ParticleVector &neighbours,
220 if(neighbours.size() > 3){
222 const double ndist=(particle.getPos()-Pi.getPos()).norm();
224 newParticle = particle;
225 if (ndist < Pi.getRad()){
227 Pi.getPos()+((particle.getPos()-Pi.getPos())*(Pi.getRad()/ndist))
230 const double dist_p = plane.sep(newParticle.getPos());
231 const double dist_3 = (newParticle.getPos()-neighbours[3]->getPos()).norm()-neighbours[3]->getRad();
232 if (dist_p > dist_3) {
233 newParticle = findAFit(newParticle, neighbours);
236 newParticle = INVALID;
240 if (newParticle.isValid() && !particleFits(newParticle)) {
241 newParticle = INVALID;
243 }
else if (newParticle.isValid()) {
258 virtual std::string getName()
const
265 const ParticleVector &particleVector
273 if (particleVector.size() > 2) {
274 const bool foundFit =
276 particleVector[0]->getPos(),
277 particleVector[1]->getPos(),
278 particleVector[2]->getPos(),
279 particleVector[0]->getRad(),
280 particleVector[1]->getRad(),
281 particleVector[2]->getRad(),
289 throw std::runtime_error(
"findAFit: particleVector argument has fewer than 3 elements.");
291 return fittedParticle;
296 const ParticleVector &neighbours,
302 if(neighbours.size() > 2){
304 const double ndist=(particle.getPos()-Pi.getPos()).norm();
306 newParticle = particle;
307 if (ndist < Pi.getRad()){
309 Pi.getPos()+((particle.getPos()-Pi.getPos())*(Pi.getRad()/ndist))
312 const double dist_p = plane.sep(newParticle.getPos());
313 const double dist_2 = (newParticle.getPos()-neighbours[2]->getPos()).norm()-neighbours[2]->getRad();
314 if (dist_p > dist_2) {
315 newParticle = findAFit(newParticle, neighbours);
318 newParticle = INVALID;
322 if (newParticle.isValid() && !particleFits(newParticle)) {
323 newParticle = INVALID;
325 }
else if (newParticle.isValid()) {
340 virtual std::string getName()
const
347 const ParticleVector &particleVector,
354 const int id = particle.getID();
356 if (particleVector.size() > 1) {
357 const bool foundFit =
359 particleVector[0]->getPos(),
360 particleVector[1]->getPos(),
363 particleVector[0]->getRad(),
364 particleVector[1]->getRad(),
372 throw std::runtime_error(
"findAFit: particleVector vector contains less than 2 particles.");
375 return fittedParticle;
380 const ParticleVector &neighbours,
386 if (neighbours.size() >= 2) {
388 const double ndist=(particle.getPos()-Pi.getPos()).norm();
393 (neighbours.size() == 2)
396 plane.sep(particle.getPos())
398 (particle.getPos()-neighbours[2]->getPos()).norm()-neighbours[2]->getRad()
402 newParticle = particle;
403 if (ndist < Pi.getRad()) {
405 Pi.getPos() + ((particle.getPos() - Pi.getPos()) * (Pi.getRad()/ndist))
408 newParticle = findAFit(newParticle, neighbours, plane);
411 if (newParticle.isValid() && !particleFits(newParticle)) {
412 newParticle = INVALID;
414 }
else if (newParticle.isValid()) {
429 virtual std::string getName()
const
436 const ParticleVector &particleVector,
443 const int id = particle.getID();
445 if (particleVector.size() > 2) {
446 const bool foundFit =
448 particleVector[0]->getPos(),
449 particleVector[1]->getPos(),
450 particleVector[2]->getPos(),
453 particleVector[0]->getRad(),
454 particleVector[1]->getRad(),
455 particleVector[2]->getRad(),
463 throw std::runtime_error(
"findAFit: particleVector vector contains less than 3 particles.");
466 return fittedParticle;
471 const ParticleVector &neighbours,
477 if (neighbours.size() >= 3) {
479 const double ndist=(particle.getPos()-Pi.getPos()).norm();
484 (neighbours.size() == 3)
487 plane.sep(particle.getPos())
489 (particle.getPos()-neighbours[3]->getPos()).norm()-neighbours[3]->getRad()
493 newParticle = particle;
494 if (ndist < Pi.getRad()) {
496 Pi.getPos() + ((particle.getPos() - Pi.getPos()) * (Pi.getRad()/ndist))
499 newParticle = findAFit(newParticle, neighbours, plane);
502 if (newParticle.isValid() && !particleFits(newParticle)) {
503 newParticle = INVALID;
505 }
else if (newParticle.isValid()) {