Jpp  19.1.0-rc.1
the software that should make you happy
JEigen3D.hh
Go to the documentation of this file.
1 #ifndef __JEIGEN3D__
2 #define __JEIGEN3D__
3 
4 #include <cmath>
5 #include <map>
6 #include <limits>
7 #include <functional>
8 
9 #include "TMatrixDSymEigen.h"
10 
11 #include "JLang/JException.hh"
12 
13 #include "JMath/JMath.hh"
14 #include "JMath/JLegendre.hh"
15 
16 #include "JGeometry3D/JVector3D.hh"
19 
20 
21 /**
22  * \author mdejong
23  */
24 
25 namespace JGEOMETRY3D {}
26 namespace JPP { using namespace JGEOMETRY3D; }
27 
28 namespace JGEOMETRY3D {
29 
30  /**
31  * Eigen values in 3D.
32  */
33  struct JEigenValues3D :
34  public std::map<double, JVector3D, std::greater<double> >
35  {
36  /**
37  * Constructor.
38  *
39  * \param __begin begin of data
40  * \param __end end of data
41  */
42  template<class T>
43  JEigenValues3D(T __begin,
44  T __end)
45  {
46  const JCenter3D center(__begin, __end);
47 
48  // RMS matrix
49 
50  TMatrixDSym A(3);
51 
52  A = 0.0;
53 
54  for (T i = __begin; i != __end; ++i) {
55 
56  const double dx = center.getX() - i->getX();
57  const double dy = center.getY() - i->getY();
58  const double dz = center.getZ() - i->getZ();
59 
60  A(0,0) += (dx * dx);
61  A(0,1) += (dx * dy);
62  A(0,2) += (dx * dz);
63 
64  A(1,1) += (dy * dy);
65  A(1,2) += (dy * dz);
66 
67  A(2,2) += (dz * dz);
68  }
69 
70  A(1,0) = A(0,1);
71  A(2,0) = A(0,2);
72  A(2,1) = A(1,2);
73 
74  const TMatrixDSymEigen result(A);
75 
76  const TVectorD& u = result.GetEigenValues();
77  const TMatrixD& V = result.GetEigenVectors();
78 
79  for (Int_t i = 0; i != u.GetNoElements(); ++i) {
80  (*this)[u[i]] = JVector3D(V(0,i),
81  V(1,i),
82  V(2,i));
83  }
84  }
85  };
86 }
87 
88 namespace JMATH {
89 
92 
93 
94  /**
95  * Template spacialisation for averaging quaternions.
96  */
97  template<>
99 
100  /**
101  * Numerical precision for eigen value evaluation.
102  */
103  static double MINIMAL_EIGEN_VALUE;
104 
105  /**
106  * Default constructor.
107  */
109  A(4),
110  W(0.0)
111  {}
112 
113 
114  /**
115  * Constructor.
116  *
117  * \param __begin begin of data
118  * \param __end end of data
119  */
120  template<class T>
121  JAverage(T __begin, T __end) :
122  A(4),
123  W(0.0)
124  {
125  for (T i = __begin; i != __end; ++i) {
126  this->put(*i);
127  }
128  }
129 
130 
131  /**
132  * Reset.
133  */
134  void reset()
135  {
136  A = 0.0;
137  W = 0.0;
138  }
139 
140 
141  /**
142  * Type conversion operator.
143  *
144  * \return quaternion
145  */
146  operator JQuaternion3D() const
147  {
148  JQuaternion3D Q;
149 
150  const TMatrixDSymEigen result(A);
151 
152  const TVectorD& u = result.GetEigenValues();
153  const TMatrixD& V = result.GetEigenVectors();
154 
155  if (u.GetNoElements() != 0 && fabs(u[0]) > MINIMAL_EIGEN_VALUE) {
156 
157  Q = JQuaternion3D(V(0,0), V(1,0), V(2,0), V(3,0));
158 
159  if (Q.getA() < 0.0) {
160  Q.negate();
161  }
162 
163  Q.pow(u[0] / W);
164 
165  } else {
166 
167  //THROW(JDivisionByZero, "Invalid eigen value.");
168  }
169 
170  return Q;
171  }
172 
173 
174  /**
175  * Put quaternion.
176  *
177  * \param Q quaternion
178  * \param w weight
179  */
180  void put(const JQuaternion3D& Q, const double w = 1.0)
181  {
182  TMatrixDSym a(4);
183 
184  a(0,0) = Q.getA()*Q.getA(); a(0,1) = Q.getA()*Q.getB(); a(0,2) = Q.getA()*Q.getC(); a(0,3) = Q.getA()*Q.getD();
185  a(1,0) = a(0,1); a(1,1) = Q.getB()*Q.getB(); a(1,2) = Q.getB()*Q.getC(); a(1,3) = Q.getB()*Q.getD();
186  a(2,0) = a(0,2); a(2,1) = a(1,2); a(2,2) = Q.getC()*Q.getC(); a(2,3) = Q.getC()*Q.getD();
187  a(3,0) = a(0,3); a(3,1) = a(1,3); a(3,2) = a(2,3); a(3,3) = Q.getD()*Q.getD();
188 
189  a *= w;
190  A += a;
191  W += w*w;
192  }
193 
194  private:
195  TMatrixDSym A;
196  double W;
197  };
198 
199 
200  /**
201  * Numerical precision for eigen value evaluation.
202  */
204 
205 
206  /**
207  * Template specialisation for function evaluation of Legendre polynome of quaternions for undefined number of degrees.
208  */
209  template<>
210  struct JLegendre<JQuaternion3D, (size_t) -1> :
211  public JLegendre_t<JQuaternion3D>
212  {
213  /**
214  * Default constructor.
215  */
217  {}
218 
219 
220  /**
221  * Constructor.
222  *
223  * \param xmin minimal abscissa value
224  * \param xmax maximal abscissa value
225  */
226  JLegendre(const double xmin,
227  const double xmax)
228  {
229  this->set(xmin, xmax);
230  }
231 
232 
233  /**
234  * Function value.
235  *
236  * \param x abscissa value
237  * \return function value
238  */
239  virtual JQuaternion3D getValue(const double x) const override
240  {
241  const double z = this->getX(x);
242  JQuaternion3D y = zero;
243 
244  for (size_t n = 0; n != this->size(); ++n) {
245  y *= pow((*this)[n], legendre(n,z));
246  }
247 
248  return y.normalise();
249  }
250  };
251 
252 
253  /**
254  * Template specialisation for function evaluation of Legendre polynome of quaternions for defined number of degrees.
255  */
256  template<size_t N>
258  JLegendre<JQuaternion3D, (size_t) -1>
259  {
261 
262 
263  /**
264  * Default constructor.
265  */
267  {
269  }
270 
271 
272  /**
273  * Constructor.
274  *
275  * \param xmin minimal abscissa value
276  * \param xmax maximal abscissa value
277  */
278  JLegendre(const double xmin,
279  const double xmax)
280  {
282  this->set(xmin, xmax);
283  }
284 
285 
286  /**
287  * Constructor.
288  *
289  * The template argument <tt>T</tt> refers to an iterator of a data structure which should have the following data members:
290  * - double %first; // abscissa
291  * - JQuaternion3D %second; // ordinate
292  *
293  * which conforms with std::pair.
294  *
295  * \param __begin begin of data
296  * \param __end end of data
297  */
298  template<class T>
299  JLegendre(T __begin, T __end)
300  {
302  this->set(__begin, __end);
303  }
304 
305 
306  /**
307  * Set Legendre polynome of quaternions.
308  *
309  * The template argument <tt>T</tt> refers to an iterator of a data structure which should have the following data members:
310  * - double %first; // abscissa
311  * - JQuaternion3D %second; // ordinate
312  *
313  * which conforms with std::pair.
314  *
315  * \param __begin begin of data
316  * \param __end end of data
317  */
318  template<class T>
319  void set(T __begin, T __end)
320  {
321  this->set(__begin, __end, __end);
322  }
323 
324 
325  /**
326  * Set Legendre polynome of quaternions.
327  *
328  * The template argument <tt>T</tt> refers to an iterator of a data structure which should have the following data members:
329  * - double %first; // abscissa
330  * - JQuaternion3D %second; // ordinate
331  *
332  * which conforms with std::pair.
333  *
334  * \param __begin begin of data
335  * \param __not not in data
336  * \param __end end of data
337  */
338  template<class T>
339  void set(T __begin, T __not, T __end)
340  {
341  // factor to be applied as power to eigen value obtained with JAverage when degree larger than zero.
342 
343  static const double factor = 1.0 / (PI * 45.0 / 180.0);
344 
345  for (size_t n = 0; n != N + 1; ++n) {
346  (*this)[n] = JQuaternion3D();
347  }
348 
349  this->xmin = std::numeric_limits<double>::max();
350  this->xmax = std::numeric_limits<double>::lowest();
351 
352  for (T i = __begin; i != __end; ++i) {
353  if (i != __not) {
354  if (i->first < this->xmin) { this->xmin = i->first; }
355  if (i->first > this->xmax) { this->xmax = i->first; }
356  }
357  }
358 
359  for (size_t n = 0; n != N + 1; ++n) {
360 
362 
363  for (T i = __begin; i != __end; ++i) {
364 
365  if (i != __not) {
366  const double x = i->first;
367  const JQuaternion3D& y = i->second;
368  const double z = this->getX(x);
369  const double w = legendre(n, z);
370  const JQuaternion3D q = this->getValue(x).getConjugate() * y;
371 
372  Q.put(q, w);
373  }
374  }
375 
376  (*this)[n] = Q;
377 
378  if (n != 0) {
379  (*this)[n].pow(factor);
380  }
381  }
382  }
383 
384 
385  /**
386  * Read Legendre polynome from input.
387  *
388  * \param in input stream
389  * \param object Legendre polynome
390  * \return input stream
391  */
392  friend inline std::istream& operator>>(std::istream& in, JLegendre& object)
393  {
394  for (typename JLegendre::iterator i = object.begin(); i != object.end(); ++i) {
395  in >> *i;
396  }
397 
398  return in;
399  }
400 
401  private:
402  void clear();
403  void resize();
404  void erase();
405  void push_back();
406  void pop_back();
407  };
408 }
409 
410 #endif
Exceptions.
Base class for data structures with artithmetic capabilities.
Data structure for unit quaternion in three dimensions.
JQuaternion3D & pow(const double y)
Raise quaternion to given power.
double getB() const
Get b value.
double getD() const
Get d value.
double getC() const
Get c value.
JQuaternion3D & negate()
Negate quaternion.
double getA() const
Get a value.
Data structure for vector in three dimensions.
Definition: JVector3D.hh:36
double getY() const
Get y position.
Definition: JVector3D.hh:104
double getZ() const
Get z position.
Definition: JVector3D.hh:115
double getX() const
Get x position.
Definition: JVector3D.hh:94
Exception for division by zero.
Definition: JException.hh:288
const double a
Definition: JQuadrature.cc:42
const double xmax
Definition: JQuadrature.cc:24
const double xmin
Definition: JQuadrature.cc:23
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
Definition: JScale.hh:47
Auxiliary classes and methods for 3D geometrical objects and operations.
Definition: JAngle3D.hh:19
Auxiliary classes and methods for mathematical operations.
Definition: JEigen3D.hh:88
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:105
T pow(const T &x, const double y)
Power .
Definition: JMath.hh:97
static const double PI
Mathematical constants.
double legendre(const size_t n, const double x)
Legendre polynome.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
data_type w[N+1][M+1]
Definition: JPolint.hh:867
const int n
Definition: JPolint.hh:786
double u[N+1]
Definition: JPolint.hh:865
Eigen values in 3D.
Definition: JEigen3D.hh:35
JEigenValues3D(T __begin, T __end)
Constructor.
Definition: JEigen3D.hh:43
Template spacialisation for averaging quaternions.
Definition: JEigen3D.hh:98
JAverage(T __begin, T __end)
Constructor.
Definition: JEigen3D.hh:121
void put(const JQuaternion3D &Q, const double w=1.0)
Put quaternion.
Definition: JEigen3D.hh:180
static double MINIMAL_EIGEN_VALUE
Numerical precision for eigen value evaluation.
Definition: JEigen3D.hh:103
JAverage()
Default constructor.
Definition: JEigen3D.hh:108
Auxiliary class to determine average of set of values.
Definition: JMath.hh:409
void put(const JValue_t &value, const double w=1.0)
Put value.
Definition: JMath.hh:467
void set(T __begin, T __not, T __end)
Set Legendre polynome of quaternions.
Definition: JEigen3D.hh:339
friend std::istream & operator>>(std::istream &in, JLegendre &object)
Read Legendre polynome from input.
Definition: JEigen3D.hh:392
JLegendre(const double xmin, const double xmax)
Constructor.
Definition: JEigen3D.hh:278
JLegendre(T __begin, T __end)
Constructor.
Definition: JEigen3D.hh:299
void set(T __begin, T __end)
Set Legendre polynome of quaternions.
Definition: JEigen3D.hh:319
JLegendre()
Default constructor.
Definition: JEigen3D.hh:266
JLegendre(const double xmin, const double xmax)
Constructor.
Definition: JEigen3D.hh:226
virtual JQuaternion3D getValue(const double x) const override
Function value.
Definition: JEigen3D.hh:239
Base class for Legendre polynome.
Definition: JLegendre.hh:28
Template definition for function evaluation of Legendre polynome.
Definition: JLegendre.hh:269
void set(T __begin, T __end)
Set Legendre polynome.
Definition: JLegendre.hh:327