OpenMEEG
analytics.h
Go to the documentation of this file.
1 /*
2 Project Name : OpenMEEG
3 
4 © INRIA and ENPC (contributors: Geoffray ADDE, Maureen CLERC, Alexandre
5 GRAMFORT, Renaud KERIVEN, Jan KYBIC, Perrine LANDREAU, Théodore PAPADOPOULO,
6 Emmanuel OLIVI
7 Maureen.Clerc.AT.inria.fr, keriven.AT.certis.enpc.fr,
8 kybic.AT.fel.cvut.cz, papadop.AT.inria.fr)
9 
10 The OpenMEEG software is a C++ package for solving the forward/inverse
11 problems of electroencephalography and magnetoencephalography.
12 
13 This software is governed by the CeCILL-B license under French law and
14 abiding by the rules of distribution of free software. You can use,
15 modify and/ or redistribute the software under the terms of the CeCILL-B
16 license as circulated by CEA, CNRS and INRIA at the following URL
17 "http://www.cecill.info".
18 
19 As a counterpart to the access to the source code and rights to copy,
20 modify and redistribute granted by the license, users are provided only
21 with a limited warranty and the software's authors, the holders of the
22 economic rights, and the successive licensors have only limited
23 liability.
24 
25 In this respect, the user's attention is drawn to the risks associated
26 with loading, using, modifying and/or developing or reproducing the
27 software by the user in light of its specific status of free software,
28 that may mean that it is complicated to manipulate, and that also
29 therefore means that it is reserved for developers and experienced
30 professionals having in-depth computer knowledge. Users are therefore
31 encouraged to load and test the software's suitability as regards their
32 requirements in conditions enabling the security of their systems and/or
33 data to be ensured and, more generally, to use and operate it in the
34 same conditions as regards security.
35 
36 The fact that you are presently reading this means that you have had
37 knowledge of the CeCILL-B license and that you accept its terms.
38 */
39 
40 #pragma once
41 
42 #include <isnormal.H>
43 #include <mesh.h>
44 
45 namespace OpenMEEG {
46 
47  inline double integral_simplified_green(const Vect3& p0x, const double norm2p0x,
48  const Vect3& p1x, const double norm2p1x,
49  const Vect3& p1p0, const double norm2p1p0)
50  {
51  // The quantity arg is normally >= 1, verifying this relates to a triangular inequality
52  // between p0, p1 and x.
53  // Consequently, there is no need of an absolute value in the first case.
54 
55  const double arg = (norm2p0x * norm2p1p0 - p0x * p1p0) / (norm2p1x * norm2p1p0 - p1x * p1p0);
56  return (std::isnormal(arg) && arg > 0.0) ? log(arg) : fabs(log(norm2p1x / norm2p0x));
57  }
58 
59  class OPENMEEG_EXPORT analyticS
60  {
61  Vect3 p0, p1, p2;
62  Vect3 p2p1, p1p0, p0p2;
63  Vect3 nu0, nu1, nu2;
65  double norm2p2p1, norm2p1p0, norm2p0p2;
66  double tanTHETA0m, tanTHETA0p, tanTHETA1m, tanTHETA1p, tanTHETA2m, tanTHETA2p;
67 
68  void init_aux() {
69  nu0 = (p1p0^n);
70  nu1 = (p2p1^n);
71  nu2 = (p0p2^n);
72  nu0.normalize();
73  nu1.normalize();
74  nu2.normalize();
75  }
76 
77  public:
80  void init(const Triangle& T)
81  {
82  // all computations needed when the first triangle of integration is changed
83  p0 = T.s1();
84  p1 = T.s2();
85  p2 = T.s3();
86 
87  p1p0 = p1-p0; p2p1 = p2-p1; p0p2 = p0-p2;
88  norm2p1p0 = p1p0.norm(); norm2p2p1 = p2p1.norm(); norm2p0p2 = p0p2.norm();
89 
90  n = T.normal(); // since triangle's normal is already normalized
91  init_aux();
92  }
93 
94  void init(const Vect3& v0, const Vect3& v1, const Vect3& v2)
95  {
96  // all computations needed when the first triangle of integration is changed
97  p0 = v0;
98  p1 = v1;
99  p2 = v2;
100  p1p0 = p1-p0; p2p1 = p2-p1; p0p2 = p0-p2;
101  norm2p1p0 = p1p0.norm(); norm2p2p1 = p2p1.norm(); norm2p0p2 = p0p2.norm();
102 
103  n = p1p0^p0p2;
104  n /= n.norm();
105  init_aux();
106  }
107 
108  inline double f(const Vect3& x) const
109  {
110  // analytical value of the internal integral of S operator at point X
111  const Vect3& p0x = p0-x;
112  const Vect3& p1x = p1-x;
113  const Vect3& p2x = p2-x;
114  const double norm2p0x = p0x.norm();
115  const double norm2p1x = p1x.norm();
116  const double norm2p2x = p2x.norm();
117 
118  const double g0 = integral_simplified_green(p0x, norm2p0x, p1x, norm2p1x, p1p0, norm2p1p0);
119  const double g1 = integral_simplified_green(p1x, norm2p1x, p2x, norm2p2x, p2p1, norm2p2p1);
120  const double g2 = integral_simplified_green(p2x, norm2p2x, p0x, norm2p0x, p0p2, norm2p0p2);
121 
122  const double alpha = p0x*n;
123 
124  return (((p0x*nu0)*g0+(p1x*nu1)*g1+(p2x*nu2)*g2)-alpha*x.solangl(p0, p1, p2));
125  }
126  };
127 
128  class OPENMEEG_EXPORT analyticD3
129  {
130  const Vect3 &v1, &v2, &v3;
131 
132  public:
133  analyticD3(const Triangle& T): v1(T.s1()), v2(T.s2()), v3(T.s3()) {}
134 
136 
137  inline Vect3 f(const Vect3& x) const {
138  //Analytical value of the inner integral in operator D. See DeMunck article for further details.
139  // for non-optimized version of operator D
140  // returns in a vector, the inner integrals of operator D on a triangle viewed as a part of the 3
141  // P1 functions it has a part in.
142  Vect3 Y1 = v1-x;
143  Vect3 Y2 = v2-x;
144  Vect3 Y3 = v3-x;
145  double y1 = Y1.norm();
146  double y2 = Y2.norm();
147  double y3 = Y3.norm();
148  double d = Y1*(Y2^Y3);
149 
150  double derr = 1e-10;
151  if ( fabs(d) < derr ) {
152  return 0.0;
153  }
154 
155  double omega = 2. * atan2(d, (y1*y2*y3 + y1*(Y2*Y3) + y2*(Y3*Y1) + y3*(Y1*Y2)));
156 
157  Vect3 Z1 = Y2^Y3;
158  Vect3 Z2 = Y3^Y1;
159  Vect3 Z3 = Y1^Y2;
160  Vect3 D1 = Y2-Y1;
161  Vect3 D2 = Y3-Y2;
162  Vect3 D3 = Y1-Y3;
163  double d1 = D1.norm();
164  double d2 = D2.norm();
165  double d3 = D3.norm();
166  double g1 = -1.0/d1*log((y1*d1+Y1*D1)/(y2*d1+Y2*D1));
167  double g2 = -1.0/d2*log((y2*d2+Y2*D2)/(y3*d2+Y3*D2));
168  double g3 = -1.0/d3*log((y3*d3+Y3*D3)/(y1*d3+Y1*D3));
169  Vect3 N = Z1+Z2+Z3;
170  double invA = 1.0/N.norm2();
171  Vect3 S = D1*g1+D2*g2+D3*g3;
172  Vect3 omega_i;
173  omega_i.x() = invA*(Z1*N*omega+d*(D2*S));
174  omega_i.y() = invA*(Z2*N*omega+d*(D3*S));
175  omega_i.z() = invA*(Z3*N*omega+d*(D1*S));
176 
177  return omega_i;
178  }
179  };
180 
181  class OPENMEEG_EXPORT analyticDipPot
182  {
185 
186  public:
189 
190  inline void init( const Vect3& _q, const Vect3& _r0) {
191  q = _q;
192  r0 = _r0;
193  }
194 
195  inline double f(const Vect3& x) const {
196  // RK: A = q.(x-r0)/||^3
197  Vect3 r = x-r0;
198  double rn = r.norm();
199  return (q*r)/(rn*rn*rn);
200  }
201  };
202 
203  class OPENMEEG_EXPORT analyticDipPotDer
204  {
205  Vect3 q, r0;
206  Vect3 H0, H1, H2;
207  Vect3 H0p0DivNorm2, H1p1DivNorm2, H2p2DivNorm2, n;
208 
209  public:
212  inline void init( const Triangle& T, const Vect3 &_q, const Vect3& _r0) {
213  q = _q;
214  r0 = _r0;
215 
216  Vect3 p0, p1, p2, p1p0, p2p1, p0p2, p1p0n, p2p1n, p0p2n, p1H0, p2H1, p0H2;
217  p0 = T.s1();
218  p1 = T.s2();
219  p2 = T.s3();
220 
221  p1p0 = p0-p1; p2p1 = p1-p2; p0p2 = p2-p0;
222  p1p0n = p1p0; p1p0n.normalize(); p2p1n = p2p1; p2p1n.normalize(); p0p2n = p0p2; p0p2n.normalize();
223 
224  p1H0 = (p1p0*p2p1n)*p2p1n; H0 = p1H0+p1; H0p0DivNorm2 = p0-H0; H0p0DivNorm2 = H0p0DivNorm2/H0p0DivNorm2.norm2();
225  p2H1 = (p2p1*p0p2n)*p0p2n; H1 = p2H1+p2; H1p1DivNorm2 = p1-H1; H1p1DivNorm2 = H1p1DivNorm2/H1p1DivNorm2.norm2();
226  p0H2 = (p0p2*p1p0n)*p1p0n; H2 = p0H2+p0; H2p2DivNorm2 = p2-H2; H2p2DivNorm2 = H2p2DivNorm2/H2p2DivNorm2.norm2();
227 
228  n = -p1p0^p0p2;
229  n.normalize();
230  }
231 
232  inline Vect3 f(const Vect3& x) const
233  {
234  Vect3 P1part(H0p0DivNorm2*(x-H0), H1p1DivNorm2*(x-H1), H2p2DivNorm2*(x-H2));
235 
236  // RK: B = n.grad_x(A) with grad_x(A)= q/||^3 - 3r(q.r)/||^5
237  Vect3 r = x-r0;
238  double rn = r.norm();
239  double EMpart = n*(q/pow(rn, 3.)-3*(q*r)*r/pow(rn, 5.));
240 
241  return -EMpart*P1part; // RK: why - sign ?
242  }
243  };
244 }
Normal & normal()
Definition: triangle.h:117
Vect3 f(const Vect3 &x) const
Definition: analytics.h:232
void init(const Vect3 &v0, const Vect3 &v1, const Vect3 &v2)
Definition: analytics.h:94
void init(const Vect3 &_q, const Vect3 &_r0)
Definition: analytics.h:190
double norm2() const
Definition: vect3.h:97
void init(const Triangle &T, const Vect3 &_q, const Vect3 &_r0)
Definition: analytics.h:212
Triangle.
Definition: triangle.h:54
double norm() const
Definition: vect3.h:96
const Vertex & s2() const
Definition: triangle.h:110
Vect3 p2
vertices of the triangle
Definition: analytics.h:61
double & x()
Definition: vect3.h:84
double f(const Vect3 &x) const
Definition: analytics.h:108
double & z()
Definition: vect3.h:90
double f(const Vect3 &x) const
Definition: analytics.h:195
void init(const Triangle &T)
Definition: analytics.h:80
Vect3.
Definition: vect3.h:62
Vect3 f(const Vect3 &x) const
Definition: analytics.h:137
analyticD3(const Triangle &T)
Definition: analytics.h:133
void normalize()
Definition: vect3.h:143
const Vect3 & v3
Definition: analytics.h:130
const Vertex & s1() const
Definition: triangle.h:109
double integral_simplified_green(const Vect3 &p0x, const double norm2p0x, const Vect3 &p1x, const double norm2p1x, const Vect3 &p1p0, const double norm2p1p0)
Definition: analytics.h:47
const Vertex & s3() const
Definition: triangle.h:111
double & y()
Definition: vect3.h:87